Boundary conditions fixes by johnomotani · Pull Request #1319 · boutproject/BOUT-dev

johnomotani · 2018-10-12T17:20:00Z

Fixes 2 bugs:

calling applyBoundary() on a field with location CELL_ZLOW did nothing, because the conditionals missed CELL_ZLOW.
Fields created on meshes that aren't the global mesh could be wrong because the methods in boundary_standard.cxx all used the global mesh.

When iterating over flags from h5cc, the HDF5_{CPPFLAGS,LDFLAGS,LIBS} variables were assembled in reverse order. This causes problems when linking against static libraries. Fixes #961

Fix order of HDF5 flags

Now consistently treats CELL_DEFAULT

If the simulation is restarted (`restart`), and the old dump files have already been squashed, dump on restart is enabled by default. This patch checks whether the two times are the same, and if so, does not duplicate that time slice.

For some reason the time index is called `x2` - and thus this condition is not matched. This fixes this, as the newly created output should be more reliable in this regard.

`glob` was previously not correctly imported, if append mode was disabled, but delete was enabled.

coordinates. Replaces throw with ASSERT

* master: AX_LIB_HDF5: fix order of flags

Advective and flux derivatives, which take v and f as inputs, had cases using yup/ydown fields for one of v and f, but not the other. These do not make sense as multiplication of a field in field-aligned coordinates with another in non-field-aligned coordinates is incorrect. Delete last of these cases and fall back to converting both v and f to field-aligned if either does not have yup/ydown fields. Also add some comments where mixed cases were removed before.

In several y-derivatives where we have to convert the input to field-aligned coordinates, the result was not transformed back to the original coordinates. Also remove uses of global 'mesh' in Vpar_Grad_par_LCtoC()

Was previously a typo that resulted in using non-field-aligned f in VDDY, which is incorrect.

The location of the result will in future need to be set before shifting it from field-aligned coordinates, if this is necessary. So it is safer to set the location to outloc as soon as 'Field3D result' is declared. Also in a couple of places just 'return apply*diff(...)' instead of creating an unneeded intermediate variable 'result'.

* next: Move handling of outloc==CELL_DEFAULT above location checks Use ASSERT1 instead of ASSERT0 for location checking in derivatives Remove unnecessary conditionals in applyXdiff() Tidy up division in D2DYDZ Simplify location checking in index_derivs.cxx Remove mixed yup/ydown-fieldaligned cases in advective/flux derivatives Move checks for no guard cells to ASSERT1 at top of functions Add ASSERT1s to check that location of Coordinates is outloc Add missing outloc argument in indexFDDY(Field2D,Field2D) Check 'vloc != diffloc' in indexFDDX Use RGN_NOY for y-derivative in D2DXDY(Field2D) Fix location of Bxy in Div_par_flux Fix location checking in VDD* and FDD* Fix staggering in SPLIT branch of FDD* Fix staggering in Div_par_K_Grad_par Don't use global mesh in bracket_location Correct dy in D2DYDZ Fix sign in Div_par_flux Fix Bxy location in Div_par

Ensure location of result set in vecops routines

d7919 · 2018-10-15T20:53:37Z

 /////////////////// ASSIGNMENT ////////////////////

 Vector2D & Vector2D::operator=(const Vector2D &rhs) {
+  fielddatamesh = rhs.fielddatamesh;


I don't recall seeing similar statements for operator= in the fields, should there be?

I think there should be... I will add.

@ZedThree could you push the changes you just made please?

While looking at Field3D::operator=(const Field2D &rhs), I have a question; should we set fieldmesh = rhs.getMesh();, or ASSERT1(fieldmesh == rhs.getMesh()); in this operator?

Since Field2D has nz=1, we can't correctly set nz for the Field3D in Field3D::operator=(const Field2D &rhs) (as is done in Field3D::operator=(const Field3D &rhs)), so I think we have to assert rather than copying the Mesh* pointer.

I think to be consistent with other places this should probably be an assert. If it is an assert we should ensure we assert for Field3D::operator=(const Field3D &rhs) etc. as well to be consistent.

More location/method fixes for derivatives

…ft-global-option Check TwistShift==true in ShiftedMetric, add ShiftWithoutTwist option to bypass

ZedThree · 2018-10-16T11:12:50Z

I had forgotten how horrifying the boundary conditions code is. Slightly regretting saying anything now...

@johnomotani I believe you've done the following in each BoundaryCondition::apply:

Turned this conditional into an ASSERT precondition:

  CELL_LOC loc = f.getLocation();	
  if(mesh->StaggerGrids && loc != CELL_CENTRE) {

Unindented the rest of the function

There are now three top-level branches in each ::apply:

CELL_XLOW
CELL_YLOW
everything else, i.e. CELL_CENTRE and CELL_ZLOW

Is that correct?

johnomotani · 2018-10-16T11:15:06Z

@ZedThree yes, exactly.

d7919 · 2018-10-16T11:36:01Z


  // Copy the data and data sizes
  fieldmesh = rhs.fieldmesh;
+  fielddatamesh = rhs.fielddatamesh;


Do we need to check that fieldmesh/fielddatamesh are identical? If so why do we need two variables?

They definitely should be the same. We should ASSERT1 in Field2D and Field3D.

There are two because there's one in each of the two base classes. I think this points towards our class hierarchy being a little weird

I was wondering if this might be an opportunity to address the class structure? We're currently introducing more opportunities to form an inconsistent state (albeit with appropriate checks to try to guard against this). All vectors and fields apart from fieldperp have fielddata and all fields have Field. I'd guess FieldPerp should also have fielddata (and if not why not?) and if so then it looks like FieldData is the generic stuff to all field like objects (which would presumably include getting the single relevant mesh) whilst Field provides the things fields have that vectors don't (coordinates, anything else?).

Or could Field inherit from FieldData?

d7919 · 2018-10-16T11:51:32Z

 }

 Vector3D & Vector3D::operator=(const Vector2D &rhs) {
+  fielddatamesh = rhs.x.getMesh();


Should this be getDataMesh() for consistency? (I know they should return the same thing)

d7919 · 2018-10-16T11:54:15Z

@@ -23,88 +23,92 @@
    lead to an out of bounds access error later but we add it here to provide a
    more explanatory message.
 */
-void verifyNumPoints(BoundaryRegion *region, int ptsRequired) {
-  TRACE("Verifying number of points available for BC");
+namespace {


Are all the changes in this file separate from the FieldData gaining a mesh pointer?

I guess the changes are separate in the sense that they could all be made without adding a mesh pointer to FieldData, but the mesh pointer in FieldData is required for the boundary conditions to work correctly when the fieldmesh is not the same as the global mesh.

in verifyNumPoints we need to get localmesh from the BoundaryRegion *region; BoundaryRegion*s are created in a method of FieldData; FieldData had to gain a mesh pointer so that BoundaryRegion*s aren't always created with the global mesh.

the BoundaryRegion* bndry of the BoundaryOps could previously be wrong if the global and local Mesh are different (e.g. wrong number of points, or wrong number of guard cells, etc.).

d7919 · 2018-10-16T11:58:52Z

+      for(; !bndry->isDone(); bndry->next1d()) {
+        BoutReal xnorm = 0.5*(   localmesh->GlobalX(bndry->x)
+            + localmesh->GlobalX(bndry->x - bndry->bx) );
+        BoutReal ynorm = localmesh->GlobalY(bndry->y);


Should we introduce Mesh::GlobalZ for symmetry?

d7919 · 2018-10-16T12:00:28Z

+      }
+    }
+  } else {
+    // CELL_CENTRE or CELL_ZLOW


Do we not care about the Z location in the CELL_ZLOW case for 3D fields?

That's a good point... If we're using a constant rather than an expression for the value of the boundary condition, then we don't care about Z location, because there are no z-boundaries. But generally, we need to change the argument to fg->generate(), so there should be a separate CELL_ZLOW case :disappointed. I'll go fix this...

Going through all the boundary operators again will be painful... There's a lot of cut-and-paste code in boundary_standard.cxx, we could probably cut it down massively by having a single set of BoundaryOp::apply() methods implemented for the BoundaryOp class, and just having the single point stencils, one staggered and one unstaggered, for each derived class.

That would change quite a lot of code (on the other hand, I've touched most of the lines in boundary_standard.cxx already, at least for whitespace) so maybe you'd rather put off until 4.3. Also, since probably nobody is using CELL_ZLOW it doesn't matter too much if it stays broken for a little while longer (at least it's less broken now, constant value boundary conditions work).

I don't like leaving things broken though, so I'll try to make a fixed version, but I'll put it in a new PR.

I think we're going to put this PR off till 4.3 -- there's a few questions about Field/FieldData and the classes that do and don't inherit from them that we really need to resolve, and that will involve some possibly major structural changes.

All the boundary stuff needs a thorough look through and rethink as well. Can we use Region, or a subclass instead of BoundaryRegion? Is there a way we can template the boundary conditions on simple kernels and stencils? @d7919 has been doing something like this for derivatives, so maybe we can extend that work to here as well.

* v4.2-rc: (58 commits) Invert conditional in Div_par and Div_par_flux Make test-yupdown work with new test for TwistShift=true Check TwistShift==true in ShiftedMetric Remove output location from Curl properly Fix whitespace in vecops header Deprecate Curl overloads that take a location Fix locations of Jacobian in Div Set yup/ydown fields for Div_par_flux() Use curly braces with if statements Fix V_dot_Grad documentation Use std::map::emplace to avoid naming complicated type with insert Add method argument to Grad2_par2 Remove const on arguments passed by value Remove bracket_location() function, replace with ASSERT1 Tidying up more location setting in derivatives Remove unnecessary interp_to in Grad2_par2 Make mesh::getRegion?? const Add const(_iterator) overload of Region::begin/end Provide `getRegion` method on Fields to return the requested Region from fieldmesh Fix some sphinx whitespace/code blocks ...

Rename Mesh::coordinates -> Mesh::getCoordinates; deprecate old name

Previously, fields at CELL_ZLOW had nothing done by any boundary condition operator. This commit passes them through to the unstaggered case, which is correct since there are no z-boundaries so no special handling is needed.

Requires adding a Mesh* pointer and getDataMesh() method to FieldData, so that FieldData::setBoundary() can use the local mesh. getDataMesh() method is needed in case the FieldData is constructed before the global mesh is created (i.e. when Field3D, etc. are declared in global scope). It must not be called getMesh() because that would clash with the method of Field.

Give default argument 'Mesh* m = nullptr' for FieldData::FieldData constructor so that it is backward compatible. Method FieldData::getDataMesh() does not need to be virtual, since it is unlikely to be overridden.

Also move FieldData ctor into header

In places where fieldmesh and fielddatamesh are set from separate pointers, check that their values are equal. Also, for consistency set fielddatamesh from *.fielddatamesh or *.getDataMesh(), rather than fieldmesh or getMesh().

bendudson · 2018-11-26T08:53:38Z

Replaced by #1334

ZedThree and others added 30 commits March 23, 2018 11:23

AX_LIB_HDF5: fix order of flags

2909334

When iterating over flags from h5cc, the HDF5_{CPPFLAGS,LDFLAGS,LIBS} variables were assembled in reverse order. This causes problems when linking against static libraries. Fixes #961

Merge pull request #962 from boutproject/bugfix-961

8389878

Fix order of HDF5 flags

Use input field/vector location for CELL_DEFAULT in vecops

bdd81a1

Now consistently treats CELL_DEFAULT

add simple test for squashoutput

a9fcbf9

Fix append mode

f7d5117

If the simulation is restarted (`restart`), and the old dump files have already been squashed, dump on restart is enabled by default. This patch checks whether the two times are the same, and if so, does not duplicate that time slice.

Fix for Field2D

de6ffad

For some reason the time index is called `x2` - and thus this condition is not matched. This fixes this, as the newly created output should be more reliable in this regard.

Fix: glob not correctly imported

25cfad4

`glob` was previously not correctly imported, if append mode was disabled, but delete was enabled.

add dump_on_restart=false sub test

5006006

Fix formatting

a555391

Ensure CELL_DEFAULT and CELL_VSHIFT are forbidden in calls to

291b9f6

coordinates. Replaces throw with ASSERT

Bump version number to 4.2.0

f02880b

Include caveats on OpenMP parallelisation in documentation

8d68443

Merge branch 'master' into v4.2-rc

a5ce343

* master: AX_LIB_HDF5: fix order of flags

Add missing fromFieldAligned calls

90648ba

In several y-derivatives where we have to convert the input to field-aligned coordinates, the result was not transformed back to the original coordinates. Also remove uses of global 'mesh' in Vpar_Grad_par_LCtoC()

Fix f->f_fa

54494d4

Was previously a typo that resulted in using non-field-aligned f in VDDY, which is incorrect.

Add ASSERTs for checking requested location is unambigous/consistent.

791ebc7

Ensure location of result set in vecops routines

Ensure Vector copy and assignment copies the location

5a079d4

Expand unit tests to check that vector copy/assignment copies location

2f39306

Handle CELL_DEFAULT in vector setLocation

b3b756e

Add unit tests for vector setLocation

04f0eca

Fix for vector non-staggered setLocation unit tests

96dc8ab

Add staggered flag check in vector setLocation

ba80a08

Remove duplicated ASSERT

546ca78

Remove setLocation calls that are not required

8c571e9

Improvements after feedback

fce307d

Adding support for CELL_VSHIFT to Grad(Field2D)

1c4f9b9

Deprecating three arg version of Grad(Field3D)

83d6b82

d7919 reviewed Oct 15, 2018

View reviewed changes

Comment thread src/mesh/boundary_standard.cxx

ZedThree added 2 commits October 16, 2018 11:00

Merge pull request #1318 from boutproject/derivatives-more-tidyup

f586fbd

More location/method fixes for derivatives

Merge pull request #1323 from boutproject/shiftedmetric-test-twistshi…

c353ef9

…ft-global-option Check TwistShift==true in ShiftedMetric, add ShiftWithoutTwist option to bypass

d7919 reviewed Oct 16, 2018

View reviewed changes

ZedThree modified the milestones: BOUT-4.2, BOUT-4.3 Oct 16, 2018

ZedThree and others added 10 commits October 16, 2018 14:22

Merge pull request #1310 from boutproject/rename-mesh-coordinates

7516b23

Rename Mesh::coordinates -> Mesh::getCoordinates; deprecate old name

Handle CELL_ZLOW case in boundary_standard.cxx

3fe02cb

Previously, fields at CELL_ZLOW had nothing done by any boundary condition operator. This commit passes them through to the unstaggered case, which is correct since there are no z-boundaries so no special handling is needed.

Tidy up location checking in boundary_standard.cxx

bfcf30f

Make FieldData() backward compatible, no virtual on getDataMesh()

ea613e1

Give default argument 'Mesh* m = nullptr' for FieldData::FieldData constructor so that it is backward compatible. Method FieldData::getDataMesh() does not need to be virtual, since it is unlikely to be overridden.

Move FieldData::fielddatamesh init into init list

38f4161

Also move FieldData ctor into header

Check Mesh*'s are the same in Field3D operator=(Field2D)

2de6e1d

Copy or check fielddatamesh in Field3D/Field2D operator=

2f36229

Check consistency of fieldmesh and fielddatamesh

df110a0

In places where fieldmesh and fielddatamesh are set from separate pointers, check that their values are equal. Also, for consistency set fielddatamesh from *.fielddatamesh or *.getDataMesh(), rather than fieldmesh or getMesh().

johnomotani force-pushed the boundary-conditions-fixes branch from 67773e9 to df110a0 Compare October 16, 2018 15:46

johnomotani changed the base branch from v4.2-rc to next October 16, 2018 15:46

johnomotani mentioned this pull request Oct 21, 2018

BoundaryOp refactor #1334

Open

bendudson closed this Nov 26, 2018

d7919 mentioned this pull request Dec 19, 2018

Use namespace for global variables 'mesh' and 'dump' #1470

Merged

ZedThree deleted the boundary-conditions-fixes branch October 13, 2023 09:32

Uh oh!

Conversation

johnomotani commented Oct 12, 2018

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

ZedThree commented Oct 16, 2018

Uh oh!

johnomotani commented Oct 16, 2018

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

bendudson commented Nov 26, 2018

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants