Peano 4
Loading...
Searching...
No Matches
RiemannSolverUtils.h
Go to the documentation of this file.
1
14#pragma once
15
16#include "tarch/Assertions.h"
17#include "tarch/la/Vector.h"
19
20#include <sstream>
21
23
24namespace kernels {
25namespace riemannsolvers {
26namespace util {
27
54template <bool useNCP,bool useFlux,int numberOfVariables, int numberOfParameters, int numQuadPoints, typename SolverType>
56 SolverType& solver,
57 const double* const qL,
58 const double* const qR,
59 const int direction,
60 double (&osherMatrix)[numberOfVariables][numberOfVariables] /*initialise with zeroes*/,
61 double (&osherNCP)[numberOfVariables] /*initialise with zeroes*/) {
62 // coordinate axes permutations (if direction is template parameter, this can be constexpr)
63 // original; keep a while for reference
64 // int in = 0; // direction == 0
65 // int is = 1;
66 // int it = 2;
67 // if ( direction == 1 ) {
68 // in=1;
69 // is=0;
70 // it=2;
71 // } else if ( direction == 2 ) {
72 // in=2;
73 // is=0;
74 // it=1;
75 // }
76 constexpr int numberOfData = numberOfVariables + numberOfParameters;
77
78 const int in = direction;
79 const int is = (direction == 0) ? 1 : 0;
80 const int it = (direction != 2) ? 2 : 1;
81
82 // integrate matrix |A| := Rs * |L| * iRs matrix over phase space (QL -> QR)
83 double s_max = -1.0;
84 for(int i=0; i<numQuadPoints; i++) {
85 const double& si = kernels::legendre::nodes[numQuadPoints-1][i];
86 double Qs[numberOfData];
87 for(int j=0; j < numberOfData; j++) {
88 Qs[j] = qL[j] + si * (qR[j] - qL[j]);
89 }
90
91 if (useFlux) {
92 // eigenstructure
93 double Ls[numberOfVariables]; // all eigenvalues need to be set by user
94 double Rs[numberOfVariables][numberOfVariables] = {0.0}; // fill with zeros as user typically performs rotation via matrix product
95 double iRs[numberOfVariables][numberOfVariables]= {0.0};
96 solver.eigenvectors(Qs/*in*/,in,is,it,Rs/*inout*/,Ls/*inout*/,iRs/*inout*/); // user solver can work internally with eigen system or eigen solvers;
97 for (int j = 0; j < numberOfVariables; j++) {
98 s_max = std::max( std::abs(Ls[j]), s_max );
99 }
100
101 #ifdef Asserts
102 double sumOfEntries = 0;
103 for(int j=0; j < numberOfVariables; j++) {
104 for(int k=0; k < numberOfVariables; k++) {
105 for(int a=0; a < numberOfVariables; a++) {
106 sumOfEntries += iRs[j][a] * Rs[a][k];
107 }
108 }
109 }
110 if ( sumOfEntries < numberOfVariables - 1e-12 || sumOfEntries > numberOfVariables + 1e-12 ) {
111 std::cerr << "Error: Left eigenvector matrix is not an inverse of the right one as sum of entries="<<sumOfEntries << "; direction="<<direction << std::endl;
112 std::terminate();
113 }
114 #endif
115
116 // compute Osher matrix
117 // scale column vectors: (|L1|*col1,|L2|*col2,...)
118 const double& wi = kernels::legendre::weights[numQuadPoints-1][i];
119 for(int j=0; j < numberOfVariables; j++) {
120 for(int k=0; k < numberOfVariables; k++) {
121 for(int a=0; a < numberOfVariables; a++) {
122 osherMatrix[j][k] += wi * Rs[j][a]*std::abs(Ls[a])*iRs[a][k];
123 }
124 }
125 }
126 }
127
128 if (useNCP) { // we directly compute the nonconservative product instead of computing the jacobian matrix
129 const double& wi = kernels::legendre::weights[numQuadPoints-1][i];
130 double gradQs[Dimensions][numberOfData] = {0.0};
131 for(int j=0; j < numberOfData; j++) {
132 gradQs[Dimensions][j] = qR[j] - qL[j];
133 }
134 double ncp[numberOfData] = {0.0};
135 solver.nonConservativeProduct(Qs, gradQs[in], ncp);
136 for(int j=0; j < numberOfData; j++) {
137 osherNCP[j] += wi * ncp[j];
138 }
139 }
140 }
141 return s_max;
142}
143
151template <int basisSize, int numberOfData>
153 const double* const Q,
154 const double* const * const weights1D,
155 double (&Qav)[numberOfData]) {
156 constexpr int order = basisSize-1;
157 #if Dimensions==2
158 idx2 idx_Q(basisSize, numberOfData);
159 for (int j = 0; j < basisSize; j++) {
160 const double wj = weights1D[order][j];
161 for (int k = 0; k < numberOfData; k++) {
162 Qav[k] += wj * Q[idx_Q(j, k)];
163 }
164 }
165 #else
166 idx3 idx_Q(basisSize, basisSize, numberOfData);
167 for (int i = 0; i < basisSize; i++) {
168 for (int j = 0; j < basisSize; j++) {
169 const double wij = weights1D[order][i]*weights1D[order][j];
170 for (int k = 0; k < numberOfData; k++) {
171 Qav[k] += wij * Q[idx_Q(i,j, k)];
172 }
173 }
174 }
175 #endif
176}
177
178
179} // namespace util
180} // namespace riemannsolvers
181} // namespace kernels
182
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$ our matrix elements will nabla phi_i dx f By this will be a sparse as these basis functions are chosen to not overlap with each other almost everywhere In other they have only local support We can read off the right hand side by taking our known right hand side f$ f f$ and integrating against an appropriate test phi_i dx f Please excuse the slight abuse of notation here There should probably be a clearer indication that we move from a continuous f$ f f$ to some discrete f$ f_i f$ We can demonstrate simply It s worth as when we discuss the discontinuous version of this it will no longer disappear We take our left hand side and discretise it
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
double computeOsherMatrix(SolverType &solver, const double *const qL, const double *const qR, const int direction, double(&osherMatrix)[numberOfVariables][numberOfVariables], double(&osherNCP)[numberOfVariables])
Computes ingredients for the generalised Osher Solomon flux [1,2].
void averageRiemannInputs(const double *const Q, const double *const *const weights1D, double(&Qav)[numberOfData])
Average Riemann input data.
This file is part of the ExaHyPE project.