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 ]
>>
>
>
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
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
