Peano 4
Loading...
Searching...
No Matches
ADERDGMGCCZ4.cpp
Go to the documentation of this file.
1#include "ADERDGMGCCZ4.h"
2#include "InitialValues.h"
4#include "Constants.h"
5#include "MGCCZ4Kernels.h"
6#include <algorithm>
7#include <limits>
8
9#include <stdio.h>
10#include <string.h>
12
13
14tarch::logging::Log examples::exahype2::mgccz4::ADERDGMGCCZ4::_log( "examples::exahype2::ADERDGMGCCZ4::MGCCZ4" );
15
16
17
19 if ( Scenario=="gaugewave-c++" || Scenario=="linearwave-c++" ) {
20 const char* name = "GaugeWave";
21 int length = strlen(name);
22 //initparameters_(&length, name);
23 }
24 else {
25 std::cerr << "initial scenario " << Scenario << " is not supported" << std::endl << std::endl << std::endl;
26 }
27}
28
30 double * __restrict__ Q,
32 double t
33) {
34 logTraceInWith2Arguments( "adjustSolution(...)", x, t);
35 if (tarch::la::equals(t,0.0) ) {
36 if ( Scenario=="gaugewave-c++" ) {
38 }
39 else if ( Scenario=="linearwave-c++" ) {
41 }
42 else {
43 logError( "adjustSolution(...)", "initial scenario " << Scenario << " is not supported" );
44 }
45
46 for (int i=0; i<NumberOfUnknowns; i++) {
47 assertion3( std::isfinite(Q[i]), x, t, i );
48 }
49
50 for (int i=NumberOfUnknowns; i<NumberOfUnknowns+NumberOfAuxiliaryVariables; i++) {
51 Q[i] = 0.0;
52 }
53 }
54 else {
56 }
57 logTraceOut( "adjustSolution(...)" );
58}
59
61 const double * __restrict__ Q,
63 double t,
64 double * __restrict__ S
65) {
66 constexpr int nVars = 59;
67 for(int i=0; i<nVars; i++){
68 assertion3( std::isfinite(Q[i]), i, x, t );
69 }
70 memset(S, 0, nVars*sizeof(double));
71 //pdesource_(S,Q); // S(Q)
72 source(S,Q); // S(Q)
73 for(int i=0; i<nVars; i++){
74 nonCriticalAssertion3( std::isfinite(S[i]), i, x, t );
75 }
76}
77
79 const double * __restrict__ Q, // Q[64+0],
81 double t,
82 int normal
83) {
84#if defined(MGCCZ4EINSTEIN)
85 const double qmin = std::min({Q[0],Q[3],Q[5]});
86 const double alpha = std::max({1.0, std::exp(Q[16])}) * std::max({1.0, std::exp(Q[54])}) / std::sqrt(qmin);
87#else
88 const double alpha = 1.0;
89#endif
90
91 constexpr double sqrtwo = 1.4142135623730951;
92 // NOTE parameters are stored in Constants.h
93 const double tempA = alpha * std::max({sqrtwo, MGCCZ4e, MGCCZ4ds, MGCCZ4GLMc/alpha, MGCCZ4GLMd/alpha});
94 const double tempB = Q[17+normal];//DOT_PRODUCT(Q(18:20),nv(:))
96 return std::max({1.0, std::abs(-tempA-tempB), std::abs(tempA-tempB)});
97 /*
98 logTraceInWith3Arguments( "eigenvalues(...)", x, t, normal );
99 // helper data structure
100 constexpr int Unknowns = 59;
101 double lambda[Unknowns];
102 for (int i=0; i<Unknowns; i++) {
103 nonCriticalAssertion4( std::isfinite(Q[i]), i, x, t, normal );
104 lambda[i] = 1.0;
105 }
106
107 // routine requires explicit normal vector
108 double normalVector[3];
109 normalVector[0] = normal % 3 == 0 ? 1.0 : 0.0;
110 normalVector[1] = normal % 3 == 1 ? 1.0 : 0.0;
111 normalVector[2] = normal % 3 == 2 ? 1.0 : 0.0;
112
113 // actual method invocation
114 pdeeigenvalues_(lambda, Q, normalVector);
115
116 // we are only interested in the maximum eigenvalue
117 double result = 0.0;
118 for (int i=0; i<Unknowns; i++) {
119 result = std::max(result,std::abs(lambda[i]));
120 }
121
122 // @todo Raus
123 assertion3( std::isfinite(result), x, t, normal );
124 //nonCriticalAssertion3( std::isfinite(result), x, t, normal );
125
126 logTraceOut( "eigenvalues(...)" );
127 return result;*/
128}
129
130
132 const double * __restrict__ Q, // Q[64+0],
133 const double * __restrict__ deltaQ, // [64+0]
135 double t,
136 int normal,
137 double * __restrict__ BgradQ // BgradQ[64]
138) {
139 logTraceInWith3Arguments( "nonconservativeProduct(...)", x, t, normal );
140
141 assertion( normal>=0 );
142 assertion( normal<Dimensions );
143 constexpr int nVars = 59;
144 double gradQSerialised[nVars*3];
145 for (int i=0; i<nVars; i++) {
146 gradQSerialised[i+0*nVars] = 0.0;
147 gradQSerialised[i+1*nVars] = 0.0;
148 gradQSerialised[i+2*nVars] = 0.0;
149
150 gradQSerialised[i+normal*nVars] = deltaQ[i];
151 }
152 //pdencp_(BgradQ, Q, gradQSerialised);
153 ncp(BgradQ, Q, gradQSerialised, normal);
154
155 for (int i=0; i<nVars; i++) {
156 nonCriticalAssertion4( std::isfinite(BgradQ[i]), i, x, t, normal );
157 }
158
159 logTraceOut( "nonconservativeProduct(...)" );
160}
161
163 const double * __restrict__ Qinside,
164 double * __restrict__ Qoutside,
166 double t,
167 int normal
168) {
169 logTraceInWith3Arguments( "boundaryConditions(...)", x, t, normal );
170 for(int i=0; i<NumberOfUnknowns+NumberOfAuxiliaryVariables; i++) {
171 assertion4( Qinside[i]==Qinside[i], x, t, normal, i );
172 Qoutside[i]=Qinside[i];
173 }
174 logTraceOut( "boundaryConditions(...)" );
175}
176
#define assertion4(expr, param0, param1, param2, param3)
#define assertion3(expr, param0, param1, param2)
#define assertion(expr)
#define logError(methodName, logMacroMessageStream)
Wrapper macro around tarch::tarch::logging::Log to improve logging.
Definition Log.h:464
#define logTraceOut(methodName)
Definition Log.h:379
#define logTraceInWith3Arguments(methodName, argument0, argument1, argument2)
Definition Log.h:372
#define logTraceInWith2Arguments(methodName, argument0, argument1)
Definition Log.h:371
#define nonCriticalAssertion3(expr, param0, param1, param2)
#define nonCriticalAssertion4(expr, param0, param1, param2, param3)
virtual void algebraicSource(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, double *__restrict__ S) override
virtual void boundaryConditions(const double *__restrict__ Qinside, double *__restrict__ Qoutside, const tarch::la::Vector< Dimensions, double > &x, double t, int normal) override
static tarch::logging::Log _log
void adjustSolution(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t) override
void nonconservativeProduct(const double *__restrict__ Q, const double *__restrict__ deltaQ, const tarch::la::Vector< Dimensions, double > &x, double t, int normal, double *__restrict__ BgradQ) override
double maxEigenvalue(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t, int normal) override
Log Device.
Definition Log.h:516
void ncp(double *BgradQ, const double *const Q, const double *const gradQSerialised, const int normal, const int MGCCZ4LapseType, const double MGCCZ4ds, const double MGCCZ4c, const double MGCCZ4e, const double MGCCZ4f, const double MGCCZ4bs, const double MGCCZ4sk, const double MGCCZ4xi, const double MGCCZ4mu)
void source(double *S, const double *const Q, const int MGCCZ4LapseType, const double MGCCZ4ds, const double MGCCZ4c, const double MGCCZ4e, const double MGCCZ4f, const double MGCCZ4bs, const double MGCCZ4sk, const double MGCCZ4xi, const double MGCCZ4itau, const double MGCCZ4eta, const double MGCCZ4k1, const double MGCCZ4k2, const double MGCCZ4k3)
void enforceMGCCZ4constraints(double *luh)
void linearWave(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &X, double t)
void gaugeWave(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t)
bool equals(const Matrix< Rows, Cols, Scalar > &lhs, const Matrix< Rows, Cols, Scalar > &rhs, const Scalar &tolerance=NUMERICAL_ZERO_DIFFERENCE)
Compares to matrices on equality by means of a numerical accuracy.
Simple vector class.
Definition Vector.h:134