Skip to content
This repository was archived by the owner on Jul 19, 2023. It is now read-only.

WIP: 2nd order pde #434

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 25 additions & 5 deletions src/MOLFiniteDifference/MOL_discretization.jl
Original file line number Diff line number Diff line change
@@ -140,7 +140,7 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret

bclocs = map(e->substitute.(indvars,e),edgevals) # location of the boundary conditions e.g. (t,0.0,y)
edgemaps = Dict(bclocs .=> [spacevals[e...] for e in edges])
initmaps = depvars

if t != nothing
initmaps = substitute.(depvars,[t=>tspan[1]])
end
@@ -187,13 +187,31 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
# Generate initial conditions and bc equations
for bc in pdesys.bcs
bcdepvar = first(get_depvars(bc.lhs, depvar_ops))
# Make sure the bcdepvar is one of the depvars
if any(u->isequal(operation(u),operation(bcdepvar)),depvars)
# Check whether the bc is an initial condition
is_ic = false
if t != nothing && operation(bc.lhs) isa Sym && !any(x -> isequal(x, t.val), arguments(bc.lhs))
# Look for u(t,0) case
depvar = bc.lhs # e.g. u(t,0)
bcop = identity
is_ic = true
elseif t != nothing && operation(bc.lhs) isa Differential
# Look for Dt(u(t,0)) case
depvar = arguments(bc.lhs)[1] # e.g. u(t,0)
if operation(depvar) isa Sym && !any(x -> isequal(x, t.val), arguments(depvar))
bcop = operation(bc.lhs) # e.g. Dt
is_ic = true
end
end
if is_ic
# initial condition
# Assume in the form `u(...) ~ ...` for now
i = findfirst(isequal(bc.lhs),initmaps)
# Assume in the form `u(...) ~ ...` or `Dt(u(...)) ~ ...` for now
# find the index of the depvar
i = findfirst(isequal(depvar),initmaps)
if i !== nothing
push!(u0,vec(depvarsdisc[i] .=> substitute.((bc.rhs,),gridvals)))
# bcop is identity for u(...) and Dt for Dt(u(...))
push!(u0,vec(bcop.(depvarsdisc[i]) .=> substitute.((bc.rhs,),gridvals)))
end
else
# Algebraic equations for BCs
@@ -207,7 +225,7 @@ function SciMLBase.symbolic_discretize(pdesys::ModelingToolkit.PDESystem,discret
# Replace symbol in the bc lhs with the spatial discretized term
lhs = substitute(lhs,depvarbcmaps[i])
rhs = substitute.((bc.rhs,),edgemaps[bcargs])
lhs = lhs isa Vector ? lhs : [lhs] # handle 1D
lhs = lhs isa Array ? lhs : [lhs] # handle 1D
push!(bceqs,lhs .~ rhs)
end
end
@@ -375,6 +393,8 @@ end

function SciMLBase.discretize(pdesys::ModelingToolkit.PDESystem,discretization::DiffEqOperators.MOLFiniteDifference)
sys, tspan = SciMLBase.symbolic_discretize(pdesys,discretization)
# Lower order for higher-order ODEs (e.g. wave equation)
sys = sys isa NonlinearSystem ? sys : ode_order_lowering(sys)
if tspan == nothing
return prob = NonlinearProblem(sys, ones(length(sys.states)))
else
55 changes: 55 additions & 0 deletions test/MOL/MOL_1D_Wave.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,55 @@
# 1D wave equation problems

# Packages and inclusions
using ModelingToolkit,DiffEqOperators,LinearAlgebra,Test,OrdinaryDiffEq, DomainSets
using ModelingToolkit: Differential

# Tests
# @testset "Test 00: Dtt(u(t,x)) ~ Dxx(u(t,x))" begin
# Method of Manufactured Solutions
u_exact = (x,t) -> exp.(-t) * cos.(x)

# Parameters, variables, and derivatives
@parameters t x
@variables u(..)
Dt = Differential(t)
Dtt = Differential(t)^2
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
eq = Dtt(u(t,x)) ~ -Dxx(u(t,x))
bcs = [u(0,x) ~ cos(x),
Dt(u(0,x)) ~ -cos(x),
u(t,0) ~ exp(-t),
u(t,Float64(π)) ~ -exp(-t)]

# Space and time domains
domains = [t ∈ Interval(0.0,1.0),
x ∈ Interval(0.0,Float64(π))]

# PDE system
pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

# Method of lines discretization
dx = range(0.0,Float64(π),length=30)
order = 2
discretization = MOLFiniteDifference([x=>dx],t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)
sys, tspan = SciMLBase.symbolic_discretize(pdesys,discretization)
# Solve ODE problem
sol = solve(prob,Tsit5(),saveat=0.1)
x_sol = dx[2:end-1]
t_sol = sol.t

# Test against exact solution
# for i in 1:length(sol)
i = 10
exact = u_exact(x_sol, t_sol[i])
# non-differential states only (ode_sys_lowering creates additional differential states)
u_approx = sol.u[i][1:28]
@test all(isapprox.(u_approx, exact, atol=0.01))
# end
# end

2 changes: 1 addition & 1 deletion test/MOL/MOL_NonlinearProblem.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# 1D diffusion problem
# 1D nonlinear problem

# Packages and inclusions
using ModelingToolkit,DiffEqOperators,LinearAlgebra,Test
1 change: 1 addition & 0 deletions test/runtests.jl
Original file line number Diff line number Diff line change
@@ -48,6 +48,7 @@ if GROUP == "All" || GROUP == "MOLFiniteDifference"
@time @safetestset "MOLFiniteDifference Interface: 1D HigherOrder" begin include("MOL/MOL_1D_HigherOrder.jl") end
@time @safetestset "MOLFiniteDifference Interface: 1D Partial DAE" begin include("MOL/MOL_1D_PDAE.jl") end
@time @safetestset "MOLFiniteDifference Interface: Stationary Nonlinear Problems" begin include("MOL/MOL_NonlinearProblem.jl") end
@time @safetestset "MOLFiniteDifference Interface: 1D Wave Equation" begin include("MOL/MOL_1D_Wave.jl") end
end

if GROUP == "All" || GROUP == "Misc"