Peano 4
Loading...
Searching...
No Matches
PatchUtils.cpp
Go to the documentation of this file.
1// This file is part of the ExaHyPE2 project. For conditions of distribution and
2// use, please see the copyright notice at www.peano-framework.org
3#include "PatchUtils.h"
4
6#include "tarch/Assertions.h"
7#include "tarch/la/la.h"
8
9#include "peano4/utils/Loop.h"
10
12
13
14
16 int unknownsPlusAuxiliaryVariables,
17 int numberOfGridCellsPerPatchPerAxis,
18 int haloSize, // same as overlap
19 int normal,
20 bool isRightLayer,
21 const double* __restrict__ srcQ,
22 double* __restrict__ destQ
23) {
24 dfore(k,numberOfGridCellsPerPatchPerAxis,normal,0) {
25 for (int i= (isRightLayer ? haloSize : 0); i< (isRightLayer ? 2*haloSize : haloSize); i++) {
27 overlapCell(normal) = i;
28 const int index = toolbox::blockstructured::serialiseVoxelIndexInOverlap(overlapCell,numberOfGridCellsPerPatchPerAxis,haloSize,normal);
29 for (int j=0; j<unknownsPlusAuxiliaryVariables; j++) {
30 destQ[index*unknownsPlusAuxiliaryVariables+j] = srcQ[index*unknownsPlusAuxiliaryVariables+j];
31 }
32 }
33 }
34
35}
36
37
40 int numberOfVolumesPerAxisInPatch
41) {
42 return h(0)/numberOfVolumesPerAxisInPatch;
43}
44
45
48 int numberOfVolumesPerAxisInPatch
49) {
50 assertion2( numberOfVolumesPerAxisInPatch>=1, h, numberOfVolumesPerAxisInPatch );
51
52 return h(0)/numberOfVolumesPerAxisInPatch;
53}
54
55
57 const double* __restrict__ Q,
58 int unknowns
59) {
60 std::string result = "(" + std::to_string(Q[0]);
61 for (int i=1; i<unknowns; i++) result += "," + std::to_string(Q[i]);
62 result += ")";
63 return result;
64}
65
66
67namespace {
73 std::string prettyFormat( double value, bool isFirst ) {
74 static int previousZeroDataEntries = 0;
75
76 // first entry equals zero
77 if ( tarch::la::equals(value,0.0,1e-8) and isFirst ) {
78 previousZeroDataEntries = 1;
79 return "0";
80 }
81 // first entry does not equal zero
82 else if ( isFirst ) {
83 previousZeroDataEntries = 0;
84 return std::to_string(value);
85 }
86 #if PeanoDebug>0
87 else {
88 return "," + std::to_string(value);
89 }
90 #else
91 else if ( not tarch::la::equals(value,0.0,1e-8) ) {
92 previousZeroDataEntries = 0;
93 return "," + std::to_string(value);
94 }
95 else if ( tarch::la::equals(value,0.0,1e-8) and previousZeroDataEntries==1 ) {
96 previousZeroDataEntries += 1;
97 return ",...";
98 }
99 else {
100 assertion(false);
101 }
102 #endif
103 return "<undef>";
104 }
105}
106
107
109 const double* __restrict__ Q,
110 int unknowns,
111 int auxiliaryVariables,
112 int numberOfVolumesPerAxisInPatch,
113 int haloSize,
114 bool prettyPrint
115) {
116 const int PatchSize = numberOfVolumesPerAxisInPatch+2*haloSize;
117
118 std::ostringstream result;
119
120 if (prettyPrint) result << std::endl;
121
122 #if Dimensions==2
123 for (int y=0; y<PatchSize; y++) {
124 for (int x=0; x<PatchSize; x++) {
125 int index = (y*PatchSize+x) * (unknowns+auxiliaryVariables);
126
127 // It is a diagonal entry if this counter is bigger than 1. If it equals
128 // 0, it is interior. If it equals 1, then this is a face-connected halo
129 // entry.
130 bool isDiagonal = (x<haloSize and y<haloSize)
131 or (x>=numberOfVolumesPerAxisInPatch+haloSize and y<haloSize)
132 or (x<haloSize and y>=numberOfVolumesPerAxisInPatch+haloSize)
133 or (x>=numberOfVolumesPerAxisInPatch+haloSize and y>=numberOfVolumesPerAxisInPatch+haloSize);
134
135 result << "(";
136 for (int i=0; i<unknowns+auxiliaryVariables; i++) {
137 const int entry = index+i;
138
139 if (not isDiagonal)
140 result << prettyFormat( Q[entry], i==0);
141 else
142 result << "x";
143 }
144 result << ")";
145 }
146 result << std::endl;
147 }
148 #else
149 dfor (k,PatchSize) {
150 int index = peano4::utils::dLinearised(k,PatchSize) * (unknowns+auxiliaryVariables);
151
152 // It is a diagonal entry if this counter is bigger than 1. If it equals
153 // 0, it is interior. If it equals 1, then this is a face-connected halo
154 // entry.
155 bool isDiagonal = (tarch::la::count(k,0) + tarch::la::count(k,PatchSize-1))>1;
156
157 result << "(";
158 for (int i=0; i<unknowns+auxiliaryVariables; i++) {
159 const int entry = index+i;
160
161 if (i!=0) result << ",";
162
163 if (haloSize==0 or not isDiagonal)
164 result << Q[entry];
165 else
166 result << "x";
167 }
168 result << ")";
169 }
170 #endif
171
172 if (prettyPrint) result << std::endl;
173
174 return result.str();
175}
176
177
178
180 const double* __restrict__ Q,
181 int unknowns,
182 int auxiliaryVariables,
183 int numberOfVolumesPerAxisInPatch,
184 int haloSize,
185 int normal,
186 bool prettyPrint
187) {
188 assertion(normal>=0);
189 assertion(normal<Dimensions);
190
191 std::ostringstream result;
192 if (prettyPrint) {
193 result << std::endl;
194
195 #if Dimensions==2
196 if (normal==0) {
197 for (int y=0; y<numberOfVolumesPerAxisInPatch; y++) {
198 for (int x=0; x<2*haloSize; x++) {
199 int index = (y*numberOfVolumesPerAxisInPatch+x) * (unknowns+auxiliaryVariables);
200 result << "(";
201 for (int i=0; i<unknowns+auxiliaryVariables; i++) {
202 const int entry = index+i;
203 result << prettyFormat( Q[entry],i==0 );
204 }
205 result << ")";
206 if (x==haloSize-1)
207 result << " | ";
208 }
209 result << std::endl;
210 }
211 }
212 else {
213 for (int y=0; y<2*haloSize; y++) {
214 for (int x=0; x<numberOfVolumesPerAxisInPatch; x++) {
215 int index = (y*numberOfVolumesPerAxisInPatch+x) * (unknowns+auxiliaryVariables);
216 result << "(";
217 for (int i=0; i<unknowns+auxiliaryVariables; i++) {
218 const int entry = index+i;
219 result << prettyFormat( Q[entry],i==0 );
220 }
221 result << ")";
222 }
223 result << std::endl;
224
225 if (y==haloSize-1) {
226 for (int x=0; x<numberOfVolumesPerAxisInPatch; x++) {
227 result << " --- ";
228 }
229 }
230 }
231 }
232 #else
233 result << "<not implemented yet>";
234 #endif
235 result << std::endl;
236 }
237 else {
238 int entry = 0;
239 for (int overlap=0; overlap<2*haloSize; overlap++)
240 #if Dimensions==3
241 for (int i=0; i<numberOfVolumesPerAxisInPatch; i++)
242 #endif
243 for (int ii=0; ii<numberOfVolumesPerAxisInPatch; ii++)
244 {
245 if (ii==0) {
246 result << "(";
247 }
248 else {
249 result << ",(";
250 }
251 for (int iii=0; iii<unknowns+auxiliaryVariables; iii++) {
252 result << prettyFormat( Q[entry],iii==0 );
253 entry++;
254 }
255 result << ")";
256 }
257 }
258
259 return result.str();
260}
261
262
264 const double *__restrict__ Q, int unknowns,
265 int auxiliaryVariables,
266 int numberOfVolumesPerAxisInPatch, int haloSize,
267 const std::string &location,
268 bool triggerNonCriticalAssertion,
269 double* minValues,
270 double* maxValues
271) {
272 #if PeanoDebug>1 && !defined(GPUOffloadingOMP)
273 static tarch::logging::Log _log( "exahype2::fv" );
274 logTraceInWith5Arguments( "validatePatch(...)", unknowns, auxiliaryVariables, numberOfVolumesPerAxisInPatch, haloSize, location );
275 const int PatchSize = numberOfVolumesPerAxisInPatch+2*haloSize;
276 dfor (k,PatchSize) {
277 int index = peano4::utils::dLinearised(k,PatchSize) * (unknowns+auxiliaryVariables);
278
279 int outsidePatchAlongCoordinateAxis = 0;
280 for (int d=0; d<Dimensions; d++) {
281 if (k(d)<haloSize) outsidePatchAlongCoordinateAxis++;
282 if (k(d)>=haloSize+numberOfVolumesPerAxisInPatch) outsidePatchAlongCoordinateAxis++;
283 }
284 bool isDiagonal = outsidePatchAlongCoordinateAxis>1;
285
286 if (haloSize==0 or not isDiagonal) {
287 for (int i=0; i<unknowns+auxiliaryVariables; i++) {
288 const int entry = index+i;
289 bool dataValid = true;
290 if (minValues!=nullptr and minValues[i]>Q[entry]) {
291 dataValid = false;
292 }
293 if (maxValues!=nullptr and maxValues[i]<Q[entry]) {
294 dataValid = false;
295 }
296 if (triggerNonCriticalAssertion) {
298 Q[entry]==Q[entry] and std::isfinite(Q[entry]) and dataValid,
299 Q[entry], unknowns,
300 auxiliaryVariables,
301 isDiagonal, haloSize, entry,
302 numberOfVolumesPerAxisInPatch, haloSize, k, i, location );
303 assertion( Q[entry]==Q[entry] and std::isfinite(Q[entry]) and dataValid );
304 }
305 else {
307 Q[entry]==Q[entry] and std::isfinite(Q[entry]) and dataValid,
308 Q[entry], unknowns,
309 auxiliaryVariables,
310 isDiagonal, haloSize, i,
311 numberOfVolumesPerAxisInPatch, haloSize, k, i, location,
312 plotPatch(Q,unknowns,auxiliaryVariables,numberOfVolumesPerAxisInPatch,haloSize)
313 );
314 }
315 }
316 }
317 }
318 logTraceOut( "validatePatch(...)" );
319 #endif
320}
#define assertion2(expr, param0, param1)
#define assertion(expr)
#define assertion12(expr, param0, param1, param2, param3, param4, param5, param6, param7, param8, param9, param10, param11)
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
#define logTraceOut(methodName)
Definition Log.h:379
#define logTraceInWith5Arguments(methodName, argument0, argument1, argument2, argument3, argument4)
Definition Log.h:374
#define dfore(counter, max, dim, value)
This is an exclusive d-dimensional for loop.
Definition Loop.h:937
#define dfor(counter, max)
d-dimensional Loop
Definition Loop.h:308
#define nonCriticalAssertion11(expr, param0, param1, param2, param3, param4, param5, param6, param7, param8, param9, param10)
Log Device.
Definition Log.h:516
void validatePatch(const double *__restrict__ Q, int unknowns, int auxiliaryVariables, int numberOfVolumesPerAxisInPatch, int haloSize, const std::string &location="", bool triggerNonCriticalAssertion=true, double *minValues=nullptr, double *maxValues=nullptr)
Just runs over the patch and ensures that no entry is non or infinite.
std::string plotPatch(const double *__restrict__ Q, int unknowns, int auxiliaryVariables, int numberOfVolumesPerAxisInPatch, int haloSize, bool prettyPrint=false)
Plot patch.
double getVolumeLength(const tarch::la::Vector< 2, double > &h, int numberOfVolumesPerAxisInPatch)
With GCC 10, it was impossible to return/copy the vector class.
void copyHalfOfHalo(int unknownsPlusAuxiliaryVariables, int numberOfGridCellsPerPatchPerAxis, int haloSize, int normal, bool isRightLayer, const double *__restrict__ srcQ, double *__restrict__ destQ)
A face always holds a left and a right overlap.
std::string plotPatchOverlap(const double *__restrict__ Q, int unknowns, int auxiliaryVariables, int numberOfGridCellsPerPatchPerAxis, int haloSize, int normal, bool prettyPrint=false)
std::string plotVolume(const double *__restrict__ Q, int unknowns)
Helper routine that I need in the log statements.
CPUGPUMethod int dLinearised(const tarch::la::Vector< Dimensions, int > &counter, int max)
Map d-dimensional vector onto integer index.
Definition Loop.cpp:106
int count(const Vector< Size, Scalar > &vector, const Scalar &value)
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.
int serialiseVoxelIndexInOverlap(const tarch::la::Vector< Dimensions, int > &overlapCell, int numberOfDoFsPerAxisInPatch, int overlap, int normal)
The volumes or elements within an overlap are always enumerated lexicographically.
tarch::logging::Log _log("examples::unittests")
Simple vector class.
Definition Vector.h:134