From 8f1a8d8e2bac0b97ec524b6da7a2cac22381ee75 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Tue, 5 Dec 2023 18:56:18 +0000 Subject: [PATCH 1/4] add complex docs --- Project.toml | 2 +- docs/pages.jl | 1 + docs/src/tutorials/schroedinger.md | 48 ++++++++++++++++++++++++++++++ 3 files changed, 50 insertions(+), 1 deletion(-) create mode 100644 docs/src/tutorials/schroedinger.md diff --git a/Project.toml b/Project.toml index 464e8a9cf..1a631d54d 100644 --- a/Project.toml +++ b/Project.toml @@ -31,7 +31,7 @@ Interpolations = "0.14" Latexify = "0.15, 0.16" ModelingToolkit = "8" OrdinaryDiffEq = "6" -PDEBase = "0.1.6" +PDEBase = "0.1.8" PrecompileTools = "1" RuntimeGeneratedFunctions = "0.5.9" SafeTestsets = "0.0.1" diff --git a/docs/pages.jl b/docs/pages.jl index 11caa32a2..cd5c5ded9 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -8,6 +8,7 @@ pages = [ "tutorials/sispde.md", "tutorials/icbc_sampled.md", "tutorials/PIDE.md", + "tutorials/schroedinger.md", ], "Performance Tips" => "performance.md", "MOLFiniteDifference" => "MOLFiniteDifference.md", diff --git a/docs/src/tutorials/schroedinger.md b/docs/src/tutorials/schroedinger.md new file mode 100644 index 000000000..72ca4b13a --- /dev/null +++ b/docs/src/tutorials/schroedinger.md @@ -0,0 +1,48 @@ +# Schrödinger Equation + +MethodOfLines can solve linear complex PDEs like the Scrödinger equation: + +```@example schro +using MethodOfLines, OrdinaryDiffEq, Plots, DomainSets, ModelingToolkit + +@parameters t, x +@variables ψ(..) + +Dt = Differential(t) +Dxx = Differential(x)^2 + +xmin = 0 +xmax = 1 + +V(x) = 0.0 + +eq = [im*Dt(ψ(t,x)) ~ Dxx(ψ(t,x)) + V(x)*ψ(t,x)] # You must enclose complex equations in a vector, even if there is only one equation + +ψ0 = x -> sin(2pi*x) + +bcs = [ψ(0,x) ~ ψ0(x), + ψ(t,xmin) ~ 0, + ψ(t,xmax) ~ 0] + +domains = [t ∈ Interval(0, 1), x ∈ Interval(xmin, xmax)] + +@named sys = PDESystem(eq, bcs, domains, [t, x], [ψ(t,x)]) + +disc = MOLFiniteDifference([x => 100], t) + +prob = discretize(sys, disc) + +sol = solve(prob, TRBDF2(), saveat = 0.01) + +discx = sol[x] +disct = sol[t] + +discψ = sol[ψ(t, x)] +anim = @animate for i in 1:length(disct) + u = discψ[i, :] + plot(discx, [real.(u), imag.(u)], ylim = (-1.5, 1.5), title = "t = $(disct[i])", xlabel = "x", ylabel = "ψ(t,x)", label = ["re(ψ)" "im(ψ)"], legend = :topleft) +end +gif(anim, "schroedinger.gif", fps = 10) +``` + +This represents the second from ground state of a particle in an infinite quantum well, try changing the potential `V(x)`, initial conditions and BCs, it is extremely interesting to see how the wave function evolves even for nonphysical combinations. Be sure to post interesting results on the discourse! \ No newline at end of file From 95f1cc844b6ff9452db7f3bcd31026537e849104 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Tue, 5 Dec 2023 19:03:36 +0000 Subject: [PATCH 2/4] remove show --- src/interface/solution/timedep.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/interface/solution/timedep.jl b/src/interface/solution/timedep.jl index ecfb21012..068df285f 100644 --- a/src/interface/solution/timedep.jl +++ b/src/interface/solution/timedep.jl @@ -53,7 +53,6 @@ function SciMLBase.PDETimeSeriesSolution(sol::SciMLBase.AbstractODESolution{T}, end |> Dict # Build Interpolations interp = build_interpolation(umap, dvs, ivs, ivgrid, sol, pdesys, discretespace.vars.replaced_vars) - @show metadata.complexmap return SciMLBase.PDETimeSeriesSolution{T,length(discretespace.ū),typeof(umap),typeof(metadata), typeof(sol),typeof(sol.errors),typeof(sol.t),typeof(ivgrid), From de97bf3ef6e104d2c5de2701bedbc5c954d341b7 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Tue, 5 Dec 2023 20:06:15 +0000 Subject: [PATCH 3/4] add cpx test --- .github/workflows/CI.yml | 1 + test/pde_systems/schroedinger.jl | 55 ++++++++++++++++++++++++++++++++ test/runtests.jl | 6 ++++ 3 files changed, 62 insertions(+) create mode 100644 test/pde_systems/schroedinger.jl diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index f87b2bd42..4345fb509 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -34,6 +34,7 @@ jobs: - Burgers - Brusselator - Mixed_Derivatives + - Complex version: - '1' - '1.6' diff --git a/test/pde_systems/schroedinger.jl b/test/pde_systems/schroedinger.jl new file mode 100644 index 000000000..dc8cb78a8 --- /dev/null +++ b/test/pde_systems/schroedinger.jl @@ -0,0 +1,55 @@ +using MethodOfLines, OrdinaryDiffEq, DomainSets, ModelingToolkit, Test + +@parameters t, x +@variables ψ(..) + +Dt = Differential(t) +Dxx = Differential(x)^2 + +xmin = 0 +xmax = 1 + +V(x) = 0.0 + +eq = [im*Dt(ψ(t,x)) ~ (Dxx(ψ(t,x)) + V(x)*ψ(t,x))] # You must enclose complex equations in a vector, even if there is only one equation + +ψ0 = x -> sin(2*pi*x) + +bcs = [ψ(0,x) ~ ψ0(x), + ψ(t,xmin) ~ 0, + ψ(t,xmax) ~ 0] + +domains = [t ∈ Interval(0, 1), x ∈ Interval(xmin, xmax)] + +@named sys = PDESystem(eq, bcs, domains, [t, x], [ψ(t,x)]) + +disc = MOLFiniteDifference([x => 100], t) + +prob = discretize(sys, disc) + +sol = solve(prob, TRBDF2(), saveat = 0.01) + +discx = sol[x] +disct = sol[t] + +discψ = sol[ψ(t, x)] + +analytic(t, x) = sqrt(2)*sin(2*pi*x)*exp(-im*4*pi^2*t) + +analψ = [analytic(t, x) for t in disct, x in discx] + +for i in 1:length(disct) + u = abs.(analψ[i, :]).^2 + u2 = abs.(discψ[i, :]).^2 + + @test u./maximum(u) ≈ u2./maximum(u2) atol=1e-3 +end + +#using Plots + +# anim = @animate for i in 1:length(disct) +# u = analψ[i, :] +# u2 = discψ[i, :] +# plot(discx, [real.(u), imag.(u)], ylim = (-1.5, 1.5), title = "t = $(disct[i])", xlabel = "x", ylabel = "ψ(t,x)", label = ["re(ψ)" "im(ψ)"], legend = :topleft) +# end +# gif(anim, "schroedinger.gif", fps = 10) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 263692a8f..368937d93 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,6 +8,12 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") # Start Test Script @time begin + if GROUP == "All" || GROUP == "Complex" + @time @safetestset "MOLFiniteDifference Interface" begin + include("pde_systems/schroedinger.jl") + end + end + if GROUP == "All" || GROUP == "Mixed_Derivatives" @time @safetestset "MOLFiniteDifference Interface: Mixed Derivatives" begin include("pde_systems/MOL_Mixed_Deriv.jl") From a06da5e7b881ea3619a09285d391d585cae759a4 Mon Sep 17 00:00:00 2001 From: Alex Jones Date: Tue, 5 Dec 2023 20:10:01 +0000 Subject: [PATCH 4/4] rebase (#4776) * encapsulating basic functionality of construct_discrete_space to make it easier to extend/specialize on grid type * adding staggered grid type * overloading transform_pdesys (and unfortunately for now symbolic_discretize, but this shouldn't be necessary in general) to get spatial staggering working. very close now, just need to tease out the details of the discretized differential operators now. * this commit should be reverted, this is debug code that I want to have available elsewhere * removing symbolic discretize specialization, think we have a better path forward specializing only low-level differential discretizer methods * implementing centered finite difference for staggered grid\! * fixes to resolve off-by-one error in staggered solve of wave equation * working staggered boundary conditions * specializing SciMLBase.discretize to put staggered discretizations into DynamicalODEProblem * tucking DynamicalODEProblem inside discretize, so removing it from example code * successfully returning DynamicalODEProblem from discretize * a bit of massaging, simplyfying, poking at generality * removing debug code * formulating guiding examples as tests * generalizing symbolic tracing operation * framing hack fixing changes * upstream pipelining just about done, need to check constructor works correctly and staggered data is available in DiscreteSpace object * DiscreteSpace now holds staggered dictionary created from staggered discretization kwarg * finished removing hacks, put in proper logic, beginning to add tests * adding staggered tests to runtests * adding additional tests...currently doesn't work * updating staggered discretize to output split ode system compatible with OTS solvers * removing duplicate/ambiguous code * correcting transparent boundary conditions, updating solution format to work with revised output from discretize * updating test to use splitode solution interface * no change, trying to trigger github actions * adding staggered wave equation tests to CI * fixing precompilation bug, vars->vars.intervals * remedying includes * adding required includes, correcting test bug * adding docs file for staggered grids * marking nonlinear laplacian tests as broken - directed by Alex for a known upstream issue * removing from heat equation example due to upstream failing code * dummy commit to trigger CI * dummy commit to trigger CI * Update Project.toml --------- Co-authored-by: chyett home Co-authored-by: Criston Hyett --- .github/workflows/CI.yml | 1 + Project.toml | 2 +- docs/src/staggered.md | 43 +++++++ docs/src/tutorials/heatss.md | 2 +- src/MethodOfLines.jl | 9 ++ src/discretization/discretize_vars.jl | 110 +++++++++++++---- src/discretization/generate_bc_eqs.jl | 52 ++++++++ .../generate_finite_difference_rules.jl | 17 +++ .../centered_difference.jl | 112 +++++++++++++++++- src/discretization/staggered_discretize.jl | 43 +++++++ src/interface/MOLFiniteDifference.jl | 5 + src/interface/grid_types.jl | 13 ++ src/interface/solution/MOLMetadata.jl | 4 + src/interface/solution/solution_utils.jl | 3 + src/interface/solution/timedep.jl | 10 +- test/components/staggered_constructors.jl | 42 +++++++ test/pde_systems/MOL_NonlinearProblem.jl | 6 +- test/pde_systems/wave_eq_staggered.jl | 108 +++++++++++++++++ test/runtests.jl | 9 ++ 19 files changed, 562 insertions(+), 29 deletions(-) create mode 100644 docs/src/staggered.md create mode 100644 src/discretization/staggered_discretize.jl create mode 100644 test/components/staggered_constructors.jl create mode 100644 test/pde_systems/wave_eq_staggered.jl diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 4345fb509..ec9cabc7d 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -34,6 +34,7 @@ jobs: - Burgers - Brusselator - Mixed_Derivatives + - Wave_Eq_Staggered - Complex version: - '1' diff --git a/Project.toml b/Project.toml index 1a631d54d..41cba5eb1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MethodOfLines" uuid = "94925ecb-adb7-4558-8ed8-f975c56a0bf4" authors = ["Alex Jones, "] -version = "0.10.2" +version = "0.10.3" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" diff --git a/docs/src/staggered.md b/docs/src/staggered.md new file mode 100644 index 000000000..ba5218a77 --- /dev/null +++ b/docs/src/staggered.md @@ -0,0 +1,43 @@ +# Staggered Grids + +For more examples of staggered grid implementations for the wave equation, see the tests `/test/pde_systems/wave_eq_staggered.jl` + +Hyperbolic PDEs are often most written as a system of first-order equations. When the original system is of second-order (e.g., wave equation), this results in a natural splitting of variables as demonstrated below. Staggered-grid numerical schemes have been developed for such equations as efficient and conservative finite difference schemes. + +To utilize this functionality, only two differences are needed to the user-interface. `grid_align=MethodOfLines.StaggeredGrid()` and defining which variable is edge aligned (vs center aligned), `edge_aligned_var=ϕ(t,x)` + +Solvers should be chosen carefully, the only officially supported solver is `SplitEuler`, but in theory any SSP solver could be used. + +Staggered grid functionality is still in its infancy. Please open issues if unexpected results occur or needed functionality is not present. + +``` +using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets +@parameters t x +@variables ρ(..) ϕ(..) +Dt = Differential(t); +Dx = Differential(x); + +a = 5.0; #wave speed +L = 8.0; +dx = 0.125; +dt = dx/a; # CFL condition +tmax = 10.0; + +initialFunction(x) = exp(-(x-L/2)^2); +eq = [Dt(ρ(t,x)) + Dx(ϕ(t,x)) ~ 0, + Dt(ϕ(t,x)) + a^2 * Dx(ρ(t,x)) ~ 0] +bcs = [ρ(0,x) ~ initialFunction(x), + ϕ(0.0,x) ~ 0.0, + ρ(t,L) ~ ρ(t,-L), + ϕ(t,-L) ~ ϕ(t,L)]; + +domains = [t in Interval(0.0, tmax), + x in Interval(-L, L)]; + +@named pdesys = PDESystem(eq, bcs, domains, [t,x], [ρ(t,x), ϕ(t,x)]); + +discretization = MOLFiniteDifference([x=>dx], t, grid_align=MethodOfLines.StaggeredGrid(), edge_aligned_var=ϕ(t,x)); +prob = discretize(pdesys, discretization); + +sol = solve(prob, SplitEuler(), dt=dt); +``` diff --git a/docs/src/tutorials/heatss.md b/docs/src/tutorials/heatss.md index ce7fbd522..8d83ada2e 100644 --- a/docs/src/tutorials/heatss.md +++ b/docs/src/tutorials/heatss.md @@ -1,7 +1,7 @@ # Steady State Heat Equation - No Time Dependence - NonlinearProblem Occasionally, it is desirable to solve an equation that has no time evolution, such as the steady state heat equation: -```@example heatss +``` using ModelingToolkit, MethodOfLines, DomainSets, NonlinearSolve @parameters x y diff --git a/src/MethodOfLines.jl b/src/MethodOfLines.jl index 5bd8e6107..314addddd 100644 --- a/src/MethodOfLines.jl +++ b/src/MethodOfLines.jl @@ -21,6 +21,11 @@ RuntimeGeneratedFunctions.init(@__MODULE__) # See here for the main `symbolic_discretize` and `generate_system` functions using PDEBase using PDEBase: unitindices, unitindex, remove, insert, sym_dot, VariableMap, depvar, x2i, d_orders, vcat!, update_varmap!, get_ops + +# staggered changes +using PDEBase: cardinalize_eqs!, make_pdesys_compatible, parse_bcs, generate_system, Interval +using PDEBase: error_analysis, add_metadata! + # To Extend import PDEBase.interface_errors import PDEBase.check_boundarymap @@ -33,6 +38,7 @@ import PDEBase.construct_differential_discretizer import PDEBase.discretize_equation! import PDEBase.generate_ic_defaults import PDEBase.generate_metadata +import PDEBase.symbolic_discretize import PDEBase.get_time import PDEBase.get_eqvar @@ -47,6 +53,8 @@ import Base.checkbounds import Base.getproperty import Base.ndims +import SciMLBase.discretize + # Interface include("interface/grid_types.jl") include("interface/scheme_types.jl") @@ -96,6 +104,7 @@ include("discretization/schemes/integral_expansion/integral_expansion.jl") include("discretization/generate_finite_difference_rules.jl") include("discretization/generate_bc_eqs.jl") include("discretization/generate_ic_defaults.jl") +include("discretization/staggered_discretize.jl") # Main include("scalar_discretization.jl") diff --git a/src/discretization/discretize_vars.jl b/src/discretization/discretize_vars.jl index 384a70135..c02163df2 100644 --- a/src/discretization/discretize_vars.jl +++ b/src/discretization/discretize_vars.jl @@ -79,6 +79,7 @@ struct DiscreteSpace{N,M,G} <: AbstractCartesianDiscreteSpace dxs Iaxies Igrid + staggeredvars end # * The move to DiscretizedVariable with a smart recursive getindex and custom dict based index type (?) will allow for sampling whole expressions at once, leading to much greater flexibility. Both Sym and Array interfaces will be implemented. Derivatives become the demarcation between different types of sampling => Derivatives are a custom subtype of DiscretizedVariable, with special subtypes for Nonlinear laplacian/spherical/ other types of derivatives with special handling. There is a pre discretized equation step that recognizes and replaces these with rules, and then the resulting equation is simply indexed into to generate the interior/BCs. @@ -90,6 +91,74 @@ function PDEBase.construct_discrete_space(vars::PDEBase.VariableMap, discretizat nspace = length(x̄) # Discretize space + axies = discretize_space(x̄, vars, discretization) + + # Define the grid on which the dependent variables will be evaluated (see #378) + # center_align is recommended for Dirichlet BCs + # edge_align is recommended for Neumann BCs (spatial discretization is conservative) + + grid = generate_grid(x̄, axies, vars.intervals, discretization) + dxs = generate_dxs(x̄, grid, vars, discretization) + + axies = Dict(axies) + grid = Dict(grid) + + # Build symbolic variables + Iaxies = [u => CartesianIndices(((axes(axies[x])[1] for x in remove(arguments(u), t))...,)) for u in depvars] + Igrid = [u => CartesianIndices(((axes(grid[x])[1] for x in remove(arguments(u), t))...,)) for u in depvars] + + depvarsdisc = discretize_dep_vars(depvars, grid, vars); + + return DiscreteSpace{nspace,length(depvars),G}(vars, Dict(depvarsdisc), axies, grid, Dict(dxs), Dict(Iaxies), Dict(Igrid), nothing) +end + +function PDEBase.construct_discrete_space(vars::PDEBase.VariableMap, discretization::MOLFiniteDifference{G}) where {G<:StaggeredGrid} + x̄ = vars.x̄ + t = vars.time + depvars = vars.ū + nspace = length(x̄) + + # Discretize space + axies = discretize_space(x̄, vars, discretization) + + # Define the grid on which the dependent variables will be evaluated (see #378) + # center_align is recommended for Dirichlet BCs + # edge_align is recommended for Neumann BCs (spatial discretization is conservative) + + grid = generate_grid(x̄, axies, vars.intervals, discretization) + dxs = generate_dxs(x̄, grid, vars, discretization) + + axies = Dict(axies) + grid = Dict(grid) + + # Build symbolic variables + Iaxies = [u => CartesianIndices(((axes(axies[x])[1] for x in remove(arguments(u), t))...,)) for u in depvars] + Igrid = [u => CartesianIndices(((axes(grid[x])[1] for x in remove(arguments(u), t))...,)) for u in depvars] + + depvarsdisc = discretize_dep_vars(depvars, grid, vars); + + # determine which variables are grid/stagger aligned + edge_aligned_var = operation(unwrap(discretization.kwargs[:edge_aligned_var])); + center_aligned_var = operation(unwrap(depvars[findfirst(u->operation(unwrap(u))!==edge_aligned_var, depvars)])); + staggered_dict = Dict(edge_aligned_var=>EdgeAlignedVar, center_aligned_var=>CenterAlignedVar); + + return DiscreteSpace{nspace,length(depvars),G}(vars, Dict(depvarsdisc), axies, grid, Dict(dxs), Dict(Iaxies), Dict(Igrid), staggered_dict) +end + + +function Base.getproperty(s::DiscreteSpace, p::Symbol) + if p in [:ū, :x̄, :ps, :time, :args, :x2i, :i2x] + getfield(s.vars, p) + else + getfield(s, p) + end +end + + +""" +Discretize space +""" +@inline function discretize_space(x̄, vars, discretization) axies = map(x̄) do x xdomain = vars.intervals[x] dx = prepare_dx(discretization.dxs[x], xdomain, discretization.grid_align) @@ -102,12 +171,13 @@ function PDEBase.construct_discrete_space(vars::PDEBase.VariableMap, discretizat end x => discx end - # Define the grid on which the dependent variables will be evaluated (see #378) - # center_align is recommended for Dirichlet BCs - # edge_align is recommended for Neumann BCs (spatial discretization is conservative) - - grid = generate_grid(x̄, axies, vars.intervals, discretization) + return axies; +end +""" +generate dxs +""" +@inline function generate_dxs(x̄, grid, vars, discretization) dxs = map(x̄) do x discx = Dict(grid)[x] if discx isa StepRangeLen @@ -119,14 +189,15 @@ function PDEBase.construct_discrete_space(vars::PDEBase.VariableMap, discretizat throw(ArgumentError("Supplied d$x is not a Number or AbstractVector, got $(typeof(discretization.dxs[x])) for $x")) end end + return dxs; +end - axies = Dict(axies) - grid = Dict(grid) - - # Build symbolic variables - Iaxies = [u => CartesianIndices(((axes(axies[x])[1] for x in remove(arguments(u), t))...,)) for u in depvars] - Igrid = [u => CartesianIndices(((axes(grid[x])[1] for x in remove(arguments(u), t))...,)) for u in depvars] - +""" +map dependent variables +""" +@inline function discretize_dep_vars(depvars, grid, vars) + x̄ = vars.x̄; + t = vars.time; depvarsdisc = map(depvars) do u op = SymbolicUtils.operation(u) if op isa SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}} @@ -144,16 +215,7 @@ function PDEBase.construct_discrete_space(vars::PDEBase.VariableMap, discretizat u => unwrap.(collect(first(@variables $sym(t)[uaxes...]))) end end - - return DiscreteSpace{nspace,length(depvars),G}(vars, Dict(depvarsdisc), axies, grid, Dict(dxs), Dict(Iaxies), Dict(Igrid)) -end - -function Base.getproperty(s::DiscreteSpace, p::Symbol) - if p in [:ū, :x̄, :ps, :time, :args, :x2i, :i2x] - getfield(s.vars, p) - else - getfield(s, p) - end + return depvarsdisc end """ @@ -231,6 +293,10 @@ map_symbolic_to_discrete(II::CartesianIndex, s::DiscreteSpace{N,M}) where {N,M} return axies end +@inline function generate_grid(x̄, axies, intervals, discretization::MOLFiniteDifference{G}) where {G<:StaggeredGrid} + return axies +end + @inline function generate_grid(x̄, axies, intervals, discretization::MOLFiniteDifference{G}) where {G<:EdgeAlignedGrid} dict = Dict(axies) return map(x̄) do x diff --git a/src/discretization/generate_bc_eqs.jl b/src/discretization/generate_bc_eqs.jl index a1b195f65..f5a2af087 100644 --- a/src/discretization/generate_bc_eqs.jl +++ b/src/discretization/generate_bc_eqs.jl @@ -116,6 +116,58 @@ function boundary_value_maps(II, s::DiscreteSpace{N,M,G}, boundary, derivweights return vcat(depvarderivbcmaps, integralbcmaps, depvarbcmaps) end +function boundary_value_maps(II, s::DiscreteSpace{N,M,G}, boundary, derivweights, indexmap) where {N,M,G<:StaggeredGrid} + u_, x_ = getvars(boundary) + ufunc(v, I, x) = s.discvars[v][I] + + depvarderivbcmaps = [] + depvarbcmaps = [] + + # * Assume that the BC is in terms of an explicit expression, not containing references to variables other than u_ at the boundary + u = depvar(u_, s) + args = ivs(u, s) + j = findfirst(isequal(x_), args) + IIold = II + # We need to construct a new index in case the value at the boundary appears in an equation one dimension lower + II = newindex(u_, II, s, indexmap) + val = filter(z -> z isa Number, arguments(u_))[1] + r = x_ => val + othervars = map(boundary.depvars) do v + substitute(v, r) + end + othervars = filter(v -> (length(arguments(v)) != 1) && any(isequal(x_), arguments(depvar(v, s))), othervars) + + depvarderivbcmaps = [(Differential(x_)^d)(u_) => central_difference(derivweights, II, s, [], (x2i(s, u, x_), x_), u, ufunc, d) for d in derivweights.orders[x_]]; + # generate_cartesian_rules(II, s, [u], derivweights, depvarbcmaps, indexmap, []); + depvarbcmaps = [v_ => s.discvars[depvar(v_, s)][II] for v_ in [u_; othervars]] + + # Only make a map if the integral will actually come out to the same number of dimensions as the boundary value + integralvs = unwrap.(filter(v -> !any(x -> safe_unwrap(x) isa Number, arguments(v)), boundary.depvars)) + + integralbcmaps = generate_whole_domain_integration_rules(IIold, s, integralvs, indexmap, nothing, x_) + + # Deal with the other relevant variables if boundary isa InterfaceBoundary + if boundary isa HigherOrderInterfaceBoundary + u__ = boundary.u2 + x__ = boundary.x2 + otheru = depvar(u__, s) + + is = [II[i] for i in setdiff(1:length(II), [j])] + j = x2i(s, otheru, x__) + + is = vcat(is[1:j-1], 1, is[j:end]) + II = CartesianIndex(is...) + + otherderivmaps = [(Differential(x__)^d)(u__) => central_difference(derivweights.map[Differential(x__)^d], II, s, [], (x2i(s, otheru, x__), x__), otheru, ufunc) for d in derivweights.orders[x__]] + otherbcmaps = [u__ => s.discvars[otheru][II]] + + depvarderivbcmaps = vcat(depvarderivbcmaps, otherderivmaps) + depvarbcmaps = vcat(depvarbcmaps, otherbcmaps) + end + + return vcat(depvarderivbcmaps, integralbcmaps, depvarbcmaps) +end + function boundary_value_maps(II, s::DiscreteSpace{N,M,G}, boundary, derivweights, indexmap) where {N,M,G<:CenterAlignedGrid} u_, x_ = getvars(boundary) ufunc(v, I, x) = s.discvars[v][I] diff --git a/src/discretization/generate_finite_difference_rules.jl b/src/discretization/generate_finite_difference_rules.jl index ac11911ec..5a527bc61 100644 --- a/src/discretization/generate_finite_difference_rules.jl +++ b/src/discretization/generate_finite_difference_rules.jl @@ -85,3 +85,20 @@ function generate_finite_difference_rules(II::CartesianIndex, s::DiscreteSpace, integration_rules = vcat(integration_rules, vec(generate_whole_domain_integration_rules(II, s, depvars, indexmap, terms))) return vcat(cb_rules, vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(mixed_deriv_rules_cartesian), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules) end + +function generate_finite_difference_rules(II::CartesianIndex, s::DiscreteSpace{W,M,G}, depvars, pde::Equation, derivweights::DifferentialDiscretizer, bmap, indexmap) where {W,M,G<:StaggeredGrid} + terms = split_terms(pde, s.x̄) + if length(II) != 0 + # Standard cartesian centered difference scheme + central_deriv_rules_cartesian = generate_cartesian_rules(II, s, depvars, derivweights, bmap, indexmap, terms) + else + central_deriv_rules_cartesian = [] + end + advection_rules = [] + nonlinlap_rules = [] + spherical_diffusion_rules = [] + integration_rules = [] + + integration_rules = vcat(integration_rules, vec(generate_whole_domain_integration_rules(II, s, depvars, indexmap, terms))) + return vcat(vec(spherical_diffusion_rules), vec(nonlinlap_rules), vec(central_deriv_rules_cartesian), vec(advection_rules), integration_rules) +end diff --git a/src/discretization/schemes/centered_difference/centered_difference.jl b/src/discretization/schemes/centered_difference/centered_difference.jl index 64dc87903..b6a60e143 100644 --- a/src/discretization/schemes/centered_difference/centered_difference.jl +++ b/src/discretization/schemes/centered_difference/centered_difference.jl @@ -32,7 +32,6 @@ function central_difference_weights_and_stencil(D::DerivativeOperator{T,N,Wind,D ndims(u, s) == 0 && return 0 # unit index in direction of the derivative I1 = unitindex(ndims(u, s), j) - # offset is important due to boundary proximity if (II[j] <= D.boundary_point_count) weights = D.low_boundary_coefs[II[j]] @@ -68,3 +67,114 @@ This is a catch all ruleset, as such it does not use @rule. Any even ordered der end )] for x in ivs(u, s)], init = []) for u in depvars], init = []) end + + +function generate_cartesian_rules(II::CartesianIndex, s::DiscreteSpace{N,M,G}, depvars, derivweights::DifferentialDiscretizer, bcmap, indexmap, terms) where {N,M,G<:StaggeredGrid} + central_ufunc(u, I, x) = s.discvars[u][I] + ufunc = central_ufunc; + xs = unique(reduce(safe_vcat, [ivs(u, s) for u in depvars], init=[])); + odd_orders = unique(filter(isodd, reduce(safe_vcat, [derivweights.orders[x] for x in xs], init=[]))); + placeholder = []; + for u in depvars + for x in xs + j = x2i(s,u,x) + jx = (j, x) + bs = filter_interfaces(bcmap[operation(u)][x]); + for d in odd_orders + ndims(u, s) == 0 && return 0 + # unit index in direction of the derivative + I1 = unitindex(ndims(u, s), j) + + # offset is important due to boundary proximity + haslower, hasupper = haslowerupper(bs, x) + boundary_point_count = derivweights.map[Differential(x)^d].boundary_point_count; + + if (II[j] <= boundary_point_count) & !haslower + if (s.staggeredvars[operation(u)] == EdgeAlignedVar)# can use centered diff + D = derivweights.windmap[1][Differential(x)^d]; + weights = derivweights.windmap[1][Differential(x)^d].stencil_coefs; + Itap = [II + (i*I1) for i in 0:1]; + else #need one-sided + D = derivweights.halfoffsetmap[1][Differential(x)^d]; + weights = D.low_boundary_coefs[II[j]] + offset = 1 - II[j] + Itap = [II + (i + offset) * I1 for i in 0:(D.boundary_stencil_length-1)] + end + elseif (II[j] > (length(s, x) - boundary_point_count)) & !hasupper + if (s.staggeredvars[operation(u)] == CenterAlignedVar) # can use centered diff + D = derivweights.windmap[1][Differential(x)^d]; + weights = derivweights.windmap[1][Differential(x)^d].stencil_coefs; + Itap = [II + (i*I1) for i in -1:0]; + else #need one-sided + D = derivweights.halfoffsetmap[1][Differential(x)^d]; + weights = D.high_boundary_coefs[length(s, x)-II[j]+1]; + offset = length(s, x) - II[j]; + Itap = [II + (i + offset) * I1 for i in (-D.boundary_stencil_length+1):1:0]; + end + else + if (s.staggeredvars[operation(u)] == CenterAlignedVar) + D = derivweights.windmap[1][Differential(x)^d]; + weights = D.stencil_coefs; + Itap = [bwrap(II + i * I1, bs, s, jx) for i in 0:1] + else + D = derivweights.windmap[1][Differential(x)^d]; + weights = D.stencil_coefs; + Itap = [bwrap(II + i * I1, bs, s, jx) for i in -1:0] + end + end + append!(placeholder, [(Differential(x)^d)(u) => sym_dot(weights, ufunc(u, Itap, x))]); + end + end + end + # Tap points of the stencil, this uses boundary_point_count as this is equal to half the stencil size, which is what we want. + return reduce(safe_vcat, placeholder, init = []) +end + +function central_difference(derivweights::DifferentialDiscretizer, II, s::DiscreteSpace{W,M,G}, bs, jx, u, ufunc, d) where {W,M,G<:StaggeredGrid} + placeholder = []; + ndims(u, s) == 0 && return 0 + j,x = jx; + # unit index in direction of the derivative + I1 = unitindex(ndims(u, s), j) + + # offset is important due to boundary proximity + haslower, hasupper = haslowerupper(bs, x) + boundary_point_count = derivweights.map[Differential(x)^d].boundary_point_count; + + if (II[j] <= boundary_point_count) & !haslower + if (s.staggeredvars[operation(u)] == EdgeAlignedVar)# can use centered diff + D = derivweights.windmap[1][Differential(x)^d]; + weights = derivweights.windmap[1][Differential(x)^d].stencil_coefs; + Itap = [II + (i*I1) for i in 0:1]; + else #need one-sided + D = derivweights.halfoffsetmap[1][Differential(x)^d]; + weights = D.low_boundary_coefs[II[j]] + offset = 1 - II[j] + Itap = [II + (i + offset) * I1 for i in 0:(D.boundary_stencil_length-1)] + end + elseif (II[j] > (length(s, x) - boundary_point_count)) & !hasupper + if (s.staggeredvars[operation(u)] == CenterAlignedVar) # can use centered diff + D = derivweights.windmap[1][Differential(x)^d]; + weights = derivweights.windmap[1][Differential(x)^d].stencil_coefs; + Itap = [II + (i*I1) for i in -1:0]; + else #need one-sided + D = derivweights.halfoffsetmap[1][Differential(x)^d]; + weights = D.high_boundary_coefs[length(s, x)-II[j]+1]; + offset = length(s, x) - II[j]; + Itap = [II + (i + offset) * I1 for i in (-D.boundary_stencil_length+1):1:0]; + end + else + if (s.staggeredvars[operation(u)] == CenterAlignedVar) + D = derivweights.windmap[1][Differential(x)^d]; + weights = D.stencil_coefs; + Itap = [bwrap(II + i * I1, bs, s, jx) for i in 0:1] + else + D = derivweights.windmap[1][Differential(x)^d]; + weights = D.stencil_coefs; + Itap = [bwrap(II + i * I1, bs, s, jx) for i in -1:0] + end + end + append!(placeholder, [sym_dot(weights, ufunc(u, Itap, x))]); + # Tap points of the stencil, this uses boundary_point_count as this is equal to half the stencil size, which is what we want. + return reduce(safe_vcat, placeholder, init = []) +end diff --git a/src/discretization/staggered_discretize.jl b/src/discretization/staggered_discretize.jl new file mode 100644 index 000000000..29d65f7b1 --- /dev/null +++ b/src/discretization/staggered_discretize.jl @@ -0,0 +1,43 @@ +function SciMLBase.discretize(pdesys::PDESystem, + discretization::MOLFiniteDifference{G}, + analytic = nothing, kwargs...) where {G<:StaggeredGrid} + sys, tspan = SciMLBase.symbolic_discretize(pdesys, discretization) + try + simpsys = structural_simplify(sys) + if tspan === nothing + add_metadata!(get_metadata(sys), sys) + return prob = NonlinearProblem(simpsys, ones(length(simpsys.states)); + discretization.kwargs..., kwargs...) + else + add_metadata!(get_metadata(simpsys), sys) + prob = ODEProblem(simpsys, Pair[], tspan; discretization.kwargs..., + kwargs...) + return symbolic_trace(prob, simpsys) + + end + catch e + error_analysis(sys, e) + end +end + +function symbolic_trace(prob, sys) + get_var_from_state(state) = operation(arguments(state)[1]); + states = get_states(sys); + u1_var = get_var_from_state(states[1]); + u2_var = get_var_from_state(states[findfirst(x->get_var_from_state(x)!=u1_var, states)]); + u1inds = findall(x->get_var_from_state(x)===u1_var, states); + u2inds = findall(x->get_var_from_state(x)===u2_var, states); + u1 = states[u1inds] + u2 = states[u2inds] + tracevec_1 = [i in u2inds ? states[i] : Num(0.0) for i in 1:length(states)] # maybe just 0.0 + tracevec_2 = [i in u1inds ? states[i] : Num(0.0) for i in 1:length(states)] + du1 = prob.f(tracevec_1, nothing, 0.0); + du2 = prob.f(tracevec_2, nothing, 0.0); + gen_du1 = eval(Symbolics.build_function(du1, states)[2]); + gen_du2 = eval(Symbolics.build_function(du2, states)[2]); + dynamical_f1(_du1,u,p,t) = gen_du1(_du1, u); + dynamical_f2(_du2,u,p,t) = gen_du2(_du2, u); + u0 = prob.u0;#[prob.u0[u1inds]; prob.u0[u2inds]]; + return SplitODEProblem(dynamical_f1, dynamical_f2, u0, prob.tspan); +#return DynamicalODEProblem{false}(dynamical_f1, dynamical_f2, u0[1:length(u1)], u0[length(u1)+1:end], prob.tspan) +end diff --git a/src/interface/MOLFiniteDifference.jl b/src/interface/MOLFiniteDifference.jl index 968e105e2..ebbe1abb1 100644 --- a/src/interface/MOLFiniteDifference.jl +++ b/src/interface/MOLFiniteDifference.jl @@ -49,6 +49,11 @@ function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, advection_sche @assert (time isa Num) | (time isa Nothing) "time must be a Num, or Nothing - got $(typeof(time)). See docs for MOLFiniteDifference." + if (grid_align == StaggeredGrid() && + !(:edge_aligned_var in keys(kwargs))) + @warn "when using StaggeredGrid(), you must set 'edge_aligned_var' keyword arg" + end + 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, useIR, callbacks, kwargs) diff --git a/src/interface/grid_types.jl b/src/interface/grid_types.jl index 84a58b16a..6489fc2ab 100644 --- a/src/interface/grid_types.jl +++ b/src/interface/grid_types.jl @@ -6,5 +6,18 @@ end struct EdgeAlignedGrid <: AbstractGrid end +struct StaggeredGrid <: AbstractGrid +end + const center_align=CenterAlignedGrid() const edge_align=EdgeAlignedGrid() +const stagger_align=StaggeredGrid() + + +abstract type AbstractVarAlign end + +struct CenterAlignedVar <: AbstractVarAlign +end + +struct EdgeAlignedVar <: AbstractVarAlign +end diff --git a/src/interface/solution/MOLMetadata.jl b/src/interface/solution/MOLMetadata.jl index 8710aecf9..c3f9953e9 100644 --- a/src/interface/solution/MOLMetadata.jl +++ b/src/interface/solution/MOLMetadata.jl @@ -45,3 +45,7 @@ end function PDEBase.generate_metadata(s::DiscreteSpace, disc::MOLFiniteDifference, pdesys::PDESystem, boundarymap, complexmap, metadata=nothing) return MOLMetadata(s, disc, pdesys, boundarymap, complexmap, metadata) end + +# function PDEBase.generate_metadata(s::DiscreteSpace, disc::MOLFiniteDifference{G,D}, pdesys::PDESystem, boundarymap, metadata=nothing) where {G<:StaggeredGrid} +# return MOLMetadata(s, disc, pdesys, boundarymap, metadata) +# end diff --git a/src/interface/solution/solution_utils.jl b/src/interface/solution/solution_utils.jl index 22628bd72..58a0d123f 100644 --- a/src/interface/solution/solution_utils.jl +++ b/src/interface/solution/solution_utils.jl @@ -25,6 +25,9 @@ function build_interpolation(umap, dvs, ivs, ivgrid, sol, pdesys, replaced_vars) end |> Dict end +# function build_interpolation(umap, dvs, ivs, ivgrid, sol{G}, pdesys, replaced_vars) where {G<:StaggeredGrid} +# end + sym_to_index(sym, syms) = findfirst(isequal(sym), syms) iscomplex(A::SciMLBase.AbstractPDESolution) = !isnothing(A.disc_data.complexmap) \ No newline at end of file diff --git a/src/interface/solution/timedep.jl b/src/interface/solution/timedep.jl index 068df285f..85a7ef987 100644 --- a/src/interface/solution/timedep.jl +++ b/src/interface/solution/timedep.jl @@ -1,4 +1,12 @@ +function generate_ivgrid(discretespace, ivs, t, metadata::MOLMetadata{G}) where {G} + return ((isequal(discretespace.time, x) ? t : discretespace.grid[x] for x in ivs)...,) +end + +function generate_ivgrid(discretespace, ivs, t, metadata::MOLMetadata{G}) where {G<:StaggeredGrid} + return #TODO ((isequal(discretespace.time, x) ? t : discretespace.grid[x] for x in ivs)...,) +end + function SciMLBase.PDETimeSeriesSolution(sol::SciMLBase.AbstractODESolution{T}, metadata::MOLMetadata) where {T} try odesys = sol.prob.f.sys @@ -6,7 +14,7 @@ function SciMLBase.PDETimeSeriesSolution(sol::SciMLBase.AbstractODESolution{T}, discretespace = metadata.discretespace ivs = [discretespace.time, discretespace.x̄...] - ivgrid = ((isequal(discretespace.time, x) ? sol.t : discretespace.grid[x] for x in ivs)...,) + ivgrid = generate_ivgrid(discretespace, ivs, sol.t, metadata) solved_states = if metadata.use_ODAE deriv_states = metadata.metadata[] diff --git a/test/components/staggered_constructors.jl b/test/components/staggered_constructors.jl new file mode 100644 index 000000000..064d7d821 --- /dev/null +++ b/test/components/staggered_constructors.jl @@ -0,0 +1,42 @@ +using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils, LinearAlgebra + +@testset "staggered constructor" begin + @parameters t x + @variables ρ(..) ϕ(..) + + Dt = Differential(t); + Dx = Differential(x); + + a = 5.0;#1.0/2.0; + L = 8.0; + dx = 1.0;#0.125; + dt = dx/a; + tmax = 1000.0; + + initialFunction(x) = exp(-(x)^2); + eq = [Dt(ρ(t,x)) + Dx(ϕ(t,x)) ~ 0, + Dt(ϕ(t,x)) + a^2 * Dx(ρ(t,x)) ~ 0] + bcs = [ρ(0,x) ~ initialFunction(x), + ϕ(0.0,x) ~ 0.0, + ρ(t,-L) ~ initialFunction(-L)*exp(-(t^2)), + ρ(t,L) ~ 0.0, + ϕ(t,-L) ~ exp((t^2)), + Dt(ϕ(t,L)) ~ -a^2*Dx(ρ(t,L))]; + + domains = [t in Interval(0.0, tmax), + x in Interval(-L, L)]; + + @named pdesys = PDESystem(eq, bcs, domains, [t,x], [ρ(t,x), ϕ(t,x)]); + + discretization = MOLFiniteDifference([x=>dx], t, grid_align=MethodOfLines.StaggeredGrid(), edge_aligned_var=ϕ(t,x)); + @test operation(Symbolics.unwrap(discretization.kwargs[:edge_aligned_var])) === operation(Symbolics.unwrap(ϕ(t,x))); + + v = MethodOfLines.VariableMap(pdesys, discretization) + disc_space = MethodOfLines.construct_discrete_space(v, discretization) + @test disc_space.staggeredvars[operation(Symbolics.unwrap(ρ(t,x)))] == MethodOfLines.CenterAlignedVar; + @test disc_space.staggeredvars[operation(Symbolics.unwrap(ϕ(t,x)))] == MethodOfLines.EdgeAlignedVar; + + orders = Dict(map(xx -> xx => [1; 2], v.x̄)) + diff_disc = MethodOfLines.construct_differential_discretizer(pdesys, disc_space, discretization, orders) + +end diff --git a/test/pde_systems/MOL_NonlinearProblem.jl b/test/pde_systems/MOL_NonlinearProblem.jl index 640d2e4aa..90c15ba47 100644 --- a/test/pde_systems/MOL_NonlinearProblem.jl +++ b/test/pde_systems/MOL_NonlinearProblem.jl @@ -23,7 +23,7 @@ using DomainSets end # Laplace's Equation, same as above but with MOL discretization -@testset "1D Laplace - constant solution" begin +@test_skip begin #@testset "1D Laplace - constant solution" begin @parameters x @variables u(..) Dxx = Differential(x)^2 @@ -47,7 +47,7 @@ end end # Laplace's Equation, linear solution -@testset "1D Laplace - linear solution" begin +@test_skip begin #@testset "1D Laplace - linear solution" begin @parameters x @variables u(..) Dxx = Differential(x)^2 @@ -73,7 +73,7 @@ end end # 2D heat -@testset "2D heat equation" begin +@test_skip begin #@testset "2D heat equation" begin @parameters x y @variables u(..) Dxx = Differential(x)^2 diff --git a/test/pde_systems/wave_eq_staggered.jl b/test/pde_systems/wave_eq_staggered.jl new file mode 100644 index 000000000..ed442ff7e --- /dev/null +++ b/test/pde_systems/wave_eq_staggered.jl @@ -0,0 +1,108 @@ +using ModelingToolkit, MethodOfLines, DomainSets, Test, Symbolics, SymbolicUtils, LinearAlgebra +using OrdinaryDiffEq + +@testset "1D wave equation, staggered grid, Mixed BC" begin + @parameters t x + @variables ρ(..) ϕ(..) + Dt = Differential(t); + Dx = Differential(x); + + a = 5.0; + L = 8.0; + dx = 0.125; + dt = (dx/a)^2; + tmax = 10.0; + + initialFunction(x) = exp(-(x)^2); + eq = [Dt(ρ(t,x)) + Dx(ϕ(t,x)) ~ 0, + Dt(ϕ(t,x)) + a^2 * Dx(ρ(t,x)) ~ 0] + bcs = [ρ(0,x) ~ initialFunction(x), + ϕ(0.0,x) ~ 0.0, + Dx(ρ(t,L)) ~ 0.0, + ϕ(t,-L) ~ 0.0];#-a^2*Dx(ρ(t,L))]; + + domains = [t in Interval(0.0, tmax), + x in Interval(-L, L)]; + + @named pdesys = PDESystem(eq, bcs, domains, [t,x], [ρ(t,x), ϕ(t,x)]); + + + discretization = MOLFiniteDifference([x=>dx], t, grid_align=MethodOfLines.StaggeredGrid(), edge_aligned_var=ϕ(t,x)); + prob = discretize(pdesys, discretization); + + sol = solve(prob, SplitEuler(), dt=dt); + + test_ind = floor(Int, (2(L-dx)/a)/(dt)) + @test maximum(sol[1:128,1] .- sol[1:128,test_ind]) < max(dx^2, dt) + @test maximum(sol[1:128,1] .- sol[1:128,2*test_ind]) < 10*max(dx^2, dt) +end + + +@testset "1D wave equation, staggered grid, Neumann BC" begin + @parameters t x + @variables ρ(..) ϕ(..) + Dt = Differential(t); + Dx = Differential(x); + + a = 5.0; + L = 8.0; + dx = 0.125; + dt = dx/a; + tmax = 10.0; + + initialFunction(x) = exp(-(x)^2); + eq = [Dt(ρ(t,x)) + Dx(ϕ(t,x)) ~ 0, + Dt(ϕ(t,x)) + a^2 * Dx(ρ(t,x)) ~ 0] + bcs = [ρ(0,x) ~ initialFunction(x), + ϕ(0.0,x) ~ 0.0, + Dt(ρ(t,L)) - (1/a)*Dt(ϕ(t,L)) ~ 0.0, + Dt(ρ(t,-L)) + (1/a)*Dt(ϕ(t,-L)) ~ 0.0] + + domains = [t in Interval(0.0, tmax), + x in Interval(-L, L)]; + + @named pdesys = PDESystem(eq, bcs, domains, [t,x], [ρ(t,x), ϕ(t,x)]); + + + discretization = MOLFiniteDifference([x=>dx], t, grid_align=MethodOfLines.StaggeredGrid(), edge_aligned_var=ϕ(t,x)); + prob = discretize(pdesys, discretization); + sol = solve(prob, SplitEuler(), dt=(dx/a)^2); + @test maximum(sol[:,end]) < 1e-3 +end + + +@testset "1D wave equation, staggered grid, Periodic BC" begin + @parameters t x + @variables ρ(..) ϕ(..) + Dt = Differential(t); + Dx = Differential(x); + + a = 5.0; + L = 8.0; + dx = 0.125; + dt = (dx/a)^2; + tmax = 10.0; + + initialFunction(x) = exp(-(x-L/2)^2); + eq = [Dt(ρ(t,x)) + Dx(ϕ(t,x)) ~ 0, + Dt(ϕ(t,x)) + a^2 * Dx(ρ(t,x)) ~ 0] + bcs = [ρ(0,x) ~ initialFunction(x), + ϕ(0.0,x) ~ 0.0, + ρ(t,L) ~ ρ(t,-L), + ϕ(t,-L) ~ ϕ(t,L)]; + + domains = [t in Interval(0.0, tmax), + x in Interval(-L, L)]; + + @named pdesys = PDESystem(eq, bcs, domains, [t,x], [ρ(t,x), ϕ(t,x)]); + + + discretization = MOLFiniteDifference([x=>dx], t, grid_align=MethodOfLines.StaggeredGrid(), edge_aligned_var=ϕ(t,x)); + prob = discretize(pdesys, discretization); + + sol = solve(prob, SplitEuler(), dt=dt); + + test_ind = round(Int, ((2*L-dx)/a)/(dt)) + @test maximum(sol[1:128,1] .- sol[1:128,test_ind]) < max(dx^2, dt) + @test maximum(sol[1:128,1] .- sol[1:128,2*test_ind]) < 10*max(dx^2, dt) +end diff --git a/test/runtests.jl b/test/runtests.jl index 368937d93..17bed13f8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -77,6 +77,9 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") @time @safetestset "ODEFunction" begin include("components/ODEFunction_test.jl") end + @time @safetestset "MOLFiniteDifference Interface: Staggered constructors" begin + include("components/staggered_constructors.jl") + end @time @safetestset "Discrete Callbacks" begin include("components/callbacks.jl") end @@ -143,4 +146,10 @@ const is_TRAVIS = haskey(ENV, "TRAVIS") include("pde_systems/MOL_1D_Linear_Convection.jl") end end + if GROUP == "All" || GROUP == "Wave_Eq_Staggered" + + @time @safetestset "MOLFiniteDifference Interface: 1D Wave Equation, Staggered" begin + include("pde_systems/wave_eq_staggered.jl") + end + end end