Use in ApproxFun

[dlfivefifty]

As I (very slowly) revamp ApproxFun.jl using InfiniteLinearAlgebra.jl and ContinuumArrays.jl I'm starting to think how solving BVPs in ApproxFun.jl will work, and likely the best option is to support problems specified via ModelingToolkit.jl. This issue is to begin discussion on how that would work.

I'm imaging for example the Airy equation would be specified as something like:

@parameters x
@variables y
D = Differential(x)
eqs = [y(-1) ~ airyai(-1), y(1) ~ airyai(1), D(D(y)) ~ x*y]
sol = solve(ContinuousNonlinearProblem(eqs, (-1,1)), ApproxFunNewton()) # returns a Fun

Here I called in NonlinearContinuousProblem as we may support, for example, integral operators. (Of course this example is actually linear but whatever, would also support nonlinear obvs.)

I suppose in ModelingToolkit we can compute the Jacobian symbolically so that I don't need to use something like the current DualFun for determining the Jacobian by auto-diff. It's not yet clear to me though how one would actually implement the solver.

[ChrisRackauckas]

I think it would be good for ApproxFun to target PDESystem. See NeuralPDE.jl as the most developed consumer of that right now:

https://neuralpde.sciml.ai/dev/pinn/system/

Where new packages tie in is by defining new discretize dispatches. I planned to setup an interface so that discretize(sys,ApproxFun()) would kick out a nonlinear system or an ODEProblem for the spectral-discretized PDE. Also note we plan to add integral operators for this fairly soon.

[dlfivefifty]

ApproxFun() won't work as its the name of the module... but can look into it.

How advanced are things on other geometries like disks and spheres? At the moment that's my main coding focus.

[ChrisRackauckas]

How advanced are things on other geometries like disks and spheres?

See https://github.com/SciML/ModelingToolkit.jl/issues/526. We just need to settle on one, or two, or 3 boundary forms that we think multiple discretization libraries can handle. I'm thinking we might want:

1. Boxes and intervals
2. Circles and spheres
3. Arbitrary combinations and unions of these
4. Level sets

But we can continue figuring that out. Slapping a sphere type down is simple, but making it so we can then take the union of a sphere with an off centered box is a bit less so.

[dlfivefifty]

Why not use DomainSets.jl which has support for UnionDomain, cartesian products, etc.?

[ChrisRackauckas]

Because I didn't know about it 😄 . That's perfect for what we need. Any way it can do level sets too?

[dlfivefifty]

Not at the moment but the "domain interface" is essentially just implementing in, via indomain. It would be pretty easy to add:

struct LevelSet{T, F} <: Domain{T}
    f::F
end

indomain(x::T, d:LevelSet{T}) where T = iszero(d.f(x))

So, for example, the unit circle would be equivalent to LevelSet{SVector{2,T}, typeof(norm)}, though we of course already have a special type UnitCircle.

Though perhaps you want f to be symbolic as opposed to a function?

@daanhb have you thought about this?

[ChrisRackauckas]

Though perhaps you want f to be symbolic as opposed to a function?

Yes, allowing it to be symbolic could potentially be helpful.

[ChrisRackauckas]

So, for example, the unit circle would be equivalent to LevelSet{SVector{2,T}, typeof(norm)}, though we of course already have a special type UnitCircle.

Yes, I think it would make sense for "most" users to try and use the more specialized forms when they can, as more discretizers can support them. But there are definitely cases like physics-informed neural networks and a discretizer wrapping Gridap.jl where we people may want to try FEM and PINNs simultaneously on some larger sets of domains, and level sets seems to be a common one. The other of course is triangularized domains, which it seems you cover via unions of simplexes.

[dlfivefifty]

Unions of simplices is what we do at the moment but it's not a good way to handle meshes. I'm sure there's another package for simplicial meshes.

Note for many applications you also need a simplicial complex, that is with knowledge of the edges and nodes in addition to the triangles. This is the case for domain decomposition methods which involve operators on the interfaces. A quick google reveals that Simplicial.jl is in active development.

[ChrisRackauckas]

Then I think it's a plan to start using DomainSets.jl everywhere.

[daanhb]

Great, this would help push development of DomainSets.jl to cover the new use cases. It won't do everything that is needed yet, but it looks like, indeed, it does at least some of it.

[daanhb]

Not at the moment but the "domain interface" is essentially just implementing in, via indomain. It would be pretty easy to add:

struct LevelSet{T, F} <: Domain{T}
    f::F
end

indomain(x::T, d:LevelSet{T}) where T = iszero(d.f(x))

So, for example, the unit circle would be equivalent to LevelSet{SVector{2,T}, typeof(norm)}, though we of course already have a special type UnitCircle.

Though perhaps you want f to be symbolic as opposed to a function?

@daanhb have you thought about this?

They have been on my mind, yes, but it was never clear what to implement. It is hard to offer generic functionality without analytical knowledge of the function f in f(x)=0. The in function like above is typically not that useful, because the domain has virtually no volume. The approx_in function is more interesting, and one could conceivably do abs(f(x)) < epsilon generically, with some assumptions on f.

In applications (an integral equation solver) I have wanted a parameterization of the level set, or at least some way of generating points in the set. Currently, this does work for a circle, because the type "knows" what it is: you can ask for its parameterization, which is a map from some other domain (in this case an interval) to the level set. That approach works, but the interface is still a bit clumsy and specific, especially with the maps involved.

Maybe if f is a symbolic function, this could give the right knowledge - I did not look into this. Perhaps with a specific use case in mind, it will also be clearer what can be done.

[ChrisRackauckas]

Let's continue the talk about domains in https://github.com/SciML/ModelingToolkit.jl/issues/526 . I think we're agreed on using DomainSets.jl now, but I think we'll need a bit more added to DomainSets.jl

[dlfivefifty]

I can close this issue for now