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 ]
