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 https://github.com/JuliaApproximation/DomainSets.jl/pull/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.

[daanhb]

Well, glad that this may be helpful! For historical reasons I tend to use elements and element(structure, i) for composite structures, because in some cases getindex already has a different meaning for the structure. Anyway, in the example above, the structure is nested:

julia> d = ProductDomain(0..1, UnitDisk(), 2..3.0)

julia> bnd = boundary(d);

julia> element(bnd, 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

julia> element(bnd, 2)
0.0..1.0 x the unit circle x 2.0..3.0

julia> element(bnd, 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)

julia> element(element(bnd, 3), 3)
the union of 2 domains:
    1.  : Point{Float64}(2.0)
    2.  : Point{Float64}(3.0)

It would probably be better to turn a product domain of unions into a union of product domains, I'll look into that. Still, for boxes and cubes the output is quite horrible at the moment with all the Point's - it would be more useful to have a union of faces or edges, but those don't exist yet.

[ChrisRackauckas]

What is getindex used for there?

[daanhb]

It is not used in this package, but the way it was designed anything that supports in can act like a domain. This includes Set's and Array's. It also felt funny to define the size of a domain as being its number of composite elements rather than, you know, its size :-) Still, there is probably no direct clash with UnionDomain and ProductDomain. Perhaps I should reconsider.

[ChrisRackauckas]

It would make a lot of nice things like broadcast work.

[daanhb]

I've opened a separate issue for the getindex behaviour (https://github.com/JuliaApproximation/DomainSets.jl/issues/79). There may be a conflict because currently we use broadcasting for domains differently, to support the easy creation of mapped domains. For example, this would be a ball with different radius and center:

julia> [0.5, 0.2, 3] .+ 2*UnitBall()
A mapped domain based on the 3-dimensional closed unit ball

Here we think of a domain as being a continuous set of points, and translation is like element-wise addition. I'm mentioning the syntax because it may be relevant in this issue, and suggestions for user-friendly creation of domains are very welcome! I've been meaning to explore something like

[ [0.5, 0.2, 3] + 2*x for x ∈ UnitBall()]

That seems closer to the kind of syntax we're looking for here. Apart from Ball(radius = 2, center = [0.5, 0.2, 3]) perhaps.

[ChrisRackauckas]

Hmm, that's an interesting use case for that. Cool!

[ChrisRackauckas]

Awesome. I think we have a general plan. We just need to get the PINNs and MOLFiniteDifference on it in order to solidify it.

[valentinsulzer]

@daanhb how would you do the equivalent of this using DomainSets.jl?

domain = IntervalDomain(0.0, 1.0)
domain.lower
domain.upper

and, how would you define a UnitDisk of radius r, then later read the radius? UnitDisk(2.5) gives an error.

[daanhb]

Hi @tinosulzer, the equivalent would be:

julia> using DomainSets

julia> d = 0..1.0
0.0..1.0

julia> infimum(d)
0.0

julia> supremum(d)
1.0

julia> DomainSets.endpoints(d)
(0.0, 1.0)

julia> 2.5 * UnitDisk()
2.5 * UnitDisk()

julia> 2.5 * UnitDisk() .+ [0.5,0.1]
2.5 * UnitDisk() .+ [0.5, 0.10000000000000002]

julia> typeof(ans)
MappedDomain{StaticArrays.SVector{2, Float64}, DomainSets.GenericAffineMap{StaticArrays.SVector{2, Float64}, Float64, StaticArrays.SVector{2, Float64}}, StaticUnitBall{StaticArrays.SVector{2, Float64}, :closed}}

Some comments:

