Nonlinear Laplacian Handling in MOLFiniteDifference (issue #354)

[emmanuellujan]

Fixes part of #354

[ChrisRackauckas]

Is there a reason why this is commented?

[emmanuellujan]

I am working with some test files in my computer that isolate the portion of code that performs this particular parsing. This code is not integrated to MOL_discretization.jl yet, first I need to solve some issues with the automatic generation of rules involving Differentials. Because I wanted to upload some progress, I added the portion of code of interest at about the place where it should be integrated, but as a comment.

[emmanuellujan]

I added code to handle cases like the following:

• Dt(u(t,x)) ~ Dx(u(t,x)^(-1) * Dx(u(t,x)))
• Dt(u(t,x)) ~ Dx(u(t,x)^2 * Dx(u(t,x)))
• Dt(u(t,x)) ~ Dx(1. / (1. + u(t,x)^2) * Dx(u(t,x)))
• Dt(u(t,x)) ~ Dx(1. / (u(t,x)^2 - 1.) * Dx(u(t,x)))
• Dt(u(t,x,y)) ~ Dx( (u(t,x,y)^2 / exp(x+y)^2 + sin(x+y+4t)^2)^0.5 * Dx(u(t,x,y))) +
                         Dy( (u(t,x,y)^2 / exp(x+y)^2 + sin(x+y+4t)^2)^0.5 * Dy(u(t,x,y)))

To do this I partially used the SymbolicUtils @rules.

Also, I added new test cases for all the examples above.

The function g  (see issue "#354") needs to be improved. Currently a basic approximation is implemented. I would like to solve this in another PR. Now I would like to address convection and convection-diffusion.

[emmanuellujan]

@ChrisRackauckas if you agree, you can merge this PR.

[ChrisRackauckas]

It's WIP and has merge conflicts?

[ChrisRackauckas]

why does this need to be done manually? @YingboMa @shashi don't the rules check the variations?

[emmanuellujan]

Unfortunately, I did not find a way to fully use SymbolicUtils to solve this issue. So I solved a part of it manually. An explanation of my issue can be found in the following document.: https://docs.google.com/document/d/1qbWokDYXuO5GV3ilI8R4fWWnvkBkrpQd3H6d7l9DrC4/edit

[ChrisRackauckas]

@YingboMa can you double check this?

[emmanuellujan]

It's WIP and has merge conflicts?

I update my fork from time to time, I estimate that the conflicts were generated because of that. Unfortunately I do not see an option in github that allows me to resolve them. If there is such an option, and you could point me to it, I would be happy to resolve them.

[valentinsulzer]

I think you need write access to resolve conflicts on github, but if you merge master into your local branch you will see them and be able to resolve them

[emmanuellujan]

I think you need write access to resolve conflicts on github, but if you merge master into your local branch you will see them and be able to resolve them

Thanks for the advice. In summary this is what I did:

git clone https://github.com/emmanuellujan/DiffEqOperators.jl/
cd DiffEqOperators.jl
git remote add upstream https://github.com/SciML/DiffEqOperators.jl/
git fetch upstream
git checkout master
git merge upstream/master
git checkout fix-issue-354
git merge master
# Solve conflict in editor 
git add src/MOLFiniteDifference/MOL_discretization.jl
git push

[emmanuellujan]

@ChrisRackauckas @tinosulzer

[valentinsulzer]

I am ok merging this as something that works for now so we can move forward. But I am not sure this is the most accurate approach: currently we do D(u_i+1/2, x_i+1/2) (i.e. function of the average) but usually I have seen D_i+1/2 (i.e. average of the function). These can be quite different if D is nonlinear. Would need to dig more into the numerical analysis to be sure which is best

[emmanuellujan]

I am ok merging this as something that works for now so we can move forward. But I am not sure this is the most accurate approach: currently we do D(u_i+1/2, x_i+1/2) (i.e. function of the average) but usually I have seen D_i+1/2 (i.e. average of the function). These can be quite different if D is nonlinear. Would need to dig more into the numerical analysis to be sure which is best

Currently a basic approximation is implemented. I would like to solve this in another PR. In the first page of this draft you can find the definition of g based on a degree 2 polynomial in a 3D case. If this is the correct definition, I would need to know which stencil nodes should I use. The "cross" gives me 7, but I need 10 to fit the coefficients, I think that the crossed derivatives could be avoided to use just 7 nodes.

[ChrisRackauckas]

Currently a basic approximation is implemented. I would like to solve this in another PR. In the first page of this draft you can find the definition of g based on a degree 2 polynomial in a 3D case. If this is the correct definition, I would need to know which stencil nodes should I use. The "cross" gives me 7, but I need 10 to fit the coefficients, I think that the crossed derivatives could be avoided to use just 7 nodes.

You split it to the 3 terms Del u = u_xx + u_yy + u_zz and then you have 2 unique coefficients per, + 1 for the middle. Otherwise you have to have the corner terms, which you could discretize to include the corners.

[ChrisRackauckas]

just remove

[emmanuellujan]

See below

[ChrisRackauckas]

? just remove if there's a new definition

[valentinsulzer]

it's back to the old definition which ignores the discretization order (I think) i.e. breaks #365 but not enough to cause tests to fail

[ChrisRackauckas]

weird.

[ChrisRackauckas]

Are you looking into this?

[valentinsulzer]

I can take a look tomorrow

[valentinsulzer]

The only thing that is broken with the stencil approach is "Test 00: Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y))" in the MOL_2D_Diffusion.jl test with orders 4 and 6. The test breaks because the solution is not accurate enough, not a syntax error. The reason it works with the old syntax is that it just ignores discretization.centered_order and pretends you're always asking for order 2. @emmanuellujan if you only test for order 2 in that test then you can use the syntax you commented out, and we should open an issue for the order 4 and 6 in 2D.

[emmanuellujan]

I understand that it is better not to remove those comments for now.

[valentinsulzer]

actually I'll just open a PR to your fork

[emmanuellujan]

Oh, we posted almost at the same time, I did not see your comment when I posted mine. I will check what you mention. Thank you :)

