99 views.m(index).template convert<real>(),
103 template<
class T,
class S>
112 size_t index)
const {
114 views.u(index).template convert<real>()
122 const ConservedVariables& input) {
123 output.rho.at(index) = input.
rho;
124 output.m(index) = input.
m;
125 output.E.at(index) = input.
E;
133 const ExtraVariables& input)
const {
134 output.p.at(index) = input.
p;
135 output.u(index) = input.
u;
144 const ConservedVariables& input)
const {
145 output.rho.at(index) += input.
rho;
146 output.m(index) += input.
m;
147 output.E.at(index) += input.
E;
172 template<
size_t direction>
174 ConservedVariables& F)
const {
175 static_assert(direction < 3,
"We only support up to three dimensions");
177 F.
rho = u.
m[direction];
178 F.
m = u.
u[direction] * u.
m;
179 F.
m[direction] += u.
p;
180 F.
E = (u.
E + u.
p) * u.
u[direction];
187 template<
size_t direction>
189 const ConservedVariables& u)
const {
191 ConservedVariables F;
192 computePointFlux<direction>(all, F);
209 v.
p = (gamma - 1) * ie;
222 primitiveVariables)
const {
237 const PrimitiveVariables& primitiveVariables)
const {
238 const vec m = primitiveVariables.
rho * primitiveVariables.
u;
241 primitiveVariables.
p / (gamma - 1)
242 + 0.5 * primitiveVariables.
rho * primitiveVariables.
u.dot(primitiveVariables.
u);
251 template<
int direction>
253 const ExtraVariables& v)
const {
254 static_assert(direction >= 0,
"Direction can not be negative");
255 static_assert(direction < 3,
"We only support dimension up to and inclusive 3");
257 return fabs(v.
u[direction]) + sqrt(gamma * v.
p / u.
rho);
270 const ExtraVariables& v)
const {
273 return u.
rho < INFINITY && (u.
rho == u.
rho) && (u.
rho > 0) && (v.
p > 0);
275 return std::isfinite(u.
rho) && (u.
rho == u.
rho) && (u.
rho > 0) && (v.
p > 0);
282 ConservedVariables conserved(rho, m, E);
287 return in.rho.at(index);
291 const ConservedVariables& conserved)
const {
292 vec u = conserved.
m / conserved.
rho;
293 real ie = conserved.
E - 0.5 * conserved.
m.dot(conserved.
m) / conserved.
rho;
295 real p = (gamma - 1) * ie;
312 const ConservedVariables& conserved)
const {
317 real(sqrt(primitiveVariables.
rho / primitiveVariables.
p)) *
318 primitiveVariables.
u,
319 real(sqrt(primitiveVariables.
rho * primitiveVariables.
p)));
339 const ConservedVariables& conserved)
const {
341 const real s = log(primitive.p) - gamma * log(conserved.
rho);
344 (primitive.u.dot(primitive.u))) / (2 * primitive.p),
345 conserved.
m / primitive.p,
346 -conserved.
rho / primitive.p);
364 template<
int direction>
367 const ConservedVariables&
370 return computeEigenVectorMatrix<direction>(conserved).transposed() *
376 template<
int direction>
378 const ConservedVariables& conserved)
const;
382 template<
int direction>
408 template<
int direction>
411 auto conserved = conservedConst;
413 if (direction == 1) {
419 }
else if (direction == 2) {
428 matrix5 matrixWithEigenVectors;
431 const real a = sqrt(gamma * primitive.p / conserved.rho);
432 const real H = (conserved.E + primitive.p) / conserved.rho;
433 matrixWithEigenVectors(0, 0) = 1;
434 matrixWithEigenVectors(0, 1) = 1;
435 matrixWithEigenVectors(0, 2) = 0;
436 matrixWithEigenVectors(0, 3) = 0;
437 matrixWithEigenVectors(0, 4) = 1;
439 matrixWithEigenVectors(1, 0) = primitive.u.x - a;
440 matrixWithEigenVectors(1, 1) = primitive.u.x;
441 matrixWithEigenVectors(1, 2) = 0;
442 matrixWithEigenVectors(1, 3) = 0;
443 matrixWithEigenVectors(1, 4) = primitive.u.x + a;
445 matrixWithEigenVectors(2, 0) = primitive.u.y;
446 matrixWithEigenVectors(2, 1) = primitive.u.y;
447 matrixWithEigenVectors(2, 2) = 1;
448 matrixWithEigenVectors(2, 3) = 0;
449 matrixWithEigenVectors(2, 4) = primitive.u.y;
451 matrixWithEigenVectors(3, 0) = primitive.u.z;
452 matrixWithEigenVectors(3, 1) = primitive.u.z;
453 matrixWithEigenVectors(3, 2) = 0;
454 matrixWithEigenVectors(3, 3) = 1;
455 matrixWithEigenVectors(3, 4) = primitive.u.z;
457 matrixWithEigenVectors(4, 0) = H - primitive.u.x * a;
458 matrixWithEigenVectors(4, 1) = 0.5 * primitive.u.dot(primitive.u);
459 matrixWithEigenVectors(4, 2) = primitive.u.y;
460 matrixWithEigenVectors(4, 3) = primitive.u.z;
461 matrixWithEigenVectors(4, 4) = H + primitive.u.x * a;
471 template<
int direction>
473 conservedConst)
const {
474 auto conserved = conservedConst;
476 if (direction == 1) {
482 }
else if (direction == 2) {
491 const real a = sqrt(gamma * primitive.p / conserved.rho);
492 return rvec5(primitive.u.x - a, primitive.u.x, primitive.u.x, primitive.u.x,
499 return rvec5{0, 0, 0};
505 template<
int direction>
512 if (direction == 1) {
519 matrix4 matrixWithEigenVectors;
525 const real a = sqrt(gamma * primitive.p / conserved.
rho);
526 const real H = (conserved.
E + primitive.p) / conserved.
rho;
527 matrixWithEigenVectors(0, 0) = 1;
528 matrixWithEigenVectors(0, 1) = 1;
529 matrixWithEigenVectors(0, 2) = 0;
530 matrixWithEigenVectors(0, 3) = 1;
532 matrixWithEigenVectors(1, 0) = primitive.u.x - a;
533 matrixWithEigenVectors(1, 1) = primitive.u.x;
534 matrixWithEigenVectors(1, 2) = 0;
535 matrixWithEigenVectors(1, 3) = primitive.u.x + a;
537 matrixWithEigenVectors(2, 0) = primitive.u.y;
538 matrixWithEigenVectors(2, 1) = primitive.u.y;
539 matrixWithEigenVectors(2, 2) = 1;
540 matrixWithEigenVectors(2, 3) = primitive.u.y;
542 matrixWithEigenVectors(3, 0) = H - primitive.u.x * a;
543 matrixWithEigenVectors(3, 1) = 0.5 * primitive.u.dot(primitive.u);
544 matrixWithEigenVectors(3, 2) = primitive.u.y;
545 matrixWithEigenVectors(3, 3) = H + primitive.u.x * a;
547 return matrixWithEigenVectors;
552 return matrixWithEigenVectors;
556 template<
int direction>
558 conservedConst)
const {
561 if (direction == 1) {
570 const real a = sqrt(gamma * primitive.p / conserved.
rho);
571 return rvec4(primitive.u.x - a, primitive.u.x, primitive.u.x,
579 return rvec4{0, 0, 0};
585 template<
int direction>
589 matrix3 matrixWithEigenVectors;
591 if (direction == 0) {
597 const real a = sqrt(gamma * primitive.p / conserved.
rho);
598 const real H = (conserved.
E + primitive.p) / conserved.
rho;
599 const real a1 = sqrt(rho / (2 * gamma));
600 const real a2 = sqrt(rho * (gamma - 1) / gamma);
605 matrixWithEigenVectors(0, 0) = a1;
606 matrixWithEigenVectors(0, 1) = a2;
607 matrixWithEigenVectors(0, 2) = a4;
609 matrixWithEigenVectors(1, 0) = a1 * (primitive.u.x - a);
610 matrixWithEigenVectors(1, 1) = a2 * primitive.u.x;
611 matrixWithEigenVectors(1, 2) = a4 * (primitive.u.x + a);
613 matrixWithEigenVectors(2, 0) = a1 * (H - primitive.u.x * a);
614 matrixWithEigenVectors(2, 1) = 0.5 * a2 * (primitive.u.dot(primitive.u));
615 matrixWithEigenVectors(2, 2) = a4 * (H + primitive.u.x * a);
618 return matrixWithEigenVectors;
624 const real a = sqrt(gamma * primitive.p / conserved.
rho);
625 const real H = (conserved.
E + primitive.p) / conserved.
rho;
626 matrixWithEigenVectors(0, 0) = 1;
627 matrixWithEigenVectors(0, 1) = 1;
628 matrixWithEigenVectors(0, 2) = 1;
630 matrixWithEigenVectors(1, 0) = primitive.u.x - a;
631 matrixWithEigenVectors(1, 1) = primitive.u.x;
632 matrixWithEigenVectors(1, 2) = primitive.u.x + a;
634 matrixWithEigenVectors(2, 0) = H - primitive.u.x * a;
635 matrixWithEigenVectors(2, 1) = 0.5 * primitive.u.dot(primitive.u);
636 matrixWithEigenVectors(2, 2) = H + primitive.u.x * a;
637 return matrixWithEigenVectors.normalized();
643 return matrixWithEigenVectors;
647 template<
int direction>
650 if (direction == 0) {
652 const real a = sqrt(gamma * primitive.p / conserved.
rho);
653 return rvec3(primitive.u.x - a, primitive.u.x, primitive.u.x + a);
658 return rvec3{0, 0, 0};
__device__ __host__ state_matrix computeEigenVectorMatrix(const ConservedVariables &conserved) const
Computes the Eigen vector matrix. See 3.2.2 for full description in http://www.springer.com/de/book/9783540252023.
real rho
rho is the density
Definition: PrimitiveVariables.hpp:58
__device__ __host__ PrimitiveVariables computePrimitiveVariables(const ConservedVariables &conserved) const
Definition: Euler.hpp:290
Definition: types.hpp:104
__device__ static __host__ void setViewAt(Views &output, size_t index, const ConservedVariables &input)
Definition: Euler.hpp:121
__device__ __host__ AllVariables fetchAllVariables(ConstViews &views, size_t index) const
Definition: Euler.hpp:96
__device__ __host__ ConservedVariables computePointFlux(const ConservedVariables &u) const
Definition: Euler.hpp:188
__device__ __host__ real getGamma() const
Definition: Euler.hpp:299
Types< nsd+2 >::matrix state_matrix
Definition: Euler.hpp:45
__device__ __host__ ConservedVariables computeConserved(const PrimitiveVariables &primitiveVariables) const
computes the extra variables from the primitive ones
Definition: Euler.hpp:236
simulator::SimulatorParameters & parameters
Definition: CellComputerFactory.cpp:60
__device__ __host__ void computePointFlux(const AllVariables &u, ConservedVariables &F) const
Definition: Euler.hpp:173
__device__ __host__ ExtraVariables fetchExtraVariables(ConstViewsExtra &views, size_t index) const
Definition: Euler.hpp:111
__device__ __host__ real computeWaveSpeed(const ConservedVariables &u, const ExtraVariables &v) const
Definition: Euler.hpp:252
#define __host__
Definition: types.hpp:46
__device__ __host__ ExtraVariables computeExtra(const ConservedVariables &u) const
Definition: Euler.hpp:204
__device__ __host__ void setExtraViewAt(ViewsExtra &output, size_t index, const ExtraVariables &input) const
Definition: Euler.hpp:132
static const std::vector< std::string > extraVariables
Definition: Euler.hpp:71
Definition: ConservedVariables.hpp:30
double real
Definition: types.hpp:65
int rho
Definition: sodshocktube.py:3
equation::euler::Views< const volume::Volume, const memory::View< const real >, nsd > ConstViews
Definition: Euler.hpp:85
__device__ __host__ real getWeight(const ConstViews &in, size_t index) const
Definition: Euler.hpp:286
Types< nsd+2 >::rvec state_vector
Definition: Euler.hpp:44
Definition: TecnoVariables.hpp:34
equation::euler::Views< volume::Volume, memory::View< real >, nsd > Views
Definition: Euler.hpp:83
__device__ __host__ vec computeEntropyPotential(const ConservedVariables &conserved) const
Definition: Euler.hpp:358
__device__ __host__ TecnoVariables computeTecnoVariables(const ConservedVariables &conserved) const
Definition: Euler.hpp:311
__device__ __host__ state_vector computeEntropyVariablesMultipliedByEigenVectorMatrix(const ConservedVariables &conserved) const
Definition: Euler.hpp:366
Definition: PrimitiveVariables.hpp:28
__device__ __host__ matrix< T, NumberOfColumns, NumberOfRows > normalized() const
Definition: mat.hpp:103
rvec m
Definition: ConservedVariables.hpp:93
euler::PrimitiveVariables< nsd > PrimitiveVariables
Definition: Euler.hpp:49
euler::EulerParameters Parameters
Definition: Euler.hpp:46
euler::ConservedVariables< nsd > ConservedVariables
Definition: Euler.hpp:47
equation::euler::ViewsExtra< volume::Volume, memory::View< real >, nsd > ViewsExtra
Definition: Euler.hpp:88
vec4< real > rvec4
Definition: types.hpp:79
__device__ __host__ void addToViewAt(Views &output, size_t index, const ConservedVariables &input) const
Definition: Euler.hpp:143
__device__ static __host__ size_t getNumberOfConservedVariables()
Definition: Euler.hpp:79
__device__ __host__ state_vector computeEntropyVariables(const ConservedVariables &conserved) const
Definition: Euler.hpp:338
#define static_assert(x, y)
Definition: types.hpp:52
__device__ __host__ bool obeysConstraints(const ConservedVariables &u, const ExtraVariables &v) const
Definition: Euler.hpp:269
rvec u
u is the velocity
Definition: PrimitiveVariables.hpp:63
vec5< real > rvec5
Definition: types.hpp:82
__device__ __host__ AllVariables makeAllVariables(real rho, vec m, real E) const
Definition: Euler.hpp:279
__device__ __host__ ExtraVariables computeExtra(const PrimitiveVariables &primitiveVariables) const
computes the extra variables from the primitive ones
Definition: Euler.hpp:221
static std::string getName()
static const std::vector< std::string > conservedVariables
Definition: Euler.hpp:61
#define __device__
Definition: types.hpp:45
static const size_t numberOfConservedVariables
Definition: Euler.hpp:77
__device__ __host__ state_vector computeEigenValues(const ConservedVariables &conserved) const
real rho
Definition: ConservedVariables.hpp:92
Definition: AllVariables.hpp:27
Euler(const EulerParameters ¶meters)
Definition: Euler.hpp:39
vec3< real > rvec3
Definition: types.hpp:76
Various utility functions to implement the tecno flux.
Definition: types.hpp:30
float p
Definition: sodshocktube.py:5
euler::TecnoVariables< nsd > TecnoVariables
Definition: Euler.hpp:51
static const std::vector< std::string > primitiveVariables
Definition: Euler.hpp:66
equation::euler::ViewsExtra< const volume::Volume, const memory::View< const real >, nsd > ConstViewsExtra
Definition: Euler.hpp:90
real p
p is the pressure
Definition: PrimitiveVariables.hpp:68
__device__ static __host__ ConservedVariables fetchConservedVariables(euler::Views< T, S, nsd > &views, size_t index)
Definition: Euler.hpp:104
euler::ExtraVariables< nsd > ExtraVariables
Definition: Euler.hpp:48
Types< nsd >::rvec vec
Definition: Euler.hpp:43
Definition: EulerParameters.hpp:24
real E
Definition: ConservedVariables.hpp:94
euler::AllVariables< nsd > AllVariables
Definition: Euler.hpp:50