Applying Two Different Diffusion Coefficients at Single FV Face

Hello All,

We have a uniformly spaced 1D mesh of length 2.0, wherein a standard elliptic
diffusion equation (\nabla. (D \nabla \phi) = f is being solved.

The diffusion coefficient, D, changes sharply at x = 1.0. Although this is
quite similar to examples.diffusion.mesh1D,  there is a key difference:

We are seeking an instantaneous change in the coefficient value at that face,
i.e. something of the effect that on one side of the face (i.e. from x = 0 up
to and including x = 1.0) diffusion coefficient D has a value D1.  Then there
is an abrupt change in diffusion coefficient from x >= 1.0 (very notable,
including the face at x = 1.0) and for the remainder of the mesh, diffusion
coefficient D has a value D2.  Currently, we have to either decide on whether
the single face at x = 1.0 gets the value D1 or D2. However, when we run with
either of these cases, the results are not accurate (when compared against a
benchmark from a commercial FV solver). We are getting a smooth interpolation
across the face , rather than the sharp change correctly predicted by the
benchmarking code.

With the code below, we have only been able to set the coefficient at that
middle face to one of either D1 or D2. This doesn't result in the desired
instantaneous change in the coefficient value, but instead, the new coefficient
value only applies from the NEXT face (which is a distance dx away) from this
critical interior boundary.

Coefficient = FaceVariable(mesh = 1D_mesh)

# Define facevariable on mesh
Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] <= 1.0))
#
Set coefficient values in 1st half of mesh
Coefficient.setValue(D2, where=((1.0 < 1D_mesh.faceCenters[0]) &
(1D_mesh.faceCenters[0] <= 2.0)))                    # Set coefficient values
in 2nd half of mesh

An alternative with the inequalities adjusted, but that still only permits a
single coefficient at that face:
Coefficient = FaceVariable(mesh = 1D_mesh)

# Define facevariable on mesh
Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] < 1.0))
#
Set coefficient values in 1st half of mesh
Coefficient.setValue(D2, where=((1.0 <= 1D_mesh.faceCenters[0]) &
(1D_mesh.faceCenters[0] <= 2.0)))                  # Set coefficient values in
2nd half of mesh

Sincerely,

-          Ian

P.S. Daniel, thank you very much for your previous reply on sweep speeds! It

_______________________________________________