Peano 4
Loading...
Searching...
No Matches
riemannSolverNonlinear.cpph
Go to the documentation of this file.
1
14// included in ../../Kernels.h
15
16#include <algorithm>
17#include <cassert>
18#include <cmath>
19#include <cstring>
20
22
23namespace kernels {
24namespace aderdg {
25namespace generic {
26namespace c {
27
37template <bool useNCP, typename SolverType>
38void riemannSolverNonlinear(SolverType& solver, double* FL, double* FR,
39 const double* const QL,
40 const double* const QR,
41 const double t,
42 const double dt,
45 const int direction) {
46 constexpr int numberOfVariables = SolverType::NumberOfVariables;
47 constexpr int numberOfParameters = SolverType::NumberOfParameters;
48 constexpr int numberOfData = numberOfVariables+numberOfParameters;
49 constexpr int order = SolverType::Order;
50 constexpr int basisSize = order+1;
51
52 // Compute the average variables and parameters from the left and the right
53 double QavL[numberOfData] = {0.0}; // ~(numberOfVariables+numberOfParameters)
54 double QavR[numberOfData] = {0.0}; // ~(numberOfVariables+numberOfParameters)
55 kernels::riemannsolvers::util::averageRiemannInputs<basisSize,numberOfData>(
57 kernels::riemannsolvers::util::averageRiemannInputs<basisSize,numberOfData>(
59
60 // Hyperbolic eigenvalues
61 const double maxHyperbolicEigenvalueL = solver.maxEigenvalue(QavL, faceCentre, dx, t, dt, direction);
62 const double maxHyperbolicEigenvalueR = solver.maxEigenvalue(QavR, faceCentre, dx, t, dt, direction);
63 const double smax = std::max(maxHyperbolicEigenvalueL, maxHyperbolicEigenvalueR);
64
65 // compute fluxes (and fluctuations for non-conservative PDEs)
66 double Qavg[numberOfData];
67 idx2 idx_gradQ(Dimensions, numberOfVariables);
68 double gradQ[numberOfVariables] = {0.0};
69 double ncp[numberOfVariables] = {0.0};
70 {
71 idx2 idx_FLR(basisSize, numberOfVariables);
72 idx2 idx_QLR(basisSize, numberOfData);
73
74 for (int j = 0; j < basisSize; j++) {
75
76 if(useNCP) { // we don't use matrixB but the NCP call here.
77 for(int l=0; l < numberOfVariables; l++) {
78 gradQ[l] = QR[idx_QLR(j, l)] - QL[idx_QLR(j, l)];
79 }
80 for(int l=0; l < numberOfData; l++) {
81 Qavg[l] = 0.5 * (QR[idx_QLR(j, l)] + QL[idx_QLR(j, l)]);
82 }
83
84 solver.nonconservativeProduct(
85 Qavg,
86 gradQ,
87 faceCentre,
88 dx,
89 t,
90 dt,
91 direction,
92 ncp
93 );
94 }
95
96 // skip parameters
97 for (int k = 0; k < numberOfVariables; k++) {
98 FL[idx_FLR(j, k)] =
99 0.5 * (FR[idx_FLR(j, k)] + FL[idx_FLR(j, k)]) -
100 0.5 * smax * (QR[idx_QLR(j, k)] - QL[idx_QLR(j, k)]);
101 //assertion4(std::isfinite(FL[idx_FLR(j, k)]), j, k, smax,
102 // maxHyperbolicEigenvalue);
103
104 if(useNCP) {
105 FR[idx_FLR(j, k)] = FL[idx_FLR(j, k)] - 0.5 * ncp[k];
106 FL[idx_FLR(j, k)] = FL[idx_FLR(j, k)] + 0.5 * ncp[k];
107 } else {
108 FR[idx_FLR(j, k)] = FL[idx_FLR(j, k)];
109 }
110 }
111 }
112 }
113}
114
115} // namespace c
116} // namespace generic
117} // namespace aderdg
118} // namespace kernels
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ so we are integrating with f$ phi_k phi_k dx
-with-ExaHyPE2-benchmark
void riemannSolverNonlinear(SolverType &solver, double *FL, double *FR, const double *const QL, const double *const QR, const double t, const double dt, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &dx, const int direction)
We implement a very simple Rusanov scheme with scalar dissipation (smax*Id).
This file is part of the ExaHyPE project.
Simple vector class.
Definition Vector.h:134