Alsvinn  0.5.3
The fast FVM simulator with UQ support
Euler.hpp
Go to the documentation of this file.
1 /* Copyright (c) 2018 ETH Zurich, Kjetil Olsen Lye
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation, either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  */
15 
16 #pragma once
22 
24 #include "alsfvm/volume/Volume.hpp"
26 
28 #include <iostream>
29 namespace alsfvm {
30 namespace equation {
31 namespace euler {
32 
35 template<int nsd>
36 class Euler {
37 public:
38 
40  : gamma(parameters.getGamma()) {
41  }
42 
43  typedef typename Types<nsd>::rvec vec;
52 
56  static std::string getName();
57 
61  static const std::vector<std::string> conservedVariables;
62 
66  static const std::vector<std::string> primitiveVariables;
67 
71  static const std::vector<std::string> extraVariables;
72 
73 
77  static const size_t numberOfConservedVariables = 2 + nsd;
78 
80  return 2 + nsd;
81  }
82 
86 
91 
92 
97  size_t index) const {
98  return makeAllVariables(views.rho.at(index),
99  views.m(index).template convert<real>(),
100  views.E.at(index));
101  }
102 
103  template<class T, class S>
104  __device__ __host__ static ConservedVariables fetchConservedVariables(
105  euler::Views<T, S, nsd>& views, size_t index) {
106  return ConservedVariables(views.rho.at(index),
107  views.m(index),
108  views.E.at(index));
109  }
110 
112  size_t index) const {
113  return ExtraVariables(views.p.at(index),
114  views.u(index).template convert<real>()
115  );
116  }
117 
121  __device__ __host__ static void setViewAt(Views& output, size_t index,
122  const ConservedVariables& input) {
123  output.rho.at(index) = input.rho;
124  output.m(index) = input.m;
125  output.E.at(index) = input.E;
126 
127  }
128 
132  __device__ __host__ void setExtraViewAt(ViewsExtra& output, size_t index,
133  const ExtraVariables& input) const {
134  output.p.at(index) = input.p;
135  output.u(index) = input.u;
136  }
137 
143  __device__ __host__ void addToViewAt(Views& output, size_t index,
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;
148 
149  }
150 
171 
172  template<size_t direction>
173  __device__ __host__ void computePointFlux(const AllVariables& u,
174  ConservedVariables& F) const {
175  static_assert(direction < 3, "We only support up to three dimensions");
176 
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];
181  }
182 
187  template<size_t direction>
189  const ConservedVariables& u) const {
190  auto all = makeAllVariables(u.rho, u.m, u.E);
191  ConservedVariables F;
192  computePointFlux<direction>(all, F);
193  return F;
194  }
195 
196 
204  __device__ __host__ ExtraVariables computeExtra(const ConservedVariables& u)
205  const {
206  ExtraVariables v;
207  real ie = u.E - 0.5 * u.m.dot(u.m) / u.rho;
208  v.u = u.m / u.rho;
209  v.p = (gamma - 1) * ie;
210  return v;
211  }
212 
221  __device__ __host__ ExtraVariables computeExtra(const PrimitiveVariables&
222  primitiveVariables) const {
223  return ExtraVariables(primitiveVariables.p,
224  primitiveVariables.u
225  );
226  }
227 
237  const PrimitiveVariables& primitiveVariables) const {
238  const vec m = primitiveVariables.rho * primitiveVariables.u;
239 
240  const real E =
241  primitiveVariables.p / (gamma - 1)
242  + 0.5 * primitiveVariables.rho * primitiveVariables.u.dot(primitiveVariables.u);
243 
244  return ConservedVariables(primitiveVariables.rho, m, E);
245  }
246 
251  template<int direction>
252  __device__ __host__ real computeWaveSpeed(const ConservedVariables& u,
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");
256 
257  return fabs(v.u[direction]) + sqrt(gamma * v.p / u.rho);
258  }
259 
269  __device__ __host__ bool obeysConstraints(const ConservedVariables& u,
270  const ExtraVariables& v) const {
271 #ifdef __CUDA_ARCH__
272 
273  return u.rho < INFINITY && (u.rho == u.rho) && (u.rho > 0) && (v.p > 0);
274 #else
275  return std::isfinite(u.rho) && (u.rho == u.rho) && (u.rho > 0) && (v.p > 0);
276 #endif
277  }
278 
280  real E) const {
281 
282  ConservedVariables conserved(rho, m, E);
283  return AllVariables(conserved, computeExtra(conserved));
284  }
285 
286  __device__ __host__ real getWeight(const ConstViews& in, size_t index) const {
287  return in.rho.at(index);
288  }
289 
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;
294 
295  real p = (gamma - 1) * ie;
296  return PrimitiveVariables(conserved.rho, u, p);
297  }
298 
300  return gamma;
301  }
302 
312  const ConservedVariables& conserved) const {
313  PrimitiveVariables primitiveVariables = computePrimitiveVariables(conserved);
314 
315 
316  return TecnoVariables(real(sqrt(primitiveVariables.rho / primitiveVariables.p)),
317  real(sqrt(primitiveVariables.rho / primitiveVariables.p)) *
318  primitiveVariables.u,
319  real(sqrt(primitiveVariables.rho * primitiveVariables.p)));
320  }
321 
322 
339  const ConservedVariables& conserved) const {
340  auto primitive = computePrimitiveVariables(conserved);
341  const real s = log(primitive.p) - gamma * log(conserved.rho);
342 
343  return state_vector((gamma - s) / (gamma - 1) - (conserved.rho *
344  (primitive.u.dot(primitive.u))) / (2 * primitive.p),
345  conserved.m / primitive.p,
346  -conserved.rho / primitive.p);
347 
348  }
349 
358  __device__ __host__ vec computeEntropyPotential(const ConservedVariables&
359  conserved) const {
360  return conserved.m;
361  }
362 
363 
364  template<int direction>
367  const ConservedVariables&
368  conserved) const {
369 
370  return computeEigenVectorMatrix<direction>(conserved).transposed() *
371  computeEntropyVariables(conserved);
372  }
373 
374 
376  template<int direction>
378  const ConservedVariables& conserved) const;
379 
382  template<int direction>
383  __device__ __host__ state_vector computeEigenValues(const ConservedVariables&
384  conserved) const;
385 private:
386  const real gamma;
387 };
388 
389 
390 template<>
391 inline std::string Euler<3>::getName() {
392  return "euler3";
393 }
394 
395 template<>
396 inline std::string Euler<2>::getName() {
397  return "euler2";
398 }
399 
400 template<>
401 inline std::string Euler<1>::getName() {
402  return "euler1";
403 }
404 
405 
406 
407 template<>
408 template<int direction>
410  const ConservedVariables& conservedConst) const {
411  auto conserved = conservedConst;
412 
413  if (direction == 1) {
414  // We use the rotation trick, see 3.2.2 and Proposition 3.19 in Toro's book
415  // http://www.springer.com/de/book/9783540252023
416  conserved = ConservedVariables(conserved.rho, rvec3{conserved.m.y, -conserved.m.x, conserved.m.z},
417  conserved.E);
418 
419  } else if (direction == 2) {
420  // We use the rotation trick, see 3.2.2 and Proposition 3.19 in Toro's book
421  // http://www.springer.com/de/book/9783540252023
422  conserved = ConservedVariables(conserved.rho, rvec3{conserved.m.z, conserved.m.y, -conserved.m.x},
423  conserved.E);
424  }
425 
426 
427  if (direction < 3) {
428  matrix5 matrixWithEigenVectors;
429 
430  auto primitive = computePrimitiveVariables(conserved);
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;
438 
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;
444 
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;
450 
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;
456 
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;
462 
463  return matrixWithEigenVectors.normalized();
464  }
465 
466 
467  assert(false);
468 }
469 
470 template<>
471 template<int direction>
473  conservedConst) const {
474  auto conserved = conservedConst;
475 
476  if (direction == 1) {
477  // We use the rotation trick, see 3.2.2 and Proposition 3.19 in Toro's book
478  // http://www.springer.com/de/book/9783540252023
479  conserved = ConservedVariables(conserved.rho, rvec3{conserved.m.y, -conserved.m.x, conserved.m.z},
480  conserved.E);
481 
482  } else if (direction == 2) {
483  // We use the rotation trick, see 3.2.2 and Proposition 3.19 in Toro's book
484  // http://www.springer.com/de/book/9783540252023
485  conserved = ConservedVariables(conserved.rho, rvec3{conserved.m.z, conserved.m.y, -conserved.m.x},
486  conserved.E);
487  }
488 
489  if (direction < 3) {
490  auto primitive = computePrimitiveVariables(conserved);
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,
493  primitive.u.x + a);
494  }
495 
496 
497  assert(false);
498 
499  return rvec5{0, 0, 0};
500 }
501 
502 
503 
504 template<>
505 template<int direction>
507  const ConservedVariables& conservedConst) const {
508 
509 
510  ConservedVariables conserved = conservedConst;
511 
512  if (direction == 1) {
513  // We use the rotation trick, see 3.2.2 and Proposition 3.19 in Toro's book
514  // http://www.springer.com/de/book/9783540252023
515  conserved = ConservedVariables(conserved.rho, rvec2{conserved.m.y, -conserved.m.x},
516  conserved.E);
517  }
518 
519  matrix4 matrixWithEigenVectors;
520 
521  if (direction < 2) {
522 
523 
524  auto primitive = computePrimitiveVariables(conserved);
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;
531 
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;
536 
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;
541 
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;
546 
547  return matrixWithEigenVectors;
548  }
549 
550  assert(false);
551 
552  return matrixWithEigenVectors;
553 }
554 
555 template<>
556 template<int direction>
558  conservedConst) const {
559  ConservedVariables conserved = conservedConst;
560 
561  if (direction == 1) {
562  // We use the rotation trick, see 3.2.2 and Proposition 3.19 in Toro's book
563  // http://www.springer.com/de/book/9783540252023
564  conserved = ConservedVariables(conserved.rho, rvec2{conserved.m.y, -conserved.m.x},
565  conserved.E);
566  }
567 
568  if (direction < 2) {
569  auto primitive = computePrimitiveVariables(conserved);
570  const real a = sqrt(gamma * primitive.p / conserved.rho);
571  return rvec4(primitive.u.x - a, primitive.u.x, primitive.u.x,
572  primitive.u.x + a);
573  }
574 
575 
576 
577  assert(false);
578 
579  return rvec4{0, 0, 0};
580 }
581 
582 
583 
584 template<>
585 template<int direction>
587  const ConservedVariables& conserved) const {
588 
589  matrix3 matrixWithEigenVectors;
590 
591  if (direction == 0) {
592 
593 #if 1
594  const real rho = conserved.rho;
595 
596  const auto primitive = computePrimitiveVariables(conserved);
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);
601  const real a4 = a1;
602 
603 
604 
605  matrixWithEigenVectors(0, 0) = a1;
606  matrixWithEigenVectors(0, 1) = a2;
607  matrixWithEigenVectors(0, 2) = a4;
608 
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);
612 
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);
616 
617 
618  return matrixWithEigenVectors;
619 
620 #else
621 
622 
623  auto primitive = computePrimitiveVariables(conserved);
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;
629 
630  matrixWithEigenVectors(1, 0) = primitive.u.x - a;
631  matrixWithEigenVectors(1, 1) = primitive.u.x;
632  matrixWithEigenVectors(1, 2) = primitive.u.x + a;
633 
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();
638 #endif
639 
640  }
641 
642  assert(false);
643  return matrixWithEigenVectors;
644 }
645 
646 template<>
647 template<int direction>
649  conserved) const {
650  if (direction == 0) {
651  auto primitive = computePrimitiveVariables(conserved);
652  const real a = sqrt(gamma * primitive.p / conserved.rho);
653  return rvec3(primitive.u.x - a, primitive.u.x, primitive.u.x + a);
654  }
655 
656  assert(false);
657 
658  return rvec3{0, 0, 0};
659 }
660 
661 
662 
663 }
664 } // namespace alsfvm
665 } // namespace equation
666 
__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: vec4.hpp:25
Definition: types.hpp:104
__device__ static __host__ void setViewAt(Views &output, size_t index, const ConservedVariables &input)
Definition: Euler.hpp:121
Definition: Views.hpp:33
__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
Definition: vec2.hpp:24
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: mat.hpp:23
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
rvec u
Definition: ExtraVariables.hpp:59
#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
Definition: vec5.hpp:25
static const size_t numberOfConservedVariables
Definition: Euler.hpp:77
Definition: Euler.hpp:36
__device__ __host__ state_vector computeEigenValues(const ConservedVariables &conserved) const
real p
Definition: ExtraVariables.hpp:58
real rho
Definition: ConservedVariables.hpp:92
Definition: AllVariables.hpp:27
Euler(const EulerParameters &parameters)
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
Definition: ViewsExtra.hpp:31
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
Definition: ExtraVariables.hpp:29
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