Feature request: CellVariable method for explicitly defining no-flux boundary conditions

[mhvwerts]

FiPy's cell-centered finite-volume method has implicit no-flux boundary conditions (BCs). These BCs are used by default when boundary conditions on a CellVariable are not specified. They work extremely well in terms of conservation of total quantity represented by the cell variable (total mass, energy...). In many cases, these implicit no-flux BCs have superior numerical behaviour compared to any hand-crafted explicit no-flux conditions based on faceGrad.

Some users do (did) not know this (see Issue #676 and Issue #674), and even those that do may have a desire to specify the default no-flux BCs explicitly, with the aim of making their code more understandable.

It would be a nice feature to have a CellVariable method called no_flux_boundary() (for instance) that allows for explicitly stating the default no-flux boundaries.

It would allow to run the following code:

u = fp.CellVariable(name = 'concentration',
                    mesh = mesh,
                    value = 1.0)
                    
u.no_flux_boundary(where = mesh.facesLeft)
u.no_flux_boundary(where = mesh.facesRight)
u.no_flux_boundary(where = mesh.facesBottom)
u.constrain(1.0, where = mesh.facesTop)  # and a Dirichlet BC on top

In essence, with a freshly created CellVariable, the method would do nothing. However, for consistence, it should probably erase any prior BCs that were defined on the relevant boundary.

Alternatively, a simpler method: no_flux_boundaries() taking no further arguments would reset all boundaries to the default implicit no-flux conditions.

[wd15]

Thanks for this feedback. We used to have boundary condition objects rather than constraints and then we went with constraints and eliminated the boundary conditions. However, it has seemed to cause confusion here and there. I think this idea would work though. The no flux would apply to the equation that solves for "u" rather than to "u" itself. There is a difference between applying a constraint on "u" and applying a boundary condition to an equation that solves for "u" although we conflate them in FiPy.

[mhvwerts]

I was not aware of this subtlety, the distinction between a boundary condition and a constraint. In the recipe we give our students for using FiPy (mainly for diffusion/advection/reaction systems), we start with defining the mesh (domain), then creating the variable(s), then applying boundary conditions (via constraints, as recommended by the documentation) and finally specifying the PDE, and solving it (typically through time-stepping). I see that we should perhaps be more precise in our vocabulary.

Does application of constraints to a variable follow a numerical scheme that is fundamentally different from applying boundary conditions to an equation? (e.g. in the ways the corresponding matrices are constructed or matrix equations solved?) With my limited knowledge of the inner works of FiPy, I can imagine that constraints on a variable are somehow more convenient especially when working with several variables and coupled equations.

[guyer]

The numerical scheme should be identical, as far as the linear equations that are solved, but there are subtleties with boundary values in either case that aren't very satisfactory. We've found that significant numbers of FiPy users are confused by the behavior of both BoundaryCondition objects and constraints. We're still trying to figure out a better treatment.

I've put together a notebook that illustrates some simple examples of the differences.

[mhvwerts]

Many thanks! Your notebook is very insightful for understanding BCs and constraints and how these work internally. I'll follow your advice and will continue using constraints for implementing boundary conditions. I know of at least one student who will be very interested in this notebook and its implications.

[guyer]

(The link to the notebook should probably be https://gist.github.com/guyer/2bd27a25bdd40b4aa26e40222cf3d046 )

Fixed, thanks

[wd15]

we give our students for using FiPy (mainly for diffusion/advection/reaction systems), we start with defining the mesh (domain), then creating the variable(s), then applying boundary conditions (via constraints, as recommended by the documentation) and finally specifying the PDE, and solving it (typically through time-stepping).

Thanks for sharing that you are using FiPy in your teaching. Any feedback you have from your teaching experience could be very useful.

[guyer]

If you haven't seen it, you and your students may find this notebook instructive, too