PDESystem boundary specification

[ChrisRackauckas]

@thomasgibson @sandreza @simonbyrne let's discuss how the PDESystem specification needs to change in order to accommodate the types of meshes that are used for Clima.

@KirillZubov @emmanuellujan should stay informed here, since we will need to make sure the PINN and finite difference methods update for any change.

The current PDE system is like:

@parameters x y θ
@variables u(..)
@derivatives Dxx''~x
@derivatives Dyy''~y

# 2D PDE
eq  = Dxx(u(x,y,θ)) + Dyy(u(x,y,θ)) ~ -sin(pi*x)*sin(pi*y)

# Boundary conditions
bcs = [u(0,y) ~ 0.f0, u(1,y) ~ -sin(pi*1)*sin(pi*y),
       u(x,0) ~ 0.f0, u(x,1) ~ -sin(pi*x)*sin(pi*1)]
# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0),
                   y ∈ IntervalDomain(0.0,1.0)]

pde_system = PDESystem(eq,bcs,domains,[x,y],[u])

One of the difficulties that I saw was that, in this form, it's "easy" to know how to identify the boundaries and write a condition to them. However, a breakdown occurs on harder domains. An even simpler example than unstructured meshes is CircleDomain, we we can no long write things like u(1,x) = ... for all x, since we now have to say things like u(x,y) = ... for all x^2 + y^2 = 1 which is a fundamentally implicit specification, whereas it can be made explicit only by directly knowing a representation of the boundary.

So it looks like we need to be able to query for boundaries of domains and specify the BCs based on them. For example, maybe something like:

domains = [x ∈ IntervalDomain(0.0,1.0),
                   y ∈ IntervalDomain(0.0,1.0)]

bcs = [boundary(domains[1])[1] => 0.f0, 
           boundary(domains[1])[2] => -sin(pi*1)*sin(pi*y),
           boundary(domains[2])[1] => 0.f0,
           boundary(domains[2])[2] => -sin(pi*x)*sin(pi*1)]

That looks a bit awful though, but gets the point across. But the real question is:

1. Does this solve the issues with BCs?
2. Is there a nicer way to get the same information?

[ChrisRackauckas]

Probably good to add @kburns here as well,

[kburns]

On the spectral side, Dedalus requires that the boundaries be coordinate level sets. So for a disk, polar coordinates are used and the boundary conditions are specified at fixed polar radius (u(r=R) = 0, etc.).

This works well for the simple domains we handle with global spectral bases (intervals, disks, annuli, spheres, spherical shells, balls, etc.) and products of such. Special attention is needed at the discretization level to deal with corners and coordinate singularities, though.

Smooth deformations of those domains can also be defined and incorporated into the PDE, so you can do things like ocean bottom topography by using stretched / terrain-following coordinates. But for anything more complicated, we usually switch to immersed boundary methods.

[emmanuellujan]

Because boundary conditions can be defined at any point in the domain, not just at the outer boundaries, I believe the best way to express them is through implicit expressions. For example:

bcs = [ (u(x,y)=x*y, x<0.5 && x^2+y^2=1.0),... ]

Unfortunately, it becomes more complex when something like this is needed:

F(x,y,z) = x^2+y^2+z^2
q(x,y,z) = -k*grad(u(x,y,z))
n(x,y,z) = grad(F(x,y,z)) / abs(grad(F(x,y,z)))
bcs = [ (dot(q(x,y,z),n(x,y,z))=0.0, x<0.5 && F(x,y,z)=1), ... ]

[emmanuellujan]

Also, if you have a point set (e.g. 3d reconstruction of an airplane wing) through which you can create a surface (or volume), you may need something like this:

points = [(138.9,177.8,10.5), (139.9,178.8,10.5), ... ]
S = calcSurface(points)
q(x,y,z) = -k*grad(u(x,y,z))
n(x,y,z) = calcNorm(S)
bcs = [ (dot(q(x,y,z),n(x,y,z))=0.0, (x,y,z) ∈ S), ... ]

[valentinsulzer]

I would suggest allowing a more "vector-like" representation, while still allowing the existing style for existing domains.
Cartesian example:

xvec = [x,y]
domain = xvec in UnitSquareDomain()
domains = [xvec in domain]
S = surface(domain)

# Existing example
# S[1], ..., S[4] pre-defined and documented for different domain geometries
bcs = [u(S[0]) ~ 0.f0, u(S[1]) ~ -sin(pi*1)*sin(pi*y),
       u(S[2]) ~ 0.f0, u(S[3]) ~ -sin(pi*x)*sin(pi*1)]
u(S[0]) == u(x,0) # true

# Dirichlet everywhere
bcs = [u(S) ~ 0.f0]

# Neumann everywhere
n = unit_normal(domain)
bcs = [grad(u(S)).n ~ 0.f0]

Spherical example

@parameters r
domain = AxisymmetricSphereDomain(0,1)
domains = [r in domain]

S = surface(domain)
n = unit_normal(domain)

# Dirichlet everywhere
bcs = [u(S) ~ 0.f0]
# or equivalently
bcs = [u(1) ~ 0.f0]

# Neumann everywhere
bcs = [grad(u(S)).n ~ 0.f0]
# or equivalently
bcs = [Dr(u(1)) ~ 0.f0]

For more complex geometries, in the python-based FEM packages they usually create the domain (and mesh) in a GUI, label the domain boundaries 1-n, then in the code specify what bc to apply on each boundary

[ChrisRackauckas]

The suggested syntax is in the WIP part of the documentation, i.e. https://mtk.sciml.ai/dev/systems/PDESystem/#Domains-(WIP)-1

t ∈ IntervalDomain(0.0,1.0)
(t,x) ∈ UnitDisk()
[v,w,x,y,z] ∈ VectorUnitBall(5)

with a suggested form for LevelSet for now as well. That seems to be equivalent to what @tinosulzer suggested. Now the difficulty here is that each of these objects will need to have an equivalent object for boundary(UnitDisk()) == UnitCircle(), which is not included in DomainSets.jl. @daanhb is there a plan for that? That seems to match what @tinosulzer wants with surface, but just calling it boundary would be simpler.

I don't think bcs = [u(S) ~ 0.f0] makes all that much sense, because you might have to mix a sphere with time and etc. It probably needs to be like:

bcs = [(u(x,y,z,t) ~ 0.f0, (x,y) ∈ boundary(d1), z ∈ boundary(d2))]

Even worse, we will need to handle:

bcs = [(u(x,y,z,t) ~ 0.f0, (x,y,z) ∈ boundary(d[1])]

and such, since there may be multiple boundaries that for a given domain.

I like the grad idea with unit_normal normal though. That's a syntax error though.

S = (x,y,z)
n = unit_normal(domain)
grad = Differential(n)
bcs = [(grad(u(S)) ~ 0.f0, S ∈ boundary(d))]

seems reasonable.

[daanhb]

There is some initial support for boundary in DomainSets, e.g. it actually is currently the case (with JuliaApproximation/DomainSets.jl#78) that

julia> boundary(UnitDisk()) == UnitCircle()
true

And it is a goal to make this work in more cases, though right now in spite of its name DomainSets does not actually define a great number of domains yet :-) It would probably need polygons and polytopes in particular.

As mentioned above, this approach will likely remain restricted to fairly simple domains and combinations of them. Still, getting all the feasible simple combinations right, with minimal assumptions on later usage, is within scope.

Right now, I'm aiming to make DomainIntegrals useful for DomainIntegrals.jl, for the evaluation of integrals on general domains. This means that I am thinking about ways of automatically generating parameterizations in terms of standard intervals and boxes, and it looks like that might work soon. I have not given any thought to friendly notations like shown above.

At one point we did have a basic PDE solver using domains from DomainSets: example notebook (probably no longer functional). Nothing here is symbolic or descriptive, it starts from a specific basis of a certain length (like a Fourier series), and everything later on are operators acting on the coefficients in that basis. Still, this solver only required generating points on the boundary, which was relatively straightforward to do. Continuous representations are a little harder.

[daanhb]

To give you an idea, right now I get something like this:

julia> d = ProductDomain(0..1, UnitDisk(), 2..3.0)
0.0..1.0 x the 2-dimensional closed unit ball x 2.0..3.0

julia> boundary(d)
the union of 3 domains:
	1.	: the union of 2 domains:
	1.	: Point{Float64}(0.0)
	2.	: Point{Float64}(1.0)
 x the 2-dimensional closed unit ball x 2.0..3.0
	2.	: 0.0..1.0 x the unit circle x 2.0..3.0
	3.	: 0.0..1.0 x the 2-dimensional closed unit ball x the union of 2 domains:
	1.	: Point{Float64}(2.0)
	2.	: Point{Float64}(3.0)

The output is messed up, but the result is some kind of union of product domains. The boundary of an interval is currently a union of Point's, that's where some of the ingredients are coming from.

[ChrisRackauckas]

Awesome, that sounds great. And getindex on a union of domains returns the domains in order? That would be the one sugar you didn't show. If that's the case, this all sounds like a plan.