Solving wave PDE with ModelingToolkit.jl and method of lines

fabiopierella · November 24, 2022, 2:36pm

Hi everyone,

I am trying to reproduce a standing wave in a 2D fluid domain. I am having troubles to go past the discretization step with the ModelingToolkit.jl - hereafter you find a detailed description of the problem and of the error.

Versions:
julia 1.8.3
ModelingToolkit 8.34.0

The problem

It’s an inviscid problem, so once the potential \varphi [m^2/s] is found, I can find the fluid velocity by taking its gradients in the x and z directions, as

u [m/s] =\frac{\partial \varphi}{\partial x}
w [m/s] =\frac{\partial \varphi}{\partial z}

\eta[m] is the position of the water free surface (meaning, the deviation of the free surface from the z=0 plane). The tilde means that the potential is evaluated at the still water level, so \tilde\varphi = \varphi(t,x,z=0).
The usefulness of having this extra variable is that it helps the definition of the time stepping algorithm. The problem can be in fact time-stepped according to the following equations:

w = \frac {\partial \tilde\varphi}{\partial z} = \frac{\partial \eta}{\partial t}

This equation means that the variation in time of the free surface is a consequence of the vertical velocity of the fluid (which makes sense physically - the free surface rate of change is due to the fluid moving upwards).

The other evolution equation for the potential at the free surface is
-g \eta = \frac {\partial \tilde\varphi}{\partial t} .

This comes from a simplified (i.e. linearized) Bernoulli equation.

The problem here is that at every time step I need to solve the Laplace equation \nabla^2 \varphi = 0 in the whole fluid domain to evaluate the vertical gradient of the potential, to find w and to be able to step the free surface elevation forward.

My domain

The domain looks like in the following sketch.
There is a simple rectangular domain, with solid impermeable walls, that therefore can be treated as a zero-gradient boundary in the normal directions.

The starting conditions are zero potential everywhere (meaning, fluid in still position), and a known free surface.

The code

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, Debugger

g = 9.81

@parameters x z t
@variables φ(..) φ̃(..) η(..)

Dt = Differential(t)
Dx = Differential(x)
Dz = Differential(z)
Dxx = Differential(x)^2
Dzz = Differential(z)^2

eqs  = [Dxx(φ(t, x, z)) + Dzz(φ(t, x, z)) ~ 0,
        Dt(φ̃(t,x)) ~ -g*η(t,x),
        Dt(η(t,x)) ~ Dz(φ(t, x, 1.))      
        ]

bcs = [
    φ(0, x, z) ~ 0,
    φ̃(0., x) ~ 0.,
    η(0., x) ~ cos(2*π*x),
    φ(t, x, 1.) ~ φ̃(t, x),
    Dx(φ(t, 0., z)) ~ 0.0,
    Dx(φ(t, 1., z)) ~ 0.0,
    Dz(φ(t, x, 0.)) ~ 0.0,
    Dx(φ̃(t, 0.)) ~ 0.0,
    Dx(φ̃(t, 1.)) ~ 0.0,
    Dx(η(t, 0.)) ~ 0.0,
    Dx(η(t, 1.)) ~ 0.0,
    ]

domains = [x ∈ Interval(0.0, 1.0),
           z ∈ Interval(0.0, 1.0),
           t ∈ Interval(0.0, 10.0)]

@named pdesys = PDESystem(eqs, bcs, domains, [t, x, z], 
    [φ(t, x, z), φ̃(t, x), η(t, x)])


dx = 0.1
dz = 0.1
order = 2

discretization = MOLFiniteDifference([x=>dx, z=>dz], t, 
    approx_order=order,
    grid_align=center_align)

println("Discretization:")
prob = discretize(pdesys, discretization)

The Error
(More details in the appendix)

I think somehow the system gets confused by the fact that \tilde\varphi and \eta do not depend on z.
It could be that this is not the right way to specify the equations - is there a way where I can specify that a variable is nothing but a pointer to a part of a larger domain?

Any help would be greatly appreciated!
Fabio

Appendix

