MethodOfLines & default vs. WENO scheme

Specific Domains Modelling & Simulations
1

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:

\frac{\partial p}{\partial t} = -v\frac{\partial p}{\partial x}, x\in \left(0,L\right] \\ v(t) = \text{known} \\ p(t,x=0) = p_0(t=0)

Of main interest is the pressure at the outlet of the pipe:

p_L(t) = p(t,x=L).

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:

image
image931×595 36.9 KB

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?
2

That epsilon is the default value, so just using WENOScheme() should be the same. I will take a look at that error though.

We’ve had best results solving with FBDF() and other implicit solvers when using WENO.

Be warned that this scheme still has some issues when used with not periodic and few BCs per boundary, relating to the interpolation used to pad the boundary which has relatively poor stability. This also appears to cause the singular exception which often arises when trying to solve with pure dirichlet boundary conditions (This is simply a bug).

Working out a better way to do this padding is an active area of research for both us, and others in the field, I have not yet seen another in the literature that is applicable to a general PDE discretizer.

Only workaround that I know works in most cases at this stage is to supply at least 2 bcs per boundary of any kind, or to use a periodic condition.

2 Likes