BoundaryOp refactor by johnomotani · Pull Request #1334 · boutproject/BOUT-dev

johnomotani · 2018-10-21T21:45:26Z

Simplify boundary_standard.cxx:

Most boundary conditions now use BoundaryOp::apply, which handles looping and staggering. The derived classes BoundaryNeumann, etc. mostly provide just methods to apply the boundary condition at a point, e.g.

BOUT-dev/src/mesh/boundary_standard.cxx

Lines 415 to 445 in cf7aa24

    
           BoundaryOp* BoundaryDirichlet_O3::clone(BoundaryRegion *region, const list<string> &args){ 
        
             verifyNumPoints(region, 2); 
        
             std::shared_ptr<FieldGenerator> newgen = nullptr; 
        
             if(!args.empty()) { 
        
               // First argument should be an expression 
        
               newgen = FieldFactory::get()->parse(args.front()); 
        
             } 
        
             return new BoundaryDirichlet_O3(region, newgen); 
        
           } 
        
           void BoundaryDirichlet_O3::applyAtPoint(Field2D &f, BoutReal val, int x, int bx, int y, int by, int z, Coordinates* UNUSED(metric)) { 
        
             f(x, y, z) = (8./3)*val - 2.*f(x - bx, y - by, z) + f(x - 2*bx, y - 2*by, z)/3.; 
        
           } 
        
           void BoundaryDirichlet_O3::applyAtPoint(Field3D &f, BoutReal val, int x, int bx, int y, int by, int z, Coordinates* UNUSED(metric)) { 
        
             f(x, y, z) = (8./3)*val - 2.*f(x - bx, y - by, z) + f(x - 2*bx, y - 2*by, z)/3.; 
        
           } 
        
           void BoundaryDirichlet_O3::applyAtPointStaggered(Field2D &f, BoutReal val, int x, int UNUSED(bx), int y, int UNUSED(by), int z, Coordinates* UNUSED(metric)) { 
        
             f(x, y, z) = val; 
        
           } 
        
           void BoundaryDirichlet_O3::applyAtPointStaggered(Field3D &f, BoutReal val, int x, int UNUSED(bx), int y, int UNUSED(by), int z, Coordinates* UNUSED(metric)) { 
        
             f(x, y, z) = val; 
        
           } 
        
           void BoundaryDirichlet_O3::extrapFurther(Field2D &f, int x, int bx, int y, int by, int z) { 
        
             extrap3rd(f, x, bx, y, by, z); 
        
           } 
        
           void BoundaryDirichlet_O3::extrapFurther(Field3D &f, int x, int bx, int y, int by, int z) { 
        
             extrap3rd(f, x, bx, y, by, z); 
        
           }

BoundaryDirichlet still has its own apply() method because it handled 2nd (and 3rd, 4th, etc.) guard cells in a non-standard way, setting them equal to the boundary value expression (evaluated at the guard cell's location); I have left this as-is. BoundaryNeumann_NonOrthogonal, BoundaryZeroLaplace, BoundaryZeroLaplace2 and BoundaryConstLaplace also implement overriding apply methods.

Remove BoundaryOpBase class. This is replaced by using templates to implement several methods for both BoundaryOp and BoundaryOpPar with a minimum of code duplication. Not totally sure this form is better/clearer. It does avoid several static_casts. Could undo this, but would take a little digging because I didn't do a good job of implementing this refactor in separate commits.

Fixes a few bugs: setting boundary conditions for fields at CELL_ZLOW (previously this case fell though and the guard cells were not changed); use the fieldmesh of the field, not the global mesh, in all boundary conditions; correct one of the coefficients in BoundaryDirichlet_4thOrder (now renamed to BoundaryDirichlet_O5; implement BoundaryNeumann_O4 for staggered grids.

Add GlobalZ method to Mesh for consistency with GlobalX, GlobalY.

Add a couple of higher order extrapolating boundary conditions, free_o4 and free_o5.

Change BoundaryRobin to take account of grid spacing, so now the derivative coefficient passed in (b) does not have to include dx or dy. Effectively derivative in boundary condition is now DDX or DDY instead of indexDDX or indexDDY.

Test the order of accuracy of Dirichlet, Neumann and free boundary conditions in MMS/derivatives3.

Add new boundary condition dirichlet_smooth, which suppresses grid-scale oscillations at the boundary enough to be used in wave-1d/wave-1d-y MMS tests. It has the same accuracy as the standard dirichlet condition, but uses two grid points to set the boundary condition. Setting the guard cell g by extrapolating from the boundary value, v, and first grid cell, f1 would give 'g1 = 2v-f1'; setting from the second grid cell would give 'g2 = 4/3v-1/3f2'; this bc takes the average of the two, setting 'g = 4/3v-1/2f1-1/6f2'.

Test both dirichlet_smooth and neumann boundary conditions for StaggerGrids=true or false in MMS/wave-1d and MMS/wave-1d-y.

Test both dirichlet and neumann boundary conditions on both boundaries in MMS/diffusion.

Also add override keywords for GlobalX, GlobalY and GlobalZ in BoutMesh declaration.

Add 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().

In Field3D::operator=(Field2D) move the checks on fieldmesh consistency to after the call to allocate(). If fields were declared in global scope and therefore initialized with a null Mesh*, fieldmesh may not be set to a sensible value until allocate() has been called.

If fieldmesh needs to be set, so will fielddatamesh. Therefore set it to global mesh in allocate() if fieldmesh is being set.

Put most of the functionality into the base BoundaryOp class, removing copy-paste code from boundary_standard.cxx. Use some templates to further shorten code.

called "free_o4" and "free_o5".

Makes it consistent with standard naming scheme. Can be called with either option "dirichlet_o5" or (now deprecated) "dirichlet_4thorder".

Previously BoundaryRobin::apply did not take account of the grid spacing dx/dy for the derivative term, so dx/dy would have to be included in its input.

Apply x- and y-boundary conditions in MMS/derivatives3 to test their accuracy and correctness. Tests Dirichlet, Neumann and Free boundary conditions.

Dirichlet boundary condition, tries to smooth out grid-scale oscillations at the boundary. Dirichlet bc using val and first grid point would be fb = 2*val - f0 using val and second grid point would be fb = 4/3*val - 1/3*f1 Here we apply the bc using the average of the two, to try and suppress grid-scale oscillations at the boundary.

Adds extra cases to runtest. Check both inner boundaries Neumann with outer boundaries Dirichlet and the opposite way round. Check true and false for mesh:staggergrids.

MMS/wave-1d-y

johnomotani · 2018-10-21T21:48:53Z

This should replace #1319.

dschwoerer · 2018-10-22T00:54:19Z

Change BoundaryRobin to take account of grid spacing, so now the derivative coefficient passed in (b) does not have to include dx or dy. Effectively derivative in boundary condition is now DDX or DDY instead of indexDDX or indexDDY.

Does it take only a uniform dx/dy into account, or also a non-uniform dx/dy?

I am asking, because to my understanding the others (I had a look at free, dirichlet, neumann) do not as far as I can tell, and that is why I added #1179

d7919 · 2018-10-22T07:33:06Z

include/boundary_standard.hxx

-  void apply(Field3D &f) override;
-  void apply(Field3D &f,BoutReal t) override;
+  void apply(Field2D &f,BoutReal t = 0.) override {
+    applyTemplate<Field2D>(f, t);


You probably don't actually need to specify Field2D or Field3D here as the compiler should be able to infer the type, although being explicit certainly doesn't hurt.

include/boundary_op.hxx

johnomotani · 2018-10-22T08:54:58Z

Does it take only a uniform dx/dy into account, or also a non-uniform dx/dy?

@dschwoerer no, the Robin bc doesn't account for variation of dx/dy, just uses the value dx/dy at the boundary grid point. Sorry, I hadn't looked at #1179 before now.

d7919 · 2018-10-22T09:34:26Z

Sorry for the duplicate comments!

BoundaryOp::applyAtPoint() and BoundaryOp::applyAtPointStaggered() should never be called, but must be implemented in the base class in case a subclass overrides BoundaryOp::apply() and so does not need these methods. This commit provides a helpful error message if these methods are ever called, instead of the previous ASSERT1(false).

These were not needed because the compiler can infer the template arguments from the type (Field2D/Field3D) of the function argument.

Also change BoundaryRegionOp from using is_same to is_convertible, in case template needs to be used with subclasses of BoundaryRegion/BoundaryRegionPar.

johnomotani · 2018-11-12T13:02:51Z

@bendudson what is the reason for having the BoundaryWidth modifier? It's being trickier than I'd anticipated to update; understanding the purpose would help make sure I don't break anything...

BoundaryOp implementations (e.g. BoundaryNeumann, etc.) that use a default 'apply()' method provided by a base class now inherit from BoundaryOpWithApply, which is templated on the derived class type. This means that its 'apply()' method can use static methods of the derived class directly, so that they can be inlined. This should improve performance when optimization is turned on. CRTP pattern: https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern

Move 'if (fg)' condition outside loops in boundary ops Introduces some code duplication, but significantly decreases time to apply boundary ops. Only calculate xnorm/ynorm if FieldGenerator fg is set Add some inline keywords

Cannot implement 'width' properly with member of BoundaryOp. Un-deprecates BoundaryWidth modifier, since this is now implemented (and needed to implement the 'width' keyword). Add 'copy' method to BoundaryRegion. Allows creating a new BoundaryRegion with a different width. Default value for 'wid' argument is -1, in which case xstart/ystart is used and the returned object is an exact copy of the original. BoundaryWidth modifier now has a member 'unique_ptr<BoundaryRegion> bndry' which is copied from the input BoundaryRegion, but with the width changed.

bendudson · 2018-11-12T17:30:14Z

@johnomotani I'm a bit unsure what the purpose was of the width modifier, and since it doesn't seem to be used in any of the examples, it can probably be removed. At least I wouldn't worry too much about trying to preserve it. I think it was either for:

Operators like C4
Staggered cells, where in the simple implementation the boundary needs to apply to an extra point. This extra point is now handled properly (?) by the boundary condition operators.

Otherwise if setBoundary() is called twice, the first set of boundary conditions will be applied and then the values overwritten by the second set. This is inefficient, and might cause unexpected behaviour, e.g. if the second set of boundary conditions specified "none" for some boundaries where the first set did not.

In Mesh/BoutMesh keep boundary objects in a vector< std::unique_ptr<BoundaryRegion> > and vector< std::unique_ptr<BoundaryRegionPar> > instead of vector<BoundaryRegion*> and vector<BoundaryRegionPar*>.

This previously did nothing, although it is not used anywhere.

johnomotani · 2018-11-12T18:27:40Z

Thanks @bendudson 👍

The tricky part about width keyword or BoundaryWidth modifier is that the BoundaryRegion needs to have its width changed. width is used in e.g. BoundaryRegionXIn::first(), so it's not enough to just replace uses of bndry->width in boundary_standard.cxx.

The neatest solution I've thought of is to use the BoundaryWidth modifier and have it make and store a copy of the BoundaryRegion with a different width. The downside is that this required adding a copy(int wid) method to BoundaryRegion. Suggestions for improvements or better designs welcome... or maybe we should just get rid of BoundaryWidth...

A few other changes:

I noticed FieldData::setBoundary(name) just added boundary conditions to the bndry_op vector without checking it was empty first. Probably it only ever gets called once, but just in case 1f588a4 makes it delete the entries and clear the vector. If it was called more than once, multiple boundary conditions could be applied on the same boundary. The last one would overwrite the others, so should be correct but inefficient.
I store the BoundaryRegion copy in BoundaryWidth as a std::unique_ptr and I thought it would be nice to have unique_ptrs for the boundary and par_boundary vectors in Mesh/BoutMesh. I've implemented in 2e17480. I guess it breaks backward compatibility in principle, although I wouldn't expect end-users to need to call any of the changed functions, Mesh::getBoundaries() or Mesh::getBoundariesPar(). Should I keep or get rid of the commit?
FieldData::setBoundary(region, op) was not implemented: in next at the moment it silently does nothing (except possibly printing "Replacing "). I've implemented what I think it was intended to do in ee22ccf but it might be better to just remove it. Obviously nobody is relying on it. Should I delete instead?

johnomotani · 2018-11-12T18:39:15Z

There was a comment in the manual about BoundaryWidth not being thread-safe because it changed the width of the BoundaryRegion while applying the BoundaryOp and then changed it back. That is no longer true, but are the boundary conditions thread-safe in general? I don't have a good understanding of what thread-safety means, but it looks to me like they are not, because they use shared BoundaryRegion objects as iterators, and BoundaryRegion::first(), BoundaryRegion::next1d(), etc. change the state of the shared object so if different threads were trying to apply boundary conditions to different fields at the same time, they might interfere with each other. Is that right? I've put a comment about the boundary conditions not being thread-safe in the manual, but someone who understands the concept better could probably write much more clearly...
See https://github.com/boutproject/BOUT-dev/blob/ebda6f316a184e419876a573a4bdd58a44b425aa/manual/sphinx/user_docs/boundary_options.rst#limitations

dschwoerer · 2018-11-14T16:43:37Z

If you are not sure, could you not mention that at all? Having wrong documentation is IMO way more confusing then no documentation. I don't think BOUT++ generally makes any sort of guarantees in that regard, thus saying one thing isn't thread safe implies the rest is - which might be the case, but I don't know.

Removes need for boundaryIsCopy flag, as well as some loops that delete members of the vectors.

dschwoerer · 2024-01-17T16:24:11Z

@johnomotani do you still intend to follow up on this?

johnomotani · 2024-01-17T17:15:10Z

@johnomotani do you still intend to follow up on this?

No. I don't remember now what exactly the purpose was, so don't know if in principle it would be worth updating and merging. In any case I don't have any time to work on it now, so I'm not going to do anything. Anyone is welcome to take it over if they think it's useful, although at this point it might be easier to start over from scratch!

johnomotani and others added 20 commits October 21, 2018 16:48

Add Mesh::GlobalZ methods for consistency with GlobalX, GlobalY

9124368

Also add override keywords for GlobalX, GlobalY and GlobalZ in BoutMesh declaration.

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

2f80b55

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

ed76776

Also move FieldData ctor into header

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

95e2ff2

Copy or check fielddatamesh in Field3D/Field2D operator=

c6adde3

Check consistency of fieldmesh and fielddatamesh

e851da6

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().

Set fielddatamesh in Field2D/Field3D::allocate()

8a07c7e

If fieldmesh needs to be set, so will fielddatamesh. Therefore set it to global mesh in allocate() if fieldmesh is being set.

Refactor boundary conditions code

af75612

Put most of the functionality into the base BoundaryOp class, removing copy-paste code from boundary_standard.cxx. Use some templates to further shorten code.

Add 4th and 5th order extrapolation boundary conditions

5fe8d3a

called "free_o4" and "free_o5".

Utility funcs for 2nd-5th order extrapolation in boundary_standard.cxx

20b6973

Rename BoundaryDirichlet_4thOrder to BoundaryDirichlet_O5

d40cd41

Makes it consistent with standard naming scheme. Can be called with either option "dirichlet_o5" or (now deprecated) "dirichlet_4thorder".

Make 'Robin' boundary condition inputs independent of grid spacing

a0c85de

Previously BoundaryRobin::apply did not take account of the grid spacing dx/dy for the derivative term, so dx/dy would have to be included in its input.

Add optional region argument to boutcore.Field*D.applyBoundary()

5e23d20

Test boundary conditions in MMS/derivatives3

e65ea0d

Apply x- and y-boundary conditions in MMS/derivatives3 to test their accuracy and correctness. Tests Dirichlet, Neumann and Free boundary conditions.

Swap boundary conditions and test staggergrids=false in wave-1d

3cdb964

Adds extra cases to runtest. Check both inner boundaries Neumann with outer boundaries Dirichlet and the opposite way round. Check true and false for mesh:staggergrids.

Test StaggerGrids=false and both dirichlet_smooth/neumann in

6e4de0f

MMS/wave-1d-y

Use dirichlet and neumann in MMS/diffusion

cf7aa24

d7919 reviewed Oct 22, 2018

View reviewed changes

include/boundary_op.hxx Outdated Show resolved Hide resolved

johnomotani added 2 commits October 22, 2018 11:25

Remove redundant template arguments

9cde433

These were not needed because the compiler can infer the template arguments from the type (Field2D/Field3D) of the function argument.

boutproject deleted a comment from d7919 Oct 22, 2018

johnomotani added 2 commits October 31, 2018 16:47

Remove unused BoundaryOpRegion template

b62bc0f

Also change BoundaryRegionOp from using is_same to is_convertible, in case template needs to be used with subclasses of BoundaryRegion/BoundaryRegionPar.

Update test_boundary_factory.cxx, compatible with refactored BoundaryOp

773cf49

johnomotani mentioned this pull request Nov 7, 2018

BoundaryOp refactor performance regression fix? #1351

Merged

johnomotani added the work in progress Not ready for merging label Nov 12, 2018

johnomotani added 6 commits November 12, 2018 17:02

Example to time application of boundary conditions

b248225

Pass (BoutReal) grid spacing instead of Coordinates*

3c0c6a3

Speed up BoundaryOps

5fcc7df

Move 'if (fg)' condition outside loops in boundary ops Introduces some code duplication, but significantly decreases time to apply boundary ops. Only calculate xnorm/ynorm if FieldGenerator fg is set Add some inline keywords

clang-format, remove commented-out code

9e16181

johnomotani force-pushed the BoundaryOp-refactor branch 2 times, most recently from 1703fa2 to ee9b86b Compare November 12, 2018 17:27

johnomotani added 2 commits November 12, 2018 17:45

Store BoundaryRegion's and BoundarRegionPar's in unique_ptr's

2e17480

In Mesh/BoutMesh keep boundary objects in a vector< std::unique_ptr<BoundaryRegion> > and vector< std::unique_ptr<BoundaryRegionPar> > instead of vector<BoundaryRegion*> and vector<BoundaryRegionPar*>.

johnomotani force-pushed the BoundaryOp-refactor branch from ee9b86b to ce640f0 Compare November 12, 2018 17:52

Implement FieldData::setBoundary(region, op)

ee22ccf

This previously did nothing, although it is not used anywhere.

johnomotani force-pushed the BoundaryOp-refactor branch from ce640f0 to ebda6f3 Compare November 12, 2018 18:16

johnomotani force-pushed the BoundaryOp-refactor branch from 48af365 to e153669 Compare November 12, 2018 19:15

bendudson mentioned this pull request Nov 26, 2018

Boundary conditions fixes #1319

Closed

johnomotani added 3 commits December 5, 2018 21:06

Update manual entry on boundary conditions

9a4b72a

Fix unit tests: use unique_ptr for BoundaryRegion/BoundaryRegionPar's

94d4e0d

Use shared_ptr for FieldData::bndry_op/bndry_op_par

d37335b

Removes need for boundaryIsCopy flag, as well as some loops that delete members of the vectors.

johnomotani force-pushed the BoundaryOp-refactor branch from e153669 to d37335b Compare December 5, 2018 21:06

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

BoundaryOp refactor#1334

BoundaryOp refactor#1334
johnomotani wants to merge 46 commits intonextfrom
BoundaryOp-refactor

johnomotani commented Oct 21, 2018

Uh oh!

johnomotani commented Oct 21, 2018

Uh oh!

dschwoerer commented Oct 22, 2018

Uh oh!

d7919 Oct 22, 2018

Uh oh!

Uh oh!

johnomotani commented Oct 22, 2018

Uh oh!

d7919 commented Oct 22, 2018

Uh oh!

johnomotani commented Nov 12, 2018

Uh oh!

bendudson commented Nov 12, 2018

Uh oh!

johnomotani commented Nov 12, 2018

Uh oh!

johnomotani commented Nov 12, 2018

Uh oh!

dschwoerer commented Nov 14, 2018 via email

Uh oh!

dschwoerer commented Jan 17, 2024

Uh oh!

johnomotani commented Jan 17, 2024

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

	BoundaryOp* BoundaryDirichlet_O3::clone(BoundaryRegion *region, const list<string> &args){
	verifyNumPoints(region, 2);

	std::shared_ptr<FieldGenerator> newgen = nullptr;
	if(!args.empty()) {
	// First argument should be an expression
	newgen = FieldFactory::get()->parse(args.front());
	}
	return new BoundaryDirichlet_O3(region, newgen);
	}

	void BoundaryDirichlet_O3::applyAtPoint(Field2D &f, BoutReal val, int x, int bx, int y, int by, int z, Coordinates* UNUSED(metric)) {
	f(x, y, z) = (8./3)val - 2.f(x - bx, y - by, z) + f(x - 2bx, y - 2by, z)/3.;
	}
	void BoundaryDirichlet_O3::applyAtPoint(Field3D &f, BoutReal val, int x, int bx, int y, int by, int z, Coordinates* UNUSED(metric)) {
	f(x, y, z) = (8./3)val - 2.f(x - bx, y - by, z) + f(x - 2bx, y - 2by, z)/3.;
	}

	void BoundaryDirichlet_O3::applyAtPointStaggered(Field2D &f, BoutReal val, int x, int UNUSED(bx), int y, int UNUSED(by), int z, Coordinates* UNUSED(metric)) {
	f(x, y, z) = val;
	}
	void BoundaryDirichlet_O3::applyAtPointStaggered(Field3D &f, BoutReal val, int x, int UNUSED(bx), int y, int UNUSED(by), int z, Coordinates* UNUSED(metric)) {
	f(x, y, z) = val;
	}

	void BoundaryDirichlet_O3::extrapFurther(Field2D &f, int x, int bx, int y, int by, int z) {
	extrap3rd(f, x, bx, y, by, z);
	}
	void BoundaryDirichlet_O3::extrapFurther(Field3D &f, int x, int bx, int y, int by, int z) {
	extrap3rd(f, x, bx, y, by, z);
	}

Conversation

johnomotani commented Oct 21, 2018

Uh oh!

johnomotani commented Oct 21, 2018

Uh oh!

dschwoerer commented Oct 22, 2018

Uh oh!

d7919 Oct 22, 2018

Choose a reason for hiding this comment

Uh oh!

Uh oh!

johnomotani commented Oct 22, 2018

Uh oh!

d7919 commented Oct 22, 2018

Uh oh!

johnomotani commented Nov 12, 2018

Uh oh!

bendudson commented Nov 12, 2018

Uh oh!

johnomotani commented Nov 12, 2018

Uh oh!

johnomotani commented Nov 12, 2018

Uh oh!

dschwoerer commented Nov 14, 2018 via email

Uh oh!

dschwoerer commented Jan 17, 2024

Uh oh!

johnomotani commented Jan 17, 2024

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants