DiffEq's PDE Story 2nd Try: Domains, BCs, Operators, and Problems #260
Comments
|
To make this complete, let me describe the highly related DiffEqOperators interface. Basically, a DiffEqOperator is the discretization of an operator to full(L) # Get a dense matrix
sparse(L) # Get a sparse matrix
banded(L) # Get a banded matrix
...it can possibly update_coefficients!(L,u,p,t,x)which updates internal coefficients of L = a(u,p,t,x)*Asome quasi-semi-linear PDE discretization operator with a matrix The reason why this is nice is that quasi-semi-time_dependent operators can all be appropriately evaluated via the sequence: update_coefficients!(L,u,p,t,x)
L*uHave a linear PDE and want the stationary solution (i.e. it's not time-dependent)? update_coefficients!(L,nothing,p,nothing,x)
# Solve L*u = bWant to write an Euler loop for a quasilinear PDE? for i in 1:10
t = t + dt
update_coefficients!(L,u,p,t,x)
u = u + dt*L*u
endEtc. You see that these codes are operator independent, and independent of the implementation of the underlying operator. All it needs to know is that We plan to support special matrix types like BlockBandedMatrices.jl as well if that can be done generically, or that can be dispatching on the DiffEqOperator concrete type. |
|
Looks like a great start. I will put up more thoughts soon but a few to start... some are on general framework while some are on the gen:
|
|
Some comments:
|
|
Sorry, just say that you answered the discretization question, but would love to see examples of |
|
It would be nice to have Then the key thing I am missing is whether the I don't have strong feelings about the separation of the DerivativeOperator and the particular discretization, but the space vs. domain would be helpful. |
Yes, that was just for space(lbc,rbc,tbc,bbc,corner1,corner2,corner3,corner4)would be what's required. Documenting this needs to be done by pattern for nD domain/space pairs, so there must be a better idea. But that crucial step is why I wanted to leave open the possibility to separate domain and space, maybe it's easier like this? What's gained? I don't really know that fully yet.
I agree, but I didn't want the lack of a good name to pause the release of the idea any more.
In case it's for a whole side of the domain. Left of a 2D grid needs a function.
It has
I'm sure most people who read this won't be using that part, just the operators (or they will be building this level of the abstraction 👍 ). But this is how we support users who say "how do I solve the Diffusion-Advection equation?" Obviously they can do it themselves by building the operators, but it's not hard with this tooling to give them a function that spits out an ODE to solve. Also, this is a great place for undergrad and GSoC projects IMO.
Yes. They will be there.
The "standard" one will be the current DiffEqOperators one. It might even good to have defaults for this. But for example, someone can write a package which computes fast Laplacians for regular girds, and so they can let
See the stuff on the
If I make it lower case, would it be easier to understand that it's supposed to be a high level "give me a discretization of this derivative to this order that's appropriate for this space using this method". |
|
It's not possible to know what boundary conditions are needed without knowing the differential operator, and even if you know the differential operator its non-trivial. |
|
I see what you mean. So then it'll just have to error at the problem/operator construction call instead of the space construction time, meaning there's no reason for |
|
Bcs should either be in the operator, or in the basis/space (for zero conditions). Putting them in the domain doesn't really make sense. Also, there's some work on a general purpose Domains.jl package: https://github.com/daanhb/Domains.jl |
|
For stochastic processes (I can't say much about physics applications) the boundaries are connected to the domain: reflecting or absorbing barriers. The problem with associating them with the operator itself is that Chris has the idea of composing the operators lazily. It is only at the level of the completely composed operator that you can apply the boundary values (for stochastic processes at least). Are we sure the operator composition approach would work (for boundaries which end up affine)? |
Yeah, this is where the idea comes from. The first approach I had in mind was: domain = RegularGrid((0,5),dx)
lbc = Dirichlet(0)
rbc = Neumann(g(t,x))
bcs = BoundarySet(lbc,rbc)
coeff(u,p,t,x) = x^2 + t^2
A = derivative_operator(domain,bcs,FiniteDifference(),coeff,2,2)
L = derivative_operator(domain,bcs,LazyUpwind(),coeff,2,2)but @jlperla asked whether there's a way to make sure the same BCs are satisfied, in which case I thought we might as well package the BCs into the "domain". Maybe that's still fine but we need a new name?
Yes. There's the odd issue that |
|
As soon as you associate the boundary conditions with the domain, it (1) breaks generality, as you can't impose constraints like \int_Ω u(x) dx = 0 and (2) doesn't make sense in the ∞ -dimension setting, where bcs are either side constraints, or incorporated in the function space. |
Wouldn't that just be
That's covered by the answer to the first, as either BCs are just a collection sent to every operator build (so it doesn't matter if there's some weird ones or none), or it needs to be at the operator instantiation call. |
|
I don't think I understand: an FD second-derivative maps |
|
Thinking less mathematically and more practically, A = derivative_operator(domain,bcs,...)
L = derivative_operator(domain,bcs,...)would I ever want those bcs (or more generally, constraints) to be different? If every operator has to always satisfy the same constraints, do they need to be kept separate? |
|
The main question I have is whether the generality to allow lazily composing operators is worth the complexity in interface complexity and short-term performance. With loop fusion, etc. it may be possible to generate an overhead free implementation, but it is going to be tricky. For economics and finance applications at least, there are very few finite difference discretizations we would need (e.g. for a diffusion process: central differences on a diffusive second derivative and upwind on the first-derivative drift). |
|
I don’t understand why the derivative takes in bcs: it seems to be conflating the application of the operator with the inversion of the operator. |
|
For example, we could have a specialized operator for a diffusion process which generates the coefficients given the passed in boundary functions: |
What can be dropped/simplified if we aren't composing operators lazily? None of the PDE story has that in it, that's the DiffEqOperator stuff which users can pretty much ignore. |
Even forward application of the operator requires the BCs. |
For me at least, I have no use for the higher level "heat equation" and "odeproblem" solution methods, so it is the What could be simplified:
Now, if in the underlying code for the operator you want to use lazily decomposed operators, that is an implementation detail, but we are discussing the higher level interface here. |
|
The derivative of a function can’t “see” a single point (or rather it’s just not defined at a point where it’s not differentiable), so I don’t agree that the notion of derivative is different for Dirichlet and Neumann. Maybe I’m missing the point of the conversation: is the point of this discussion to make a system for doing finite differences (and only finite differences), or is the point to make something general purpose where you can swap out the Discretization method? |
The stencil matrix still has to be modified at the edges to satisfy the BCs. That's what the naive loop is essentially doing.
Where do you have to worry about that? Show me an example of where you have to worry about the details of loop fusion. Just get and sparse the operator if all you want is the sprase operator?
Is this about the lazy implementation of DiffEqOperators via Fornberg's algorithm? That has nothing to do with the interface... but anyways, even that has been mostly built by undergrads so it's at least at that level. Most of the code is just a loop over a |
The latter. For example, this should be generic enough that I can have a wrapper over an ApproxFun domain, specify some BCs in the DiffEq style,
But the differential operator is the derivative at each point in the discretization. At some point that operator has to be instantiated, and in order to know its action it has to have the BCs in it. |
|
But in FD it shouldn’t be each point: each derivative drops one point. So D maps values at n points to values at n-1 points, D^2 maps values at n points to n-2 points and so on. |
|
Only in the periodic case (which really is a different domain since the topology is a 1-torus) do you get a map from n points to n points. |
|
When I said
What I meant was that if the user only directly works with fully composed operators in their higher level interface, then we don't have the current issue we are discussing with the boundary values being connected to the derivative. If the user works directly with the fully composed differential operator, then it makes sense to group that with the appropriate boundaries. |
Not necessarily. If you know the BCs then you can incorporate them into the operator to make the map not drop points. The periodic case is well known where you just put stencil coefficients in the upper right and lower left corner of the matrix (thinking about it as a matrix). But for Dirichlet/Neumann, in order for the BCs to be satisfied it defines what the values have to be at the boundaries, giving you an n -> n map by filling in the ends via the BCs. I'm trying to think of whether there's a case that can't be true for. |
|
That incorporation is exactly equivalent to a step of Gaussian elimination ... I don’t see why you would want to try doing Gaussian elimination before you’ve setup the problem .. |
|
For a user in my field, what they really want is to (1) specify the stochastic process and (2) specify what happens at the boundaries of that stochastic process . There are only two types of boundaries that matter: (a) a reflecting barrier and (b) an absorbing barrier. I would love to specify these things in an abstract way independent of the discretization, but it isn't necessary. In fact, there are practical considerations (i.e. in solving HJBE the drift changes all the time as part of the algorithm, and you solve for the new drifts at grid points) that make sticking with finite differences reasonable for this exact problem. So, in my way of thinking you are having me specify the https://en.wikipedia.org/wiki/Infinitesimal_generator_(stochastic_processes) for the stochastic process directly, which is fine. If you want me to specify it in parts and have lazy composition then I think that is fine as well, but somehow you are going to have to attach the boundary types to it when you generate the full matrix that puts things together. The general algebra behind that is not at all obvious to me, let alone the numeric implementation or what happens when I want to take the adjoint (as you do when moving between solving the KBE and KFE). |
|
So here's the gist of it. I think I have most of it clear now. Copying my notes for future reference. The domain is discretized and
Now let's say we want to solve A[:,bidxs]*(B[:,bidxs]\[bc1;bc2]) = A*Qband R*(A - A[:,bidxs*(B[:,bidxs]\B)) = A*QaI don't have a good way showing why that should be true though, so maybe @dlfivefifty can chime in on how that definition actually works. So that can be used to derive But once we're at this point, to solve the actual differential equation we only want to solve the interior, since the boundaries are given directly by the interior v_t = A*Qa*v + A*Qbis how the operator is applied. For a steady state problem, it's the same idea. Discretize some This shows that what's important and non-trivial is Dirichlet0. This case is easy. If you have any Dirichlet. The boundary is not dependent on any points in the interior, and therefore Neumann. You have to pick a discretization of Robin. Left as an exercise for the reader :P. It's just the same ideas from before to write out the affine relation. It needs a choice of And that's all of the FDM discretizations. Let's take a quick look at the spectral discretizations then. In this case, let We can actually make clear where the affine-ness comes from from this. If you let Of course this extends to finite element methods since it's usually discussed in an operator form much like this. So this covers the different discretizations as long as you can find |
|
I wanted to finish up with the discussion of how exponentiation of the operator works. Using the notation of the previous post, we have that our diffeq is given via: v_t = A*Q*vwhere v_t = L*vand it's clear that we just use But now let's take v_t = A*Qa*v + A*Qb + f(t,v)but then if we define v_t = L*v + f̃(t,v)which is the semilinear ODE with a square I think that covers the questions for how to get it right. Now we need to implement it all.... Reference: the Chebfun paper https://academic.oup.com/imajna/article/36/1/108/2363563 |
|
The abstract definition of the operator, along with domains, has been punted to SciCompDSL.jl. Operator instantiation will be the dispatches discussed above for DiffEqOperators.jl and DiffEqApproxFun.jl. |
|
At some point we should talk about unifying |
|
I think we're going to have to make such an interface in SciCompDSL.jl, so that's where it will likely happen. It's getting hard to predict though, I think it just needs to get done and then refactored after we have something working. |
|
Maybe the bests steps are:
|
|
Yup, sounds good. |
|
FYI I think this approach to boundary conditions works for eigenvalue problems too. |
|
I am not sure whether this approach works for general PDEs and discretisations. As far as I've understood the discussion above, the boundary conditions are imposed in a strong sense, i.e. the numerical solution has to fulfil the boundary conditions exactly. For hyperbolic conservation laws (linear advection, wave equation, Burgers' equation, Euler, ...), the imposition of boundary conditions might be a little bit different. Such a strong enforcement of boundary conditions can result in unstable methods. Therefore, summation-by-parts (SBP) operators using simultaneos approximation terms (SATs) have been designed to construct provably stable semidiscretisations (https://doi.org/10.1016/j.jcp.2014.02.031). The basic idea is to mimic integration by parts discretely by SBP operators, i.e. (finite difference) derivative operators that are defined on the complete grid, including the boundary nodes. In the interior, these derivative operators can be just the classical centered derivative operators. Near the boundary, special boundary stencils are chosen to mimic integration by parts discretely. Thus, the derivative is computed without using boundary conditions. To enforce boundary conditions, the SATs are added at the boundary nodes. In the simplest case, these terms contain some difference of the value of the numerical solution at the boundary and the desired boundary value. Thus, the numerical solution does not necessarily satisfy the boundary conditions exactly and the numerical solution has to be computed at all nodes (including the boundary nodes). However, this procedure allows the creation of provably stable semidiscretisations for hyperbolic equations. If I understood the discussion above correctly, I don't know how such an approach can be represented in the framework discussed above. For hyperbolic equations, well-posed boundary conditions might be specified only at certain boundaries and not everywhere (inflow or outflow). Moreover, whether it is possible to enforce (Dirichlet) boundary conditions can also depend on the numerical solution itself. But this might be handled by the corresponding semidiscretisation. I can see how the approach discussed above works for parabolic PDEs. Please, correct me if I'm wrong, but I don't see immediately how this approach can be extended to weakly imposed boundary conditions using SATs. |
|
@ranocha brings up a good point. In fact, what it basically points out is that this version of boundary condition handling we have been discussing is only compatible with linear boundary conditions, basically boundary conditions of the form
since then you can discretize the derivative terms and then the expansion operator has to be a linear operator. But if you have a global boundary condition or something which is nonlinear, like But I think it highlights that what really needs to be discussed is what we want out of the PDE story. (Note: I'm mostly using this to prepare my thoughts for today's workshop). Essentially, we want two things:
What we tried to layout was a formulation for (2), which ran into multiple issues which we could get different answers for, but I think that it's ultimately doomed to failure because the developers need the flexibility in order to handle whatever wonky PDEs and boundary conditions that they are interested in. So I think it makes sense to abandon (2) as seemingly impossible without giving up both performance and feature-completeness. FEM, FDM, etc. libraries are going to have to act differently in order to exploit their unique features, and anyone working at that level will need to make tough but isolated decisions. DiffEq can hold a variety of tools like this, with DiffEqOperators being the FEM operators with linear boundary conditions which is core to many discretizations, while others like ApproxFun play a central role in the surrounding ecosystem as well. But (1) is the goal we should be trying to attain. The question is stated as: how can someone specify a high level what PDE they want to solve without simultaneously specifying a discretization and solution method? Mathematically the PDE doesn't have a discretization, that's just the numerical method, so it should be possible to setup an interface where a user writes down a PDE, chooses a method, and then a solver method runs. There's two possible levels this can be written at. First of all, it can be written at a level like The second option is to use a DSL like ModelingToolkit.jl to abstractly specify the PDEs. It can define So I think the intermediate solution is to do both. Setup the From there a nice ecosystem could evolve. More I think we will want to separate out the PDE docs though since I imagine this will grow quickly. I see it sectioning out like the DiffEq docs: each |
|
Integral constraints like Are there actually cases where people have non-linear constraints that are reducible to explicit solvers in an adhoc way? |
|
I don't think so, but why would be looking at an explicit solver for a problem with nonlinear constraints? Seems like it would be natural to use a discretization to a DAE. |
|
For nonlinear conservation laws, useful boundary conditions are often nonlinear. Nevertheless, they can be solved using spatial semidiscretisations such as finite volume or finite difference methods and explicit Runge-Kutta schemes. |
|
In economics there are lots of nonlinear problems (some of which you saw in SciML/PDERoadmap#22 post) but... I am not sure that this is possible to solve in general, and it is usually the interior themselves that are nonlinear rather than the boundaries. Many of the approaches to those use an iterative apporach where you define a linear differential equation, solve it, change the parameters, repeat. Of course, you could also define these as generating nonlinear systems of equations instead... but I think we should wait on that sort of generality. When Capstan.jl is mature and reverse-mode works well, the "bigass nonlinear system" approach might work better. So, my feeling is that this library should stay focused on affine and linear operators and boundaries, and that the fully nonlinear stuff should wait until this works. As an aside note, one of the first nonlinear problems I need to solve are differential quasi-variational inequalities, that is, writing down the stationary problem to solve for This isn't a linear problem, but it is a linear-complementarity problem when discretized. That is |
|
This issue is on DifferentialEquations.jl, not DiffEqOperators.jl. There was never a plan to limit it to affine boundary conditions. |
|
My mistake, I was confused. I don't have much to add on a general interface, as problems in economics tend to be of a very particular set of sub-types. What I will say, though, from my limited knowledge is that I think it may be too ambitious to try to come up with design for arbitrary PDEs at this point, since the methods for solving them are so radically different. Having a shared set of types/interfaces/etc. for common things (e.g. domains, grids, spaces, etc.) is a different question. If ApproxFun, DiffEqOperators, etc. share those things, then it is a step forward. |
Solving and specification are two completely different things. We don't have to commit to solve every representable PDE, just have a language to specify what problem to solve. What I already showed in the slides can do that if you allow the AbstractDomain to be something that can be package defined as well, so it's not very ambitious. |
|
Supersceded by being partially done, and the higher level part is #469 |
DiffEq's first PDE attempt (which started the package) was a failure. The issue was that it was tied directly to an FEM subpackage I made internally, but then it couldn't scale to all of the problems it needed to face. The rest of DiffEq was then built up around decentralization. The core idea is that if you separate the data that specifies a problem and leave the solver support to dispatch, many different approaches and packages can step in to offer unique (and fast) ways to solve the problems. That has worked splendidly, so the FEM support was dropped in DiffEq 4.0 since it didn't conform to this.
What I want to return with is a PDE story worthy of the rest of DiffEq. So the goals it needs to meet is:
Here is a prototype for a design that could hopefully handle this.
General Flow: Heat Equation FDM
Let's talk about the Heat Equation. First we start with a domain type. Domain types are specified by types
<:AbstractDomain. For example,is a regular grid of
[0,5]with spacingdx. But now we have to start talking about boundaries. It's the domain that knows what boundaries it has, and thus what has to be specified. So to the domain we need to add boundary functions.Boundary conditions are specified by types
<:AbstractBoundaryCondition. So we could have:Now we'd package that together:
Those domains are compatible with those boundary conditions, so it works without error. Once we have a space, our Heat Equation is specified via:
so we specify:
This is just a bundle of information that describes what we want to solve. But how do we want to solve it? We need to do a call that converts it into something with usable operators. Let's do it method of lines via a conversion call:
where
FiniteDifferencedispatches this conversion to DiffEqOperators.jl with default order.FiniteDifference(order=4)can set the spatial discretization order, or we can choose other packages with adaptive space operators etc. This we can take to the ODE solver and solve.Details
Now let's talk about details.
Domains
Domains are a specification of the domain to solve on. The basic ones will be provided by DiffEqPDEBase.jl. Wrapper libraries such as DiffEqApproxFun.jl or FEniCS.jl could provide extras, like
SpectralGridorFEMMesh. It can be as lean as possible, though in some cases like FEM it will need to hold node/mesh pairings.Boundary Conditions
These are defined in DiffEqPDEBase.jl and have the minimal information describing a specific kind of boundary condition. Either holding constants, arrays, or functions, they just encapsulate the idea of the BC choice.
Space
Space is how domains and BCs come together. For a given domain, it will make sure the BCs are compatible with the domain. Can this be merged with Domains?
Operators
This part wasn't shown in the simple example. One we have a space, that is enough information to define operators. The high level interface is something like
DerivativeOperator(space,alg,derivative_order,approximation_order)etc, and then via
algit dispatches to the appropriate package to spit outDiffEqOperators (this stuff is already defined) appropriate for that space. So like DiffEqOperators would just need to be changed to implement that dispatch.Problems
Problems are just an abstraction over operators which basically knows what operators it should be building given the space, alg, and problem. They don't do anything until the conversion call, in which case they build these operators to then build something to solve.
Conversions
The obvious conversions are to ODE/SDE/DAE/DDEProblems and keep on going. We need two more problem types though. LinearProblem and NonlinearProblem. Then we cover it all. Everything becomes instantiates operators through some conversion.
High-Level Conversions
Then there can also be some standard conversions like
DiffusionAdvection(prob,space,tspan)whereprob=SDEProblem. This is just high level sugar to make everything boil down to conversions on PDE problems to then give linear, nonlinear, or diffeqs out in the end.Big Questions
Documentation
I've learned by now to plan documentation levels at the same time as the interface. Where are things documented? I think the general workflow is documented, and then boundary conditions are documented generally, some general domains, the operator calls, and problems are all documented package-wide generic. Then unique domains, space, operator algorithms, and conversions are documented per-package. DiffEqOperators describes no unique domains, it describes what BCs it can support, choices for operator construction, and problem conversions it can perform. DiffEqApproxFun links over to ApproxFun for choices of domains, has a single
Spectral()for generating the spectral operators (with choices in there), and then can generate ODEs for the Heat Equation, or a fully implicit discretization to a LinearProblem via a 2D discretization (can it?).Representation Information?
When you get this ODEProblem... what does it actually represent? We need a way to describe back to the user what
u0means. Coefficients of ...?u0[i,j,k]means what(x,y,z)? Not sure how to handle this.Time?
Time is kept separate because it usually acts differently when it exists. For algorithms like a fully implicit discretization,
LinearProblem(prob,FullyFiniteDifference(dt))or something.Passing space
This is just minor bikeshedding, but instead of like making the Heat Equation require the user give
f(u,p,t,x,y,z)for 3D, should it bef(u,p,t,x)forxa tuple/number? That keeps it the same regardless of dimension. Minor bikeshedding.Conclusion
This allows us to specify our problem at a high level, interface with different packages in an extendable way to generate appropriate FDM/FVM/FEM/spectral discretizations, and spit out linear/nonlinear/diffeq problems to take to solvers. The hope is that this give enough places to inject information and package-specific components that it can be efficient, yet gives a good level of automation and transferrability to other setups.
Thanks. Pulling in some people to bikeshed this high level story. @dextorious @dlfivefifty @alanedelman @jlperla
The text was updated successfully, but these errors were encountered: