# Solving wave PDE with ModelingToolkit.jl and method of lines

**URL:** <https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750>\
**Category:** Modelling & Simulations\
**Tags:** question\
**Created:** [November 24, 2022, 2:36pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750 "2022-11-24T14:36:19Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![fabiopierella](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fabiopierella/32/26642_2.png) [@fabiopierella](https://discourse.julialang.org/u/fabiopierella)\
**Post date:** [November 24, 2022, 2:36pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/1 "2022-11-24T14:36:19Z")

</div>

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.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/2/225b45a80f02156cb04fe0c0402fce7c09bd7cc0.png)

**The code**

```julia
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)

 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/c/fc479580e3bc02f811a891a99d1405faffc26abb.png)

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**

```julia
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

```

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [November 24, 2022, 3:40pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/2 "2022-11-24T15:40:14Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![fabiopierella](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fabiopierella/32/26642_2.png) [@fabiopierella](https://discourse.julialang.org/u/fabiopierella)\
**Post date:** [November 24, 2022, 6:22pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/3 "2022-11-24T18:22:10Z")

</div>

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.

![image](https://global.discourse-cdn.com/julialang/original/3X/6/d/6d38c76d8f5f517d93cf48f7014c056f01d5bdd2.png)

 ![image](https://global.discourse-cdn.com/julialang/original/3X/8/0/80fbcc3e8c892b9d7ad1e151b0b12c7eb44af732.png)

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

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [November 24, 2022, 6:37pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/4 "2022-11-24T18:37:11Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![fabiopierella](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fabiopierella/32/26642_2.png) [@fabiopierella](https://discourse.julialang.org/u/fabiopierella)\
**Post date:** [November 24, 2022, 10:31pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/5 "2022-11-24T22:31:43Z")

</div>

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](https://www.sciencedirect.com/science/article/pii/S0021999108006190)  
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 😉

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [November 28, 2022, 4:58pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/6 "2022-11-28T16:58:05Z")

</div>

> [@fabiopierella](#):
>
> ```julia
> 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)
> 
> ```

> **[GitHub - xtalax/MethodOfLines.jl at interfaces](https://github.com/xtalax/MethodOfLines.jl/tree/interfaces)**
>
> interfaces

This branch has the fix 😁

---

<div class="post-metadata">

**Author:** ![fabiopierella](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fabiopierella/32/26642_2.png) [@fabiopierella](https://discourse.julialang.org/u/fabiopierella)\
**Post date:** [November 29, 2022, 12:08am UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/7 "2022-11-29T00:08:01Z")

</div>

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
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:

```julia
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:

```julia
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:

```julia
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

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [November 29, 2022, 12:12am UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/8 "2022-11-29T00:12:24Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [November 29, 2022, 7:48pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/9 "2022-11-29T19:48:10Z")

</div>

> [@fabiopierella](#):
>
> 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](https://www.sciencedirect.com/science/article/pii/S0021999108006190)  
> 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.

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.

---

<div class="post-metadata">

**Author:** ![fabiopierella](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fabiopierella/32/26642_2.png) [@fabiopierella](https://discourse.julialang.org/u/fabiopierella)\
**Post date:** [November 29, 2022, 11:51pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/10 "2022-11-29T23:51:20Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/d/f/df54f72af5094d529a7f823e28cb3924844143e0.gif)

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.

```julia
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! 🙏

---

<div class="post-metadata">

**Author:** ![fabiopierella](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fabiopierella/32/26642_2.png) [@fabiopierella](https://discourse.julialang.org/u/fabiopierella)\
**Post date:** [November 30, 2022, 12:03am UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/11 "2022-11-30T00:03:25Z")

</div>

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.

> **[GitHub - apengsigkarup/OceanWave3D-Fortran90: OceanWave3D - a very efficient...](https://github.com/apengsigkarup/OceanWave3D-Fortran90)**
>
> OceanWave3D - a very efficient coastal engineering research tool used worldwide for simulation of nonlinear and dispersive free surface waves in varying bathymetries from very deep to shallow water...

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [November 30, 2022, 2:59pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/12 "2022-11-30T14:59:15Z")

</div>

> [@fabiopierella](#):
>
> I am not so sure about the constant zeros at the four corners

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.
> 
> ```julia
> 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.

---

<div class="post-metadata">

**Author:** ![fabiopierella](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fabiopierella/32/26642_2.png) [@fabiopierella](https://discourse.julialang.org/u/fabiopierella)\
**Post date:** [November 30, 2022, 10:11pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/13 "2022-11-30T22:11:23Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/4/7/477a904023ccd8aa666e7d15fd00bc43bc4d27b3.gif)

List of packages:

```julia
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

```julia
# 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)

```

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [December 1, 2022, 1:09pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/14 "2022-12-01T13:09:56Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![fabiopierella](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fabiopierella/32/26642_2.png) [@fabiopierella](https://discourse.julialang.org/u/fabiopierella)\
**Post date:** [December 2, 2022, 9:01am UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/15 "2022-12-02T09:01:03Z")

</div>

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/](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.

 ![dae_saveat](https://global.discourse-cdn.com/julialang/original/3X/a/0/a07c3089fee425b6a6b57377ab76611d5511486f.png)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 2, 2022, 10:09am UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/16 "2022-12-02T10:09:06Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![fabiopierella](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fabiopierella/32/26642_2.png) [@fabiopierella](https://discourse.julialang.org/u/fabiopierella)\
**Post date:** [December 2, 2022, 12:45pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/17 "2022-12-02T12:45:26Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 2, 2022, 12:58pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/18 "2022-12-02T12:58:29Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [December 2, 2022, 1:04pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/19 "2022-12-02T13:04:59Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 2, 2022, 2:11pm UTC](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750/20 "2022-12-02T14:11:08Z")

</div>

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.

[Next page](https://discourse.julialang.org/t/solving-wave-pde-with-modelingtoolkit-jl-and-method-of-lines/90750.md?page=2)