RE: Applying Two Different Diffusion Coefficients at Single FV Face

Hi Raymond,

Thanks a lot for your inputs. We have been working along the same lines that
you have suggested.

There is a slightly more complication that wasn’t fully discussed
(intentionally, to keep it simple, and to find out if any other options
existed).

Assuming that we are working on a normalised co-ordinate and the two bulk
regions within the domain have non-uniform properties.  In this situation, we
have to normalise the PDEs too.

Thus, the diffusion coefficient becomes D_region/L_region  , where region can
be either region 1 or region 2. , and L_region denotes the original length of
the region.

Now, although it’s a perfectly straightforward idea to define D in the cell
centers and let fiPy handle the harmonic facevalue computation,  I don’t know
if it’s meaningful to compute harmonicfacevalue of a length (which appears in
the denominator).

What do you suggest for L_region ?  Define that at the cell-centers again, and
use arithmeticfacevalue , or stick with harmonicfacevalue ?

Best Regards

Krishna & Ian.

From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of Raymond
Smith
Sent: Friday, September 16, 2016 7:18 PM
To: fipy@nist.gov
Subject: Re: Applying Two Different Diffusion Coefficients at Single FV Face

A side note, here it may actually be more convenient to think about the D in
terms of the volumes, so you could also define it as a cell variable with
values specified as D1 when x<=1 and D2 when x>1, then use the
harmonicFaceValue attribute when you put it into the governing equation to have
FiPy do this calculation for you and avoid awkward use of the mesh spacing in
the specification of the variable values.

On Fri, Sep 16, 2016 at 11:12 AM, Raymond Smith
<smit...@mit.edu<mailto:smit...@mit.edu>> wrote:
Hi, Ian.
I don't think there is such a thing as having two different flux coefficients
at the same face for the same governing PDE. The flux through a given face is
calculated by the coefficient at that face times some approximation of the
gradient in a field variable at that face, like
both evaluated at the face, so there is only one coefficient in that
expression. And in the FiPy FV approach D doesn't come into play at all except
at the faces, so it can't change immediately to the right or left of a face but
only from one face to the next.
That said, perhaps you could try following a common approach used when a FV
interface is directly between two bulk regions with different continuum
properties -- use some sort of a mean of the continuum transport coefficient in
the adjacent volumes as the approximation for the coefficient at the face. I
would suggest the harmonic mean (one simple reason: it preserves zero flux
through the face if either of the adjacent cells has a zero coefficient). You
could try:
Dint = 2*D1*D2/(D1+D2)
Coefficient.setValue(Dint, where=((1.0-dx/2 < 1D_mesh.faceCenters[0]) &
(1D_mesh.faceCenters[0] < 1.0+dx/2)))
If you do that, you should see a change in the slope of the field variable
right near the face between the two bulk regions, which seems to be what you're
missing. Of course, ideally your mesh is fine enough that the dx is
insignificant in terms of the solution, but this is a common approach to get a
slightly more reasonable estimate for what happens between two regions.

Best,
Ray

On Fri, Sep 16, 2016 at 10:44 AM, Campbell, Ian
<i.campbel...@imperial.ac.uk<mailto:i.campbel...@imperial.ac.uk>> wrote:
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

_______________________________________________
fipy mailing list
fipy@nist.gov<mailto:fipy@nist.gov>
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]


_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]