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 ]
