2D stationary diffusion problem & (?) MethodOfLines

Specific Domains Modelling & Simulations
,
1

I’d like to model a mass exchanger which can be represented as diffusion in a fluid flow with given stationary flow field and given coordinate-dependent diffusion coefficient. 2D, rectangular domain, rectangular grid, steady state, all linear. Dirichlet bc on part of borders (inlets), Neumann otherwise.

MethodOfLines.jl apeared to be what I need, however the first attempt was unsuccessful, the example from the manual don’t work. I’ve reformulated the same toy problem as SteadyStateProblem, it appears to work, and now I have several questions.

Code 
using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets, Plots

@parameters t x y
@variables u(..)
Dt = Differential(t)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

eq  = Dt(u(t, x, y)) ~ Dxx(u(t, x, y)) + Dyy(u(t, x, y))

bcs = [u(0, x, y) ~ x * y,
    u(t, 0, y) ~ x * y,
    u(t, 1, y) ~ x * y,
    u(t, x, 0) ~ x * y,
    u(t, x, 1) ~ x * y]

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

@named pdesys = PDESystem(eq, bcs, domains, [t, x, y], [u(t, x, y)])

g = 20 
dx = dy = 1/g

order = 2
discretization = MOLFiniteDifference([x => dx, y => dy], t)

t0 = @elapsed prob = discretize(pdesys,discretization)

steadystateprob = SteadyStateProblem(prob)

t1 = @elapsed sol = solve(steadystateprob, DynamicSS(Tsit5()))
t2 = @elapsed sol = solve(steadystateprob, DynamicSS(Tsit5()))

s = reshape(copy(sol.u), g-1, g-1)
heatmap(s)

The solution is returned as 1D array. It is not a problem to reshape it, but why?

The time to compute increases quite steep with grid density. Most of the time is spent in the discretize step, and solve appears to spend most time by compiling, not computing:
For the grid 100x100 the timings are

t0 = 925  # discretize
t1 = 181 # 1st solve
t2 = 0.358 # 2nd solve

The total time about 6 min for 75x75 and 18.5 min for 100x100: it looks like the time is proportional to the 4th power of grid side size (i.e. 2nd power of number of nodes).

My system is a bit more complex than the toy example above - should I expect much longer compute time?

I’d like to test the model with different parameter, e.g. different diffusion coefficients. Would it be possible to reuse discretization for some situations? My feeling, the grid like 100x100 or 50x200 should be (just) enough for realistic modelling in my case, but I have little experience, it could do with less.

One alternative would be to write out the finite differences by hand (a bit tedious). Still another possibility is some finite elements package (a lot to learn). Any more options? My purpose is to solve a specific technical problem with spending minimum time.

2

Just index the solution via the symbolic, i.e. sol[u], to get the shaped version. The reason is because if you have multiple dependent variables, this symbolic form will be able to generate each one and know the interpretation.

This is something being worked on by JuliaSimCompiler and symbolic array tooling.

Just do remake(prob, p = [myparameter => 1.0)) etc.

Right now manual ODE definition code for PDEs would give you faster compile time, though we’re working to remove that gap.

Finite element does not make that much sense if you have a square domains though.

1 Like
3

Chris, first thanks a lot for your support!

Now, following issues:

@parameters t x y
@variables u(..)
###
julia> sol[u]
ERROR: ArgumentError: invalid index: u⋆ of type Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}}

This one is actually just a minor inconvenience in my case.

I’ve first tried to run the Tutorial example, which errored - I’ve opened an issue.

Then tried to apply parameter to my code, with the same error.

My code

using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets, Plots

@parameters t x y
@parameters Px, Py
@variables u(…)
Dt = Differential(t)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

eq = Dt(u(t, x, y)) ~ Dxx(u(t, x, y)) * (1 - Px + Px * y^2) + Dyy(u(t, x, y)) * (1 - Py + Py * x^2)

bcs = [u(0, x, y) ~ 1,
u(t, 0, y) ~ 0,
u(t, 1, y) ~ 0,
u(t, x, 0) ~ x, # y = 0
u(t, x, 1) ~ 1-x, # y = 1
]

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

@named pdesys = PDESystem(eq, bcs, domains, [t, x, y], [u(t, x, y)], [Px=>0, Py=>0])

I hope it was some change in the syntax, not yet reflected in the manual, and not an actual bug.

4

sol[u(t,x,y)]

Something happened with parameters in the recent change. I need to handle that.

1 Like
5

For parameters, the Tutorial example is working with MethodOfLines@0.11.0 and erroring with MethodOfLines@0.11.1. I’ve commented in the issue.

So now I’ve downgraded to 11.0, and it works. I hope now I can start to work on the actual model :slight_smile:

julia> sol[u(t,x,y)]
ERROR: ArgumentError: u(t, x, y) is neither an observed nor an unknown variable.
6

The next problem: Interfaces (“You may want to connect regions with differing dynamics”)

Tried to run the Tutorial example. As the example contains a problem definition only, added the discretization and solving.

code 
using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets

@parameters t x1 x2
@variables c1(..)
@variables c2(..)
Dt = Differential(t)

Dx1 = Differential(x1)
Dxx1 = Dx1^2

Dx2 = Differential(x2)
Dxx2 = Dx2^2

D1(c) = 1 + c / 10
D2(c) = 1 / 10 + c / 10

eqs = [Dt(c1(t, x1)) ~ Dx1(D1(c1(t, x1)) * Dx1(c1(t, x1))),
    Dt(c2(t, x2)) ~ Dx2(D2(c2(t, x2)) * Dx2(c2(t, x2)))]

bcs = [c1(0, x1) ~ 1 + cospi(2 * x1),
    c2(0, x2) ~ 1 + cospi(2 * x2),
    Dx1(c1(t, 0)) ~ 0,
    c1(t, 0.5) ~ c2(t, 0.5), # Relevant interface boundary condition
    -D1(c1(t, 0.5)) * Dx1(c1(t, 0.5)) ~ -D2(c2(t, 0.5)) * Dx2(c2(t, 0.5)), # Higher order interface condition
    Dx2(c2(t, 1)) ~ 0]

domains = [t ∈ Interval(0.0, 0.15),
    x1 ∈ Interval(0.0, 0.5),
    x2 ∈ Interval(0.5, 1.0)]

@named pdesys = PDESystem(eqs, bcs, domains,
    [t, x1, x2], [c1(t, x1), c2(t, x2)])

# from this point my code
order = 2
discretization = MOLFiniteDifference([x1 => 0.05, x2 => 0.05], t)
;
prob = MethodOfLines.discretize(pdesys, discretization);
sol = solve(prob, NewtonRaphson(), saveat = 0.1);

Getting a warning

prob = MethodOfLines.discretize(pdesys, discretization);
 Warning: The system contains interface boundaries, which are not compatible with system transformation. The system will not be transformed. Please post an issue if you need this feature.

and then an error on calling solve.

ERROR: MethodError: anyeltypedual(::ArrayPartition{Union{}, Tuple{}}) is ambiguous.

MethodOfLines@0.11.0
Anything wrong with my code, or with the Tutorial example?

7

We need to update some tutorials. See the discussion in MethodOfLines complains about my model... · Issue #249 · SciML/MethodOfLines.jl · GitHub

8

Now updated to MethodOfLines@0.11.3 (i.e. yesterday’s release), adapted my code according to updated manual - it works now :smiling_face:

Also got the Tutorial example with an interface working (just had to change the algorithm in my code to Rosenbrock23).

Now trying to adapt that code with interface to my 2D problem. The code below describes a simplified case of a heat/mass exchanger with a counterflow, and semi-permeable wall in the middle.

Code counterflow exchanger 
using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets 

@parameters t x y1 y2
@variables c1(..)
@variables c2(..)

Dt = Differential(t)

Dx = Differential(x)
Dxx = Dx^2

Dy1 = Differential(y1)
Dyy1 = Dy1^2

Dy2 = Differential(y2)
Dyy2 = Dy2^2

li_vel(y) = 2 * (0.5 - y)

eq1  = Dt(c1(t, x, y1)) ~ Dxx(c1(t, x, y1)) + Dyy1(c1(t, x, y1)) - Dx(c1(t, x, y1))*li_vel(y1)
eq2  = Dt(c2(t, x, y2)) ~ Dxx(c2(t, x, y2)) + Dyy2(c2(t, x, y2)) - Dx(c2(t, x, y2))*li_vel(y2)

eqs = [eq1, eq2]

domains = [t ∈ Interval(0.0, 0.15),
    x ∈ Interval(0.0, 1.0),
    y1 ∈ Interval(0.0, 0.5),
    y2 ∈ Interval(0.5, 1.0)]

bcs = [
    c1(0, x, y1) ~ 0.5, # init
    c2(0, x, y2) ~ 0.5,
    c1(t, 0, y1) ~ 0, # flow in
    c2(t, 1, y2) ~ 1,
    Dx(c1(t, 1, y1)) ~ 0, # flow out
    Dx(c2(t, 0, y1)) ~ 0, 
    Dy1(c1(t, x, 0)) ~ 0, # bottom & top walls
    Dy2(c2(t, x, 1)) ~ 0,
    Dy1(c1(t, x, 0.5)) + c1(t, x, 0.5) ~ c2(t, x, 0.5), # membrane
    Dy2(c2(t, x, 0.5)) + c2(t, x, 0.5) ~ c1(t, x, 0.5), 
]

domains = [t ∈ Interval(0.0, 0.15),
    y1 ∈ Interval(0.0, 0.5),
    y2 ∈ Interval(0.5, 1.0)]

@named pdesys = PDESystem(eqs, bcs, domains,
    [t, x, y1, y2], [c1(t, x, y1), c2(t, x, y2)])

discretization = MOLFiniteDifference([x => 0.1, y1 => 0.1, y2 => 0.1], t)

prob = MethodOfLines.discretize(pdesys, discretization) 
# error on previous line
alg = Rosenbrock23() 
sol = solve(prob, alg, saveat = 0.1);

Unfortunately I get an error on MethodOfLines.discretize(pdesys, discretization):

ERROR: ArgumentError: invalid index: nothing of type Nothing
Stack trace 
ERROR: ArgumentError: invalid index: nothing of type Nothing
Stacktrace:
  [1] to_index(i::Nothing)
    @ Base ./indices.jl:300
  [2] to_index(A::Vector{Symbolics.VarDomainPairing}, i::Nothing)
    @ Base ./indices.jl:277
  [3] _to_indices1(A::Vector{Symbolics.VarDomainPairing}, inds::Tuple{Base.OneTo{Int64}}, I1::Nothing)
    @ Base ./indices.jl:359
  [4] to_indices
    @ ./indices.jl:354 [inlined]
  [5] to_indices
    @ ./indices.jl:345 [inlined]
  [6] getindex
    @ ./abstractarray.jl:1291 [inlined]
  [7] (::PDEBase.var"#54#62"{Vector{Symbolics.VarDomainPairing}})(x::SymbolicUtils.BasicSymbolic{Real})
    @ PDEBase ~/.julia/packages/PDEBase/o7LGc/src/variable_map.jl:34
  [8] iterate
    @ ./generator.jl:47 [inlined]
  [9] collect_to!(dest::Vector{Pair{…}}, itr::Base.Generator{Vector{…}, PDEBase.var"#54#62"{…}}, offs::Int64, st::Int64)
    @ Base ./array.jl:892
 [10] collect_to_with_first!(dest::Vector{…}, v1::Pair{…}, itr::Base.Generator{…}, st::Int64)
    @ Base ./array.jl:870
 [11] _collect(c::Vector{…}, itr::Base.Generator{…}, ::Base.EltypeUnknown, isz::Base.HasShape{…})
    @ Base ./array.jl:864
 [12] collect_similar(cont::Vector{Any}, itr::Base.Generator{Vector{Any}, PDEBase.var"#54#62"{Vector{…}}})
    @ Base ./array.jl:763
 [13] map(f::Function, A::Vector{Any})
    @ Base ./abstractarray.jl:3285
 [14] PDEBase.VariableMap(pdesys::PDESystem, disc::MOLFiniteDifference{…}; replaced_vars::Dict{…})
    @ PDEBase ~/.julia/packages/PDEBase/o7LGc/src/variable_map.jl:33
 [15] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
    @ PDEBase ~/.julia/packages/PDEBase/o7LGc/src/symbolic_discretize.jl:17
 [16] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…}; analytic::Nothing, kwargs::@Kwargs{})
    @ PDEBase ~/.julia/packages/PDEBase/o7LGc/src/discretization_state.jl:57
 [17] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
    @ PDEBase ~/.julia/packages/PDEBase/o7LGc/src/discretization_state.jl:54
 [18] top-level scope
    @ ~/Julia/MembH2Otransport.jl/src/PDE/m20_2D_middle-interf_4discourse.jl:52
Some type information was truncated. Use `show(err)` to see complete types.

Any ideas?
The manual says

Note that if you want to use a higher order interface condition, this may not work if you have no simple condition of the form c1(t, 0.5) ~ c2(t, 0.5)

I do not have that kind of conditions. However if I add a simple condition,

bcs with a simple condition 
bcs = [
    c1(0, x, y1) ~ 0.5, # init
    c2(0, x, y2) ~ 0.5,
    c1(t, 0, y1) ~ 0, # flow in
    c2(t, 1, y2) ~ 1,
    Dx(c1(t, 1, y1)) ~ 0, # flow out
    Dx(c2(t, 0, y1)) ~ 0, 
    Dy1(c1(t, x, 0)) ~ 0, # bottom & top walls
    Dy2(c2(t, x, 1)) ~ 0,
    Dy1(c1(t, x, 0.5)) + c1(t, x, 0.5) ~ c2(t, x, 0.5), # membrane
    # Dy2(c2(t, x, 0.5)) + c2(t, x, 0.5) ~ c1(t, x, 0.5), 
    c2(t, x, 0.5) ~ c1(t, x, 0.5), # simple condition
]

I still get the same error, so it is probably not the error cause, but the next problem.

9

Okay that’s an odd error :sweat_smile:. Open an issue on that.

10

Issue opened