Namespaces
namespace	tests

Functions
void	applyDirichletBoundaryConditions (std::function< void(const tarch::la::Vector< Dimensions, double > &x, double t, double Q)> boundaryCondition, const peano4::datamanagement::FaceMarker &marker, const tarch::la::Vector< Dimensions, double > &domainOffset, const tarch::la::Vector< Dimensions, double > &domainSize, const std::bitset< Dimensions > &periodicBC, int numberOfGridCellsPerAxis, int overlap, int unknownsPlusAuxiliaryVariables, double t, double faceData)
	Set the ghost cells of a face along the domain boundary through a Dirichlet boundary condition.

void	applySommerfeldBoundaryConditions (std::function< double(const double Q, const tarch::la::Vector< Dimensions, double > &x, int normal)> waveSpeed, std::function< void(const tarch::la::Vector< Dimensions, double > &x, double Q)> farFieldSolution, const peano4::datamanagement::FaceMarker &marker, const tarch::la::Vector< Dimensions, double > &domainOffset, const tarch::la::Vector< Dimensions, double > &domainSize, const std::bitset< Dimensions > &periodicBC, int numberOfGridCellsPerAxis, int overlap, int numberOfUnknowns, int numberOfAuxiliaryVariables, const tarch::la::Vector< Dimensions, double > &origin, double timeStepSize, const double faceDataOld, double faceDataNew)
	Advance the ghost cells of a face along the domain boundary through a Sommerfeld (radiative, outgoing) boundary condition.

Function Documentation

◆ applyDirichletBoundaryConditions()

void exahype2::blockstructured::applyDirichletBoundaryConditions	(	std::function< void(const tarch::la::Vector< Dimensions, double > &x, double t, double *Q)>	boundaryCondition,
		const peano4::datamanagement::FaceMarker &	marker,
		const tarch::la::Vector< Dimensions, double > &	domainOffset,
		const tarch::la::Vector< Dimensions, double > &	domainSize,
		const std::bitset< Dimensions > &	periodicBC,
		int	numberOfGridCellsPerAxis,
		int	overlap,
		int	unknownsPlusAuxiliaryVariables,
		double	t,
		double *	faceData
	)

Set the ghost cells of a face along the domain boundary through a Dirichlet boundary condition.

A face of a block-structured patch carries an overlap of 2*overlap grid cells (Finite Volumes or Finite Difference grid cells) along its normal. The left half belongs to the cell left of the face, the right half to the cell right of it. On the domain boundary, one of these two cells lies outside of the domain, so nobody ever writes its half of the face data. This routine fills exactly this half: for each of the outside's grid cells, it invokes

boundaryCondition( x, t, Q );

with x being the centre of the (ghost) grid cell, t the time stamp handed in, and Q pointing to the grid cell's unknowns and auxiliary variables within faceData. Layer 0 is the ghost grid cell adjacent to the face.

Whether the face is on the domain boundary is derived from its position, i.e. compared against domainOffset and domainOffset+domainSize. Faces along periodic axes and faces which are not on the domain boundary are left unchanged. Hence, the routine can be called for any face.

The routine replaces the code which the action sets exahype2.api.fv.actionsets.ApplyDirichletBoundaryConditions and exahype2.api.rkfd4.actionsets.ApplyDirichletBoundaryConditions used to generate. Both generate a call to it now. The geometry of a Finite Volume and of a Finite Difference grid cell is the same, so the two schemes share this routine.

Parameters

boundaryCondition	Functor setting the unknowns and auxiliary variables Q at position x and time t. Typically a lambda which forwards to the solver's dirichletBoundaryConditions().
marker	Marker of the face, i.e. the face's centre, the size of the adjacent cells and the number of the face.
domainOffset	Lower left corner of the computational domain.
domainSize	Size of the computational domain.
periodicBC	Periodic boundary conditions per axis. Faces along periodic axes are left unchanged.
numberOfGridCellsPerAxis	Number of grid cells per axis within one patch.
overlap	Number of grid cell layers per side of the face, i.e. the halo size.
unknownsPlusAuxiliaryVariables	Number of doubles per grid cell.
t	Time stamp handed over to boundaryCondition.
faceData	The face's overlap data. Its layout is the one of toolbox::blockstructured::serialiseVoxelIndexInOverlap().

◆ applySommerfeldBoundaryConditions()

