exagrype::ccz4 Namespace Reference

Namespaces
namespace	internal
namespace	tests

Functions
void	ncp (double *BgradQ, const double *const Q, const double *const gradQSerialised, const int normal, const int CCZ4LapseType, const double CCZ4ds, const double CCZ4c, const double CCZ4e, const double CCZ4f, const double CCZ4bs, const double CCZ4sk, const double CCZ4xi, const double CCZ4mu, const double CCZ4SO)
	Non-conservative product of the first-order CCZ4 formulation.
void	source (double *S, const double *const Q, const int CCZ4LapseType, const double CCZ4ds, const double CCZ4c, const double CCZ4e, const double CCZ4f, const double CCZ4bs, const double CCZ4sk, const double CCZ4xi, const double CCZ4itau, const double CCZ4eta, const double CCZ4k1, const double CCZ4k2, const double CCZ4k3, const double CCZ4SO)
	The source term is one out of two terms that we use in our CCZ4 formulation.
void	maxEigenvalue (const double *const Q, int normal, const double CCZ4e, const double CCZ4ds, const double CCZ4GLMc, const double CCZ4GLMd, double *maxEigenvalue)
void	admconstraints (double *constraints, const double *const Q, const double *const gradQSerialised)
	This is a postprocessing routine to monitor if the physical constraints are fulfilled.
void	ThetaOutputNCP (double *NCPterm, const double *const Q, const double *const gradQSerialised, int normal)
	A temporary test function, to output Hamilton constraint related term in theta, 1 terms: RPlusTwoNablaZNCP.
void	TestingOutput (double *terms, const double *const Q, const double *const gradQSerialised)
	A temporary test function, to output some testing values 0,1 entries: Hamilton constraint related term in theta, 2 terms: RPlusTwoNablaZSrc, pure Src 2-10 entry: no-symmetry version of R_ij.
void	Psi4Calc (double *Psi4, const double *const Q, const double *const gradQSerialised, double *coor)
	This function is for the calculation of psi4, a quantity related to gravitional wave.
void	enforceCCZ4constraints (double *__restrict__ newQ)
	A postprocessing routine which pushes the volume solution back into the area of the CCZ4 constraints.
void	addKreissOligerDissipation (const double *__restrict__ QWithHalo, double *__restrict__ rightHandSide, const tarch::la::Vector< 3, double > &patchSize, int numberOfGridCellsPerPatchPerAxis, int haloSize, int unknowns, int auxiliaryVariables, double sigma, const int *quantities, int numberOfQuantities)
	Add the Kreiss-Oliger dissipation of a patch to its right-hand side.
void	recomputeAuxiliaryVariablesFD4_4thOrder (::exahype2::CellData< double, double > &patchData, int numberOfGridCellsPerPatchPerAxis, int haloSize, int unknowns, int auxiliaryVariables)
	Recompute auxiliary variables for FD4 scheme with a 4th order scheme.
void	recomputeAuxiliaryVariablesFD4_centralDifferences (::exahype2::CellData< double, double > &patchData, int numberOfGridCellsPerPatchPerAxis, int haloSize, int unknowns, int auxiliaryVariables)
void	recomputeAuxiliaryVariablesFD4_leftDifferences (::exahype2::CellData< double, double > &patchData, int numberOfGridCellsPerPatchPerAxis, int haloSize, int unknowns, int auxiliaryVariables)
void	recomputeAuxiliaryVariablesFD4_rightDifferences (::exahype2::CellData< double, double > &patchData, int numberOfGridCellsPerPatchPerAxis, int haloSize, int unknowns, int auxiliaryVariables)

Function Documentation

addKreissOligerDissipation()

void exagrype::ccz4::addKreissOligerDissipation	(	const double *__restrict__	QWithHalo,
		double *__restrict__	rightHandSide,
		const tarch::la::Vector< 3, double > &	patchSize,
		int	numberOfGridCellsPerPatchPerAxis,
		int	haloSize,
		int	unknowns,
		int	auxiliaryVariables,
		double	sigma,
		const int *	quantities,
		int	numberOfQuantities
	)

Add the Kreiss-Oliger dissipation of a patch to its right-hand side.

The routine adds

\( \frac{\sigma}{64 \Delta x_d} (Q_{i-3} - 6 Q_{i-2} + 15 Q_{i-1} - 20 Q_i + 15 Q_{i+1} - 6 Q_{i+2} + Q_{i+3}) \)

