[WIP] First steps in PDE autodiscretization. #250

emmanuellujan · 2020-06-11T09:18:10Z

One-dimensional PDE auto-discretization using finite-difference (work in progress)

Current PR is associated to issue: Automated PDE solving fit for DifferentialEquations.jl DifferentialEquations.jl#469
Auto-discretization code is located in src/MOL_discretization.jl. Tests are located in test/ directory.
The following 1D CFD cases are supported, nonetheless more testing is needed.
- Typical cases of the 1D diffusion equation:
  - Dt(u(t,x)) ~ Dxx(u(t,x))
  - Dt(u(t,x)) ~ D*Dxx(u(t,x))
  - Dt(u(t,x)) ~ Dx(D(t,x))*Dx(u(t,x))+D(t,x)*Dxx(u(t,x))
- Typical cases of the 1D convection equation:
  - Dt(u(t,x)) ~ -Dx(u(t,x))
  - Dt(u(t,x)) ~ -Dx(u(t,x)) + S, S is constant
  - Dt(u(t,x)) ~ -v*Dx(u(t,x)), v is constant
  - Dt(u(t,x)) ~ -Dx(v(t,x))*u(t,x)-v(t,x)*Dx(u(t,x)), v(t,x)=1
  - Dt(u(t,x)) ~ -Dx(v(t,x))*u(t,x)-v(t,x)Dx(u(t,x)), v(t,x)=0.999 + 0.001 * t * x
    - Typical cases of 1D diffusion-convection equation:
    - [x] Dt(u(t,x)) ~ Dxx(u(t,x))-vDx(u(t,x))
    - Source/sink term
    - [x] All cases above support a constant source/sink term in rhs
    - [x] All cases above support a source/sink term which depends on (t,x) in rhs
    - Boundary conditions
    - [x] All cases above support Dirichlet boundary conditions

…rc/MOL_discretization.jl). 1D toy examples added (test/): linear convection and diffusion.

test/MOL_1D_Linear_Convection.jl

src/MOL_discretization.jl

ChrisRackauckas · 2020-06-11T19:36:46Z

Looking great. @YingboMa comments?

From ChrisRackauckas: you don't need to allocate the array here. This will use a lazy concatenation to decrease the memory cost. Co-authored-by: Christopher Rackauckas <accounts@chrisrackauckas.com>

…tion.jl) to runtests.jl. Minor changes in these tests.

ChrisRackauckas · 2020-06-13T15:38:42Z

I'll look into the unrelated test failures on master.

ChrisRackauckas · 2020-06-13T20:57:18Z

Retriggering build due to upstream fixes

YingboMa

I commented on a few stylistic things. I don't quite understand what extract_bc is doing.

src/MOL_discretization.jl

ChrisRackauckas · 2020-06-15T01:42:34Z

Rebase since master now passes tests

Co-authored-by: Yingbo Ma <mayingbo5@gmail.com>

…uellujan/DiffEqOperators.jl into wip_pde_autodiscretization

ChrisRackauckas · 2020-09-02T05:57:53Z

src/MOL_discretization.jl

+function discretize_2(input,grade,order,dx,lhs_deriv_depvars,lhs_nonderiv_depvars)
+    if input isa ModelingToolkit.Constant
+        return :($input.value)
+    elseif input isa Operation
+        if input.op isa Variable
+            if haskey(lhs_nonderiv_depvars,input.op)
+                x = lhs_nonderiv_depvars[input.op]
+                if x isa ModelingToolkit.Constant
+                    expr = :($x.value)
+                else
+                    expr = Expr(x)
+                    expr = :(x=i*$dx;eval($expr))
+                end
+            elseif grade == 1
+                # TODO: the discretization order should not be the same for
+                #       first derivatives and second derivarives
+                j = findfirst(x->x==input.op, lhs_deriv_depvars)
+                expr = :((u[$j][i]-u[$j][i-1])/$dx)
+            elseif grade == 2
+                j = findfirst(x->x==input.op, lhs_deriv_depvars)
+                expr = :((u[$j][i+1]-2.0*u[$j][i]+u[$j][i-1])/($dx*$dx))
+            else
+                j = findfirst(x->x==input.op, lhs_deriv_depvars)
+                expr = :(u[$j][i])
+            end
+            return expr
+        elseif input.op isa Differential
+            grade += 1
+            return discretize_2(input.args[1],grade,order,dx,lhs_deriv_depvars,lhs_nonderiv_depvars)
+        elseif input.op isa typeof(*)
+            expr1 = discretize_2(input.args[1],grade,order,dx,lhs_deriv_depvars,lhs_nonderiv_depvars)
+            expr2 = discretize_2(input.args[2],grade,order,dx,lhs_deriv_depvars,lhs_nonderiv_depvars)
+            return Expr(:call,:*,expr1,expr2)
+        elseif input.op isa typeof(/)
+            expr1 = discretize_2(input.args[1],grade,order,dx,lhs_deriv_depvars,lhs_nonderiv_depvars)
+            expr2 = discretize_2(input.args[2],grade,order,dx,lhs_deriv_depvars,lhs_nonderiv_depvars)
+            return Expr(:call,:/,expr1,expr2)
+        elseif input.op isa typeof(-)
+            if size(input.args,1) == 2
+                expr1 = discretize_2(input.args[1],grade,order,dx,lhs_deriv_depvars,lhs_nonderiv_depvars)
+                expr2 = discretize_2(input.args[2],grade,order,dx,lhs_deriv_depvars,lhs_nonderiv_depvars)
+                return Expr(:call,:-,expr1,expr2)
+            else #if size(input.args,1) == 1
+                expr1 = discretize_2(input.args[1],grade,order,dx,lhs_deriv_depvars,lhs_nonderiv_depvars)
+                return Expr(:call,:*,:(-1),expr1)
+            end
+        elseif input.op isa typeof(+)
+            if size(input.args,1) == 2
+                expr1 = discretize_2(input.args[1],grade,order,dx,lhs_deriv_depvars,lhs_nonderiv_depvars)
+                expr2 = discretize_2(input.args[2],grade,order,dx,lhs_deriv_depvars,lhs_nonderiv_depvars)
+                return Expr(:call,:+,expr1,expr2)
+            else #if size(input.args,1) == 1
+                expr1 = discretize_2(input.args[1],grade,order,dx,lhs_deriv_depvars,lhs_nonderiv_depvars)
+                return Expr(expr1)
+            end
+        end
+    end
+end


All of this code is quite unnecessary. It should just be a call to CenteredDifference and UpwindDifference

The aim of this algorithm is to traduce an rhs such as:

ln( v(t,x)^Dx(phi(t,x)) / (u(t,x)*(1-u(t,x))) ) where, v(t,x) = 1+2*t*x

to an expression like this:

ln( (1+2*t*i*0.025)^((u[1][i]-u[1][i-1])/0.025) / (u[1][i]*(1-u[1][i])) )

where x=i * dx=i * 0.025, and "u" stores each variable (phi has index 1).

The obtained expression is then evaluated at each time step in the function "f(du,u,p,t)". Using this strategy I avoid using matrices, which is convenient in non-linear problems, and I store every variable within "u", which is consistent with the "f" interface. How would you translate the rhs that I gave using CenteredDifference?

On the other hand, last month I have been thinking how to reimplement this prototype code in a simpler way. Here is an alternative:

Assuming lhs has expressions of this form: Dx(phi(t,x)) or v(t,x)

For each rhs (e.g. Dx(v(t,x)*phi(t,x))) associated to an lhs variable (e.g. phi)

expand derivatives (does mtk provide this feature?)

e.g.: Dx(v(t,x)*phi(t,x)) => Dx(v(t,x))*phi(t,x)+v(t,x)*Dx(phi(t,x))

replace derivatives and variables with its discretizations

e.g.: Dx(v(t,x))*phi(t,x)+v(t,x)*Dx(phi(t,x)) => ((u[1][i]-u[1][i-1])/0.025)*u[2][i]+u[1][i]*((u[2][i]-u[2][i-1])/0.025)

simplify obtained expression (does mtk provide this feature?)

solve the ODE problem

to an expression like this:

You don't want the expression though, since with that form you cannot hit the stencil compiler. This is the reason for the operator form since otherwise you will not have speed or GPU-compatibility.

How would you translate the rhs that I gave using CenteredDifference?

That's an upwind operator call and then a broadcast expression.

You don't want the expression though, since with that form you cannot hit the stencil compiler. This is the reason for the operator form since otherwise you will not have speed or GPU-compatibility.

The whole code is based on building the discretizations as expressions that can be evaluated within f(du,u,p,t), avoiding expressions implies reimplementing many things.

That's an upwind operator call and then a broadcast expression.

Could you please explicitly write the outcome you expect to get when an Operation like the one below is discretized?:

ln( v(t,x)^Dx(phi(t,x)) / (u(t,x)*(1-u(t,x))) )

Dx(phi(t,x)) would be a forward difference operator, and the rest would be a broadcast expression

Dx(phi(t,x)) would be a forward difference operator,

That I had understood.

and the rest would be a broadcast expression

Here is where I do not follow you, would you please write a single line of code with the discretized rhs, using a broadcast expression and the UpwindDifference operator?

rhs = ln( v(t,x)^Dx(phi(t,x)) / (u(t,x)*(1-u(t,x))) ) discretized_rhs = ...

This code should discretize the derivative, with UpwindDifference, but also the other variables: v and u. Furthermore, this discretized code (discretized_rhs) should be evaluated w.r.t. t at each time iteration, that is one of the reasons I used expressions to build the discretization.

ChrisRackauckas · 2020-09-02T05:58:09Z

src/MOL_discretization.jl

+function get_bcs(bcs,tdomain,domain)
+    lhs_deriv_depvars_bcs = Dict()
+    n = size(bcs,1)
+    for i = 1:n
+        var = bcs[i].lhs.op
+        if var isa Variable
+            if !haskey(lhs_deriv_depvars_bcs,var)
+                lhs_deriv_depvars_bcs[var] = Array{Expr}(undef,3)
+            end
+            if isequal(bcs[i].lhs.args[1],tdomain.lower) # u(t=0,x) 
+                lhs_deriv_depvars_bcs[var][1] = Expr(bcs[i].rhs)
+            elseif isequal(bcs[i].lhs.args[2],domain.lower) # u(t,x=x_init)
+                lhs_deriv_depvars_bcs[var][2] = :(var=$(bcs[i].rhs.value))
+            elseif isequal(bcs[i].lhs.args[2],domain.upper) # u(t,x=x_final)
+                lhs_deriv_depvars_bcs[var][3] = :(var=$(bcs[i].rhs.value))
+            end
+        end
+    end
+    return lhs_deriv_depvars_bcs
+end


This should be building a Q operator?

It is a possibility, we should take into account that the BCs should be evaluated w.r.t. t at each time step within "f(du,u,p,t)".

Yes, Q is lazy so you'd just create the type on the fly before evaluation.

Ok. I assume this is planned for an upcoming PR specializing in BCs, is that right?

That should already work

Note: the function get_bcs works, the Q operator is created below in the code:

### Get boundary conditions ################################ # TODO: generalize to N equations lhs_deriv_depvars_bcs = get_bcs(pdesys.bcs,tdomain,domain[1]) t = 0.0 interior = domain[1].lower+dx[1]:dx[1]:domain[1].upper-dx[1] u0 = Array{Float64}(undef,length(interior),length(discretization)) Q = Array{RobinBC}(undef,length(discretization)) i = 1 for var in keys(discretization) bcs = lhs_deriv_depvars_bcs[var] g = eval(:((x,t) -> $(bcs[1]))) u0[:,i] = @eval $g.($interior,$t) u_x0 = eval(bcs[2]) u_x1 = eval(bcs[3]) Q[i] = DirichletBC(u_x0,u_x1) i = i+1 end

ChrisRackauckas · 2020-09-02T05:59:14Z

It looks like tests are failing?

ChrisRackauckas · 2020-09-02T06:01:41Z

Let's have a call to get in sync.

emmanuellujan · 2020-09-02T16:14:58Z

Let's have a call to get in sync.

I download the github code, at least in my notebook, the tests associated with this code are working fine.

ChrisRackauckas · 2020-09-02T16:17:35Z

I download the github code, at least in my notebook, the tests associated with this code are working fine.

https://travis-ci.com/github/SciML/DiffEqOperators.jl/builds/180792060#L602-L613

The failure is in the new code so I think it must be related.

emmanuellujan · 2020-09-02T20:09:17Z

The failure is in the new code so I think it must be related.

This is the method of lines discretization scheme

struct MOLFiniteDifference{T} <: DiffEqBase.AbstractDiscretization
  dxs::T
  order::Int
end
MOLFiniteDifference(args...;order=2) = MOLFiniteDifference(args,order)

If you invoke MOLFiniteDifference this way:

MOLFiniteDifference(0.1) 
or 
dx = 0.1
MOLFiniteDifference(dx)

Then typeof(dxs) is a Tuple{Float64}, and dxs=(0.1,)

But if you invoke MOLFiniteDifference this way:

dx = 0.1
order = 2 
MOLFiniteDifference(dx,order)

Then typeof(dxs) is a Float64

This produced the error when parsing dxs. Any idea about how to modify MOLFiniteDifference to avoid the Tuple case?

ChrisRackauckas · 2020-09-02T20:37:30Z

Oh, that needs to be an inner constructor then:

struct MOLFiniteDifference{T} <: DiffEqBase.AbstractDiscretization
  dxs::T
  order::Int
  MOLFiniteDifference(args...;order=2) = new{typeof(args)}(args,order)
end

emmanuellujan · 2020-09-07T09:15:00Z

Oh, that needs to be an inner constructor then:

This way makes dxs=(0.1,2) and order=2. I added a small change so dxs is not a tuple:

struct MOLFiniteDifference{T} <: DiffEqBase.AbstractDiscretization
    dxs::T
    order::Int
    MOLFiniteDifference(args...;order=2) = new{typeof(args[1])}(args[1],order)
end

ChrisRackauckas · 2020-10-02T01:22:38Z

src/MOL_discretization.jl

+                elseif deriv_order == 1
+                    # TODO: approx_order and forward/backward should be 
+                    #       input parameters of each derivative
+                    approx_order = 1
+                    L = UpwindDifference(deriv_order,approx_order,dx[i],len[i]-2,-1)
+                    expr = :(-1*($L*Q[$j]*u[:,$j]))
+                elseif deriv_order == 2


The difference here isn't 1 or 2, but more generally, when it's odd you want to use upwinding and when it's even you want to use centered difference.

ChrisRackauckas · 2020-10-02T01:24:19Z

src/MOL_discretization.jl

+    # TODO: discretize the following cases
+    #
+    #   1) PDE System
+    #        1.a) Transient
+    #                There is more than one indep. variable, including  't'
+    #                E.g. du/dt = d2u/dx2 + d2u/dy2 + f(t,x,y)
+    #        1.b) Stationary
+    #                There is more than one indep. variable, 't' is not included
+    #                E.g. 0 = d2u/dx2 + d2u/dy2 + f(x,y)
+    #   2) ODE System
+    #        't' is the only independent variable
+    #        The ODESystem is packed inside a PDESystem
+    #        E.g. du/dt = f(t)
+    #
+    #   Note: regarding input format, lhs must be "du/dt" or "0".


Note that it shouldn't always be an ODEProblem. If there isn't a "preferred dimension", then it should become a NonlinearProblem. We should discuss this a bit more.

ChrisRackauckas · 2020-10-02T01:25:16Z

src/MOL_discretization.jl

+                    #       input parameters of each derivative
+                    approx_order = 1
+                    L = UpwindDifference(deriv_order,approx_order,dx[i],len[i]-2,-1)
+                    expr = :(-1*($L*Q[$j]*u[:,$j]))


yup that's what I had in mind, basically. Though the way you're writing it is going to allocate and it's not going to extend to higher dimensions. We should talk next about how to get to higher dimensions.

ChrisRackauckas · 2020-10-02T01:26:08Z

src/MOL_discretization.jl

+        bcs = lhs_deriv_depvars_bcs[var]
+
+        g = eval(:((x,t) -> $(bcs[1])))
+        u_t0[:,j] = @eval $g.($(X[2][2:len[2]-1]),$t)


This should make use of https://github.com/SciML/RuntimeGeneratedFunctions.jl

ChrisRackauckas · 2020-10-02T01:26:57Z

This looks good to merge. I made a few comments on things I saw that need to change, but they can come in follow ups and some of them are very simple. Thanks!

ChrisRackauckas · 2020-10-02T01:27:06Z

🎉

emmanuellujan · 2020-10-02T01:48:56Z

Great!

Thanks for the comments, I will take them into account for the next PR :-)

ChrisRackauckas · 2020-10-02T01:54:55Z

I opened an issue per comment to keep track of them.

emmanuellujan · 2020-10-02T01:59:20Z

I opened an issue per comment to keep track of them.

Excellent

First steps in PDE autodiscretization. discretize function updated (s…

41ff011

…rc/MOL_discretization.jl). 1D toy examples added (test/): linear convection and diffusion.

emmanuellujan changed the title ~~First steps in PDE autodiscretization.~~ [WIP] First steps in PDE autodiscretization. Jun 11, 2020

ChrisRackauckas reviewed Jun 11, 2020

View reviewed changes

test/MOL_1D_Linear_Convection.jl Show resolved Hide resolved

ChrisRackauckas reviewed Jun 11, 2020

View reviewed changes

src/MOL_discretization.jl Outdated Show resolved Hide resolved

ChrisRackauckas reviewed Jun 11, 2020

View reviewed changes

src/MOL_discretization.jl Outdated Show resolved Hide resolved

ChrisRackauckas mentioned this pull request Jun 11, 2020

General PDE solving with NNs SciML/NeuralPDE.jl#69