[ChrisRackauckas]

I am ok merging this as something that works for now so we can move forward. But I am not sure this is the most accurate approach: currently we do D(u_i+1/2, x_i+1/2) (i.e. function of the average) but usually I have seen D_i+1/2 (i.e. average of the function). These can be quite different if D is nonlinear. Would need to dig more into the numerical analysis to be sure which is best

That would make the approximation O(dx) instead of O(dx^2), but that can be improved in a subsequent PR. Let's get an issue open for this.

[emmanuellujan]

Currently a basic approximation is implemented. I would like to solve this in another PR. In the first page of this draft you can find the definition of g based on a degree 2 polynomial in a 3D case. If this is the correct definition, I would need to know which stencil nodes should I use. The "cross" gives me 7, but I need 10 to fit the coefficients, I think that the crossed derivatives could be avoided to use just 7 nodes.

You split it to the 3 terms Del u = u_xx + u_yy + u_zz and then you have 2 unique coefficients per, + 1 for the middle. Otherwise you have to have the corner terms, which you could discretize to include the corners.

I know, but if you want to fit a degree 2 polynomial in 3D you need 10 coefficients, ergo 10 nodes. If you use a typical stencil (without the corners), you have 7 nodes, ergo you need to neglect some terms of the polynomial. I think the cross derivative terms could be removed to this end. Another option is to use the corners. In this case, you have 15 nodes available, but you only need 10, so you need to chose which corners to use. Probably any choice in which the nodes are sufficiently distributed is enough.

[emmanuellujan]

I am ok merging this as something that works for now so we can move forward. But I am not sure this is the most accurate approach: currently we do D(u_i+1/2, x_i+1/2) (i.e. function of the average) but usually I have seen D_i+1/2 (i.e. average of the function). These can be quite different if D is nonlinear. Would need to dig more into the numerical analysis to be sure which is best

That would make the approximation O(dx) instead of O(dx^2), but that can be improved in a subsequent PR. Let's get an issue open for this.

Okay, we need to wait until this PR is merged to post the new issue.

Note: currently issue #371 contains this "sub-issue".

[ChrisRackauckas]

I know, but if you want to fit a degree 2 polynomial in 3D you need 10 coefficients, ergo 10 nodes. If you use a typical stencil (without the corners), you have 7 nodes, ergo you need to neglect some terms of the polynomial. I think the cross derivative terms could be removed to this end. Another option is to use the corners. In this case, you have 15 nodes available, but you only need 10, so you need to chose which corners to use. Probably any choice in which the nodes are sufficiently distributed is enough.

Yes, that is what I'm saying. The typical stencil is not a 2nd order 2D polynomial: it's the cross product of 2 1D polynomials. That's the derivation we should probably use.

Though the 9 coefficient 2D polynomial is used though. See the examples here:

https://en.wikipedia.org/wiki/Discrete_Laplace_operator#Implementation_via_operator_discretization

Generating the higher precision stencils wouldn't be a bad idea, given that it's easy (and doesn't impact the cache efficiency)