Re: Mailing list archive

[Raymond Smith]

Oh, thanks! You may want to update the pointer on the FiPy website
<http://www.ctcms.nist.gov/fipy/documentation/MAIL.html> when you get a
chance to avoid having to deal with more of these sorts of requests.

On Wed, Nov 23, 2016 at 7:03 PM, Daniel Wheeler <daniel.wheel...@gmail.com>
wrote:

> Hi Ray,
>
> Gmain is dead unfortunately. Try this instead,
> https://www.mail-archive.com/fipy@nist.gov/.
>
> Cheers,
>
> Daniel
>
> On Wed, Nov 23, 2016 at 9:54 AM, Raymond Smith <smit...@mit.edu> wrote:
> > Hi, FiPy.
> >
> > I'm having trouble accessing the mailing list archive,
> > http://dir.gmane.org/gmane.comp.python.fipy
> > I can get to that page but I can't search or browse through the archive.
> I'm
> > wondering if it's just me or if it's down.
> >
> > Ray
> >
> > ___
> > fipy mailing list
> > fipy@nist.gov
> > http://www.ctcms.nist.gov/fipy
> >   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
> >
>
>
>
> --
> Daniel Wheeler
> ___
> 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 ]

Mailing list archive

[Raymond Smith]

Hi, FiPy.

I'm having trouble accessing the mailing list archive,
http://dir.gmane.org/gmane.comp.python.fipy
I can get to that page but I can't search or browse through the archive.
I'm wondering if it's just me or if it's down.

Ray
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: problem with my initial conditions

[Raymond Smith]

Hi, Damien.

Try using u.setValue(0.5 + 0.5*sin(pi*x)*sin(pi*y)) and similarly for v.
These lines could naturally be placed right after the definition of u and v
(certainly before the time stepping begins).

See
http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.mesh1D.html

Ray

On Wed, Oct 5, 2016 at 5:32 AM, Damien CHUDEAU 
wrote:

> Hello everybody,
>
> I am a French student and I am trying to model the emergence of patterns
> on the fur of animals using the Turing equations (involving the
> concentrations u and v) and to code it thanks to the FiPy library. A
> pattern appearing only if the system is unstable, it is necessary that
> the initial conditions on u and v can be changed. Here is my program
> where everything is initialized to 0.5 ... which does not suit me so. I
> have tried to find help on google and on French forums but no one could
> help me. Can anyone explain me how to solve my problem ?
>
> from fipy import *
> from math import *
>
> a = 0.25
> b = 0.75
> d = 20
> nx = 20
> ny = nx
> dx = 1.
> dy = dx
> L = dx * nx
>
> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
>
> u = CellVariable(name = "concentration u", mesh = mesh, hasOld=True,
> value = 0.5)
> v = CellVariable(name = "concentration v", mesh = mesh, hasOld=True,
> value = 0.5)
>
> D = 1.
>
> eqn0 = TransientTerm(var=u) == a - u + u*u*v + DiffusionTerm(1, var=u)
> eqn1 = TransientTerm(var=v) == b - u*u*v + DiffusionTerm(D, var=v)
>
> eqn = eqn0 & eqn1
>
> #pas d'échange avec l'extérieur -> conditions de Neumann
>
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesTop)
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesBottom)
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesRight)
> u.faceGrad.constrain(((0,),(0,)), where=mesh.facesLeft)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesTop)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesBottom)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesRight)
> v.faceGrad.constrain(((0,),(0,)), where=mesh.facesLeft)
>
> vi = Viewer((u, v))
>
> for t in range(1):
> u.updateOld()
> v.updateOld()
> eqn.solve(dt=1.e-3)
> vi.plot()
>
> I tried adding the following lines of code ... but it does not work.
>
> x, y = mesh.cellCenters
> u.value = 0.5 + 0.5*sin(pi*x)*sin(pi*y)
> v.value = 1 + sin(pi*x)*sin(pi*y)
>
> Thank you in advance
>
> Damien
> ___
> 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 ]

Re: FiPy sweep convergence bottoms out

[Raymond Smith]

I am confused as well now. Comparing with the plots on this Wiki
<https://en.wikipedia.org/wiki/Rate_of_convergence> page which are also
semi-log, it looks to me like Krishna is seeing linear convergence.

On Tue, Sep 27, 2016 at 3:57 PM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> I don't understand what you mean by "supra-linear trend in the semiology
> plot". You show clear 2nd order convergence, which is what I would expect.
>
> > On Sep 27, 2016, at 4:37 PM, Krishna <krishnaku...@imperial.ac.uk>
> wrote:
> >
> > As you can see,  we need supra-linear trend in the semiology plot, such
> as to continue with the  linear drops achieved in the first few sweeps.
> >
> > I.e. the solver is effectively bottoming out. Under-relaxation factors,
> or solver-changes don't seem to work.
> >
> > In fact, for certain under-relaxation factors (including 1.0 for 2 of
> the variables), it breaks the simulation, and produces NaNs right from the
> first sweep.
> >
> > Krishna
> >
> >
> >
> >  Original Message 
> > From: "Gopalakrishnan, Krishnakumar" <krishnaku...@imperial.ac.uk>
> > Sent: Tuesday, September 27, 2016 09:02 PM
> > To: fipy@nist.gov
> > Subject: RE: FiPy sweep convergence bottoms out
> >
> > Thank you Ray, Thanks for pointing that out.
> >
> >
> >
> > Here’s the link to Semilog plot. It takes nearly 22 sweeps to achieve a
> tolerance of 10^-4 for \phi_e and \phi_s_neg.
> >
> >
> >
> > Furthermore, the time spent in sweeping (within each time-step)
> increases as time progresses.
> >
> >
> >
> > https://imperialcollegelondon.box.com/s/4ix6pozs1h9syt1r3fbkw2pi05ooicmy
> >
> >
> >
> > Krishna
> >
> >
> >
> > From: fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] On Behalf Of
> Raymond Smith
> > Sent: Tuesday, September 27, 2016 8:51 PM
> > To: fipy@nist.gov
> > Subject: Re: FiPy sweep convergence bottoms out
> >
> >
> >
> > Hi, Krishna.
> >
> > It would be more clear to plot the residuals on a semi-log plot (or
> equivalently plot the log of residual vs sweep number) to more clearly show
> the value of the small residuals, as the plots in that link make it look to
> me like the residuals all go to zero.
> >
> > Ray
> >
> >
> >
> > On Tue, Sep 27, 2016 at 12:42 PM, Gopalakrishnan, Krishnakumar <
> krishnaku...@imperial.ac.uk> wrote:
> >
> >
> >
> > Hi,
> >
> >
> >
> > We are solving a system of 5 coupled non-linear PDEs.
> >
> >
> >
> > As shown in this plot of residuals vs. sweep count
> https://imperialcollegelondon.box.com/s/9davbq2gq5eani98xuuj2cw9tmz4mbu3
> ,  our residuals die down very slowly, i.e. the solver bottoms out. The
> drop in all the residuals is linear at first, and then asymptotically
> bottoms out to a value.
> >
> >
> >
> > How do we get our residuals to drop faster, i.e. with lesser sweeps and
> faster convergence ? I tried changing solvers and tolerances, but curiously
> enough the results remain identical.
> >
> >
> >
> > Any pointers on this will be much appreciated.
> >
> >
> >
> >
> >
> > Krishna
> >
> >
> >
> >
> >
> >
> >
> >
> > ___
> > 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 ]

Re: FiPy sweep convergence bottoms out

[Raymond Smith]

Hi, Krishna.

It would be more clear to plot the residuals on a semi-log plot (or
equivalently plot the log of residual vs sweep number) to more clearly show
the value of the small residuals, as the plots in that link make it look to
me like the residuals all go to zero.

Ray

On Tue, Sep 27, 2016 at 12:42 PM, Gopalakrishnan, Krishnakumar <
krishnaku...@imperial.ac.uk> wrote:

> Hi,
>
>
>
> We are solving a system of 5 coupled non-linear PDEs.
>
>
>
> As shown in this plot of residuals vs. sweep count
> https://imperialcollegelondon.box.com/s/9davbq2gq5eani98xuuj2cw9tmz4mbu3
> ,  our residuals die down very slowly, i.e. the solver bottoms out. The
> drop in all the residuals is linear at first, and then asymptotically
> bottoms out to a value.
>
>
>
> How do we get our residuals to drop faster, i.e. with lesser sweeps and
> faster convergence ? I tried changing solvers and tolerances, but curiously
> enough the results remain identical.
>
>
>
> Any pointers on this will be much appreciated.
>
>
>
>
>
> Krishna
>
>
>
>
>
>
>
> ___
> 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 ]

Re: Applying Two Different Diffusion Coefficients at Single FV Face

[Raymond Smith]

In general, if you're going to scale things in your problem, I tend to
recommend fully non-dimensionalizing a problem with reference lengths
scales, time scales, etc. and not using different references for different
regions described by the same PDE. In other words, if you have some sort of
relevant time scale, t_ref and a chosen length scale L_ref (which could be,
e.g., L1, L2, L1+L2), I would scale everything with those values, such that
D* = D/D_ref
D_ref = L_ref^2/t_ref
and so on for other quantities in the PDE's including differential
operators and field variables.
In that case, the harmonic mean of D1* and D2* is the same as the harmonic
mean of D1 and D2 divided by D_ref. If you change the references with
position along the length of your PDE, you run into issues with things like
this.

Ray

On Fri, Sep 16, 2016 at 11:36 AM, Gopalakrishnan, Krishnakumar <
krishnaku...@imperial.ac.uk> wrote:

> Hi Raymond,
>
>
>
> Thanks a lot for your inputs. We have been working along the same lines
> that you have suggested.
>
>
>
> There is a slightly more complication that wasn’t fully discussed
> (intentionally, to keep it simple, and to find out if any other options
> existed).
>
>
>
> Assuming that we are working on a normalised co-ordinate and the two bulk
> regions within the domain have non-uniform properties.  In this situation,
> we have to normalise the PDEs too.
>
>
>
> Thus, the diffusion coefficient becomes D_region/L_region  , where region
> can be either region 1 or region 2. , and L_region denotes the original
> length of the region.
>
>
>
> Now, although it’s a perfectly straightforward idea to define D in the
> cell centers and let fiPy handle the harmonic facevalue computation,  I
> don’t know if it’s meaningful to compute harmonicfacevalue of a length
> (which appears in the denominator).
>
>
>
> What do you suggest for L_region ?  Define that at the cell-centers again,
> and use arithmeticfacevalue , or stick with harmonicfacevalue ?
>
>
>
> Best Regards
>
>
>
> Krishna & Ian.
>
>
>
>
>
>
>
> *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf
> Of *Raymond Smith
> *Sent:* Friday, September 16, 2016 7:18 PM
> *To:* fipy@nist.gov
> *Subject:* Re: Applying Two Different Diffusion Coefficients at Single FV
> Face
>
>
>
> A side note, here it may actually be more convenient to think about the D
> in terms of the volumes, so you could also define it as a cell variable
> with values specified as D1 when x<=1 and D2 when x>1, then use the
> harmonicFaceValue attribute when you put it into the governing equation to
> have FiPy do this calculation for you and avoid awkward use of the mesh
> spacing in the specification of the variable values.
>
>
>
> On Fri, Sep 16, 2016 at 11:12 AM, Raymond Smith <smit...@mit.edu> wrote:
>
> Hi, Ian.
>
> I don't think there is such a thing as having two different flux
> coefficients at the same face for the same governing PDE. The flux through
> a given face is calculated by the coefficient at that face times some
> approximation of the gradient in a field variable at that face, like
> D * grad(c),
> both evaluated at the face, so there is only one coefficient in that
> expression. And in the FiPy FV approach D doesn't come into play at all
> except at the faces, so it can't change immediately to the right or left of
> a face but only from one face to the next.
>
> That said, perhaps you could try following a common approach used when a
> FV interface is directly between two bulk regions with different continuum
> properties -- use some sort of a mean of the continuum transport
> coefficient in the adjacent volumes as the approximation for the
> coefficient at the face. I would suggest the harmonic mean (one simple
> reason: it preserves zero flux through the face if either of the adjacent
> cells has a zero coefficient). You could try:
> Dint = 2*D1*D2/(D1+D2)
> Coefficient.setValue(Dint, where=((1.0-dx/2 < 1D_mesh.faceCenters[0]) &
> (1D_mesh.faceCenters[0] < 1.0+dx/2)))
> If you do that, you should see a change in the slope of the field variable
> right near the face between the two bulk regions, which seems to be what
> you're missing. Of course, ideally your mesh is fine enough that the dx is
> insignificant in terms of the solution, but this is a common approach to
> get a slightly more reasonable estimate for what happens between two
> regions.
>
>
>
> Best,
>
> Ray
>
>
>
> On Fri, Sep 16, 2016 at 10:44 AM, Campbell, Ian <
> i.campbel...@imperial.ac.uk> wrote:
>
> Hello All,
>
>
>
> We have a uniformly spaced 1D mesh of length 2.0, wherein a standard
> el

Re: Applying Two Different Diffusion Coefficients at Single FV Face

[Raymond Smith]

A side note, here it may actually be more convenient to think about the D
in terms of the volumes, so you could also define it as a cell variable
with values specified as D1 when x<=1 and D2 when x>1, then use the
harmonicFaceValue attribute when you put it into the governing equation to
have FiPy do this calculation for you and avoid awkward use of the mesh
spacing in the specification of the variable values.

On Fri, Sep 16, 2016 at 11:12 AM, Raymond Smith <smit...@mit.edu> wrote:

> Hi, Ian.
>
> I don't think there is such a thing as having two different flux
> coefficients at the same face for the same governing PDE. The flux through
> a given face is calculated by the coefficient at that face times some
> approximation of the gradient in a field variable at that face, like
> D * grad(c),
> both evaluated at the face, so there is only one coefficient in that
> expression. And in the FiPy FV approach D doesn't come into play at all
> except at the faces, so it can't change immediately to the right or left of
> a face but only from one face to the next.
>
> That said, perhaps you could try following a common approach used when a
> FV interface is directly between two bulk regions with different continuum
> properties -- use some sort of a mean of the continuum transport
> coefficient in the adjacent volumes as the approximation for the
> coefficient at the face. I would suggest the harmonic mean (one simple
> reason: it preserves zero flux through the face if either of the adjacent
> cells has a zero coefficient). You could try:
> Dint = 2*D1*D2/(D1+D2)
> Coefficient.setValue(Dint, where=((1.0-dx/2 < 1D_mesh.faceCenters[0]) &
> (1D_mesh.faceCenters[0] < 1.0+dx/2)))
> If you do that, you should see a change in the slope of the field variable
> right near the face between the two bulk regions, which seems to be what
> you're missing. Of course, ideally your mesh is fine enough that the dx is
> insignificant in terms of the solution, but this is a common approach to
> get a slightly more reasonable estimate for what happens between two
> regions.
>
> Best,
> Ray
>
> On Fri, Sep 16, 2016 at 10:44 AM, Campbell, Ian <
> i.campbel...@imperial.ac.uk> wrote:
>
>> Hello All,
>>
>>
>>
>> We have a uniformly spaced 1D mesh of length 2.0, wherein a standard
>> elliptic diffusion equation (\nabla. (D \nabla \phi) = f is being solved.
>>
>>
>>
>> The diffusion coefficient, D, changes sharply at x = 1.0. Although this
>> is quite similar to examples.diffusion.mesh1D,  there is a key difference:
>>
>>
>>
>> We are seeking an instantaneous change in the coefficient value at that
>> face, i.e. something of the effect that on one side of the face (i.e.
>> from x = 0 up to and including x = 1.0) diffusion coefficient D has a value
>> D1.  Then there is an abrupt change in diffusion coefficient from x >= 1.0
>> (very notable, including the face at x = 1.0) and for the remainder of the
>> mesh, diffusion coefficient D has a value D2.  Currently, we have to either
>> decide on whether the single face at x = 1.0 gets the value D1 or
>> D2. However, when we run with either of these cases, the results are not
>> accurate (when compared against a benchmark from a commercial FV solver).
>> We are getting a smooth interpolation across the face , rather than the
>> sharp change correctly predicted by the benchmarking code.
>>
>>
>>
>> With the code below, we have only been able to set the coefficient at
>> that middle face to one of either D1 or D2. This doesn’t result in the
>> desired instantaneous change in the coefficient value, but instead, the new
>> coefficient value only applies from the NEXT face (which is a distance dx
>> away) from this critical interior boundary.
>>
>>
>>
>> Coefficient = FaceVariable(mesh = 1D_mesh)
>>
>>  #
>> Define facevariable on mesh
>>
>> Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] <= 1.0))
>>
>> # Set coefficient values in 1st half of mesh
>>
>> Coefficient.setValue(D2, where=((1.0 < 1D_mesh.faceCenters[0]) &
>> (1D_mesh.faceCenters[0] <= 2.0)))# Set coefficient
>> values in 2nd half of mesh
>>
>>
>>
>> An alternative with the inequalities adjusted, but that still only
>> permits a single coefficient at that face:
>>
>> Coefficient = FaceVariable(mesh = 1D_mesh)
>>
>>  # Define facevariable
>> on mesh
>>
>> Coefficient.setValue(D1, wh

Re: Applying Two Different Diffusion Coefficients at Single FV Face

[Raymond Smith]

Hi, Ian.

I don't think there is such a thing as having two different flux
coefficients at the same face for the same governing PDE. The flux through
a given face is calculated by the coefficient at that face times some
approximation of the gradient in a field variable at that face, like
D * grad(c),
both evaluated at the face, so there is only one coefficient in that
expression. And in the FiPy FV approach D doesn't come into play at all
except at the faces, so it can't change immediately to the right or left of
a face but only from one face to the next.

That said, perhaps you could try following a common approach used when a FV
interface is directly between two bulk regions with different continuum
properties -- use some sort of a mean of the continuum transport
coefficient in the adjacent volumes as the approximation for the
coefficient at the face. I would suggest the harmonic mean (one simple
reason: it preserves zero flux through the face if either of the adjacent
cells has a zero coefficient). You could try:
Dint = 2*D1*D2/(D1+D2)
Coefficient.setValue(Dint, where=((1.0-dx/2 < 1D_mesh.faceCenters[0]) &
(1D_mesh.faceCenters[0] < 1.0+dx/2)))
If you do that, you should see a change in the slope of the field variable
right near the face between the two bulk regions, which seems to be what
you're missing. Of course, ideally your mesh is fine enough that the dx is
insignificant in terms of the solution, but this is a common approach to
get a slightly more reasonable estimate for what happens between two
regions.

Best,
Ray

On Fri, Sep 16, 2016 at 10:44 AM, Campbell, Ian  wrote:

> Hello All,
>
>
>
> We have a uniformly spaced 1D mesh of length 2.0, wherein a standard
> elliptic diffusion equation (\nabla. (D \nabla \phi) = f is being solved.
>
>
>
> The diffusion coefficient, D, changes sharply at x = 1.0. Although this is
> quite similar to examples.diffusion.mesh1D,  there is a key difference:
>
>
>
> We are seeking an instantaneous change in the coefficient value at that
> face, i.e. something of the effect that on one side of the face (i.e.
> from x = 0 up to and including x = 1.0) diffusion coefficient D has a value
> D1.  Then there is an abrupt change in diffusion coefficient from x >= 1.0
> (very notable, including the face at x = 1.0) and for the remainder of the
> mesh, diffusion coefficient D has a value D2.  Currently, we have to either
> decide on whether the single face at x = 1.0 gets the value D1 or
> D2. However, when we run with either of these cases, the results are not
> accurate (when compared against a benchmark from a commercial FV solver).
> We are getting a smooth interpolation across the face , rather than the
> sharp change correctly predicted by the benchmarking code.
>
>
>
> With the code below, we have only been able to set the coefficient at that
> middle face to one of either D1 or D2. This doesn’t result in the desired
> instantaneous change in the coefficient value, but instead, the new
> coefficient value only applies from the NEXT face (which is a distance dx
> away) from this critical interior boundary.
>
>
>
> Coefficient = FaceVariable(mesh = 1D_mesh)
>
>  #
> Define facevariable on mesh
>
> Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] <= 1.0))
>
> # Set coefficient values in 1st half of mesh
>
> Coefficient.setValue(D2, where=((1.0 < 1D_mesh.faceCenters[0]) &
> (1D_mesh.faceCenters[0] <= 2.0)))# Set coefficient
> values in 2nd half of mesh
>
>
>
> An alternative with the inequalities adjusted, but that still only permits
> a single coefficient at that face:
>
> Coefficient = FaceVariable(mesh = 1D_mesh)
>
>  # Define facevariable on
> mesh
>
> Coefficient.setValue(D1, where=(1D_mesh.faceCenters[0] < 1.0))
>
> # Set coefficient values in 1st half of mesh
>
> Coefficient.setValue(D2, where=((1.0 <= 1D_mesh.faceCenters[0]) &
> (1D_mesh.faceCenters[0] <= 2.0)))  # Set coefficient
> values in 2nd half of mesh
>
>
>
> Sincerely,
>
>
>
> -  Ian
>
>
>
> P.S. Daniel, thank you very much for your previous reply on sweep speeds!
> It was very helpful.
>
>
>
> ___
> 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 ]

Re: Accuracy of DiffusionTerm for non-uniform mesh

[Raymond Smith]

Thanks.

And no, I'm not sure about the normalization for grid spacing. I very well
could have calculated the error incorrectly. I just reported the root mean
square error of the points and didn't weight by volume or anything like
that. I can change that, but I'm not sure of the correct approach.

On Wed, Jul 20, 2016 at 4:25 PM, Daniel Wheeler <daniel.wheel...@gmail.com>
wrote:

> On Wed, Jul 20, 2016 at 1:30 PM, Raymond Smith <smit...@mit.edu> wrote:
> > Hi, FiPy.
> >
> > I was looking over the diffusion term documentation,
> >
> http://www.ctcms.nist.gov/fipy/documentation/numerical/discret.html#diffusion-term
> > and I was wondering, do we lose second order spatial accuracy as soon as
> we
> > introduce any non-uniform spacing (anywhere) into our mesh? I think the
> > equation right after (3) for the normal component of the flux is only
> second
> > order if the face is half-way between cell centers. If this does lead to
> > loss of second order accuracy, is there a standard way to retain 2nd
> order
> > accuracy for non-uniform meshes?
>
> This is a different issue than the non-orthogonality issue, my mistake
> in the previous reply.
>
> > I was playing around with this question here:
> > https://gist.github.com/raybsmith/e57f6f4739e24ff9c97039ad573a3621
> > with output attached, and I couldn't explain why I got the trends I saw.
> > The goal was to look at convergence -- using various meshes -- of a
> simple
> > diffusion equation with a solution both analytical and non-trivial, so I
> > picked a case in which the transport coefficient varies with position
> such
> > that the solution variable is an arcsinh(x). I used three different
> styles
> > of mesh spacing:
> > * When I use a uniform mesh, I see second order convergence, as I'd
> expect.
> > * When I use a non-uniform mesh with three segments and different dx in
> each
> > segment, I still see 2nd order convergence. In my experience, even
> having a
> > single mesh point with 1st order accuracy can drop the overall accuracy
> of
> > the solution, but I'm not seeing that here.
> > * When I use a mesh with exponentially decreasing dx (dx_i = 0.96^i *
> dx0),
> > I see 0.5-order convergence.
>
> That's strange. Are you sure that all the normalization for grid
> spacing is correct when calculation the norms in that last case?
>
> > I can't really explain either of the non-uniform mesh cases, and was
> curious
> > if anyone here had some insight.
>
> I don't have any immediate insight, but certainly needs to addressed.
>
> --
> Daniel Wheeler
> ___
> 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 ]

Accuracy of DiffusionTerm for non-uniform mesh

[Raymond Smith]

Hi, FiPy.

I was looking over the diffusion term documentation,
http://www.ctcms.nist.gov/fipy/documentation/numerical/discret.html#diffusion-term
and I was wondering, do we lose second order spatial accuracy as soon as we
introduce any non-uniform spacing (anywhere) into our mesh? I think the
equation right after (3) for the normal component of the flux is only
second order if the face is half-way between cell centers. If this does
lead to loss of second order accuracy, is there a standard way to retain
2nd order accuracy for non-uniform meshes?

I was playing around with this question here:
https://gist.github.com/raybsmith/e57f6f4739e24ff9c97039ad573a3621
with output attached, and I couldn't explain why I got the trends I saw.
The goal was to look at convergence -- using various meshes -- of a simple
diffusion equation with a solution both analytical and non-trivial, so I
picked a case in which the transport coefficient varies with position such
that the solution variable is an arcsinh(x). I used three different styles
of mesh spacing:
* When I use a uniform mesh, I see second order convergence, as I'd expect.
* When I use a non-uniform mesh with three segments and different dx in
each segment, I still see 2nd order convergence. In my experience, even
having a single mesh point with 1st order accuracy can drop the overall
accuracy of the solution, but I'm not seeing that here.
* When I use a mesh with exponentially decreasing dx (dx_i = 0.96^i * dx0),
I see 0.5-order convergence.

I can't really explain either of the non-uniform mesh cases, and was
curious if anyone here had some insight.

Thanks,
Ray


fipy_convergence.pdf
Description: Adobe PDF document
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: Possibility of setting NaN constraints for CellVariable ?

[Raymond Smith]

Hi, Krishna.

This strikes me as a strange way to go about implementing something, so it
might be helpful to have a bit more physical details about the system in
order to see if the list has proposals for alternative ways to address the
problem. The worry about 0's propagating into the rest of the system makes
me wonder how you expect the rest of the PDE to interact with the nan's. If
values of one PDE are undefined over a particular region, then perhaps the
mesh they are defined on should actually end, and boundary conditions
should be applied. If only one of the PDE's is undefined over that region,
then it would mean having multiple partially overlapping meshes, and that
leads to implementation details, but I give this just as an example
interpretation of your situation.

Ray

On Wed, Jul 13, 2016 at 9:49 AM, Gopalakrishnan, Krishnakumar <
krishnaku...@imperial.ac.uk> wrote:

> Hello,
>
>
>
> I have a regular 2D Cartesian uniformly spaced mesh.
>
>
>
> Due to some peculiarities of the physical system that I am currently
> working with, the variable of interest in one of my PDEs, (fiPy
> CellVariable), is undefined in certain interior node locations, i.e. we
> need to constrain them to a value of float(‘nan’), when solving this PDE.
>
>
>
> It is not satisfactory to merely constrain them to have a value of 0.0 in
> these interior co-ordinates, because a) It is physically inaccurate since
> the values are not defined here and b) Because we need to couple this PDE
> to other PDEs, the incorrect 0.0 values will propagate to other PDEs,
> thereby polluting/corrupting their solutions too.  And hence, the full
> system of
>
>
>
>
>
> However, the following code snippet was deemed to be illegal when I tried
> it out, since FiPy was checking for cellvariables to be strictly numeric
> (got a value < eps error)
>
>
>
> phi.constrain(float(‘nan’),mesh.y<=1.5)
>
>
>
>
>
>
>
> How else can we enforce NaNs in these interior nodes in FiPy ?
>
>
>
>
>
>
>
> Best Regards,
>
>
>
> Krishna
>
>
>
>
>
> ___
> 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 ]

Re: casting implicit Boundary Conditions in FiPy

[Raymond Smith]

Okay, I see your dilemma more clearly now; thanks for clarifying and yes, I
was missing a few things, so this makes more sense now. Anyway, at this
point, we're stepping beyond of my FiPy-fu, so I'll be curious to learn as
well.

Best,
Ray

On Fri, Jun 10, 2016 at 7:59 AM, Gopalakrishnan, Krishnakumar <
k.gopalakrishna...@imperial.ac.uk> wrote:

> Hi Ray,
>
>
>
> Thanks a lot for your further inputs.
>
>
>
> I realise that perhaps I have not been articulating the problem correctly,
> and there is a miscommunication.
>
>
>
> In your code, you have the k3 term as constant. More importantly,  the
> boundary constraint code is set outside the loop , i.e. before beginning
> time-stepping.
>
>
>
> In this case, a single time-step difference will not make a difference,
> i.e . an explicit use of ‘j’ from previous time-step will not be too bad
> (we’ll have to think a bit more about stability, due to the time-delay
> introduced)
>
>
>
> However, my problem is quite different.  The linearization of my
> (non-linear) RHS of the boundary condition is performed around the
> operating point c*.  The operating point changes at every time-step.  Thus,
> the line of code implementing the BC constraint must be inside the loop.
>   The problem now becomes the following:
>
>
>
> 1.   The code print c.faceValue return a numerical result. That means
> that this expression is evaluated numerically, rather than being left as
> part of an implicit exprerssion.  The only way a numerical solution can be
> obtained is to simply return the existing solution, i.e. c.faceValue
> simply uses solution values from the previous time-step.
>
> 2.   After computing the solution, the c* point also needs to be
> changed to the new operating point, i.e. c* = c.  (You always linearise
> around your latest operating point – principle of small-signal perturbation
> and linearization)
>
> 3.   Thus, in the next iteration and all subsequent time-steps,
> c.faceValue – c*  term evaluates to zero again. Thus, we are solving the
> PDE only with the DC component of the Taylor series expansion (i.e.
> effectively ignoring the first-order derivatives too !).
>
>
>
> Note that this sharply contrasts with the problem with the implicit PDE
> set-up, i.e. one could have an implicit source term, which is a linear
> function of solution variable, and fiPy has a built-in method to handle
> this in terms of ImplicitSourceTerm(coeff=…).  You simply obtain the
> linearised mathematical expression of the source term in paper, and feed it
> to the coeff term in the ImplicitSourceTerm method.   c.faceValue  must not
> get evaluated numerically, and must be used implicitly as part of obtaining
> the solution at this time-step.
>
>
>
> I guess describing the problem in words are falling a bit short, since the
> idea is complex enough.  I have drawn a simple flow-chart of the workflow
> for solving this Implicit BC problem.  Can you perhaps kindly take a look
> at it ?
> https://onedrive.live.com/redir?resid=618ADC32B6C2B152!19186=!AHrlNkAcvENEZ0M=3=photo%2cjpg
>
>
>
> Meanwhile, may I also ask if other Fipy users or developers had to deal
> with non-linear Implicit Neumann boundary conditions  in their problems?
>
>
>
>
>
> Krishna
>
>
>
> *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf
> Of *Raymond Smith
> *Sent:* 10 June 2016 00:16
>
> *To:* fipy@nist.gov
> *Subject:* Re: casting implicit Boundary Conditions in FiPy
>
>
>
> Hi, Krisnha.
>
> I don't know the details about how FiPy handles the boundary conditions,
> so I can't confirm whether or not it's solving it with the current
> (implicit) or previous (explicit) time step's value when constraining the
> gradient. I think that will have to come from one of the devs or someone
> more involved than I. However, I did note that in the simple example code
> below, I put it in the same form as you have with the extra terms,
> n*grad(c) = k1 + k2*(c - k3)
> where k1, k2, and k3 are all constants.
> It also seems to work fine. It may be treated explicitly, but it does at
> least seem to update in time correctly. Is there a significant problem with
> it being explicit (using previous time step value) anyway? The overall
> order of accuracy should be the same as implicit (current time step),
> although I guess stability could be worse?
>
> In other words, for the first time step, the (c - cstar) term _should_ be
> very close to zero (exactly zero when treated explicitly), because c is
> initially equal to cstar at the surface. If your time step is large enough
> that c at the surface deviates "significantly" from cstar within a single
&g

Re: casting implicit Boundary Conditions in FiPy

[Raymond Smith]

Hi, Krishna.

Perhaps I'm misunderstanding something, but I'm still not convinced the
second version you suggested -- c.faceGrad.constrain([-(j_at_c_star +
partial_j_at_op_point*(c.faceValue - c_star))], mesh.facesTop) -- isn't
working like you want. Could you look at the example I suggested to see if
that behaves differently than you expect?

Here's the code I used. To me it looks very similar in form to c constraint
above and at first glance it seems to behave exactly like we want -- that
is, throughout the time stepping, n*grad(phi) is constrained to the value,
-phi at the surface. Correct me if I'm wrong, but my impression is that
this is the behavior you desire.

from fipy import *

nx = 50
dx = 1.
mesh = Grid1D(nx=nx, dx=dx)

phi = CellVariable(name="field variable", mesh=mesh, value=1.0)
D = 1.

valueLeft  = 0.0
phi.constrain(valueLeft, mesh.facesLeft)
#phi.faceGrad.constrain([phi], mesh.facesRight) # This is the problematic BC
#phi.faceGrad.constrain(phi.harmonicFaceValue, mesh.facesRight) # This is
the problematic BC
phi.faceGrad.constrain([-phi.harmonicFaceValue], mesh.facesRight) # This is
the problematic BC

eq = TransientTerm() == DiffusionTerm(coeff=D)

timeStep = 0.9 * dx**2 / (2 * D)
steps = 800

viewer = Viewer(vars=phi, datamax=1., datamin=0.)

for step in range(steps):
eq.solve(var=phi, dt=timeStep)
viewer.plot()

Cheers,
Ray

On Thu, Jun 9, 2016 at 11:44 AM, Gopalakrishnan, Krishnakumar <
k.gopalakrishna...@imperial.ac.uk> wrote:

> Hi Ray,
>
>
>
> Yes. You make a good point.  I see that the analytical solution to the
> particular problem I have posted is also zero.
>
>
>
> The reason I posted is because I wanted to present an (oversimplified)
> analogous problem when posting to the group, retaining the generality,
> since many other subject experts might have faced similar situations.
>
>
>
> The actual problem I am solving is the solid diffusion PDE (only 1
> equation) in a Li-ion battery.   I am solving this PDE in a pseudo-2D
> domain.  i.e. I have defined a Cartesian 2D space, wherein the y-coordinate
> corresponds to the radial direction. The bottom face corresponds to
> particle centres, and the top face corresponding to surface of each
> spherical particle.   The x-axis co-ordinate corresponds to particles along
> the width (or thickness) of the positive electrode domain.  Diffusion of Li
> is restricted to be within the solid particle (i.e. y-direction only), by
> defining a suitable tensor diffusion coefficient as described in the
> Anisotropic diffusion example and FAQ in FiPy.  I have normalised my x and
> y dimensions to have a length of unity.
>
>
>
> Now, the boundary condition along the top face is
>
>
>
> Now, j is non-linear (Butler-Volmer), and I am using a Taylor-expanded
> linear version for this boundary condition. All other field variables are
> assumed as constants.   The idea is to set up the infrastructure and solve
> this problem independently, before worrying upon the rubrics of setting up
> the coupled system.  In a similar fashion, I have built up and solved the
> solid phase potential PDE (thanks to your help for pointing out about the
> implicit source term). Thus, the idea is to build up the coupled P2D Newman
> model piecemeal.
>
>
>
> The linearised version of my BC’s RHS at a given operating point ()is
>
>
>
> As you can see, the linearised Boundary condition, is cast in terms of the
> field variable, .  Hence, we need it in an implicit form corresponding
> to  (pseudocode: c.faceGrad.constrain([-( (j_at_c_star -
> partial_j_at_op_point*c_star) + coeff =
> partial_j_at_op_point)],mesh.facesTop) , or something of this
> form/meaning.   (just like the very useful ImplicitSourceterm method)
>
>
>
> If I instead apply the c.faceValue method,i.e. using it in setting the BC
> as
>
>
>
> c.faceGrad.constrain([-(j_at_c_star + partial_j_at_op_point*(c.faceValue -
> c_star))], mesh.facesTop),  then c.faceValue gets immediately evaluated
> at the operating point, c_star,  and we are left with 0 multiplying the
> first-order derivative.
>
>
>
> ie.  the Boundary conditions becomes,
>
>
>
> Leading to huge loss of accuracy.
>
>
>
> Is there any hope at all in this situation ? J . Cheers and thanks for
> your help thus far.
>
>
>
>
>
> Krishna
>
>
>
> *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf
> Of *Raymond Smith
> *Sent:* 09 June 2016 16:06
>
> *To:* fipy@nist.gov
> *Subject:* Re: casting implicit Boundary Conditions in FiPy
>
>
>
> Hi, Krishna.
>
> Could you give a bit more detail and/or an example about how you know it's
> doing the wrong thing throughout the solution process? In the example you
> sent, the

Re: Pure Neumann Boundary Conditions; Elliptic PDE

[Raymond Smith]

Ah, you make a good point: I was using terms incorrectly. I don't actually
know if anyone has studied that PDE with the non-linear Butler-Volmer and
show existence/uniqueness. However, it is true that without phi dependence
on the RHS, it admits infinite solutions. This is demonstrated by noting
that phi only appears as grad(phi) within the equation. Thus, the system
could only ever be determined to an additive constant (the electric _field_
is the actually relevant quantity). However, because the BV equation is
written in terms of phi directly, then there is explicit dependence not on
the field but on the actual value of the electric potential. This is done
by introducing a reference state somewhere else within the equations such
that the absolute value of phi can be given some meaning. To be honest,
I've never done simulations of the solid phase without coupling it to the
electrolyte phase, so I can't comment much more about that. The reaction
term should of course also depend on various other field variables in the
system such as the electric potential in the electrolyte, but I don't
immediately see that assuming that's a uniform constant would cause
troubles.

In terms of implementation, I'm not sure exactly what you're asking.
Specifying a flux can be done indirectly by specifying a gradient
phi.faceGrad.constrain([value], mesh.facesRight)

Ray

On Thu, May 26, 2016 at 3:43 PM, Campbell, Ian <i.campbel...@imperial.ac.uk>
wrote:

> Thank you for your descriptive reply Ray.
>
>
>
> Can you please provide an explanation (or links to an explanation) of the
> connection between the implicit nature of the equation (phi_s on both sides
> of the PDE) and the well-posedness / existence and uniqueness of the
> solution?
>
>
>
> For clarity for the rest of the users on the list, the RHS of the PDE in
> question is strongly non-linear (the j term expands to the Butler-Volmer
> equation).
>
>
>
> Since the right hand side of the PDE is non-linear, theoretically there is
> no guarantee of a unique solution. In this context, we would like to better
> understand your concluding statement that the solution gets automatically
> pinned uniquely for this PDE. The implementation strategy – i.e. how to
> represent these boundary conditions in FiPy – still isn’t clear to us – can
> you kindly elaborate?
>
>
>
> With best regards,
>
>
>
> -  Ian & Krishna
>
>
>
> *From:* fipy-boun...@nist.gov [mailto:fipy-boun...@nist.gov] *On Behalf
> Of *Raymond Smith
> *Sent:* 26 May 2016 19:42
> *To:* fipy@nist.gov
>
> *Subject:* Re: Pure Neumann Boundary Conditions; Elliptic PDE
>
>
>
> At least for now, why don't we keep this on the list.
>
> A couple of thoughts. First, I have actually already implemented a
> Newman-like model in Python. I didn't use Fipy, but I did use the finite
> volume method. It's missing a few features compared to dualfoil, but has
> several extra things. If you are interested in using it, please be in touch
> off the list. I plan to publish it, but I haven't quite gotten around to
> that yet. Hopefully quite soon.
>
> To your questions/comments about the PDE, you raise a couple of critical
> differences. First, your source term, j, is a function of phi (which you
> didn't actually state but is true) and second, the specified flux boundary
> condition at the positive current collector changes the well-posed-ness
> entirely. In this case, the integral of the source term has to match the
> outlet flux into the current collector. That allows for a steady state.
> Then, the fact that j depends on phi is what actually "pins" the solution
> rather than allowing any vertical shift.
>
> Best,
>
> Ray
>
>
>
> On Thu, May 26, 2016 at 2:25 PM, Gopalakrishnan, Krishnakumar <
> k.gopalakrishna...@imperial.ac.uk> wrote:
>
> Hi Ray,
>
>
>
> Many thanks for your quick reply.  I am Krishna, Ian’s colleague here at
> Imperial.  Hope you remember me from our earlier conversations about access
> to Martin’s 10.626 materials.
>
>
>
> Ian & I are currently working on implementing the basic pseudo-2D porous
> electrode Newman model using FiPy.
>
>
>
> The issue is that all of the PDE Boundary conditions are cast as Neumann.
> However, when we try to apply your recommendation in FiPy  (i.e. omit the
> BC definitions/revert to default), it still does not converge.
>
>
>
> ·The PDE that we are currently solving is the solid-potential in
> the positive electrode. Sorry, there was an typo in the last email. What we
> have is:
>
> o   at the positive/separator boundary, there is no potential flux. (a
> no-flux BC)
>
> o at the positive current collector, we have a fixed flux  ( )  BC
>
&

Re: IndexError: index is out of bounds for axis 1 with size

[Raymond Smith]

I think Daniel's point is that Pf is not defined yet when it is first
called in the script. You could fix this by moving the line using Pf that
Daniel pointed to after the definition that you refer to. I think the
general point, however, is that it makes it much easier to help if you try
to make sure that the script you send gives the errors that you see.

Best,
Ray


Virus-free.
www.avast.com

<#m_-68387843283325897_DDB4FAA8-2DD7-40BB-A1B8-4E2AA1F9FDF2>

On Fri, May 6, 2016 at 1:00 PM, Mohammad Kazemi  wrote:

> Thanks Daniel for your response. I thought defining variables as
> cellvariable should be enough.
>
> Pf = CellVariable(mesh=mesh, value=Pdi,hasOld=True)
>
>
> Should I define it in some other way?
>
> Thanks
>
> On Fri, May 6, 2016 at 12:10 PM, Daniel Wheeler  > wrote:
>
>> On Fri, May 6, 2016 at 12:09 PM, Daniel Wheeler
>>  wrote:
>> > Hi Mo,
>>
>> Apologies for that.
>>
>> --
>> Daniel Wheeler
>> ___
>> fipy mailing list
>> fipy@nist.gov
>> http://www.ctcms.nist.gov/fipy
>>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>>
>
>
>
> --
> Mohammad Kazemi
> West Virginia University
> Department of Petroleum and Natural Gas Engineering
>
> ___
> 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 ]

Re: Steady-Source solution of PDE with fixed sources and sinks

[Raymond Smith]

I posted a reply focusing on the parts of this thread which are in the
original question (and thus on the SO site). Feel free to edit/suggest
edits.

Ray

On Mon, Apr 4, 2016 at 7:19 PM, Raymond Smith <smit...@mit.edu> wrote:

> Thanks, Jon.
>
> I'll post something to StackOverflow in the next few days.
>
> Also, that makes sense about the steady state ignoring the initial
> condition and finding the zero solution. I'd forgotten seeing that in the
> example.
>
> Ray
>
> On Mon, Apr 4, 2016 at 10:06 AM, Guyer, Jonathan E. Dr. (Fed) <
> jonathan.gu...@nist.gov> wrote:
>
>> Ray -
>>
>> Thanks for your very complete and very correct answers to Dario.
>>
>> If you wouldn't mind transcribing your answer to Dario's question on
>> StackOverflow, I'd be happy to up-vote it.
>>
>>
>> In answer to *your* questions:
>>
>> > On Apr 3, 2016, at 11:09 AM, Raymond Smith <smit...@mit.edu> wrote:
>> >
>> > Actually, I'm not sure how FiPy treats the steady-state initial guess
>> for Laplace's equation with no flux boundary conditions like yours here.
>> The governing equation + BC's without an initial condition admits any
>> uniform profile as a solution.
>>
>>
>> > I'm still unsure about the treatment of the initial phi values in this
>> sense as mentioned above.
>> > However, the direct-to-steady approach merits a few words of caution.
>> First, there isn't really a good numerical way of directly computing steady
>> state solutions for general systems. Often, your best bet is actually to
>> solve the transient equation by time stepping from some initial condition
>> until you reach steady state, as that's actually probably the most robust
>> algorithm for solving for steady state profiles.
>>
>> The situation is as you expect. We talk about this toward the end of
>>
>>
>> http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.mesh1D.html
>>
>> (search for "Fully implicit solutions are not without their pitfalls").
>> Basically, a steady-state diffusive problem will "lose" its initial
>> condition and should instead be solved by relaxation. The timestep can be
>> made very large; there just needs to be a TransientTerm.
>>
>>
>> Note, also, that it's not necessary to constrain the gradient to zero to
>> get no-flux. FiPy is a cell-centered finite volume code, so no-flux is the
>> natural boundary condition if nothing else is specified.
>>
>> - Jon
>> ___
>> 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 ]

Re: Steady-Source solution of PDE with fixed sources and sinks

[Raymond Smith]

It seems there's a bit of confusion here about initial condition vs. source
term.
Some comments in line.
Hope it helps.

Ray

On Sun, Apr 3, 2016 at 8:38 AM, Dario Panada  wrote:

> Hi,
>
> I originally asked this on
> http://stackoverflow.com/questions/36385367/solving-pde-with-sources-in-python
> but am copying it to this list in case somebody has any advice :) I've
> already read through
> https://www.mail-archive.com/fipy@nist.gov/msg02439.html but didn't solve
> my problem.
>
The list archives can be a good source, and I'd probably even more highly
recommend the examples on the fipy webpage in terms of learning per time.

>
> I'm using FiPy to address a problem inspired by Biology.
>
> Essentially, I want to represent a 2D plane where at different points I
> have sources and sinks. Sources emit substrate at a fixed rate (different
> sources can have different rates) and sinks consume substrate at a fixed
> rate (different sinks can have different rates). My code:
>
> import numpy.matlib
> from fipy import CellVariable, Grid2D, Viewer, TransientTerm,
> DiffusionTerm, ImplicitSourceTerm, ExplicitDiffusionTerm
> from fipy.tools import numerix
> from time import *
>
> nx = 10
> ny = nx
> dx = 1.
> dy = dx
> L = dx*nx
> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
>
> arr_grid = numerix.array((
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,),'d')
>
> arr_source = numerix.array((
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0.5,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,),'d')
>
> arr_sink = numerix.array((
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0,
> 0,0,0,0,0,0,0,0,0,0.5,),'d')
>
> source = CellVariable(mesh=mesh, value = arr_source)
> sink = CellVariable(mesh=mesh, value = arr_sink)
> phi = CellVariable(name = "solution variable", mesh = mesh, value =
> arr_grid)
>
As arr_grid is simply zeros, you could write
phi = CellVariable(name="solution variable", mesh=mesh, value=0.)
and avoid arr_grid entirely.

> X,Y = mesh.cellCenters
> phi.setValue(3.0, where=(X < 2.0) & (X > 1.0))
> phi.setValue(-1.0, where=(X < 6.0) & (X > 5.0))
>
You might consider setting up your source/sink CellVariables with this same
sort of "where" logic so that you avoid writing full arrays out explicitly.
What you've done with this line is establish an initial condition for the
field variable (or an initial guess if you're trying to do a steady-state
solve directly). Also, because this is an initial condition -- assuming phi
is representing, e.g., concentration -- then, it's strange that phi is set
to negative. I think this gets back to the issue that this is an initial
condition and not a source/sink term.

> D = 1.
> eq = TransientTerm() == DiffusionTerm(coeff=D)
>
There is no source/sink in this equation, so the steady profile will be a
perfectly uniform system with the average system concentration defined by
the initial conditions. Actually, I'm not sure how FiPy treats the
steady-state initial guess for Laplace's equation with no flux boundary
conditions like yours here. The governing equation + BC's without an
initial condition admits any uniform profile as a solution. Anyway, the
approach you took below does have the source/sink included in the actual
equations, so that seems like the direction to pursue.

>
>
>
> viewer = Viewer(vars=phi, datamin=0., datamax=1.)
>
> steadyState = False
>
> if(steadyState):
> print("SteadyState")
> DiffusionTerm().solve(var=phi)
> viewer.plot()
> sleep(20)
> else:
> print("ByStep")
> timeStepDuration = 10 * 0.9 * dx**2 / (2 * D)
> steps = 500
> for step in range(steps):
> print(step)
> eq.solve(var=phi,
> dt=timeStepDuration)
> if __name__ == '__main__':
> viewer.plot()
>
> This works well, but FiPy treats the sources as "non-renewable" and
> eventually I get an homogeneous concentration throughout the space as would
> be expected. The alternative was to delete:
>
Here, I'm not sure what you mean by "non-renewable." As written above,
there isn't any source/sink at all.

>
> X,Y = mesh.cellCenters
> phi.setValue(3.0, where=(X < 2.0) & (X > 1.0))
> phi.setValue(-1.0, where=(X < 6.0) & (X > 5.0))
>
> And change the equation to:
>
> eq = 

Re: unsubscribe

[Raymond Smith]

Chris, please try following the instructions on the website and putting
unsubscribe in the body and sending the mail to fipy-requ...@nist.gov ...

http://www.ctcms.nist.gov/fipy/documentation/MAIL.html

On Thu, Mar 24, 2016 at 12:15 PM, Christopher Jones 
wrote:

>
>
> Chris Jones
>
> *Project Electrical Engineer*
>
> cjo...@tantaya.com | 603-601-3340
>
>
>
>
> *ANTAYA* *SCIENCE & TECHNOLOGY*
>
> 7A Merrill Industrial Drive
>
> Hampton,  NH  03842
>
> 603-601-7474
>
>
> This message may contain information that is confidential, privileged or
> otherwise protected from disclosure. If you are not the intended recipient,
> you are hereby notified that any use, disclosure, dissemination,
> distribution, or copying of this message, or any attachments, is strictly
> prohibited. If you have received this message in error, please advise the
> sender by reply e-mail, and delete the message and any attachments. Thank
> you.
>
> ___
> 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 ]

Re: simple convection

[Raymond Smith]

First, I hope you're considering the broad set of FiPy examples
 as part of the
documentation. They can serve as a great starting point for many problems.
Second, solving convection problems without any diffusion is notoriously a
bit difficult with standard finite volume techniques. So, despite the
apparent simplicity of your starting example, you've happened to pick one
that's likely to disagree with the analytical solution because of numerical
challenges. That said, you may try using different forms of the convection
term (mentioned here
)
and letting the simulation run for convection beyond a single grid point.
Anyway, it's usually a good idea to start with simple examples like this,
but here, you might consider letting FiPy's examples guide your choice in
getting started unless you're specifically looking at problems of the form
you've begun with here.

If you really need to solve convection problems without the numerically
smoothing diffusion, then you might consider different numerical schemes
designed for hyperbolic equations, such as those implemented in clawpack:
http://www.clawpack.org/
Very similar to your problem:
http://www.clawpack.org/pyclaw/gallery/gallery_all.html#dimensional-advection

Best,
Ray

On Wed, Mar 23, 2016 at 7:54 AM, Tyler ABBOT 
wrote:

> Hello!
>
> I'm new to FiPy and am having a lot of trouble finding my way using only
> the documentation...
>
> I am trying to understand how FiPy works by working an example, in
> particular I would like to solve the following simple convection equation
> with periodic boundary:
>
> $$\partial_t u + \partial_x u = 0$$
>
> If initial data is given by $u(x, 0) = F(x)$, then the analytical solution
> is $u(x, t) = F(x - t)$. I do get a solution, but it is not correct.
>
> What am I missing? Is there a better resource for understanding FiPy than
> the documentation? It is very sparse...
>
> Here is my attempt
>
> from fipy import *
> import numpy as np
>
> # Generate mesh
> nx = 20
> dx = 2*np.pi/nx
> mesh = PeriodicGrid1D(nx=nx, dx=dx)
>
> # Generate solution object with initial discontinuity
> phi = CellVariable(name="solution variable", mesh=mesh)
> phiAnalytical = CellVariable(name="analytical value", mesh=mesh)
> phi.setValue(1.)
> phi.setValue(0., where=x > 1.)
>
> # Define the pde
> D = [[-1.]]
> eq = TransientTerm() == ConvectionTerm(coeff=D)
>
> # Set discretization so analytical solution is exactly one cell translation
> dt = 0.01*dx
> steps = 2*int(dx/dt)
>
> # Set the analytical value at the end of simulation
> phiAnalytical.setValue(np.roll(phi.value, 1))
>
> for step in range(steps):
> eq.solve(var=phi, dt=dt)
>
> print(phi.allclose(phiAnalytical, atol=1e-1))
>
> Thanks for any guidance you could give!
>
> Tyler Abbot
> Phd Student
> Department of Economics
> Sciences Po, Paris
> tyler.ab...@sciencespo.fr
>
> ___
> 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 ]

Re: Ideal Free Energy: Computation of x*log(x) at x -> 0

[Raymond Smith]

Hi, Aniruddha.

I certainly don't know the details of your system, but my first thought
would be that this term in the free energy should (physically, at least)
prevent the concentration from actually ever reaching zero, as the chemical
potential diverges when x->0. So I don't understand why you would
initialize the system with zero concentrations. From my not-super-fresh
memory of deriving the ideal solution model from a lattice model in Thermo,
I didn't think it could really apply until you could at least apply
Stirling's approximation, implying a moderately large number of the species
in the system. Anyway, is there a strong reason that you need the
concentrations to be zero? In my experience with such models, when run
dynamically, they keep concentrations away from zero.

Ray

On Tue, Feb 9, 2016 at 8:11 PM, Aniruddha Jana  wrote:

> Hello Everyone,
>
> I am trying to use an ideal chemical free energy formulation for a
> problem, where I have terms of the form x*log(x). Computing x*log(x) at x =
> 0 gives the error "RuntimeWarning: divide by zero encountered in log".
> Even if I try to prevent computation of x*log(x) at x = 0 by doing (for
> example):
>
> epsilon = 1.0e-3
> G.setValue((x>epsilon)*log(x>epsilon))
>
> or
>
> G.setValue(x*log(x), where=(x>0))
>
> the error persists.
>
> A test script is attached with this email.
>
> Looking up online, I read that "log" is computed before "where". Another
> suggestion says to catch and handle the exception (links below).
>
>
> http://stackoverflow.com/questions/25087769/runtimewarning-divide-by-zero-error-how-to-avoid-python-numpy
>
>
> http://stackoverflow.com/questions/13497891/python-getting-around-division-by-zero
>
> Any thoughts on how to get rid of the problem, especially for FiPy, and/or
> the ideal free energy formulation?
>
> Thank You!
>
> Best regards,
> Aniruddha
>
>
>
> ___
> 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 ]

Re: Ideal Free Energy: Computation of x*log(x) at x -> 0

[Raymond Smith]

Hi, Aniruddha.

Thanks for the clarification. I've also run into issues with concentrations
hitting zero in ideal-solution-like simulations in the past.

One thing worth mentioning, if you're imposing a flux that's too large for
the physical system to handle, (e.g. above a diffusion limited current),
then the system will certainly "explode" at some finite time when the
concentration at the surface reaches zero. The phenomenon of Sand's Time in
electrochemistry is related to this issue.

So, assuming you're not trying to ask the system to do something
unphysical, I've _still_ run into issues in the past often associated with
my poor numerical choices. The first simple point is to make sure your time
steps aren't too large, such that any given cell isn't overstepping in the
direction toward zero. Also, one issue I had at one point was that I was
(in an electroneutral electrolyte system with electromigration flux terms)
using a linear mean for faceValue's instead of a harmonic mean. That can
lead to situations in which flux leaves one side of an "empty" cell with no
inlets.

Anyway, my first thought is that the issue might be avoided with more
careful numerics, but perhaps that is not your issue. I've also run into
this issue for more recalcitrant systems where I couldn't figure out how to
avoid hitting zero with my ideal-solution-like concentrations. What I ended
up doing to "cheat" is not an elegant solution, but it kind-of works. I
linearly extended the log terms beyond zero with a piecewise function. That
is:
f(x) = log(x) if x > epsilon
f(x) = log(epsilon) + 1/epsilon * (x - epsilon)
which should preserve continuity in the function and first derivative but
nothing above that. The result is that concentrations in the system can
actually go negative, but they really don't like being there (as long as
epsilon is small enough), and at least don't go very negative. Again, not a
great solution, but it may be something to try.

Ray

On Wed, Feb 10, 2016 at 2:30 PM, Aniruddha Jana <aj...@purdue.edu> wrote:

> Hello Ray,
>
> Thank you for the reply. I should have explained my problem better, sorry
> for that.
>
> Yes, I do not initialize the concentration to zero, unlike in the test
> script I have shared.  But due to a flux at the boundary in my actual
> problem, the concentration in some regions can decrease to very low values
> close to zero.
>
> Any thoughts on how to address the log(x>epsilon) problem that gives the
> Runtime Warning, that I mentioned in my previous email? Even though it is a
> "warning", the solutions to my equations (not in test script) result in a
> "nan" because of undesired values computed by the log function, so I cannot
> ignore the Runtime Warning.
>
> Thanks,
> Aniruddha
>
> On Wed, Feb 10, 2016 at 9:13 AM, Raymond Smith <smit...@mit.edu> wrote:
>
>> Hi, Aniruddha.
>>
>> I certainly don't know the details of your system, but my first thought
>> would be that this term in the free energy should (physically, at least)
>> prevent the concentration from actually ever reaching zero, as the chemical
>> potential diverges when x->0. So I don't understand why you would
>> initialize the system with zero concentrations. From my not-super-fresh
>> memory of deriving the ideal solution model from a lattice model in Thermo,
>> I didn't think it could really apply until you could at least apply
>> Stirling's approximation, implying a moderately large number of the species
>> in the system. Anyway, is there a strong reason that you need the
>> concentrations to be zero? In my experience with such models, when run
>> dynamically, they keep concentrations away from zero.
>>
>> Ray
>>
>> On Tue, Feb 9, 2016 at 8:11 PM, Aniruddha Jana <aj...@purdue.edu> wrote:
>>
>>> Hello Everyone,
>>>
>>> I am trying to use an ideal chemical free energy formulation for a
>>> problem, where I have terms of the form x*log(x). Computing x*log(x) at x =
>>> 0 gives the error "RuntimeWarning: divide by zero encountered in log".
>>> Even if I try to prevent computation of x*log(x) at x = 0 by doing (for
>>> example):
>>>
>>> epsilon = 1.0e-3
>>> G.setValue((x>epsilon)*log(x>epsilon))
>>>
>>> or
>>>
>>> G.setValue(x*log(x), where=(x>0))
>>>
>>> the error persists.
>>>
>>> A test script is attached with this email.
>>>
>>> Looking up online, I read that "log" is computed before "where". Another
>>> suggestion says to catch and handle the exception (links below).
>>>
>>>
>>> http://stackoverflow.com/questions/25087769/run

Re: trouble installing Fipy3.1

[Raymond Smith]

Looks to me like pip thinks you already have fipy installed (Requirement
already satisfied). I guess you've already tried running the test and it
failed? What was the output of that?

Ray

On Sat, Feb 28, 2015 at 2:25 PM, Franck Kalala franckkal...@googlemail.com
wrote:

 I have install ez_setup package, now I get this

 divine@franck-laptop:~$ sudo pip install fipy
 Requirement already satisfied (use --upgrade to upgrade): fipy in
 /usr/local/lib/python2.7/dist-packages/FiPy-3.1-py2.7.egg
 Cleaning up...


 cheers

 2015-02-28 0:26 GMT+02:00 Raymond Smith smit...@mit.edu:

 It means you don't have a python module called ez_setup. See
 https://pypi.python.org/pypi/ez_setup. Maybe try
 $ sudo pip install ez_setup

 Ray

 On Fri, Feb 27, 2015 at 5:23 PM, Franck Kalala 
 franckkal...@googlemail.com wrote:

 I am trying to install Fipy, but I am getting this

 divine@franck-laptop:~/Downloads/FiPy-3.1$ sudo python setup.py install
 [sudo] password for divine:
 Traceback (most recent call last):
   File setup.py, line 44, in module
 import ez_setup
 ImportError: No module named ez_setup

 Any help?

 --
 Franck Kalala Mutombo, Ph.D.
 https://sites.google.com/a/aims.ac.za/dr-franck-kalalaf-mutombo/
 Tel: +243 824018007





 ___
 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 ]




 --
 Franck Kalala Mutombo, Ph.D.
 https://sites.google.com/a/aims.ac.za/dr-franck-kalalaf-mutombo/
 Tel: +243 824018007





 ___
 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 ]

Re: trouble installing Fipy3.1

[Raymond Smith]

You should install packages through the package manager when possible. Try
installing pip through the repositories (and you should first try to
uninstall what you did manually). Then try to install fipy using pip.

On Sat, Feb 28, 2015 at 2:51 PM, Franck Kalala franckkal...@googlemail.com
wrote:

 UBUNTU 14.04

 no fipy in the repositories.
 intalles pip manually


 2015-02-28 21:48 GMT+02:00 Raymond Smith smit...@mit.edu:

 What Linux distribution are you using? Does it have FiPy in the
 repositories? If so, you should probably install FiPy from there (unless
 you need a newer version that your distribution provides). Did you install
 pip from the distribution's repositories?

 On Sat, Feb 28, 2015 at 2:45 PM, Franck Kalala 
 franckkal...@googlemail.com wrote:

 divine@franck-laptop:~$ python -c from distutils.sysconfig import
 get_python_lib; print(get_python_lib())
 /usr/lib/python2.7/dist-packages
 divine@franck-laptop:~$  ls -l `which python`
 lrwxrwxrwx 1 root root 9 Feb 26 10:05 /usr/bin/python - python2.7
 divine@franck-laptop:~$


 2015-02-28 21:39 GMT+02:00 Raymond Smith smit...@mit.edu:

 It looks like pip is installing to
 /usr/local/lib/python2.7/dist-packages/ (very normal). Perhaps that's not
 where you're python is looking for packages. What is the output of

 $  python -c from distutils.sysconfig import get_python_lib;
 print(get_python_lib())

 (from
 https://stackoverflow.com/questions/122327/how-do-i-find-the-location-of-my-python-site-packages-directory
 )
 Also, is there an issue with python multiple versions and, e.g.
 python pointing to python3? What is output from

 $ ls -l `which python`

 Ray

 On Sat, Feb 28, 2015 at 2:33 PM, Franck Kalala 
 franckkal...@googlemail.com wrote:

 test give this

 divine@franck-laptop:~$ python -c import fipy; fipy.test()
 Traceback (most recent call last):
   File string, line 1, in module
 ImportError: No module named fipy
 divine@franck-laptop:~$


 2015-02-28 21:27 GMT+02:00 Raymond Smith smit...@mit.edu:

 Looks to me like pip thinks you already have fipy installed
 (Requirement already satisfied). I guess you've already tried running 
 the
 test and it failed? What was the output of that?

 Ray


 On Sat, Feb 28, 2015 at 2:25 PM, Franck Kalala 
 franckkal...@googlemail.com wrote:

 I have install ez_setup package, now I get this

 divine@franck-laptop:~$ sudo pip install fipy
 Requirement already satisfied (use --upgrade to upgrade): fipy in
 /usr/local/lib/python2.7/dist-packages/FiPy-3.1-py2.7.egg
 Cleaning up...


 cheers

 2015-02-28 0:26 GMT+02:00 Raymond Smith smit...@mit.edu:

 It means you don't have a python module called ez_setup. See
 https://pypi.python.org/pypi/ez_setup. Maybe try
 $ sudo pip install ez_setup

 Ray

 On Fri, Feb 27, 2015 at 5:23 PM, Franck Kalala 
 franckkal...@googlemail.com wrote:

 I am trying to install Fipy, but I am getting this

 divine@franck-laptop:~/Downloads/FiPy-3.1$ sudo python setup.py
 install
 [sudo] password for divine:
 Traceback (most recent call last):
   File setup.py, line 44, in module
 import ez_setup
 ImportError: No module named ez_setup

 Any help?

 --
 Franck Kalala Mutombo, Ph.D.
 https://sites.google.com/a/aims.ac.za/dr-franck-kalalaf-mutombo/
 Tel: +243 824018007





 ___
 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 ]




 --
 Franck Kalala Mutombo, Ph.D.
 https://sites.google.com/a/aims.ac.za/dr-franck-kalalaf-mutombo/
 Tel: +243 824018007





 ___
 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
 ]




 --
 Franck Kalala Mutombo, Ph.D.
 https://sites.google.com/a/aims.ac.za/dr-franck-kalalaf-mutombo/
 Tel: +243 824018007





 ___
 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 ]




 --
 Franck Kalala Mutombo, Ph.D.
 https://sites.google.com/a/aims.ac.za/dr-franck-kalalaf-mutombo/
 Tel: +243 824018007





 ___
 fipy mailing list
 fipy@nist.gov
 http://www.ctcms.nist.gov/fipy
   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy

Re: trouble installing Fipy3.1

[Raymond Smith]

It looks like pip is installing to /usr/local/lib/python2.7/dist-packages/
(very normal). Perhaps that's not where you're python is looking for
packages. What is the output of

$  python -c from distutils.sysconfig import get_python_lib;
print(get_python_lib())

(from
https://stackoverflow.com/questions/122327/how-do-i-find-the-location-of-my-python-site-packages-directory
)
Also, is there an issue with python multiple versions and, e.g. python
pointing to python3? What is output from

$ ls -l `which python`

Ray

On Sat, Feb 28, 2015 at 2:33 PM, Franck Kalala franckkal...@googlemail.com
wrote:

 test give this

 divine@franck-laptop:~$ python -c import fipy; fipy.test()
 Traceback (most recent call last):
   File string, line 1, in module
 ImportError: No module named fipy
 divine@franck-laptop:~$


 2015-02-28 21:27 GMT+02:00 Raymond Smith smit...@mit.edu:

 Looks to me like pip thinks you already have fipy installed (Requirement
 already satisfied). I guess you've already tried running the test and it
 failed? What was the output of that?

 Ray


 On Sat, Feb 28, 2015 at 2:25 PM, Franck Kalala 
 franckkal...@googlemail.com wrote:

 I have install ez_setup package, now I get this

 divine@franck-laptop:~$ sudo pip install fipy
 Requirement already satisfied (use --upgrade to upgrade): fipy in
 /usr/local/lib/python2.7/dist-packages/FiPy-3.1-py2.7.egg
 Cleaning up...


 cheers

 2015-02-28 0:26 GMT+02:00 Raymond Smith smit...@mit.edu:

 It means you don't have a python module called ez_setup. See
 https://pypi.python.org/pypi/ez_setup. Maybe try
 $ sudo pip install ez_setup

 Ray

 On Fri, Feb 27, 2015 at 5:23 PM, Franck Kalala 
 franckkal...@googlemail.com wrote:

 I am trying to install Fipy, but I am getting this

 divine@franck-laptop:~/Downloads/FiPy-3.1$ sudo python setup.py
 install
 [sudo] password for divine:
 Traceback (most recent call last):
   File setup.py, line 44, in module
 import ez_setup
 ImportError: No module named ez_setup

 Any help?

 --
 Franck Kalala Mutombo, Ph.D.
 https://sites.google.com/a/aims.ac.za/dr-franck-kalalaf-mutombo/
 Tel: +243 824018007





 ___
 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 ]




 --
 Franck Kalala Mutombo, Ph.D.
 https://sites.google.com/a/aims.ac.za/dr-franck-kalalaf-mutombo/
 Tel: +243 824018007





 ___
 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 ]




 --
 Franck Kalala Mutombo, Ph.D.
 https://sites.google.com/a/aims.ac.za/dr-franck-kalalaf-mutombo/
 Tel: +243 824018007





 ___
 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 ]

Re: hints on creating a Source in simple convection diffusion problem

[Raymond Smith]

Have you tried printing scalarSource to make sure it is at least 0.0 at
some locations and 100.0 at other locations? This might be easier to see if
you make a small, dx=6 x dy=2 mesh.
And you do want it to be a source term like
dcdt = diffusion + convection + 100.0 * c
right? As in, a first order reaction, not zero order. Otherwise, I could
see doing a simple test in which you start out with c=0 everywhere and
having nothing being generated. If you want zero order reaction, you should
add the term as
dcdt = diffusion + convection + scalarSource

Ray

On Fri, Jan 23, 2015 at 3:49 PM, Richard Gillilan r...@cornell.edu wrote:

 Doing simple convection diffusion in a narrow channel: simulating the
 accumulation of damaged solute in a dilute solution as a result of x-ray
 exposure. I’m having trouble getting a Source term to work the way I think
 it should. I want a strip of volume elements in the middle of the channel
 (the x-ray illuminated volume) to act as a source of “damaged” solute. So
 any solution that convects into those pixels should have the concentration
 of damaged material increased at a constant rate. Note vel is simple
 Poiseuille flow, but set to zero for testing purposes.

 Can anyone spot any obvious errors? Here’s what I tried:

 I’m only showing some of the code here for brevity

 print volume element size: ,dx,x,dy

 Lx = dx*nx
 Ly = dy*ny

 mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

 print dimensions of capillary: ,Lx,Ly
 ...
 xc,yx = mesh.cellCenters
 ...
 scalarSource = CellVariable(mesh, value=0.0, rank=0)
 strip = ((xcLx/3)  (xc2*Lx/3))
 scalarSource.setValue([100.0], where=strip)

 eq = TransientTerm() ==
 (DiffusionTerm(coeff=D)+VanLeerConvectionTerm(coeff=vel))+ImplicitSourceTerm(scalarSource)

 ___
 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 ]

Re: Sweep causes solution to zero out.

[Raymond Smith]

This seems to be a common issue here... I think the diffusion 1D example on
the website may soon comment on this. The correct equation for conservation
of thermal energy is that the time rate of change of the enthalpy is
related to the divergence of the flux. Enthalpy is related to rho * cp * T,
and flux of enthalpy is related to k * grad(T). Thus, the correct way to do
this is something more like

TransientTerm(coeff=rho * cp) == DiffusionTerm(coeff=k) + Source

Best,
Ray

On Thu, Jan 8, 2015 at 4:40 PM, Christopher Jones cjo...@tantaya.com
wrote:

 Thank you for the insight.

 The coefficient for the diffusion term is

 thermal conductivity / (density * specific heat)

 where thermal conductivity and specific heat are both non linearly
 temperature dependent. Thermal conductivity is approximated by
 10^(a(u)/b(u)) where a and b are polynomials in u.

 Specific heat is approximated by a piece-wise defined cubic spline in
 u/100.

 So, just to be sure, would you recommend using the chain rule as given in
 your first reply? (DiffusionTerm(coeff=A) + A.grad.dot(u.grad))

 Thanks again,
 Chris




 Chris Jones

 *Project Electrical Engineer*

 cjo...@tantaya.com | 603-601-3340




 *ANTAYA* *SCIENCE  TECHNOLOGY*

 7A Merrill Industrial Drive

 Hampton,  NH  03842

 603-601-7474


 This message may contain information that is confidential, privileged or
 otherwise protected from disclosure. If you are not the intended recipient,
 you are hereby notified that any use, disclosure, dissemination,
 distribution, or copying of this message, or any attachments, is strictly
 prohibited. If you have received this message in error, please advise the
 sender by reply e-mail, and delete the message and any attachments. Thank
 you.

 On Thu, Jan 8, 2015 at 2:05 PM, Guyer, Jonathan E. Dr. 
 jonathan.gu...@nist.gov wrote:


 On Jan 8, 2015, at 10:30 AM, Christopher Jones cjo...@tantaya.com
 wrote:

  I am modelling nonlinear heat conduction using time steps and sweeps.
 
  Basic equation is:
 
  du/dt = A(u) * d2u/dx2 + S(x, t)
 
  where u(x, t)
 
  I have functions written for A(u), the DiffusionTerm coefficient, and
 S(x, t), the ImplicitSourceTerm. I implement them with the following
 snippet:
 
  eq = TransientTerm() == DiffusionTerm(coeff=A) +
 ImplicitSourceTerm(coeff=S)


 DiffusionTerm represents d/dx(A(u) * du/dx) (or, more rigorously,
 nabla\cdot(A(u) \nabla u)).

 If your functional form *really* is A(u) * d2u/dx2, then you need to run
 the chain rule to obtain \nabla\cdot(A(u) \nabla u) - \nabla A(u) \cdot
 \nabla u, which equates to DiffusionTerm(coeff=A) + A.grad.dot(u.grad). It
 is rare, however for this to be needed. The form of DiffusionTerm (with A
 inside the first derivative) is the one almost invariably encountered.

  S is really the whole source term, not the coefficient, so have I coded
 this correctly?

 Raymond has already pointed out the linearity of ImplicitSourceTerm. If
 your source is linear in u, then using +ImplicitSourceTerm(coeff=S1) is
 more efficient than using +S1*u. If your source isn't linear in u, then
 don't worry about it and write the source explicitly.

  I have declared A as a FaceVariable (because diffusion happens from
 cell to cell, and S as a CellVariable, is this correct?

 Yes.

  Since S depends on t as well, I created a Variable t. I then created a
 small function and included it with S with this snippet:
 
  B = 50 * numerix.exp(-t)
 
  S.setValue(B)

 
  After sweeping, I increment the elapsed time and then t, with
 
  t.setValue(elapsedTime)
 
  I can print B's value to see that it does decay.

 .setValue() sets the instantaneous value of S. S won't change with
 further changes in t.

 Skip B and just write

 S = 50 * numerix.exp(-t)



 ___
 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 ]

Re: Sweep causes solution to zero out.

[Raymond Smith]

As a comment, if S is the whole source, you don't want to use
ImplicitSourceTerm, because that assumes that the source is linear in the
variable with coefficient as given. This is apparently usually the
recommended way to do it if there is a way to write your source that way.
However, if you have a source that can't be written that way, you can just
type the source term into the equation directly: documentation
http://www.ctcms.nist.gov/fipy/documentation/numerical/discret.html#source-term
(and an example here
http://www.ctcms.nist.gov/fipy/examples/convection/generated/examples.convection.exponential1DSource.mesh1D.html#module-examples.convection.exponential1DSource.mesh1D
)

Ray

On Thu, Jan 8, 2015 at 10:30 AM, Christopher Jones cjo...@tantaya.com
wrote:

 Hello,

 I am modelling nonlinear heat conduction using time steps and sweeps.

 Basic equation is:

 du/dt = A(u) * d2u/dx2 + S(x, t)

 where u(x, t)

 I have functions written for A(u), the DiffusionTerm coefficient, and S(x,
 t), the ImplicitSourceTerm. I implement them with the following snippet:

 eq = TransientTerm() == DiffusionTerm(coeff=A) +
 ImplicitSourceTerm(coeff=S)

 S is really the whole source term, not the coefficient, so have I coded
 this correctly?

 I have declared A as a FaceVariable (because diffusion happens from cell
 to cell, and S as a CellVariable, is this correct?

 Since S depends on t as well, I created a Variable t. I then created a
 small function and included it with S with this snippet:

 B = 50 * numerix.exp(-t)

 S.setValue(B)

 After sweeping, I increment the elapsed time and then t, with

 t.setValue(elapsedTime)

 I can print B's value to see that it does decay.


 I provide an initial distribution: u(x, 0). For simplicity, I will choose
 no boundary conditions. The first set of sweeps causes u to zero out.

 What is happening? Have I not set this up correctly? What examples could I
 use to iron this out?

 Thanks,
 Chris

 ___
 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 ]

Re: Convergence suggestions?

[Raymond Smith]

Hi, Kyle.

I'll go ahead and keep this on the list for posterity (and also to allow
other people to correct me if I say incorrect things).

On Sun, Nov 23, 2014 at 12:29 AM, Kyle Briton Lawlor klawlor...@gmail.com
wrote:

 Hi, Ray.

 I’ve been taking a look at your code to learn some stuff. I have a few
 questions. They are not particularly related to the questions you are
 asking, so I thought I would message you separate from the FiPy list.

 Anyways, I am working on a problem where I need to invoke a boundary
 condition involving a gradient and I saw in this code you had a similar
 thing.

 phi.constrain(0., mesh.facesLeft)
 phi.constrain(V, mesh.facesRight)
 c_m.constrain(c0, mesh.facesLeft)
 c_m.faceGrad.constrain(c_m.faceValue*phi.faceGrad, mesh.facesRight)

 My question is does the value in parentheses here,
 (c_m.faceValue)*phi.faceGrad, change as the solution is computed for c_m? I
 suppose it is fixed here, c_m.faceValue=1. Is this correct?


The value of c_m.faceValue does indeed change as the solution is computed
for c_m. So, no c_m.faceValue is not fixed, nor is phi.faceGrad. Both are
terms which are dependent on the governing equations as well as the other
boundary conditions, so that the computed answer is the solution such that
the normal gradient of c_m at the boundary is equal to the value of c_m at
the boundary times the normal gradient of phi at the boundary. You can
verify this by post-computing n*grad(c_m), c_m,and n*grad(phi) at the
boundary. Then, you can see that n*grad(c_m) = c_m * n*grad(phi), where all
of those terms come directly from the solution variables. (here, n* is the
external surface normal vector dot operation). This may be made a bit more
clear in the way the analytical solution is worked out in the original
attachment, and agreement between analytical and numerical results confirms
that the non-linear, coupled boundary condition is treated properly in FiPy.

This is one of the very nice aspects of FiPy in my opinion. Although
convergence is perpetually an issue with coupled non-linear problems, the
ease with which you can implement even particularly nasty-looking boundary
conditions is extremely convenient.


 Best,
 Kyle


 On Nov 19, 2014, at 5:12 PM, Raymond Smith smit...@mit.edu wrote:

 Hi, FiPy.

 I'm trying to solve a number of problems involving (quasi-)neutral
 electrolyte solutions. I've attached some notes I made a while ago on one
 of the simplest example problems, for which there is an analytical
 solution. The equations are non-linear, and my typical boundary conditions
 are also often non-linear, as in the attached example.

 I'm curious if there are suggestions for improving convergence for this
 type of problem in FiPy. In the simple example (code here
 https://gist.github.com/raybsmith/6e724e6f1ecc98f9303b and below), we
 have the analytical solution for any applied voltage. For moderate positive
 applied voltages, I can get convergence -- up to around 10 thermal volts.
 However, for negative applied voltages (which corresponds to depletion of
 the salt in the electrolyte), I only get convergence down to about -0.5
 thermal volts, which is notably less than I'd like to simulate. I realize
 that sometimes solving non-linear problems numerically is simply difficult,
 but perhaps there are some things I could try that someone here could
 recommend?

 I've thought about doing the transient version of the problem and stepping
 toward a steady state with more extreme applied potentials. However,
 because the electrostatic part is static, one of the equations will always
 be stationary, even in the dynamic case. (avoiding electrodynamics here) So
 in some similar problems, I've still found I can't apply large potentials.

 Thanks,
 Ray

 import fipy as fp

 Lx = 1.
 nx = 1000
 dx = Lx/nx
 mesh = fp.Grid1D(nx=nx, dx=dx)
 phi = fp.CellVariable(name=elec. potential, mesh=mesh, value=0.)
 c_m = fp.CellVariable(name=conc., mesh=mesh, value=1.0)

 #V = -6e-1
 V = -4e-1
 #V = 4.0e0
 c0 = 1.
 c_m.setValue(c0)
 D_m = 2.03e-9 # m^2/s
 D_p = 1.33e-9 # m^2/s
 D_p = D_p/D_m
 D_m = 1.

 eq1 = (0 == fp.DiffusionTerm(coeff=-D_m, var=c_m) +
 fp.DiffusionTerm(coeff=D_m*c_m, var=phi))
 eq2 = (0 == fp.DiffusionTerm(coeff=-(D_p - D_m), var=c_m) +
 fp.DiffusionTerm(coeff=-(D_p + D_m)*c_m, var=phi))

 phi.constrain(0., mesh.facesLeft)
 phi.constrain(V, mesh.facesRight)
 c_m.constrain(c0, mesh.facesLeft)
 c_m.faceGrad.constrain(c_m.faceValue*phi.faceGrad, mesh.facesRight)

 Xcc = mesh.cellCenters[0]
 J = 1 - fp.numerix.exp(V)
 c_m_analyt = fp.CellVariable(name=analyt. conc., mesh=mesh)
 c_m_analyt.setValue(1 - J*Xcc)
 phi_analyt = fp.CellVariable(name=analyt. phi, mesh=mesh)
 phi_analyt.setValue(fp.numerix.log(1-J*Xcc))

 eq = eq1  eq2

 res = 1e10
 nIt = 0
 while res  1e-6:
 res = eq.sweep()
 print Iteration:, nIt, Residual:, res
 nIt += 1

 anstol = 1e-3
 print c_m.allclose(c_m_analyt, rtol=anstol, atol=anstol)
 print phi.allclose(phi_analyt

Re: Convergence suggestions?

[Raymond Smith]

Well, I think I spoke a bit too soon. By doing the transient version of the
problem and by switching to GMRESSolver, I was able to approach steady
state for considerably higher applied voltages (about an order of magnitude
higher), which is considerably better and getting quite close to full salt
depletion on one side. I still have trouble sweeping to convergence on the
time steps for large applied voltages, but perhaps by taking smaller and
smaller time steps I'll be able to continue to make progress. Tips for
solving the steady state directly (apart from simply beginning with great
initial guesses, which isn't feasible for me once I move to 2D and have
more complicated boundary conditions) are certainly still welcome, but I
think this at least gives me a way to move forward.

Best,
Ray

On Wed, Nov 19, 2014 at 5:12 PM, Raymond Smith smit...@mit.edu wrote:

 Hi, FiPy.

 I'm trying to solve a number of problems involving (quasi-)neutral
 electrolyte solutions. I've attached some notes I made a while ago on one
 of the simplest example problems, for which there is an analytical
 solution. The equations are non-linear, and my typical boundary conditions
 are also often non-linear, as in the attached example.

 I'm curious if there are suggestions for improving convergence for this
 type of problem in FiPy. In the simple example (code here
 https://gist.github.com/raybsmith/6e724e6f1ecc98f9303b and below), we
 have the analytical solution for any applied voltage. For moderate positive
 applied voltages, I can get convergence -- up to around 10 thermal volts.
 However, for negative applied voltages (which corresponds to depletion of
 the salt in the electrolyte), I only get convergence down to about -0.5
 thermal volts, which is notably less than I'd like to simulate. I realize
 that sometimes solving non-linear problems numerically is simply difficult,
 but perhaps there are some things I could try that someone here could
 recommend?

 I've thought about doing the transient version of the problem and stepping
 toward a steady state with more extreme applied potentials. However,
 because the electrostatic part is static, one of the equations will always
 be stationary, even in the dynamic case. (avoiding electrodynamics here) So
 in some similar problems, I've still found I can't apply large potentials.

 Thanks,
 Ray

 import fipy as fp

 Lx = 1.
 nx = 1000
 dx = Lx/nx
 mesh = fp.Grid1D(nx=nx, dx=dx)
 phi = fp.CellVariable(name=elec. potential, mesh=mesh, value=0.)
 c_m = fp.CellVariable(name=conc., mesh=mesh, value=1.0)

 #V = -6e-1
 V = -4e-1
 #V = 4.0e0
 c0 = 1.
 c_m.setValue(c0)
 D_m = 2.03e-9 # m^2/s
 D_p = 1.33e-9 # m^2/s
 D_p = D_p/D_m
 D_m = 1.

 eq1 = (0 == fp.DiffusionTerm(coeff=-D_m, var=c_m) +
 fp.DiffusionTerm(coeff=D_m*c_m, var=phi))
 eq2 = (0 == fp.DiffusionTerm(coeff=-(D_p - D_m), var=c_m) +
 fp.DiffusionTerm(coeff=-(D_p + D_m)*c_m, var=phi))

 phi.constrain(0., mesh.facesLeft)
 phi.constrain(V, mesh.facesRight)
 c_m.constrain(c0, mesh.facesLeft)
 c_m.faceGrad.constrain(c_m.faceValue*phi.faceGrad, mesh.facesRight)

 Xcc = mesh.cellCenters[0]
 J = 1 - fp.numerix.exp(V)
 c_m_analyt = fp.CellVariable(name=analyt. conc., mesh=mesh)
 c_m_analyt.setValue(1 - J*Xcc)
 phi_analyt = fp.CellVariable(name=analyt. phi, mesh=mesh)
 phi_analyt.setValue(fp.numerix.log(1-J*Xcc))

 eq = eq1  eq2

 res = 1e10
 nIt = 0
 while res  1e-6:
 res = eq.sweep()
 print Iteration:, nIt, Residual:, res
 nIt += 1

 anstol = 1e-3
 print c_m.allclose(c_m_analyt, rtol=anstol, atol=anstol)
 print phi.allclose(phi_analyt, rtol=anstol, atol=anstol)



___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: Diffusion in a 2D Annulus.

[Raymond Smith]

Hey, Kyle.

Given inner radius rI and outer radius rO, perhaps you could do something
like

rMid = 0.5 * (rI + rO)
Xfc, Yfc = mesh.faceCenters
innerFaces = (mesh.exteriorFaces  (Xfc**2 + Yfc**2  rMid))
outerFaces = (mesh.exteriorFaces  (Xfc**2 + Yfc**2  rMid))

Then constrain using where=innerFaces and where=outerFaces.

Hopefully that might help.

Ray

On Wed, Nov 19, 2014 at 5:48 PM, Kyle Briton Lawlor klawlor...@gmail.com
wrote:

 Hi again FiPy,

 What is the easiest/most efficient way to make the distinction between the
 “outer” exterior faces and the “inner” exterior faces of the annulus?
 I would like to enforce to different boundary conditions at the outer and
 inner radii.

 Can I use the where=() function? Or something like it within the
 phi.constrain() function?
 I’ve seen this used to specify different coefficient regions in the
 1D-Diffusion example.
 Which was used as D.setValue(0.1, where=(specify some ranges)).

 Look forward to hearing more. :)

 Best, Kyle.

 On Nov 19, 2014, at 3:43 PM, Kyle Briton Lawlor klawlor...@gmail.com
 wrote:

  Thanks Daniel and Jonathan,
  All of this information helps very much.
  It is now working properly, so here is the code.
  In case it is useful to anyone.
  It is basic, but I figure its nice to share.
 
  cellSize=0.05
  radius=1.
  from fipy import *
  mesh = Gmsh2D('''cellSize = %(cellSize)g;
 radius = %(radius)g;
 Point(1) = {0, 0, 0, cellSize};
 Point(2) = {radius, 0, 0, cellSize};
 Point(3) = {0, radius, 0, cellSize};
 Point(5) = {-radius, 0, 0, cellSize};
 Point(6) = {0, -radius, 0, cellSize};
 Circle(1) = {2, 1, 3};
 Circle(2) = {3, 1, 5};
 Circle(3) = {5, 1, 6};
 Circle(4) = {6, 1, 2};
 Line Loop(1) = {1, 2, 3, 4};
 Point(200) = {radius/3, 0, 0, cellSize};
 Point(300) = {0, radius/3, 0, cellSize};
 Point(500) = {-radius/3, 0, 0, cellSize};
 Point(600) = {0, -radius/3, 0, cellSize};
 Circle(100) = {200, 1, 300};
 Circle(200) = {300, 1, 500};
 Circle(300) = {500, 1, 600};
 Circle(400) = {600, 1, 200};
 Line Loop(100) = {100, 200, 300, 400};
 Plane Surface(1)={1,100};
 ''' % locals())
  phi=CellVariable(name=solution,mesh=mesh,value=0.)
  D=1
  eq=TransientTerm()==DiffusionTerm(coeff=D)
  X,Y=mesh.faceCenters
  phi.constrain(1,mesh.exteriorFaces)
  timeStepDuration=10*0.9*cellSize**2/(2*D)
  steps=10
  for step in range(steps):
 eq.solve(var=phi,dt=timeStepDuration)
 
 TSVViewer(vars=(phi)).plot(filename=diffusion.annulus.+str(step)+.tsv”)
 
  Thanks!
  Kyle.
 
  On Nov 19, 2014, at 2:16 PM, Guyer, Jonathan E. Dr. 
 jonathan.gu...@nist.gov wrote:
 
  Just to clarify: as far as FiPy is concerned, interiorFaces have a
 cell on either side of them, whereas exteriorFaces have a cell on only
 one side; the other side is outside the mesh. It doesn't matter whether
 the mesh has simple topological connectivity (square, sphere, etc.) or
 complex topological connectivity (torus). interiorFaces separate mesh
 from more mesh; exteriorFaces separate mesh from not mesh.
 
 
  On Nov 19, 2014, at 10:53 AM, Daniel Wheeler daniel.wheel...@gmail.com
 wrote:
 
  Hi Kyle,
 
  #Boundary Conditions
  phi.constrain(X,mesh.exteriorFaces)
  phi.constrain(X,mesh.interiorFaces)
 
  The above is the problem. You are constraining the internal faces,
 which makes no sense in FiPy. I am not even sure how FiPy behaves when that
 constraint is added. However, I assume that is not what you want to do. If
 you remove that constraint the result seems to look nice. BTW,
 exteriorFaces refers to both the inner and outer external faces of the
 annulus. Perhaps you thought interiorFaces referred to the interior
 annulus faces that are external to the mesh.
 
  Also, you don't need to import the results into Mathematica. You can
 just use viewer = fipy.Viewer(phi); viewer.plot() to look at the results.
 
  Cheers,
 
  Daniel
 
  --
  Daniel Wheeler
  ___
  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 ]

Re: Unexpected result, possibly wrong, result solving 1D unsteady heat equation with spatially-varying diffusion coefficient

[Raymond Smith]

Thanks, Daniel. Congrats on making the switch! I submitted a pull request,
but I'm not sure everything is done correctly; feel free to follow up. In
particular, I didn't know how to test to make sure all my equations/images,
etc. would be rendered correctly by the mesh1D.py-browser process.

Ray

On Thu, Sep 25, 2014 at 12:16 PM, Daniel Wheeler daniel.wheel...@gmail.com
wrote:

 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 smit...@mit.edu 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 smit...@mit.edu 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 raymondbarrettsm...@gmail.com
 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. 
 jonathan.gu...@nist.gov 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

Mailing list archive link

[Raymond Smith]

Hi, Fipy.

I just noticed that when I click on the mailing list archive links at

http://www.ctcms.nist.gov/fipy/documentation/MAIL.html

I get a redirection error:

You have selected a link that connects to a redirection script on this
server, but the provided target URL is either malformed or is not approved
for redirection.

The provided url is in the location bar of this page, after url=. If it
makes sense to you, you are welcome to try it out. You are hereby advised
that the target site is not on this server.

Please contact the author of the referring page to advise them of the
problem.

Is this just me or is there a dead link that needs to be fixed?

Cheers,
Ray
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: Accuracy problem on results

[Raymond Smith]

Hi, John.

I'm not completely sure if this is related, but I do want to make sure it's
not something like that discussed in the recent thread about the heat
equation on this list:

http://thread.gmane.org/gmane.comp.python.fipy/3556/focus=3565

Is your diffusion coefficient (typically called alpha, thermal
diffusivity, in heat transfer) a function of position? If so, it's
important to keep only the thermal conductivity, k, as the coefficient for
DiffusionTerm (more detail in the linked thread). The density and heat
capacity should be coefficients for the TransientTerm.

Cheers,
Ray

On Tue, Sep 9, 2014 at 10:00 AM, John Assael iass...@gmail.com wrote:

 Hi,
 I keep getting inaccurate results on my heat conduction experiment, 0.1
 missing accuracy except if the shape is extremely and abnormally dense.
 In the beginning I though its the time steps but they are already very
 small and make no difference if I shorten the intervals.

 I am trying to replicate a simulation from COMSOL and if I take the mesh
 and the timesteps I cannot replicate the results accurately.

 I believe the problem is that all my variables are cell variables.
 However, I read in the FAQ that diffusion terms have to be FaceVariables.

 More specifically, my mesh is 2-D and I have the following equation and
 variables:

 *dT = CellVariable(name=r'$\Delta T$', mesh=mesh, value=0.)*
 *Q = CellVariable(name=r'$Q$', mesh=mesh, value=0.)*
 *T = CellVariable (name=r'$T$', mesh=mesh, value=0.)*
 *D = CellVariable (name=r'$D$', mesh=mesh, value=0.)*

 *heat = (TransientTerm(coeff=D) == DiffusionTerm(coeff=T) + Q)*

 And then I go to the areas I have defined and set the values. Here is an
 example of two of them.


 *x, y = mesh.cellCentersmeshW = (x ** 2 + y ** 2 =
 properties['r'] ** 2)*
 *Q.setValue(properties['Q'], where= meshW)*
 *T.setValue(properties['T'], where= meshW)*
 *D.setValue(properties['D'], where= meshW)*

 *meshW1 = (x ** 2 + y ** 2  properties['r'] ** 2)*
 *T.setValue(properties['T1'], where= **meshW1**)*
 *D.setValue(properties['D2'], where= **meshW1**)*

 However, when I do that I get an error for the face variable for every
 single one I try to change.
 Would sweeps instead of solve help me solve the problem?
 Do you have any suggestions?

 Thank you very much!

 Best regards,
 John


 --
 ..:: ic3man.gr ::..
 software | design | development
 _/_/_/ email : iass...@gmail.com
 _/_/_/ www   : http://www.ic3man.gr
 --
 This e-mail and any attachments are confidential. You may not copy
 or disseminate any information  contained in them  to anyone other
 than the intended recipient. If you are not the intended recipient
 please contact the sender by reply  e-mail and destroy all  copies
 of the original message immediately.
 --

 ___
 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 ]

Re: Accuracy problem on results

[Raymond Smith]

Perhaps your errors in assigning the values are coming from the initial
creation of the coefficients. I do recommend looking over the 1D diffusion
example
http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.mesh1D.html
carefully. Note that when declaring the coefficient for a diffusion term,
it should be declared as a FaceVariable, not a CellVariable, as you've
noted.

Ray

On Tue, Sep 9, 2014 at 10:50 AM, John Assael iass...@gmail.com wrote:

 Dear Raymond,

 Thank you for your fast reply.

 No I don't have that problem as my T variable represents the thermal
 conductivity k,  and it is not described in terms of position, it just
 changes from material to material.

 Exactly rho * cp are the transient terms.

 Finally the equation results dT the temperature difference.

 What else could it be?

 Best regards,
 John
 On 9 Sep 2014 17:16, Raymond Smith smit...@mit.edu wrote:

 Hi, John.

 I'm not completely sure if this is related, but I do want to make sure
 it's not something like that discussed in the recent thread about the heat
 equation on this list:

 http://thread.gmane.org/gmane.comp.python.fipy/3556/focus=3565

 Is your diffusion coefficient (typically called alpha, thermal
 diffusivity, in heat transfer) a function of position? If so, it's
 important to keep only the thermal conductivity, k, as the coefficient for
 DiffusionTerm (more detail in the linked thread). The density and heat
 capacity should be coefficients for the TransientTerm.

 Cheers,
 Ray

 On Tue, Sep 9, 2014 at 10:00 AM, John Assael iass...@gmail.com wrote:

 Hi,
 I keep getting inaccurate results on my heat conduction experiment, 0.1
 missing accuracy except if the shape is extremely and abnormally dense.
 In the beginning I though its the time steps but they are already very
 small and make no difference if I shorten the intervals.

 I am trying to replicate a simulation from COMSOL and if I take the mesh
 and the timesteps I cannot replicate the results accurately.

 I believe the problem is that all my variables are cell variables.
 However, I read in the FAQ that diffusion terms have to be FaceVariables.

 More specifically, my mesh is 2-D and I have the following equation and
 variables:

 *dT = CellVariable(name=r'$\Delta T$', mesh=mesh, value=0.)*
 *Q = CellVariable(name=r'$Q$', mesh=mesh, value=0.)*
 *T = CellVariable (name=r'$T$', mesh=mesh, value=0.)*
 *D = CellVariable (name=r'$D$', mesh=mesh, value=0.)*

 *heat = (TransientTerm(coeff=D) == DiffusionTerm(coeff=T) + Q)*

 And then I go to the areas I have defined and set the values. Here is an
 example of two of them.


 *x, y = mesh.cellCentersmeshW = (x ** 2 + y ** 2 =
 properties['r'] ** 2)*
 *Q.setValue(properties['Q'], where= meshW)*
 *T.setValue(properties['T'], where= meshW)*
 *D.setValue(properties['D'], where= meshW)*

 *meshW1 = (x ** 2 + y ** 2  properties['r'] ** 2)*
 *T.setValue(properties['T1'], where= **meshW1**)*
 *D.setValue(properties['D2'], where= **meshW1**)*

 However, when I do that I get an error for the face variable for every
 single one I try to change.
 Would sweeps instead of solve help me solve the problem?
 Do you have any suggestions?

 Thank you very much!

 Best regards,
 John


 --
 ..:: ic3man.gr ::..
 software | design | development
 _/_/_/ email : iass...@gmail.com
 _/_/_/ www   : http://www.ic3man.gr
 --
 This e-mail and any attachments are confidential. You may not copy
 or disseminate any information  contained in them  to anyone other
 than the intended recipient. If you are not the intended recipient
 please contact the sender by reply  e-mail and destroy all  copies
 of the original message immediately.
 --

 ___
 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 ]

Re: Unexpected result, possibly wrong, result solving 1D unsteady heat equation with spatially-varying diffusion coefficient

[Raymond Smith]

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 conor.flem...@eng.ox.ac.uk
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:* fipy-boun...@nist.gov [fipy-boun...@nist.gov] on behalf of Conor
 Fleming [conor.flem...@eng.ox.ac.uk]
 *Sent:* 08 August 2014 17:24
 *To:* fipy@nist.gov
 *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:* fipy-boun...@nist.gov [fipy-boun...@nist.gov] on behalf of Kris
 Kuhlman [kristopher.kuhl...@gmail.com]
 *Sent:* 08 August 2014 17:04
 *To:* fipy@nist.gov
 *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 conor.flem...@eng.ox.ac.uk
 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 - 

Re: Unexpected result, possibly wrong, result solving 1D unsteady heat equation with spatially-varying diffusion coefficient

[Raymond Smith]

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. 
jonathan.gu...@nist.gov 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 smit...@mit.edu 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 
 conor.flem...@eng.ox.ac.uk 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

Re: Unexpected result, possibly wrong, result solving 1D unsteady heat equation with spatially-varying diffusion coefficient

[Raymond Smith]

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 smit...@mit.edu 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 raymondbarrettsm...@gmail.com
 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. 
 jonathan.gu...@nist.gov 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 smit...@mit.edu wrote:

  Hi, Jonathan.
 
  I didn't know how

Re: 2 moving boundaries problem

[Raymond Smith]

Hey Olivier,

It seems like this is a phase separation problem. One option is to try to
model infinitesimal moving interfaces as you seem to be thinking now, with
some sort of Stefan condition to dictate their motion. However, another
approach which is useful in modeling phase separating systems is to start
with a Cahn-Hilliard style equation to model the dynamics of phase
separation. Have a look at the Phase Field examples, which show some phase
separation using what are referred to as non-conserved order parameters
(e.g. a parameter that goes from 0 if liquid to 1 if solid). Then, consider
the Cahn-Hilliard examples, including 
http://www.ctcms.nist.gov/fipy/examples/cahnHilliard/generated/examples.cahnHilliard.mesh2D.html;,
which is an alternate implementation of the first Cahn-Hilliard example
(it's linked to via an image on the main FiPy page, but seems not to be
included in the main Examples page...?). Cahn-Hilliard equations are
designed to model a conserved order parameter, such as concentration of a
species, which is more what you're considering.

Best,
Ray


On Wed, Dec 11, 2013 at 10:50 AM, Olivier DEZELLUS 
olivier.dezel...@univ-lyon1.fr wrote:

  Hello,

 I am a completely newbie in the use of FiPy but I tried to understand by
 using the different examples and particularly the ones on diffusion
 problems. However without any success sinces several weeks... this is the
 reason of my email to the list to find some help from more qualified
 users...

 My diffusion problem is summarized in the attached pdf file.

- An alpha(mostly A with a limited solubility of B)+beta(pure B) alloy
is in contact with a third element and a reaction occurs.
- The thickening of the tau reaction layer, that is rich in element B,
induces a depletion in B of the two phase layer
- Depletion in B leads to the formation of a pure alpha layer
- As the process continues, tau and alpha layers thickens

 From experimental study I have the rate of growth (constants Ktau and Kas)
 of the Tau and alpha layers.

 I would like to know if this diffusion process could be simulated by using
 FiPy ?  The quantity I am interested in is the concentration in B at the
 left hand side of the alpha layer, at interface with the tau layer.

 I tried to simplify the problem by using only one moving boundary (the
 alpha/alpha+beta) and giving a flux condition on the tau/alpha interface bu
 without any success.

 Does anyone have an idea on how to model this in FiPy ?

 Regards,
 Olivier


  --

 -
 La France est en 5e position pour le nombre de publications et en 26e pour le 
 financement de la recherche académique.
 Notre système de recherche est-il inefficace ou insuffisamment financé ?
 http://recherche-en-danger.apinc.org/

 --
 Dr Olivier DEZELLUS
 Laboratoire des Multimatériaux et Interfaces (UMR 5615)http://lmi.cnrs.fr/
 43, Bd du 11 Novembre 1918
 69622 Villeurbanne Cedex
 Tel: 04 72 44 83 86
 Fax: 04 72 44 06 18



 ___
 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 ]

Re: non-linear boundary conditions

[Raymond Smith]

Hi Fridolin,

I've used FiPy with non-linear boundary conditions before, e.g. (quite like
yours)
n\cdot\nabla c = c/(c+k)
in which k is a constant and c is an independent variable.

You can probably take care of the geometry with Gmsh, which is referenced
in the install documentation for FiPy (and also in the comprehensive
tutorials).

Best of luck,
Ray


On Tue, Dec 10, 2013 at 7:36 AM, Fridolin Gross fridolin.gr...@hu-berlin.de
 wrote:

 Hello,

 I’m looking for a PDE tool for a biological model describing a population
 of cells that interact via a diffusive messenger molecule.
 Before learning to use fipy, I’d like to know if it is the right tool for
 my purposes.

 I have a rectangular domain, intersected with circular regions (=cells).
 The interior of the cells is not part of the domain. In the intercellular
 region there is diffusion of the messenger molecule,
 while reactions take place on the circular boundaries (=cell membranes)
 that bind or degrade the molecules.
 The model requires a non-trivial Neumann boundary condition on those
 circular boundaries, with flux proportional to u/(k+u), where u is the
 concentration of messenger at the boundary.
 I was disappointed to find out that the matlab pdetoolbox cannot deal with
 parabolic problems with boundary conditions that depend non-linearly on the
 solution.
 Would this be possible in fipy?

 Thanks a lot for any kind of help!

 Fridolin
 ___
 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 ]

Re: Convection terms and algebraic equations

[Raymond Smith]

For the algebraic thing, could you do something like
eq = ( fp.TransientTerm(coeff=0., var=whatever) == \sum{z_i*C_i} )?


On Wed, Jul 24, 2013 at 3:35 PM, Edwin Sze Lun Khoo edwin...@mit.eduwrote:

 Hello,

 I have a couple of quick questions about the use of convection terms and
 defining algebraic equations of dependent variables in FiPy:

 1) What is meant by an inlet or an outlet boundary condition in
 http://www.ctcms.nist.gov/fipy/documentation/USAGE.html#applying-outlet-or-inlet-boundary-conditions
 ?

 2) What is the best way to write algebraic equations that relate
 CellVariables in FiPy? For instance, the electroneutrality assumption in
 electrochemistry, \sum{z_{i}C_{i} = 0, relates ion concentrations with an
 algebraic constraint.

 Best,
 Edwin



 ___
 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 ]

Re: divergence operator

[Raymond Smith]

Hi Jonathan,

Thanks for the reply, and you're right. Most of the simulations I'd been
doing involved only positive quantities like concentration, so I was safe
using harmonicFaceValue, but my displacement variables u1 and u2 can be
positive or negative. I just re-ran it using divergence operators on
arithmetic face values of displacements to confirm, and with that it works
fine.

Thanks,
Ray


On Tue, Jun 25, 2013 at 12:16 PM, Jonathan Guyer gu...@nist.gov wrote:


 On Jun 24, 2013, at 4:03 PM, Raymond Smith smit...@mit.edu wrote:

 I'm having confusion about what the divergence operator is doing in my
 situation. I ran two simulations and got very different results using two
 similar scripts, attached. The scripts are relatively long, but I'm
 including the full scripts in case there's anything else that someone may
 end up wanting to reference. My focus now is on the diff of the two
 scripts (lines 145-174 and lines 196-225). Basically, I'm not sure why
 they're different. I'm in 2 dimensions with a cellVariable called c, and
 the term I want in my equation is
 elem1 * dc/dx + elem2 * dc/dy,
 where elem1 and elem2 are known, constant parameters.


 My guess would have been that the treatments should be equivalent (except
 perhaps in the orders of accuracy?), but the results (movie_DIV.mpg for
 *divop.py and movie_DOT.mpg for the *dotop.py file) are very dissimilar.


 Is there something I'm missing about how the divergence operator is acting
 on the FaceVariables?


 It appears to me that the difference is not DIV vs. DOT, per se, but
 rather your use of the harmonicFaceValue for fields that are not uniformly
 positive. The harmonic average is only meaningful for quantities that are
 greater than zero.


  import fipy as fp
  from fipy import numerix as nx
  mesh = fp.Grid2D(nx=100, Lx=4 * nx.pi, ny=100, Ly=4 * nx.pi)
  x, y = mesh.cellCenters[0], mesh.cellCenters[1]
  c = y * nx.sin(y) * nx.cos(x)
  fp.Viewer(c).plot(filename=field.png)

  elem1 = 3.
  elem2 = 5.

 Note, there's no reason to make the nhats and dot them. Just ask for the
 scalar component you want:

  dotty = elem1 * c.grad[0] + elem2 * c.grad[1]
  fp.Viewer(dotty).plot(filename=dotty.png)


  a1 = fp.Variable((elem1, elem2))

 The harmonic mean produces some surprising artifacts:

  divvy_harmonic = (a1 * c.harmonicFaceValue).divergence
  fp.Viewer(divvy_harmonic).plot(filename=divvy_harmonic.png)

 The arithmetic mean does not:

  divvy_arithmetic = (a1 * c.arithmeticFaceValue).divergence
  fp.Viewer(divvy_arithmetic).plot(filename=divvy_arithmetic.png)

 The difference between the divergence of the arithmetic mean and the sum
 of the components is quite small:

  fp.Viewer(divvy_arithmetic - dotty).plot(filename=diff.png)


 Also, out of curiosity, examples.diffusion.anisotropy uses
 fp.numerix.NUMERIX.dot for some dotting. Is there a reason for the
 preference of that over fp.numerix.dot, and if so, how is it different?


 fp.numerix.dot performs dot products on each position of a field of
 vectors.
 fp.numerix.NUMERIX.dot (an alias for numpy.dot) does not understand how
 FiPy treats fields of vectors, but tries to dot the entire field as a
 single vector (or tensor).

 This is why we always stress to use the functions from fipy.numerix,
 rather than from numpy directly.

 numpy.dot is, nonetheless, a useful function for FiPy's internals.




 ___
 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 ]

Re: Fipy install test failures

[Raymond Smith]

Hi Johnathan,

I went ahead and re-updated my scipy and matplotlib as well, so a test
output with
numpy-1.7.0
pysparse-1.1.1
scipy-0.11.0
matplotlib-1.2.0
I got 6 errors
http://pastebin.com/BHadSKBE

A number of them seem to be related to gmsh, but I also do get the two
errors you mentioned.

I can rebuild gmsh in case there was something amiss having updated
numpy/scipy/matplotlib, but I'd be a bit surprised by that because gmsh
doesn't depend on any of those. I'm also available to try other things you
think may be helpful.

Cheers,
Ray




On Wed, Mar 20, 2013 at 4:54 PM, Jonathan Guyer gu...@nist.gov wrote:


 On Feb 22, 2013, at 2:27 PM, Raymond Smith wrote:

  So as far as I can tell, you were exactly right, and numpy-1.7.0 was
 causing the issues. If there's anything else I can do that would help
 clarify what's going on or help ease the transition to numpy-1.7.0, let me
 know.

 Can you try again with the released version of NumPy 1.7.0?

 I see from your log at http://pastebin.com/grn7Maki that you were running
 numpy version 1.7.0b2 when you got all those errors. I have not yet tested
 on Ubuntu, but on my Mac OS X machine I get only two of those errors
 (IndexError: index 34417792 is out of bounds for axis 0 with size 2 and
 ImportError: cannot import name _formatInteger), plus a third error that
 seems unique to my machine (and which I also get with NumPy 1.6.1).

 The first two issues are filed as http://matforge.org/fipy/ticket/556 and
 http://matforge.org/fipy/ticket/557.
 ___
 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 ]

Re: Time dependent spatially varying 2D source term representation

[Raymond Smith]

Hi Serbulent,

I was curious about the same thing, so I gave it a shot. I'm not sure if
what I've done is a good way to do things, but I think I managed to get a
source term that's a function of x, y. I modified your script in a number
of ways, and wrote it so that myFunc takes positions as arguments:

1) dx should be a float (dx = 1., not dx =1)
2) As far as I know, it's better to use fipy's numerix package than numpy.
3) I think for a source term that's dependent on position, you can use a
CellVariable. So I declared a CellVariable, which I then modify values of
at each time step.
4) Rather than loop through x and y and then trying to find the indexes
that correspond to the a value close to that, you can just loop through all
the indexes of the grid points in the problem. You can get the X and Y
arrays which are indexed by this total number of points from
mesh.CellCenters, then loop through all your points (nx*ny in the
rectangular 2D case).
5) For your time dependence, when you're calling myFunc, you need to have
calculated what the current time is before you call it (t =
step*timeStepDuration)
6) I made an arbitrary function of x, y, t for myFunc. The way I cast it,
it takes _values_ of position and time as arguments. The way you wrote it,
x and y need to be indexes of your tmp array.

File at http://pastebin.com/sZLHBTRV

Cheers,
Ray


On Tue, Mar 19, 2013 at 9:43 AM, Serbulent UNSAL serbule...@gmail.comwrote:

 Hi,

 First thanks to FiPy community for such a practical tool. I'm a newbie
 in FiPy and try to solve a problem based on source term.

 My equation is a 2D diffusion equation which has a source term as a
 function depends on x,y cordinates and time t

\frac{ \partial c( \vec{x},t)}{\partial t}  =
 D_{c}\bigtriangledown ^{2} c( \vec{x},t)} - f_{c}( \vec{x},t)

 I am googleing about 3 days found some clues [1] but they could'nt
 solve my problem. My code is working with a constant source term but I
 still couldn't find how can I represent it as a function depends on
 x,y and t.

 I'm trying to get my source term values from a mock up matrix for sake
 of simplicity. I'll be appriciate for any ideas/solutions. My code can
 be seen at http://pastie.org/6605575

 Thanks,

 Serbulent


 [1]

 http://bb10.com/python-fipy/2012-12/msg00016.html
 http://osdir.com/ml/python.fipy/2007-10/msg9.html
 http://bb10.com/python-fipy/2011-05/msg00030.html
 ___
 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 ]

Re: Point Source

[Raymond Smith]

Hi Adrian,

I'm not sure if this is what you're going for, but perhaps you could use an
implicit source to apply the internal boundary condition as outlined in
http://www.ctcms.nist.gov/fipy/documentation/USAGE.html#applying-internal-boundary-conditions.
If you're just trying to keep a particular point at a fixed value, maybe
that will work?

Cheers,
Ray


On Fri, Mar 1, 2013 at 12:11 PM, Adrian Jacobo ajac...@mail.rockefeller.edu
 wrote:

  Hi All,

  I'm trying to solve the diffusion equation with a point source at an
 arbitrary location inside a circle. My idea is to adapt the example
 examples.diffusion.circlehttp://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.circle.html#module-examples.diffusion.circlebut
  I'm not sure how to proceed.
  I think the idea should be create a variable which is zero everywhere
 except for where the point source is located, and add this variable as the
 source. But I don't know how to do this.

 Any ideas or help will be greatly appreciated.

 Best,
 Adrian.


 ___
 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 install test failures

[Raymond Smith]

Hi all,

I've just installed fipy, and I got a number of errors while running the
test. My output is located at http://pastebin.com/grn7Maki.
I'm using Arch Linux, 64 bit, with fipy-3.0, and I've installed pysparse
(1.1.1), gmsh (2.6.1), numpy (1.7.0), scipy (0.11.0), matplotlib (1.2.0),
mayavi (4.2.0).

I looked through some of the tests an couldn't find a clear pattern for
which tests succeeded or failed. I tried going to an interactive python
prompt and executing the code in the docstrings which is being executed in
the tests. I got the same errors, but I also got errors trying to do that
in tests that did not fail, so I didn't seem to be copying the tests
accurately at the prompt.

In particular, in the modules,
* models.levelSet.electroChem
* meshes.mesh.Mesh
* fipy.meshes.gmshMesh.MSHFile
* fipy.meshes.abstractMesh.AbstractMesh
I got
ValueError: operands could not be broadcast together with shapes (2)
(2,640) (2)

And I also got a few other errors, including a couple of
IndexError: axis 1 out of bounds [0, 1)
in a few seemingly random other places.

Any thoughts?

Thanks,
Raymond
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]