Peano 4
Loading...
Searching...
No Matches
GPRDR.cpp
Go to the documentation of this file.
1#include "GPRDR.h"
2#include "Constants.h"
3#include "InitialData.h"
4#include "PDE.h"
5#include "ODE.h"
6
8#include "peano4/utils/Loop.h"
9
10
11
12
13#include <algorithm>
15
16
17tarch::logging::Log examples::exahype2::gprdr::GPRDR::_log( "examples::exahype2::gprdr::GPRDR" );
18
19
20
22 const int length=Scenario.length();
23 initparameters_(&length,&Scenario[0]);
24
25 /*
26 if (constants.isValueValidString("cgfile")) {
27 std::string cgfile = constants.getValueAsString("cgfile");
28 if (constants.isValueValidDouble("cx") && constants.isValueValidDouble("cy") && constants.isValueValidInt("isbinary") ) {
29 const double cx = constants.getValueAsDouble("cx");
30 const double cy = constants.getValueAsDouble("cy");
31 const int isbinary = constants.getValueAsInt("isbinary");
32 const int lengthcg=cgfile.length();
33
34 logInfo("init(...)","CG File:"<<cgfile);
35 logInfo("init(...)","CG center:("<<cx<<","<<cy<<")");
36 logInfo("init(...)","CG binary:"<<isbinary);
37 std::cout << cgfile << std::endl;
38 //loadcgfile_(&_domainOffset[0],&_domainSize[0],&lengthcg,&cgfile[0],&cx,&cy,&isbinary);
39 loadcgfile_(&_domainOffset[0],&_domainSize[0],&cx,&cy,&isbinary);
40 }else{
41 logInfo("init(...)","CG data CX is not correct.");
42 std::abort();
43 }
44}
45*/
46}
47
48
49
50// @todo I think this is underspecified, and we need dt here on the long run, too
52 double * __restrict__ Q,
55 double t
56) {
57 logTraceInWith3Arguments( "adjustSolution(...)", volumeX, volumeH, t );
58
59 double x[3];
60 #if Dimensions==2
61 x[0] = volumeX(0);
62 x[1] = volumeX(1);
63 x[2] = 0.0;
64 #else
65 x[0] = volumeX(0);
66 x[1] = volumeX(1);
67 x[2] = volumeX(2);
68 #endif
69
70 if (tarch::la::equals(t,0.0) ) {
71 //
72 // Initial conditions
73 //
74 int md = 0;
75 double cms = getMaxMeshSize();
76 const int order = _NumberOfFiniteVolumesPerAxisPerPatch;
77
78 std::fill_n(Q,NumberOfUnknowns,0.0);
79
80 initialdata_(x, &t, Q);
81 }
82 else {
83 //
84 // Point-source
85 //
86 double source[NumberOfUnknowns];
87 pdesource_(source, Q);
88 for (int i=0; i<NumberOfUnknowns; i++) {
89 Q[i] += getMinTimeStepSize() * source[i];
90 }
91
92 //
93 // Guess this is the actual rupture part
94 //
95 #if Dimensions==2
96 double slp = std::abs(Q[24]) ;
97 #endif
98
99 dynamicrupture_(x,&t,Q,&slp); //Duo April 28
100
101 double Q_cell_new[NumberOfUnknowns];
102 double dt = getMinTimeStepSize(); // not nice; see comment above
103 updatesolutionode_(Q_cell_new,Q,&dt);
104 std::copy_n(Q_cell_new,NumberOfUnknowns,Q);
105 }
106
107 //
108 // Sanity check
109 //
110 for(int i = 0; i< NumberOfUnknowns; i++){
111 // If check fails, we write one snapshot and then terminate. A real assertion
112 // would abbort immediately.
113 nonCriticalAssertion(std::isfinite(Q[i]));
114 }
115
116 logTraceOut( "adjustSolution(...)" );
117}
118
119
121 const double * __restrict__ QL,
122 const double * __restrict__ QR,
123 const tarch::la::Vector<Dimensions,double>& faceCentre,
124 double volumeH,
125 double t,
126 double dt,
127 int normal,
128 double * __restrict__ FL,
129 double * __restrict__ FR
130 ) {
131 hllemfluxfv_(FL, FR, QL, QR, &normal);
132}
133
134
135/*double examples::exahype2::gprdr::GPRDR::maxEigenvalue(
136 const double * __restrict__ Q,
137 const tarch::la::Vector<Dimensions,double>& faceCentre,
138 const tarch::la::Vector<Dimensions,double>& volumeH,
139 double t,
140 int normal
141) {
142 logTraceInWith4Arguments( "eigenvalues(...)", faceCentre, volumeH, t, normal );
143
144 double nv[3] = {0.0};
145 nv[normal] = 1.0;
146 double lambda[NumberOfUnknowns];
147 pdeeigenvalues_(lambda, Q, nv);
148
149 double result = 0.0;
150 for (int i=0; i<NumberOfUnknowns; i++) {
151 result = std::max( result, std::abs(lambda[i]) );
152 }
153
154 logTraceOutWith1Argument( "eigenvalues(...)", result );
155 return result;
156}
157
158
159
160void examples::exahype2::gprdr::GPRDR::flux(
161 const double * __restrict__ Q,
162 const tarch::la::Vector<Dimensions,double>& faceCentre,
163 const tarch::la::Vector<Dimensions,double>& volumeH,
164 double t,
165 int normal,
166 double * __restrict__ F
167) {
168 logTraceInWith4Arguments( "flux(...)", faceCentre, volumeH, t, normal );
169
170 double F_tmp[3][NumberOfUnknowns];
171 pdeflux_(F_tmp[0], F_tmp[1], F_tmp[2], Q);
172
173 for (int i=0; i<NumberOfUnknowns; i++) {
174 F[i] = F_tmp[normal][i];
175 }
176
177 logTraceOut( "flux(...)" );
178}
179
180
181
182
183void examples::exahype2::gprdr::GPRDR::nonconservativeProduct(
184 const double * __restrict__ Q,
185 const double gradQ[27+3][Dimensions],
186 const tarch::la::Vector<Dimensions,double>& faceCentre,
187 const tarch::la::Vector<Dimensions,double>& volumeH,
188 double t,
189 int normal,
190 double * __restrict__ BgradQ // BgradQ[27]
191) {
192 logTraceInWith4Arguments( "nonconservativeProduct(...)", faceCentre, volumeH, t, normal );
193
194 double gradQ_tmp[NumberOfUnknowns*3];
195
196 for (int i=0; i<NumberOfUnknowns; i++) {
197 gradQ_tmp[i+0*NumberOfUnknowns] = gradQ[i][0];
198 gradQ_tmp[i+1*NumberOfUnknowns] = gradQ[i][1];
199 #if Dimensions==2
200 gradQ_tmp[i+2*NumberOfUnknowns] = 0.0;
201 #else
202 gradQ_tmp[i+2*NumberOfUnknowns] = gradQ[i][2];
203 #endif
204 }
205
206 pdencp_(BgradQ, Q, gradQ_tmp);
207
208 logTraceOut( "nonconservativeProduct(...)" );
209}*/
210
#define logTraceOut(methodName)
Definition Log.h:379
#define logTraceInWith3Arguments(methodName, argument0, argument1, argument2)
Definition Log.h:372
#define nonCriticalAssertion(expr)
void updatesolutionode_(double *Qnew, const double *const Q0, const double *const dt)
void initparameters_(const int *length, const char *parsetup)
void solveRiemannProblem(const double *__restrict__ QL, const double *__restrict__ QR, const tarch::la::Vector< Dimensions, double > &faceCentre, double volumeH, double t, double dt, int normal, double *__restrict__ FL, double *__restrict__ FR) override
Definition GPRDR.cpp:120
static tarch::logging::Log _log
Definition GPRDR.h:26
void adjustSolution(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &volumeCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t) override
Definition GPRDR.cpp:51
Log Device.
Definition Log.h:516
void initialdata_(const double *x, const double *const t, double *Q)
void pdesource_(double *S, const double *const Q)
void dynamicrupture_(const double *x, const double *const t, double *Q, const double *slp)
void hllemfluxfv_(double *FL, double *FR, const double *const QL, const double *const QR, const int *normalNonZeroIndex)
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