void exahype2::blockstructured::applySommerfeldBoundaryConditions	(	std::function< double(const double *Q, const tarch::la::Vector< Dimensions, double > &x, int normal)>	waveSpeed,
		std::function< void(const tarch::la::Vector< Dimensions, double > &x, double *Q)>	farFieldSolution,
		const peano4::datamanagement::FaceMarker &	marker,
		const tarch::la::Vector< Dimensions, double > &	domainOffset,
		const tarch::la::Vector< Dimensions, double > &	domainSize,
		const std::bitset< Dimensions > &	periodicBC,
		int	numberOfGridCellsPerAxis,
		int	overlap,
		int	numberOfUnknowns,
		int	numberOfAuxiliaryVariables,
		const tarch::la::Vector< Dimensions, double > &	origin,
		double	timeStepSize,
		const double *	faceDataOld,
		double *	faceDataNew
	)

Advance the ghost cells of a face along the domain boundary through a Sommerfeld (radiative, outgoing) boundary condition.

What Sommerfeld conditions are

Many simulations cut a finite box out of an infinite domain. The classic example is a black hole spacetime: The solution is defined everywhere, but we can only afford to simulate \( [-L,L]^d \). Waves that are emitted in the interior should leave this box as if the box was not there. Dirichlet conditions, which hold the boundary at a prescribed value, do not achieve this: Everything that travels outward is reflected back into the domain. It is like a wave hitting a wall.

So instead, we apply Sommerfeld conditions. To be precise, we apply outgoing Sommerfeld conditions: Normally, these conditions handle incoming waves (there should be none) different to outgoing waves (they should not reflect). The calculations here always assume that the waves travel outwards, i.e. the eigenvalues along the normal are positive.

We observe that we cannot use plain homogeneous Neumann conditions either: If we set the voxels outside of the domain (in the halo) to the same value as their interior counterpart, they adopt to the solution too quickly. A wave front approaches the boundary and the voxel outside jumps immediately to the wave value in the domain. What it has to do instead, it has to lag a little bit behind. Just enough to mimick the behaviour as if the wave travelled outside since the previous time stamp.

More mathematical: If we derive the solution along the radius r, so look how it changes along the outgoing normal, and then derive by time, so see how things change as we walk out of the domain over time, then the result should be exactly as if we took our first-order PDE along the normal, which also describes how the solution changes over time. We assume that this outgoing beam describes a characteristic.

Now let's assume that we are far away and the PDE does not really do anything anymore. It only propagates the solution similar to a pure advection. It totally has deteriorated. In this case, the values outside of the domain (our ghost values) change depending on the the distance from their inner counterpart and the signal speed given by the maximum eigenvalue.

Let's make one further assumption: They assume that it is totally sufficient to know the maximum eigenvalue. For systems of PDEs, we would have to know the eigenvalue per equation, but we just use the maximum eigenvalue as stand-in approximation, ignoring that we look along the radius r and not along the Cartesian axes.

With all the considerations, we assume that, close to the boundary, every unknown \( u \) behaves like an outgoing spherical wave which decays towards a far-field value \( u_\infty \):

\begin{eqnarray*} u(r,t) = u_\infty + \frac{f(r - \lambda _{\text{max}} t)}{r}. \end{eqnarray*}

Here, \( \lambda _{\text{max}} \) is the speed at which the field propagates, \( f \) an arbitrary profile, and \( r = |x-x_0| \) the distance from the centre \( x_0 \) (origin) of the system, i.e. the point that emits the waves. As the equation depends on r only, the ansatz is radially symmetric. The argument in \( f \) models how quickly a signal propagates through the boundary. It almost serves as a countdown.

Differentiating the ansatz with respect to \( t \) and \( r \) eliminates the unknown profile \( f \), and we obtain a partial differential equation which the solution has to satisfy along the boundary:

\begin{eqnarray*} \partial_t u + \lambda _{\text{max}} \partial_r u + \frac{ \lambda _{\text{max}} }{r} ( u - u_\infty ) = 0 \qquad \text{with} \qquad \partial_r u = \frac{x_i - x_{0,i}}{r} \partial_i u. \end{eqnarray*}

The condition is exact for a spherical wave with constant speed. Everything else (non-spherical waves, waves that do not originate from \( x_0 \), several speeds) is reflected partially. The reflections become smaller the further the boundary is away from the sources. Different to a Dirichlet condition, the boundary value is not known a priori: It follows from the solution inside the domain, and it has its own time evolution. A Sommerfeld boundary condition hence is an evolution equation for the ghost cells.

Discretisation

This routine works on the grid cells (Finite Volumes or Finite Difference points) of a block-structured patch. Along a boundary face with normal \( n \), there are overlap layers of ghost cells outside of the domain. We only discretise the derivative along the normal, as a ghost cell has no neighbours outside of the domain along the tangential directions which we could use. For a field which depends on \( r \) only, \( \partial_n u = \frac{x_n - x_{0,n}}{r} \partial_r u \). We invert this relation, i.e. we reconstruct the radial derivative from the normal one:

\begin{eqnarray*} \partial_r u \approx \frac{r}{|x_n - x_{0,n}|} \partial_n u. \end{eqnarray*}

Along the coordinate axes, \( \partial_r u = \partial_n u \). In the edges and corners of the domain, the factor grows up to \( \sqrt{d} \). The alternative, i.e. to take the definition of \( \partial_r u \) above and to drop the tangential derivatives, would yield the factor \( |x_n - x_{0,n}|/r \leq 1 \) instead, and would underestimate the radial derivative of a spherical wave.

Let \( o \) be a ghost cell and \( i \) its neighbour one step further inside, i.e. the first cell inside the domain for ghost layer 0, and ghost layer \( l-1 \) for ghost layer \( l \). We discretise the condition with a box scheme: It is centred between \( o \) and \( i \) in space, and between \( t^n \) and \( t^{n+1} = t^n + \Delta t \) in time:

\begin{eqnarray*} \frac{o^{n+1}-o^n+i^{n+1}-i^n}{\Delta t} + \sigma \left( o^{n+1}-i^{n+1}+o^n-i^n \right) + \frac{ \lambda _{\text{max}} }{r} \left( o^n+i^n-2u_\infty \right) = 0 \qquad \text{with} \qquad \sigma = \frac{ \lambda _{\text{max}} }{h} \frac{r}{|x_n - x_{0,n}|}. \end{eqnarray*}

\( h \) is the grid cell size along the normal. The damping term is evaluated explicitly at the old time. Solved for the new ghost value, this yields

\begin{eqnarray*} o^{n+1} = \frac{ (1 - \sigma \Delta t) o^n + (1 + \sigma \Delta t) i^n - (1 - \sigma \Delta t) i^{n+1} - \Delta t \frac{ \lambda _{\text{max}} }{r} \left( o^n+i^n-2u_\infty \right) }{ 1 + \sigma \Delta t }. \end{eqnarray*}

The routine evaluates \( \lambda _{\text{max}} \) via waveSpeed() for \( o^n \), \( u_\infty \) via farFieldSolution(), and \( r \) and \( x_n \) at the centre of the ghost cell. It runs through the ghost layers from the face outward. Therefore, \( i^{n+1} \) is the new interior solution for layer 0, and the ghost value which the routine has just computed for all further layers. The scheme is consistent for any overlap, i.e. it serves the one ghost layer of a first-order Finite Volume scheme as well as the three layers of fourth-order Finite Differences.

A state which equals \( u_\infty \) everywhere is preserved exactly. If \( \Delta t = 0 \), the routine copies the old ghost values, i.e. it sets \( o^{n+1} = o^n \). This is not what the formula above yields in the limit: there, the jump of the interior solution would be passed on to the ghost cell. See the discussion of Runge-Kutta schemes below for why zero time step sizes arise.

The auxiliary variables are not evolved. Each ghost cell takes the auxiliary variables of \( i^{n+1} \), i.e. they are extrapolated from the new interior solution.

How to use the routine within a block-structured scheme

Block-structured schemes such as exahype2.api.fv and exahype2.api.rkfd4 hold the halo data of their patches on the faces: A face carries 2*overlap layers of grid cells along its normal. The left half belongs to the cell left of the face, the right half to the cell right of it. On the domain boundary, one of the two halves lies outside of the domain. These are the ghost cells, and nobody but the boundary treatment writes them. The layout of the face data is the one of toolbox::blockstructured::serialiseVoxelIndexInOverlap().

Different to exahype2::blockstructured::applyDirichletBoundaryConditions(), which only needs a position and a time stamp, the Sommerfeld update needs two versions of the face data:

faceDataOld at \( t^n \): both halves have to be valid, i.e. the interior half holds \( i^n \) and the ghost half \( o^n \).
faceDataNew at \( t^{n+1} \): the interior half holds the new interior solution. The routine writes the ghost half.

