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