per axis d to the right-hand side of each grid cell i of a patch, i.e. a sixth-order dissipation term which damps oscillations on the scale of the grid cells. \( \Delta x_d \) is the patch size along d divided by the number of grid cells per axis. The term is added to the selected quantities only, all other entries of the right-hand side remain unchanged.

The routine is the dissipation which the legacy fourth-order Finite Differences of ExaHyPE built into their kernels. These kernels have been removed. The Finite Differences of exahype2.api.rkfd4 have no dissipation. exagrype.api.FD4WithKODissipation invokes this routine after their right-hand side kernel, within each Runge-Kutta stage.

Both arrays are stored lexicographically in AoS format, i.e. all quantities of one grid cell are stored together, and the first axis runs fastest. The stencil reaches three grid cells into the halo, i.e. haloSize has to be at least three.

Parameters

QWithHalo	Input of the Runge-Kutta stage with its halo, i.e. (numberOfGridCellsPerPatchPerAxis+2*haloSize)^3 grid cells with unknowns+auxiliaryVariables entries each.
rightHandSide	Right-hand side of the stage without halo, i.e. numberOfGridCellsPerPatchPerAxis^3 grid cells with unknowns entries each. The routine adds the dissipation to it.
patchSize	Size of the patch, i.e. of the mesh cell.
quantities	Indices of the unknowns which are subject to the dissipation. Each index has to be smaller than unknowns.
numberOfQuantities	Number of entries of quantities.

admconstraints()

void exagrype::ccz4::admconstraints	(	double *	constraints,
		const double *const	Q,
		const double *const	gradQSerialised
	)

This is a postprocessing routine to monitor if the physical constraints are fulfilled.

Used in python script.

enforceCCZ4constraints()

void exagrype::ccz4::enforceCCZ4constraints

(

double *__restrict__

newQ

)

A postprocessing routine which pushes the volume solution back into the area of the CCZ4 constraints.

This is a local operation that we can invoke after each time step per volume.

If works on the full set of unknowns of the first-order formulation

maxEigenvalue()

void exagrype::ccz4::maxEigenvalue	(	const double *const	Q,
		int	normal,
		const double	CCZ4e,
		const double	CCZ4ds,
		const double	CCZ4GLMc,
		const double	CCZ4GLMd,
		double *	maxEigenvalue
	)

ncp()

void exagrype::ccz4::ncp	(	double *	BgradQ,
		const double *const	Q,
		const double *const	gradQSerialised,
		const int	normal,
		const int	CCZ4LapseType,
		const double	CCZ4ds,
		const double	CCZ4c,
		const double	CCZ4e,
		const double	CCZ4f,
		const double	CCZ4bs,
		const double	CCZ4sk,
		const double	CCZ4xi,
		const double	CCZ4mu,
		const double	CCZ4SO
	)

Non-conservative product of the first-order CCZ4 formulation.

ncp(), source() and maxEigenvalue() have one implementation in Kernels.cpp which serves both the CPU and OpenMP offloading. With GPUOffloadingOMP, they are device functions (omp declare target), declared here and defined in Kernels.cpp, i.e. they are not inlined into the GPU kernels.

Older LLVM OpenMP implementations failed on local arrays in target functions, so Peano used to ship a separate array-free GPU variant of these three functions (KernelsOpenMP.cpph). That copy had drifted from the CPU version and has been removed.

source() discusses the remaining annotations.

Psi4Calc()

void exagrype::ccz4::Psi4Calc	(	double *	Psi4,
		const double *const	Q,
		const double *const	gradQSerialised,
		double *	coor
	)

This function is for the calculation of psi4, a quantity related to gravitional wave.

Not used yet.

recomputeAuxiliaryVariablesFD4_4thOrder()

void exagrype::ccz4::recomputeAuxiliaryVariablesFD4_4thOrder	(	::exahype2::CellData< double, double > &	patchData,
		int	numberOfGridCellsPerPatchPerAxis,
		int	haloSize,
		int	unknowns,
		int	auxiliaryVariables
	)

Recompute auxiliary variables for FD4 scheme with a 4th order scheme.

We are given a patch as well as the number of mesh cells per axis within this patch. I also parameterise over the unknowns and the auxiliaryVariables. By default, we'd expect 59 and 0 here, but we want to routine to work for augmented systems, too.

The reconstruction uses a 4th order Finite Difference scheme. The implementation in exagrype::ccz4::internal::recomputeAuxiliaryVariablesFD4_4thOrder_LoopBody() realises the normal stencil

[1 -8 0 8 -1]

scaled with \( \frac{1}{12h} \). There is a fundamental problem with this variant which is worth mentioning: To evaluate this stencil, we need a halo layer of at least two. For a given haloSize, we can never update the outermost two layers of the code.

Shortcomings

As we use a fourth order scheme, we need two neighbours around each point to construct the outcome. That means, when we are given a patch with size 3 (left illustation above), we effectively can only calculate the auxiliary variables within the patch (green) and the first halo layer around it (brownish). This is illustrated to the right.

There's one further complexity: we usually do not have the digonal values in ExaHyPE. What we really get is not the illustration above to the left but the one in the centre. The diagonal blocks hold garbage. As a consequence, the auxiliary data in the brown data layer to the right is not properly computed. In the left brown layer, only the x-derivatives are properly reconstructed. The y-derivatives contain garbage.

Parameters

patchData	Host the actual patch data. For this function, we only have to ensure that the QIn (data plus halo) are properly set, and that the patch size in patchData is correct. All the other (meta data) properties have no influence and can hold any data.
numberOfGridCellsPerPatchPerAxis	Has to be at least 3, as the halo is also at least three. More general, has to be bigger or equal to haloSize.
haloSize	Has to be at least 3.
unknowns	Typically 59, but can be bigger if you hold additional quantities within the PDE.
auxiliaryVariables	Typically 0.

See also

exagrype::ccz4::internal::recomputeAuxiliaryVariablesFD4_4thOrder_LoopBody()

recomputeAuxiliaryVariablesFD4_centralDifferences()

void exagrype::ccz4::recomputeAuxiliaryVariablesFD4_centralDifferences	(	::exahype2::CellData< double, double > &	patchData,
		int	numberOfGridCellsPerPatchPerAxis,
		int	haloSize,
		int	unknowns,
		int	auxiliaryVariables
	)

recomputeAuxiliaryVariablesFD4_leftDifferences()

void exagrype::ccz4::recomputeAuxiliaryVariablesFD4_leftDifferences	(	::exahype2::CellData< double, double > &	patchData,
		int	numberOfGridCellsPerPatchPerAxis,
		int	haloSize,
		int	unknowns,
		int	auxiliaryVariables
	)

recomputeAuxiliaryVariablesFD4_rightDifferences()

void exagrype::ccz4::recomputeAuxiliaryVariablesFD4_rightDifferences	(	::exahype2::CellData< double, double > &	patchData,
		int	numberOfGridCellsPerPatchPerAxis,
		int	haloSize,
		int	unknowns,
		int	auxiliaryVariables
	)

source()

void exagrype::ccz4::source	(	double *	S,
		const double *const	Q,
		const int	CCZ4LapseType,
		const double	CCZ4ds,
		const double	CCZ4c,
		const double	CCZ4e,
		const double	CCZ4f,
		const double	CCZ4bs,
		const double	CCZ4sk,
		const double	CCZ4xi,
		const double	CCZ4itau,
		const double	CCZ4eta,
		const double	CCZ4k1,
		const double	CCZ4k2,
		const double	CCZ4k3,
		const double	CCZ4SO
	)

The source term is one out of two terms that we use in our CCZ4 formulation.

As it is used within the compute kernel, we also want to use it on the GPU, and we want to vectorise over it aggressively. To make this work, we have to inline the function plus ensure that the interna of the function are not vectorised before we inline the routine.

So we first declare it as simd and then we also say explicitly that it should be inlined. In my opinion, the inlining should imply that we wanna use it as SIMD, but my guess is that the simd statement ensures that the compiler doesn't prematurely vectorise.

See also

ncp()

Referenced by applications::exahype2::ccz4::CCZ4::sourceTerm(), and benchmarks::exahype2::ccz4::CCZ4::sourceTerm().

Here is the caller graph for this function:

TestingOutput()

void exagrype::ccz4::TestingOutput	(	double *	terms,
		const double *const	Q,
		const double *const	gradQSerialised
	)

A temporary test function, to output some testing values 0,1 entries: Hamilton constraint related term in theta, 2 terms: RPlusTwoNablaZSrc, pure Src 2-10 entry: no-symmetry version of R_ij.

ThetaOutputNCP()

void exagrype::ccz4::ThetaOutputNCP	(	double *	NCPterm,
		const double *const	Q,
		const double *const	gradQSerialised,
		int	normal
	)

A temporary test function, to output Hamilton constraint related term in theta, 1 terms: RPlusTwoNablaZNCP.