Merged

Update src/MOL_discretization.jl

3cd458a

From ChrisRackauckas: you don't need to allocate the array here. This will use a lazy concatenation to decrease the memory cost. Co-authored-by: Christopher Rackauckas <accounts@chrisrackauckas.com>

ChrisRackauckas requested review from ChrisRackauckas and YingboMa June 12, 2020 09:56

Adding new tests (MOL_1D_Linear_Diffusion.jl and MOL_1D_Linear_Convec…

134184e

…tion.jl) to runtests.jl. Minor changes in these tests.

ChrisRackauckas closed this Jun 13, 2020

ChrisRackauckas reopened this Jun 13, 2020

YingboMa reviewed Jun 13, 2020

View reviewed changes

emmanuellujan mentioned this pull request Jun 15, 2020

Replaced eval function with a proper substitution #254

Closed

emmanuellujan and others added 12 commits June 15, 2020 04:43

Update src/MOL_discretization.jl

6d05988

Co-authored-by: Yingbo Ma <mayingbo5@gmail.com>

Update src/MOL_discretization.jl

2a981f4

Co-authored-by: Yingbo Ma <mayingbo5@gmail.com>

Update src/MOL_discretization.jl

30bb9cd

Co-authored-by: Yingbo Ma <mayingbo5@gmail.com>

Update src/MOL_discretization.jl

35f385c

Co-authored-by: Yingbo Ma <mayingbo5@gmail.com>

Minor changes in src/MOL_discretization.jl

fea539c

Co-authored-by: Yingbo Ma <mayingbo5@gmail.com>

Minor changes

c5ae0ac

Minor changes

1ed256d

Minor changes in src/MOL_discretization.jl

25fa3aa

Co-authored-by: Yingbo Ma <mayingbo5@gmail.com>

New tests

898c4f7

Merge branch 'wip_pde_autodiscretization' of https://github.com/emman…

67e64b1

…uellujan/DiffEqOperators.jl into wip_pde_autodiscretization

Minor fix

96fc98d

Advances in coupled systems.

b29f1c2

emmanuellujan mentioned this pull request Aug 31, 2020

Automated PDE solving fit for DifferentialEquations.jl SciML/DifferentialEquations.jl#469

Closed

8 tasks

ChrisRackauckas reviewed Sep 2, 2020

View reviewed changes

Minor fix

f277d38

Minor fix.

7b01f33

emmanuellujan added 4 commits September 7, 2020 06:28

Small improvement in parser (^,abs)

637b4ac

Reimplementing using broadcast expressions

7b4f12b

Removing plots from tests

beddc6c

first steps in multidimensions

dc7b601

ChrisRackauckas reviewed Oct 2, 2020

View reviewed changes

ChrisRackauckas merged commit a08a6cd into SciML:master Oct 2, 2020

This was referenced Oct 2, 2020

MOLFiniteDifference: Upwind vs central difference based on odd vs even #281

Closed

MOLFiniteDifference: NonlinearProblem #282

Closed

emmanuellujan deleted the wip_pde_autodiscretization branch October 2, 2020 01:49

ChrisRackauckas mentioned this pull request Oct 2, 2020

Use RuntimeGeneratedFunctions for the evals in the MOLFiniteDifference #283

Closed

Uh oh!

[WIP] First steps in PDE autodiscretization. #250

[WIP] First steps in PDE autodiscretization. #250

Uh oh!

Conversation

emmanuellujan commented Jun 11, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

One-dimensional PDE auto-discretization using finite-difference (work in progress)

Uh oh!

Uh oh!

Uh oh!

Uh oh!

ChrisRackauckas commented Jun 11, 2020

Uh oh!

ChrisRackauckas commented Jun 13, 2020

Uh oh!

ChrisRackauckas commented Jun 13, 2020

Uh oh!

YingboMa left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

ChrisRackauckas commented Jun 15, 2020

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ChrisRackauckas commented Sep 2, 2020

Uh oh!

ChrisRackauckas commented Sep 2, 2020

Uh oh!

emmanuellujan commented Sep 2, 2020

Uh oh!

ChrisRackauckas commented Sep 2, 2020

Uh oh!

emmanuellujan commented Sep 2, 2020

Uh oh!

ChrisRackauckas commented Sep 2, 2020

Uh oh!

emmanuellujan commented Sep 7, 2020

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ChrisRackauckas commented Oct 2, 2020

Uh oh!

ChrisRackauckas commented Oct 2, 2020

emmanuellujan commented Jun 11, 2020 •

edited

Loading