11.08., 9:00 - 11:00: Due to updates GitLab will be unavailable for some minutes between 09:00 and 11:00.

Commit 954bc81d authored by Jean-Matthieu Gallard's avatar Jean-Matthieu Gallard

Mixed Euler - update test case to use predictor recompute + vect pde

parent bdaa2308
......@@ -8,7 +8,7 @@
"architecture": "hsw",
"computational_domain": {
"dimension": 2,
"time_steps": 205,
"time_steps": 50,
"offset": [
0.0,
0.0,
......@@ -58,10 +58,10 @@
"terms": [
"flux", "ncp"
],
"space_time_predictor": {},
"space_time_predictor": {"predictor_recompute": true, "vectorise_terms":true},
"optimised_terms": [],
"optimised_kernel_debugging": [],
"implementation": "generic"
"implementation": "optimised"
},
"variables": [
{
......@@ -85,17 +85,9 @@
"type": "vtk::Cartesian::vertices::ascii",
"name": "ErrorPlotter",
"time": 0.0,
"repeat": 0.01,
"repeat": 0.001,
"output": "./errors",
"variables": 15
},
{
"type": "vtk::Legendre::cells::ascii",
"name": "ConservedQuantitiesWriter",
"time": 0.0,
"repeat": 0.01,
"output": "./conserved",
"variables": 5
}
]
}
......
......@@ -34,6 +34,10 @@ 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) {
}
/*******
Scalar default PDEs
*/
void Euler::EulerSolver_ADERDG::flux(const double* const Q, double** const F) {
constexpr double gamma = 1.4;
......@@ -78,6 +82,7 @@ void Euler::EulerSolver_ADERDG::flux(const double* const Q, double** const F) {
#endif
}
void Euler::EulerSolver_ADERDG::nonConservativeProduct(const double* const Q, const double* const gradQ, double* const BgradQ) {
//ncp is div(F(Q))
......@@ -117,6 +122,152 @@ void Euler::EulerSolver_ADERDG::nonConservativeProduct(const double* const Q, co
}
/*****
Recompute Predictor PDEs
*/
//Flux in header (delegate to default)
void Euler::EulerSolver_ADERDG::nonConservativeProduct2(const double* const Q,const double* const P, const double* const* const gradQ, double* const BgradQ) {
//ncp is div(F(Q))
std::memset(BgradQ,0,NumberOfVariables*sizeof(double));
constexpr double gamma = 1.4;
const double irho = 1./Q[0];
// x
// done in flux
// y
const double* const Q_dy = gradQ[1];
const double irho_dy = -Q_dy[0]*irho*irho;
#if DIMENSIONS==3
const double j2 = Q[1]*Q[1] + Q[2]*Q[2] + Q[3]*Q[3];
const double j2_dy = 2*Q[1]*Q_dy[1] + 2*Q[2]*Q_dy[2] + 2*Q[3]*Q_dy[3];
#else
const double j2 = Q[1]*Q[1] + Q[2]*Q[2];
const double j2_dy = 2*Q[1]*Q_dy[1] + 2*Q[2]*Q_dy[2];
#endif
const double p = (gamma-1) * (Q[4] - 0.5*irho*j2);
const double p_dy = (gamma-1) * (Q_dy[4] - 0.5*(j2_dy*irho+j2*irho_dy));
BgradQ[0] += Q_dy[2];
BgradQ[1] += irho_dy*Q[1]*Q[2] + irho*Q_dy[1]*Q[2] + irho*Q[1]*Q_dy[2];
BgradQ[2] += irho_dy*Q[2]*Q[2] + irho*Q_dy[2]*Q[2]*2 + p_dy;
BgradQ[3] += irho_dy*Q[3]*Q[2] + irho*Q_dy[3]*Q[2] + irho*Q[3]*Q_dy[2];
BgradQ[4] += irho_dy*Q[4]*Q[2] + irho*Q_dy[4]*Q[2] + irho*Q[4]*Q_dy[2] + irho_dy*p*Q[2] + irho*p_dy*Q[2] + irho*p*Q_dy[2];
#if DIMENSIONS==3
// z
// done in flux
#endif
}
/*****
Recompute Predictor vectorized PDEs
*/
void Euler::EulerSolver_ADERDG::flux_vect(const double* const Q,const double* const P,double** const F) {
constexpr double gamma = 1.4;
for(int i=0; i<VectLength; i++) {
const double irho = 1./Q[0*VectStride+i];
#if DIMENSIONS==3
const double j2 = Q[1*VectStride+i]*Q[1*VectStride+i]+Q[2*VectStride+i]*Q[2*VectStride+i]+Q[3*VectStride+i]*Q[3*VectStride+i];
#else
const double j2 = Q[1*VectStride+i]*Q[1*VectStride+i]+Q[2*VectStride+i]*Q[2*VectStride+i];
#endif
const double p = (gamma-1) * (Q[4*VectStride+i] - 0.5*irho*j2);
// col 1
F[0][0*VectStride+i] = Q[1*VectStride+i];
F[0][1*VectStride+i] = irho*Q[1*VectStride+i]*Q[1*VectStride+i] + p;
F[0][2*VectStride+i] = irho*Q[2*VectStride+i]*Q[1*VectStride+i];
F[0][3*VectStride+i] = irho*Q[3*VectStride+i]*Q[1*VectStride+i];
F[0][4*VectStride+i] = irho*(Q[4*VectStride+i]+p)*Q[1*VectStride+i];
// 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];
*/
// done in ncp
F[1][0*VectStride+i] = 0.;//Q[2];
F[1][1*VectStride+i] = 0.;//irho*Q[1]*Q[2];
F[1][2*VectStride+i] = 0.;//irho*Q[2]*Q[2] + p;
F[1][3*VectStride+i] = 0.;//irho*Q[3]*Q[2];
F[1][4*VectStride+i] = 0.;//irho*Q[4]*Q[2]+irho*p*Q[2];
#if DIMENSIONS==3
// col 3
F[2][0*VectStride+i] = Q[3*VectStride+i];
F[2][1*VectStride+i] = irho*Q[1*VectStride+i]*Q[3*VectStride+i];
F[2][2*VectStride+i] = irho*Q[2*VectStride+i]*Q[3*VectStride+i];
F[2][3*VectStride+i] = irho*Q[3*VectStride+i]*Q[3*VectStride+i] + p;
F[2][4*VectStride+i] = irho*(Q[4*VectStride+i]+p)*Q[3*VectStride+i];
#endif
}
}
void Euler::EulerSolver_ADERDG::nonConservativeProduct_vect(const double* const Q, const double* const P, const double* const * const gradQ, double* const BgradQ) {
//ncp is div(F(Q))
std::memset(BgradQ,0,NumberOfVariables*VectStride*sizeof(double));
constexpr double gamma = 1.4;
const double* const Q_dy = gradQ[1];
for(int i=0; i<VectLength; i++) {
const double irho = 1./Q[0*VectStride+i];
// x
// done in flux
// y
const double irho_dy = -Q_dy[0*VectStride+i]*irho*irho;
#if DIMENSIONS==3
const double j2 = Q[1*VectStride+i]*Q[1*VectStride+i] + Q[2*VectStride+i]*Q[2*VectStride+i] + Q[3*VectStride+i]*Q[3*VectStride+i];
const double j2_dy = 2*Q[1*VectStride+i]*Q_dy[1*VectStride+i] + 2*Q[2*VectStride+i]*Q_dy[2*VectStride+i] + 2*Q[3*VectStride+i]*Q_dy[3*VectStride+i];
#else
const double j2 = Q[1*VectStride+i]*Q[1*VectStride+i] + Q[2*VectStride+i]*Q[2*VectStride+i];
const double j2_dy = 2*Q[1*VectStride+i]*Q_dy[1*VectStride+i] + 2*Q[2*VectStride+i]*Q_dy[2*VectStride+i];
#endif
const double p = (gamma-1) * (Q[4*VectStride+i] - 0.5*irho*j2);
const double p_dy = (gamma-1) * (Q_dy[4*VectStride+i] - 0.5*(j2_dy*irho+j2*irho_dy));
BgradQ[0*VectStride+i] += Q_dy[2*VectStride+i];
BgradQ[1*VectStride+i] += irho_dy*Q[1*VectStride+i]*Q[2*VectStride+i] + irho*Q_dy[1*VectStride+i]*Q[2*VectStride+i] + irho*Q[1*VectStride+i]*Q_dy[2*VectStride+i];
BgradQ[2*VectStride+i] += irho_dy*Q[2*VectStride+i]*Q[2*VectStride+i] + irho*Q_dy[2*VectStride+i]*Q[2*VectStride+i]*2 + p_dy;
BgradQ[3*VectStride+i] += irho_dy*Q[3*VectStride+i]*Q[2*VectStride+i] + irho*Q_dy[3*VectStride+i]*Q[2*VectStride+i] + irho*Q[3*VectStride+i]*Q_dy[2*VectStride+i];
BgradQ[4*VectStride+i] += irho_dy*Q[4*VectStride+i]*Q[2*VectStride+i] + irho*Q_dy[4*VectStride+i]*Q[2*VectStride+i] + irho*Q[4*VectStride+i]*Q_dy[2*VectStride+i] + irho_dy*p*Q[2*VectStride+i] + irho*p_dy*Q[2*VectStride+i] + irho*p*Q_dy[2*VectStride+i];
#if DIMENSIONS==3
// z
// done in flux
#endif
}
}
/*****
Other
*/
void Euler::EulerSolver_ADERDG::eigenvalues(const double* const Q,
const int direction,
double* const lambda) {
......@@ -138,43 +289,13 @@ void Euler::EulerSolver_ADERDG::eigenvalues(const double* const Q,
lambda[4] += c;
}
// Sin waves
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::sinWave(const double* const x,double t, double* const Q) {
const double gamma = 1.4;
constexpr double width = 0.3;
constexpr double amplitude = 1.5;
#if DIMENSIONS==2
tarch::la::Vector<DIMENSIONS,double> xVec(x[0],x[1]);
tarch::la::Vector<DIMENSIONS,double> v0(2.,1.);
#else
tarch::la::Vector<DIMENSIONS,double> xVec(x[0],x[1],x[2]);
tarch::la::Vector<DIMENSIONS,double> v0(2.,1.,0.0);
#endif
constexpr double v0[3] = {2., 1., 0.};
Q[0] = 2 + amplitude/4 * sin(10*2*M_PI*(x[0]-v0[0]*t)) * (3+sin(5*2*M_PI*(x[1]-v0[1]*t)))/4;
Q[1] = Q[0] * v0[0];
......@@ -183,17 +304,18 @@ void Euler::EulerSolver_ADERDG::sinWave(const double* const x,double t, double*
// 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);
sinWave(x,t,Q);
/*
// if nVar == 8
Q[5] = 0;
Q[6] = 0;
Q[7] = 0;
*/
}
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);
entropyWave(x,0.0,Q);
}
}
......@@ -230,7 +352,7 @@ void Euler::EulerSolver_ADERDG::boundaryValues(const double* const x,const doubl
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);
entropyWave(x,ti,Q);
flux(Q,F);
for (int v=0; v<NumberOfVariables; v++) {
stateOut[v] += Q[v] * kernels::legendre::weights[Order][i];
......
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