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 ]
