add Piecewise Parabolic Method discretization method

src/MOL_discretization.jl

@@ -4,8 +4,8 @@ function interface_errors(depvars, indvars, discretization)
     for x in indvars
         @assert haskey(discretization.dxs, Num(x)) || haskey(discretization.dxs, x) "Variable $x has no step size"
     end
-    if !(typeof(discretization.advection_scheme) ∈ [UpwindScheme, WENOScheme])
-        throw(ArgumentError("Only `UpwindScheme()` and `WENOScheme()` are supported advection schemes."))
+    if !(typeof(discretization.advection_scheme) ∈ [UpwindScheme, WENOScheme, PPMScheme])
+        throw(ArgumentError("Only `UpwindScheme()`, `PPMScheme()`, and `WENOScheme()` are supported advection schemes."))
     end
     if !(typeof(discretization.disc_strategy) ∈ [ScalarizedDiscretization])
         throw(ArgumentError("Only `ScalarizedDiscretization()` are supported discretization strategies."))
@@ -112,7 +112,7 @@ function SciMLBase.symbolic_discretize(pdesys::PDESystem, discretization::Method
         if disc_strategy isa ScalarizedDiscretization
             # Generate the equations for the interior points
             discretize_equation!(alleqs, bceqs, pde, interiormap, eqvar, bcmap,
-                depvars, s, derivweights, indexmap)
+                depvars, s, derivweights, indexmap, get(discretization.kwargs, :show_symbolic, false))
         else
             throw(ArgumentError("Only ScalarizedDiscretization is currently supported"))
         end

src/MethodOfLines.jl

@@ -66,6 +66,7 @@ include("discretization/schemes/half_offset_centred_difference.jl")
 include("discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl")
 include("discretization/schemes/spherical_laplacian/spherical_laplacian.jl")
 include("discretization/schemes/WENO/WENO.jl")
+include("discretization/schemes/PPM/PPM.jl")
 include("discretization/schemes/integral_expansion/integral_expansion.jl")
 
 # System Discretization
@@ -79,6 +80,6 @@ include("scalar_discretization.jl")
 include("MOL_discretization.jl")
 
 export MOLFiniteDifference, discretize, symbolic_discretize, ODEFunctionExpr, generate_code, grid_align, edge_align, center_align, get_discrete
-export UpwindScheme, WENOScheme
+export UpwindScheme, WENOScheme, PPMScheme
 
 end

src/discretization/generate_finite_difference_rules.jl

@@ -58,6 +58,10 @@ function generate_finite_difference_rules(II::CartesianIndex, s::DiscreteSpace,
             advection_rules = generate_WENO_rules(II, s, depvars, derivweights, bmap, indexmap, terms)
             advection_rules = vcat(advection_rules,
                                 generate_winding_rules(II, s, depvars, derivweights, bmap, indexmap, terms; skip = [1]))
+        elseif derivweights.advection_scheme isa PPMScheme
+            advection_rules = generate_PPM_rules(II, s, depvars, derivweights, bmap, indexmap, terms)
+            advection_rules = vcat(advection_rules,
+                                generate_winding_rules(II, s, depvars, derivweights, bmap, indexmap, terms; skip = [1]))
         else
             error("Unsupported advection scheme $(derivweights.advection_scheme) encountered.")
         end

src/discretization/schemes/PPM/PPM.jl

@@ -0,0 +1,131 @@
+function ppm_interp(Im2, Im1, II, Ip1, Ip2, udisc, dx::Number) 
+    u_m2 = udisc[Im2]
+    u_m1 = udisc[Im1]
+    u_0 = udisc[II]
+    u_p1 = udisc[Ip1]
+    u_p2 = udisc[Ip2]
+
+    ap = (7/12) * (u_0 + u_p1) - (1/12) * (u_p2 - u_m1)
+    am = (7/12) * (u_m1 + u_0) - (1/12) * (u_p1 - u_m2)
+    am, ap, u_0
+end
+
+function ppm_interp(Im2, Im1, II, Ip1, Ip2, udisc, dx::AbstractVector) 
+    u_m2 = udisc[Im2]
+    u_m1 = udisc[Im1]
+    u_0 = udisc[II]
+    u_p1 = udisc[Ip1]
+    u_p2 = udisc[Ip2]
+
+    # CW 1.7
+    d_a0 = (dx[II] / (dx[Im1] + dx[II] + dx[Ip1])) * ((2 * dx[Im1] + dx[II]) * (u_p1 - u_0) / (dx[Ip1] + dx[II]) + (dx[II] + 2 * dx[Ip1]) * (u_0 - u_m1) / (dx[Im1] + dx[II]))
+    d_am1 = (dx[Im1] / (dx[Im2] + dx[Im1] + dx[II])) * ((2 * dx[Im2] + dx[Im1]) * (u_0 - u_m1) / (dx[II] + dx[Im1]) + (dx[Im1] + 2 * dx[II]) * (u_m1 - u_m2) / (dx[Im2] + dx[Im1]))
+
+    # CW 1.8
+    ud_p1 = u_p1 - u_0
+    ud_0 = u_0 - u_m1
+    ud_m1 = u_m1 - u_m2
+    if sign(ud_p1) == sign(ud_0)
+        d_a0 = min(abs(d_a0), min(2 * abs(ud_p1), 2 * abs(ud_0))) * sign(d_a0)
+    else
+        d_a0 = 0
+    end
+    if sign(ud_0) == sign(ud_m1)
+        d_am1 = min(abs(d_am1), min(2 * abs(ud_0), 2 * abs(ud_m1))) * sign(d_am1)
+    else
+        d_am1 = 0
+    end
+
+    coeffp1 = dx[II] / (dx[II] + dx[Ip1])
+    coeffp2 = 1 / (dx[Im1] + dx[II] + dx[Ip1] + dx[Ip2])
+    coeffp3 = (2 * dx[Ip1] * dx[II]) / (dx[II] + dx[Ip1])
+    coeffp4 = (dx[Im1] + dx[II]) / (2 * dx[II] + dx[Ip1])
+    coeffp5 = (dx[Ip2] + dx[Ip1]) / (2 * dx[Ip1] + dx[II])
+    coeffp6 = dx[II] * (dx[Im1] + dx[II]) / (2 * dx[II] + dx[Ip1])
+    coeffp7 = dx[Ip1] * (dx[Ip1] + dx[Ip2]) / (dx[II] + 2 * dx[Ip1])
+
+    coeffm1 = dx[Im1] / (dx[Im1] + dx[II])
+    coeffm2 = 1 / (dx[Im2] + dx[Im1] + dx[II] + dx[Ip1])
+    coeffm3 = (2 * dx[II] * dx[Im1]) / (dx[Im1] + dx[II])
+    coeffm4 = (dx[Im2] + dx[Im1]) / (2 * dx[Im1] + dx[II])
+    coeffm5 = (dx[Ip1] + dx[II]) / (2 * dx[II] + dx[Im1])
+    coeffm6 = dx[Im1] * (dx[Im2] + dx[Im1]) / (2 * dx[Im1] + dx[II])
+    coeffm7 = dx[II] * (dx[II] + dx[Ip1]) / (dx[Im1] + 2 * dx[II])
+
+    ap = u_0 + coeffp1 * (u_p1 - u_0) + coeffp2 * (coeffp3 * (coeffp4 - coeffp5) * (u_p1 - u_0) - coeffp6 * ud_p1 + coeffp7 * ud_0)
+    am = u_m1 + coeffm1 * (u_0 - u_m1) + coeffm2 * (coeffm3 * (coeffm4 - coeffm5) * (u_0 - u_m1) - coeffm6 * ud_0 + coeffm7 * ud_m1)
+    am, ap, u_0
+end
+
+"""
+Corrects values between zone boundaries as per CW84 1.10
+Slightly inefficient in the case that cond is true, because ifelse can't handle symbolic tuples very well
+"""
+function ppm_zone_boundary_correct(am, ap, a0)
+    C1 = (ap - am) * (a0 - (am + ap)/2)
+    C2 = (ap - am)^2 / 6
+    cond = sign(ap - a0) != sign(a0 - am)
+    aL = ifelse(cond, a0, ifelse(C1 > C2, 3*a0 - 2*ap, am))
+    aR = ifelse(cond, a0, ifelse(-C2 > C1, 3*a0 - 2*am, ap))
+    aL, aR
+end
+
+"""
+Computes the spatial derivative as per CW84 1.11-1.13.
+True PPM requires knowledge of u and Δt, which we don't have;
+instead, we use the Courant number as an interpolation hyperparameter
+and always use the positive-advection version of PPM.
+"""
+function ppm_du(aL, aR, a0, courant)
+    a6 = 6 * (a0 - (aL + aR) / 2)
+    return aR - (courant / 2) * (aR - aL - (1 - 2 * courant / 3) * a6)
+end
+
+function ppm_dudx(am, ap, ap2, a0, a1, Cx, dx::Number, II, Ip1)
+    return (ppm_du(ap, ap2, a1, Cx / dx) - ppm_du(am, ap, a0, Cx / dx)) / dx
+end
+
+function ppm_dudx(am, ap, ap2, a0, a1, Cx, dx::AbstractVector, II, Ip1)
+    return (ppm_du(ap, ap2, a1, Cx / dx[Ip1]) - ppm_du(am, ap, a0, Cx / dx[II])) / dx[II]
+end
+
+"""
+Implements the piecewise parabolic method for fixed dx.
+Uses Equation 1.9 from Collela and Woodward, 1984. 
+"""
+function ppm(II::CartesianIndex, s::DiscreteSpace, ppmscheme::PPMScheme, bs, jx, u, dx::Union{Number,AbstractVector})
+    j, x = jx
+    I1 = unitindex(ndims(u, s), j)
+
+    Im2 = bwrap(II - 2I1, bs, s, jx)
+    Im1 = bwrap(II - I1, bs, s, jx)
+    Ip1 = bwrap(II + I1, bs, s, jx)
+    Ip2 = bwrap(II + 2I1, bs, s, jx)
+    Ip3 = bwrap(II + 3I1, bs, s, jx)
+    is = map(I -> I[j], [Im1, II, Ip1, Ip2])
+    for i in is
+        if i < 1
+            return nothing
+        elseif i > length(s, x)
+            return nothing
+        end
+    end
+
+    if Ip3[j] > length(s, x)
+        Ip3 = Ip2
+    end
+    if Im2[j] < 1
+        Im2 = Im1
+    end
+
+    am, ap, a0 = ppm_interp(Im2, Im1, II, Ip1, Ip2, s.discvars[u], dx) # dispatches to number or abstractvector
+    _, ap2, a1 = ppm_interp(Im1, II, Ip1, Ip2, Ip3, s.discvars[u], dx)
+    am, ap = ppm_zone_boundary_correct(am, ap, a0)
+    _, ap2 = ppm_zone_boundary_correct(ap, ap2, a1)
+
+    ppm_dudx(am, ap, ap2, a0, a1, ppmscheme.Cx, dx, II, Ip1)
+end 
+
+function generate_PPM_rules(II::CartesianIndex, s::DiscreteSpace, depvars, derivweights::DifferentialDiscretizer, bcmap, indexmap, terms)
+    return reduce(safe_vcat, [[(Differential(x))(u) => ppm(Idx(II, s, u, indexmap), s, derivweights.advection_scheme, filter_interfaces(bcmap[operation(u)][x]), (x2i(s, u, x), x), u, s.dxs[x]) for x in params(u, s)] for u in depvars], init = [])
+end

src/discretization/schemes/centered_difference/centered_diff_weights.jl

@@ -2,7 +2,7 @@
 A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the centered scheme.
 """
 function CompleteCenteredDifference(derivative_order::Int,
-                                    approximation_order::Int, dx::T) where {T <: Real}
+                                    approximation_order::Int, dx::T) where T
     @assert approximation_order>1 "approximation_order must be greater than 1."
     stencil_length = derivative_order + approximation_order - 1 +
                      (derivative_order + approximation_order) % 2
@@ -54,7 +54,7 @@ end
 
 function CompleteCenteredDifference(derivative_order::Int,
                                     approximation_order::Int,
-                                    x::AbstractVector{T}) where {T <: Real}
+                                    x::AbstractVector{T}) where T
     stencil_length = derivative_order + approximation_order - 1 +
                      (derivative_order + approximation_order) % 2
     boundary_stencil_length = derivative_order + approximation_order

src/discretization/schemes/fornberg_calculate_weights.jl

@@ -17,7 +17,7 @@ Inputs:
                                                  derivative values respectively.                                             
 =#
 
-function calculate_weights(order::Int, x0::T, x::AbstractVector; dfdx::Bool = false) where T<:Real
+function calculate_weights(order::Int, x0::T, x::AbstractVector; dfdx::Bool = false) where T
     N = length(x)
     @assert order < N "Not enough points for the requested order."
     M = order

src/discretization/schemes/half_offset_weights.jl

@@ -3,7 +3,7 @@ A helper function to compute the coefficients of a derivative operator including
 """
 function CompleteHalfCenteredDifference(derivative_order::Int,
                                         approximation_order::Int,
-                                        dx::T) where {T <: Real}
+                                        dx::T) where T
     @assert approximation_order>1 "approximation_order must be greater than 1."
     centered_stencil_length = approximation_order + 2 * Int(floor(derivative_order / 2)) +
                               (approximation_order % 2)
