Add a FDM partial operator with Dirichlet boundary conditions

[protaisM]

We may need to have a centred FDM operator of order two which can handle boundary conditions, starting with Dirichlet ones

[EmilyBourne]

Are you sure you mean "Dirichlet boundary conditions"? (i.e. the value of the derivative at the boundary is known).

When we discussed this I thought you were talking about adding an extrapolation condition, e.g. ConstantExtrapolation to allow the derivative to be calculated from approximations of the function outside the domain?

[protaisM]

I think I mean Dirichlet boundary conditions, i.e. the value of the function at the boundary is known. When the value of the derivative is known we call this Neumann boundary condition (see BC examples).

However, numerically the Dirichlet BC is often taken to be outside of the domain instead of in the last mesh point. This is what is called in gyselalib extrapolation rules, which are very specific to splines because they are used to interpolate. For FDM (of order 2) we only need one point outside of the domain of derivation, whose value we call Dirichlet BC. The boundary condition are however not part of the FDM operator but of the function. For instance :

We have a physics problem where the heat equation has to be solved on a grid $x_1, ..., x_n$, simulating a rod with fixed temperature $T$ at its extremities. To solve this we need to compute the laplacian of the temperature. Numerically this is given by $\Delta T_i = (T_{i-1}-2T_i+T_{i+1})/h^2$
but what do we do at the extremities ? This is given by the Dirichlet BC, $\Delta T_n = (T_{n-1}-2T_n+T)/ h^2$
The discretized laplacian does not contain the boundary condition. It is always the same, for all functions, but take as input the Dirichlet BC it needs to compute the derivative according to the physics assumption.
I hope it is clear, and that the latex will be interpreted...

[EmilyBourne]

I think the fact that we are talking about an FDM confuses this somewhat.

It is not clear to me if $f$ is the function being differentiated or if it is the result of the operation. Wikipedia thinks it is the latter:

However, numerically the Dirichlet BC is often taken to be outside of the domain instead of in the last mesh point.

We took care when naming boundary conditions/extrapolation conditions to avoid confusing these two. A boundary condition is a condition imposed at the boundary (and nowhere else). I think it is always a closure condition for an operator. An extrapolation condition is an approximation of a function outside the simulation domain. So a boundary condition is something like $f(x_0) = 0$, while an extrapolation condition is something like $f(x) = f(x_0) \forall x<x_0$. It completes the definition of a function.
What you are describing:

one point outside of the domain of derivation

not part of the FDM operator but of the function

certainly sounds like an extrapolation condition and not a boundary condition.

[EmilyBourne]

Suppose we want to find $f(x) = \partial F(x)$ on the grid $x_0,...,x_n$.

For me a Dirichlet boundary condition on the solution would be:
$f(x_0) = \alpha$
$f(x_n) = \beta$
with $\alpha$ and $\beta$ known and not calculated from $F(x)$ as part of the FDM operator.

An extrapolation operator would be $F(x) = \alpha \forall x<x_0$ or $F(x) = F(x_0) \forall x<x_0$, then the FDM operator which is numerically defined as:
$\Delta F_i = (F_{i-1} - 2 F_i + F_{i+1}) / h^2$ with $F_i = F(x_i)$ can be defined at $i=0$:
$$\Delta F_0 = (F_{-1} - 2 F_0 + F_{1}) / h^2 = (F_0 - 2 F_0 + F_{1}) / h^2$$

The discretised laplacian does not contain the boundary condition, it works in the same way at all points

[protaisM]

Okay, let's take a step back. I have a parabolic PDE (technically it is not uniformly parabolic, because the coefficients are evolving in time, but it does not change my point) of the form
$\partial_t n - \partial_x \Gamma_n = f(t,x,n)$
$\Gamma_n = g(t,x)\partial_x n - D \Delta n$
where $f$ and $g$ are functions we do not care too much about (actually they depend on other quantities that are not relevant here), and $D$ is a time-dependent coefficient.
What is the solution of this PDE? It is a couple $(n, \Gamma_n)$ depending on time. Fine, but this is not well determined. We need to know what happens at the boundary (if for example we want to solve it on $[0,1]$). Physically, we have no condition on $n$: no fixed value, no fixed derivative of any order, nothing. However, we need to have a boundary condition on $\Gamma_n$. Let's take the easiest case: for $n$ to be a conserved quantity (its integral over all $x$ to be constant in time) we need to have $\Gamma_n$ equals $0$ at the boundary. This is what we call a Dirichlet boundary condition (very precisely a Dirichlet boundary condition on $\Gamma_n$). With this (and an initial value for $n$ and $\Gamma_n$), our problem goes from a parabolic PDE to a Dirichlet problem for a parabolic PDE.
Now we need to discretise this Dirichlet problem to solve it numerically. Where will this Dirichlet boundary condition come into play? The only place it will appear is in the evaluation of $\partial_x \Gamma_n$.

Little mathematical parenthesis here to explain why the boundary condition will not come into play in the second equation, namely
$\Gamma_n = g(t,x)\partial_x n - D \Delta n$.
This is because PDE, in mathematics, are solved in open domains, i.e. if we want to solve a Dirichlet problem on $[0,1]$ we actually solve the PDE on $(0,1)$ and enforce the boundary condition on $0$ and $1$. Thus, theoretically we should not need our last mesh points to be at the boundary of the domain. However, it makes no sense to remove the first and last points of the dicretised domain, so instead for a discretisation $x_1, ..., x_n$ we just consider the domain of the PDE to be $[x_1-\epsilon, x_n+\epsilon]$, which is just simpler numerically, and does not change anything to the result.

