Nonlinear Laplacian Handling in MOLFiniteDifference

[ChrisRackauckas]

This requires a rule application so that the substitution can look for an arbitrary D(a * D(u)) and rewrite with a . Some of the tricks for this are shown in the upwinding section of the code, but this requires an additional feature in the SymbolicUtils rules for interpolation into rules.

@shashi @YingboMa

[valentinsulzer]

Can you elaborate more on this? How would you rewrite something like

D(4 * exp(- 3 * u) * D(u))

?

[ChrisRackauckas]

The rewrite rule is essentially a regex, so it would catch D(4 * exp(- 3 * u) * D(u)) as D(a * D(u)) with a = 4 * exp(- 3 * u). The rewrite would then be to use the same rules as @mjsheikh 's PR #327, which is the generalized divergence form (scheme 2 of https://math.stackexchange.com/questions/1949795/explicit-finite-difference-scheme-for-nonlinear-diffusion/1950816 ) which is a bit more stable than reducing to two terms. So that would then be:

(a(x+1/2) * (u[i+1] - u[i]) - a(x-1/2) * (u[i] - u[i-1]) / dx^2

which you can then generalize as in that PR with Fornberg. Now if it's just a(x,t,u) == a(x,t) then you're done. If it's dependent on u then you have another complication. But the proper derivation of finite difference schemes comes from differentiation of interpolating polynomials:

https://mitmath.github.io/18337/lecture14/pdes_and_convolutions

So implicitly to the scheme there us a polynomial that is fit and would thus give the evaluation of u(x+1/2). For example, at second order it is:

Evaluating that at u(x+1/2) thus gives an approximation to u that retains the 2nd order accuracy in the equation and uses precisely the same points.

Because of finite difference stability issues and errors, it might only be useful to hardcode the 2nd order version, and maybe the 4th order version, since the 6th order version isn't quite useful with standard Float64 numbers. But if there's a nice way to generate these then 🤷 let's go for it. In fact, I show exactly how you could do this symbolically, because the coefficients of the polynomial can be generated from a linear solve against a matrix with a repeated structure. So you could use Symbolics.jl to generate the matrix, use symbolic \ to derive the polynomial, then use simplify.(substitute.(polys,(x->x+1/2,)) and you'd get the polynomial expression at higher order in terms of the same points as the discretization, evaluated at the half grid point. I'm not sure if doing all of that is necessary though.

[valentinsulzer]

@emmanuellujan will take the lead on this one and I will do spherical domain (#367 )

[emmanuellujan]

Hello @ChrisRackauckas, in a 3D case, do you use a polynomial like this?
g(i, j, k) = a1 x(i)^2 + a2 x(i) + a3 y(j)^2 + a4 y(j) + a5 z(k)^2 + a6 z(k) + a7 x(i) y(j) + a8 x(i) z(k) + a9 y(j) z(k) + a10
If so, which stencil nodes should I use? The cross gives me 7, but I need 10 to fit the coefficients.

[shashi]

@Suavesito-Olimpiada If you want bigger fish to fry, this is the issue to look at!

[shashi]

@Suavesito-Olimpiada Does @numrule help here at all?