Hi Raymond, We moved FiPy to Github very recently and I was wondering if you would like to submit this patch as a pull-request so you can get credit for your changes.
Thanks, Daniel On Mon, Aug 18, 2014 at 4:48 PM, Raymond Smith <[email protected]> wrote: > Also, I'm not sure if this is what you meant by the doctests, but I have > added something that may be along the lines of what you were talking about. > > I attached a diff to the ticket in case the git path doesn't work,. > > Ray > > diff --git a/examples/diffusion/mesh1D.py b/examples/diffusion/mesh1D.py > index 71f2da4..9e6e705 100755 > --- a/examples/diffusion/mesh1D.py > +++ b/examples/diffusion/mesh1D.py > @@ -450,6 +450,81 @@ And finally, we can plot the result > > > ------ > > +Note that for problems involving heat transfer and other similar > +conservation equations, it is important to ensure that we begin with > +the correct form of the equation. For example, for heat transfer with > +:math:`\phi` representing the temperature, > + > +.. math:: > + \frac{\partial \rho \hat{C}_p \phi}{\partial t} = \nabla \cdot [ k > \nabla \phi ]. > + > +With constant and uniform density :math:`\rho`, heat capacity > :math:`\hat{C}_p` > +and thermal conductivity :math:`k`, this is often written like Eq. > +:eq:`eq:diffusion:mesh1D:constantD`, but replacing :math:`D` with > :math:`\alpha = > +\frac{k}{\rho \hat{C}_p}`. However, when these parameters vary either in > position > +or time, it is important to be careful with the form of the equation > used. For > +example, if :math:`k = 1` and > + > +.. math:: > + \rho \hat{C}_p = \begin{cases} > + 1& \text{for \( 0 < x < L / 4 \),} \\ > + 10& \text{for \( L / 4 \le x < 3 L / 4 \),} \\ > + 1& \text{for \( 3 L / 4 \le x < L \),} > + \end{cases}, > + > +then we have > + > +.. math:: > + \alpha = \begin{cases} > + 1& \text{for \( 0 < x < L / 4 \),} \\ > + 0.1& \text{for \( L / 4 \le x < 3 L / 4 \),} \\ > + 1& \text{for \( 3 L / 4 \le x < L \),} > + \end{cases}. > + > +However, using a ``DiffusionTerm`` with the same coefficient as that in > the > +section above is incorrect, as the steady state governing equation > reduces to > +:math:`0 = \nabla^2\phi`, which results in a linear profile in 1D, unlike > that > +for the case above with spatially varying diffusivity. Similar care must > be > +taken if there is time dependence in the parameters in transient problems. > + > +We can illustrate the differences with an example. We define field > +variables for the correct and incorrect solution > + > +>>> phiT = CellVariable(name="correct", mesh=mesh) > +>>> phiF = CellVariable(name="incorrect", mesh=mesh) > +>>> phiT.faceGrad.constrain([fluxRight], mesh.facesRight) > +>>> phiF.faceGrad.constrain([fluxRight], mesh.facesRight) > +>>> phiT.constrain(valueLeft, mesh.facesLeft) > +>>> phiF.constrain(valueLeft, mesh.facesLeft) > +>>> phiT.setValue(0) > +>>> phiF.setValue(0) > + > +The relevant parameters are > + > +>>> k = 1. > +>>> alpha_false = FaceVariable(mesh=mesh, value=1.0) > +>>> X = mesh.faceCenters[0] > +>>> alpha_false.setValue(0.1, where=(L / 4. <= X) & (X < 3. * L / 4.)) > +>>> eqF = 0 == DiffusionTerm(coeff=alpha_false) > +>>> eqT = 0 == DiffusionTerm(coeff=k) > +>>> eqF.solve(var=phiF) > +>>> eqT.solve(var=phiT) > + > +Comparing to the correct analytical solution, :math:`\phi = x` > + > +>>> x = mesh.cellCenters[0] > +>>> phiAnalytical.setValue(x) > +>>> print phiT.allclose(phiAnalytical, atol = 1e-8, rtol = 1e-8) # > doctest: +SCIPY > +1 > + > +and finally, plot > + > +>>> if __name__ == '__main__': > +... Viewer(vars=(phiT, phiF)).plot() > +... raw_input("Non-uniform thermal conductivity. Press <return> to > proceed...") > + > +------ > + > Often, the diffusivity is not only non-uniform, but also depends on > the value of the variable, such that > > > > > On Mon, Aug 18, 2014 at 4:35 PM, Raymond Smith <[email protected]> wrote: > >> Hi Jonathan, >> >> I pulled, branched, committed, and fetched/merged develop according to >> the instructions on the linked site, but my repo isn't publicly available >> online. When I do a >> >> $ git format-patch >> >> I get no output. My git log shows >> >> commit 6d571b83dc5dfdeb16cae51b016b69cb279d5450 >> Author: Raymond Smith <[email protected]> >> Date: Mon Aug 18 16:23:15 2014 -0400 >> >> Add heat transfer example to mesh1D. >> >> andgit status shows >> >> On branch >> ticket670-Remind_users_of_different_types_of_conservation_equations >> nothing to commit, working directory clean >> >> Sorry, I must be missing something. >> >> Ray >> >> >> On Mon, Aug 18, 2014 at 3:14 PM, Guyer, Jonathan E. Dr. < >> [email protected]> wrote: >> >>> Your patch looks a good start to me and I think after the non-uniform >>> diffusivity is the right place to put it. Can you add doctests of both >>> desired and undesired behavior? >>> >>> We expect to be migrating to github "soon", when pull requests will be >>> easier, but in the meantime an email patch is OK. Somewhat better would be >>> what we do in our own development: >>> >>> >>> http://www.ctcms.nist.gov/fipy/documentation/ADMINISTRATA.html#submit-branch-for-code-review >>> >>> In short, file a ticket on matforge and then attach either your `git >>> request-pull` or your `git format-patch` results to the ticket. This makes >>> it a bit easier to keep track and make sure we don't permanently lose any >>> patches (for instance, I'd test and merge your patch right now, but I just >>> moved to a new laptop and everything is broken). >>> >>> On Aug 18, 2014, at 2:00 PM, Raymond Smith <[email protected]> wrote: >>> >>> > Hi, Jonathan. >>> > >>> > I didn't know how to submit a pull request on MatForge, so here >>> (attached and pasted below) is a git diff of mesh1D.py. I added a small >>> section discussing the import of keeping the original governing equation in >>> mind when parameters vary with a specific example from heat transfer. I >>> didn't know where it would fit best in the example; I added it after the >>> non-uniform diffusivity section. If you have suggestions, I can edit. >>> > >>> > Ray >>> > >>> > >>> > >>> > diff --git a/examples/diffusion/mesh1D.py >>> b/examples/diffusion/mesh1D.py >>> > index 71f2da4..3c5209f 100755 >>> > --- a/examples/diffusion/mesh1D.py >>> > +++ b/examples/diffusion/mesh1D.py >>> > @@ -450,6 +450,45 @@ And finally, we can plot the result >>> > >>> > ------ >>> > >>> > +Note that for problems involving heat transfer and other similar >>> > +conservation equations, it is important to ensure that we begin with >>> > +the correct form of the equation. For example, for heat transfer with >>> > +:math:`\phi` representing the temperature, >>> > + >>> > +.. math:: >>> > + \frac{\partial \rho \hat{C}_p \phi}{\partial t} = \nabla \cdot [ >>> k \nabla \phi ]. >>> > + >>> > +With constant and uniform density :math:`\rho`, heat capacity >>> :math:`\hat{C}_p` >>> > +and thermal conductivity :math:`k`, this is often written like Eq. >>> > +:eq:`eq:diffusion:mesh1D:constantD`, but replacing :math:`D` with >>> :math:`\alpha = >>> > +\frac{k}{\rho \hat{C}_p}`. However, when these parameters vary either >>> in position >>> > +or time, it is important to be careful with the form of the equation >>> used. For >>> > +example, if :math:`k = 1` and >>> > + >>> > +.. math:: >>> > + \rho \hat{C}_p = \begin{cases} >>> > + 1& \text{for \( 0 < x < L / 4 \),} \\ >>> > + 10& \text{for \( L / 4 \le x < 3 L / 4 \),} \\ >>> > + 1& \text{for \( 3 L / 4 \le x < L \),} >>> > + \end{cases}, >>> > + >>> > +then we have >>> > + >>> > +.. math:: >>> > + \alpha = \begin{cases} >>> > + 1& \text{for \( 0 < x < L / 4 \),} \\ >>> > + 0.1& \text{for \( L / 4 \le x < 3 L / 4 \),} \\ >>> > + 1& \text{for \( 3 L / 4 \le x < L \),} >>> > + \end{cases}. >>> > + >>> > +However, using a ``DiffusionTerm`` with the same coefficient as that >>> in the >>> > +section above is incorrect, as the steady state governing equation >>> reduces to >>> > +:math:`0 = \nabla^2\phi`, which results in a linear profile in 1D, >>> unlike that >>> > +for the case above with spatially varying diffusivity. Similar care >>> must be >>> > +taken if there is time dependence in the parameters in transient >>> problems. >>> > + >>> > +------ >>> > + >>> > Often, the diffusivity is not only non-uniform, but also depends on >>> > the value of the variable, such that >>> > >>> > >>> > >>> > On Mon, Aug 18, 2014 at 10:54 AM, Guyer, Jonathan E. Dr. < >>> [email protected]> wrote: >>> > Conor and Raymond - >>> > >>> > Thank you both for posting your findings and interpretation of the >>> issue. I agree that this would be a useful issue to make clear in the >>> documentation. We would welcome a patch or pull request from either of you >>> to illustrate this situation in 'examples.diffusion.mesh1D'. >>> > >>> > >>> > On Aug 18, 2014, at 9:36 AM, Raymond Smith <[email protected]> wrote: >>> > >>> > > Hi Conor, >>> > > >>> > > Just to add to your observations, I'm guessing FiPy is "correct" >>> here in both situations, and as you noticed, the original input with a >>> thermal diffusivity is simply not the correct representation of your >>> physical situation. The actual conservation equation for thermal energy is >>> (neglecting convection and source terms) >>> > > >>> > > \frac{\partial rho * c_p * T}{\partial t} = \nabla\cdot(k\nabla T) >>> > > >>> > > So if you have spatially varying rho and c_p, you cannot sweep them >>> into the coefficient for the flux term, else they will appear inside FiPy's >>> divergence operator. They are often combined with k to form alpha for the >>> convenience of writing it like the mass conservation (diffusion) equation >>> (when chemical diffusivity is constant), but it's always important to start >>> from the correct conserved quantities. I would guess that if you run it >>> with any rho, c_p that are the same in the concrete and insulator, you will >>> get the correct result with the "alpha form" because, being both constant >>> and uniform, they can be swept into and out of either type of partial >>> derivative. >>> > > >>> > > Cheers, >>> > > Ray >>> > > >>> > > >>> > > On Mon, Aug 18, 2014 at 8:51 AM, Conor Fleming < >>> [email protected]> wrote: >>> > > Hi Kris, >>> > > >>> > > >>> > > I have identified the cause of this difference - I had not specified >>> the governing equation correctly in FiPy. Originally, I had prescribed the >>> heat equation as follows: >>> > > >>> > > >>> > > eqn = TransientTerm() == DiffusionTerm(coeff=alpha) >>> > > >>> > > >>> > > where the thermal diffusivity alpha is >>> > > >>> > > >>> > > alpha = k / (rho * c_p), >>> > > >>> > > >>> > > and FiPY gives an erroneous answer. Additionally, solving the steady >>> equation 'DiffusionTerm(coeff=alpha)==0' also yields an incorrect result >>> where (rho*c_p) has a value other than unity. I have realised that the >>> equation should be specified as >>> > > >>> > > >>> > > eqn = TransientTerm(coeff=C) == DiffusionTerm(coeff=k) >>> > > >>> > > >>> > > where the volumetric heat capacity is C = rho * c_p and k is the >>> thermal conductivity. For a steady case, this equation reduces to >>> 'DiffusionTerm(coeff=k)==0' and gives a correct result. >>> > > >>> > > >>> > > I have attached a figure with the updated comparison of FiPy and >>> TEMP/W for an insulated concrete slab, showing perfect agreement. >>> > > >>> > > >>> > > I hope this result is useful to other FiPy users. It may be helpful >>> to add a note to the documentation page 'examples.diffusion.mesh1D', >>> explaining that for some applications, e.g. the heat equation, it is >>> appropriate to separate the (thermal) diffusivity into two portions, which >>> act on the TransientTerm and DiffusionTerm respectively. >>> > > >>> > > >>> > > Many thanks, >>> > > >>> > > Conor >>> > > >>> > > >>> > > From: [email protected] [[email protected]] on behalf of >>> Conor Fleming [[email protected]] >>> > > Sent: 08 August 2014 17:24 >>> > > To: [email protected] >>> > > Subject: RE: Unexpected result, possibly wrong, result solving 1D >>> unsteady heat equation with spatially-varying diffusion coefficient >>> > > >>> > > Hi Kris, >>> > > >>> > > >>> > > Thank you for the prompt response. You are right - altering the >>> insulation conductivity in the FiPy model to k_ins=0.1 improves the >>> agreement. I have just checked TEMP/W again, and can confirm that >>> k_ins=0.212 in that model. Specific heat capacity and density match as >>> well. I will read up on FE and FD implementations generally to see if I can >>> spot any issues on that side. >>> > > >>> > > >>> > > Many thanks, >>> > > >>> > > Conor >>> > > >>> > > >>> > > From: [email protected] [[email protected]] on behalf of >>> Kris Kuhlman [[email protected]] >>> > > Sent: 08 August 2014 17:04 >>> > > To: [email protected] >>> > > Subject: Re: Unexpected result, possibly wrong, result solving 1D >>> unsteady heat equation with spatially-varying diffusion coefficient >>> > > >>> > > Conor, >>> > > >>> > > if you reduce the thermal conductivity in the insulation to about >>> 0.1, the fipy solution looks about like the other model (the knee in T is >>> about at 400 degrees C). Is there an issue with how your compute or >>> specify this in fipy or the other model? >>> > > >>> > > Kris >>> > > >>> > > >>> > > On Fri, Aug 8, 2014 at 9:35 AM, Conor Fleming < >>> [email protected]> wrote: >>> > > Hi, >>> > > >>> > > I am using FiPy to determine the depth of heat penetration into >>> concrete structures due to fire over a certain period of time. I am solving >>> the unsteady heat equation on a 1D grid, and modelling various scenarios, >>> e.g. time-dependent temperature boundary condition, temperature-dependent >>> diffusion coefficient. For these cases, the model compares very well to >>> results from other solvers (for example the commercial finite element >>> solver TEMP/W). However I am having trouble modelling problems with a >>> spatially-varying diffusion coefficient, in particular, a two layer model >>> representing a concrete wall covered with a layer of insulating material. >>> > > >>> > > I have attached a number of images to clarify the issue. The first, >>> 'FiPy_TEMPW_insulation_only.png' shows the temperature distribution in a >>> single material (the insulation) with constant diffusivity, D_ins, when a >>> constant temperature of 1200 C is applied at one boundary for 3 hours. The >>> FiPy result agrees very well with the analytical solution >>> > > >>> > > 1 - erf(x/2(sqrt(D_ins t))*1200 +20, >>> > > >>> > > taken from the examples.diffusion.mesh1D example (scaled and shifted >>> appropriately) and with a numerical solution calculated using TEMP/W (using >>> the same spatial and temporal discretisation, material properties and >>> boundary conditions). The three results agree well, showing that the FiPy >>> model is performing as expected. >>> > > >>> > > The second figure, 'FiPy_TEMPW_insulated_concrete.png', presents the >>> temperature distribution through an insulated concrete wall (where for >>> simplicity the concrete is also modelled with constant diffusivity, D_con) >>> for the same surface temperature and period. There is now a considerable >>> difference between the FiPy and TEMP/W predictions. The FiPy model predicts >>> a lower temperature gradient in the insulation layer, which leads to higher >>> temperatures throughout the domain. >>> > > >>> > > I am confident that the TEMP/W result is accurate, as it agrees >>> perfectly with a simple explicit finite difference solution coded in >>> FORTRAN. I have tried to identify any coding errors I have made in my FiPy >>> script. I am aware that diffusion terms are solved at the grid faces, so >>> when the diffusion coefficient is a function of a CellVariable, an >>> appropriate interpolation scheme must be used to obtain sensible face >>> values. However, in my case the diffusion coefficients, D_ins and D_conc, >>> are created as FaceVariables and assigned constant values. I have also >>> examined the effects of space and time discretisation, implicit and >>> explicit DiffusionTerm, multiple sweeps per timestep, but these have made >>> no significant difference. >>> > > >>> > > I would be very interested to hear anyone's opinion on what I might >>> be doing wrong here. Also, does anyone think it is possible for FiPy to >>> deliver an accurate result, and for the finite difference and finite volume >>> solvers to be wrong? Below this email I have written out a minimal working >>> example, which reproduces the 'FiPy' curve from figure >>> 'FiPy_TEMPW_insulated_concrete.png'. >>> > > >>> > > Many thanks, >>> > > Conor >>> > > >>> > > -- >>> > > Conor Fleming >>> > > Research student, Civil Engineering Group >>> > > Dept. of Engineering Science, University of Oxford >>> > > >>> > > ################################################################### >>> > > # FiPy script to solve 1D heat equation for two-layer material >>> > > # >>> > > import fipy as fp >>> > > >>> > > >>> > > nx = 45 >>> > > dx = 0.009 # grid spacing, m >>> > > dt = 20. # timestep, s >>> > > >>> > > >>> > > mesh = fp.Grid1D(nx=nx, dx=dx) # define 1D grid >>> > > >>> > > >>> > > # define temperature variable, phi >>> > > phi = fp.CellVariable(name="Fipy", mesh=mesh, value=20.) >>> > > >>> > > >>> > > # insulation thermal properties >>> > > thick_ins = 0.027 # insulation thickness >>> > > k_ins = 0.212 >>> > > rho_ins = 900. >>> > > cp_ins = 1000. >>> > > D_ins = k_ins / (rho_ins * cp_ins) # insulation diffusivity >>> > > >>> > > >>> > > # concrete thermal properties >>> > > k_con = 1.5 >>> > > rho_con = 2300. >>> > > cp_con = 1100. >>> > > D_con = k_con / (rho_con * cp_con) # concrete diffusivity >>> > > >>> > > >>> > > valueLeft = 1200. # set temperature at edge of domain >>> > > phi.constrain(valueLeft, mesh.facesLeft) # apply boundary condition >>> > > >>> > > >>> > > # create diffusion coefficient as a FaceVariable >>> > > D = fp.FaceVariable(mesh=mesh, value=D_ins) >>> > > X = mesh.faceCenters.value[0] >>> > > D.setValue(D_con, where=X>thick_ins) # change diffusivity in >>> concrete region >>> > > >>> > > >>> > > # unsteady heat equation >>> > > eqn = fp.TransientTerm() == fp.DiffusionTerm(coeff=D) >>> > > >>> > > >>> > > # specify viewer >>> > > viewer = fp.Viewer(vars=phi, datamin=0., datamax=1200.) >>> > > >>> > > >>> > > # solve equation >>> > > t = 0. >>> > > while t < 10800: # simulate for 3 hours >>> > > t += dt # advance time >>> > > eqn.solve(var=phi, dt=dt) # solve equation >>> > > viewer.plot() # plot result >>> > > _______________________________________________ >>> > > fipy mailing list >>> > > [email protected] >>> > > http://www.ctcms.nist.gov/fipy >>> > > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy >>> ] >>> > > >>> > > >>> > > >>> > > _______________________________________________ >>> > > fipy mailing list >>> > > [email protected] >>> > > http://www.ctcms.nist.gov/fipy >>> > > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy >>> ] >>> > > >>> > > >>> > > _______________________________________________ >>> > > fipy mailing list >>> > > [email protected] >>> > > http://www.ctcms.nist.gov/fipy >>> > > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy >>> ] >>> > >>> > >>> > _______________________________________________ >>> > fipy mailing list >>> > [email protected] >>> > http://www.ctcms.nist.gov/fipy >>> > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] >>> > >>> > <mesh1D.diff>_______________________________________________ >>> > fipy mailing list >>> > [email protected] >>> > http://www.ctcms.nist.gov/fipy >>> > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] >>> >>> >>> _______________________________________________ >>> fipy mailing list >>> [email protected] >>> http://www.ctcms.nist.gov/fipy >>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] >>> >> >> > > _______________________________________________ > fipy mailing list > [email protected] > http://www.ctcms.nist.gov/fipy > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ] > > -- Daniel Wheeler
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