@@ -21,9 +21,10 @@ function CompleteHalfCenteredDifference(derivative_order::Int,
     #deriv_spots             = (-div(stencil_length,2)+1) : -1  # unused
     L_boundary_deriv_spots = xoffset[1:div(centered_stencil_length, 2)]
 
+    half = convert(T, 0.5)
     stencil_coefs = convert(SVector{centered_stencil_length, T},
                             (1 / dx^derivative_order) *
-                            calculate_weights(derivative_order, convert(T, 0.5), dummy_x))
+                            calculate_weights(derivative_order, half, dummy_x))
     # For each boundary point, for each tappoint in the half offset central difference stencil, we need to calculate the coefficients for the stencil.
 
     _low_boundary_coefs = [convert(SVector{boundary_stencil_length, T},
@@ -55,7 +56,7 @@ end
 
 function CompleteHalfCenteredDifference(derivative_order::Int,
                                         approximation_order::Int,
-                                        x::T) where {T <: AbstractVector{<:Real}}
+                                        x::T) where {T <: AbstractVector}
     @assert approximation_order>1 "approximation_order must be greater than 1."
     centered_stencil_length = approximation_order + 2 * Int(floor(derivative_order / 2)) +
                               (approximation_order % 2)

src/interface/MOLFiniteDifference.jl

@@ -49,5 +49,5 @@ function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, advection_sche
 
     dxs = dxs isa Dict ? dxs : Dict(dxs)
 
-    return MOLFiniteDifference{typeof(grid_align), typeof(discretization_strategy)}(dxs, time, approx_order, advection_scheme, grid_align, should_transform, use_ODAE, discretization_strategy, kwargs)
+    return MOLFiniteDifference{typeof(grid_align), typeof(discretization_strategy)}(dxs, time, approx_order, advection_scheme, grid_align, should_transform, use_ODAE, discretization_strategy, Dict(kwargs))
 end

src/interface/scheme_types.jl

@@ -26,4 +26,13 @@ function extent(::WENOScheme, dorder)
     return 2
 end
 
+struct PPMScheme <: AbstractScheme 
+    Cx
+end
+
+function extent(::PPMScheme, dorder)
+    @assert dorder == 1 "PPM requires first-order derivatives." # I think
+    return 2
+end
+
 # Note: This type and its subtypes will become important later with the stencil interfaces as we will need to dispatch on derivative order and approximation order

src/scalar_discretization.jl

@@ -1,10 +1,18 @@
-function discretize_equation!(alleqs, bceqs, pde, interiormap, eqvar, bcmap, depvars, s, derivweights, indexmap)
+function discretize_equation!(alleqs, bceqs, pde, interiormap, eqvar, bcmap, depvars, s, derivweights, indexmap, show_symbolic)
     # Handle boundary values appearing in the equation by creating functions that map each point on the interior to the correct replacement rule
     # Generate replacement rule gen closures for the boundary values like u(t, 1)
     boundaryvalfuncs = generate_boundary_val_funcs(s, depvars, bcmap, indexmap, derivweights)
     # Find boundaries for this equation
     eqvarbcs = mapreduce(x -> bcmap[operation(eqvar)][x], vcat, s.x̄)
     # Generate the boundary conditions for the correct variable
+
+    if show_symbolic
+        isyms = @. Symbol("i_" * string(unwrap(s.vars.x̄)) * "::Int")
+        symindices = map(sym -> Sym{Int, Nothing}(sym, nothing), isyms)
+        idx = CartesianIndex(symindices...) |> eval
+        println(idx, typeof(idx))
+        println(discretize_equation_at_point(idx, s, depvars, pde, derivweights, bcmap, eqvar, indexmap, boundaryvalfuncs))
+    end
     for boundary in eqvarbcs
         generate_bc_eqs!(bceqs, s, boundaryvalfuncs, interiormap, boundary)
     end
@@ -14,6 +22,11 @@ function discretize_equation!(alleqs, bceqs, pde, interiormap, eqvar, bcmap, dep
     # Set invalid corner points to zero
     generate_corner_eqs!(bceqs, s, interiormap, ndims(s.discvars[eqvar]), eqvar)
     # Extract Interior
+    generate_interior_eqs!(alleqs, interiormap, s, depvars, pde, derivweights, bcmap, eqvar, indexmap, boundaryvalfuncs)
+end
+
+# this should live in discretization/generate_interior_eqs.jl for consistency
+function generate_interior_eqs!(alleqs, interiormap, s, depvars, pde, derivweights, bcmap, eqvar, indexmap, boundaryvalfuncs)
     interior = interiormap.I[pde]
     # Generate the discrete form ODEs for the interior
     eqs = if length(interior) == 0

test/demo_parker_singular.jl

@@ -0,0 +1,56 @@
+using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets
+using Symbolics: @register_symbolic
+
+@parameters x t
+@variables ρ(..) v(..)
+Dt = Differential(t)
+Dx = Differential(x)
+
+GM = 4e20
+cs = 2.87e5
+cs2 = 8.24e10
+
+eq = [
+        Dt(ρ(x,t)) ~ -exp(-x) * (2*ρ(x,t)*v(x,t) + v(x,t)*Dx(ρ(x,t)) + ρ(x,t)*Dx(v(x,t))),
+        Dt(v(x,t)) ~ -v(x,t)*exp(-x) * Dx(v(x,t)) - cs2/(ρ(x,t)) * exp(-x) * Dx(ρ(x,t)) - GM*exp(-2*x)
+    ]
+
+x_min, x_max = 20.2, 21.8
+t_min, t_max = 0.0, 30_000.0
+
+domains = [x ∈ Interval(x_min, x_max), t ∈ Interval(t_min, t_max)]
+
+ρ_init(x) = 5e-15 * exp(4.84e9 * exp(-x))
+function v_init(x)
+    if (x > 21.6)
+        return 2 * cs
+    else
+        return 0.0
+    end
+end
+
+@register_symbolic v_init(x)
+
+bcs = [
+    ρ(x,t_min) ~ ρ_init(x),
+    ρ(x_min,t) ~ ρ_init(x_min),
+    v(x,t_min) ~ v_init(x),
+    v(x_min,t) ~ v_init(x_min),
+    Dx(v(x_min,t)) ~ 0.0,
+    Dx(ρ(x_min,t)) ~ 0.0,
+    Dx(v(x_max,t)) ~ 0.0,
+    Dx(ρ(x_max,t)) ~ 0.0,
+    #Dt(v(x,t_min)) ~ 0.0,
+    #Dt(ρ(x,t_min)) ~ 0.0,
+] 
+
+@named pdesys = PDESystem(eq, bcs, domains, [x,t], [ρ(x,t),v(x,t)])
+
+N = 20
+dx = (x_max-x_min)/N
+order = 2
+discretization = MOLFiniteDifference([x=>dx], t, advection_scheme=WENOScheme())
+
+prob = discretize(pdesys, discretization)
+
+sol = solve(prob, ImplicitEuler())

test/demo_ppm_runs.jl

@@ -0,0 +1,31 @@
+using MethodOfLines, DomainSets, OrdinaryDiffEq, ModelingToolkit
+
+@parameters x t
+@variables u(..)
+Dx = Differential(x)
+Dt = Differential(t)
+x_min = 0.0
+x_max = 1.0
+t_min = 0.0
+t_max = 6.0
+
+analytic_u(t, x) = x / (t + 1)
+
+eq = Dt(u(t, x)) ~ -u(t, x) * Dx(u(t, x))
+
+bcs = [u(0, x) ~ x,
+    u(t, x_min) ~ analytic_u(t, x_min),
+    u(t, x_max) ~ analytic_u(t, x_max)]
+
+domains = [t ∈ Interval(t_min, t_max),
+    x ∈ Interval(x_min, x_max)]
+
+dx = 0.05
+
+@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x)])
+
+disc = MOLFiniteDifference([x => dx], t, advection_scheme=PPMScheme(0.1))
+
+prob = discretize(pdesys, disc)
+
+sol = solve(prob, Euler(), dt=0.1)