Boundary conditions as a property of a CellVariable

[mhvwerts]

Closes #36.

Following the discussion in #36, the definition of the boundary conditions has been embedded in the cell variable object CellVariable, without breaking any existing code. While boundary conditions may still be defined separately before association with a cell variable, a new and more readable workflow now exists (see the example below, which was taken from the new README.md).

The basic code changes implementing this new feature are surprisingly modest, which is probably a good sign.

Perhaps more tests should be adapted to reflect the new and improved workflow, so that testing becomes more complete. Also, the example scripts should some day be converted to tests executed by pytest. Anyways, everything seems to still work fine.

import pyfvtool as pf

# Solving a 1D diffusion equation with a fixed concentration 
# at the left boundary and a closed boundary on the right side

# Calculation parameters
Nx = 20 # number of finite volume cells
Lx = 1.0 # [m] length of the domain 
c_left = 1.0 # left boundary concentration
c_init = 0.0 # initial concentration
D_val = 1e-5 # diffusion coefficient (gas phase)
t_simulation = 7200.0 # [s] simulation time
dt = 60.0 # [s] time step
Nskip = 10 # plot every Nskip-th profile

# Define mesh
mesh = pf.Grid1D(Nx, Lx)

# Create a cell variable with initial concentration
# By default, 'no flux' boundary conditions are applied
c = pf.CellVariable(mesh, c_init)

# Switch the left boundary to Dirichlet: fixed concentration
c.BCs.left.a[:] = 0.0
c.BCs.left.b[:] = 1.0
c.BCs.left.c[:] = c_left
c.apply_BCs()

# Assign diffusivity to cells
D_cell = pf.CellVariable(mesh, D_val)
D_face = pf.geometricMean(D_cell) # average value of diffusivity at the interfaces between cells

# Time loop
t = 0
nplot = 0
while t<t_simulation:
    # Compose discretized terms for matrix equation
    eqn = [ pf.transientTerm(c, dt, 1.0),
           -pf.diffusionTerm(D_face)]

    # Solve PDE
    pf.solvePDE(c, eqn)
    t+=dt

    if (nplot % Nskip == 0):
        pf.visualizeCells(c)
    nplot+=1

[mhvwerts]

Converting this to draft, since there are several places in the code where a new CellVariable is created to take the place of an input CellVariable, basically creating a modified copy. Here, the embedded BCs should also be copied between the input and the new output CellVariable.

This is in particular the case for CellVariable.__add__(), CellVariable.__sub__() and all the arithmetic dunder methods of CellVariable.

All current tests and examples work perfectly, though. This issue, however, might create unexpected behavior in more elaborate code.

[mhvwerts]

We should perform a search in the code for all calls to CellVariable and check. In particular, the cases where we have = CellVariable and return CellVariable and variants.

The copying is not necessarily needed in all cases (e.g. where a CellVariable is a non-solution variable).

[mhvwerts]

An initial survey (multi-file searches) suggests that the need for copying BCs only exists in cell.py (aforementioned arithmetic dunder methods and funceval).

[mhvwerts]

Perhaps CellVariable.value should become an internal (private) property CellVariable._value that the user is not supposed to touch, because it becomes confusing with the ghost cells included.

Most user cases are better handled via the CellVariable.innerCellValues property, followed by an explicit call to CellVariable.apply_BCs() where necessary.

[mhvwerts]

The call to CellVariable.apply_BCs() might perhaps be done automatically when setting innerCellValues. This would systematically keep the ghost cells 'in sync' with the latest inner cell values and the current boundary conditions.

[mhvwerts]

For now, let's already start with avoiding the use of CellVariable.value in user code (i.e. the supplied example scripts and Notebooks). In most cases, replacement with innerCellValues actually corrects small, undetected"off by one" indexing errors.

[mhvwerts]

@simulkade I have the feeling that in most cases, the ghost cells should not be set/copied explicitly, but rather be (re)calculated from the inner cell values and the governing BCs. Am I right?

(This is actually what instantiation of a CellVariable has been doing since the beginning: take the inner cell values and compute the ghost cell values from the supplied BCs. If no BCs were supplied when creating a CellVariable, default no-flux / zero Neumann BCs were already created, even when the CellVariable is not a solution variable and no BCs are needed)

[mhvwerts]

Only the inner cells should be involved in the arithmetic of CellVariable objects. Inclusion of the ghost cells may lead to 'divide by zero' warnings/errors which cannot be controlled by the user. Ghost cells should be re-calculated after doing the the arithmetic on the inner cells.

[mhvwerts]

Personally, I will probably continue to do most CellVariable-related arithmetic operations via the .innerCellValues property, with - for now - explicit updating of the ghost cells via apply_BCs() if necessary. We might use direct CellVariable arithmetic in the future.

[mhvwerts]

All pytest tests (scripts, and Notebooks still via pytest-notebook) complete successfully, so this may be merged into main as far as I am concerned.

This PR request does not break any existing code that runs with 0.3.0, except for the renaming of internalCellValues property to innerCellValues. It only adds automated handling of BC information by CellVariable objects.

[simulkade]

Perhaps CellVariable.value should become an internal (private) property CellVariable._value that the user is not supposed to touch, because it becomes confusing with the ghost cells included.

Most user cases are better handled via the CellVariable.innerCellValues property, followed by an explicit call to CellVariable.apply_BCs() where necessary.

I agree for most cases. However, when the user, for example, needs to assign spatially varying initial values to a cell, it can become useful. Although, it can always be done by creating a new CellVariable. Thinking about more advanced problems I solve with PyFVTool, e.g., sequentially non-iterative reactive transport, I almost always create a new cell variable by assigning updated inner cell values!
I agree with the innerCellValues property. As you said, it has two advantages of solving the unintuitive indexing and bringing more consistency to ghost cell values.

[simulkade]

Only the inner cells should be involved in the arithmetic of CellVariable objects. Inclusion of the ghost cells may lead to 'divide by zero' warnings/errors which cannot be controlled by the user. Ghost cells should be re-calculated after doing the the arithmetic on the inner cells.

Totally agree. The previous implementation has been out of laziness rather than a good design.

[simulkade]

Personally, I will probably continue to do most CellVariable-related arithmetic operations via the .innerCellValues property, with - for now - explicit updating of the ghost cells via apply_BCs() if necessary. We might use direct CellVariable arithmetic in the future.

I'm a middle of a switch to another department, which has made my schedule way too busy. I will push myself during next holidays to bring some of the advanced cases to the new PyFVTool. The Buckley-Leverett problem is one of the best ones that takes advantage of almost every functionality of the solver.

[mhvwerts]

Thanks for the swift reaction. Good luck with the switch. Such a project can indeed be time and energy consuming. Take care!

Concerning PyFVTool, I think we have a arrived at a stage where indeed it is best to add more and more examples of increasing complexity. We should use the 'new' workflow with embedded BCs and the new solvePDE() where possible, only addressing innerCellValues (or simply create a new CellVariable at each step as you do, which is equivalent to calling apply_BCs() after an update of the innerCellValues).

There are many minor or cosmetic changes that may still be considered (e.g., rename innerCellValues once again to something simpler..., or automatically calling apply_BCs() when updating the innerCellValues), but I do not expect that big code changes will be needed. It will mainly be changes to the API naming and adding convenient short-cuts (e.g. add a CellVariable.centers shortcut to CellVariable.domain.cellcenters and the like).

For now, let's enjoy the new PyFVTool version 0.3.1!