From c8ee6ea8e1568bfb78f728097fe88f4f93d3bbd0 Mon Sep 17 00:00:00 2001 From: aditya-sengupta Date: Sat, 4 Feb 2023 11:36:59 -0800 Subject: [PATCH 1/2] PPM implementation --- src/MOL_discretization.jl | 6 +- src/MethodOfLines.jl | 3 +- .../generate_finite_difference_rules.jl | 4 + src/discretization/schemes/PPM/PPM.jl | 102 ++++++++++++++++++ .../centered_diff_weights.jl | 4 +- .../schemes/fornberg_calculate_weights.jl | 2 +- .../schemes/half_offset_weights.jl | 13 ++- src/interface/MOLFiniteDifference.jl | 2 +- src/interface/scheme_types.jl | 9 ++ src/scalar_discretization.jl | 15 ++- test/demo_parker_singular.jl | 56 ++++++++++ test/demo_ppm_runs.jl | 31 ++++++ 12 files changed, 235 insertions(+), 12 deletions(-) create mode 100644 src/discretization/schemes/PPM/PPM.jl create mode 100644 test/demo_parker_singular.jl create mode 100644 test/demo_ppm_runs.jl diff --git a/src/MOL_discretization.jl b/src/MOL_discretization.jl index 196c8b6d9..6369f9535 100644 --- a/src/MOL_discretization.jl +++ b/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.")) @@ -116,7 +116,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 diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl index 1f85a1eb4..382fb163d 100644 --- a/src/MethodOfLines.jl +++ b/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 diff --git a/src/discretization/generate_finite_difference_rules.jl b/src/discretization/generate_finite_difference_rules.jl index 379be5a89..047ac07d9 100644 --- a/src/discretization/generate_finite_difference_rules.jl +++ b/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 diff --git a/src/discretization/schemes/PPM/PPM.jl b/src/discretization/schemes/PPM/PPM.jl new file mode 100644 index 000000000..2c382e818 --- /dev/null +++ b/src/discretization/schemes/PPM/PPM.jl @@ -0,0 +1,102 @@ +""" +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::Number) + j, x = jx + I1 = unitindex(ndims(u, s), j) + + udisc = s.discvars[u] + + 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) + is = map(I -> I[j], [Im2, Im1, Ip1, Ip2]) + for i in is + if i < 1 + return nothing + elseif i > length(s, x) + return nothing + end + end + + 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) + return (ap - am) / dx +end + +function ppm(II::CartesianIndex, s::DiscreteSpace, b, jx, u, dx::AbstractVector) + j, x = jx + I1 = unitindex(ndims(u, s), j) + + udisc = s.discvars[u] + + 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) + is = map(I -> I[j], [Im2, Im1, Ip1, Ip2]) + for i in is + if i < 1 + return nothing + elseif i > length(s, x) + return nothing + end + end + + 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) + return (ap - am) / dx[II] +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 \ No newline at end of file diff --git a/src/discretization/schemes/centered_difference/centered_diff_weights.jl b/src/discretization/schemes/centered_difference/centered_diff_weights.jl index 2b6ad7526..ec3e90d29 100644 --- a/src/discretization/schemes/centered_difference/centered_diff_weights.jl +++ b/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 diff --git a/src/discretization/schemes/fornberg_calculate_weights.jl b/src/discretization/schemes/fornberg_calculate_weights.jl index 3aac12116..dcf8472af 100644 --- a/src/discretization/schemes/fornberg_calculate_weights.jl +++ b/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 diff --git a/src/discretization/schemes/half_offset_weights.jl b/src/discretization/schemes/half_offset_weights.jl index 86e4125e7..8b8727379 100644 --- a/src/discretization/schemes/half_offset_weights.jl +++ b/src/discretization/schemes/half_offset_weights.jl @@ -1,9 +1,10 @@ +using Unitful """ A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the half offset centered scheme. See table 2 in https://web.njit.edu/~jiang/math712/fornberg.pdf """ 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 +22,15 @@ 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 = nothing + if T <: Unitful.Quantity + half = dx / dx.val / 2.0 + else + half = convert(T, 0.5) + end 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 +62,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) diff --git a/src/interface/MOLFiniteDifference.jl b/src/interface/MOLFiniteDifference.jl index 5f354fdce..1562d4bfe 100644 --- a/src/interface/MOLFiniteDifference.jl +++ b/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 diff --git a/src/interface/scheme_types.jl b/src/interface/scheme_types.jl index 3d7c9035f..9bd8b424e 100644 --- a/src/interface/scheme_types.jl +++ b/src/interface/scheme_types.jl @@ -26,4 +26,13 @@ function extent(::WENOScheme, dorder) return 2 end +struct PPMScheme <: AbstractScheme + dt +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 diff --git a/src/scalar_discretization.jl b/src/scalar_discretization.jl index 2ecb1f840..f4e3df4c2 100644 --- a/src/scalar_discretization.jl +++ b/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 diff --git a/test/demo_parker_singular.jl b/test/demo_parker_singular.jl new file mode 100644 index 000000000..8ade431a6 --- /dev/null +++ b/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()) \ No newline at end of file diff --git a/test/demo_ppm_runs.jl b/test/demo_ppm_runs.jl new file mode 100644 index 000000000..1f8826188 --- /dev/null +++ b/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) From 6a688b7d84bb4c53d842b037b5430632f640b69c Mon Sep 17 00:00:00 2001 From: aditya-sengupta Date: Wed, 8 Feb 2023 15:39:53 -0800 Subject: [PATCH 2/2] Use full version of PPM with mid-cell interpolation --- src/discretization/schemes/PPM/PPM.jl | 119 +++++++++++------- .../schemes/half_offset_weights.jl | 8 +- src/interface/scheme_types.jl | 4 +- 3 files changed, 77 insertions(+), 54 deletions(-) diff --git a/src/discretization/schemes/PPM/PPM.jl b/src/discretization/schemes/PPM/PPM.jl index 2c382e818..b8c5aa864 100644 --- a/src/discretization/schemes/PPM/PPM.jl +++ b/src/discretization/schemes/PPM/PPM.jl @@ -1,26 +1,4 @@ -""" -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::Number) - j, x = jx - I1 = unitindex(ndims(u, s), j) - - udisc = s.discvars[u] - - 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) - is = map(I -> I[j], [Im2, Im1, Ip1, Ip2]) - for i in is - if i < 1 - return nothing - elseif i > length(s, x) - return nothing - end - end - +function ppm_interp(Im2, Im1, II, Ip1, Ip2, udisc, dx::Number) u_m2 = udisc[Im2] u_m1 = udisc[Im1] u_0 = udisc[II] @@ -29,28 +7,10 @@ function ppm(II::CartesianIndex, s::DiscreteSpace, ppmscheme::PPMScheme, bs, jx, 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) - return (ap - am) / dx -end - -function ppm(II::CartesianIndex, s::DiscreteSpace, b, jx, u, dx::AbstractVector) - j, x = jx - I1 = unitindex(ndims(u, s), j) - - udisc = s.discvars[u] - - 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) - is = map(I -> I[j], [Im2, Im1, Ip1, Ip2]) - for i in is - if i < 1 - return nothing - elseif i > length(s, x) - return nothing - end - end + 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] @@ -94,9 +54,78 @@ function ppm(II::CartesianIndex, s::DiscreteSpace, b, jx, u, dx::AbstractVector) 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) - return (ap - am) / dx[II] + 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 \ No newline at end of file diff --git a/src/discretization/schemes/half_offset_weights.jl b/src/discretization/schemes/half_offset_weights.jl index 8b8727379..462791aba 100644 --- a/src/discretization/schemes/half_offset_weights.jl +++ b/src/discretization/schemes/half_offset_weights.jl @@ -1,4 +1,3 @@ -using Unitful """ A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the half offset centered scheme. See table 2 in https://web.njit.edu/~jiang/math712/fornberg.pdf """ @@ -22,12 +21,7 @@ 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 = nothing - if T <: Unitful.Quantity - half = dx / dx.val / 2.0 - else - half = convert(T, 0.5) - end + half = convert(T, 0.5) stencil_coefs = convert(SVector{centered_stencil_length, T}, (1 / dx^derivative_order) * calculate_weights(derivative_order, half, dummy_x)) diff --git a/src/interface/scheme_types.jl b/src/interface/scheme_types.jl index 9bd8b424e..4866f0ade 100644 --- a/src/interface/scheme_types.jl +++ b/src/interface/scheme_types.jl @@ -26,8 +26,8 @@ function extent(::WENOScheme, dorder) return 2 end -struct PPMScheme <: AbstractScheme - dt +struct PPMScheme <: AbstractScheme + Cx end function extent(::PPMScheme, dorder)