The routine is part of the roll over

Block-structured schemes project their patches onto the faces into a temporary buffer, and then copy (roll over) this buffer into the persistent face data, which the reconstruction of the halo reads in the next sweep. Directly before the roll over, the persistent face data still holds the old solution, and the temporary buffer already holds the new interior solution. This is exactly what the routine needs, so it does not need any additional storage.

Therefore, the Sommerfeld update replaces the plain roll over of a face: There is no separate boundary task in the scheme's task graph, which sets the ghost cells before the reconstruction. Instead, the roll over first calls the routine with faceDataOld being the persistent face data and faceDataNew the temporary buffer, and then copies both halves of the face. After the roll over, the persistent ghost half holds \( o^{n+1} \), i.e. it is ready for the next time step or Runge-Kutta stage. A typical roll over within the generated touchFaceLastTime() reads

::exahype2::blockstructured::applySommerfeldBoundaryConditions(
  [](
    const double*                                Q,
    const tarch::la::Vector<Dimensions, double>& x,
    int                                          normal
  ) -> double {
    double lambdaMax;
    repositories::instanceOfMySolver.maxEigenvalue(Q, x, normal, &lambdaMax);
    return lambdaMax;
  },
  [](const tarch::la::Vector<Dimensions, double>& x, double* Q) -> void {
    repositories::instanceOfMySolver.sommerfeldFarFieldSolution(x, Q);
  },
  marker,
  DomainOffset,
  DomainSize,
  PeriodicBC,
  NumberOfGridCellsPerAxis,
  Overlap,
  NumberOfUnknowns,
  NumberOfAuxiliaryVariables,
  {0.0, 0.0, 0.0},
  repositories::instanceOfMySolver.getMinTimeStepSize(),
  fineGridFaceQCurrentTimeStep.value,
  fineGridFaceQUpdate.value
);
// copy both halves of fineGridFaceQUpdate into fineGridFaceQCurrentTimeStep

The roll over runs in touchFaceLastTime() of the sweep that projects the new solution. The solver's time stamp and time step size are still the ones of the time step which has just been computed, as the solver only advances them in unprepareTraversal(). The roll over does not need to know whether a face is on the domain boundary: It may call the routine for all faces. The routine leaves faces which are not on the domain boundary or which lie along a periodic axis unchanged, and it never writes into faceDataOld.

Within the Python API, users add an exahype2.api.SommerfeldBoundaryConditions term to their solver, and the schemes of exahype2.api.fv and exahype2.api.rkfd4 then replace their roll over with RollOverUpdatedFaceWithSommerfeldBoundaryConditions. The schemes do not add a separate boundary task to their time step in this case (exahype2.api.blockstructured.BoundaryConditions).

A scheme has to take care of four further points:

Initial data. The ghost half of the persistent face data has to hold valid data before the first time step, as it is the \( o^n \) of the first update. A natural choice is to evaluate the initial condition at the centres of the ghost cells within the roll over of the grid construction. This is exactly what exahype2::blockstructured::applyDirichletBoundaryConditions() does if it is handed the initial condition as boundary condition.
Single-step schemes such as explicit Euler in time call the routine once per time step with the time step size \( \Delta t \).
Runge-Kutta schemes project the input \( Y_k \) of each stage onto the faces. The ghost cells of stage \( k \) should represent the boundary at the stage's time stamp \( t^n + c_k \Delta t \). Therefore, the roll over after each stage advances the ghost cells from one stage's time stamp to the next one, i.e. by \( (c_{k+1}-c_k) \Delta t \), and by \( (1-c_{s-1}) \Delta t \) after the last stage. Over a time step, these increments add up to \( \Delta t \). Stages with the same time stamp (classic RK4 has \( c = 0, \frac{1}{2}, \frac{1}{2}, 1 \)) yield an increment of zero. The stage inputs then differ, but they represent the same time, so the ghost cells remain unchanged. This is the reason why the routine copies the ghost cells for \( \Delta t = 0 \).
The old face data has to survive until the roll over. A scheme with Dirichlet conditions may discard the persistent face data once it has reconstructed the halo, and does not have to load it again before the projection, as the roll over overwrites it completely. With Sommerfeld conditions, the roll over reads it as faceDataOld. Within a task graph, every task of the time step which uses the persistent face data hence has to declare it as READ_AND_WRITE (exahype2.api.blockstructured.BoundaryConditions.get_access_to_persistent_face_data()). Otherwise, the generated code neither stores the face data after the reconstruction nor loads it before the projection, and the routine reads garbage.

