ExaHyPE-Engine issueshttps://gitlab.lrz.de/exahype/ExaHyPE-Engine/-/issues2018-06-15T15:13:43+02:00https://gitlab.lrz.de/exahype/ExaHyPE-Engine/-/issues/62Non cubic simulation area no more working2018-06-15T15:13:43+02:00Ghost UserNon cubic simulation area no more workingWhen I start the [SRHD_FV.exahype](https://gitlab.lrz.de/gi26det/ExaHyPE/blob/master/Code/ApplicationExamples/SRHD_FV.exahype) with non-cubic simulation domain, ie.
```
computational-domain
dimension = 2
width ...When I start the [SRHD_FV.exahype](https://gitlab.lrz.de/gi26det/ExaHyPE/blob/master/Code/ApplicationExamples/SRHD_FV.exahype) with non-cubic simulation domain, ie.
```
computational-domain
dimension = 2
width = 1.0, 0.2
offset = 0.0, 0.0
end-time = 1.0
end computational-domain
```
then I get this funny error while Peano tries to setup the initial grid:
```
...
0.0305133 info peano::utils::UserInterface::writeHeader() Application based upon the PDE framework Peano - 3rd Generation
0.0305371 info peano::utils::UserInterface::writeHeader() build: dim=2
0.030556 info peano::utils::UserInterface::writeHeader() optimisations: d-loop persistent-attributes packed opt-static-subtrees recursion-unrolling
0.030575 info peano::utils::UserInterface::writeHeader() (C) 2005 - 2016 www.peano-framework.org
0.0305951 info peano::utils::UserInterface::writeHeader() processes: 1, threads: 1
0.579099 info exahype::mappings::LoadBalancing::endIteration(State) memoryUsage =39 MB
0.579301 info exahype::runners::Runner::createGrid() grid setup iteration #1, max-level=4, state=(mergeMode:MergeNothing,sendMode:ReduceAndMergeTimeStepData,reinitTimeStepData:0,fuseADERDGPhases:0,stabilityConditionOfOneSolverWasViolated:114,timeStepSizeWeightForPredictionRerun:1.14324e+243,minMeshWidth:[0.037037,0.037037],maxMeshWidth:[1,1],numberOfInnerVertices:208,numberOfBoundaryVertices:154,numberOfOuterVertices:352,numberOfInnerCells:144,numberOfOuterCells:306,numberOfInnerLeafVertices:104,numberOfBoundaryLeafVertices:64,numberOfOuterLeafVertices:176,numberOfInnerLeafCells:135,numberOfOuterLeafCells:266,maxLevel:4,hasRefined:1,hasTriggeredRefinementForNextIteration:1,hasErased:0,hasTriggeredEraseForNextIteration:0,hasChangedVertexOrCellState:1,hasModifiedGridInPreviousIteration:1,isTraversalInverted:1), idle-nodes=1
1.21773 info exahype::mappings::LoadBalancing::endIteration(State) memoryUsage =39 MB
1.21792 info exahype::runners::Runner::createGrid() grid setup iteration #2, max-level=4, state=(mergeMode:MergeNothing,sendMode:ReduceAndMergeTimeStepData,reinitTimeStepData:0,fuseADERDGPhases:0,stabilityConditionOfOneSolverWasViolated:114,timeStepSizeWeightForPredictionRerun:1.14324e+243,minMeshWidth:[0.037037,0.037037],maxMeshWidth:[1,1],numberOfInnerVertices:208,numberOfBoundaryVertices:154,numberOfOuterVertices:548,numberOfInnerCells:144,numberOfOuterCells:423,numberOfInnerLeafVertices:104,numberOfBoundaryLeafVertices:64,numberOfOuterLeafVertices:270,numberOfInnerLeafCells:135,numberOfOuterLeafCells:370,maxLevel:4,hasRefined:1,hasTriggeredRefinementForNextIteration:0,hasErased:0,hasTriggeredEraseForNextIteration:0,hasChangedVertexOrCellState:0,hasModifiedGridInPreviousIteration:1,isTraversalInverted:0), idle-nodes=1
assertion in file /home/sven/numrel/exahype/master/Code/./Peano/peano/parallel/loadbalancing/Oracle.cpp, line 227 failed: _currentOracle>=0
ExaHyPE-SRHD: /home/sven/numrel/exahype/master/Code/./Peano/peano/parallel/loadbalancing/Oracle.cpp:227: int peano::parallel::loadbalancing::Oracle::getRegularLevelAlongBoundary() const: Assertion `false' failed.
Abgebrochen (Speicherabzug geschrieben)
```
This is bad when trying to run effective 1D simulations. :-(Tobias WeinzierlTobias Weinzierlhttps://gitlab.lrz.de/exahype/ExaHyPE-Engine/-/issues/61Compile time constants should be available as static const2018-06-15T15:13:43+02:00Ghost UserCompile time constants should be available as static constCurrently, we do have
```c
int exahype::solvers::Solver::getNumberOfVariables() const {
return _numberOfVariables;
}
```
as a method in the Solvers. However, this is not accessible at user kernels as the Solver object is no more acce...Currently, we do have
```c
int exahype::solvers::Solver::getNumberOfVariables() const {
return _numberOfVariables;
}
```
as a method in the Solvers. However, this is not accessible at user kernels as the Solver object is no more accessible. As these are compile time constants, they should be generated by the code, ie. in a file like [GeneratedConstants.h](https://gitlab.lrz.de/gi26det/ExaHyPE/blob/85cc9213e8fa71aea55e7063067285d68fae965d/Code/ApplicationExamples/SRHD/GeneratedConstants.h):
```c
#ifndef __MY_GENERATED_CONSTANTS__
#define __MY_GENERATED_CONSTANTS__
/**
* These constants should be created by the toolkit instead
* of scattering numbers around in the code. The practice to
* write naked numbers somewhere, as in
* ADERDGSolver("SRHDSolver", 3, 2, ...)
* is called "magic numbers" and they are accepted as bad
* coding practice.
*
* As we currently cannot run the toolkit for SRHD,
* these constants have to be always kept equal to the
* toolkit.
*
**/
static const int MY_POLYNOMIAL_DEGREE = 1;
static const int MY_NUMBER_OF_VARIABLES = 5;
static const int MY_NUMBER_OF_PARAMETERS = 0;
#endif /* __MY_GENERATED_CONSTANTS__ */
```https://gitlab.lrz.de/exahype/ExaHyPE-Engine/-/issues/60Spec file parsing errors in comments2018-06-15T15:13:43+02:00Ghost UserSpec file parsing errors in commentsThis specification file does not compile:
```
....
optimisation
fuse-algorithmic-steps = on
fuse-algorithmic-steps-factor = 0.99
end optimisation
/*
Q0 = das skalare Feld "dens", die relativistische ruh...This specification file does not compile:
```
....
optimisation
fuse-algorithmic-steps = on
fuse-algorithmic-steps-factor = 0.99
end optimisation
/*
Q0 = das skalare Feld "dens", die relativistische ruhemasse,
(Q1,Q2,Q3) = den relativistischen Geschwindigkeitsvektor (Impuls),
Q4 = das skalare Feld der inneren Energie (Hydrodynamik),
(Q5,Q6,Q7) = das Magnetfeld (Vektorfeld, für SRHD=0).
Q8 = Die Divergenz des Magnetfelds, die verschwinden muss und ein Mass fuer die numerische Exaktheit ist
*/
solver ADER-DG MHDSolver
variables = 9
parameters = 0
order = 3
...
```
This is due to the comment. **Comments in the specification file are not working like comments but they are parsed**. The comment was actually added by Tobias himself in https://gitlab.lrz.de/gi26det/ExaHyPE/commit/3f7028f1989ce5367d92336a09bb25fc640eca43 and when doing this, the toolkit fails with the unhelpful message `IOException: "Pushback buffer overflow"`, without mentioning any line number or so.
This is a real drawback of the proprietary configuration file format, as already mentioned by Fabian in https://gitlab.lrz.de/gi26det/ExaHyPE/issues/56#note_16127. If we had instead a simple configuration format as INI, JSON, YAML or so, such errors would not occur and instead, we could have fail proof comments like `/* comment */`, `// comment` or `# comment`.https://gitlab.lrz.de/exahype/ExaHyPE-Engine/-/issues/58Reduce memory footprint2019-09-20T15:46:51+02:00Ghost UserReduce memory footprint# Open issues
* We have consecutive heap indices for volume data and face data. We thus need to store one index for each
and can get the others by incrementing the index up to the bound we know of as developers.
We distinguish b...# Open issues
* We have consecutive heap indices for volume data and face data. We thus need to store one index for each
and can get the others by incrementing the index up to the bound we know of as developers.
We distinguish between cell and face data since we have helper cells that do not allocate cell data but face data.
* We can remove prediction and volumeFlux fields completely if we perform also the time integration
in the volume integral and the boundary extrapolation routines.
This would further make it easier to switch between global and local time stepping.
Here, we would just load a different kernel for the boundary extrapolation and allocate space-time face data
if the user switches local time stepping on.
* Allocate all the temporary arrays like the rhs, lQi_old, etc. only once per Thread and not dynamically
during the kernel calls. **(This could be done easily now in each solver!)**
* Create one "big" ADERDGTimeStep function in kernels/solver. This might help the compiler/is more Cache-friendly.
# Done
I think the following has been Tobias' idea originally;
It is not necessary to store temporary data on the heap for every cell description.
We need to analyse which ADER-DG fields are temporary and which
need to be stored persistently on the heap.
From my point of view, the following variables are temporary:
* spaceTimePredictor
* predictor
* spaceTimeVolumeFlux (includes sources)
* volumeFlux (includes sources)
The spacetime fields have a massive memory footprint.
They scale with (N+1)^{d+1} and d*(N+1)^{d+1}.
I thus propose that we assign each thread its own spaceTimePredictor spaceTimeVolumeFlux, predictor, and volume flux fields
and remove the fields from the heap cell descriptions.
This would reduce the memory footprint of the ADER-DG method dramatically (and might further lead to more cache-friendly code ?).
In a second step, we should kick out the volumeFlux field completely, don't do the time integration of the spaceTimeVolumeFlux,
and directly perform the volume integral with the spaceTimeVolumeFlux.
Implementation details
* Allocate arrays in Prediction mapping for each threadhttps://gitlab.lrz.de/exahype/ExaHyPE-Engine/-/issues/57Finite volumes solver / Limiter2019-09-20T15:46:56+02:00Ghost UserFinite volumes solver / Limiter# Treating NaNs in the update vector
If NaNs occur in the update vector, a rollback
using the update vector is not possible anymore.
We thus need to introduce a "previousSolution" field
to the ADER-DG solution.
We might be able to g...# Treating NaNs in the update vector
If NaNs occur in the update vector, a rollback
using the update vector is not possible anymore.
We thus need to introduce a "previousSolution" field
to the ADER-DG solution.
We might be able to get rid of the "update" field if we directly
add a weighted (by dt and quad weight) update to the "solution" field
from the volume integral and surface integral evaluation.
# Status of the implementation of the limiting ADER-DG scheme (unordered notes)
* The limiter workflow works now for uniform meshes with Intel's TBBs.
[sod-shock-tube_limiting-aderdg-P_3_Godunov-N_7_dmp-only.avi](/uploads/b71d24c7b82f805c604f31c1b322fcae/sod-shock-tube_limiting-aderdg-P_3_Godunov-N_7_dmp-only.avi)
* Found and fixed some more bugs. Now everything looks fine. The limiter has a local characteristic and only
requires a rollback/reallocation of memory in every fifth time step (49 out of 255) in the experiment
shown [sod-shock-tube_limiting-aderdg-P_3_Godunov-N_7_dmp-only.avi](/uploads/1e065c7c5c2a52aefb551ef0a03f0ac5/sod-shock-tube_limiting-aderdg-P_3_Godunov-N_7_dmp-only.avi), where we ran an experiment with P=3 order ADER-DG approx, \delta_0=1e-2, \epsilon=1e-3.
As long as the limiter domain does not change, we do not need to reallocate new memory and do a
recomputation of certain cells in our novel implementation.
* Find another simulation using a 9-th order ADER-DG approximation and parameters \delta_0=10^-4 , \epsilon=10^-3 here:
[sod-shock-tube_limiting-aderdg-P_9_Godunov-N_19_dmp_pad.avi](/uploads/fb466322f7df5339790aac4afca9ce6c/sod-shock-tube_limiting-aderdg-P_9_Godunov-N_19_dmp_pad.avi)
~~We observee a strange x velocity in this benchmark. Can't be seen in the video. The both profile of the Shock-Tube differs signifantly
in vicinity of the contact discontinuity.~~This is actually the momentum density. Everything is fine.
* Explosion problem with P=5 and delta_0=1e-4 and epsilon=1e-3:
[explosion_limiting-aderdg-P_5_Godunov-N_11_dmp_pad.avi](/uploads/4ace50cdcbfcbc1706432d0cc180fcea/explosion_limiting-aderdg-P_5_Godunov-N_11_dmp_pad.avi). Here the limiter domain must be adjusted after most of the time steps
* I performed some further optimisations of the DMP evaluation. I found that we can easily perform a
"loop" fusion of the min and max computation and the DMP.
* Implemented the physical admissibility detection (PAD) now as well and use it to
detect non-physical oscillations in the initial conditions. Here, I found that we can directly
pass it the solution min and max we computed as part of the DMP calculation. No need
to loop over the whole degrees of freedom and the Lobatto nodes again.
* ~~There is still an issue with the discrete maximum principle which detects too many
cells to be troubled.~~ Was resolved by tuning the DMP parameters.
* AMR and Limiting need to be combined. This is in principle "a simple" wiring of FV patches
along the ADER-DG tree. We further need spatial averaging and interpolation operators.
Then, we can follow the AMR timestepping implementation of the ADER-DG scheme.
# Writing a Godunov type first order method
Due to the issues I encountered while [rewriting](#issueswithrewritingthefinitevolumessolver) the
pseudo-TVD donor cell type finite volumes solver, Tobias and me decided
to write a simple first order Godunov type finite volumes solver that only exchanges volume averages/fluxes
between direct (face) neighbours.
Replacement of this simple method by a more complex one will be tackled in later stages of the project -- if necessary.
# Issues with higher order finite volumes solvers
* Higher-order methods have a reconstruction stencil which is larger than the
Riemann solver stencil (simple star).
* Reconstruction and time evolution of boundary extrapolated values at one face of a patch can only be performed
after all volume averages from the direct neighbour and corner neighbour patches are available.
I have identified the following phases of the FVM solver now:
* Gather: Just get the neighbour values (arithmetic intensity = 0). This will move into the Merging mapping.
* Spatial reconstruction: Compute spatial slopes in the cells or something similar. This will move into the SolutionUpdate mapping.
* Temporal evolution of boundary-extrapolated values: This will move into the SolutionUpdate mapping.
* Riemann solve: This will move into the SolutionUpdate mapping.
* Solution update; This will move into SolutionUpdate mapping.
* Extrapolation of volume averages: This will send layers of the volume averages to the neighbours.
We send layers instead of a single layer in order to only have one data exchange instead of two.
With the above strategy, we can merge FVM solvers into the current framework. The difference to the ADER-DG solver
is that the computational intensity of the neighbour merging operation is zero.
We will further need to consider neighbour data exchange over edges (3-d) and corners in our
merging scheme. Neighbour data exchange over edges in 3-d will require additional falgs.
Exchange over corners does not require additional flags.
# Issues with rewriting the finite volumes solver
I tried to decompose our pseudo-TVD donor cell type finite volumes
solver into a Riemann solve and a solution update part,
and ran into the following issues.
## Open issues with the donor cell type pseudo-TVD finite volumes solver implementation:
* The current solver uses updated solution values in the solution update.
We have to distinguish between old values and new values or
have to introduce an update.
* We need to consider two layers of the neighbour to compute all extrapolated boundary
values (wLx,wRx,wLy,wRy). Currently only one layer ist considered. This means
that the boundary extrapolated values and thus the face fluxes at the outermost faces
are computed wrongly.
* To tackle the above problem, we either need to exchange a single cell of
the diagonal neighbours or we need to rely on two data exchanges.
Tobias told me there exists however a trick to circumvent this:
Reconstructing the diagonal neighbours' contributions with values from
the direct neighbours.
## Limitations of the donor cell type pseudo-TVD finite volumes solver:
* The solver is not TVD in multiple dimensions.
* We do neither perform corner correction nor (dimensional) operator splitting. We should
thus observe loss of mass conservation, large dispersion errors, and a reduced CFL stability limit.
The reduced CFL stability limit was taken into account in our implementation.
This is also a problem of the multi-dim. Godunov method in the unsplit form.
## Follow up: Low-order Euler-DG patch with direct coupling to ADER-DG method?
## Dissemination
For the following benchmarks I always used copy boundary conditions which are equal to
outflow/inflow boundary conditions as long as everything flows out of the domain.
* Euler (Compressible Hydrodynamics)
* Explosion Problem with 27^2 grid, P=3, and N=2*3+1 FV subgrid:
![lim-aderdg_explosion](/uploads/6459f9b4a775a482cdfec24e129c56ae/lim-aderdg_explosion.png)
[lim-aderdg_explosion.avi](/uploads/6a38fa00a6bf41fe3aefbee380a03396/lim-aderdg_explosion.avi)
* Sod Shock Tube with 27^2 grid, P=3 and N=2*3+1 FV subgrid:
![lim-aderdg_euler_sod-shock-tube](/uploads/eb2a5e57ca8f37665c2dffe48dbf1e22/lim-aderdg_euler_sod-shock-tube.png)
[lim-aderdg_euler_sod-shock-tube.avi](/uploads/61dfedf7db4ead996d7fd6b9796ba5c0/lim-aderdg_euler_sod-shock-tube.avi)
![hires_antialias](/uploads/26843a0f50e52ded919cfc12b777360b/draft2_hires_antialias.png)
* SRMHD (Special Relativistic Magnetohydrodynamics; basically: Special Relativistic Euler+Maxwell)
* Blast Wave setup as in 10.1016/j.cpc.2014.03.018 with 27^2 grid, P=5, and N=2*5+1 FV subgrid:
![lim-aderg_mhd__blast-wave](/uploads/feb026167c7c2d882c132f3b775c37ac/lim-aderg_mhd__blast-wave.png)
[lim-aderdg_mhd_blast-wave.avi] (/uploads/4d715c47a03bd3850d3c0398e39da58c/lim-aderdg_mhd_blast-wave.avi)
* Rotor setup as in http://adsabs.harvard.edu/abs/2004MSAIS...4...36D with 9^2 grid, P=9, and N=2*9+1 FV subgrid:
![lim-aderdg_mhd_rotor_9x9_P9](/uploads/b35e7dc2dc137068e7416a1952d02a1a/lim-aderdg_mhd_rotor_9x9_P9.png)
[lim-aderdg_mhd_rotor_9x9_P9.avi](/uploads/c5ee8ff29f18a7201be66a1a26b664dc/lim-aderdg_mhd_rotor_9x9_P9.avi)
https://gitlab.lrz.de/exahype/ExaHyPE-Engine/-/issues/51Clean up Cell.cpp and mappings2016-12-21T16:01:26+01:00Ghost UserClean up Cell.cpp and mappingsExplicit separation ADER-DG and FV solvers. They are not distinguished by a CellDescription field "solverType" anymore.Explicit separation ADER-DG and FV solvers. They are not distinguished by a CellDescription field "solverType" anymore.https://gitlab.lrz.de/exahype/ExaHyPE-Engine/-/issues/49More control over initial data creation2017-01-27T12:27:52+01:00Ghost UserMore control over initial data creationThere are plenty of cases where I have the strong whish to break out of Peano's *Hollywood principle*. Another case is that I really need more control when setting up the initial data. I know some codes which work, in pseudo code, like
...There are plenty of cases where I have the strong whish to break out of Peano's *Hollywood principle*. Another case is that I really need more control when setting up the initial data. I know some codes which work, in pseudo code, like
```
Begin Program
Startup: Populate grid with initial data.
Load files, etc, then loop over grid and set data.
Evolution:
do loop over timesteps
WaveToyC: Evolution of 3D wave equation
t = t+dt
Do the online analysis part
Write out data if neccessary
enddo
End Program.
```
This allows typical use cases for the inital data as opening files where to read initial data from (for instance, data created with the [LORENE](http://www.lorene.obspm.fr/) code), or checking parameters and precomputing quantities how to compute the intitial data. For instance, in Astrophysics there are quite easy codes to create simple space times (TOV, discs) but they cannot *just return* the wanted quantities in a
```c++
void getInitialDataAt(const double* pointCenter, double* theQplease);
```
fashion. Instead, they precompute grid transformations, potentials, etc. There has to be a place where they can do that, naturally on a per-process but not per-thread basis.
Just to turn the tables, how idiotic the following pseudo code would be:
>
> FUNCTION getIntialDataAt(a point, returning the Q):
> * somehow retrieve the initial data properties from environment variables or open and parse another config file
> * 300 lines of code computing global derived quantities from the properties, solving several integrals and sums etc.
> * Interpolation of data on a rectangular grid with dx spacing (I have programs where this interpolation takes 40 minutes)
> * Throwing away everything and just returning the quantities for the requested point.
>
> END FUNCTION
>
We would never see ExaHyPE starting to evolve something well before we get retired.https://gitlab.lrz.de/exahype/ExaHyPE-Engine/-/issues/47Peano FileWriter should create directories or quit program2017-01-27T12:28:32+01:00Ghost UserPeano FileWriter should create directories or quit programIn order to preserve a clean directory structure, my ExaHyPE specification files write the solutions into a **subdirectory**,
```
plot vtk::binary
time = 0.0
repeat = 0.004
output = ./output/solution
...In order to preserve a clean directory structure, my ExaHyPE specification files write the solutions into a **subdirectory**,
```
plot vtk::binary
time = 0.0
repeat = 0.004
output = ./output/solution
select = {all}
end plot
```
However, if this subdirectory `output/` does not exist, `tarch` (component of Peano) writes an error everytime it is supposed to write something out:
```
1265.35 error tarch::plotter::griddata::unstructured::vtk::VTKBinaryFileWriter::close() unable to write output file ./output/solution-49.vtk
```
However, the simulation does **not** stop but just runs until the end. Afterwards, there are no results anywhere, nothing has been written out. This is not helpful but wastes computing power if something went wront.
* Simple solution: An error of tarch should stop the simulation.
* Nicer solution: The code should create the directory structure for the output files. This is not complicated at all, eg. there is [Boosts `create_directories`](http://www.boost.org/doc/libs/1_57_0/libs/filesystem/doc/reference.html#create_directories). (Yes I know, Tobias does not want any dependency like boost... You know what? You could just copy the neccessary functions)Tobias WeinzierlTobias Weinzierlhttps://gitlab.lrz.de/exahype/ExaHyPE-Engine/-/issues/37Precision issues2019-04-15T12:37:23+02:00Ghost UserPrecision issuesI observed that we can only apply the restriction and prolongation operations
a certain number of times until the projected constant values in my test are not equal anymore
to the reference values (with respect to a tolerance of 1e-12)....I observed that we can only apply the restriction and prolongation operations
a certain number of times until the projected constant values in my test are not equal anymore
to the reference values (with respect to a tolerance of 1e-12).
For more infos see the todos in 2D and 3D definitions of
GenericEulerKernelTest::testFaceUnknownsProjection(),
and
GenericEulerKernelTest::testVolumeUnknownsProjection()https://gitlab.lrz.de/exahype/ExaHyPE-Engine/-/issues/35Naming convention for arguments of Solver.h functions vs. template functions ...2018-03-02T13:07:20+01:00Ghost UserNaming convention for arguments of Solver.h functions vs. template functions for generic kernelsThe current naming convention is confusing and should be rethought.
<p>
We need to emphasise that the user flux, eigenvalue and initial value functions for generic kernels
do pointwise evaluations while the actual Solver.h functions lik...The current naming convention is confusing and should be rethought.
<p>
We need to emphasise that the user flux, eigenvalue and initial value functions for generic kernels
do pointwise evaluations while the actual Solver.h functions like the refinement criterion
work with all basis function coefficients of a cell.
<p/>
<p>
This issue led to a lot of confusion on the user side.
<p/>