04.06., 9:00 - 12:00: GitLab will be migrated to a new server environment and upgraded to Enterprise Edition Ultimate. The estimated downtime will be 2-3 hours. Please see https://doku.lrz.de/display/PUBLIC/GitLab+Ultimate+Migration for more details about changes related to the migration.

Commit 909b848c authored by Jean-Matthieu's avatar Jean-Matthieu

Add mixed Euler testcase

parent 078d829d
{
"project_name": "Euler",
"paths": {
"peano_kernel_path": "./Peano",
"exahype_path": "./ExaHyPE",
"output_directory": "./ApplicationExamples/TestCases/mixed_Euler"
},
"architecture": "hsw",
"computational_domain": {
"dimension": 2,
"time_steps": 200,
"offset": [
0.0,
0.0,
0.0
],
"width": [
1.0,
1.0,
1.0
]
},
"_shared_memory": {
"cores": 4,
"properties_file": "sharedmemory.properties",
"autotuning_strategy": "dummy"
},
"_distributed_memory": {
"timeout": 60,
"load_balancing_type": "static",
"buffer_size": 64,
"load_balancing_strategy": "hotspot",
"node_pool_strategy": "fair",
"ranks_per_node": 1
},
"optimisation": {
"fuse_algorithmic_steps": "none",
"fuse_algorithmic_steps_rerun_factor": 0.99,
"fuse_algorithmic_steps_diffusion_factor": 0.99,
"spawn_predictor_as_background_thread": false,
"spawn_amr_background_threads": false,
"disable_vertex_exchange_in_time_steps": true,
"time_step_batch_factor": 0.5,
"disable_metadata_exchange_in_batched_time_steps": true,
"double_compression": 0.0,
"spawn_double_compression_as_background_thread": false
},
"solvers": [
{
"type": "ADER-DG",
"name": "EulerSolver_ADERDG",
"order": 7,
"maximum_mesh_size": 0.05,
"time_stepping": "globalfixed",
"aderdg_kernel": {
"language": "C",
"nonlinear": true,
"terms": [
"flux", "ncp"
],
"space_time_predictor": {},
"optimised_terms": [],
"optimised_kernel_debugging": [],
"implementation": "generic"
},
"variables": [
{
"name": "rho",
"multiplicity": 1
},
{
"name": "j",
"multiplicity": 3
},
{
"name": "E",
"multiplicity": 1
}
],
"parameters": {
"reference": "entropywave"
},
"plotters": [
{
"type": "vtk::Cartesian::vertices::ascii",
"name": "ErrorPlotter",
"time": 0.0,
"repeat": 0.01,
"output": "./errors",
"variables": 15
},
{
"type": "vtk::Legendre::cells::ascii",
"name": "ConservedQuantitiesWriter",
"time": 0.0,
"repeat": 0.01,
"output": "./conserved",
"variables": 5
}
]
}
]
}
\ No newline at end of file
#include "ErrorPlotter.h"
#include "EulerSolver_ADERDG.h"
Euler::ErrorPlotter::ErrorPlotter(EulerSolver_ADERDG& solver) {
// @todo Please insert your code here
}
Euler::ErrorPlotter::~ErrorPlotter() {
// @todo Please insert your code here
}
void Euler::ErrorPlotter::startPlotting(double time) {
// @todo Please insert your code here
}
void Euler::ErrorPlotter::finishPlotting() {
// @todo Please insert your code here
}
void Euler::ErrorPlotter::mapQuantities(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& x,
const tarch::la::Vector<DIMENSIONS, int>& pos,
double* Q,
double* outputQuantities,
double timeStamp
) {
constexpr int numberOfVariables = AbstractEulerSolver_ADERDG::NumberOfVariables;
double QAna[numberOfVariables];
EulerSolver_ADERDG::referenceSolution(x.data(),timeStamp,QAna);
for (int i=0; i<numberOfVariables; i++){
outputQuantities[i] =std::abs( QAna[i] - Q[i] );
}
for (int i=0; i<numberOfVariables; i++){
outputQuantities[i+numberOfVariables] = Q[i];
}
for (int i=0; i<numberOfVariables; i++){
outputQuantities[i+2*numberOfVariables] = QAna[i];
}
}
// This file is generated by the ExaHyPE toolkit.
// Please do not modify - it will be overwritten by the next
// ExaHyPE toolkit call.
//
// ========================
// www.exahype.eu
// ========================
#include "exahype/plotters/Plotter.h"
namespace Euler{
class ErrorPlotter;
/**
* Forward declaration
*/
class EulerSolver_ADERDG;
}
class Euler::ErrorPlotter: public exahype::plotters::Plotter::UserOnTheFlyPostProcessing{
public:
ErrorPlotter(EulerSolver_ADERDG& solver);
virtual ~ErrorPlotter();
void startPlotting(double time) override;
void finishPlotting() override;
void mapQuantities(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& x,
const tarch::la::Vector<DIMENSIONS, int>& pos,
double* Q,
double* outputQuantities,
double timeStamp) override;
};
/**
* This file is part of the ExaHyPE project.
* Copyright (c) 2016 http://exahype.eu
* All rights reserved.
*
* The project has received funding from the European Union's Horizon
* 2020 research and innovation programme under grant agreement
* No 671698. For copyrights and licensing, please consult the webpage.
*
* Released under the BSD 3 Open Source License.
* For the full license text, see LICENSE.txt
**/
#include "EulerSolver_ADERDG.h"
#include "EulerSolver_ADERDG_Variables.h"
#include "tarch/la/MatrixVectorOperations.h"
#include <algorithm>
#include <string>
#include <math.h>
#include "peano/utils/Loop.h"
#include "kernels/KernelUtils.h"
#include "kernels/GaussLegendreBasis.h"
tarch::logging::Log Euler::EulerSolver_ADERDG::_log("Euler::EulerSolver_ADERDG");
void Euler::EulerSolver_ADERDG::init(const std::vector<std::string>& cmdlineargs,const exahype::parser::ParserView& constants) {
}
void Euler::EulerSolver_ADERDG::flux(const double* const Q, double** const F) {
constexpr double gamma = 1.4;
const double irho = 1./Q[0];
#if DIMENSIONS==3
const double j2 = Q[1]*Q[1]+Q[2]*Q[2]+Q[3]*Q[3];
#else
const double j2 = Q[1]*Q[1]+Q[2]*Q[2];
#endif
const double p = (gamma-1) * (Q[4] - 0.5*irho*j2);
// col 1
F[0][0] = Q[1];
F[0][1] = irho*Q[1]*Q[1] + p;
F[0][2] = irho*Q[2]*Q[1];
F[0][3] = irho*Q[3]*Q[1];
F[0][4] = irho*(Q[4]+p)*Q[1];
// col 2
F[1][0] = Q[2];
F[1][1] = irho*Q[1]*Q[2];
F[1][2] = irho*Q[2]*Q[2] + p;
F[1][3] = irho*Q[3]*Q[2];
F[1][4] = irho*(Q[4]+p)*Q[2];
#if DIMENSIONS==3
// col 3
F[2][0] = Q[3];
F[2][1] = irho*Q[1]*Q[3];
F[2][2] = irho*Q[2]*Q[3];
F[2][3] = irho*Q[3]*Q[3] + p;
F[2][4] = irho*(Q[4]+p)*Q[3];
#endif
}
void Euler::EulerSolver_ADERDG::nonConservativeProduct(const double* const Q, const double* const gradQ, double* const BgradQ) {
BgradQ[0] = 0.;
BgradQ[1] = 0.;
BgradQ[2] = 0.;
BgradQ[3] = 0.;
BgradQ[4] = 0.;
}
void Euler::EulerSolver_ADERDG::eigenvalues(const double* const Q,
const int direction,
double* const lambda) {
constexpr double gamma = 1.4;
const double irho = 1./Q[0];
#if DIMENSIONS==3
const double j2 = Q[1]*Q[1]+Q[2]*Q[2]+Q[3]*Q[3];
#else
const double j2 = Q[1]*Q[1]+Q[2]*Q[2];
#endif
const double p = (gamma-1) * (Q[4] - 0.5*irho*j2);
const double u_n = Q[direction + 1] * irho;
const double c = std::sqrt(gamma * p * irho);
std::fill_n(lambda,5,u_n);
lambda[0] -= c;
lambda[4] += c;
}
void Euler::EulerSolver_ADERDG::entropyWave(const double* const x,double t, double* const Q) {
const double gamma = 1.4;
constexpr double width = 0.3;
constexpr double amplitude = 0.3;
#if DIMENSIONS==2
tarch::la::Vector<DIMENSIONS,double> xVec(x[0],x[1]);
tarch::la::Vector<DIMENSIONS,double> v0(0.5,0.5);
tarch::la::Vector<DIMENSIONS,double> x0(0.5,0.5);
#else
tarch::la::Vector<DIMENSIONS,double> xVec(x[0],x[1],x[2]);
tarch::la::Vector<DIMENSIONS,double> v0(0.5,0.5,0.0);
tarch::la::Vector<DIMENSIONS,double> x0(0.5,0.5,0.5);
#endif
const double distance = tarch::la::norm2( xVec - x0 - v0 * t );
Q[0] = 0.5 + amplitude * std::exp(-distance / std::pow(width, DIMENSIONS));
Q[1] = Q[0] * v0[0];
Q[2] = Q[0] * v0[1];
Q[3] = 0.0;
// total energy = internal energy + kinetic energy
const double p = 1.;
Q[4] = p / (gamma-1) + 0.5*Q[0] * (v0[0]*v0[0]+v0[1]*v0[1]); // v*v; assumes: v0[2]=0
}
void Euler::EulerSolver_ADERDG::referenceSolution(const double* const x,double t, double* const Q) {
entropyWave(x,t,Q);
}
void Euler::EulerSolver_ADERDG::adjustPointSolution(const double* const x,const double t,const double dt, double* const Q) {
if (tarch::la::equals(t, 0.0)) {
referenceSolution(x,0.0,Q);
}
}
exahype::solvers::Solver::RefinementControl
Euler::EulerSolver_ADERDG::refinementCriterion(
const double* const luh, const tarch::la::Vector<DIMENSIONS, double>& center,
const tarch::la::Vector<DIMENSIONS, double>& dx, double t,
const int level) {
double maxE = -std::numeric_limits<double>::max();
double minE = +std::numeric_limits<double>::max();
const int nodes = std::pow((Order+1), DIMENSIONS);
for(int i=0;i<nodes;i++) {
maxE = std::max(maxE, luh[i*NumberOfVariables+NumberOfVariables-1]);
minE = std::min(minE, luh[i*NumberOfVariables+NumberOfVariables-1]);
}
if ( maxE/minE > 1.05 ) {
return RefinementControl::Refine;
} else if ( level > getCoarsestMeshLevel() ) {
return RefinementControl::Erase;
} else {
return RefinementControl::Keep;
}
}
void Euler::EulerSolver_ADERDG::boundaryValues(const double* const x,const double t,const double dt,const int faceIndex,const int direction,const double* const fluxIn,const double* const stateIn,const double* const gradStateIn,double* const fluxOut,double* const stateOut) {
double Q[NumberOfVariables] = {0.0};
double _F[3][NumberOfVariables] = {0.0};
double* F[3] = {_F[0],_F[1],_F[2]};
// initialise
std::fill_n(stateOut, NumberOfVariables, 0.0);
std::fill_n(fluxOut, NumberOfVariables, 0.0);
for (int i=0; i<Order+1; i++) {
const double ti = t + dt * kernels::legendre::nodes[Order][i];
referenceSolution(x,ti,Q);
flux(Q,F);
for (int v=0; v<NumberOfVariables; v++) {
stateOut[v] += Q[v] * kernels::legendre::weights[Order][i];
fluxOut[v] += F[direction][v] * kernels::legendre::weights[Order][i];
}
}
}
void Euler::EulerSolver_ADERDG::mapDiscreteMaximumPrincipleObservables(double* const observables, const double* const Q) const {
std::copy_n(Q, NumberOfVariables, observables);
}
bool Euler::EulerSolver_ADERDG::isPhysicallyAdmissible(
const double* const solution,
const double* const observablesMin,const double* const observablesMax,
const bool wasTroubledInPreviousTimeStep,
const tarch::la::Vector<DIMENSIONS,double>& center,
const tarch::la::Vector<DIMENSIONS,double>& dx,
const double t) const {
if (observablesMin[0] <= 0.0) return false;
if (observablesMin[4] < 0.0) return false;
return true;
}
#ifndef __EulerSolver_ADERDG_CLASS_HEADER__
#define __EulerSolver_ADERDG_CLASS_HEADER__
// This file is generated by the ExaHyPE toolkit.
// Please do not modify - it will be overwritten by the next
// ExaHyPE toolkit call.
//
// ========================
// www.exahype.eu
// ========================
#include <ostream>
#include "AbstractEulerSolver_ADERDG.h"
#include "exahype/parser/ParserView.h"
namespace Euler{
class EulerSolver_ADERDG;
}
class Euler::EulerSolver_ADERDG: public Euler::AbstractEulerSolver_ADERDG {
private:
/**
* (Smooth solution)
*
* Entropy wave is a moving Gaussian matter distribution where it is simple
* to give an analytic result.
*
* See also chapter 7.13.2 in "I do like CFD, VOL.1" by Katate Masatsuka.
*/
static void entropyWave(const double* const x,double t, double* const Q);
/**
* Log device
*/
static tarch::logging::Log _log;
public:
EulerSolver_ADERDG(const double maximumMeshSize,int maximumAdaptiveMeshDepth,const int haloCells,const int haloBufferCells,const int limiterBufferCells,const int regularisedFineGridLevels,const exahype::solvers::Solver::TimeStepping timeStepping,const int DMPObservables);
/**
* Initialise the solver.
*
* \param[in] cmdlineargs the command line arguments.
*/
void init(const std::vector<std::string>& cmdlineargs,const exahype::parser::ParserView& constants) final override;
/**
* calls ::entropyWave
* ErrorWriter, ErrorPlotter write errors of numerical solution.
*/
static void referenceSolution(const double* const x, const double t, double* const Q);
/**
* Adjust the conserved variables and parameters (together: Q) at a given time t at the (quadrature) point x.
*
* \note Use this function and ::useAdjustSolution to set initial conditions.
*
* \param[in] x the physical coordinate on the face.
* \param[in] w (deprecated) the quadrature weight corresponding to the quadrature point w.
* \param[in] t the start of the time interval.
* \param[in] dt the width of the time interval.
* \param[inout] Q the conserved variables (and parameters) associated with a quadrature point
* as C array (already allocated).
*/
void adjustPointSolution(const double* const x,const double t,const double dt,double* const Q) override;
/**
* Compute the flux tensor.
*
* \param[in] Q the conserved variables (and parameters) associated with a quadrature point
* as C array (already allocated).
* \param[inout] F the fluxes at that point as C array (already allocated).
*/
virtual void flux(const double* const Q,double** const F);
virtual void nonConservativeProduct(const double* const Q, const double* const gradQ, double* const BgradQ);
/**
* Compute the eigenvalues of the flux tensor per coordinate direction \p d.
*
* \param[in] Q the conserved variables associated with a quadrature node
* as C array (already allocated).
* \param[in] d the column of the flux vector (d=0,1,...,DIMENSIONS).
* \param[inout] lambda the eigenvalues as C array (already allocated).
*/
void eigenvalues(const double* const Q,const int d,double* const lambda);
/**
* Impose boundary conditions at a point on a boundary face
* within the time interval [t,t+dt].
*
* \param[in] x the physical coordinate on the face.
* \param[in] t the start of the time interval.
* \param[in] dt the width of the time interval.
* \param[in] faceIndex indexing of the face (0 -- {x[0]=xmin}, 1 -- {x[1]=xmax}, 2 -- {x[1]=ymin}, 3 -- {x[2]=ymax}, and so on,
* where xmin,xmax,ymin,ymax are the bounds of the cell containing point x.
* \param[in] d the coordinate direction the face normal is pointing to.
* \param[in] QIn the conserved variables at point x from inside of the domain
* and time-averaged (over [t,t+dt]) as C array (already allocated).
* \param[in] FIn the normal fluxes at point x from inside of the domain
* and time-averaged (over [t,t+dt]) as C array (already allocated).
* \param[inout] QOut the conserved variables at point x from outside of the domain
* and time-averaged (over [t,t+dt]) as C array (already allocated).
* \param[inout] FOut the normal fluxes at point x from outside of the domain
* and time-averaged (over [t,t+dt]) as C array (already allocated).
*/
void boundaryValues(const double* const x,const double t,const double dt,const int faceIndex,const int normalNonZero,const double* const fluxIn,const double* const stateIn,const double* const gradStateIn,double* const fluxOut,double* const stateOut);
/**
* Evaluate the refinement criterion within a cell.
*
* \note Instead of a variables array at a single quadrature point we give
* you all NumberOfVariables*(Order+1)^DIMENSIONS solution degrees of freedom.
*
* \note Use this function and ::adjustSolution to set initial conditions.
*
* \param[in] centre The centre of the cell.
* \param[in] dx The extent of the cell.
* \param[in] t the start of the time interval.
* \param[in] dt the width of the time interval.
* \return One of exahype::solvers::Solver::RefinementControl::{Erase,Keep,Refine}.
*/
exahype::solvers::Solver::RefinementControl refinementCriterion(const double* const luh,const tarch::la::Vector<DIMENSIONS,double>& centre,const tarch::la::Vector<DIMENSIONS,double>& dx,double t,const int level) override;
void mapDiscreteMaximumPrincipleObservables(double* const observables, const double* const Q) const override;
bool isPhysicallyAdmissible(
const double* const solution,
const double* const observablesMin,const double* const observablesMax,
const bool wasTroubledInPreviousTimeStep,
const tarch::la::Vector<DIMENSIONS,double>& center,
const tarch::la::Vector<DIMENSIONS,double>& dx,
const double t) const override;
};
#endif // __EulerSolver_ADERDG_CLASS_HEADER__
# MPI Logfilter.
# black or white list entry
# -1 means all ranks
# default entries
debug tarch -1 black
debug peano -1 black
info tarch -1 black
info tarch::logging::CommandLineLogger -1 white
info peano -1 black
info peano::utils::UserInterface -1 white
info peano::performanceanalysis -1 white
info exahype -1 white
# on first node
info mpibalancing -1 black
# too many messages
info peano::parallel::SendReceiveBufferAbstractImplementation::releaseSentMessages -1 black
info sharedmemoryoracles -1 black
info Z4Solver -1 black
info exahype::solvers::Plotter -1 black
info exahype::runner::Runner::runAsWorker -1 black
info exahype::mappings::LoadBalancing -1 white
info exahype::runners::Runner::runAsWorker -1 black
info exahype::runners::Runner::init -1 black
info exahype::runners::Runner::createRepository -1 white
info exahype::mappings::LoadBalancing::endIteration -1 black
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment