Skip to content

Commit 8c5750d

Browse files
nicolasaunainicolasaunai
committed
fine_boundary_represents_var for copy **before** refinement
- PatchData NaN initialized on construction - fix tests failing as result of above - comment a field refinement test (useless, wrong refinement op for E,B) - debug plots for advance field overlap test - copy done before refinement (boolean false in variable) - overwrite_interior false also for refinement is default for FieldFillPattern - J manually init to zero in model init, fine init and regrid init (Jx unused in ampere but used in Ohm with its now NaN values) - Grid/NdarrayVector take default value overrides (for test) - UsableTensorField is default constructed with zero init.
1 parent 66f1230 commit 8c5750d

23 files changed

+422
-139
lines changed

‎src/amr/data/field/field_data.hpp‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -398,6 +398,11 @@ namespace amr
398398
SAMRAI::hier::Box const intersectionBox{box * transformedSource
399399
* destinationBox};
400400

401+
std::cout << "copy_ debug all boxes: "
402+
<< "sourceBox: " << sourceBox
403+
<< ", destinationBox: " << destinationBox
404+
<< ", transformedSource: " << transformedSource
405+
<< ", intersectionBox: " << intersectionBox << std::endl;
401406
if (!intersectionBox.empty())
402407
copy_<Operator>(intersectionBox, transformedSource, destinationBox,
403408
source, dst);

‎src/amr/data/field/field_geometry.hpp‎

Lines changed: 1 addition & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -28,8 +28,6 @@ namespace amr
2828
// generic BoxGeometry into the specific geometry but cannot cast into
2929
// the FieldGeometry below because it does not have the GridLayoutT and
3030
// PhysicalQuantity for template arguments.
31-
// this class is thus used instead and provide the method pureInteriorFieldBox()
32-
// used in FieldFillPattern::calculateOverlap()
3331
template<std::size_t dimension>
3432
class FieldGeometryBase : public SAMRAI::hier::BoxGeometry
3533
{
@@ -43,34 +41,17 @@ namespace amr
4341
, ghostFieldBox_{ghostFieldBox}
4442
, interiorFieldBox_{interiorFieldBox}
4543
, centerings_{centerings}
46-
, pureInteriorFieldBox_{pureInteriorBox_(interiorFieldBox, centerings)}
4744
{
4845
}
4946

50-
auto const& pureInteriorFieldBox() const { return pureInteriorFieldBox_; }
47+
auto const& interiorFieldBox() const { return interiorFieldBox_; }
5148

5249
SAMRAI::hier::Box const patchBox;
5350

5451
protected:
5552
SAMRAI::hier::Box const ghostFieldBox_;
5653
SAMRAI::hier::Box const interiorFieldBox_;
5754
std::array<core::QtyCentering, dimension> const centerings_;
58-
SAMRAI::hier::Box const pureInteriorFieldBox_;
59-
60-
private:
61-
static SAMRAI::hier::Box
62-
pureInteriorBox_(SAMRAI::hier::Box const& interiorFieldBox,
63-
std::array<core::QtyCentering, dimension> const& centerings)
64-
{
65-
auto noSharedNodeBox{interiorFieldBox};
66-
SAMRAI::hier::IntVector growth(SAMRAI::tbox::Dimension{dimension});
67-
for (auto dir = 0u; dir < dimension; ++dir)
68-
{
69-
growth[dir] = (centerings[dir] == core::QtyCentering::primal) ? -1 : 0;
70-
}
71-
noSharedNodeBox.grow(growth);
72-
return noSharedNodeBox;
73-
}
7455
};
7556

7657
template<typename GridLayoutT, typename PhysicalQuantity>

‎src/amr/data/field/field_variable.hpp‎

