OK – I solved a “stupid” gas pipe example with manual semi-discretization + ModelingToolkit and lots of help from @baggepinnen, “MTK & vector equations…”.
Now, I have figured out the basic solution of the problem using MethodOfLines
– and wonder about how to solve it using the WENOScheme
…
Here is the problem I’m considering:
Of main interest is the pressure at the outlet of the pipe:
Here is the basic solution using MethodOfLines
.
# Packages
using ModelingToolkit
using DifferentialEquations
using Plots, LaTeXStrings
using MethodOfLines, DomainSets
# Unknowns, etc.
@parameters A M R T
@variables t x p(..) v(t)
Dt = Differential(t)
Dx = Differential(x)
# -- registering functions
@register_symbolic u_p(t)
@register_symbolic u_v(t)
# -- specifying driving functions
p_a = 1.03e5
p0_step(t) = t<10 ? 10*p_a : 12*p_a
v_const(t) = 0.2 # m/s
# Model
eqs = [Dt(p(t,x)) ~ -v*Dx(p(t,x)),
v ~ u_v(t)]
# -- domains
L = 10.0
doms = [t ∈ Interval(0.0, 150), x ∈ Interval(0, L)]
# -- boundary conditions, including Initial conditions
bcs = [ p(0,x) ~ 10*p_a,
p(t,0) ~ u_p(t)]
@named pde = PDESystem(eqs,bcs,doms,[t,x],[p(t,x), v])
# -- setting registered functions equal to driving functions
u_p(t) = p0_step(t)
u_v(t) = v_const(t)
# Discretization
Nmol = 10
order = 2 # This may be increased to improve accuracy of some schemes
# -- MOL scheme
disc = MOLFiniteDifference([x=>Nmol], t, approx_order=order)
# -- discretizing model
prob = discretize(pde,disc)
# Solving model
solmol = solve(prob, Tsit5(), reltol=1e-6)
The result makes sense:
With Nmol
in [10,300]
(dashed lines), the result is very similar to what I previously got with manual semi-discretization.
Anyways, I now want to try to use the WENOScheme
istead of the default scheme for MethodOfLines
. From the documentation, it looks like I just need to replace the above statement with the following:
disc = MOLFiniteDifference([x=>Nmol], t, approx_order=order, advection_scheme = WENOScheme(epsilon = 1e-6))
However, this gives an error message. I then try to remove the argument of WENOScheme
and just set:
...advection_scheme = WENOScheme()
Now, I don’t get the error message.
However, I’m told I cannot use the Tsit5()
solver because that solver doesn’t support a mass matrix. So I remove specifying the solver and hope that the default solver will work.
OK – at first, it seems to work, but I then get an error message about singularity.
- Question: what is the proper way to switch to the WENO scheme for the above very simple case?