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 <[email protected]> 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 [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