• it is called infimum and supremum because the intervals may be open or closed. If that added complexity is not needed, you can use endpoints (which is currently not exported - but that can change, please file issues at will)
• a UnitDisk always has radius 1. We focused on a few primitives when building the framework, and the idea is that domains do not always correspond to a type but rather to behaviour. A disk is a disk if in(x, domain) is true when norm(x) <= 1, regardless of the type of the domain. A UnitDisk has no state, so it is a very cheap type to use.
• Right now, 2.5 * UnitDisk() makes a mapped domain. But there could easily be a function disk(...) function with arguments for radius and center, and there could be a more general Disk type. That could be useful if you want to dispatch on a domain being a disk, but in that case you may not catch all domains that are actually a disk. As you can see, the type of the mapped domain is not very pleasant to look at.

The release of DomainSets 0.5 is imminent and it will be a lot more robust. The set of domains it actually supports right now is still pretty limited, but it can be extended either in the package or elsewhere - again, feel free to file issues!

[daanhb]

What other domains are there to include? I see a circle and product domain at least in domains.jl.

The AbstractDomain{N,T} in that file would correspond to a Domain{SVector{N,T}}. The abstract type Domain{T} is a bit more general than that, and univariate domains typically have T=Float64. For example, an interval is usually a Domain{Float64} and not a Domain{SVector{1,Float64}} so it may be safer to work with Domain{T} from the start.

The interval notation 0..1 is actually defined in IntervalSets.jl. The .. operator does not have a strong precedence however, so you often have to write (0..1) instead when used in an expression.

Another pitfall is that 0..1 actually gives a Domain{Int}, and its endpoints will be integers. But the interval contains all numbers x that satisfy 0 <= x <= 1, including floats. So the T in Domain{T} is just some indication of what the elements of the domain are, there is no type strictness. The true definition of a domain is whatever in(x,d::Domain) returns true for, which is closer to the mathematical definition. Typing 0..1.0 does give an Interval{Float64}. You can also write Interval(0.0, 1.0).

There is a UnitSimplex for now, but no general triangles or polytopes. It seems a useful library for that is Meshes.jl, so it may be worth making sure that DomainSets and Meshes interoperate well rather than duplicating efforts.

Coming back to your question, recovering the radius of a circle is not actually possible, weirdly enough. Currently you'd have to check the mapped domain to see whether the domain being mapped is a circle, whether the map is linear and diagonal, and whether all diagonal entries are equal... :-)

Is there anything I can do here to help moving forward?

[valentinsulzer]

Thanks! For now I am happy with just cuboids, disks and spheres. The main thing to figure out here is what we want the user-facing syntax to be, then actual implementation is very easy (just replacing what's in domains.jl)

I like just being able to use Interval(0.0, 1.0) and those who want to be fancy can use 0..1.0. It is important that the type be Float64 and not Int when defining the time domain. And using endpoints instead of lower and upper should be fine too.

For the circle, I'm feeling stupid for asking how to change the radius of a unit disk haha. But it would definitely be useful to have a more general Disk, and to be able to read its radius.

but in that case you may not catch all domains that are actually a disk

why not?

[daanhb]

I like just being able to use Interval(0.0, 1.0) and those who want to be fancy can use 0..1.0. It is important that the type be Float64 and not Int when defining the time domain. And using endpoints instead of lower and upper should be fine too.

Yeah, 0..1 returning an interval of integers is annoying, but automatic conversion in IntervalSets.jl is not ideal either because some users might actually want to create intervals of integers. Perhaps it is possible to catch this case in post-processing? The statement convert(Domain{Float64}, interval) should do the trick.

but in that case you may not catch all domains that are actually a disk

why not?

It might not be a problem in practice, depending on what interface you provide here. But DomainSets is ridiculously general and so it allows many different ways to create a domain that acts like a disk. Even if you have a Disk type, someone might map a unit disk with a linear map, or an affine map, or... Or in 1D, a disk is really just an interval and vice-versa. You can't dispatch on all those things based on type alone if you want to catch a disk, that was my thinking.

I'll add a ball/sphere with radius and center in DomainSets 0.5.1 because it is such a common case (v0.5 is out today). Perhaps the type can stay the same after scaling and translation, because that would preserve the shape.

[mforets]

Hi there :wave: I'm coming late to the party but.. Note that some of the groundwork in LazySets.jl could be applied to this issue (and vice-versa: we're interested in PDE modeling combined with set-based techniques in Julia so I'm interested to follow and maybe at some point participate in this issue). There is an upcoming JuliaCon workshop for our JuliaReach stack [1]. We've been using Symbolics.jl variables to create lazysets types for a while now. I haven't looked at depth into DomainSets.jl but we could create new convert methods so it's easier for newcomers and to reuse more code.

[1] https://github.com/JuliaReach/JuliaCon-2021-Workshop-Its-All-Set

[ChrisRackauckas]

Perhaps it is possible to catch this case in post-processing?

There are some PDEs with integer-valued independent variables which are interesting, like the chemical master equation.

Note that some of the groundwork in LazySets.jl could be applied to this issue (and vice-versa: we're interested in PDE modeling combined with set-based techniques in Julia so I'm interested to follow and maybe at some point participate in this issue).

Where should I be looking for the high level overview of the domain interface?

[mforets]

Where should I be looking for the high level overview of the domain interface?

You could check [1-3].

Now let me state upfront that specialized algorithms to efficiently iterate over complex domains (say, a finite element mesh) are not the focus of LazySets.jl -- the focus is on writing new algorithms for set computations, conversions, approximations, and so on (especially in dimensions higher than 3). The library Meshes.jl will be much better suited for the purpose of mesh generation, exploration, and computational geometry à la CGAL.

[1] https://juliareach.github.io/LazySets.jl/dev/man/tour/#A-Tour-of-LazySets [2] https://juliareach.github.io/LazySets.jl/dev/man/tour/#Exploring-the-type-hierarchy [3] https://arxiv.org/abs/1901.10736

[valentinsulzer]

We're using DomainSets.jl for now after #976 was merged. Thanks @daanhb for adding the features for circle domains. Anything else in this issue or should we close it?

[ChrisRackauckas]

I don't think we're completely set until we're settled and using the boundary specifications as well: intervals and boxes can be easy on the boundary specifications but things like circles, not so much.

[daanhb]

I think the suggested syntax above should work, in which a variable is declared to be associated with a domain, but indeed it may require some experimenting to get the non-product domains right. It seems easier to use v ∈ UnitDisk() than (x,y) ∈ UnitDisk() in the definition. But when passing v to a function later on, you somehow need to keep track of whether that function expects one argument (v) or two (x and y).

I found that you have some control over this when using the native generator syntax, and I can use this to show what I mean. The generator does not associate a variable with a domain, but it associates a function with a domain. There are several scenarios:

julia> using DomainSets

julia> gen = (exp(x) for x in Interval(2,3))
Base.Generator{ClosedInterval{Int64}, typeof(exp)}(exp, 2..3)

julia> gen.f
exp (generic function with 15 methods)

julia> gen.iter
2..3

Here, the generator has two fields f and iter which give you precisely the function and the domain. (The variable x is not bound to anything and does not exist outside of the generator expression, so it could not be a ModelingToolkit parameter.)

A function on a domain with product structure could be specified with a repeated for:

julia> gen2 = (exp(x+y) for x in Interval(2,3.5), y in Interval(4, 5.5))
Base.Generator{Base.Iterators.ProductIterator{Tuple{ClosedInterval{Float64}, ClosedInterval{Float64}}}, var"#11#12"}(var"#11#12"(), Base.Iterators.ProductIterator{Tuple{ClosedInterval{Float64}, ClosedInterval{Float64}}}((2.0..3.5, 4.0..5.5)))

julia> gen2.f
#7 (generic function with 1 method)

julia> gen2.f( (0.2, 0.3))
1.6487212707001282

julia> exp(0.2+0.3)
1.6487212707001282

julia> gen2.iter
Base.Iterators.ProductIterator{Tuple{ClosedInterval{Float64}, ClosedInterval{Float64}}}((2.0..3.5, 4.0..5.5))

julia> D = ProductDomain(gen2.iter.iterators...)
(2.0..3.5) × (4.0..5.5)

julia> p = point_in_domain(D)
2-element StaticArrays.SVector{2, Float64} with indices SOneTo(2):
 2.75
 4.75

julia> gen2.f(p)
1808.0424144560632

Here, the function gen2.f expects a tuple of two arguments (or any t that defines t[1] and t[2] I think), but they are named x and y in the expression. The iter field stores a tuple of iterators, which correspond exactly to the respective domains, and hence you can make a ProductDomain out of it. The eltype of that domain would be an SVector{2}, which you can pass to gen2.f.

Alternatively, you can work with vectors to support non-product domains:

julia> gen3 = (cos(x[1]+x[2]+x[3]) for x in UnitBall())
Base.Generator{EuclideanUnitBall{3, Float64, :closed}, var"#13#14"}(var"#13#14"(), UnitBall())

julia> gen3.f
#13 (generic function with 1 method)

julia> gen3.iter
UnitBall()

Here, the function gen3.f again expects a single argument, but the generator expression also treats x as a vector. The domain is simple again, gen3.iter literally is the unit ball.

Finally, the trickiest setting is where you combine both structures, which probably arises e.g. for time-dependent problems:

julia> gen4 = (exp(x[1]+x[2]+y) for x in UnitDisk(), y in Interval(2..3.0))
Base.Generator{Base.Iterators.ProductIterator{Tuple{EuclideanUnitBall{2, Float64, :closed}, ClosedInterval{Float64}}}, var"#15#16"}(var"#15#16"(), Base.Iterators.ProductIterator{Tuple{EuclideanUnitBall{2, Float64, :closed}, ClosedInterval{Float64}}}((UnitDisk(), 2.0..3.0)))

julia> gen4.f
#15 (generic function with 1 method)

julia> gen4.iter.iterators
(UnitDisk(), 2.0..3.0)

julia> D1 = ProductDomain(gen4.iter.iterators)
UnitDisk() × (2.0..3.0)

julia> eltype(D1)
StaticArrays.SVector{3, Float64} (alias for StaticArrays.SArray{Tuple{3}, Float64, 1, 3})

julia> D2 = TupleProductDomain(gen4.iter.iterators)
UnitDisk() × (2.0..3.0)

julia> eltype(D2)
Tuple{StaticArrays.SVector{2, Float64}, Float64}

julia> p = point_in_domain(D2)
([0.0, 0.0], 2.5)

julia> gen4.f(p)
12.182493960703473

Now, the function gen4.f expects a nested tuple of the form ( (x[1],x[2]), y). As before, gen4.iter.iterators completely specifies the domain. If you pass it to the ProductDomain constructor, it will make a 3D Euclidean domain with eltype SVector{3}, because that is what you usually want in geometry. However, you can explicitly make it a TupleProductDomain instead: this constructor does not concatenate vectors and scalars, but simply joins the element types of the respective domains into a tuple. The result has the right structure to pass to gen4.f.

The code for product domains in DomainSets is a little tricky because I wanted to support all possibilities. The good news is that probably makes it flexible enough to get everything to work, the bad news is that it becomes, well, tricky.

[daanhb]

I think the boundary function in DomainSets still needs some improvement. Some more work will be needed to support gradients and directional derivatives. One concern is that domains in DomainSets are not oriented, though implicitly they do have an orientation. I noticed that ApproxFun defines an oriented domain to get around this, with which you can wrap a domain and change its (default) orientation.

[ChrisRackauckas]

Gradients in the normal direction will be required for Neumann BCs.

I kind of think the generator use is a little too cute, but I still don't know.

[daanhb]

True, cute, but probably too limited. I wanted to show what considerations arise when treating functions defined on non-product domains - these issues or similar ones will show up one way or another, I think.

[xtalax]

Can we move this issue to PDEBase? @ChrisRackauckas Looks like a lot of stuff that will need to be thought about eventually is documented here.

Also PS what is Dedalus?