In January 2021 we will introduce a 10 GB quota for project repositories. Higher limits for individual projects will be available on request. Please see https://doku.lrz.de/display/PUBLIC/GitLab for more information.

Commit 74557523 authored by Jean-Matthieu Gallard's avatar Jean-Matthieu Gallard

Merge branch 'jm/splitPicard' into 'master'

Jm/split picard

See merge request exahype/ExaHyPE-Engine!37
parents 15bf23fe 3d2bafda
......@@ -19,8 +19,8 @@
#include "kernels/KernelUtils.h"
#ifdef OPT_KERNELS
#include "kernels/NavierStokes_NavierStokesSolver_ADERDG/Kernels.h"
#include "kernels/NavierStokes_NavierStokesSolver_ADERDG/converter.h"
#include "kernels/NavierStokes_NavierStokesSolver_ADERDG/aderdg/Kernels.h"
#include "kernels/NavierStokes_NavierStokesSolver_ADERDG/aderdg/converter.h"
#endif
......@@ -57,6 +57,25 @@ void NavierStokes::NavierStokesSolver_ADERDG::algebraicSource(const tarch::la::V
scenario->source(x, t, ns, Q, S);
}
void NavierStokes::NavierStokesSolver_ADERDG::algebraicSource_vect(const tarch::la::Vector<DIMENSIONS, double>& x, double t, const double *const Q, const double *const P, double *S) {
// naive fallback to scalar
double Q2[NumberOfVariables + NumberOfParameters];
double S2[NumberOfVariables];
for(int i=0; i<VectLength; i++) {
for(int j=0 ; j<NumberOfVariables; j++) {
Q2[j] = Q[j*VectStride+i];
}
for(int j=0; j<NumberOfParameters; j++) {
Q2[j+NumberOfVariables] = P[j*VectStride+i];
}
algebraicSource(x,t,Q2,S2);
for(int j=0 ; j<NumberOfVariables; j++) {
S[j*VectStride+i] = S2[j];
}
}
}
void NavierStokes::NavierStokesSolver_ADERDG::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) {
constexpr auto basisSize = Order + 1;
constexpr auto NumberOfData = NumberOfVariables + NumberOfParameters;
......@@ -385,6 +404,47 @@ void NavierStokes::NavierStokesSolver_ADERDG::viscousFlux(const double *const Q,
ns.evaluateFlux(Q, gradQ, F, true);
}
void NavierStokes::NavierStokesSolver_ADERDG::viscousFlux2(const double *const Q, const double *const P, const double *const *const gradQ, double **F) {
constexpr auto gradQ2Size = NumberOfVariables*DIMENSIONS;
double gradQ2[gradQ2Size];
for(int d = 0; d < DIMENSIONS; d++) {
for(int i=0; i < NumberOfVariables; i++) {
gradQ2[i+NumberOfVariables*d] = gradQ[d][i];
}
}
ns.evaluateFlux(Q, gradQ2, F, true);
}
void NavierStokes::NavierStokesSolver_ADERDG::viscousFlux_vect(const double *const Q, const double *const P, const double *const *const gradQ, double **F) {
// naive fallback to scalar
double Q2[NumberOfVariables + NumberOfParameters];
double gradQ2[NumberOfVariables*DIMENSIONS];
double Fx[NumberOfVariables];
double Fy[NumberOfVariables];
double Fz[NumberOfVariables];
double* F2[3] = {Fx,Fy,Fz};
for(int i=0; i<VectLength; i++) {
for(int j=0 ; j<NumberOfVariables; j++) {
Q2[j] = Q[j*VectStride+i];
}
for(int j=0; j<NumberOfParameters; j++) {
Q2[j+NumberOfVariables] = P[j*VectStride+i];
}
for(int d = 0; d < DIMENSIONS; d++) {
for(int j=0; j < NumberOfVariables; j++) {
gradQ2[j+NumberOfVariables*d] = gradQ[d][j*VectStride+i];
}
}
viscousFlux(Q2,gradQ2,F2);
for(int d=0 ; d<DIMENSIONS; d++) {
for(int j=0; j < NumberOfVariables; j++) {
F[d][j*VectStride+i] = F2[d][j];
}
}
}
}
void NavierStokes::NavierStokesSolver_ADERDG::boundaryConditions( double* const fluxIn, const double* const stateIn, const double* const gradStateIn, const double* const luh, const tarch::la::Vector<DIMENSIONS, double>& cellCentre, const tarch::la::Vector<DIMENSIONS,double>& cellSize, const double t,const double dt, const int direction, const int orientation) {
#ifdef OPT_KERNELS
......
......@@ -156,9 +156,13 @@ class NavierStokes::NavierStokesSolver_ADERDG : public NavierStokes::AbstractNav
*/
//void flux(const double* const Q,double** F) final override;
void viscousFlux(const double* const Q, const double* const gradQ, double** F) final override;
void viscousFlux2(const double* const Q, const double* const P, const double* const * const gradQ, double** F);
void viscousFlux_vect(const double* const Q, const double* const P, const double* const * const gradQ, double** F);
// TODO(Lukas) Add to toolkit!
void algebraicSource(const tarch::la::Vector<DIMENSIONS, double>& x, double t, const double *const Q, double *S) override;
inline void algebraicSource2(const tarch::la::Vector<DIMENSIONS, double>& x, double t, const double *const Q, const double *const P, double *S) {algebraicSource(x,t,Q,S);}
void algebraicSource_vect(const tarch::la::Vector<DIMENSIONS, double>& x, double t, const double *const Q, const double *const P, double *S);
/* nonConservativeProduct() function is not included, as requested in the specification file */
......
......@@ -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
}
]
}
......
......@@ -21,6 +21,7 @@
#include <string>
#include <math.h>
#include <cmath> //sin
#include "peano/utils/Loop.h"
......@@ -33,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;
......@@ -77,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))
......@@ -116,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) {
......@@ -137,39 +289,33 @@ 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;
constexpr double amplitude = 1.5;
constexpr double v0[3] = {2., 1., 0.};
#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[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];
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);
/*
// 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);
}
}
......@@ -206,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];
......
......@@ -50,7 +50,7 @@ public:
* calls ::entropyWave
* ErrorWriter, ErrorPlotter write errors of numerical solution.
*/
static void referenceSolution(const double* const x, const double t, double* const Q);
inline static void referenceSolution(const double* const x, const double t, double* const Q) {entropyWave(x,t,Q);}
/**
* Adjust the conserved variables and parameters (together: Q) at a given time t at the (quadrature) point x.
......@@ -74,8 +74,12 @@ public:
* \param[inout] F the fluxes at that point as C array (already allocated).
*/
virtual void flux(const double* const Q,double** const F);
inline void flux2(const double* const Q,const double* const P,double** const F){flux(Q,F);} //add P array
void flux_vect(const double* const Q,const double* const P,double** const F); //add P array
virtual void nonConservativeProduct(const double* const Q, const double* const gradQ, double* const BgradQ);
void nonConservativeProduct2(const double* const Q, const double* const P, const double* const * const gradQ, double* const BgradQ); // add P array and gradQ like F
void nonConservativeProduct_vect(const double* const Q, const double* const P, const double* const * const gradQ, double* const BgradQ); // add P array and gradQ like F
/**
* Compute the eigenvalues of the flux tensor per coordinate direction \p d.
......
......@@ -115,9 +115,9 @@ optional arguments:
* --useMaterialParamVect enable vectorized material parameters
* --usePointSources numberOfPointSources enable numberOfPointSources point sources
* --useCERKGuess use CERK for SpaceTimePredictor inital guess (nonlinear only)
* --useSplitCKScalar use split Cauchy–Kowalevski formulation (linear only)
* --useSplitCKVect use split Cauchy–Kowalevski formulation with vect PDE (linear only)
* --useSplitCK use split Cauchy–Kowalevski formulation (linear only)
* --useGaussLobatto use Gauss Lobatto Quadrature instead of Gauss Legendre
* --useVectPDE use vectorized PDE terms (applies when present to: Flux, NCP, Source, FusedSource and MaterialParam)
* --tempVarsOnStack put the big scratch arrays on the stack instead of the
heap (you can use ulimit -s to increase the stack size)
......@@ -194,7 +194,7 @@ generated converter
Using the C index order convention with index in italic being absent in dim 2
Advanced options like ``--useSplitCKVect`` may use other data layout, see the
Advanced options like ``--useSplitCK`` or ``--useVectPDE`` may use other data layout, see the
respective templates
......
......@@ -62,21 +62,22 @@ class ArgumentParser:
("architecture", ArgType.MandatoryString, "the microarchitecture of the target device"),
# optional arguments
("useFlux", ArgType.OptionalBool, "enable flux"),
("useFluxVect", ArgType.OptionalBool, "enable vectorized flux (include useFlux)"),
("useFluxVect", ArgType.OptionalBool, "enable vectorized flux (include useFlux)"), #TODO JMG: legacy, use usePDEVect instead
("useViscousFlux", ArgType.OptionalBool, "enable viscous flux"),
("useNCP", ArgType.OptionalBool, "enable non conservative product"),
("useNCPVect", ArgType.OptionalBool, "enable vectorized non conservative product (include useNCP)"),
("useNCPVect", ArgType.OptionalBool, "enable vectorized non conservative product (include useNCP)"), #TODO JMG: legacy, use usePDEVect instead
("useSource", ArgType.OptionalBool, "enable source terms"),
("useSourceVect", ArgType.OptionalBool, "enable vectorized source terms (include useSource)"),
("useSourceVect", ArgType.OptionalBool, "enable vectorized source terms (include useSource)"), #TODO JMG: legacy, use usePDEVect instead
("useFusedSource", ArgType.OptionalBool, "enable fused source terms (include useSource)"),
("useFusedSourceVect", ArgType.OptionalBool, "enable vectorized fused source terms (include useFusedSource and useSourceVect)"),
("useFusedSourceVect", ArgType.OptionalBool, "enable vectorized fused source terms (include useFusedSource and useSourceVect)"), #TODO JMG: legacy, use usePDEVect instead
("useMaterialParam", ArgType.OptionalBool, "enable material parameters"),
("useMaterialParamVect",ArgType.OptionalBool, "enable vectorized material parameters"),
("useMaterialParamVect",ArgType.OptionalBool, "enable vectorized material parameters"), #TODO JMG: legacy, use usePDEVect instead
("usePointSources", ArgType.OptionalInt , "enable numberOfPointSources point sources", -1, "numberOfPointSources"),
("useCERKGuess", ArgType.OptionalBool, "use CERK for SpaceTimePredictor inital guess (nonlinear only)"),
("useSplitCKScalar", ArgType.OptionalBool, "use split Cauchy–Kowalevski formulation (linear only)"),
("useSplitCKVect", ArgType.OptionalBool, "use split Cauchy–Kowalevski formulation with vect PDE (linear only)"),
("useSplitCK", ArgType.OptionalBool, "use split Cauchy–Kowalevski formulation (linear only)"),
("useGaussLobatto", ArgType.OptionalBool, "use Gauss Lobatto Quadrature instead of Gauss Legendre"),
("predictorRecompute", ArgType.OptionalBool, "predictor step will recompute the PDE instead of relying on stored values from the picard loop (nonlinear only)"),
("useVectPDE", ArgType.OptionalBool, "use vectorized PDE terms (applies when present to: Flux, NCP, Source, FusedSource and MaterialParam)"),
("tempVarsOnStack", ArgType.OptionalBool, "put the big scratch arrays on the stack instead of the heap (you can use ulimit -s to increase the stack size)")
]
......
......@@ -95,8 +95,9 @@ class Controller:
"useMaterialParamVect" : args["useMaterialParamVect"],
"quadratureType" : ("Gauss-Lobatto" if args["useGaussLobatto"] else "Gauss-Legendre"),
"useCERKGuess" : args["useCERKGuess"],
"useSplitCKScalar" : args["useSplitCKScalar"],
"useSplitCKVect" : args["useSplitCKVect"]
"useSplitCK" : args["useSplitCK"],
"useVectPDE" : args["useVectPDE"],
"predictorRecompute" : args["predictorRecompute"]
})
self.config["useSourceOrNCP"] = self.config["useSource"] or self.config["useNCP"]
elif self.config["kernelType"] == "limiter":
......
......@@ -63,24 +63,25 @@ class ConfigurationParametersModel(AbstractModelBaseClass):
self.context["PSiSize"] = -1
self.context["PSiDerivativeSize"] = -1
if self.context["isLinear"]:
if(self.context["useSplitCKScalar"]):
# Linear + split CK
self.context["lQiSize"] = nVarPad*(nDof**nDim)
self.context["lQiNextSize"] = nVarPad*(nDof**nDim)
self.context["lPiSize"] = nParPad*(nDof**nDim)
self.context["lQhiSize"] = nVarPad*(nDof**nDim)
self.context["lFhiSize"] = nVarPad*(nDof**nDim)
self.context["gradQSize"] = nVarPad*(nDof**nDim)
self.context["PSiSize"] = nDof*(nDof**nDim)*nVarPad
elif(self.context["useSplitCKVect"]):
# Linear + split CK
self.context["lQiSize"] = nDof3D*nDof*nVar*nDofPad
self.context["lQiNextSize"] = nDof3D*nDof*nVar*nDofPad
self.context["lPiSize"] = nDof3D*nDof*nPar*nDofPad
self.context["lQhiSize"] = nDof3D*nDof*nVar*nDofPad
self.context["lFhiSize"] = nDof3D*nDof*nVar*nDofPad
self.context["gradQSize"] = nDof3D*nDof*nVar*nDofPad
self.context["PSiSize"] = nDof*nDof3D*nDof*nVar*nDofPad
if(self.context["useSplitCK"]):
if self.context["useVectPDE"]:
# Linear + split CK vect
self.context["lQiSize"] = nDof3D*nDof*nVar*nDofPad
self.context["lQiNextSize"] = nDof3D*nDof*nVar*nDofPad
self.context["lPiSize"] = nDof3D*nDof*nPar*nDofPad
self.context["lQhiSize"] = nDof3D*nDof*nVar*nDofPad
self.context["lFhiSize"] = nDof3D*nDof*nVar*nDofPad
self.context["gradQSize"] = nDof3D*nDof*nVar*nDofPad
self.context["PSiSize"] = nDof*nDof3D*nDof*nVar*nDofPad
else:
# Linear + split CK scalar
self.context["lQiSize"] = nVarPad*(nDof**nDim)
self.context["lQiNextSize"] = nVarPad*(nDof**nDim)
self.context["lPiSize"] = nParPad*(nDof**nDim)
self.context["lQhiSize"] = nVarPad*(nDof**nDim)
self.context["lFhiSize"] = nVarPad*(nDof**nDim)
self.context["gradQSize"] = nVarPad*(nDof**nDim)
self.context["PSiSize"] = nDof*(nDof**nDim)*nVarPad
else:
# default linear
self.context["lQiSize"] = nDataPad*(nDof**nDim)*(1+nDof)
......@@ -96,17 +97,41 @@ class ConfigurationParametersModel(AbstractModelBaseClass):
self.context["PSiSize"] = (nDof+1)*(nDof**nDim)*nVarPad
self.context["PSiDerivativeSize"] = self.context["PSiSize"]
else:
# nonlinear
self.context["lQiSize"] = nDataPad*(nDof**(nDim+1))
self.context["lQhiSize"] = nDataPad*(nDof**nDim)
if self.context["useFlux"]:
self.context["lFiSize"] = nVarPad*(nDof**(nDim+1))*nDim
self.context["lFhiSize"] = nVarPad*(nDof**nDim)*nDim
if self.context["useSource"] or self.context["useNCP"]:
self.context["lSiSize"] = nVarPad*(nDof**(nDim+1))
self.context["lShiSize"] = nVarPad*(nDof**nDim)
if self.context["useNCP"]:
self.context["gradQSize"] = nVarPad*(nDof**nDim)*nDim
# nonlinear
if self.context["predictorRecompute"]:
if self.context["useVectPDE"]:
self.context["lQiSize"] = nDofPad*nVar*(nDof**nDim)
self.context["lQhiSize"] = nDofPad*nVar*nDof*nDof3D
if nPar > 0:
self.context["lPiSize"] = nDofPad*nPar*nDof*nDof3D
if self.context["useFlux"]:
self.context["lFhiSize"] = nDofPad*nVar*nDof*nDof3D*nDim
if self.context["useSource"] or self.context["useNCP"]:
self.context["lShiSize"] = nDofPad*nVar*nDof*nDof3D
if self.context["useNCP"] or self.context["useViscousFlux"]:
self.context["gradQSize"] = nDofPad*nVar*nDof*nDof3D*nDim
else: #scalar predictorRecompute
self.context["lQiSize"] = nVarPad*(nDof**(nDim+1))
self.context["lQhiSize"] = nVarPad*(nDof**nDim)
if nPar > 0:
self.context["lPiSize"] = nParPad*(nDof**nDim)
if self.context["useFlux"]:
self.context["lFhiSize"] = nVarPad*(nDof**nDim)*nDim
if self.context["useSource"] or self.context["useNCP"]:
self.context["lShiSize"] = nVarPad*(nDof**nDim)
if self.context["useNCP"] or self.context["useViscousFlux"]:
self.context["gradQSize"] = nVarPad*(nDof**nDim)*nDim
else: # default nonlinear
self.context["lQiSize"] = nDataPad*(nDof**(nDim+1))
self.context["lQhiSize"] = nDataPad*(nDof**nDim)
if self.context["useFlux"]:
self.context["lFiSize"] = nVarPad*(nDof**(nDim+1))*nDim
self.context["lFhiSize"] = nVarPad*(nDof**nDim)*nDim
if self.context["useSource"] or self.context["useNCP"]:
self.context["lSiSize"] = nVarPad*(nDof**(nDim+1))
self.context["lShiSize"] = nVarPad*(nDof**nDim)
if self.context["useNCP"] or self.context["useViscousFlux"]:
self.context["gradQSize"] = nVarPad*(nDof**nDim)*nDim