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?