General boundary conditions

Question

I’d appreciate help on setting general boundary conditions, -grad(y) + g(y) = 0 where g is some function of the unknown y. Here’s a simple 1D example that I can’t get to work:

N=3
h=1./(float(N)-1.)

mesh = Grid1D(nx=N, dx=h)

c=CellVariable(mesh=mesh,value=0.5)

## Dirichlet boundary conditions
#c.constrain(2., mesh.facesLeft)
#c.constrain(1., mesh.facesRight)

## Neumann boundary conditions
c.faceGrad.constrain(-1, where=mesh.facesLeft)
c.faceGrad.constrain( -c.faceValue , where=mesh.facesRight)

Eq = DiffusionTerm(coeff=1.0)
Eq.cacheMatrix()
Eq.cacheRHSvector()
Eq.solve(var=c)
m = Eq.matrix.numpyArray
b = Eq.RHSvector

This code won’t solve but I do get to see the matrix and RHS:

m= array([[-2.,  2.,  0.],
          [ 2., -4.,  2.],
          [ 0.,  2., -2.]])

b= array([-1. ,  0. ,  0.5])

The matrix, m, is clearly singular because the source term was not included in the last line. Any suggestions on how to include it?

Answer 1

[edit: add derivation and demonstration of 1st order implementation]

There are known issues with general boundary conditions.

It is possible to implement such a boundary condition as a source. Using the discretization of the DiffusionTerm $\sum_f D_f (n\cdot\nabla(y))_f A_f$ we treat $f=R$ as a special case and substitute of the desired boundary condition -n\cdot\nabla(y) - y = 0. We accomplish this by zeroing D_(f=R) in the DiffusionTerm

c.faceGrad.constrain([-1], where=mesh.facesLeft)

D = 1.
Dface = FaceVariable(mesh=mesh, value=D)
Dface.setValue(0., where=mesh.facesRight)

and then adding an implicit source D_(f=R) (n\cdot\nabla(y))_(f=R) A_(f=R) or D_(f=R) (-y)_(f=R) A_(f=R). Sources are defined at cell centers, so we locate the cell adjacent to $f=R$

mask_coeff = (mesh.facesRight * mesh.faceNormals).divergence

and then add the source

Af = mesh._faceAreas[mesh.facesRight.value][0]
Eq = DiffusionTerm(coeff=Dface) - ImplicitSourceTerm(coeff=mask_coeff * D * Af)

This treatment is only 0th order accurate as the ImplicitSourceTerm operates on the value at the cell center, whereas the boundary condition is defined at the adjacent face center.

We can make the boundary condition 1st order accurate in space by projecting the cell value to the face along the gradient from the boundary condition: y_(f=R) ~ y_P + n\cdot\nabla(y)_(f=R) dPR where y_P is the value of y at the cell center closest to f=R and dPR is the distance from point P to face R.

Thus the boundary condition -n\cdot\nabla(y)_(f=R) - y_(f=R) = 0 becomes -n\cdot\nabla(y)_(f=R) - (y_P + n\cdot\nabla(y)_(f=R) dPR) = 0, which we can solve for n\cdot\nabla(y)_(f=R) = -y_P / (1 + dPR). The implicit source corresponding to the boundary for the DiffusionTerm is thus D_(f=R) (-y_P / (1 + dPR)) A_(f=R)

dPR = mesh._cellDistances[mesh.facesRight.value][0]
Af = mesh._faceAreas[mesh.facesRight.value][0]
Eq = DiffusionTerm(coeff=Dface) - ImplicitSourceTerm(coeff=mask_coeff * D * Af / (1 + dPR))

This was discussed last summer on the FiPy mailing list, starting at https://www.mail-archive.com/fipy@nist.gov/msg03671.html. Yes, it's all pretty clumsy right now.