Interface for PDE solver on more general domains

[termi-official]

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

[xtalax]

linking for relevance
https://github.com/JuliaPDE/SurveyofPDEPackages

[xtalax]

PDEBase.construct_discrete_space is for creating an object that represents the variables in both your pure mathematical pde representation, and their equivalents in their discretized pde representation, and has the means to translate between the two.

Bear in mind that the present state of the package is more of a "throwing things at the wall to see what sticks" idea, than something that is intended to be complete.

The real task of this issue is working out what needs to be in here to support discretization of all types - this will evolve as we come to know the full representation space.

[xtalax]

Let's think about how one would represent a boundary condition on a fully unstructured (read: arbitrary) mesh, and what it would take to recognize and discretize this...

I think that we would have to allow domains (as long as they were subdomains of the full domain) to be passed to bcs in the position of usual symbolic variables - do you agree?

[xtalax]

Getting this in here early: We are going to need traits of some form to do this properly - i.e whether a certain discretization can handle a certain type of domain, and to throw early if it is unsupported.

[termi-official]

For fully unstructured domains we need to cook boundary markers into the domain description and might want to pick up on the marked boundaries for our boundary conditions. At least this is the standard workflow (to keep sane, because writing down functions to describe the boundary can become very complicated).

Neumann boundaries are usually described with the integrals popping out of the weak formulation after applying e.g. Divergence theorem. This is a bit different from what we do in FDM.

[xtalax]

How would that look? Please write some preferred syntax.

[xtalax]

Regarding what you wrote in the first post WRT syntax, a user should supply

disc = FerricsDiscretization(RitzGalerkin(...), TessellationInfo(...), PDEFrameworkInfo(...))

prob = discretize(pdesys, disc, kwargs_to_prob...)

And perhaps there is some system to work out defaults where a user has not supplied the required information.
This keeps the syntax equivalent and general enough not to be confusing

[xtalax]

The argument signature of discretize, symbolic_discretize, and their return type is somewhat fixed; this is our only limitation.

[termi-official]

