# Re: Applying Two Different Diffusion Coefficients at Single FV Face

Yes!

> On Sep 16, 2016, at 2:18 PM, Raymond Smith <smit...@mit.edu> wrote:
>
> 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> 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
> D * grad(c),
> 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>
> 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
> was very helpful.
>
>
>
>
> _______________________________________________
> fipy mailing list
> 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 ]


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