Hi Dr. Guyer, I updated my code with your changes. It now matches the analytical solution pretty well, after several grid refinements. This is of course not a comprehensive test.
The Python script: https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/ConvectionTestProblem2D_01_20181215Version.py A report showing agreement between FiPy and analytical solution: https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/Report20181216.pdf The README.md of my repo has been updated to provide a better explanation of what is going on. I apologize for not learning git/github well enough to have github properly track changes between us as we went along. Instead, I put comments at the top of my scripts to document how I was looking at your code and documentation. Thanks On Fri, Dec 14, 2018 at 8:20 PM Guyer, Jonathan E. Dr. (Fed) via fipy < fipy@nist.gov> wrote: > Drew - > > I had found the same sign error in the denominator. You're also correct > that the units didn't agree. Multiplying by the face areas is redundant to > taking the divergence. > > I've updated our notes: > > > https://github.com/usnistgov/fipy/blob/dd3420fb71884d74850051ad2280bff525301824/documentation/USAGE.rst > > and your latest code: > > > https://github.com/guyer/Convection2DFiPyExample01/commit/ca59e8955319545f9965705c8fadfcbb5abd5951 > > With these changes, the results look much more plausible. > > - Jon > > > On Dec 12, 2018, at 1:22 PM, Drew Davidson <davidson...@gmail.com> > wrote: > > > > Hi Dr. Guyer, > > > > I looked over your more recent documentation: > > > https://github.com/usnistgov/fipy/blob/70b72b7abb267c85ada47886b8b3573e1819fffc/documentation/USAGE.rst > > > > I am embarrassed I did not understand how to choose values for a and b > in the Robin condition before. > > > > I also cloned your repo to my local computer and ran it: > > https://github.com/guyer/Convection2DFiPyExample01 > > > > When I run the code in your repo, the viewer does show a variation in T > with respect to r, but the size of that variation is extremely tiny. > Essentially, the solution is pinned at the initial temperature. The > variation also has the wrong sign (solid object is warming up), as the > ambient air is colder than the solid object (T_infinity < T_initial), and > the solid object should be cooling down. Maybe this is due to confusion I > sowed by having var be T minus T_infinity; sorry. > > > > I tried another script in my repo: > > > > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/ConvectionTestProblem2D_01_20181211Version.py > > > > This script is upgraded in that now, the log file contains information > about the solution (max and min temperatures). I obtain max and min > temperatures using var.value.min() and var.value.max() in order to rule out > goof-ups in my get_solution_along_ray function. > > > > The solution variable var is now exactly the temperature, in order to > get rid of the previous source of confusion in which var was T minus > T_infinity. > > > > The solution in my script is still pinned at the initial temperature. > This persists even if I crank up the convection coefficient to a huge value > (keep multiplying it by 1000 over and over and still seeing solution pinned > at initial condition). > > > > I made handwritten notes on top of a pdf I made via Kile Editor of your > more recent documentation: > > > > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/RobinBoundaryConditionsFiPyUsage20181211_Annot_CTTC.pdf > > > > It seems to me that the units of the source terms should match the units > of the terms in the partial differential equation for heat conduction. I > thought that is how it worked from the FiPy examples, but I haven’t gone > back and made sure. I find that for the first source term RobinCoeff*g, the > units don’t seem to match the units of the terms in the PDE for heat > conduction (I didn’t check the 2nd source term since it should be the same > situation). I did this in kind of a hurry so maybe I am goofing up. I have > not studied the finite volume method, and am just going on a hunch. > > > > It also seemed to me that there might be a minus sign missing in the > RobinCoeff. The solution I get is still pinned at initial condition > regardless of a minus sign or no minus sign there. > > > > Thanks > > > > On Fri, Dec 7, 2018 at 6:05 PM Guyer, Jonathan E. Dr. (Fed) via fipy < > fipy@nist.gov> wrote: > > Drew - > > > > > > Thanks to your notes, I found a couple of errors in our Robin > discussion. As you noted, the units don't work for the conversion of the > divergence of gamma grad phi to the sum over faces (midway on page 1 of > your notes). The left hand side should be integrated over volume. That line > is expressing the discretization of the divergence theorem. > > > > I also discovered that the .divergence operator in FiPy doesn't work > reliably on scalars. There's no real reason it should, but I thought I had > convinced myself that it did. As a result, everything in the Robin terms > needs to be multiplied by the surface normal > > > > Those changes to the FiPy documentation are [here]( > https://github.com/usnistgov/fipy/pull/615) for now. > > > > Beyond addition that, there are some changes necessary to put your > boundary condition into Robin form. The point of the Robin condition is > that it ties the gradient of the field to the value of the field. > > So, g isn't h(T - T_inf) / (-k); it's just -h T_inf / (-k). > > Likewise, \vec{a} isn't zero. Rather, \hat{n}\cdot\vec{a} = h/(-k). > > > > Here are the changes I made to your script to reflect these changes: > > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/compare/master...guyer:master > > > > I think I got your normalized temperature accounted for properly, and > put the thermal diffusivity and convection coefficient in the proper > places, but it's worth checking through. > > > > With these changes, the temperature field is now rotationally symmetric > (it wasn't before, which is why I had to multiply by the normal before > taking the divergence). > > > > Heat seems to be fluxing in from the outside, so I probably have the > sign wrong. > > > > I don't have Octave, so I have no idea how this compares with your > analytical solution. > > > > - Jon > > > > > On Dec 4, 2018, at 2:34 PM, Drew Davidson <davidson...@gmail.com> > wrote: > > > > > > Hi Dr. Guyer, > > > > > > First I tried getting rid of the square brackets in > ConvectionTestProblem2D_01.py (commit ‘changed square brackets to > parenthesis in convection BC, but get same…’), but results are the same > (still wrong). > > > > > > Next: > > > As you directed, I took a look at: > > > > https://github.com/usnistgov/fipy/blob/develop/documentation/USAGE.rst#applying-robin-boundary-conditions > > > > > > Firefox was showing me raw latex rather than human-readable equations > at that web address, so I made a version I could read using the kile editor: > > > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/RobinBoundaryConditionsFiPyUsage.pdf > > > > > > That material is rather challenging for me. I took a stab at it, > resulting in: > > > commit: ‘made a first attempt at a formulation in view of suggestions > by fipy …’ > > > > > > ConvectionTestProblem2D_01_2ndTry.py > > > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/fb47960548650c994eb0c6f990e0db297566e174/ConvectionTestProblem2D_01_2ndTry.py > > > > > > a few handwritten notes indicating what I was thinking: > > > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/blob/master/RobinBoundaryConditionsFiPyUsageConvectionBC.pdf > > > > > > Here is a bit of the code: > > > > > > #ref: > https://github.com/usnistgov/fipy/blob/develop/documentation/USAGE.rst#applying-robin-boundary-conditions > > > > > > #warning: avoid confusion between convectionCoeff in that fipy > documentation, which refers to terms involving "a", and convectionCoeff > here, which refers to a heat transfer convection coefficient at a boundary > > > > > > Gamma0=D_thermal > > > > > > Gamma = FaceVariable(mesh=mesh, value=Gamma0) > > > > > > mask=surfaceFaces > > > > > > Gamma.setValue(0., where=mask) > > > > > > dPf = FaceVariable(mesh=mesh, value=mesh._faceToCellDistanceRatio * > mesh.cellDistanceVectors) > > > > > > Af = FaceVariable(mesh=mesh, value=mesh._faceAreas) > > > > > > #RobinCoeff = (mask * Gamma0 * Af / (dPf.dot(a) + b)).divergence #a is > zero in our case > > > > > > b=1. > > > > > > RobinCoeff = (mask * Gamma0 * Af / b).divergence #a is zero in our case > > > > > > #eq = (TransientTerm() == DiffusionTerm(coeff=Gamma) + RobinCoeff * g > - ImplicitSourceTerm(coeff=RobinCoeff * mesh.faceNormals.dot(a))) #a is > zero in our case > > > > > > # g in this formulation is -convectionCoeff/k*var, where > var=T-T_infinity > > > > > > eq = (TransientTerm() == DiffusionTerm(coeff=Gamma) + > ImplicitSourceTerm(RobinCoeff * -convectionCoeff/k)) > > > > > > > > > Now the solution variable remains stuck at the initial condition, as > if the boundary condition is not being applied. I am sort of out of my > depth at this point. I was guessing that an ImplicitSourceTerm was the > right thing to do, since 'g' in the Robin condition depends on the solution > variable. I did mess around a little in IPython seeing if any terms are > coming out all zeros. > > > > > > Again, I put everything at > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01. > > > > > > Thanks > > > > > > On Mon, Dec 3, 2018 at 4:21 PM Guyer, Jonathan E. Dr. (Fed) via fipy < > fipy@nist.gov> wrote: > > > Drew - > > > > > > Apologies for the delayed reply. > > > > > > There are a couple of issues here: > > > > > > `-convectionCoeff/k*(var.faceValue-T_infinity)` describes a FiPy > FaceVariable. > > > > > > `[-convectionCoeff/k*(var.faceValue-T_infinity)]` is a single element > Python list that contains a FiPy FaceVariable. > > > > > > Multiplying that list by other elements has rather unpredictable > results. > > > > > > In short, get rid of the square brackets. We do our best to treat > Python lists like FiPy vector fields, but there's only so much we can do. A > list that holds a numpy array is not a vector field. > > > > > > > > > Beyond that, what you've described looks like a Robin boundary > condition to me. Our best recommendation for Robin conditions is covered at: > > > > > > > https://github.com/usnistgov/fipy/blob/develop/documentation/USAGE.rst#applying-robin-boundary-conditions > > > > > > Please don't hesitate to ask for more clarification if this doesn't > get you where you need. > > > > > > - Jon > > > > > > > > > > > > > On Nov 16, 2018, at 11:55 PM, Drew Davidson <davidson...@gmail.com> > wrote: > > > > > > > > Hello, > > > > > > > > I am stuck in how to correctly apply a simple convection boundary > condition in FiPy, in the context of simple transient heat conduction in 2D. > > > > > > > > Is this correct: > > > > > var.faceGrad.constrain([-convectionCoeff/k*(var.faceValue-T_infinity)]*mesh.faceNormals,where=surfaceFaces) > > > > > > > > I have a 2D mesh generated in gmsh. The convection boundary > condition is applied to a curved boundary. > > > > > > > > The code and results are found at: > > > > > https://github.com/cashTangoTangoCash/Convection2DFiPyExample01/tree/master > > > > > > > > The project appears to have a heap of files, but It’s a really just > a simple 2D problem with a comparison to an analytical solution. The basic > script is ConvectionTestProblem2D_01.py. Report20181116.pdf shows current > results, which are clearly wrong. > > > > > > > > Am I correctly applying the convection boundary condition in terms > of the FiPy syntax/language for this 2D problem with a gmsh mesh and a > curved boundary? > > > > > > > > Thanks > > > > _______________________________________________ > > > > 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 ] > > > _______________________________________________ > 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 ]