With this, I think I have everything for my problem. I want to implement a "partial derivative" operator that I will use to evaluate $\partial_x \Gamma_n$. I should give it my $\Gamma_n$ and the boundary condition needed to evaluate the derivative on the last mesh point. Now back to our questions on vocabulary.

It is not clear to me if f is the function being differentiated or if it is the result of the operation.

Here is the tricky point. In the case of this operator, $f$ (so $\Gamma_n$) is the function being differentiated and is also the solution of the PDE, on which we have boundary conditions. These two are completely uncorrelated.

We took care when naming boundary conditions/extrapolation conditions to avoid confusing these two.

The only time in the code I am seeing boundary conditions are for splines, and they are used to know what should the spline looks like near the boundary. The only time I'm seeing extrapolation conditions is to precise what we should take for the value of an interpolation when we are outside of the interpolation domain. And when I search in the code base for boundary condition I very often see things like that:
ddc::ConstantExtrapolationRule<R, Theta> boundary_condition_r_left
so I thing these two are very confused in the code.

So a boundary condition is something like $f(x_0) = 0$ , while an extrapolation condition is something like $f(x) = f (x_0) \forall x < x_0$ .

Just for you to know, mathematically speaking the first one is called a Dirichlet boundary condition and the second one is called an Neumann boundary condition.

To conclude (otherwise I can continue to write indefinitely...) I think the main problem is that we try to use the same code object for different mathematical objects.

• an ExtrapolationRule is defined in the context of interpolation: the interpolation defines how we recover the (approximated) value of a function inside of a domain and the extrapolation defines what we should return when we are outside of the domain
• a boundary condition is defined in the context of PDEs, and is very specific to the PDE we are studying. It can be of several types and applied to various quantities. It seems hard to come up with a general class implementing "boundary conditions" for all types of PDEs
• in the context of partial operators we need to give a rule on "how to compute the derivative at the boundary". This is given by the boundary condition of the Dirichlet problem (in our case) and serves the same purpose as a ExtrapolationRule (namely: what do I do when I am outside of the interior of the domain). So asking whether it is a boundary condition or an ExtrapolationRule we should implement is maybe not the right question to ask.

What I want to implement is the possibility to call the operator partial_x with a DFieldXY (for instance) as f, and two "boundary values" i.e two DFieldY bdy_value_left and bdy_value_right that will act as value taken by f when we want to evaluate the derivative at the boundary. Maybe calling them "boundary values" will be clearer for everybody.

I hope this is clear enough ! And thank's for the reading... 😄

[Philipp137]

I understand Matthieu's point. In the future, we may use approaches beyond just semi-Lagrangian methods, so we need an infrastructure that can pass boundary conditions (BCs) in any form to the discretization scheme—whether it's FEM (Finite Element Method), FD (Finite Difference), or BSL (backward Semi-Lagrangian).
Having a centralized infrastructure for this would be beneficial.

From my experience in other projects, boundary information is often already encoded in the meshing file. For example:

In 2D, an edge between two nodes may need to specify:

• Whether it’s a boundary edge. (or node)
• What type of condition applies (e.g., edge_type="inflow_static_potential").

This doesn’t directly translate to standard BCs (like Robin or Dirichlet), because the same edge or node could enforce different conditions for different quantities. For instance:

"inflow_static_potential" could mean:

• A time-dependent Dirichlet BC for density.
• A Neumann BC for electric potential.

Thus, different solvers (for kinetic part or poisson part) may interpret the same edge/node differently. Additionally, each solver might handle BCs in its own way.

We might want to look at the codes that we want to couple in the future. Like the FEM code in garchingen, how do they tag the boundary conditions at different positions.
Also, we want to communicate this with the codes now written for gysela (Mesher).

[PaulineVidal]

I think the main problem is that we try to use the same code object for different mathematical objects.

I am not sure about it. Currently the ddc::BoundCond refers to the "closure conditions" used to build a spline representation. To evaluate the spline, other objects were defined to defined the extrapolation type (ddc::ConstantExtrapolationRule, ddc::NullExtrapolationRule, ...).
There is no tag for Dirichlet, Neumann, Robin boundary conditions. I think the operators (Poisson, ... ) were defined for a specific one which should be precised in the doc. 🤔

[alex-m-h]

Dirichlet boundary condition are the values of the solution of a PDE which are fixed on the boundary. I agree with Emily that if values outside the domain are used in a discrete operator, on ghost points for a FD stencil for example, they should not be called "boundary condition".

For the future naming, I have another suggestion. The boundary condition is a constraint. "Closure condition" would actually only be appropriate, when this constraint leads to a closed system of equations. There are situations where this is not the case, e.g. one can find approximate solution using some norm minimization techniques even when the system of equations is underdetermined. So if we want to stay general, calling them "constraint" would be accurate. If they are exclusively used for interpolation and always lead to a closed system of equations, just calling them "interpolation type" or something like that would be the easiest one to understand.