Peano 4
Loading...
Searching...
No Matches
BoundaryConditions.cpp
Go to the documentation of this file.
2
4#include "peano4/utils/Loop.h"
5
6#include "tarch/logging/Log.h"
8
11
12
14 std::function< void(
15 const double* __restrict__ Qinside,
16 double * __restrict__ Qoutside,
19 double t,
20 double dt,
21 int normal
22 ) > boundaryCondition,
25 double t,
26 double dt,
27 int numberOfGridCellsPerPatchPerAxis,
28 int overlap,
29 int unknownsPlusAuxiliaryVariables,
30 int faceNumber,
31 double* __restrict__ Q
32) {
33 static tarch::logging::Log _log( "exahype2::fd" );
34 logTraceInWith4Arguments( "applyBoundaryConditions(...)", faceCentre, patchSize, numberOfGridCellsPerPatchPerAxis, faceNumber);
35
37 faceNumber,
38 numberOfGridCellsPerPatchPerAxis,
39 overlap,
40 unknownsPlusAuxiliaryVariables,
41 0
42 );
43
44 tarch::la::Vector<Dimensions,double> volumeH = exahype2::fd::getGridCellSize(patchSize, numberOfGridCellsPerPatchPerAxis);
45 tarch::la::Vector<Dimensions,double> faceOffset = faceCentre - 0.5 * patchSize;
46 faceOffset(faceNumber%Dimensions) += 0.5 * patchSize(faceNumber%Dimensions);
47
48 dfore(volume,numberOfGridCellsPerPatchPerAxis,faceNumber % Dimensions,0) {
49 tarch::la::Vector<Dimensions,int> insideGridCell = volume;
50 tarch::la::Vector<Dimensions,int> outsideGridCell = volume;
51 tarch::la::Vector<Dimensions,double> x = faceOffset + tarch::la::multiplyComponents( tarch::la::convertScalar<double>(volume)+tarch::la::Vector<Dimensions,double>(0.5), volumeH);
52
53 x(faceNumber%Dimensions) -= 0.5 * volumeH(faceNumber%Dimensions);
54
55 for (int layer=0; layer<overlap; layer++) {
56 if (faceNumber<Dimensions) {
57 insideGridCell(faceNumber % Dimensions) = overlap-layer;
58 outsideGridCell(faceNumber % Dimensions) = overlap-1-layer;
59 }
60 else {
61 insideGridCell(faceNumber % Dimensions) = layer+overlap-1;
62 outsideGridCell(faceNumber % Dimensions) = layer+overlap;
63 }
64
66 "applyBoundaryConditions(...)",
67 faceCentre << " x " << faceNumber << ": " <<
68 insideGridCell << "->" << outsideGridCell << " (" << enumerator(insideGridCell,0) << "->" << enumerator(outsideGridCell,0) << "): " <<
69 *(Q + enumerator(insideGridCell,0))
70 );
71
72 boundaryCondition(
73 Q + enumerator(insideGridCell,0),
74 Q + enumerator(outsideGridCell,0),
75 x, volumeH, t, dt, faceNumber
76 );
77 }
78 }
79
80 logTraceOut( "applyBoundaryConditions(...)" );
81}
82
83
85 std::function< double(
86 const double* __restrict__ Q,
89 double t,
90 double dt,
91 int normal
92 ) > maxEigenvalue,
93 std::function< void(
94 double* __restrict__ Q,
97 ) > farFieldSolution,
100 double t,
101 double dt,
102 int numberOfGridCellsPerPatchPerAxis,
103 int overlap,
104 int numberOfUnknowns,
105 int numberOfAuxiliaryVariables,
106 int faceNumber,
107 const tarch::la::Vector<Dimensions,double>& systemOrigin,
108 double* __restrict__ Qold,
109 double* __restrict__ Qnew
110) {
111 static tarch::logging::Log _log( "exahype2::fd" );
112 logTraceInWith4Arguments( "applySommerfeldConditions(...)", faceCentre, patchSize, numberOfGridCellsPerPatchPerAxis, faceNumber);
113
115 faceNumber,
116 numberOfGridCellsPerPatchPerAxis,
117 overlap,
118 numberOfUnknowns+numberOfAuxiliaryVariables,
119 0
120 );
121
122 double* QInf = tarch::allocateMemory<double>( numberOfUnknowns+numberOfAuxiliaryVariables, tarch::MemoryLocation::Heap );
123
124 tarch::la::Vector<Dimensions,double> volumeH = exahype2::fd::getGridCellSize(patchSize, numberOfGridCellsPerPatchPerAxis);
125 tarch::la::Vector<Dimensions,double> faceOffset = faceCentre - 0.5 * patchSize;
126 faceOffset(faceNumber%Dimensions) += 0.5 * patchSize(faceNumber%Dimensions);
127
128 dfore(volume,numberOfGridCellsPerPatchPerAxis,faceNumber % Dimensions,0) {
129 tarch::la::Vector<Dimensions,int> insideGridCell = volume;
130 tarch::la::Vector<Dimensions,int> outsideGridCell = volume;
131 tarch::la::Vector<Dimensions,double> x = faceOffset + tarch::la::multiplyComponents( tarch::la::convertScalar<double>(volume)+tarch::la::Vector<Dimensions,double>(0.5), volumeH);
132
134
135 for (int layer=0; layer<overlap; layer++) {
136 if (faceNumber<Dimensions) {
137 insideGridCell(faceNumber % Dimensions) = overlap-layer;
138 outsideGridCell(faceNumber % Dimensions) = overlap-1-layer;
139 xOnLayer(faceNumber%Dimensions) = x(faceNumber%Dimensions) - layer * volumeH(faceNumber%Dimensions); //for each layer the x is different
140 }
141 else {
142 insideGridCell(faceNumber % Dimensions) = layer+overlap-1;
143 outsideGridCell(faceNumber % Dimensions) = layer+overlap;
144 xOnLayer(faceNumber%Dimensions) = x(faceNumber%Dimensions) + layer * volumeH(faceNumber%Dimensions); //for each layer the x is different
145 }
146
147 if (tarch::la::equals(t,0.0)){ //for the initial timestep, we need to fill the outside volume for old data as well.
148 for(int i=0; i<numberOfUnknowns+numberOfAuxiliaryVariables; i++) {
149 (Qold + enumerator(outsideGridCell,0))[i]=(Qold+enumerator(insideGridCell,0))[i];
150 }
151 }
152
153 const double rOnLayer=tarch::la::norm2(xOnLayer-systemOrigin);
154 const double waveSpeed=maxEigenvalue( Qold + enumerator(outsideGridCell,0), faceCentre, volumeH, t, dt, faceNumber % Dimensions);
155
156 logDebug(
157 "applySommerfeldConditions(...)",
158 faceCentre << " x " << faceNumber << ": " <<
159 insideGridCell << "->" << outsideGridCell << " (" << enumerator(insideGridCell,0) << "->" << enumerator(outsideGridCell,0) << "): " <<
160 *(Qnew + enumerator(insideGridCell,0))
161 );
162
163 farFieldSolution(
164 QInf,
165 xOnLayer,
166 volumeH
167 );
168
169 //do something here, we have Qold, Qnew, waveSpeed, xOnLayer, rOnLayer, dt, volumeH, Qinf
170 //we try the first approach first to see if it is already sufficient.
171 double dtInverse=1/dt;
172 double spaceFactor=waveSpeed*rOnLayer/volumeH(faceNumber%Dimensions)/xOnLayer(faceNumber%Dimensions);
173 //logInfo("text","factor "<< spaceFactor <<" waveSpeed "<<waveSpeed<<" rOnLayer "<<rOnLayer<<" volumeH "<<volumeH(faceNumber%Dimensions)
174 // <<" xOnLayer "<<xOnLayer(faceNumber%Dimensions));
175
176 for(int i=0; i<numberOfUnknowns; i++) {
177 double numerator=-(waveSpeed/rOnLayer)*((Qold + enumerator(outsideGridCell,0))[i]+(Qold+enumerator(insideGridCell,0))[i]-2*QInf[i])
178 +(+dtInverse-spaceFactor)*(Qold + enumerator(outsideGridCell,0))[i]
179 +(-dtInverse+spaceFactor)*(Qnew + enumerator( insideGridCell,0))[i]
180 +(+dtInverse+spaceFactor)*(Qold + enumerator( insideGridCell,0))[i];
181 //logInfo("numerator","numerator "<<numerator );
182 (Qnew + enumerator(outsideGridCell,0))[i]=numerator/(dtInverse+spaceFactor);
183 }
184 for(int i=numberOfUnknowns; i<numberOfUnknowns+numberOfAuxiliaryVariables; i++) {
185 (Qnew + enumerator(outsideGridCell,0))[i]=(Qold+enumerator(insideGridCell,0))[i];
186 } //in principle we do not need to assign auxiliary variables, but I do this here for safety
187
188 }
189 }
190
192
193 logTraceOut( "applySommerfeldConditions(...)" );
194}
195
196
198 std::function< double(
199 const double* __restrict__ Q,
200 const tarch::la::Vector<Dimensions,double>& faceCentre,
202 double t,
203 double dt,
204 int normal
205 ) > maxEigenvalue,
206 const tarch::la::Vector<Dimensions,double>& faceCentre,
208 double t,
209 double dt,
210 int numberOfGridCellsPerPatchPerAxis,
211 int overlap,
212 int numberOfUnknowns,
213 int numberOfAuxiliaryVariables,
214 int faceNumber,
215 const tarch::la::Vector<Dimensions,double>& systemOrigin,
216 double* __restrict__ Qold,
217 double* __restrict__ Qnew
218) {
220 maxEigenvalue,
221 [&] (
222 double* __restrict__ Q,
223 const tarch::la::Vector<Dimensions,double>& faceCentre,
225 ) -> void {
226 for (int i=0; i<numberOfUnknowns+numberOfAuxiliaryVariables; i++) {
227 Q[i] = 0.0;
228 }
229 },
230 faceCentre,
231 patchSize,
232 t,
233 dt,
234 numberOfGridCellsPerPatchPerAxis,
235 overlap,
236 numberOfUnknowns,
237 numberOfAuxiliaryVariables,
238 faceNumber,
239 systemOrigin,
240 Qold,
241 Qnew
242 );
243}
244
245
#define logDebug(methodName, logMacroMessageStream)
Definition Log.h:50
#define logTraceOut(methodName)
Definition Log.h:379
#define logTraceInWith4Arguments(methodName, argument0, argument1, argument2, argument3)
Definition Log.h:373
#define dfore(counter, max, dim, value)
This is an exclusive d-dimensional for loop.
Definition Loop.h:937
Log Device.
Definition Log.h:516
void applyBoundaryConditions(std::function< void(const double *__restrict__ Qinside, double *__restrict__ Qoutside, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, double dt, int normal) > boundaryCondition, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &patchSize, double t, double dt, int numberOfGridCellsPerPatchPerAxis, int overlap, int unknownsPlusAuxiliaryVariables, int faceNumber, double *__restrict__ Q)
Apply (generic) boundary conditions.
void applySommerfeldConditions(std::function< double(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &gridCellH, double t, double dt, int normal) > maxEigenvalue, std::function< void(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &gridCellH) > farFieldSolution, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &patchSize, double t, double dt, int numberOfGridCellsPerPatchPerAxis, int overlap, int numberOfUnknowns, int numberOfAuxiliaryVariables, int faceNumber, const tarch::la::Vector< Dimensions, double > &systemOrigin, double *__restrict__ Qold, double *__restrict__ Qnew)
static tarch::la::Vector< 2, double > getGridCellSize(const tarch::la::Vector< 2, double > &h, int numberOfGridCellsPerPatchPerAxis)
Definition PatchUtils.h:22
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.
Matrix< Rows, Cols, Scalar > multiplyComponents(const Matrix< Rows, X, Scalar > &lhs, const Matrix< X, Cols, Scalar > &rhs)
Scalar norm2(const Vector< Size, Scalar > &vector)
Computes the 2-norm of the vector, i.e.
void freeMemory(void *data, MemoryLocation location)
@ Heap
Create data on the heap of the local device.
tarch::logging::Log _log("examples::unittests")
Simple vector class.
Definition Vector.h:134