Generalizing MOLFiniteDifference to N-order PDEs

[akashkgarg]

This is an attempt to resolve #353. The current solution determines the order of each term in the equation and enumerates those orders to create symbolic substitution rules for centered finite differences.

[akashkgarg]

This is still a work in progress. Although all current tests pass, I ran into some issues when testing with the Kuramoto–Sivashinsky (test/MOL/MOL_1D_HigherOrder.jl). Currently this fails with this error:

InvalidSystemException: The system is unbalanced. There are 51 highest order derivative variables and 53 equations.

I believe this is due to the presence of the Neumann boundary conditions. Not sure if this is something in my setup or a bug in the way Neumann conditions are handled; do we need to add some "ghost cells" or if additional independent variables need to be removed from the centered-difference equations. I noticed that "interior" points are hard-coded to ignore the left-most and right-most points:

    interior = indices[[2:length(s)-1 for s in space]...]

[YingboMa]

Try to list all the discretized variables and equations on a piece of paper (with no simplification whatsoever), then try to make them balanced. My guess is that you subconsciously reduced the boundary cells by hand without noticing. I don't think we need "ghost cells" at all, they arise automatically from tearing.

[akashkgarg]

The error was in fact from my defining both Dirichlet and Neumann boundary conditions together on the same boundary variables. The following set seems to work now, except that I'm getting relatively high error against the analytical solution. I'm sure I'm doing something incorrect in the way I've set things up here, so a pair of eyes with those familiar with this would be greatly appreciated.

Here is test snippet that sets up the KS system:

    @parameters x, t
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dx2 = Differential(x)^2
    Dx3 = Differential(x)^3
    Dx4 = Differential(x)^4

    α = 1
    β = 4
    γ = 1
    eq = Dt(u(x,t)) ~ -u(x,t)*Dx(u(x,t)) - α*Dx2(u(x,t)) - β*Dx3(u(x,t)) - γ*Dx4(u(x,t))

    u_analytic(x,t;z = -x/2+t) = 11 + 15*tanh(z) -15*tanh(z)^2 - 15*tanh(z)^3
    du(x,t;z = -x/2+t) = 15/2*(tanh(z) + 1)*(3*tanh(z) - 1)*sech(z)^2

    bcs = [u(x,0) ~ u_analytic(x,0),
           Dx(u(-10,t)) ~ du(-10,t),
           Dx(u(10,t)) ~ du(10,t)]

    # Space and time domains
    domains = [x ∈ IntervalDomain(-10.0,10.0),
               t ∈ IntervalDomain(0.0,0.5)]
    # Discretization
    dx = 0.4; dt = 0.2

    discretization = MOLFiniteDifference([x=>dx],t;centered_order=4)
    pdesys = PDESystem(eq,bcs,domains,[x,t],[u(x,t)])
    prob = discretize(pdesys,discretization)

    sol = solve(prob,Tsit5(),saveat=0.1,dt=dt)
    xs = domains[1].domain.lower+dx:dx:domains[1].domain.upper-dx
    ts = sol.t

    u_predict = sol.u
    u_real = [[u_analytic(x, t) for x in xs] for t in ts]
    u_diff = u_real - u_predict
    plot(xs, u_diff)

The predicted u per time-step is shown here:

The difference against the analytical solution:

[ChrisRackauckas]

BC error on the right hand side, or something to do with how the coefficients are dropped off near the edge?

Instead of analytic, what happens if you solve it with the discretization written out by hand? Does that have the same properties?

[akashkgarg]

My hand-written discretization blows up pretty quickly:

I did find that setting both Dirichlet and Neumann boundary conditions is required to get a decent solution to this problem:

    bcs = [u(x,0) ~ u_analytic(x,0),
           u(-10,t) ~ u_analytic(-10,t),
           u(10,t) ~ u_analytic(10,t),
           Dx(u(-10,t)) ~ du(-10,t),
           Dx(u(10,t)) ~ du(10,t)]

To accomplish this though I had to manually change the interior discrete vars to be

    interior = indices[[3:length(s)-2 for s in space]...]

currently, interior is hard-coded to be indices[[2:length(s)-1 for s in space]...]. One simple fix would be to look at the number of boundary condition equations available and alter the interior indices based on that; although that seems precarious and prone to bugs. Another alternative perhaps is to analyze the boundary-condition equations and remove the necessary equations from the interior?

[ChrisRackauckas]

structural_simplify will fail if it's over-determined. But I think the KS might need periodic BCs in order to be well-defined? https://encyclopediaofmath.org/wiki/Kuramoto-Sivashinsky_equation At least, now that I think about it, I haven't seen other BCs on it. I never thought too deeply about that though... but maybe I chose a bad test equation.

I'd be willing to have this marked as @test_broken and revisit as something to fix when doing periodic BCs.

[ChrisRackauckas]

The new test files need to get added to runtests.jl to be ran.

[bgroenks96]

I am very excited for this PR 😀

[valentinsulzer]

I did find that setting both Dirichlet and Neumann boundary conditions is required to get a decent solution to this problem

There are exceptions but as a rule of thumb you need as many bcs as the order of the highest-order derivative (only 1 - upstream - for convection, 2 for diffusion, etc)

[ChrisRackauckas]

