sciml / pdebase.jl Goto Github PK

Calling all lurkers, yes you 👀.

Chances are if you are reading this you are interested in solving PDEs - we want to hear your ideas for a common base either here or in an issue.

What information do you need to discretize a system of PDEs, mathematically defined with ModelingToolkit.jl, in your package of choice, returning a SciMLSystem with symbolic_discretize, or SciMLProblem with discretize, or both?

We have a system for recognizing location and type of boundary condition, and extracting independent and dependent variables from the equations. In MethodOfLines.jl and (soon) NeuralPDE.jl we are using Symbolics and SymbolicUtils as an IR, see here for information on how to use rewrite rules, and here for how a simple discretizer is implemented in MOL. We want to expand this common parsing to get enough information to make it easy to discretize symbolic systems in the package of your choice.

We are hoping to collect a case study of how to write a discretizer to include in the docs for others to learn. If we can align on a common frontend language and parser (optional in stages), and common benchmarking based on PDESystemLibrary, we can massively improve the state of solving PDEs in Julia.

But what IR/form do you need for the main equations? What about domain types? What is important to you? Any ideas please let us know.

Any questions about the current state of the SciML PDEs ecosystem also welcome.

At the moment whole equations are recognized as belonging to a particular boundary, or being a type of interface between regions.
We want a method for seeing whether they are dirichlet, robin, neumann, periodic or interface as many subpackages will rely on these primatives.

This should be robust to different forms, i.e. u(0, t) - 1 ~ 0 and 0 ~ 1-u(0, t) should be recognized as dirichlet BCs, not robin, and equivalent to u(0, t) ~ 1

Should have the answer stored in a few types of struct such as:

struct DirichletBC <: AbstractAtomicBC
    """
    The dependent variable of the BC (e.g. u(0, t))
    """
    u
    """
    The independent variable of the BC (e.g. x), where this is the one held constant
    """
    x
    """
    Where this bc is defined, unused for now but will be important for exotically shaped domains
    """
    support
    """
    The RHS
    """
    expr
end

Document the interface and key assumptions about the pdesystem

I saw this update to PDEBase was pushed (#9) , so I tried using it and it looks like the array variables are coming out a bit odd.

I used the same code as in the discourse thread (https://discourse.julialang.org/t/vector-equations-with-modelingtoolkit-and-methodoflines-updated-syntax/98254/4) also pasted below, along with outputs.

using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets, PDEBase

# Parameters, variables, and derivatives
n_comp = 2
@parameters t, x, p[1:n_comp], q[1:n_comp]
@variables u(..)[1:n_comp]
#@variables u[1:n_comp](..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
params = Symbolics.scalarize(reduce(vcat,[p .=> [1.5, 2.0], q .=> [1.2, 1.8]]))
params = reduce(vcat,[p .=> [1.5, 2.0], q .=> [1.2, 1.8]])
# 1D PDE and boundary conditions

eqs  = [Dt(u(t, x)[i]) ~ p[i] * Dxx(u(t, x)[i]) for i in 1:n_comp]

bcs = [[u(0, x)[i] ~ q[i] * cos(x),
        u(t, 0)[i] ~ sin(t),
        u(t, 1)[i] ~ exp(-t) * cos(1),
        Dx(u(t,0)[i]) ~ 0.0] for i in 1:n_comp]
bcs_collected = reduce(vcat, bcs)

# Space and time domains
domains = [t ∈ Interval(0.0, 1.0),
           x ∈ Interval(0.0, 1.0)]

# PDE system

@named pdesys = PDESystem(eqs, bcs_collected, domains, [t,x], Symbolics.scalarize(reduce(vcat,[u(..)[i] for i in 1:n_comp])), Symbolics.scalarize(params))


# Method of lines discretization
dx = 0.1
order = 2
discretization = MOLFiniteDifference([x => dx], t; approx_order = order)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

# Solve ODE problem
    sol = solve(prob, Tsit5(), saveat=0.2)
pdesys.dvs
prob.problem_type.pdesys.dvs
sol.dvs

dvs from PDEsystem:

2-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 (u(..))[1]
 (u(..))[2]

dvs from problem:

2-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 var"u_Any[1]"(..)
 var"u_Any[2]"(..)

dvs from solution:

2-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 var"u_Any[1]"(..)
 var"u_Any[2]"(..)

I've had some trouble figuring out how (if it's possible) to index into the solution array variables. Is this the intended output, or is there something else I am missing.

edited to remove unused packages from "using" line

To improve generality, bcs should be recognized to belong to some lower dimensional subdomain of arbitrary shape, whether they are lower or upper or internal or interface can be recognized by comparing the bcs domain with the pdesystem's

we need to include a case study for how to implement a basic discretizer to illustrate the process

In nonlinear continuum mechanics one has to seperate the reference domain $\Omega$ from the deformed domain $\Omega_t$. Thereby, one has to distinguish the gradients w.r.t. reference and deformed ("current") coordinates. For a short introduction see https://ferrite-fem.github.io/Ferrite.jl/stable/examples/hyperelasticity/#Introduction. In the example below the gradients are distinguished by Grad and grad where the capitalized version stands for the reference domain. The first example uses a formulation based on the reference domain (in the sense of integrating over the reference domain, total Lagrangian approach) where P corresponds to the first Piola Kirchhoff stress.

using ModelingToolkit, SciMLBase
import Symbolics: ClosedInterval

@parameters X₁ x₁ f
@variables u(..) # should be u(X₁,x₁) ~ x₁ - X₁, X₁ in Ω and x₁ in deformed domain
Ω = ClosedInterval(0,1)

# TODO correct definitions.
grad(var) = Differential(x₁)(var)
Grad(var) = Differential(X₁)(var)
div(var) = Differential(x₁)(var)
Div(var) = Differential(X₁)(var)
P(Gradu) = .. # P is nonlinear stress (nonlinear generalized fluxes)

eqs = [Div(P(Grad(u)) ~ f]
bcs = [u(0.0) ~ 0., u(1.0) ~ 0.]

@named pde = PDESystem(eqs,bcs,[Ω],[X₁],[u(X₁,x₁)])

disc = MyFrameworkDiscretization(RitzGalerkin(...), TessellationInfo(...), MyFrameworkInfo(...))
prob = discretize(pdesys, disc, kwargs_to_prob...)

Another equivalent example uses the actual configuration, i.e. an updated Lagrangian approach

σ(gradu) = .. # Cauchy stress

eqs = [div(σ(grad(u)) ~ f]
bcs = [u(0.0) ~ 0., u(1.0) ~ 0.]

@named pde = PDESystem(eqs,bcs,[Ω],[x₁],[u(X₁,x₁)]) # Now, Ω needs to hold x₁, i.e. x₁ = X₁ + u

disc = MyFrameworkDiscretization(RitzGalerkin(...), TessellationInfo(...), MyFrameworkInfo(...))
prob = discretize(pdesys, disc, kwargs_to_prob...)

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

This package has potential to be a solid foundation as being the interface between ModelingToolkit.jl and several FEM/FVM/DG/... packages, which can solve PDE-based models on unstructured domains. We also might want to use some PDE discretizers to derive a system of ODEs or DAEs, which can then be fed into DifferentialEquations.jl for time integration.

A good first start might be to figure out how the interface might look like in 1D before moving forward to higher dimensions (which also might require interaction with some geometric discretizer). As a first step let us setup a steady state heat problem via ModelingToolkit.jl in strong form. From my understanding this can be done as follows (with some shortcuts taken for now):

using ModelingToolkit, SciMLBase
import Symbolics: ClosedInterval

@parameters x₁ f
@variables u(..)
Ω = ClosedInterval(0,1)

# TODO correct definitions.
grad(var) = Differential(x₁)(var)
div(var) = Differential(x₁)(var)

eqs = [div(grad(u)) ~ f]
bcs = [u(0.0) ~ 0., u(1.0) ~ 0.]

@named pde = PDESystem(eqs,bcs,[Ω],[x₁],[u(x₁)])

disc = MyFrameworkDiscretization(RitzGalerkin(...), TessellationInfo(...), MyFrameworkInfo(...))
prob = discretize(pdesys, disc, kwargs_to_prob...)

where RitzGalerkin is a type to describe general finite element discretizations of a PDE problem. This step is likely independent of the chosen framework. The symbolic discretization should hold information required for building the system of equation on some finite element. (Let us ignore trace spaces and the fancy stuff for now.) Further TessellationInfo is information regarding the geometric discretization of the domain (assuming that $\Omega$ is not discretized at this point) and PDEFrameworkInfo carries all framework specific information to assemble the linear problem. The resulting linear problem can then be solved by some solver from LinearSolve.jl.

For finite element discretizations we start with rewriting the problem into some weak form. With the current interface this should be done by dispatching should_transform and transform_pde_system. This will also introduces new variables for the test spaces, which have to be introduced without collision with existing symbols. The resulting system for the standard formulation for the heat problem using H1 ansatz space should be equivalent to:

using ModelingToolkit, LinearAlgebra
import Symbolics: ClosedInterval

@parameters x₁ f
@variables u(..) v(..)
Ω = (x₁ ∈ ClosedInterval(0,1))
∫ = Integral(Ω)
# TODO correct definitions.
grad(var) = Differential(x₁)(var)

# Bilinear form
a(u,v) = ∫(grad(u) ⋅ grad(v))
# Linear form
b(v) = ∫(f*v)
# System
eqs_weak = [a(u(x₁),v(x₁)) ~ b(v(x₁))]
bcs = [u(0.0) ~ 0., u(1.0) ~ 0.]

@named pde = PDESystem(eqs_weak,bcs,[Ω],[x₁],[u(x₁)])

where v is the introduced test variable. Automatizing the transformations should be the smaller issue here. It should also be possible to directly feed the weak problem into the PDESystem and discretize it. Four more functions remain to be dispatched, with my understanding of what they should do as comments following them:

• PDEBase.construct_discrete_space I am not sure what exactly this one should do
• PDEBase.construct_disc_state for example, if we have a quadrilateral as our element (e.g. from the RitzGalerkin struct), linear Lagrange ansatz space and the heat problem as stated above, then this should return 4 discrete variables for $u$ ($\hat{u}_1,...\hat{u}_4$) and an additional 4 variables for $v$ ($\hat{v}_1,...,\hat{v}_4$).
• PDEBase.discretize_equation! I think this should replace $u(x_1)$ with $\sum_i \hat{u}_i \phi_i(x_1)$, where $\phi_i$ is the local ansatz funciton, and $v$ with $\phi_1,...,\phi_4$, such that we have 4 discrete equations. Further, the the unknown is taken out of the integral (coming from above) and the remaining integral is either resolved analytically or replaced with numerical quadrature.

The detailed process is written down here https://ferrite-fem.github.io/Ferrite.jl/stable/manual/fe_intro/ where I tried to keep the symbols used in this issue consistent with this reference.

cc @koehlerson