Not confident yet about how this could be approached. The first issue here would be to describe the (nonlinear) geometry itself. Maybe a more concrete example could help here. Let us assume something standard: vortex shedding in 2D, where our domain is a rectangle minus some hole. Currently we generate the grid via GMSH (see e.g. https://github.com/Ferrite-FEM/Ferrite.jl/blob/mk/nsgmsh/docs/src/literate/ns_vs_diffeq.jl#L126-L155) and just mark stuff along when we describe the geometry. It should be possible to describe the operations here on a higher level.

To give another use-case, my data comes from some imaging pipeline where we start with some voxel image, followed by some remeshing procedure to generate smoother surfaces (and coarser grids). Here the resulting surfaces usually do not come with an analytical description. Such workflows are also not that uncommon. Some meshes are also hand-optimized. This should showcase that it is still useful to allow a discretized domain as input.

Regarding the syntax, maybe we can have some symbol ∂Ω₁ = Boundary(Ω, "inflow") and then in bcs [u(∂Ω₁) ~ 1.0]? Or via some selector ∂Ω₁ = Boundary(Ω, [x₁ ~ 0.])

[xtalax]

I like it, I'd just really like to keep things as interoperable as possible, leaving any discretization specific stuff out of the system, but realise that this may be infeasible for more complex domains. Perhaps we have it as an opt in thing, where we encourage users to use domainsets where possible (especially for models going in the library), but allow GMSH objects where necessary. As for this particular vortex shedding example this domain could be constructed with DomainSets using setdiff, but it is a toy example.

Perhaps we also have utilities for automated meshing of domainsets domains with parameters specified at the discretization level.

Can you elaborate on what ∂Ω₁ = Boundary(Ω, "inflow") means, what are the arguments? Is "inflow" just a name? There also needs to be a way to split up boundaries in to sub sections which may have different conditions - perhaps this can be specified by face for product domains, and by submeshes with GMSH domains.
As for having "inflow" boundaries specifically, at the PDESystem level, I know this naming convention is common in CFD circles but it isn't universal - I'd like for any particular property of boundaries to be equational where possible.

[termi-official]

GMSH was just an example. We might want to think about some interface to take some general mesh object (Note: such an interface does not exist in Julia, yet). I mean, what infromation do we need to query from the domain? Btw having CSG capabilities in DomainSets is great!

Regarding ∂Ω₁ = Boundary(Ω, "inflow") my idea was that "inflow" is some string identifier for the subdomain. Here the (possibly discretized) domain Ω somehow needs to carry these identifiers associated to parts of the boundary. These identifiers will be specific to the mesh anyway. Taking this to the initial problem above, one example could be

@parameters x₁ f
@variables u(..) v(..)
Ω = (x₁ ∈ ClosedInterval(0,1)) # Assuming for this example that it has the subdomains "left" and "right".
∂Ωl = Boundary(Ω, "left") # This is just the point x₁=0.0
∂Ωr = Boundary(Ω, "right")  # This is just the point x₁=1.0
∫ = Integral(Ω)

# Bilinear form
a(u,v) = ∫(∇(u) ⋅ ∇(v))
# Linear form
b(v) = ∫(f*v)
# System
eqs_weak = [a(u(x₁),v(x₁)) ~ b(v(x₁))]
bcs = [u(∂Ωl) ~ 0., u(∂Ωr) ~ 0.] # Now we use the boundary object instead of describing the the boundary analytically

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

Alternatively we could think about choosing the boundary via functions.

∂Ωl = Boundary(Ω, [x₁ ~ 0.])
∂Ωr = Boundary(Ω, [x₁ ~ 1.])

Or marking adding the meta information to the domain object by some API like

mark_boundary!(Ω, [x₁ ~ 0.], "left")
mark_boundary!(Ω, [x₁ ~ 1.], "right")

[xtalax]

We might want to think about some interface to take some general mesh object

Agreed, we probably want a wrapper type around some julia meshing library (any suggestions?) which is robust to different types of input, but we can start simple.

what infromation do we need to query from the domain?

We're going to need it's boundary at the very least, dimensionality, and sooner or later we're going to need the whole mesh itself. Again I support doing this with a wrapper type so we can use the DomainSets interface to give it the symbols corresponding to each dimension, but I think that's it. Any thoughts here?

struct MeshDomain{T}
     mesh
end
Base.eltype(::MeshDomain{T}) where T = T

function Base.in(a::AbstractVector, d::MeshDomain)
     # Perhaps we implement this by tracing rays and counting intersections with the boundary?
end

What do you mean by CSG?

I like your first syntax for the boundaries, do we need more than this?

struct Boundary{P}
    domain::P
    name::Union{Symbol, String}
    Boundary(domain, name) = new{typeof(domain)}(domain, name)
end

[termi-official]

Agreed, we probably want a wrapper type around some julia meshing library (any suggestions?) which is robust to different types of input, but we can start simple.

I think our best bet for now would be GeometryBasics.jl, but we do not have a general interface for nonlinear geometry yet (see JuliaGeometry/GeometryBasics.jl#161 may need to cycle back to this issue sooner or later).

Regarding the first code example, I am not wure what the role of Base.in is. Also, I am not sure if we really need the whole mesh itself inside of PDEBase - I think this is an object on that we query some information about boundaries and discretization details like for example which geometry has been used for the discretization (e.g. quadrilaterals) to setup the equations. Maybe we should have a more open interface like abstract type DiscretizedDomain{T} end with some open dispatches for the queries which we need internally when setting up everything.

What do you mean by CSG?

The ability to describe the domain analytically with union, setdiff and intersection operations on primitives.

For the second example I think this is fine, we should be just careful with the wording. Maybe something like "descriptor" instead of "domain" is more accurrate.

[xtalax]

If we can have something that describes nolinear geometries in a pure math form I'm all for going that way.

Base.in is just needed to give the domains their symbols, like

domains = [x₁ ∈ Interval(0.0, 1.0),
           x₂ ∈ Interval(0.0, 1.0)]

or

domains = [(x_1, x_2) \in ProductDomain(Interval(0, 1), Interval(0, 1))]

If you have some ideas for how you think this should work, please write a mock up - we're getting close I think to opening the PR so running code is becoming important.

For the second example I think this is fine, we should be just careful with the wording. Maybe something like "descriptor" instead of "domain" is more accurrate.

I agree here

[termi-official]

If we can have something that describes nolinear geometries in a pure math form I'm all for going that way.

Tried to tackle this several times, but had not yet found a good angle on this. Basically, what I wanted to do in GeometryBasics.jl is to derive some general interface, such that we can grab domain (and solution field) information for visualization purposes, so not everyone needs to write their own visualization tool (and same for standard pre/post-processing tools). The problem is that we have different kinds of non-linearities which may need slightly different interfaces to be really efficient. I will try to cycle back to this if I find some spare time.

Regarding the overload of Base.in I am still a bit lost. Let me sketch what I have in mind

import IntervalSets: Domain # This is the type used in DomainSets.jl and transiently Symbolics.jl

abstract type DiscretizedDomain{T} <: Domain{T} end # This should allow the correct dispatches via https://github.com/JuliaSymbolics/Symbolics.jl/blob/561e917aa7417aa24e3f6522b8b978f2f6ec3a1c/src/domains.jl#L4-L7 

struct MeshDomain{T} <: DiscretizedDomain
    mesh
end

With this in and `eltype should be correctly dispatched, or am I missing something?

[xtalax]

import IntervalSets: Domain # This is the type used in DomainSets.jl and transiently Symbolics.jl

abstract type DiscretizedDomain{T} <: Domain{T} end # This should allow the correct dispatches via https://github.com/JuliaSymbolics/Symbolics.jl/blob/561e917aa7417aa24e3f6522b8b978f2f6ec3a1c/src/domains.jl#L4-L7 

struct MeshDomain{T} <: DiscretizedDomain{T}
    mesh
end

Yes, you are right, you just dropped a T, I had gotten confused about the requirements.

[vpuri3]

@xtalax check this out for general code structure. Adding mesh domains is a great idea. You also want to be able to add metadata, boundary tags, etc

https://github.com/vpuri3/AbstractPDEInterfaces.jl/tree/master/src/Domains

[xtalax]

@vpuri3 can you expand on what you mean by boundary tags?

[koehlerson]

I would like to add something very important for solid mechanics. In solid mechanics there are very elaborate constitutive laws where a certain part of the gradient is obeyed by an evolution law. So,

# 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.]

is not true anymore since your grad has a more elaborate structure, in the simplest case it is

grad(var, rem) = Differential(x₁)(var) + rem

and another equation needs to be added

eqs = [div(grad(u,rem)) ~ f,dt(rem) ~ law(iv)]

+- syntax that I used wrong. Here, rem corresponds to something that is usually called internal variable iv which is a variable of the system that does not live at the nodes of the discretization but rather only at the integration points. Often they depend again on the gradient. The hello world example for this is probably scalar valued isotropic damage as e.g. https://gridap.github.io/Tutorials/dev/pages/t010_isotropic_damage/

Here, your gradient remains the same but the flux takes the form (1-D)* (materialparameter * grad(u)) and your equations are basically

eqs = [div((1-D)*materialparameter*grad(u)) ~ f, dt(D) ~ something(grad(u))]

[xtalax]

Would it be possible to have this more elaborate grad defined as a closure around the base grad? I favour keeping the base operators as simple as possible.

[termi-official]

Thanks for the reminder here Maxi. Let me try to elaborate on the issue further.

In solid mechanics we define the unknown displacement $u$ with the associated deformation of the domain. Assuming we have some undeformed reference configuration with coordinate $X$ and some deformed configuration with coordinate $x$, then our unknown displacement can be written as $u = x - X$. However, this induces an ambiguity. The gradient can be taken w.r.t to $x$ or w.r.t $X$. Both are used, depending on the specific formulation of the problem. Does this help?

Regarding the second point things are actually a bit more complicated. The "placement" of the "internal problems" heavily depends on the chosen discretization scheme. The most common one is to associate the internal problems with the chosen quadrature space (used to numerically integrate remaining integrals), because it is the canonical choice in the a large class of discretizations schemes (simply because we do not need to interpolate these values). For more details I leave a paper by Krishnamoorthi here (Numerical quadrature and operator splitting in finite element methods for cardiac electrophysiology) showing two different possibilities.

[koehlerson]

The ambiguity that Dennis mentioned is yet another problem in solid mechanics.
The more complicated gradient definition (ignoring the problem with which coordinate you take the gradient) will always remain a problem dependent on yet another equation that in the worst case lives internally in the integration points. Typical examples are plasticity, viscoelasticity, growth modelling.

[xtalax]

@koehlerson Perhaps we pencil in handling this for a subsequent iteration, can you create a seperate issue detailing this problem on deforming grids?

[xtalax]

Hey guys! Have just been shown this thread, very relevant to our discussion here: https://discourse.julialang.org/t/advanced-tutorials-pinns/80595/4

Particularly important points:

• Mesh format must be compatible with gridap
• Looks like there is an Nvidia mesher that supports .stl files, it would be good to support this too. I have actually already written an STL import function in my first ever julia project, so I could adapt that for this purpose
• GMSH compatibility for sure

[termi-official]

I would rather argue that instead of targeting a specific mesh format we should think about what the interface to a mesh should look like, so we can target compatibility with multiple frameworks (and not only GridAP). Note that the Gmsh compatibility is usually tied to the specific framework (see e.g. https://github.com/gridap/GridapGmsh.jl , https://github.com/Ferrite-FEM/FerriteGmsh.jl , Trixi.jl via P4est IIRC,and there are more - just follow the link here https://github.com/JuliaPDE/SurveyofPDEPackages which you referred to above). Gmsh can also handle meshing via .stl files.

[xtalax]

Great to know, thanks for educating me. I agree, it would be better to be as flexible as possible.