Lines changed: 10 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -29,13 +29,18 @@ namespace amr
2929
*
3030
* FieldVariable represent a data on a patch, it does not contain the data itself,
3131
* after creation, one need to register it with a context : see registerVariableAndContext.
32+
*
33+
*
34+
* Note that `fineBoundaryRepresentsVariable` is set to false so that
35+
* coarse-fine interfaces are handled such that copy happens **before**
36+
* refining. See https://github.com/LLNL/SAMRAI/issues/292
3237
*/
3338
FieldVariable(std::string const& name, PhysicalQuantity qty,
34-
bool fineBoundaryRepresentsVariable = true)
35-
: SAMRAI::hier::Variable(
36-
name,
37-
std::make_shared<FieldDataFactory<GridLayoutT, FieldImpl>>(
38-
fineBoundaryRepresentsVariable, computeDataLivesOnPatchBorder_(qty), name, qty))
39+
bool fineBoundaryRepresentsVariable = false)
40+
: SAMRAI::hier::Variable(name,
41+
std::make_shared<FieldDataFactory<GridLayoutT, FieldImpl>>(
42+
fineBoundaryRepresentsVariable,
43+
computeDataLivesOnPatchBorder_(qty), name, qty))
3944
, fineBoundaryRepresentsVariable_{fineBoundaryRepresentsVariable}
4045
, dataLivesOnPatchBorder_{computeDataLivesOnPatchBorder_(qty)}
4146
{

‎src/amr/data/field/field_variable_fill_pattern.hpp‎

Lines changed: 16 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -95,6 +95,7 @@ class FieldFillPattern : public SAMRAI::xfer::VariableFillPattern
9595
{
9696
NULL_USE(node_fill_boxes);
9797

98+
9899
/*
99100
* For this (default) case, the overlap is simply the intersection of
100101
* fill_boxes and data_box.
@@ -104,7 +105,21 @@ class FieldFillPattern : public SAMRAI::xfer::VariableFillPattern
104105

105106
SAMRAI::hier::BoxContainer overlap_boxes(fill_boxes);
106107
overlap_boxes.intersectBoxes(data_box);
107-
return pdf.getBoxGeometry(patch_box)->setUpOverlap(overlap_boxes, transformation);
108+
109+
auto geom = pdf.getBoxGeometry(patch_box);
110+
auto basic_overlap
111+
= pdf.getBoxGeometry(patch_box)->setUpOverlap(overlap_boxes, transformation);
112+
113+
if (overwrite_interior_)
114+
// if (true)
115+
return basic_overlap;
116+
117+
auto& overlap = dynamic_cast<FieldOverlap const&>(*basic_overlap);
118+
auto destinationBoxes = overlap.getDestinationBoxContainer();
119+
auto& casted = dynamic_cast<FieldGeometryBase<dimension> const&>(*geom);
120+
destinationBoxes.removeIntersections(casted.interiorFieldBox());
121+
122+
return std::make_shared<FieldOverlap>(destinationBoxes, overlap.getTransformation());
108123
}
109124

110125
bool overwrite_interior_;

‎src/amr/data/field/refine/electric_field_refiner.hpp‎

Lines changed: 73 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -94,7 +94,8 @@ class ElectricFieldRefiner
9494
//
9595
// therefore in all cases in 1D we just copy the coarse value
9696
//
97-
fineField(locFineIdx[dirX]) = coarseField(locCoarseIdx[dirX]);
97+
if (std::isnan(fineField(locFineIdx[dirX])))
98+
fineField(locFineIdx[dirX]) = coarseField(locCoarseIdx[dirX]);
9899
}
99100

100101
template<typename FieldT>
@@ -119,14 +120,16 @@ class ElectricFieldRefiner
119120
{
120121
// we're on a fine edge shared with coarse mesh
121122
// take the coarse face value
122-
fineField(ilfx, ilfy) = coarseField(ilcx, ilcy);
123+
if (std::isnan(fineField(ilfx, ilfy)))
124+
fineField(ilfx, ilfy) = coarseField(ilcx, ilcy);
123125
}
124126
else
125127
{
126128
// we're on a fine edge in between two coarse edges
127129
// we take the average
128-
fineField(ilfx, ilfy)
129-
= 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx, ilcy + 1));
130+
if (std::isnan(fineField(ilfx, ilfy)))
131+
fineField(ilfx, ilfy)
132+
= 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx, ilcy + 1));
130133
}
131134
}
132135
// Ey
@@ -140,14 +143,16 @@ class ElectricFieldRefiner
140143
// both fine Ey e.g. at j=100 and 101 will take j=50 on coarse
141144
// so no need to look at whether jfine is even or odd
142145
// just take the value at the local coarse index
143-
fineField(ilfx, ilfy) = coarseField(ilcx, ilcy);
146+
if (std::isnan(fineField(ilfx, ilfy)))
147+
fineField(ilfx, ilfy) = coarseField(ilcx, ilcy);
144148
}
145149
else
146150
{
147151
// we're on a fine edge in between two coarse ones
148152
// we take the average
149-
fineField(ilfx, ilfy)
150-
= 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx + 1, ilcy));
153+
if (std::isnan(fineField(ilfx, ilfy)))
154+
fineField(ilfx, ilfy)
155+
= 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx + 1, ilcy));
151156
}
152157
}
153158
// and this is now Ez
@@ -156,19 +161,29 @@ class ElectricFieldRefiner
156161
{
157162
if (onCoarseXFace_(fineIndex) and onCoarseYFace_(fineIndex))
158163
{
159-
fineField(ilfx, ilfy) = coarseField(ilcx, ilcy);
164+
if (std::isnan(fineField(ilfx, ilfy)))
165+
fineField(ilfx, ilfy) = coarseField(ilcx, ilcy);
160166
}
161167
else if (onCoarseXFace_(fineIndex))
162-
fineField(ilfx, ilfy)
163-
= 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx, ilcy + 1));
168+
{
169+
if (std::isnan(fineField(ilfx, ilfy)))
170+
fineField(ilfx, ilfy)
171+
= 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx, ilcy + 1));
172+
}
164173
else if (onCoarseYFace_(fineIndex))
165-
fineField(ilfx, ilfy)
166-
= 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx + 1, ilcy));
174+
{
175+
if (std::isnan(fineField(ilfx, ilfy)))
176+
fineField(ilfx, ilfy)
177+
= 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx + 1, ilcy));
178+
}
167179
else
168-
fineField(ilfx, ilfy)
169-
= 0.25
170-
* (coarseField(ilcx, ilcy) + coarseField(ilcx + 1, ilcy)
171-
+ coarseField(ilcx, ilcy + 1) + coarseField(ilcx + 1, ilcy + 1));
180+
{
181+
if (std::isnan(fineField(ilfx, ilfy)))
182+
fineField(ilfx, ilfy)
183+
= 0.25
184+
* (coarseField(ilcx, ilcy) + coarseField(ilcx + 1, ilcy)
185+
+ coarseField(ilcx, ilcy + 1) + coarseField(ilcx + 1, ilcy + 1));
186+
}
172187
}
173188
}
174189

@@ -197,33 +212,37 @@ class ElectricFieldRefiner
197212
// just copy the coarse value
198213
if (onCoarseYFace_(fineIndex) and onCoarseZFace_(fineIndex))
199214
{
200-
fineField(ilfx, ilfy, ilfz) = coarseField(ilcx, ilcy, ilcz);
215+
if (std::isnan(fineField(ilfx, ilfy, ilfz)))
216+
fineField(ilfx, ilfy, ilfz) = coarseField(ilcx, ilcy, ilcz);
201217
}
202218
// we share the Y face but not the Z face
203219
// we must be one of the 2 X fine edges on a Y face
204220
// thus we take the average of the two surrounding edges at Z and Z+DZ
205221
else if (onCoarseYFace_(fineIndex))
206222
{
207-
fineField(ilfx, ilfy, ilfz)
208-
= 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy, ilcz + 1));
223+
if (std::isnan(fineField(ilfx, ilfy, ilfz)))
224+
fineField(ilfx, ilfy, ilfz)
225+
= 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy, ilcz + 1));
209226
}
210227
// we share a Z face but not the Y face
211228
// we must be one of the 2 X fine edges on a Z face
212229
// we thus take the average of the two X edges at y and y+dy
213230
else if (onCoarseZFace_(fineIndex))
214231
{
215-
fineField(ilfx, ilfy, ilfz)
216-
= 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz));
232+
if (std::isnan(fineField(ilfx, ilfy, ilfz)))
233+
fineField(ilfx, ilfy, ilfz)
234+
= 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz));
217235
}
218236
else
219237
{
220238
// we don't share any face thus we're on one of the 2 middle X edges
221239
// we take the average of the 4 surrounding X averages
222-
fineField(ilfx, ilfy, ilfz)
223-
= 0.25 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz))
224-
+ 0.25
225-
* (coarseField(ilcx, ilcy, ilcz + 1)
226-
+ coarseField(ilcx, ilcy + 1, ilcz + 1));
240+
if (std::isnan(fineField(ilfx, ilfy, ilfz)))
241+
fineField(ilfx, ilfy, ilfz)
242+
= 0.25 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz))
243+
+ 0.25
244+
* (coarseField(ilcx, ilcy, ilcz + 1)
245+
+ coarseField(ilcx, ilcy + 1, ilcz + 1));
227246
}
228247
}
229248
// now this is Ey
@@ -235,7 +254,8 @@ class ElectricFieldRefiner
235254
if (onCoarseXFace_(fineIndex) and onCoarseZFace_(fineIndex))
236255
{
237256
// we thus just copy the coarse value
238-
fineField(ilfx, ilfy, ilfz) = coarseField(ilcx, ilcy, ilcz);
257+
if (std::isnan(fineField(ilfx, ilfy, ilfz)))
258+
fineField(ilfx, ilfy, ilfz) = coarseField(ilcx, ilcy, ilcz);
239259
}
240260
// now we only have same X face, but not (else) the Z face
241261
// so we're a new fine Y edge in between two coarse Y edges
@@ -247,27 +267,30 @@ class ElectricFieldRefiner
247267
// this means we are on a Y edge that lies in between 2 coarse edges
248268
// at z and z+dz
249269
// take the average of these 2 coarse value
250-
fineField(ilfx, ilfy, ilfz)
251-
= 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy, ilcz + 1));
270+
if (std::isnan(fineField(ilfx, ilfy, ilfz)))
271+
fineField(ilfx, ilfy, ilfz)
272+
= 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy, ilcz + 1));
252273
}
253274
// we're on a Z coarse face, but not on a X coarse face
254275
// we thus must be one of the 2 Y edges on a Z face
255276
// and thus we take the average of the 2 Y edges at X and X+dX
256277
else if (onCoarseZFace_(fineIndex))
257278
{
258-
fineField(ilfx, ilfy, ilfz)
259-
= 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz));
279+
if (std::isnan(fineField(ilfx, ilfy, ilfz)))
280+
fineField(ilfx, ilfy, ilfz)
281+
= 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz));
260282
}
261283
// now we're not on any of the coarse faces
262284
// so we must be one of the two Y edge in the middle of the cell
263285
// we thus average over the 4 Y edges of the coarse cell
264286
else
265287
{
266-
fineField(ilfx, ilfy, ilfz)
267-
= 0.25
268-
* (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz)
269-
+ coarseField(ilcx, ilcy, ilcz + 1)
270-
+ coarseField(ilcx + 1, ilcy, ilcz + 1));
288+
if (std::isnan(fineField(ilfx, ilfy, ilfz)))
289+
fineField(ilfx, ilfy, ilfz)
290+
= 0.25
291+
* (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz)
292+
+ coarseField(ilcx, ilcy, ilcz + 1)
293+
+ coarseField(ilcx + 1, ilcy, ilcz + 1));
271294
}
272295
}
273296
// now let's do Ez
@@ -279,34 +302,38 @@ class ElectricFieldRefiner
279302
// we thus copy the coarse value
280303
if (onCoarseXFace_(fineIndex) and onCoarseYFace_(fineIndex))
281304
{
282-
fineField(ilfx, ilfy, ilfz) = coarseField(ilcx, ilcy, ilcz);
305+
if (std::isnan(fineField(ilfx, ilfy, ilfz)))
306+
fineField(ilfx, ilfy, ilfz) = coarseField(ilcx, ilcy, ilcz);
283307
}
284308
// here we're on a coarse X face, but not a Y face
285309
// we must be 1 of the 2 Z edges on a X face
286310
// thus we average the 2 surrounding Z coarse edges at Y and Y+dY
287311
else if (onCoarseXFace_(fineIndex))
288312
{
289-
fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ])
290-
= 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz));
313+
if (std::isnan(fineField(ilfx, ilfy, ilfz)))
314+
fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ])
315+
= 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz));
291316
}
292317
// here we're on a coarse Y face, but not a X face
293318
// we must be 1 of the 2 Z edges on a Y face
294319
// thus we average the 2 surrounding Z coarse edges at X and X+dX
295320
else if (onCoarseYFace_(fineIndex))
296321
{
297-
fineField(ilfx, ilfy, ilfz)
298-
= 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz));
322+
if (std::isnan(fineField(ilfx, ilfy, ilfz)))
323+
fineField(ilfx, ilfy, ilfz)
324+
= 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz));
299325
}
300326
// we're not on any coarse face thus must be one of the 2 Z edges
301327
// in the middle of the coarse cell
302328
// we therefore take the average of the 4 surrounding Z edges
303329
else
304330
{
305-
fineField(ilfx, ilfy, ilfz)
306-
= 0.25
307-
* (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz)
308-
+ coarseField(ilcx, ilcy + 1, ilcz + 1)
309-
+ coarseField(ilcx + 1, ilcy + 1, ilcz));
331+
if (std::isnan(fineField(ilfx, ilfy, ilfz)))
332+
fineField(ilfx, ilfy, ilfz)
333+
= 0.25
334+
* (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz)
335+
+ coarseField(ilcx, ilcy + 1, ilcz + 1)
336+
+ coarseField(ilcx + 1, ilcy + 1, ilcz));
310337
}
311338
}
312339
}

‎src/amr/data/field/refine/field_refine_operator.hpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -108,6 +108,10 @@ class FieldRefineOperator : public SAMRAI::hier::RefineOperator
108108
auto intersectionBox = destFieldBox * box;
109109

110110

111+
std::cout << "debug refine operator: "
112+
<< "destinationFieldBox: " << destFieldBox
113+
<< ", sourceFieldBox: " << sourceFieldBox << ", box: " << box
114+
<< ", intersectionBox: " << intersectionBox << std::endl;
111115

112116

113117
if constexpr (dimension == 1)

‎src/amr/data/field/refine/field_refiner.hpp‎

Lines changed: 7 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -91,7 +91,8 @@ namespace amr
9191
{
9292
fieldValue += sourceField(xStartIndex + iShiftX) * leftRightWeights[iShiftX];
9393
}
94-
destinationField(fineIndex[dirX]) = fieldValue;
94+
if (std::isnan(destinationField(fineIndex[dirX])))
95+
destinationField(fineIndex[dirX]) = fieldValue;
9596
}
9697

9798

@@ -119,7 +120,8 @@ namespace amr
119120
fieldValue += Yinterp * xLeftRightWeights[iShiftX];
120121
}
121122

122-
destinationField(fineIndex[dirX], fineIndex[dirY]) = fieldValue;
123+
if (std::isnan(destinationField(fineIndex[dirX], fineIndex[dirY])))
124+
destinationField(fineIndex[dirX], fineIndex[dirY]) = fieldValue;
123125
}
124126

125127

@@ -157,7 +159,9 @@ namespace amr
157159
fieldValue += Yinterp * xLeftRightWeights[iShiftX];
158160
}
159161

160-
destinationField(fineIndex[dirX], fineIndex[dirY], fineIndex[dirZ]) = fieldValue;
162+
if (std::isnan(destinationField(fineIndex[dirX], fineIndex[dirY], fineIndex[dirZ])))
163+
destinationField(fineIndex[dirX], fineIndex[dirY], fineIndex[dirZ])
164+
= fieldValue;
161165
}
162166
}
163167

0 commit comments

Comments
 (0)