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.