Hi Jonathan, I'm about to submit a pull request that changes (hopefully fixes...) the analytical solution in the Robin BC example. I'm attaching/pasting <https://gist.github.com/raybsmith/4b89f247cf1f056043fa> a Sage <http://www.sagemath.org/> script verifying the calculus/algebra done in the commit, for your reference. It can be run by calling "sage fipy_robin.sage" at the command line. I'm emailing because GitHub didn't seem to want to let me attach anything but an image file to the comment. Hope it may be helpful.
Best, Ray On Fri, Oct 24, 2014 at 1:55 PM, Guyer, Jonathan E. Dr. < jonathan.gu...@nist.gov> wrote: > Roberto - > > Thank you very much for your contributions. Presumably you're not new to > numerical solutions to PDEs, but just to FiPy? > > On Oct 24, 2014, at 8:32 AM, Roberto Rocca <roberto.ro...@cern.ch> wrote: > > > Dear Evan, dear all, > > > > This is my first post here, so first of all I'd like to thank the > developers for sharing their great work. > > > > I've only been using FiPy for a few weeks, I'm not an expert, but I'd > like to share my experience which I think is relevant to this post. > > > > As a preliminary step before tackling more complex problems, I've been > testing a simple model of thermal diffusion with internal heat generation, > both in steady-state and transient, to be benchmarked with analytic > solutions or with a commercial FEM software. Like the original poster, I > also had difficulties implementing a convection boundary condition, or in > general a BC that is a function of the variable solved for. These are my > main remarks: > > > > 1. It seems to me that when such a BC is implemented using for example > phi.faceGrad.constrain(-hh * (phi.faceValue - T0)/k, mesh.facesRight), the > value of the variable is not updated properly throughout the calculation. > At least in the cases I was testing. The workaround I came up with, after > much trial and error, is to place this call within a sweeping loop, which > seems to work and yield the expected results. This could be possibly a FiPy > issue. > > There is something wrong here, for sure. There are some subtleties that > can cause this problem (but aren't the whole story). > > Evan's code says > > phi.faceGrad.constrain(-hTopSurface*(phi.faceValue()-tinf)/k, topSurface) > > There are a couple of problems with this: > > * The constraint is not dynamic > > `phi.faceValue` is the property of `phi` that interpolates from cell > centers to face centers. > It is a FaceVariable that updates as `phi` changes. > > `phi.getFaceValue()` is the old-style method of `phi` that does the same > thing. > > `phi.faceValue()` (with parentheses) is the instantaneous value of > `phi.faceValue`. > This is equivalent to writing `phi.faceValue.value`. > It is a FaceVariable, but its value is static and does not change with > `phi`. > > * The constraint assigns a scalar to a vector. I *think* what this does is > assign the same scalar value to the x, y, z components of the faceGrad, > which is certainly not what's desired. > > To properly assign the Robin value to the normal component of > `phi.faceGrad`, use > > phi.faceGrad.constrain((-hTopSurface*(phi.faceValue-tinf)/k) * > mesh2.faceNormals, where=topSurface) > > This change doesn't fix things, though. As you've found, the value of the > constraint still doesn't update for some reason. > > > Our examples/convection/robin.py attempts to solve the second problem by > writing > > C.faceGrad.constrain([-P + P * C.faceValue], mesh.facesLeft) > > in an attempt to make the constraint value a vector, but even though `-P + > P * C.faceValue` is dynamic, the list `[-P + P * C.faceValue]` is not. When > I "fix" this, though, the numerical and analytical solutions completely > disagree. > > > > 2. Even in presence of a Dirichlet condition on part of the boundary, > when a convection boundary is implemented the solution could diverge. I've > been using the standard solver, which with pysparse is a PCG solver. This > can be dealt with with a proper choice of a coefficient of under relaxation. > > Good suggestion. What underRelaxation did you use? I had to set it to > 0.999 to get it to converge at any reasonable rate at all. > > > 3. In Evan's case there are no source or transient terms and in this > case I am not so sure of the proper choice for the coefficient of the > diffusion term; as the coefficient is constant, it could be canceled out of > the equation which reduces to a Laplace equation. I would personally put > either 1.0 or the thermal conductivity k, rather then the thermal > diffusivity k/rho.cp. Maybe someone can comment on this. > > My remark is that with Evan's figures, D ~7e-5 and a tolerance 1e-5, the > calculation would stop at the first iteration without reaching convergence. > It seems to work with a tolerance 1e-9. > > Also a very good observation. > > > 4. With the above changes to the code I get a global temperature range > from 30.8C to 31.4C, which is still lower than the reported expected result. > > The values I get depend on what I choose for underRelaxation and for > restol, but the range I get is similar to what you find. > > > I hope this helps. > > Very much. Thanks for your contributions. > > I filed an issue about all of this at > https://github.com/usnistgov/fipy/issues/426. I don't have any further > suggestions for the moment and have run out of time to look at it for > awhile. Hopefully Daniel can take a look when he returns next week(?). > > > _______________________________________________ > fipy mailing list > fipy@nist.gov > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] >
fipy_robin.sage
Description: Binary data
_______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
