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 ]