ERROR: KeyError: key z not found
Stacktrace:
  [1] getindex
    @ .\dict.jl:498 [inlined]
  [2] (::MethodOfLines.var"#426#429"{Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}})(x::Sym{Real, Base.ImmutableDict{DataType, Any}})
    @ MethodOfLines .\none:0
  [3] iterate
    @ .\generator.jl:47 [inlined]
  [4] collect(itr::Base.Generator{Vector{Sym{Real, Base.ImmutableDict{DataType, Any}}}, MethodOfLines.var"#426#429"{Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}}})     
    @ Base .\array.jl:787
  [5] boundary_value_maps(II::CartesianIndex{1}, s::MethodOfLines.DiscreteSpace{2, 3, MethodOfLines.CenterAlignedGrid}, boundary::MethodOfLines.LowerBoundary, derivweights::MethodOfLines.DifferentialDiscretizer{Pair{Sym{Real, Base.ImmutableDict{DataType, Any}}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Float64, StaticArraysCore.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, indexmap::Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64})
    @ MethodOfLines C:\Users\fabpi\.julia\packages\MethodOfLines\kko34\src\discretization\generate_bc_eqs.jl:82
  [6] #413
    @ C:\Users\fabpi\.julia\packages\MethodOfLines\kko34\src\discretization\generate_bc_eqs.jl:34 [inlined]
  [7] (::MethodOfLines.var"#432#434"{CartesianIndex{1}})(f::MethodOfLines.var"#413#418"{MethodOfLines.LowerBoundary, MethodOfLines.DiscreteSpace{2, 3, MethodOfLines.CenterAlignedGrid}, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, MethodOfLines.DifferentialDiscretizer{Pair{Sym{Real, Base.ImmutableDict{DataType, Any}}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Float64, StaticArraysCore.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}})
    @ MethodOfLines C:\Users\fabpi\.julia\packages\MethodOfLines\kko34\src\discretization\generate_bc_eqs.jl:98
  [8] _mapreduce(f::MethodOfLines.var"#432#434"{CartesianIndex{1}}, op::typeof(vcat), #unused#::IndexLinear, A::Vector{MethodOfLines.var"#413#418"{b, MethodOfLines.DiscreteSpace{2, 
3, MethodOfLines.CenterAlignedGrid}, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, MethodOfLines.DifferentialDiscretizer{Pair{Sym{Real, Base.ImmutableDict{DataType, Any}}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Float64, StaticArraysCore.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}} where b})
    @ Base .\reduce.jl:438
  [9] _mapreduce_dim
    @ .\reducedim.jl:365 [inlined]
 [10] #mapreduce#765
    @ .\reducedim.jl:357 [inlined]
 [11] mapreduce
    @ .\reducedim.jl:357 [inlined]
 [12] (::MethodOfLines.var"#431#433"{MethodOfLines.DiscreteSpace{2, 3, MethodOfLines.CenterAlignedGrid}, Vector{MethodOfLines.var"#413#418"{b, MethodOfLines.DiscreteSpace{2, 3, MethodOfLines.CenterAlignedGrid}, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, MethodOfLines.DifferentialDiscretizer{Pair{Sym{Real, Base.ImmutableDict{DataType, Any}}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Float64, StaticArraysCore.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}} where b}, MethodOfLines.LowerBoundary, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, Equation})(II::CartesianIndex{1})
    @ MethodOfLines C:\Users\fabpi\.julia\packages\MethodOfLines\kko34\src\discretization\generate_bc_eqs.jl:98
 [13] iterate
    @ .\generator.jl:47 [inlined]
 [14] _collect(c::Vector{CartesianIndex{1}}, itr::Base.Generator{Vector{CartesianIndex{1}}, MethodOfLines.var"#431#433"{MethodOfLines.DiscreteSpace{2, 3, MethodOfLines.CenterAlignedGrid}, Vector{MethodOfLines.var"#413#418"{b, MethodOfLines.DiscreteSpace{2, 3, MethodOfLines.CenterAlignedGrid}, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, MethodOfLines.DifferentialDiscretizer{Pair{Sym{Real, Base.ImmutableDict{DataType, Any}}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Float64, StaticArraysCore.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}} where b}, MethodOfLines.LowerBoundary, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, Equation}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base .\array.jl:807
 [15] collect_similar(cont::Vector{CartesianIndex{1}}, itr::Base.Generator{Vector{CartesianIndex{1}}, MethodOfLines.var"#431#433"{MethodOfLines.DiscreteSpace{2, 3, MethodOfLines.CenterAlignedGrid}, Vector{MethodOfLines.var"#413#418"{b, MethodOfLines.DiscreteSpace{2, 3, MethodOfLines.CenterAlignedGrid}, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, MethodOfLines.DifferentialDiscretizer{Pair{Sym{Real, Base.ImmutableDict{DataType, Any}}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Float64, StaticArraysCore.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}} where b}, MethodOfLines.LowerBoundary, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, Equation}})
    @ Base .\array.jl:716
 [16] map(f::Function, A::Vector{CartesianIndex{1}})
    @ Base .\abstractarray.jl:2933
 [17] generate_bc_eqs(s::MethodOfLines.DiscreteSpace{2, 3, MethodOfLines.CenterAlignedGrid}, boundaryvalfuncs::Vector{MethodOfLines.var"#413#418"{b, MethodOfLines.DiscreteSpace{2, 3, MethodOfLines.CenterAlignedGrid}, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, MethodOfLines.DifferentialDiscretizer{Pair{Sym{Real, Base.ImmutableDict{DataType, Any}}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Float64, StaticArraysCore.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}} where b}, boundary::MethodOfLines.LowerBoundary, interiormap::MethodOfLines.InteriorMap, indexmap::Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64})        
    @ MethodOfLines C:\Users\fabpi\.julia\packages\MethodOfLines\kko34\src\discretization\generate_bc_eqs.jl:97
 [18] generate_bc_eqs!(bceqs::Vector{Any}, s::MethodOfLines.DiscreteSpace{2, 3, MethodOfLines.CenterAlignedGrid}, boundaryvalfuncs::Vector{MethodOfLines.var"#413#418"{b, MethodOfLines.DiscreteSpace{2, 3, MethodOfLines.CenterAlignedGrid}, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, MethodOfLines.DifferentialDiscretizer{Pair{Sym{Real, Base.ImmutableDict{DataType, Any}}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Float64, StaticArraysCore.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}} where b}, interiormap::MethodOfLines.InteriorMap, boundary::MethodOfLines.LowerBoundary)
    @ MethodOfLines C:\Users\fabpi\.julia\packages\MethodOfLines\kko34\src\discretization\generate_bc_eqs.jl:4
 [19] discretize_equation!(alleqs::Vector{Any}, bceqs::Vector{Any}, pde::Equation, interiormap::MethodOfLines.InteriorMap, eqvar::Term{Real, Base.ImmutableDict{DataType, Any}}, bcmap::Dict{Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing}, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Vector{Any}}}, depvars::Vector{Any}, s::MethodOfLines.DiscreteSpace{2, 3, MethodOfLines.CenterAlignedGrid}, derivweights::MethodOfLines.DifferentialDiscretizer{Pair{Sym{Real, Base.ImmutableDict{DataType, Any}}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Float64, StaticArraysCore.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, indexmap::Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Int64}, pmap::MethodOfLines.PeriodicMap{Val{false}()})
    @ MethodOfLines C:\Users\fabpi\.julia\packages\MethodOfLines\kko34\src\scalar_discretization.jl:9
 [20] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ MethodOfLines C:\Users\fabpi\.julia\packages\MethodOfLines\kko34\src\MOL_discretization.jl:119
 [21] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ MethodOfLines C:\Users\fabpi\.julia\packages\MethodOfLines\kko34\src\MOL_discretization.jl:138
 [22] top-level scope
    @ Julia\matrixTests\linearWaveSolver.jl:51

xtalax · November 24, 2022, 3:40pm

Hi, I develop MethodofLines.

I think I know what’s going on here - this is caused by the boundary value term appearing in Dt(η(t,x)) ~ Dz(φ(t, x, 1.)). Because the Dt term in this equation is of a dimension lower, the boundary handling is having difficulty discretizing this. Since you have that boundary value defined as φ̃(t, x), could you try substituting that in instead? This may not work, since the derivative is taken in z on that line.

I will use this code to debug the boundary handling, thanks for pointing this out! I’m currently working on improving this handling, this will provide a useful test case

fabiopierella · November 24, 2022, 6:22pm

Hi @xtalax

thanks a lot for your reply! And congrats, amazing work you are doing with the MethodOfLines package I have tried your suggestion, but I get the same error.

Indeed the reason for using Dt(η(t,x)) ~ Dz(φ(t, x,1.)) over Dt(η(t,x)) ~ Dz(φ̃(t, x)) was to make sure that the z-derivative was defined.

So as far as I understand from your answer, you don’t see any particular problem with defining the additional Laplacian inside the equations that are to be solved?

Let me know if you need further information on the problem, and thanks once again!

Fabio

xtalax · November 24, 2022, 6:37pm

Some boundary values in equations are supported, and some boundary derivatives, this looks like an uncaught case. I expect this is due to an error in index mapping, where we map variables in the argument signature of a variable to the indices they appear at so we can successfully reconstitute a cartesian index for one variable in its discrete form in to an index for another of potentially differing dimension. This is likely a sub 10 line fix, I’ll keep you posted.

As for the laplacian, I suppose you are referring to the equation Dxx(φ(t, x, z)) + Dzz(φ(t, x, z)) ~ 0, yes, that is completely fine, this steady state equation will be solved and its boundary value used to condition the other variables automatically when the above error is fixed.

Thanks for your congratulations by the way, I’m glad you find the package useful.

fabiopierella · November 24, 2022, 10:31pm

Yes, very useful indeed. If we manage to have this working, I’d also like to try to implement a more advanced nonlinear wave solver with the recipe described here: An efficient flexible-order model for 3D nonlinear water waves - ScienceDirect
The possibility to write down the system of equations in a high-level mathematical language, and to quickly get a discretized version is super useful in my view.

And yes, I meant the steady state equation, as you say

xtalax · November 28, 2022, 4:58pm

fabiopierella:

g = 9.81

@parameters x z t
@variables φ(..) φ̃(..) η(..)

Dt = Differential(t)
Dx = Differential(x)
Dz = Differential(z)
Dxx = Differential(x)^2
Dzz = Differential(z)^2

eqs  = [Dxx(φ(t, x, z)) + Dzz(φ(t, x, z)) ~ 0,
        Dt(φ̃(t,x)) ~ -g*η(t,x),
        Dt(η(t,x)) ~ Dz(φ(t, x, 1.))      
        ]

bcs = [
    φ(0, x, z) ~ 0,
    φ̃(0., x) ~ 0.,
    η(0., x) ~ cos(2*π*x),
    φ(t, x, 1.) ~ φ̃(t, x),
    Dx(φ(t, 0., z)) ~ 0.0,
    Dx(φ(t, 1., z)) ~ 0.0,
    Dz(φ(t, x, 0.)) ~ 0.0,
    Dx(φ̃(t, 0.)) ~ 0.0,
    Dx(φ̃(t, 1.)) ~ 0.0,
    Dx(η(t, 0.)) ~ 0.0,
    Dx(η(t, 1.)) ~ 0.0,
    ]

domains = [x ∈ Interval(0.0, 1.0),
           z ∈ Interval(0.0, 1.0),
           t ∈ Interval(0.0, 10.0)]

@named pdesys = PDESystem(eqs, bcs, domains, [t, x, z], 
    [φ(t, x, z), φ̃(t, x), η(t, x)])


dx = 0.1
dz = 0.1
order = 2

discretization = MOLFiniteDifference([x=>dx, z=>dz], t, 
    approx_order=order,
    grid_align=center_align)

println("Discretization:")
prob = discretize(pdesys, discretization)

This branch has the fix

fabiopierella · November 29, 2022, 12:08am

Hi @xtalax, thanks so much!
I have installed your branch and the discretization step runs perfectly now!

However, shortly after I have found this other hiccup.
If I run the solve (prob, Tsit(5)) I encounter this error:

julia> solve(prob, Tsit5())
ERROR: This solver is not able to use mass matrices.

After some reading, I seem to understand I need to solve a DAE system and not a plain PDE system.
The reason behing is it that as there is an implicit step to be taken in between each solve to get to write down the time derivatives \partial \varphi / \partial t and \partial\eta/\partial t, namely the solution of the Laplace equation \nabla ^2 \varphi = 0 ( Dxx(φ(t, x, z)) + Dzz(φ(t, x, z)) ~ 0). This means I need a different problem formulation and family of solvers (e.g. Sundials.jl).

Therefore, I have added this flag to the finite difference discretization:

discretization = MOLFiniteDifference([x=>dx, z=>dz], t, 
   approx_order=order,
   grid_align=center_align,
   use_ODAE=true)

The discretization step now fails with the stacktrace below. it seems to me that the DAE problem constructor expects du0map, u0map, tspan in addition to the system:

function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan,
                                    parammap = DiffEqBase.NullParameters();
                                    check_length = true, kwargs...) where {iip}

However, they are not passed by the discretize routine:

ERROR: LoadError: MethodError: no method matching (DAEProblem{true})(::ODESystem)
Closest candidates are:
  (DAEProblem{iip})(::ModelingToolkit.AbstractODESystem, ::Any, ::Any, ::Any) where iip at C:\Users\fabpi\.julia\packages\ModelingToolkit\EUBAR\src\systems\diffeqs\abstractodesystem.jl:761
  (DAEProblem{iip})(::ModelingToolkit.AbstractODESystem, ::Any, ::Any, ::Any, ::Any; check_length, kwargs...) where iip at C:\Users\fabpi\.julia\packages\ModelingToolkit\EUBAR\src\systems\diffeqs\abstractodesystem.jl:761
  (DAEProblem{iip})(::Any, ::Any, ::Any, ::Any) where iip at C:\Users\fabpi\.julia\packages\SciMLBase\MVhUm\src\problems\dae_problems.jl:98
  ...
Stacktrace:
 [1] DAEProblem(::ODESystem; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit C:\Users\fabpi\.julia\packages\ModelingToolkit\EUBAR\src\systems\diffeqs\abstractodesystem.jl:758
 [2] DAEProblem(::ODESystem)
   @ ModelingToolkit C:\Users\fabpi\.julia\packages\ModelingToolkit\EUBAR\src\systems\diffeqs\abstractodesystem.jl:757
 [3] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
   @ MethodOfLines C:\Users\fabpi\.julia\packages\MethodOfLines\YbYqN\src\MOL_discretization.jl:150
 [4] top-level scope
   @ C:\Users\fabpi\OneDoc\Julia\matrixTests\linearWaveSolver.jl:53
 [5] include(fname::String)
   @ Base.MainInclude .\client.jl:476
 [6] top-level scope
   @ REPL[25]:1
in expression starting at C:\Users\fabpi\OneDoc\Julia\matrixTests\linearWaveSolver.jl:53

Any hints here?
Thanks in advance for your patience

Fabio

xtalax · November 29, 2022, 12:12am

Hey, that flag use_ODAE doesn’t really work yet .

However, without the flag it is sufficient to just use a Mass Matrix solver like FBDF(), I was able to get your problem to solve with this - didn’t check if the result was sane though.

xtalax · November 29, 2022, 7:48pm

Could you mock up what a PDESystem specification of this system would look like in devectorized form? Then we can take a look at what we need to discretize it, and any interface changes we might need to make to support this type of boundary stencil.

This is just a particular way of handling a boundary, correct? I don’t see an interior discretization in the paper.

fabiopierella · November 29, 2022, 11:51pm

Hi @xtalax

so I have run the solution and it looks reasonable therefore I have marked the topic as solved - I hope this is ok, and that it is ok that we continue developing the topic a bit in this post.

heatmap

This is an animation of the potential \varphi(x,y) and you can see the periodic oscillations in the potential that are expected when you have a standing wave
I am not so sure about the constant zeros at the four corners, but it may be related with the way I specify the boundary conditions. I think I need to find a way to ultimately get rid of the \tilde \varphi variable and only solve for \varphi, it may help the solver.

Just for the record, the gif runs at a non-constant time step (so it may seem like it’s accelerating at some point). This is not physical, and the reason is that I haven’t been able to use something like sol = solve(prob, FBDF(), saveat=0.1) to have a constant time step in the output. I get this error, from which it seems that the interpolator is trying to act on a different object than what it is expecting - may ba a minor fix I suppose.

ERROR: MethodError: no method matching interp_summary(::Dict{Num, Interpolations.GriddedInterpolation{Float64, N, TCoefs, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}} where {N, TCoefs}})
Closest candidates are:
  interp_summary(!Matched::Union{DESolution, SciMLSolution}) at C:\Users\fabpi\.julia\packages\SciMLBase\MVhUm\src\interpolation.jl:42
  interp_summary(!Matched::SciMLBase.LinearInterpolation) at C:\Users\fabpi\.julia\packages\SciMLBase\MVhUm\src\interpolation.jl:36
  interp_summary(!Matched::OrdinaryDiffEq.OrdinaryDiffEqInterpolation{cacheType}) where cacheType at C:\Users\fabpi\.julia\packages\OrdinaryDiffEq\psStX\src\interp_func.jl:25
  ...
Stacktrace:

Thanks a lot for your help, very much appreciated!

fabiopierella · November 30, 2022, 12:03am

Hi @xtalax

A bit of background: when solving the nonlinear version of the linear problem that I have used for this code, then the free surface boundary conditions become more complex. Moreover, they need to be applied at the instantaneous true position of the free surface, and not at z=0 as in the above linear problem.

With the paper formulation, it is possible to rewrite the equations by substituting z with a derived variable \sigma \in [0,1] where

\sigma = \frac{z+h}{\eta+h} where \eta(x,t) is the instantaneous free surface elevation and h is the water depth. When this is done, the nonlinear problem can still be solved in a time-invariant grid function of the new variables (x,y,\sigma). This substitution though leads to a different formulation of the Laplace operator, which becomes a bit more complicated (as you can see in the paper). Hence my interest in the ModelingToolkit + MethodOfLines which would allow to write the operator out and use the magic of the package I will try to write down a ModelingToolkit system of equations and post it here.

As for the boundary stencil, I have to confess I have not studied very much on that part, as I haven’t coded this approach myself. However, the code from the paper’s author where the the boundary handling is implemented, is available at the below repository (in FORTRAN).
I have the feeling it should be definitely possible to express it in the ModelingToolkit framework.

xtalax · November 30, 2022, 2:59pm

The corners are held at dirichlet 0, this is an artifact of the boundary handling - I am exploring using an interpolation instead but this has led to difficulties so far.

fabiopierella:

I get this error, from which it seems that the interpolator is trying to act on a different object than what it is expecting - may ba a minor fix I suppose.

ERROR: MethodError: no method matching interp_summary(::Dict{Num, Interpolations.GriddedInterpolation{Float64, N, TCoefs, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}} where {N, TCoefs}})
Closest candidates are:
  interp_summary(!Matched::Union{DESolution, SciMLSolution}) at C:\Users\fabpi\.julia\packages\SciMLBase\MVhUm\src\interpolation.jl:42
  interp_summary(!Matched::SciMLBase.LinearInterpolation) at C:\Users\fabpi\.julia\packages\SciMLBase\MVhUm\src\interpolation.jl:36
  interp_summary(!Matched::OrdinaryDiffEq.OrdinaryDiffEqInterpolation{cacheType}) where cacheType at C:\Users\fabpi\.julia\packages\OrdinaryDiffEq\psStX\src\interp_func.jl:25
  ...
Stacktrace:

Thanks a lot for your help, very much appreciated!

I thought this had been fixed, are you running MOL v0.7.0? In any case, suppressing the output with a ; at the end of the line should avoid this. If you are running 0.7.0, please post full code in an issue on github and I’ll tackle it in the next bug fixing round.

WRT the nonlinear version, we have some pretty good nonlinear laplacian handling already baked in to the code, so I expect just the boundary handling will need work - those boundary discretizations can then go on to help anyone else with similar boundary conditions . Once we have an MTK PDESystem we can take the next steps.

fabiopierella · November 30, 2022, 10:11pm

Hi @xtalax

you were right - I was using an older package from your “Interfaces” branch. I did not notice you had merged it back into the main MethodOfLines repository

Now I have another problem apparently in the interpolation / saveat command: while the heatmap shows a correct behavior when there is no “saveat” keyword in the solve, it seems to go a bit crazy when I add it. It looks like the interpolation adds some noise on top of the real solution that is still distinguishable in the background.

heatmap

List of packages:

Status `C:\Users\fabpi\.julia\environments\v1.8\Project.toml`
  [31a5f54b] Debugger v0.7.6
  [0c46a032] DifferentialEquations v7.6.0
  [5b8099bc] DomainSets v0.5.14
  [28b8d3ca] GR v0.71.1
  [42fd0dbc] IterativeSolvers v0.9.2
  [94925ecb] MethodOfLines v0.7.0
  [961ee093] ModelingToolkit v8.36.0
  [8913a72c] NonlinearSolve v1.0.1
  [09606e27] ODEInterfaceDiffEq v3.11.0
  [1dea7af3] OrdinaryDiffEq v6.33.3
  [91a5bcdd] Plots v1.36.6
  [d330b81b] PyPlot v2.11.0
  [c3572dad] Sundials v4.11.1

Code

# For each step of the MOL, we need to solve another problem
using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

g = 9.81

@parameters x z t
@variables φ(..) φ̃(..) η(..)

Dt = Differential(t)
Dx = Differential(x)
Dz = Differential(z)
Dxx = Differential(x)^2
Dzz = Differential(z)^2

eqs  = [Dxx(φ(t, x, z)) + Dzz(φ(t, x, z)) ~ 0,
        Dt(φ̃(t,x)) ~ -g*η(t,x),
        Dt(η(t,x)) ~ Dz(φ(t, x, 1.))      
        ]

bcs = [
    φ(0, x, z) ~ 0,
    φ̃(0., x) ~ 0.,
    η(0., x) ~ 0.1*cos(2*π*x),
    φ(t, x, 1.) ~ φ̃(t, x),
    Dx(φ(t, 0., z)) ~ 0.0,
    Dx(φ(t, 1., z)) ~ 0.0,
    Dz(φ(t, x, 0.)) ~ 0.0,
    Dx(φ̃(t, 0.)) ~ 0.0,
    Dx(φ̃(t, 1.)) ~ 0.0,
    Dx(η(t, 0.)) ~ 0.0,
    Dx(η(t, 1.)) ~ 0.0,
    ]

domains = [x ∈ Interval(0.0, 1.0),
           z ∈ Interval(0.0, 1.0),
           t ∈ Interval(0.0, 5.0)]

### Works - xtalax "interfaces" branch
### https://github.com/xtalax/MethodOfLines.jl/tree/interfaces
@named pdesys = PDESystem(eqs, bcs, domains, [t, x, z], 
   [φ(t, x, z), φ̃(t, x), η(t, x)])

dx = 0.1
dz = 0.1
order = 2

discretization = MOLFiniteDifference([x=>dx, z=>dz], t, 
   approx_order=order,
   grid_align=center_align)

println("Discretization:")
prob = discretize(pdesys, discretization)

sol = solve(prob, FBDF(), saveat=0.025);
eta  = sol[η(t,x)]
phi = sol[φ(t,x,z)]

# import PyPlot as plt
# plt.plot(eta)

using Plots

n = size(phi, 1)
anim = @animate for t ∈ 1:n
    heatmap(phi[t, :, :]', clims=(-0.25,0.25))
end


gif(anim, "./heatmap.gif", fps=15)

xtalax · December 1, 2022, 1:09pm

Ok this is pretty strange, Did this occur also on the older interfaces branch? I’m sorry that you’ve run in to so many difficulties.

If this is not present without saveat, there is another option to get equally spaced results by calling the sol’s interpolation rather than using the sover’s with sol(0:0.025:5, :, :; dv = φ(t,x,z)),
actually you can supply whichever axes you desire for x and z too and these will be interpolated accordingly. This interpolation is unfortunately limited to linear at the moment, pending improvements in Interpolations.jl.

In this case I expect that some dodgey interpolation is the culprit.

fabiopierella · December 2, 2022, 9:01am

It works with your suggestion @xtalax and no need to be sorry - I am so grateful for the support you are giving me here.

I have found out this note here https://diffeq.sciml.ai/stable/solvers/dae_solve/ so I guess it could actually be a know issue - sorry I did not notice it before. I think it may go away if I use saveat with a linear interpolation and not with the default method.

ChrisRackauckas · December 2, 2022, 10:09am

Which method are you using? The standard recommended methods have an interpolation overload and thus that note doesn’t apply to them.

fabiopierella · December 2, 2022, 12:45pm

Hi @ChrisRackauckas

Thanks for your reply - I am using the FBDF() solver for the ODE.
As for the interpolation method, in my comment I assumed one can specify it somehow, but I haven’t looked into how one can do it. I am currently using saveat with the default options.

Fabio

ChrisRackauckas · December 2, 2022, 12:58pm

FBDF is one of the methods with Hermite and thus it’s not recommended for saveat on algebraic variables as the warning mentions.

xtalax · December 2, 2022, 1:04pm

Can you recommend an implicit mass matrix solver without this limitation?

ChrisRackauckas · December 2, 2022, 2:11pm

Rodas5P. Set it to GMRES mode for larger sizes like this and it should get rather close to FBDF. Use an iLU or AMG preconditioner.