Constraint: the stage time stamps must not decrease

As the roll over chains the ghost cells from one stage to the next, a Runge-Kutta scheme may only use the routine if its stage time stamps never decrease, i.e. if

\begin{eqnarray*} 0 = c_0 \leq c_1 \leq \dots \leq c_{s-1} \leq 1. \end{eqnarray*}

A scheme which uses Sommerfeld conditions has to check this constraint for its Butcher tableau and has to refuse tableaus which violate it. exahype2.api.rkfd4.CentralDifferencesWithRungeKutta raises a ValueError in this case. All tableaus of exahype2.api.rungekutta.ButcherTableau fulfil it: Their stage time stamps are 0 (explicit Euler), \( 0, \frac{1}{2} \) (explicit midpoint rule), \( 0, \frac{1}{2}, 1 \) (Kutta's third-order scheme), and \( 0, \frac{1}{3}, \frac{2}{3}, 1 \) (Kutta's 3/8 rule).

The order conditions of Runge-Kutta schemes do not require non-decreasing time stamps, and several popular schemes violate the constraint. The strong-stability-preserving third-order scheme of Shu and Osher has \( c = 0, 1, \frac{1}{2} \), the Runge-Kutta-Fehlberg scheme of order 4(5) has \( c = 0, \frac{1}{4}, \frac{3}{8}, \frac{12}{13}, 1, \frac{1}{2} \). For such a scheme, an increment \( c_{k+1}-c_k \) is negative. The increments still add up to \( \Delta t \), but the roll over then integrates the damped outgoing wave equation backward in time: Instead of damping the ghost cells, it amplifies them, and the denominator \( 1 + \sigma \Delta t \) of the update shrinks.

Schemes with decreasing stage time stamps have to start the update of every stage from the beginning of the time step instead, i.e. advance the ghost cells of stage \( k \) from \( t^n \) by \( c_k \Delta t \). This requires the ghost cells at \( t^n \) throughout all stages of a time step, i.e. an additional persistent buffer per face which the roll over does not overwrite.

Choice of the wave speed and of the far field

The wave speed is the speed at which the outgoing waves travel. The routine evaluates it per ghost cell and uses one speed for all unknowns. The maximum eigenvalue of the PDE along the normal is a safe choice: It is the speed of the fastest wave. Some codes use a known physical speed instead, such as the speed of light for most fields of a spacetime. Both choices are up to the waveSpeed() functor.

The far-field solution \( u_\infty \) is the value the solution approaches at infinity, e.g. the flat spacetime for a black hole. It may depend on the position, but not on time.

Limitations

The scheme only uses derivatives along the normal. Ghost cells in the edges and corners of the domain are never set, as the faces do not hold them.
The condition assumes spherical waves around origin. All other waves are partially reflected.
The time stepping of the ghost cells is first order for Runge-Kutta stages of higher order, as the box scheme only sees the interior solution at the stages' time stamps.
Runge-Kutta schemes may only use the routine if their stage time stamps do not decrease (see the constraint above).
The routine works on one resolution level. It is up to the scheme to provide valid data on all boundary faces it reads from.

Differences to the legacy implementation

The routine replaces the legacy exahype2::fd::applySommerfeldConditions() (removed), which the legacy Python API (python/exahype2) injected through the TemplateHandleBoundary_KernelCalls attribute of its solvers. The legacy routine realised the same box scheme, but differed in a few details:

It used the signed coordinate \( x_n \) within \( \sigma \). On lower boundary faces of a domain around the origin, \( \sigma \) hence was negative, i.e. the condition modelled an incoming wave. The update of the ghost cell then has an amplification factor of \( (1+|\sigma|\Delta t)/(1-|\sigma|\Delta t) > 1 \). We use the distance \( |x_n - x_{0,n}| \) instead.
It measured \( r \) relative to the origin, but took \( x_n \) from the absolute coordinate.
It evaluated the wave speed at the face centre rather than at the ghost cell.
It copied the inner data into the ghost cells of the old face data if the time stamp was zero or equalled a checkpoint time stamp, as the ghost cells had not been initialised. Here, it is the scheme's responsibility to initialise them.
The auxiliary variables were taken from \( i^n \).