Wow, thank you for both of those responses... they were extremely helpful. And thanks for wading through my beginner-level code. I have been experimenting with things to see if I could get FiPy to update faceValue during solve which is why I called it with faceValue() - but I forgot to change it back before posting. I was confused about faceGrad orientation as well, thanks for correcting me.
I saw someone else having trouble with this boundary condition in an older thread. There was a mention to try using a divergence term to set the robin condition. This now works outside of the sweeping loop for me as long as I am using underRelaxation. Here is the equation: eq = DiffusionTerm(coeff=D) + (-hTopSurface*(phi.faceValue-tinf)/k * mesh2.faceNormals * mesh2.facesBack).divergence After following your advice I was getting similar numbers, ~31.32C maximum value. This is only ~6-9% error from my accepted value. The error was because of the way I set up the source term by selecting cell centers within the area I specified. If I set the mesh in x and y to numbers divisible by 4 I get closer to the correct answer. Is there anyway to make the source "bounds" register as key mesh points in the mesh using grid3D? I haven't found anything like that in the documentation. Suggestions are welcome too. I can try to implement something. I do want to eventually parametrize this with multiple sources. My goal is to create a larger python tool that uses this calculation as a subroutine. Both the neumann and the robin BCs will be functions of space (and eventually time.) And they will iteratively update according to a simpler iterative routine (which calculates flow, then htc from nusselt correlations.) I do want to be able to solve this problem relatively fast as I will want to sweep parameters to optimize a design. I am also only able to use high values of underRelaxation (0.999,0.9999). The problem I have now is that i'm not sure what the optimal underRelaxation and tolerances are to get a fast (second or two) convergence. I found one way is to change the increase (lower) the underrelaxation as I get res closer to restol. Suggestions are welcome again! It is still relatively slow... but it works now!!! Thank you very much!!! Evan On Fri, Oct 24, 2014 at 10:55 AM, 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 mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