Yeah, and

interior is hard-coded to be indices[[2:length(s)-1 for s in space]...]. One simple fix would be to look at the number of boundary condition equations available and alter the interior indices based on that; although that seems precarious and prone to bugs. Another alternative perhaps is to analyze the boundary-condition equations and remove the necessary equations from the interior?

I think it does need to get fancy and reduce the size of the interior based on the number of boundary conditions in order to ensure that the equation isn't overdetermined.

This one got much harder than intended really fast though 😅 , so I think adding the files to runtests.jl and marking it @test_broken is good, with following up in proceeding PRs and documenting what's left from this in an issue.

[valentinsulzer]

A simpler system than KS to start with would be the 1d beam equation https://en.wikipedia.org/wiki/Euler%E2%80%93Bernoulli_beam_theory

Plus a good stress test for the boundary conditions as you can have 3 on one side, 1 on the other, and up to 3rd order BCs http://www.geom.uiuc.edu/education/calc-init/static-beam/boundary.html

[ChrisRackauckas]

Oh that's a much better idea.

[akashkgarg]

There are exceptions but as a rule of thumb you need as many bcs as the order of the highest-order derivative (only 1 - upstream - for convection, 2 for diffusion, etc)

Phew. Glad to hear I wasn't crazy for thinking that.

A simpler system than KS to start with would be the 1d beam equation

I suspect we'll run into the same issue we are running into with KS; without the boundary conditions being handled properly, the system will be "over-determined" and structural simplification will fail. It seems that without the proper handling of the boundary conditions and reducing the interior this is unlikely to work for most PDEs. One simple idea might be to simply count the boundary conditions on either side and remove those variables from the interior in order to make the system balanced; although I guess a fancier approach that actually looks at the system structure is probably required for a more robust solution.

[ChrisRackauckas]

although I guess a fancier approach that actually looks at the system structure is probably required for a more robust solution.

That won't be possible here. This kind of overdetermined doesn't lead to redundant equations but instead alternative equations. For example, Dirchlet on the derivative and the value at the end determines both boundary point values. An extra interior equation is not equal to either statement u[1] = a; u[2] = b. So the extra equation is not structurally satisfied and you cannot purely symbolically eliminate u[3] - u[2] + u[1] = c knowing those other two, since in light of the BCs this becomes an equation for determining u[3] instead of being a key determiner for u[2].

[bgroenks96]

So I checked out this branch and tried it out with a simple non-linear heat diffusion problem:

@parameters t z
@variables T(..) k(..)
Dz = Differential(z)
Dt = Differential(t)
eqs = [
    Dt(T) ~ Dz(k(t,z)*Dz(T(t,z)))
]
# Space and time domains
domains = [t ∈ IntervalDomain(0.0,24*3600.0),
           z ∈ IntervalDomain(0.0,1000.0)]
bcs = [
    T(0,z) ~ sin(π*z),
    T(t,0.0) ~ 0.0,
    T(t,1000.0) ~ 10.2,
    k(t,z) ~ 1 + 0.5*sin(z)
]
pdesys = PDESystem(eqs, bcs, domains, [t,z], [T,k])
dz = range(0.0,1000.0,length=100)
discretization = MOLFiniteDifference([z=>dz],t)
sys = discretize(pdesys, discretization)

And I get the following error from discretize:

ERROR: LoadError: MethodError: no method matching operation(::Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing})

Maybe I set something up wrong?

[akashkgarg]

@bgroenks96 sorry, this is still a work in progress. I believe the source is missing a using SymbolicUtils declaration in utils.jl. I'll do another push later today with this and other fixes. Thanks for trying it out and letting me know.

[valentinsulzer]

This isn't the branch that's fixing non-linear diffusion, that's #371

[bgroenks96]

Ah, I didn't realize they were separate PRs.

[akashkgarg]

I added some work in progress changes to address this issue. I still need to make sure to test all the corner cases and that it works well with multi-dimensional problems, problems with multiple dependent and independent variables, etc. As suggested, I also need to add a beam PDE to test higher order boundary conditions.

Currently, the underlying assumption is that the provided BC equations is the same in count for all independent/dependent vars; this is probably not a safe assumption to make but we'll have to restructure much of the interior FD equations in order to handle the more generic case I think.

[akashkgarg]

I added some additional tests for the Beam equation and found a couple of additional issues, which I'll tackle next.

1. The boundary conditions currently only handle first order derivatives. Higher order derivatives aren't handled.
2. Re-writing the beam equation as a first order ODE with displacement and velocity, results in un-equal boundary conditions on the dependent vars u and v. Currently the code assumes that all dependent vars have the same number of boundary conditions.

[ChrisRackauckas]

The boundary conditions currently only handle first order derivatives. Higher order derivatives aren't handled.

I think that can be left as a future PR. Open an issue on that?

Re-writing the beam equation as a first order ODE with displacement and velocity, results in un-equal boundary conditions on the dependent vars u and v. Currently the code assumes that all dependent vars have the same number of boundary conditions.

Oh interesting. Hmm.. I'd say we merge and do some next PRs.

[ChrisRackauckas]

instead of commenting, mark @test_broken.

[akashkgarg]

I pushed the required fixes and opened a new issue for handling higher order BCs here: #385