Re: bug in solvers/scipy/linearLUSolver.py

[James Pringle]

update -- as of scipy 0.19, linalg.splu had a depreciated option "drop_tol"
which did nothing:
https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/generated/scipy.sparse.linalg.splu.html

As of version 1.0, the keyword seems to have been dropped entirely, and
should be removed from linearLUSolver.py. It had not done anything, and so
can be removed.

This is the only place this keyword exists.

Jamie

On Thu, Nov 9, 2017 at 6:16 PM, James Pringle  wrote:

> Dear all --
>
> I just installed fipy from the conda channel with a fresh install of
> anaconda. It is based on python 2.7.14. When I run with the scipy solver it
> crashes with
>
>   File "/home/pringle/anaconda3/envs/MYFIPYENV27/lib/python2.7/
> site-packages/fipy/solvers/scipy/linearLUSolver.py", line 64, in _solve_
> permc_spec=3)
> TypeError: splu() got an unexpected keyword argument 'drop_tol'
>
>
> splu is loaded from scipy.sparse.linalg, and indeed it does not have the
> drop_tol keyword. When I remove the drop_tol keyword, it appears to run as
> intended, and produce good solutions for my problem.
>
> However, I am not sure that is the optimal solution...
>
> Jamie Pringle
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

bug in solvers/scipy/linearLUSolver.py

[James Pringle]

Dear all --

I just installed fipy from the conda channel with a fresh install of
anaconda. It is based on python 2.7.14. When I run with the scipy solver it
crashes with

  File
"/home/pringle/anaconda3/envs/MYFIPYENV27/lib/python2.7/site-packages/fipy/solvers/scipy/linearLUSolver.py",
line 64, in _solve_
permc_spec=3)
TypeError: splu() got an unexpected keyword argument 'drop_tol'


splu is loaded from scipy.sparse.linalg, and indeed it does not have the
drop_tol keyword. When I remove the drop_tol keyword, it appears to run as
intended, and produce good solutions for my problem.

However, I am not sure that is the optimal solution...

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

Re: fipy python 2.7

[James Pringle]

Worked fine.
Thank you, and apologize for the noise
Jamie

On Thu, Nov 9, 2017 at 3:44 AM, Benny Malengier 
wrote:

> I suppose the normal way adding  python=2.7
> to the conda arguments was tried and did not work?
>
> 2017-11-08 23:01 GMT+01:00 James Pringle :
>
>> Dear all -- (but mainly Daniel Wheeler, I guess)
>>
>>I am in the midst of finishing up a paper, but am trying to transition
>> to a new workstation. I have always installed fipy with
>>
>> conda create --name  --channel guyer --channel conda-forge fipy 
>> nomkl
>>
>>
>> and now, under anaconda, it installs fipy 3.1.3 with python 3.6.3. Is
>> there anyway to use conda to install it under python 2.7? I have nothing
>> against upgrading to python 3.6 -- but not as a finish a manuscript
>>
>> This is not urgent -- I have my old setup. But if there was any easy way
>> to use conda to install a new copy of fipy under python 2.7 until the paper
>> is out the door, that would be great.
>>
>> Thanks,
>> Jamie
>>
>> ___
>> fipy mailing list
>> fipy@nist.gov
>> http://www.ctcms.nist.gov/fipy
>> <https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=DwMFaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=vdnNYo9_ORW2mYc9Ku2lOxZGtzmS_sJzxpso3kcVYrw&s=YhNAIreEcwuy7sPlFRIzS1MCB7ZDlM5-e_GlVvhgHys&e=>
>>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy
>> <https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=DwMFaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=vdnNYo9_ORW2mYc9Ku2lOxZGtzmS_sJzxpso3kcVYrw&s=B2eHJ9lsS-wIo0YXhISTIEJ0UT_pdLbs3R6sMKn5OX4&e=>
>> ]
>>
>>
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

fipy python 2.7

[James Pringle]

Dear all -- (but mainly Daniel Wheeler, I guess)

   I am in the midst of finishing up a paper, but am trying to transition
to a new workstation. I have always installed fipy with

conda create --name  --channel guyer --channel conda-forge fipy nomkl


and now, under anaconda, it installs fipy 3.1.3 with python 3.6.3. Is there
anyway to use conda to install it under python 2.7? I have nothing against
upgrading to python 3.6 -- but not as a finish a manuscript

This is not urgent -- I have my old setup. But if there was any easy way to
use conda to install a new copy of fipy under python 2.7 until the paper is
out the door, that would be great.

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

Re: solvers for elliptic equations with dominant lower order terms

[James Pringle]

For anyone else who finds this thread with google. A couple of points that
helped me get over the issues I was having:

   1. Of the non-trilinos solvers, scipy was most robust in solving
   problems where the advective terms locally dominated over diffusive terms.
   (I could not get trilinos to run).
   2. I found that decomposing the domains into squares was more robust
   than decomposing it into triangles. This may be because I was computing the
   coefficients of equations with the gradient and divergence operators in
   fipy.
   3. However, *BY FAR *the most important thing I had to get right was
   making sure the spatial variation in the coefficients was well resolved by
   the grid/mesh. All of the things you learned when writing your own codes
   remains true when using someone else's black box solver If your problem
   is not well resolved, your solver will do something weird or fail. In my
   case I had issues with steep topography in some small parts of the in the
   realistic oceanographic problem I was solving; this will be familiar to
   ocean modelers. I blame the glaciers that carved the Laurentian Channel for
   my problems...

To find problem areas, I found it useful to add a transient term to the
equation (making, essentially, a heat equation) and solving it in time with
a simplified boundary condition of 1 on the boundaries and starting with
zero in the interior. As the solution evolved, it became easy to spot
trouble areas, since the "true" solution should stay bounded between 0 and
1 and should approach a steady solution of  1 everywhere eventually (as
long as diffusion exists everywhere; if this is not true your problem may
be under-constrained and certainly will be hard to solve). Where the
solution starts becoming >>1 or << 0 are good places to start thinking
about... especially if those locations are right next to each other. It is
also reassuring if the solution to the transient problem converges to the
steady solution in time.

Jamie

(p.s. how to solve the above heat equation problem with a BC of eta=1 on
all boundaries at the long time limit? Note that 1 is a solution of the
original elliptic equation above (as is any constant), and that it meets
the boundary condition. The elliptic equation is the steady state of the
corresponding heat equation. Thus 1 is a solution to the problem at the
limit of long time.)

On Thu, Jun 1, 2017 at 8:34 PM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> Unfortunately, James notified me offline of a breakage in PyTrilinos on
> Macs. I need to find time to debug.
>
> - Jon
>
> > On Jun 1, 2017, at 11:37 AM, Daniel Wheeler 
> wrote:
> >
> > Hi James,
> >
> > I don't have any hints for preconditioners, but using Trilinos's GMRES
> > and then systematically changing the GMRES parameters and
> > preconditioners might help you find one that works. I have done this
> > in the past, for example,
> >
> > https://urldefense.proofpoint.com/v2/url?u=https-3A__gist.
> github.com_wd15_c28ab796cb3d9781482b01fb67a7ec2d-23file-
> 2Dunsegregatedequation-2Dpy-2DL187&d=DwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=
> 7HJI3EH6RqKWf16fbYIYKw&m=GzwUdgPQ01r5uE6757k5b_5G60EzT_ebxic5VzxCPLw&s=
> ecDSV0HowRiX7RxqKQIM8k2pN6MWtkHeJVE0FsRUHKA&e=
> >
> > You can also get information about how well the solver is converging.
> >
> > Cheers,
> >
> > Daniel
> >
> > On Wed, May 31, 2017 at 3:31 PM, James Pringle  wrote:
> >>
> >> Dear all --
> >>
> >> I need to solve a second order PDE in 2D of the form
> >>
> >> A(x,y)*(eta_xx+eta_yy)+B(x,y)*eta_x+C(x,y)*eta_y = 0
> >>
> >> This is an advective diffusive balance, with B and C representing the
> >> advective velocities.
> >>
> >>
> >> Over much of the domain the low order terms dominate. For some
> geometries,
> >> the default solver I am using (SciPy LinearLUSolver) is having trouble
> find
> >> an answer (especially when the characteristics defined by B and C form
> >> loops.
> >>
> >> Does anyone have recommendations for how to choose a better solver? (I
> have
> >> blindly tried all the ones in SciPy, the ones besides LinearLUSolver do
> >> worse).
> >>
> >> And are there any hints on how to use preconditioners? There is
> discussion
> >> of preconditioners and Trilinnos, but I can't find hints on how one
> would
> >> choose to use them and when they might help.
> >>
> >> Pointers to the broader literature are welcome!
> >>
> >> Thank you,
> >> Jamie Pringle
> >>
> >> ___

Re: solvers for elliptic equations with dominant lower order terms

[James Pringle]

I may have been unclear. This is a single PDE in the variable eta. A(x,y),
B(x,y) and C(x,y) are coefficients of the PDE. They are indeed separate,
and fipy does a nice job of solving this sort of equation in many limits. I
am having problems when the lower order parts of the equation
(B(x,y)*eta_x+C(x,y)*eta_y)
dominate the solution and the characteristics defined by B and C form
closed loops.

To understand this better, imagine the equation as a steady state fluid
flow problem, with the A(x,y)*(eta_xx+eta_yy) representing diffusion (and
to be correct, it should be written "Del . (A(x,y) Grad eta)"), and the low
order terms representing advection, with B and C being the velocities doing
the advection. When the diffusion is negligible, eta should be constant
along the flow lines (characteristics) defined by B and C.  But if those
flow line form a closed loop, the diffusive term must determine the overall
solution inside the loop...

Thanks,
Jamie

On Wed, May 31, 2017 at 3:32 PM, Sergio Manzetti <
sergio.manze...@fjordforsk.no> wrote:

> James, are A(x,y), B(x,y) and C(x,y) three different variable functions
> for the PDE?
>
> FIPY can solve at least when these three are equal. If it solves when they
> differ, I am not sure.
>
> Sergio
>
>
> Sergio Manzetti
>
>
> <https://urldefense.proofpoint.com/v2/url?u=http-3A__www.fjordforsk.no_logo-5Fhr2.jpg&d=DwMFaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=INzNz0EhAfOz2J4uGWxrKjo-qBQYnUHjXOF9rSxEXTQ&s=U0xwV0Da6Bipu3Yzfcw9cPhJIwPF4QPA_zdG-Gms5uc&e=>
>
> Fjordforsk AS
> <https://urldefense.proofpoint.com/v2/url?u=http-3A__www.fjordforsk.no&d=DwMFaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=INzNz0EhAfOz2J4uGWxrKjo-qBQYnUHjXOF9rSxEXTQ&s=U6qdDuDFx7m_FBcXq8wrKgktDC12I3ar-kaEGfdQ4kA&e=>
>
> <https://urldefense.proofpoint.com/v2/url?u=http-3A__www.fjordforsk.no&d=DwMFaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=INzNz0EhAfOz2J4uGWxrKjo-qBQYnUHjXOF9rSxEXTQ&s=U6qdDuDFx7m_FBcXq8wrKgktDC12I3ar-kaEGfdQ4kA&e=>
> Midtun
> 6894 Vangsnes
> Norge
> Org.nr. 911 659 654
> Tlf: +47 57695621 <+47%2057%2069%2056%2021>
> Økolab
> <https://urldefense.proofpoint.com/v2/url?u=http-3A__www.oekolab.com&d=DwMFaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=INzNz0EhAfOz2J4uGWxrKjo-qBQYnUHjXOF9rSxEXTQ&s=K6jUoxb_DR3t0-8jNTHjwDpAu7N5bdMUKatMC-pun-U&e=>
> |  Nanofactory
> <https://urldefense.proofpoint.com/v2/url?u=http-3A__www.nanofact.no&d=DwMFaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=INzNz0EhAfOz2J4uGWxrKjo-qBQYnUHjXOF9rSxEXTQ&s=i8VNBHNUWx6dxM6qgwo0hBV4EVmqK-NBY5bRDJBmfAM&e=>
> |  AQ-Lab
> <https://urldefense.proofpoint.com/v2/url?u=http-3A__www.aq-2Dlab.no&d=DwMFaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=INzNz0EhAfOz2J4uGWxrKjo-qBQYnUHjXOF9rSxEXTQ&s=i8G56fQKIsGv9dk8ev8ozkaCLIsd-LOEHrqpX9cY5yI&e=>
> |  FAP
> <https://urldefense.proofpoint.com/v2/url?u=http-3A__www.phap.no&d=DwMFaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=INzNz0EhAfOz2J4uGWxrKjo-qBQYnUHjXOF9rSxEXTQ&s=m0sz0QVLjdTnHTiVYdJ0GZj1VHf1oMsqyQyjdWiq-_o&e=>
>
>
> --
> *From: *"James Pringle" 
> *To: *"fipy" 
> *Sent: *Wednesday, May 31, 2017 9:31:51 PM
> *Subject: *solvers for elliptic equations with dominant lower order terms
>
>
> Dear all --
>
> I need to solve a second order PDE in 2D of the form
>
> A(x,y)*(eta_xx+eta_yy)+B(x,y)*eta_x+C(x,y)*eta_y = 0
>
> This is an advective diffusive balance, with B and C representing the
> advective velocities.
>
>
> Over much of the domain the low order terms dominate. For some geometries,
> the default solver I am using (SciPy LinearLUSolver) is having trouble
> find an answer (especially when the characteristics defined by B and C form
> loops.
>
> Does anyone have recommendations for how to choose a better solver? (I
> have blindly tried all the ones in SciPy, the ones besides LinearLUSolver
> do worse).
>
> And are there any hints on how to use preconditioners? There is discussion
> of preconditioners and Trilinnos, but I can't find hints on how one would
> choose to use them and when they might help.
>
> Pointers to the broader literature are welcome!
>
> Thank you,
> Jamie Pringle
>
> ___
> 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 ]

solvers for elliptic equations with dominant lower order terms

[James Pringle]

Dear all --

I need to solve a second order PDE in 2D of the form

A(x,y)*(eta_xx+eta_yy)+B(x,y)*eta_x+C(x,y)*eta_y = 0

This is an advective diffusive balance, with B and C representing the
advective velocities.


Over much of the domain the low order terms dominate. For some geometries,
the default solver I am using (SciPy LinearLUSolver) is having trouble find
an answer (especially when the characteristics defined by B and C form loops
.

Does anyone have recommendations for how to choose a better solver? (I have
blindly tried all the ones in SciPy, the ones besides LinearLUSolver do
worse).

And are there any hints on how to use preconditioners? There is discussion
of preconditioners and Trilinnos, but I can't find hints on how one would
choose to use them and when they might help.

Pointers to the broader literature are welcome!

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

Re: an-isotropic mixing

[James Pringle]

Dear all --

   Please ignore the last message... for those who found the message by
googling, the answer is the first question in the fipy FAQ at
http://www.ctcms.nist.gov/fipy/documentation/FAQ.html .

   Well, that was a silly way to waste a day.

Jamie

On Mon, May 29, 2017 at 12:38 PM, Pringle, James 
wrote:

> Dear all --
>
> I am just double checking my sanity. If I have a diffusion problem
> with an-isotropic and spatially varying diffusion:
>
>   (d/dx)(A(x,y) dPsi/dx)+d/dy(B(x,y) dPsi/dy) = something
>
> The correct way to set up the diffusion term is code of the form
>
> diffCoef=fipy.CellVariable(mesh=mesh,rank=1) #define on center
> diffCoef.value[0,:]=A
> diffCoef.value[1,:]=B
> 
> eq=(fipy.DiffusionTerm(var=psiPart,coeff=diffCoef)+sourceTerm)
> eq.solve(var=psiPart,solver=TheSolver())
>
>
> Where A and B are vectors whose length is equal to the number of centers.
> I don't have to set up a complete tensor with 0's on the diagonal, do I?
> And if I do, How would I do so?
>
> I am just sanity checking some odd results.
>
> Thanks,
> Jamie
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

an-isotropic mixing

[James Pringle]

Dear all --

I am just double checking my sanity. If I have a diffusion problem with
an-isotropic and spatially varying diffusion:

  (d/dx)(A(x,y) dPsi/dx)+d/dy(B(x,y) dPsi/dy) = something

The correct way to set up the diffusion term is code of the form

diffCoef=fipy.CellVariable(mesh=mesh,rank=1) #define on center
diffCoef.value[0,:]=A
diffCoef.value[1,:]=B

eq=(fipy.DiffusionTerm(var=psiPart,coeff=diffCoef)+sourceTerm)
eq.solve(var=psiPart,solver=TheSolver())


Where A and B are vectors whose length is equal to the number of centers. I
don't have to set up a complete tensor with 0's on the diagonal, do I? And
if I do, How would I do so?

I am just sanity checking some odd results.

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

Re: bug in Neumann boundary condition to Laplacian

[James Pringle]

Daniel  --

Thank you, that makes a great deal of sense. On the boundary conditions you
are write about faceNorm; I had misinterpreted the text from
http://www.ctcms.nist.gov/fipy/documentation/USAGE.html

If the gradient normal to the boundary (*e.g.*, [image:
\hat{n}\cdot\nabla\phi]) is to be set to a scalar value of 2, use

>>> var.faceGrad.constrain(2 * mesh.faceNormals, where=mesh.exteriorFaces)


to mean that the gradient boundary conditions was defined normal and
outward of the solution domain. But it is clear from the text above that
this is not so; my bad.

Jamie





On Tue, May 9, 2017 at 1:55 PM, Daniel Wheeler 
wrote:

> On Fri, May 5, 2017 at 5:37 PM, James Pringle  wrote:
> > Jon --
> >
> > Thanks for this reply. There is something in it I don't understand.
> Why
> > is it necessary in fipy constrain the tangential derivative to the
> solution
> > to the Poisson equation in fipy?
>
> I don't think that the tangential derivative matters. I think that in
> James's version of "DemoBug.py" multiplying by the "faceNorms" is
> giving the wrong sign in some places from the normal component. By
> removing the "faceNorms", the normal component is correct. I think
> that you can zero out the tangential component and it has no impact.
> For example, the following also works,
>
> 
> faceNorm_x = faceNorm.copy()
> faceNorm_y = faceNorm.copy()
> faceNorm_x[:] = 0.
> faceNorm_y[:] = 0.
> faceNorm_x[0, :] = 1.
> faceNorm_y[1, :] = 1.
> psiInvert.faceGrad.constrain(psiTrue.faceGrad * faceNorm_y,
> where=yeq0bound)
> psiInvert.constrain(psiTrue.faceValue,where=yeq1bound)
> psiInvert.faceGrad.constrain(psiTrue.faceGrad * faceNorm_x,
> where=xeq0bound)
> psiInvert.faceGrad.constrain(psiTrue.faceGrad * faceNorm_x,
> where=xeq1bound)
> 
>
> In the above the tangential component is zeroed out and the normal
> component is unity, so the sign of the gradient is unaffected.
>
> --
> Daniel Wheeler
> ___
> fipy mailing list
> fipy@nist.gov
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.
> ctcms.nist.gov_fipy&d=DwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=
> 7HJI3EH6RqKWf16fbYIYKw&m=xq7jzzUAUm5bpCMpllA6E_eUfqWzoAEZQTfsDZ33klc&s=
> NdwgYlVTn0mfJ7DeAyliZl5Tk3V6adv7fSpuszW9Hls&e=
>   [ NIST internal ONLY: https://urldefense.proofpoint.
> com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_
> fipy&d=DwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=
> xq7jzzUAUm5bpCMpllA6E_eUfqWzoAEZQTfsDZ33klc&s=
> E958ChgiyZsgdzhMDnik1d5wRFV3zQVlGDphpMvParY&e=  ]
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: bug in Neumann boundary condition to Laplacian

[James Pringle]

Jon --

Thanks for this reply. There is something in it I don't understand. Why
is it necessary in fipy constrain the tangential derivative to the solution
to the Poisson equation in fipy?  In analytical solutions of Poisson's
equation one can solve the problem either by specifying the normal
derivative of the solution on the boundary OR the value of the function on
the boundary. To do both generally over-specifies the solution.

Specifying the tangential gradients is equivalent to specifying the
value of the solution on the boundary (to within a constant), for one can
integrate the tangential gradient along the boundary to find the solution
along the boundary.

Why is it in general consistent to specify both the normal and
tangential boundary conditions?

Thanks,
Jamie




On Fri, May 5, 2017 at 3:33 PM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> Jamie -
>
>
>
> Sorry it took so long to get back to this. I'd actually figured out what
> the issue was several days ago, but I haven't had time to verify and write
> anything up until now.
>
>
>
> The short answer is that you need to constrain the gradient of PsiInvert
> to be equal to the gradient of PsiTruth, but what you're doing (aside from
> being a lot of unnecessary typing and I'm pretty sure having low-order
> spacial accuracy) is only constraining the normal component of the boundary
> gradient. You need the tangential component, too.
>
>
>
> I've updated your script below with the proper boundary gradients, e.g.,
>
>
>
>   psiInvert.faceGrad.constrain(psiTrue.faceGrad,where=yeq0bound)
>
>
>
> and the recommended way to set omega to the Laplacian:
>
>
>
>   omega = psiTrue.faceGrad.divergence
>
>
>
> I get errors of O(1e-10) for all three cases.
>
>
>
> In the course of working through this, I realized that I'd not finished
> converting .faceGrad to math in what I posted before, such that I also
> neglected the tangential term. I've updated the notebook at
>
>
>
>   https://urldefense.proofpoint.com/v2/url?u=https-3A__gist.
> github.com_guyer_6b7699531f75b353f49a0cf36d683a
> a4&d=DwIGaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=
> mjWx1u6TsXKmtow-qItsKS4VE7r1nIMW-c395mi9wzM&s=EQ1vHADeK-
> 0hEbWdEYmXbmGl2msCVM8sAPo2PX5wb7w&e=
>
>
>
> to correct this and also to illustrate a *very hacky* implementation of
> boundary upwinding.
>
>
>
> - Jon
>
>
>
> > On Apr 26, 2017, at 9:02 AM, James Pringle  wrote:
>
> >
>
> > Dear all --
>
> >
>
> > I think I have discovered a bug in the implementation of Neumann
> boundary conditions in fipy for the solution of Poisson's equation. If it
> is not a bug, it exposes a significant inefficiency in the numerics.
> Attached is a test code (DemoBug.py) which illustrates the problem. The
> test code takes a given function PsiTruth, takes the the Laplacian of it,
> and then solves Possion's equation to recover Psi (we call this recovered
> version of Psi "PsiInvert," and it should ideally be exactly equal to
> PsiTruth):
>
> >
>
> > \Del^2 PsiInvert = \Del^2 PsiTruth
>
> >
>
> > The problem is solved in a rectangular domain of size (Lx,Ly) with the
> following boundary condition's:
>
> >   • PsiInvert=PsiTruth(x,Ly) on y=Ly (a Dirchelet BC)
>
> >   • The normal derivative of PsiInvert = The normal derivative of
> PsiTruth on all other boundaries.
>
> > The normal gradient of PsiTruth for the boundary condition can be solved
> either directly from PsiTruth or by using the fipy.faceGradAverage() on
> PsiTruth on the mesh. This choice is made in the code after the comment
> "#Choose Boundary Condition". I have three choices of the PsiTruth to
> examine. The error is defined as the STD(PsiInvert-PsiTruth). All have a
> peak amplitude of about 1.0, so the standard deviation of the error in the
> solution should be <<1.
>
> >   • An isolated gaussian that does not extend to the boundaries;
> this works very well with a STD(error) of 3.7e-4
>
> >   • cos(2*pi/wavelength*x)*cos(2*pi/wavelength*y), where wavelength
> is 160 the resolution of the grid (so is very well resolved). This works
> less well but ok with an error of 5.0e-2 (two orders of magnitude larger
> than above, about 5% of the solution magnitude).
>
> >   • abs(cos(2*pi/wavelength*x)*cos(2*pi/wavelength*y)), same wave
> length as above. Error is 4.8e-1; the error in the solution is 3 orders of
> magnitude greater than case 1 and an order of magnitude worse than case 2.
> It is 50% of the solution magnitude.
>

Re: gradients on faces

[James Pringle]

Jon --

Thank you for that clear explanation. What remains unclear to me is why
the .faceGrad.constrain() is not constraining the solution to Poisson's
equation properly. I have given a clearer example of that in a follow on
email titled "bug in Neumann boundary condition to Laplacian"

Thanks,
Jamie

On Wed, Apr 26, 2017 at 12:06 PM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> James -
>
> Thank you for your test script. It helped me understand what you were
> seeing and what behavior you were looking for. I've put together a notebook
> that I think concisely illustrates the source of the behavior you're
> encountering.
>
>   https://urldefense.proofpoint.com/v2/url?u=https-3A__gist.
> github.com_guyer_6b7699531f75b353f49a0cf36d683a
> a4&d=DwIGaQ&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=
> 8Otm5er2tke4qe_pfoFqjD7xv-xYoMTBdUD7HW01qxg&s=mIlZMfmcm_
> vIrOBVFcGZFNwOGxucRHOzR7wyx2D946U&e=
>
> The bottom line is that FiPy assumes zero gradient on all boundaries
> unless some constraint is applied. FiPy does not do any upwinding to
> project interior values and slopes out to the boundary.
>
> Further, Dirichlet constraints are reflected in boundary gradients, but
> Neumann constraints are not reflected in boundary values. Doing the latter
> (while desirable) would require unwinding (expensive) and care to avoid
> circular definitions (expensive). We'll think about how to achieve this and
> we'd certainly welcome code contributions from anybody who wants to tackle
> it.
>
> Finally, take care with your parentheses. .grad, .faceGrad, etc. are
> properties, not methods. Appending () to the end causes them to be
> evaluated to static numeric values (phi.grad() is the same as
> phi.grad.value).
>
> - Jon
>
> > On Apr 21, 2017, at 6:33 PM, James Pringle  wrote:
> >
> > To clarify and answer my question, and to give an analysis code for
> those interested in figuring these out, here is some text. (I would love
> comments from those who know more!)
> >
> > First a restatement of the question. To test the solvers, create a field
> PsiTruth, and solve for PsiInvert with
> >
> > \Del^2 PsiInvert = \Del^2 PsiTruth
> >
> > With various boundary conditions. The attached code does this. (To
> understand the notation in the code, see my prior email). My conclusions
> are that for maximum accuracy:
> >   • Compute \Del^2 PsiTruth in the interior on center points with
> the grad() operator. It is significantly more consistent with the solver
> numerics than the leastSquaresGrad() operator. This results in smaller
> error.
> >   • On the boundaries, use .faceGradAverage() to calculate the
> gradients of PsiTruth; .faceGrad() is incorrect on the boundaries.
> >   • At least part of one boundary must be of the form
> PsiInvert.constrain(), for otherwise the solution is arbitrary to within an
> additive constant, and that constant can be so large that it compromises
> the accuracy of the solution due to the finite number of significant digits
> in floating point numbers.
> >   • When gradient boundary conditions are specified, accuracy
> increases linearly with 1/resolution, and is much less than when the
> boundary values are fixed.
> > I hope this helps someone, and I would like to hear from others about
> the choice of gradient operators and how to achieve greater accuracy in
> calculating and specifying gradients on the boundary.
> >
> > In particular, I would love to have a gradient calculation that is
> exactly consistent with how the gradient is specified in the boundary
> condition constraint.
> >
> > It would also be nice to have a Laplacian operator that is consistent
> with the solver in the model -- or do I get that by using the .grad()
> operator twice?
> >
> > Thanks all,
> > Jamie
> >
> >
> > On Fri, Apr 21, 2017 at 11:28 AM, Pringle, James 
> wrote:
> > Dear all --
> >
> > I am using fipy to solve a series of fluid dynamics problems. As part of
> this, I need to compute a stream function Psi from the vorticity of flow
> omega=V_x-U_y.  U and V are calculated from the gradient of the solution to
> another set of elliptic equations, and are available on the same mesh as
> Psi will be calculated on.
> >
> > The equation relating the stream function Psi to the vorticity omega is
> >
> > \Del^2 Psi = omega
> >
> > And the boundary conditions on Psi are that the gradient normal to the
> boundary is equal to the flow parallel to the boundary. All very standard.
> >
> > I calculate velocities U and V from the gradient of another

Re: gradients on faces

[James Pringle]

To clarify and answer my question, and to give an analysis code for those
interested in figuring these out, here is some text. (I would love comments
from those who know more!)

First a restatement of the question. To test the solvers, create a field
PsiTruth, and solve for PsiInvert with

\Del^2 PsiInvert = \Del^2 PsiTruth


With various boundary conditions. The attached code does this. (To
understand the notation in the code, see my prior email). My conclusions
are that for maximum accuracy:

   1. Compute \Del^2 PsiTruth in the interior on center points with the
   grad() operator. It is significantly more consistent with the solver
   numerics than the leastSquaresGrad() operator. This results in smaller
   error.
   2. On the boundaries, use .faceGradAverage() to calculate the gradients
   of PsiTruth; .faceGrad() is incorrect on the boundaries.
   3. At least part of one boundary must be of the form
   PsiInvert.constrain(), for otherwise the solution is arbitrary to within an
   additive constant, and that constant can be so large that it compromises
   the accuracy of the solution due to the finite number of significant digits
   in floating point numbers.
   4. When gradient boundary conditions are specified, accuracy increases
   linearly with 1/resolution, and is much less than when the boundary values
   are fixed.

I hope this helps someone, and I would like to hear from others about the
choice of gradient operators and how to achieve greater accuracy in
calculating and specifying gradients on the boundary.

In particular, I would love to have a gradient calculation that is exactly
consistent with how the gradient is specified in the boundary condition
constraint.

It would also be nice to have a Laplacian operator that is consistent with
the solver in the model -- or do I get that by using the .grad() operator
twice?

Thanks all,
Jamie


On Fri, Apr 21, 2017 at 11:28 AM, Pringle, James 
wrote:

> Dear all --
>
> I am using fipy to solve a series of fluid dynamics problems. As part of
> this, I need to compute a stream function Psi from the vorticity of flow
> omega=V_x-U_y.  U and V are calculated from the gradient of the solution to
> another set of elliptic equations, and are available on the same mesh as
> Psi will be calculated on.
>
> The equation relating the stream function Psi to the vorticity omega is
>
> \Del^2 Psi = omega
>
>
> And the boundary conditions on Psi are that the gradient normal to the
> boundary is equal to the flow parallel to the boundary. All very standard.
>
> I calculate velocities U and V from the gradient of another variable (a
> pressure like term). My question is this: which of the routines that
> calculate the gradient of a fipy variable on the faces is most consistent
> with how the normal derivative is specified on the boundary? Which of
>  .faceGrad() or faceGradAverage() is most consistent with how the gradient
> is implemented in  Psi.faceGrad.constrain()? or some other gradient
> calculation? I ask because .faceGrad() does not seem to calculate the
> gradient on the boundary faces...
>
> Thanks,
> Jamie Pringle
>
> Director Ocean Process Analysis Laboratory in the Insitute for Earth,
> Ocean and Space
>
> Associate Professor of Earth Sciences
>
> University of New Hampshire
>
> http://oxbow.sr.unh.edu
>
from pylab import *
from numpy import *
import time
import matplotlib.tri as tri
import fipy

#this code creates and arbitrary streamfunction psi, finds velocity
#from it, and then inverts the velocity back to the stream function
#using the vorticity equation. The result should be the same as the
#original streamfunction, to within an additive constant.

#==
#define domain
dx=1.0
Lx=125.0
nx=int(round(Lx/dx)); ny=int(round(Lx/dx))
mesh=fipy.Grid2D(dx=dx,dy=dx,nx=nx,ny=ny)

assert max(mesh.faceCenters.value[0,:])==Lx,'dx must be multiple of Lx, or the BCs will fail'

#define streamfunction as cos(2*pi/wavelength*mesh.x)*cos(2*pi/wavelength*mesh.y) and plot
wavelength=80.0
psiTrue=fipy.CellVariable(name='PsiTrue',mesh=mesh,
  value=cos(2*pi/wavelength*mesh.x)*cos(2*pi/wavelength*mesh.y))
print 'wavelength/dx=',wavelength/dx

figure(1)
clf()
tricontourf(mesh.x,mesh.y,psiTrue.value,30)
colorbar()
title('True Psi')

#==
#calculate U and V from streamfunction with U=-Psi_y and V=Psi_x,
#and plot with quiver.
if False:
#calculate gradient with leastSquaresGrad
print 'Calculating gradient in interior with leastSquaresGrad'
jnk=psiTrue.leastSquaresGrad()
Vcenter=jnk[0,:]
Ucenter=-jnk[1,:]
#==
#Now calculate the vorticity omega=V_x-U_y
Ufipy=fipy.CellVariable(name='U',mesh=mesh,value=Ucenter); Ugrad=Ufipy.leastSquaresGrad()
Vfipy=fipy.CellVariable(name='V',mesh=mesh,value=Vcenter); Vgrad=Vfipy.leastSquaresGrad()

gradients on faces

[James Pringle]

Dear all --

I am using fipy to solve a series of fluid dynamics problems. As part of
this, I need to compute a stream function Psi from the vorticity of flow
omega=V_x-U_y.  U and V are calculated from the gradient of the solution to
another set of elliptic equations, and are available on the same mesh as
Psi will be calculated on.

The equation relating the stream function Psi to the vorticity omega is

\Del^2 Psi = omega


And the boundary conditions on Psi are that the gradient normal to the
boundary is equal to the flow parallel to the boundary. All very standard.

I calculate velocities U and V from the gradient of another variable (a
pressure like term). My question is this: which of the routines that
calculate the gradient of a fipy variable on the faces is most consistent
with how the normal derivative is specified on the boundary? Which of
 .faceGrad() or faceGradAverage() is most consistent with how the gradient
is implemented in  Psi.faceGrad.constrain()? or some other gradient
calculation? I ask because .faceGrad() does not seem to calculate the
gradient on the boundary faces...

Thanks,
Jamie Pringle

Director Ocean Process Analysis Laboratory in the Insitute for Earth, Ocean
and Space

Associate Professor of Earth Sciences

University of New Hampshire

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

Re: cauchy boundary condition

[James Pringle]

Dear all --

I would like to withdraw my question. My blithe assertion that there
had to be a solution was naive. It turns out in going from a stream
function to a pressure formulation, my boundary conditions had become
non-local.

Thus fipy was not giving me a solution because there was no solution --
which tends to give any solver fits. My apologies for wasted time.

Jamie

On Thu, Apr 13, 2017 at 3:10 PM, Pringle, James 
wrote:

> Dear all --
>
> I need to solve a second order PDE in 2D of the form
>
> A(x,y)*(eta_xx+eta_yy)+B(x,y)*eta_x+C(x,y)*eta_y = 0
>
> on a rectangular domain. On three sides there is a BC of no normal
> gradient. On the fourth side their is a Cauchy boundary condition:
>
> @x=0, eta_x=D(y)  *AND* eta=0 @x=0
>
> Lets assume that there is a solution that meets these boundary conditions
> (i.e. I have not screwed up). Should fipy be able to solve this by
> specifying (*both at the same place, *x0side)
>
> var.constrain(0., where=x0side)
> var.faceGrad.constrain(D_rank2, where=x0side)
>
> Or will this break the numerics? My actual problem is considerably more
> complex, and I want to check that this is even plausible before starting
> down this path. I can make this boundary simpler, but at the cost of
> considerable complexity elsewhere.
>
> Thanks,
> Jamie Pringle
>
> Director Ocean Process Analysis Laboratory in the Insitute for Earth,
> Ocean and Space
>
> Associate Professor of Earth Sciences
>
> University of New Hampshire
>
> http://oxbow.sr.unh.edu
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

cauchy boundary condition

[James Pringle]

Dear all --

I need to solve a second order PDE in 2D of the form

A(x,y)*(eta_xx+eta_yy)+B(x,y)*eta_x+C(x,y)*eta_y = 0

on a rectangular domain. On three sides there is a BC of no normal
gradient. On the fourth side their is a Cauchy boundary condition:

@x=0, eta_x=D(y)  *AND* eta=0 @x=0

Lets assume that there is a solution that meets these boundary conditions
(i.e. I have not screwed up). Should fipy be able to solve this by
specifying (*both at the same place, *x0side)

var.constrain(0., where=x0side)
var.faceGrad.constrain(D_rank2, where=x0side)

Or will this break the numerics? My actual problem is considerably more
complex, and I want to check that this is even plausible before starting
down this path. I can make this boundary simpler, but at the cost of
considerable complexity elsewhere.

Thanks,
Jamie Pringle

Director Ocean Process Analysis Laboratory in the Insitute for Earth, Ocean
and Space

Associate Professor of Earth Sciences

University of New Hampshire

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

Re: sources errors in advection/diffusion problems, and one solution

[James Pringle]

Thank you.
I am busy for the week, but then will get a chance to think about this :-)

Cheers,
Jamie

On Tue, Oct 4, 2016 at 1:16 PM, Daniel Wheeler 
wrote:

> Hi James,
>
> Thanks for this example. I think I understand what is happening. Let's
> start with the code for solving the problem,
>
> 
> import fipy
> import numpy
>
> mesh = fipy.Grid1D(nx=10)
>
> var = fipy.CellVariable(mesh=mesh)
> var[:] = 2 * numpy.random.random(mesh.numberOfCells)
> var.constrain(1., where=mesh.facesRight)
> var.faceGrad.constrain(0, where=mesh.facesLeft)
>
> eqn = fipy.ConvectionTerm([-1]) == fipy.DiffusionTerm(1.0)
>
> for _ in range(100):
> eqn.sweep(var)
> print var
> ~~~
>
> Note that the "var.faceGrad.constrain" is required with convection
> terms to make sure that there is an outlet flow. This is a weird
> decision in many ways, but was made for backwards compatibility with
> previous FiPy versions, but this is not the main issue here.
>
> Notice that when the diffusion coefficient is small or large then this
> system solves in one sweep or close to one sweep (set it to 1000.0 or
> 0.0001 for example). The number of sweeps required to solve the system
> is at its worst when the diffiusion and convection coefficients are
> close in magnitude. So what is going on? Basically, the system is not
> fully implicit.
>
> The left hand boundary condition is dependent on the left hand cell
> value but gets added to the b vector in the case of non upwind
> convection terms. In the case of fully upwind convection terms this
> contribution is zero so as the diffusion coefficient goes to zero the
> CovectionTerm becomes fully upwind (or if a UpwindConvectionTerm is
> selected). As the diffusion coefficient gets large the contribution of
> the convection term decreases and thus the contribution from the
> explicit boundary condition is reduced.
>
> The discretization for the left most cell for a central difference
> scheme (we're using power law above, but it has some central
> difference contribution at D=1 and u=-1) is
>
>   D * phi0 - (u / 2 + D) * phi1 = -u * phiB / 2
>
> where phi0 is the value at the left most cell and phi1 is the cell to
> the right. In FiPy, phiB (the boundary value of phiB) is not
> calculated implicitly, but is explicitly set for the convection term
> boundary condition. Basically phiB is the previous value of phi0 in
> the FiPy code. That may not be the best way to deal with this. It may
> be that it needs to be explicit for stability reasons.
>
> A possible way to deal with this issue is to use terms for the outlet
> condition that are upwind at the outflow, but power law everywhere
> else.
>
> Anyway, I hope that helps with understanding what's going on even if
> the implementation is sub-optimal.
>
> Cheers,
>
> Daniel
>
> On Thu, Sep 15, 2016 at 5:13 PM, James Pringle  wrote:
> > Dear FiPy users --
> >
> >This is a simple example of how and why fipy may fail to solve a
> >steady advection diffusion problem, and how solving the transient
> >problem can reduce the error. I also found something that was a
> >surprise -- the "initial" condition of a steady problem can affect
> >the solution for some solvers.
> >
> >At the end are two interesting questions for those who want to
> >understand what FiPy is actually doing I admit to being a bit
> >lost
> >
> >The equation I am solving is
> >
> > \Del\dot (D\del psi + u*psi) =0
> >
> >Where the diffusion D is 1, and u is a vector (1,0) -- so there is
> >only a flow of speed -1 in the x direction.  I solve this equation
> >on a 10 by 10 grid. The boundary conditions are no normal gradient
> >on the y=0 and y=10 boundary:
> >
> > psi_y =0 at y=0 and y=10
> >
> >For the x boundary, I impose a value of x=1 on the inflow boundary at
> > x=10
> >(this is a little tricky -- the way the equation is written, u is
> >the negative of velocity).
> >
> > psi=1 at x=10
> >
> >and a no-normal-gradient condition at x=0.
> >
> > psi_x=0 at x=0
> >
> >since all of the domain and boundary is symmetrical with respect to
> >y, we can assume psi_y=0 is zero everywhere. This reduces the
> equation to
> >
> > psi_xx + psi_x =0
> >
> >The general solution to this equation is
> >
> > psi=C1+C2*exp(-x)
> >
> >Where C1 and C2 are constants. For these boundary conditions, C1=1
> >and C2=0, so psi=1 everywhere.
> >
> &g

Re: sources errors in advection/diffusion problems, and one solution

[James Pringle]

Interesting. One addition to Zhekai's comment. When I read it, I thought
Upwind might be doing better because it is adding additional numerical
diffusion (as upwind usually does, at least in finite difference codes).
However this is not the case. If you keep the advection u and diffusion D
the 1D equation is

 D*psi_xx + u*psi_x =0

and the exponential term of the solution (the one that appears in the
error) has the form exp((u/D)*x).  An increase in diffusion D would make
the error decay more slowly away from the boundary (you can test that this
is so by increasing "DiffCoef" in the code). This increase in decay scale
is not what appears when UpwindConvectionTerm is used for the convection
term. Instead, the solution becomes perfect (Note -- in this context
advection and convection are the same thing.)

But it does give a hint as to the source of the error. The upwind schemes only
use the gradient calculated in the direction the advection the current is
coming from -- in this case u is from the positive x direction. This means
an advection scheme in a finite-difference code (which is all that I know,
I am new to finite element codes) would not usually use the gradient
calculate on the x=0 boundary since the current is coming from +x).

But the upwindConvection scheme is sensitive to the boundary condition on
the x=0 boundary. (change the line
"psi.faceGrad.constrain(0.0*faceNorm,where=xeq0bound)" to see this.) And it
appears to be doing the right thing when the gradient is set to some other
value than 0.

So how does the upwind scheme implement the normal gradient boundary
condition differently from the other schemes? At this point I should get
off my lazy behind and dig through the code, but I have a number of
deadlines for the rest of the week as I switch from being an incompetent
applied mathematician to an incompetent writer and administrator.

But if no one else has an idea, I might look next week.

Jamie


On Tue, Sep 20, 2016 at 6:52 PM, Zhekai Deng <
zhekaideng2...@u.northwestern.edu> wrote:

> Hi James,
>
> Thanks for providing this demo to illustrate the problem. I don't have any
> particular ideas exactly why the initial value of psi helps to reduce the
> error and why transient problem  "advect" out. However, I have some
> findings that may help on this.
>
> Finding # 1: I noticed you have tried ExponentialConvectionTerm,
> PowerLawConvectionTerm, CentralDifferenceConvectionTerm, all of these
> give exponential error. However, if you tried UpwindConvectionTerm, this
> will give the right result on the steady state solution. Thus, maybe using
> upwind method, the convection does not require the value from downstream,
> consequently, the error from the downstream will not excite the error
> toward the upstream. However, I am still very surprise with the magnitude
> of the error from other methods, and how similar
>
> Finding # 2: In the transient state problem, if I increase the mesh size
> from 50*50 to 100*100. The error actually grows larger for the
> ExponentialConvectionTerm, PowerLawConvectionTerm,
> CentralDifferenceConvectionTerm. To show this, if I change the mesh size
> to 100*100, the max psi value I have is around 1.1 . However, if I change
> the mesh size to 50*50, the max psi is 1.005366, which is several orders of
> magnitude lower in terms of difference to the exact solution. This is also
> the case for the UpwindConvectionTerm, however, the error for both mesh
> size are very small (max(psi) = 1e-13 or 1e-14  + 1). So even in the
> transient state, the mesh size appears to somehow amplify the error if we
> use finer mesh. I am confused by this.
>
> To now, it seems UpWindConvectionTerm appears to the the work around
> method to this issue. Of course, I will be willing to learn more about what
> other people think on this.
>
> Best,
>
> Zhekai
>
>
>
>
> On Fri, Sep 16, 2016 at 8:53 AM, James Pringle  wrote:
>
>> No worries -- I had to do it to figure out the problem in my more complex
>> domain and equation... The issue which surprised me was that the value the
>> variable was initialized to had an effect on the steady solution.
>>
>> Jamie
>>
>> On Fri, Sep 16, 2016 at 8:14 AM, Guyer, Jonathan E. Dr. (Fed) <
>> jonathan.gu...@nist.gov> wrote:
>>
>>> James -
>>>
>>> I think Daniel will have more insight into why this is happening and if
>>> there is anything to be done about it besides artificial relaxation, but I
>>> just want to say how much I appreciate your putting this together. This is
>>> a very lucid illustration.
>>>
>>> - Jon
>>>
>>> > On Sep 15, 2016, at 5:13 PM, James Pringle  wrote:
>>> >
>>> > Dear FiPy users

Re: sources errors in advection/diffusion problems, and one solution

[James Pringle]

No worries -- I had to do it to figure out the problem in my more complex
domain and equation... The issue which surprised me was that the value the
variable was initialized to had an effect on the steady solution.

Jamie

On Fri, Sep 16, 2016 at 8:14 AM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> James -
>
> I think Daniel will have more insight into why this is happening and if
> there is anything to be done about it besides artificial relaxation, but I
> just want to say how much I appreciate your putting this together. This is
> a very lucid illustration.
>
> - Jon
>
> > On Sep 15, 2016, at 5:13 PM, James Pringle  wrote:
> >
> > Dear FiPy users --
> >
> >This is a simple example of how and why fipy may fail to solve a
> >steady advection diffusion problem, and how solving the transient
> >problem can reduce the error. I also found something that was a
> >surprise -- the "initial" condition of a steady problem can affect
> >the solution for some solvers.
> >
> >At the end are two interesting questions for those who want to
> >understand what FiPy is actually doing I admit to being a bit
> >lost
> >
> >The equation I am solving is
> >
> > \Del\dot (D\del psi + u*psi) =0
> >
> >Where the diffusion D is 1, and u is a vector (1,0) -- so there is
> >only a flow of speed -1 in the x direction.  I solve this equation
> >on a 10 by 10 grid. The boundary conditions are no normal gradient
> >on the y=0 and y=10 boundary:
> >
> > psi_y =0 at y=0 and y=10
> >
> >For the x boundary, I impose a value of x=1 on the inflow boundary at
> x=10
> >(this is a little tricky -- the way the equation is written, u is
> >the negative of velocity).
> >
> > psi=1 at x=10
> >
> >and a no-normal-gradient condition at x=0.
> >
> > psi_x=0 at x=0
> >
> >since all of the domain and boundary is symmetrical with respect to
> >y, we can assume psi_y=0 is zero everywhere. This reduces the
> equation to
> >
> > psi_xx + psi_x =0
> >
> >The general solution to this equation is
> >
> > psi=C1+C2*exp(-x)
> >
> >Where C1 and C2 are constants. For these boundary conditions, C1=1
> >and C2=0, so psi=1 everywhere.
> >
> >Now run the code SquareGrid_HomemadeDelaunay and look at figure(3)
> >-- this is the plot of psi versus x, and you can see that it does
> >not match the true solution of psi=1 everywhere! Instead, it
> >appears to be decaying exponential. The blue line is actually just
> >(1+exp(-x)). What is going on? It appears that small errors in the
> >boundary condition at x=0 are allowing C2 to not be exactly 0, and
> >this error is this exponential mode. The error is the artificial
> >exiting of a correct mode of the interior equation, albeit one that
> >should not be excited by these BC's.
> >
> >Interestingly, the size of this error depends on the value the psi
> >is initially set to. If the line
> >
> >psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0)
> >
> >is changed so psi is initially 1, the error goes away entirely; if
> >it is set to some other value, you get different errors. I do not
> >know if this is a bug, or just numerical error in a well programmed
> >solver.
> >
> >Now if you run SquareGrid_HomemadeDelaunay_transient  which
> implements
> >
> >  psi_t = \Del\dot (D\del psi + u*psi)
> >
> >you can see that the error in the numerical solution is advected
> >out of the x=0 boundary, and the solution approaches the true
> >solution of psi=1 rapidly.
> >
> >The interesting question is if the formulation of the boundary
> >condition at x=0 could be altered to less excite the spurious mode?
> >
> >Also, why does the "initial" condition have any effect on the
> >steady equation?
> >
> >Cheers,
> >Jamie
> >
> >  transient.py>__
> _
> > fipy mailing list
> > fipy@nist.gov
> > https://urldefense.proofpoint.com/v2/url?u=http-3A__www.
> ctcms.nist.gov_fipy&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=
> 7HJI3EH6RqKWf16fbYIYKw&m=IaTL4nEetwjm4G4qfj8cwLzvztKzui
> Fsw_Nhksv_oWQ&s=vG-nxTf76KxE_CqEHTjt2jIkoy0l9M6X8bm01ypXaBQ&e=
> >  [ NIST internal ONLY: https://urldefens

sources errors in advection/diffusion problems, and one solution

[James Pringle]

Dear FiPy users --

   This is a simple example of how and why fipy may fail to solve a
   steady advection diffusion problem, and how solving the transient
   problem can reduce the error. I also found something that was a
   surprise -- the "initial" condition of a steady problem can affect
   the solution for some solvers.

   At the end are two interesting questions for those who want to
   understand what FiPy is actually doing I admit to being a bit
   lost

   The equation I am solving is

\Del\dot (D\del psi + u*psi) =0

   Where the diffusion D is 1, and u is a vector (1,0) -- so there is
   only a flow of speed -1 in the x direction.  I solve this equation
   on a 10 by 10 grid. The boundary conditions are no normal gradient
   on the y=0 and y=10 boundary:

psi_y =0 at y=0 and y=10

   For the x boundary, I impose a value of x=1 on the inflow boundary at
x=10
   (this is a little tricky -- the way the equation is written, u is
   the negative of velocity).

psi=1 at x=10

   and a no-normal-gradient condition at x=0.

psi_x=0 at x=0

   since all of the domain and boundary is symmetrical with respect to
   y, we can assume psi_y=0 is zero everywhere. This reduces the equation to

psi_xx + psi_x =0

   The general solution to this equation is

psi=C1+C2*exp(-x)

   Where C1 and C2 are constants. For these boundary conditions, C1=1
   and C2=0, so psi=1 everywhere.

   Now run the code SquareGrid_HomemadeDelaunay and look at figure(3)
   -- this is the plot of psi versus x, and you can see that it does
   not match the true solution of psi=1 everywhere! Instead, it
   appears to be decaying exponential. The blue line is actually just
   (1+exp(-x)). What is going on? It appears that small errors in the
   boundary condition at x=0 are allowing C2 to not be exactly 0, and
   this error is this exponential mode. The error is the artificial
   exiting of a correct mode of the interior equation, albeit one that
   should not be excited by these BC's.

   Interestingly, the size of this error depends on the value the psi
   is initially set to. If the line

   psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0)

   is changed so psi is initially 1, the error goes away entirely; if
   it is set to some other value, you get different errors. I do not
   know if this is a bug, or just numerical error in a well programmed
   solver.

   Now if you run SquareGrid_HomemadeDelaunay_transient  which implements

 psi_t = \Del\dot (D\del psi + u*psi)

   you can see that the error in the numerical solution is advected
   out of the x=0 boundary, and the solution approaches the true
   solution of psi=1 rapidly.

   The interesting question is if the formulation of the boundary
   condition at x=0 could be altered to less excite the spurious mode?

   Also, why does the "initial" condition have any effect on the
   steady equation?

   Cheers,
   Jamie
from pylab import *
from numpy import *
import fipy
from collections import defaultdict
import copy as Copy

def ScipyDelaunay2mesh(points,simplices,debug=False):
'''ScipyDelaunay2mesh(points,simplices):

Creates a FiPy 2Dmesh object from the output of a
scipy.spatial.Delaunay triangulation, 

Also returns a list of border faces for the application of
boundary conditions. The list of border conditions is computed
within this routine; it will handle the case where
triangles/simplices have been removed from the output of the
Delaunay() triangulation.


Parameters
--
points : array

   an (N,2) array of the location of N points that made the
   triangulation. If Delaunay is called with
   tri=scipy.spatial.Delaunay(), points is tri.points

simplices : array

   an (M,3) array of the location of M triangles returned by
   Delaunay. If Delaunay is called with
   tri=scipy.spatial.Delaunay(), simplices is tri.simplices

debug=False : Boolean

   if this is set to true, the mesh will be drawn in blue,
   boundaries will be drawn in red, boundary point numbers will be
   written at vertices, and numbers for triangles will be drawn in
   the triangle. This output is useful for finding where triangles
   have been removed incorrectly. 

Returns
---

 mesh: a FiPy mesh2D object

'''

#Create the arrays that fipy.mesh.mesh2D needs.  vertexCoords is
#of shape (2,N), where N is the number of points, and contains the
#coordinates of the vertices. It is equivalent to points in
#content, but reshaped.
print 'Entering ScipyDelaunay2Mesh; Making vertexCoords'
vertexCoords=points.T

#faceVertexID is of shape (2,M) where M is the number of faces/edges
#(faces is the fipy terminology, edges is the scipy.spatial.Delaunay
#terminology). faceVertexID contains the points (as indices into
#vertexCoords) which make up each face. This data can be extracted
#from s

Re: Solving advection, segregation and diffusion equation

[James Pringle]

Dear all -- in a bit I am going to post a little idealized
advection/diffusion problem to give an idea of how some of these errors in
steady solutions arise, and how iterating in time helps. It will be much
more basic then this work, but you might find it helpful. It was very
helpful for me in figuring out my more complex system to reproduce the
error in a simple way. It also reveals some surprising issues in the fipy
code -- for example, the error in a steady solution can depend on the
"initial" condition you feed the solver.

It will be under the title "sources errors in advection/diffusion problems,
and one solution." I hope to post it tonight.

Jamie

On Thu, Sep 15, 2016 at 1:27 PM, Zhekai Deng <
zhekaideng2...@u.northwestern.edu> wrote:

> Thanks for the reply. Those are very useful suggestions. I have some new
> findings which may help shedding some light on this problem.
>
> I have added TransientTerm and tried to see how the solution is developed
> over time. I also remove the datamin and datamax input on the viewer, so
> that I can see the full range of the solution variable data.
>
> When I see the temporal development of the solution (attached in the
> email), It looks stable. However, I noticed that the right outlet boundary
> has extreme high concentration and did not outflow the flow domain. I think
> It has due to the improper setup of my right outlet condition. Since the
> materials keep building up on the right outlet boundary, it appears to
> affect the solution in the actual flow domain(also attached in email) and
> drives the solution unstable. This could be the main reason that solving
> steady state equation directly never converges.
>
> I then begin to look into how to properly set up "open outlet" boundary,
> the closest I could find is following:
>
> https://www.mail-archive.com/fipy@nist.gov/msg02797.html
> 
> and in the section of "Applying fixed flux boundary conditions" of
> official website
> http://www.ctcms.nist.gov/fipy/documentation/USAGE.html
> 
>
> I tried to adapt the official method listed on the official website, but
> no success is achieved.  Thus, I wonder is it possible to point out some
> examples that used this method to define fixed flux outlet boundary
> condition ? You may find my current progress on this part in the code I
> attached.
>
> Another thing I noticed is that,  there is a reported issue on robin
> boundary condition:
>
> https://github.com/usnistgov/fipy/issues/426
> 
>
> Would it be possible to explain more about this ? I don't think the up &
> bottom boundary condition in my case causes the non convergent problem,
> because from the Transient equation, the segregation of the material starts
> to develop at the top&bottom boundary, which is consistent with my
> expectation.
>
> I have attached my current progress of the code, and the temporal
> development of the solution.
>
> Best,
>
> Zhekai
>
>
>
>
>
>
> On Thu, Sep 15, 2016 at 9:33 AM, Daniel Wheeler  > wrote:
>
>> On Wed, Sep 14, 2016 at 7:52 PM, Guyer, Jonathan E. Dr. (Fed)
>>  wrote:
>> >
>> > I don't know offhand. With a Péclet number of 100, your problem is
>> almost completely hyperbolic, which FiPy (and cell-centered Finite Volume)
>> isn't very good at. Daniel knows more about this and may have some
>> suggestions.
>> >
>> > You might consider adding a TransientTerm for artificial relaxation and
>> trying the VanLeerConvectionTerm.
>>
>> Definitely use a TransientTerm and see the time step you can get away
>> with. Don't sweep, just use time steps. Using the VanLeerConvetionTerm
>> isn't necessary as the accuracy is not important in a steady state
>> problem.
>>
>> --
>> Daniel Wheeler
>>
>> ___
>> 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 convection term

[James Pringle]

Have not looked at your code. But a quick note -- Burger's equation is
prone to forming shocks, whose width is set by the diffusive terms. In
numerical solutions, the solution then is often determined by numerical
diffusivity, which is non-optimal. Play around with analytical solutions
with shocks and make sure your code is doing something sensible. It is a
hard one.

Jamie

On Mon, Aug 29, 2016 at 9:09 AM, Sebastian Glane  wrote:

> Dear fipy-users and developers,
>
> I am new to fipy and would like to solve Burgers' equation as a first
> step. This problem involves a non-linear convection term. According to
> the FAQs, for the diffusion term it is to program the nonlinearity using
> coeff. I did the same for the convection term:
> fipy.VanLeerConvectionTerm( coeff = (0.5 * phi, ), var = phi ). But
> somehow it does not work out. I also tried with phi.faceValue.
>
> Best regards
>
> Sebastian
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: error in gradient even with orthogonal grid

[James Pringle]

Thank you -- since my problem solves for a stream function, this was a very
helpful tip, and has made my solutions appear much better.

Jamie

On Mon, Aug 1, 2016 at 4:49 PM, Daniel Wheeler 
wrote:

> Hi Jamie,
>
> Sorry for the delayed reply.
>
> On Thu, Jul 21, 2016 at 11:31 AM, James Pringle  wrote:
> >
> > However, I am still getting large errors in estimates of the gradient
> even
> > when the mesh is nearly perfectly orthogonal. These errors DO NOT become
> > smaller as resolution is increased. I have two question:
> >
> > Are these errors to be expected in unstructured meshes?
>
> Cell centered FV has two issues with unstructured meshes. You know
> about the orthogonality issue, but there is a second issue. The
> docstring from the following code might help explain,
>
>
> https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_usnistgov_fipy_blob_develop_fipy_variables_cellVariable.py-23L257&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=F8ThJ6NxKtDlYIWkDUqGq2ql_J0shUHQMI1lQCB1Wvo&s=RM15FSI_eegFR5VfTb9il1PZ2B1_yaQ6eKmZHq-8Hx0&e=
>
> See, `\frac{1}{V_P} \sum_f \vec{n} \phi_f A_f`. This is the
> discretized gradient calculation. In this formulation using the
> triangular mesh, the calculated \phi_f is not the correct average for
> the face as the line between the cell centers doesn't cross the face
> in the middle of the face. I think this is called non-conjunctionality
> in finite volume speak. See equation 9.8 in the following link for
> further explanation (we should probably implement this approach).
>
>
> https://urldefense.proofpoint.com/v2/url?u=http-3A__bit.ly_2aqRiHt&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=F8ThJ6NxKtDlYIWkDUqGq2ql_J0shUHQMI1lQCB1Wvo&s=EJ34X0r43RF50EgJ2p_clCS-H5yY8nSTitKGqtbINts&e=
>
> I believe that this is the source of the error in the gradient
> calculation. It's a consistent error independent of the mesh
> refinement I think (or maybe first order). An alternative gradient
> calculation is the least squares gradient, see
>
>
> https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_usnistgov_fipy_blob_develop_fipy_variables_cellVariable.py-23L270&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=F8ThJ6NxKtDlYIWkDUqGq2ql_J0shUHQMI1lQCB1Wvo&s=l6F4maMJnB3bfM6--FXXLzrUKyaH1i7KHuo1atLYJqM&e=
>
> Using that
>
> 
> import fipy as fp
> nx = 25
> dx = 1. / nx
> mesh = fp.Tri2D(nx=nx, ny=nx, dx=dx, dy=dx)
> psi = fp.CellVariable(mesh=mesh)
> psi.constrain(1.0, where=mesh.facesLeft)
> psi.faceGrad.constrain(1.0, where=mesh.facesRight)
> eq = fp.DiffusionTerm().solve(psi)
> print psi.leastSquaresGrad[:, :100]
> 
>
> seems to give lots of 1s and 0s which seems correct. The cells next to
> the boundary are not correct as they aren't using the boundary
> conditions in the calculation. I'm not sure about the order of
> accuracy of the least squares or why it happens to work for this mesh.
>
> > Is it acceptable to average the gradient over space to reduce the error?
>
> I don't think so unless it can be shown to improve the accuracy with
> analysis.
>
> --
> Daniel Wheeler
> ___
> fipy mailing list
> fipy@nist.gov
>
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=F8ThJ6NxKtDlYIWkDUqGq2ql_J0shUHQMI1lQCB1Wvo&s=A7D0acjE011uAk-Vl-5OhrTNVdH9kgmq0rshSN6tbsg&e=
>   [ NIST internal ONLY:
> https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=DQICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=F8ThJ6NxKtDlYIWkDUqGq2ql_J0shUHQMI1lQCB1Wvo&s=V8y_Fr-VcNq0fTtFnfiyabj8MXL8mmGa5LFCMv0Sx_Q&e=
> ]
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: scipy's Delaunay output to fipy mesh2D object

[James Pringle]

Dear all --

A warning about using the above code. Since the choice of triangles
formed by a Delaunay triangulation can be somewhat arbitrary in the case of
a somewhat regular distribution of points, it will choose some
triangulations that lead to severe orthogonality errors in fipy, even if
that was not necessary for the points it was given.

   So use with caution and monitor non-orthogonality with the
mesh._nonOrthogonal property. For my mostly regular grid I am resorting to
triangulation-by-hand.

Cheers,
Jamie Pringle
University of New Hampshire

On Fri, Jun 17, 2016 at 10:14 AM, James Pringle  wrote:

> Thanks --
>
>Once I have pounded on this a bit more, I will write the function. It
> should be straightforward. I concentrated on making it simple and easy to
> understand; since it is run only to make the grid, efficiency was not a
> great issue for me.
>
> Jamie
>
> On Fri, Jun 17, 2016 at 10:03 AM, Daniel Wheeler <
> daniel.wheel...@gmail.com> wrote:
>
>> James, this is awesome, thanks for posting. It would be a very good
>> idea to have a DelaunayMesh in FiPy that used Scipy for the
>> triangulation as it would give another way to make and test
>> unstructured meshes without the Gmsh dependency. Something like,
>>
>> mesh = DelaunayMesh(points_in_2D)
>>
>> In the long run, it would be good to change your code so that it
>> doesn'y loop in Python. I suspect that your code is quite inefficient
>> due to these loops. However, you have something that works, which is
>> great and it is only a one time cost at the beginning of the
>> simulation. Is efficiency an issue for you with this?
>>
>> On Thu, Jun 16, 2016 at 4:29 PM, James Pringle  wrote:
>> > For the use of others --
>> >
>> >  Sometimes when dealing with arbitrary grids, it is useful to use
>> the
>> > output of routines other than Gmsh to create the fipy mesh. In my case,
>> I
>> > need to create a complex grid from ocean bathymetry, and I ended up
>> using
>> > scipy.spatial.Delaunay to do so. With help from Daniel and Jonathan, I
>> > created a code to create a mesh2D object from the output of Delaunay.
>> It is
>> > written to emphasize clarity over speed. I hope others find it useful.
>> >
>> > James Pringle
>> >
>> > #==
>> > #This code creates a fipy grid from the output of a Delaunay
>> > #triangulation. The code is written to prioratize clarity over speed.
>> >
>> > from pylab import *
>> > from numpy import *
>> > import fipy
>> > from scipy.spatial import Delaunay
>> >
>> >
>> > #make a simple example set of points to triangulate with delaunay
>> > points=array([[ 0.,  0.],
>> >[ 0.,  1.],
>> >[ 1.,  0.],
>> >[ 1.,  1.]])
>> > tri=Delaunay(points)
>> >
>> > #when this code is run as is, it should produce
>> > #  faceVertexIDs=
>> > #   [[3, 1, 0, 2, 0
>> > #[1, 0, 3, 3, 2]]
>> > #and
>> > #  cellFaceIds=
>> > #[[0, 3],
>> > # [1, 2],
>> > # [2, 4]]
>> >
>> >
>> > #now create the arrays that fipy.mesh.mesh2D needs.  vertexCoords is
>> > #of shape (2,N), where N is the number of points, and contains the
>> > #coordinates of the vertices. It is equivalent to tri.points in
>> > #content, but reshaped.
>> > print 'Making vertexCoords'
>> > vertexCoords=tri.points.T
>> >
>> > #faceVertexID is of shape (2,M) where M is the number of faces/edges
>> > #(faces is the fipy terminology, edges is the scipy.spatial.Delaunay
>> > #terminology). faceVertexID contains the points (as indices into
>> > #vertexCoords) which make up each face. This data can be extracted
>> > #from simplices from Delaunay, which contains the faces of each
>> > #triangle. HOWEVER, simplices contains many repeated faces, since each
>> > #triangle will border on others. Only one copy of each face should be
>> > #inserted into faceVertexID.
>> >
>> > print 'Making faceVertexIDs'
>> > faceIn={} #this is a dictionary of all faces that have already been
>> >   #inserted into the faceVertexID, with the key a tuple of
>> >   #indices, sorted
>> > faceVertexIDs_list=[] #structure to build
>> >
>> > for ntri in xrange(tri.simplices.shape[0]):
>> > for nface in xrange(3):
>> > if nface==0:
>> > 

Re: faceGradAverage broken or I am confused?

[James Pringle]

Sorry for the multiple emails, but I spoke in haste. The psi.grad and
psi.faceGrad appear to have systemic errors that are larger than I
anticipated. The attached code shows the error in the solution (small, and
it decreases as the grid spacing decreases), and the error in gradient as
calculated by both psi.faceGrad and psi.grad. For both gradient
calculations, the the error of the gradient does not appear to decrease as
resolution becomes finer, and remains about 40% in both.

Resolution is controlled by the line "xvec=linspace(0.0,1.0,50)"

As before, the code is straightforward if you skip to the "if __main__"
section.  The analytical solution is psi = x.

   1. Figure 1 shows the grid and solution
   2. Figure 2 shows the solution minus the analytical solution
   3. Figure 3 shows the gradient from psi.faceGrad
   4. Figure 4 shows the gradient from psi.grad


Cheers,
Jamie Pringle

On Mon, Jul 18, 2016 at 4:42 PM, Pringle, James 
wrote:

> Ok, Then there is also a subtle issue with psi.faceGrad or something else.
> As above, I am solving the problem
>
> eq=(fipy.DiffusionTerm(var=psi,coeff=1.0))
>
>
> With the BC's
>
> \Del\Psi \dot \nhat = -1 at x=0
>
> \Del\Psi \dot \nhat =0 at y=0 and y=1
>
> Psi =1 at x=1
>
>
> With the code
>
> #y=0 boundary
> yeq0bound=yFace==0.0
> psi.faceGrad.constrain(0.0*faceNorm,where=yeq0bound)
>
> #y=1.0 boundary
> yeq1bound=yFace==1.0
> psi.faceGrad.constrain(0.0*faceNorm,where=yeq1bound)
>
> #x=0 boundary
> xeq0bound=xFace==0.0
> psi.faceGrad.constrain(-1.0*faceNorm,where=xeq0bound)
>
> #x=1 boundary
> xeq1bound=xFace==1.0
> psi.constrain(0.0,where=xeq1bound)
>
> The solution on the cell variables matches the analytical solution of
> psi=x, but when I get the gradient with psi.faceGrad it is incorrect on the
> y=0 and y=1 boundary (it is correct everywhere else, including the x
> boundaries). On the y boundaries, psi.faceGrad returns (0,0) at each point,
> instead of the correct (1,0).  This is not such a big deal, but should be
> documented.
>
> Code illustrating the error is attached.
>
> Thanks,
> Jamie
>
>
>
>
>
> On Mon, Jul 18, 2016 at 3:17 PM, Guyer, Jonathan E. Dr. (Fed) <
> jonathan.gu...@nist.gov> wrote:
>
>> They return \nabla\psi, the gradient, not the projection onto the normal
>> (or anything else).
>>
>> > On Jul 18, 2016, at 1:47 PM, James Pringle  wrote:
>> >
>> > Dear Jonathan --
>> >
>> >I have to admit I picked faceGradAverage over faceGrad randomly.
>> One comment on my test code -- it shuffles the order of the points that
>> Delaunay uses to make the mesh each run.  I did this as a debugging step of
>> my Delaunay code.  If the issue is the order of the triangle points, you
>> should get different answers each time, and indeed you do.   Just search
>> for "shuffle" in my test code.
>> >
>> >Also, the documentation for faceGrad and faceGradAverage (
>> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy_fipy_generated_fipy.variables.html-23fipy.variables.cellVariable.CellVariable.faceGrad&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=gtWqG3GBoMw65_GwtEbeSFgCLT0l_wMGfV7MCrPo89k&s=vHBVDYOH9hZ7yTR_To6QH5i-5bO4QbHKmXpEShgMJ54&e=
>> ) is unclear. Do these routines return \Del\psi, or \Del\psi\dot\nhat?
>> I.e. do they return the gradient or the gradient projected onto the unit
>> normal to the face?
>> >
>> >If it is the former, there is also a problem in faceGrad.
>> >
>> > Thanks!
>> > Jamie
>> >
>> > On Mon, Jul 18, 2016 at 1:23 PM, Guyer, Jonathan E. Dr. (Fed) <
>> jonathan.gu...@nist.gov> wrote:
>> > James -
>> >
>> > I didn't even know we had a .faceGradAverage!
>> >
>> > The issue turns out to be with .grad (.faceGradAverage is calculated
>> from .grad and .grad is calculated from .faceValue, whereas .faceGrad is
>> calculated from .value (at cell centers)). Neither .grad nor
>> .faceGradAverage is used in calculating the solution to the Laplace
>> problem, which is why the solution looks OK.
>> >
>> > My guess is that some of the triangles you send to Mesh2D are
>> "upside-down". We don't stipulate an ordering and that upside-downess
>> should impact the calculation of the diffusion stencil, so our best guess
>> is that we correct for the vertex-face-cell ordering elsewhere and just
>> need to apply the same correction to .grad. We'll try to sort this out in
>> the next few days.
>&g

Re: faceGradAverage broken or I am confused?

[James Pringle]

Ok, Then there is also a subtle issue with psi.faceGrad or something else.
As above, I am solving the problem

eq=(fipy.DiffusionTerm(var=psi,coeff=1.0))


With the BC's

\Del\Psi \dot \nhat = -1 at x=0

\Del\Psi \dot \nhat =0 at y=0 and y=1

Psi =1 at x=1


With the code

#y=0 boundary
yeq0bound=yFace==0.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq0bound)

#y=1.0 boundary
yeq1bound=yFace==1.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq1bound)

#x=0 boundary
xeq0bound=xFace==0.0
psi.faceGrad.constrain(-1.0*faceNorm,where=xeq0bound)

#x=1 boundary
xeq1bound=xFace==1.0
psi.constrain(0.0,where=xeq1bound)

The solution on the cell variables matches the analytical solution of
psi=x, but when I get the gradient with psi.faceGrad it is incorrect on the
y=0 and y=1 boundary (it is correct everywhere else, including the x
boundaries). On the y boundaries, psi.faceGrad returns (0,0) at each point,
instead of the correct (1,0).  This is not such a big deal, but should be
documented.

Code illustrating the error is attached.

Thanks,
Jamie





On Mon, Jul 18, 2016 at 3:17 PM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> They return \nabla\psi, the gradient, not the projection onto the normal
> (or anything else).
>
> > On Jul 18, 2016, at 1:47 PM, James Pringle  wrote:
> >
> > Dear Jonathan --
> >
> >I have to admit I picked faceGradAverage over faceGrad randomly.  One
> comment on my test code -- it shuffles the order of the points that
> Delaunay uses to make the mesh each run.  I did this as a debugging step of
> my Delaunay code.  If the issue is the order of the triangle points, you
> should get different answers each time, and indeed you do.   Just search
> for "shuffle" in my test code.
> >
> >Also, the documentation for faceGrad and faceGradAverage (
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy_fipy_generated_fipy.variables.html-23fipy.variables.cellVariable.CellVariable.faceGrad&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=gtWqG3GBoMw65_GwtEbeSFgCLT0l_wMGfV7MCrPo89k&s=vHBVDYOH9hZ7yTR_To6QH5i-5bO4QbHKmXpEShgMJ54&e=
> ) is unclear. Do these routines return \Del\psi, or \Del\psi\dot\nhat?
> I.e. do they return the gradient or the gradient projected onto the unit
> normal to the face?
> >
> >If it is the former, there is also a problem in faceGrad.
> >
> > Thanks!
> > Jamie
> >
> > On Mon, Jul 18, 2016 at 1:23 PM, Guyer, Jonathan E. Dr. (Fed) <
> jonathan.gu...@nist.gov> wrote:
> > James -
> >
> > I didn't even know we had a .faceGradAverage!
> >
> > The issue turns out to be with .grad (.faceGradAverage is calculated
> from .grad and .grad is calculated from .faceValue, whereas .faceGrad is
> calculated from .value (at cell centers)). Neither .grad nor
> .faceGradAverage is used in calculating the solution to the Laplace
> problem, which is why the solution looks OK.
> >
> > My guess is that some of the triangles you send to Mesh2D are
> "upside-down". We don't stipulate an ordering and that upside-downess
> should impact the calculation of the diffusion stencil, so our best guess
> is that we correct for the vertex-face-cell ordering elsewhere and just
> need to apply the same correction to .grad. We'll try to sort this out in
> the next few days.
> >
> > - Jon
> >
> > > On Jul 15, 2016, at 3:20 PM, James Pringle  wrote:
> > >
> > > I am trying to debug some boundary conditions in a complex problem. I
> have run into a problem with faceGradAverage or my understanding of it. It
> is illustrated in a simple problem
> > >
> > >  \Del^2 Psi =0 on a unit square
> > >  normal derivative = 1 on x=0
> > >  normal derivative = 0 on y=0
> > >  normal derivative = 0 on y=1
> > >  Psi =1 on x=1
> > >
> > > The analytical solution to this problem is
> > >
> > >  Psi = x
> > >
> > > The code is in the attached file; skip until "if __name__ ==
> "__main__":" , the earlier code converts a Delaunay triangulation to a fipy
> mesh. I am doing this to match my more complex case.
> > >
> > > Anyway, when the code is run, I get the analytical solution, and so
> gradient \Del\Psi should be (1,0) everywhere. It is not. Note also that
> since \Del\Psi is spatially uniform, how the gradient is averaged in space
> should not make a difference.
> > >
> > > How am I confused?
> > >
> > > So that you don't have to open the attachment, the code that solves
> th

Re: faceGradAverage broken or I am confused?

[James Pringle]

Dear Jonathan --

   I have to admit I picked faceGradAverage over faceGrad randomly.  One
comment on my test code -- it shuffles the order of the points that
Delaunay uses to make the mesh each run.  I did this as a debugging step of
my Delaunay code.  If the issue is the order of the triangle points, you
should get different answers each time, and indeed you do.   Just search
for "shuffle" in my test code.

   Also, the documentation for faceGrad and faceGradAverage (
http://www.ctcms.nist.gov/fipy/fipy/generated/fipy.variables.html#fipy.variables.cellVariable.CellVariable.faceGrad)
is unclear. Do these routines return \Del\psi, or \Del\psi\dot\nhat?  I.e.
do they return the gradient or the gradient projected onto the unit normal
to the face?

   If it is the former, there is also a problem in faceGrad.

Thanks!
Jamie

On Mon, Jul 18, 2016 at 1:23 PM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> James -
>
> I didn't even know we had a .faceGradAverage!
>
> The issue turns out to be with .grad (.faceGradAverage is calculated from
> .grad and .grad is calculated from .faceValue, whereas .faceGrad is
> calculated from .value (at cell centers)). Neither .grad nor
> .faceGradAverage is used in calculating the solution to the Laplace
> problem, which is why the solution looks OK.
>
> My guess is that some of the triangles you send to Mesh2D are
> "upside-down". We don't stipulate an ordering and that upside-downess
> should impact the calculation of the diffusion stencil, so our best guess
> is that we correct for the vertex-face-cell ordering elsewhere and just
> need to apply the same correction to .grad. We'll try to sort this out in
> the next few days.
>
> - Jon
>
> > On Jul 15, 2016, at 3:20 PM, James Pringle  wrote:
> >
> > I am trying to debug some boundary conditions in a complex problem. I
> have run into a problem with faceGradAverage or my understanding of it. It
> is illustrated in a simple problem
> >
> >  \Del^2 Psi =0 on a unit square
> >  normal derivative = 1 on x=0
> >  normal derivative = 0 on y=0
> >  normal derivative = 0 on y=1
> >  Psi =1 on x=1
> >
> > The analytical solution to this problem is
> >
> >  Psi = x
> >
> > The code is in the attached file; skip until "if __name__ ==
> "__main__":" , the earlier code converts a Delaunay triangulation to a fipy
> mesh. I am doing this to match my more complex case.
> >
> > Anyway, when the code is run, I get the analytical solution, and so
> gradient \Del\Psi should be (1,0) everywhere. It is not. Note also that
> since \Del\Psi is spatially uniform, how the gradient is averaged in space
> should not make a difference.
> >
> > How am I confused?
> >
> > So that you don't have to open the attachment, the code that solves the
> problem and prints the face gradients is included as text here:
> >
> > #setup fipy solution for Del^2 Psi = 1 with BC of 0 everywhere
> > #first, make grid variable
> > psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0)
> >
> > #apply boundary of Del\Psi \dot \nhat=-1 on x=0 boundary (equivalent
> to d Psi/dx=1)
> > #  Del\Psi \dot \nhat=0 on y=0,1 boundary
> > #and Psi=1 on x=1
> >
> > xFace,yFace=mesh.faceCenters
> > faceNorm=mesh.faceNormals
> >
> > #y=0 boundary
> > yeq0bound=yFace==0.0
> > psi.faceGrad.constrain(0.0*faceNorm,where=yeq0bound)
> >
> > #y=1.0 boundary
> > yeq1bound=yFace==1.0
> > psi.faceGrad.constrain(0.0*faceNorm,where=yeq1bound)
> >
> > #x=0 boundary
> > xeq0bound=xFace==0.0
> > psi.faceGrad.constrain(-1.0*faceNorm,where=xeq0bound)
> >
> > #x=1 boundary
> > xeq1bound=xFace==1.0
> > psi.constrain(0.0,where=xeq1bound)
> >
> >
> > #make equatiion and solve
> > eq=(fipy.DiffusionTerm(var=psi,coeff=1.0))
> > eq.solve(var=psi)
> >
> > #now plot solution
> > subplot(2,1,2)
> > tripcolor(psi.mesh.vertexCoords[0,:],psi.mesh.vertexCoords[1,:],
> >   tri.simplices,psi.numericValue)
> > colorbar()
> > title('solution')
> >
> > #now print faceGradients
> > print 'Face Gradients on x=0, should be
> 1,0',psi.faceGradAverage[:,xeq0bound]
> > print ' '
> > print 'Face Gradients on y=0, should be
> 1,0',psi.faceGradAverage[:,yeq0bound]
> >
> > Thanks,
> > Jamie
> > _

faceGradAverage broken or I am confused?

[James Pringle]

I am trying to debug some boundary conditions in a complex problem. I have
run into a problem with faceGradAverage or my understanding of it. It is
illustrated in a simple problem

 \Del^2 Psi =0 on a unit square
 normal derivative = 1 on x=0
 normal derivative = 0 on y=0
 normal derivative = 0 on y=1
 Psi =1 on x=1

The analytical solution to this problem is

 Psi = x

The code is in the attached file; skip until "if __name__ == "__main__":" ,
the earlier code converts a Delaunay triangulation to a fipy mesh. I am
doing this to match my more complex case.

Anyway, when the code is run, I get the analytical solution, and so
gradient \Del\Psi should be (1,0) everywhere. It is not. Note also that
since \Del\Psi is spatially uniform, how the gradient is averaged in space
should not make a difference.

How am I confused?

So that you don't have to open the attachment, the code that solves the
problem and prints the face gradients is included as text here:

#setup fipy solution for Del^2 Psi = 1 with BC of 0 everywhere
#first, make grid variable
psi=fipy.CellVariable(name='psi',mesh=mesh,value=0.0)

#apply boundary of Del\Psi \dot \nhat=-1 on x=0 boundary (equivalent to
d Psi/dx=1)
#  Del\Psi \dot \nhat=0 on y=0,1 boundary
#and Psi=1 on x=1

xFace,yFace=mesh.faceCenters
faceNorm=mesh.faceNormals

#y=0 boundary
yeq0bound=yFace==0.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq0bound)

#y=1.0 boundary
yeq1bound=yFace==1.0
psi.faceGrad.constrain(0.0*faceNorm,where=yeq1bound)

#x=0 boundary
xeq0bound=xFace==0.0
psi.faceGrad.constrain(-1.0*faceNorm,where=xeq0bound)

#x=1 boundary
xeq1bound=xFace==1.0
psi.constrain(0.0,where=xeq1bound)


#make equatiion and solve
eq=(fipy.DiffusionTerm(var=psi,coeff=1.0))
eq.solve(var=psi)

#now plot solution
subplot(2,1,2)
tripcolor(psi.mesh.vertexCoords[0,:],psi.mesh.vertexCoords[1,:],
  tri.simplices,psi.numericValue)
colorbar()
title('solution')

#now print faceGradients
print 'Face Gradients on x=0, should be
1,0',psi.faceGradAverage[:,xeq0bound]
print ' '
print 'Face Gradients on y=0, should be
1,0',psi.faceGradAverage[:,yeq0bound]

Thanks,
Jamie
# These modules takes the output of scipy.spatial.Delaunay's
# triangulation of data points, and returns both the FiPy mesh for the
# points, and a list of the boundary edges so that boundary conditions
# can be set for the points.
#
# This code handles the case where triangles/simplices have been
# removed from output of the Delaunay triangulation after the
# triangulation.  This triangle removal can be useful to create
# interior holes in the domain, or to handle concave boundaries of the
# domain.
#
# The documentation for the routines can be found in there docstrings,
# and examples of there use can be found in __main__

from numpy import *
import fipy
from collections import defaultdict
import copy as Copy

def ScipyDelaunay2mesh(points,simplices,debug=False):
'''ScipyDelaunay2mesh(points,simplices):

Creates a FiPy 2Dmesh object from the output of a
scipy.spatial.Delaunay triangulation, 

Also returns a list of border faces for the application of
boundary conditions. The list of border conditions is computed
within this routine; it will handle the case where
triangles/simplices have been removed from the output of the
Delaunay() triangulation.


Parameters
--
points : array

   an (N,2) array of the location of N points that made the
   triangulation. If Delaunay is called with
   tri=scipy.spatial.Delaunay(), points is tri.points

simplices : array

   an (M,3) array of the location of M triangles returned by
   Delaunay. If Delaunay is called with
   tri=scipy.spatial.Delaunay(), simplices is tri.simplices

debug=False : Boolean

   if this is set to true, the mesh will be drawn in blue,
   boundaries will be drawn in red, boundary point numbers will be
   written at vertices, and numbers for triangles will be drawn in
   the triangle. This output is useful for finding where triangles
   have been removed incorrectly. 

Returns
---

 mesh: a FiPy mesh2D object

 borderList: a list of lists
  
   Each list in borderList is a list of the indices of the faces
   as indexed in mesh.cellFaceIDs. Each list is for a single
   contiguous border of the domain. These will be used to set the
   boundary conditions on the mesh.

'''

#Create the arrays that fipy.mesh.mesh2D needs.  vertexCoords is
#of shape (2,N), where N is the number of points, and contains the
#coordinates of the vertices. It is equivalent to points in
#content, but reshaped.
print 'Entering ScipyDelaunay2Mesh; Making vertexCoords'
vertexCoords=points.T

#faceVertexID is of shape (2,M) where M is the number of faces/edge

Re: Gradient boundary conditions for spherical shell

[James Pringle]

Please ignore -- I just realized I am being a dolt. For those of you who
stumble across this thread and would like avoid making the same mistake, I
believe that

 psi.faceGrad.constrain(1.0*mesh.faceNormals,faceMask)

constrains the faceGrad to some value, here 1.0*mesh.faceNormals, where
faceMask is true. If I want to set it the gradient to some value which is
normal in my coordinate system, I just have to provide a (2,numberOfFaces)
array with the appropriate values of the derivatives of psi in it. There is
nothing magical about the mesh.faceNormals other than that it gives you
useful information about the orientation of the face in your variable
space.

Sorry to waste everyones time.
Jamie

On Fri, Jul 1, 2016 at 10:35 AM, Pringle, James 
wrote:

> Dear fipy folks --
>
> I am solving a problem of ocean circulation on a thin spherical shell; the
> domain is too large to ignore the spherical nature of the earth. The
> problem is of the form
>
>eq = (DiffusionTerm(var=Psi,coeff=DiffCoeff)+
> ExponentialConvectionTerm(var=Psi,coeff=convCoeff))
>
> And I have incorporated all the of the appropriate metric terms for a
> spherical shell into the coefficients above. My problem comes as I strive
> to make the solution more realistic by applying Neumann BC's to some parts
> of the problem. The usual formulation for setting \Del\dot\Psi \dot normal
> = 1
>
>psi.faceGrad.constrain(1.0*mesh.faceNormals,faceMask)
>
> will not work because on a spherical shell the gradient is equal to
>
>\Del f = (1/r) dF/d(lat)*(unit vector latitude) +
> (1/(r*cos(lat)))*dF/d(long)*(unit vector longitude)
>
> Is there any easy way to constrain the faceGrad to something like
>
> constant = A*dPsi/dx+B*dPsi/dy
>
> where x and y are the spatial coordinates, and A and B are arbitrary
> factors?  That way I could easily include the metrics...
>
> Thanks,
> Jamie Pringle
> University of New Hampshire
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: Gradient boundary conditions for spherical shell

[James Pringle]

Just to be clear -- this problem is on a shell -- there is no radial
dependence, and the mesh is two-dimensional.

Jamie

On Fri, Jul 1, 2016 at 10:35 AM, Pringle, James 
wrote:

> Dear fipy folks --
>
> I am solving a problem of ocean circulation on a thin spherical shell; the
> domain is too large to ignore the spherical nature of the earth. The
> problem is of the form
>
>eq = (DiffusionTerm(var=Psi,coeff=DiffCoeff)+
> ExponentialConvectionTerm(var=Psi,coeff=convCoeff))
>
> And I have incorporated all the of the appropriate metric terms for a
> spherical shell into the coefficients above. My problem comes as I strive
> to make the solution more realistic by applying Neumann BC's to some parts
> of the problem. The usual formulation for setting \Del\dot\Psi \dot normal
> = 1
>
>psi.faceGrad.constrain(1.0*mesh.faceNormals,faceMask)
>
> will not work because on a spherical shell the gradient is equal to
>
>\Del f = (1/r) dF/d(lat)*(unit vector latitude) +
> (1/(r*cos(lat)))*dF/d(long)*(unit vector longitude)
>
> Is there any easy way to constrain the faceGrad to something like
>
> constant = A*dPsi/dx+B*dPsi/dy
>
> where x and y are the spatial coordinates, and A and B are arbitrary
> factors?  That way I could easily include the metrics...
>
> Thanks,
> Jamie Pringle
> University of New Hampshire
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Gradient boundary conditions for spherical shell

[James Pringle]

Dear fipy folks --

I am solving a problem of ocean circulation on a thin spherical shell; the
domain is too large to ignore the spherical nature of the earth. The
problem is of the form

   eq = (DiffusionTerm(var=Psi,coeff=DiffCoeff)+
ExponentialConvectionTerm(var=Psi,coeff=convCoeff))

And I have incorporated all the of the appropriate metric terms for a
spherical shell into the coefficients above. My problem comes as I strive
to make the solution more realistic by applying Neumann BC's to some parts
of the problem. The usual formulation for setting \Del\dot\Psi \dot normal
= 1

   psi.faceGrad.constrain(1.0*mesh.faceNormals,faceMask)

will not work because on a spherical shell the gradient is equal to

   \Del f = (1/r) dF/d(lat)*(unit vector latitude) +
(1/(r*cos(lat)))*dF/d(long)*(unit vector longitude)

Is there any easy way to constrain the faceGrad to something like

constant = A*dPsi/dx+B*dPsi/dy

where x and y are the spatial coordinates, and A and B are arbitrary
factors?  That way I could easily include the metrics...

Thanks,
Jamie Pringle
University of New Hampshire
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: determine cell neighbors.

[James Pringle]

Thanks -- unfortunately, I have a complex real world mesh with interior
points (islands) which must be treated in an iterative fashion, so I need
to identify them in the grid.

On Fri, Jul 1, 2016 at 9:30 AM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> You want `mesh._cellToCellIDs`.
>
> This isn't a basic question at all. We don't encourage thinking about the
> discretized mesh and don't make it easy. _cellToCellIDs isn't documented
> anywhere as far as I can tell.
>
> > On Jun 29, 2016, at 4:07 PM, James Pringle  wrote:
> >
> > I have to find a group of contagious cells which match a certain
> criteria, which involves values on the cell center. I can figure out which
> cell centers meet the criteria.
> >
> > However, from the mesh2D variable, how do I figure out which cells are
> contiguous to a given cell? Must I build a data structure from
> variable.mesh.cellFaceIDs? Is not that connectivity stored somewhere?
> >
> > I am sorry for what must be a very basic question, but I have been
> searching for a while.
> >
> > Thanks,
> > Jamie Pringle
> > ___
> > fipy mailing list
> > fipy@nist.gov
> >
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=j21xGXOI7t5RRCHtFFIovyygu3Lofr9ph6qEsGEmLc8&s=uWPzbXrhm_AJViYOrZmt-fspuUiMLR9EhlUF-7pFpsU&e=
> >  [ NIST internal ONLY:
> https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=j21xGXOI7t5RRCHtFFIovyygu3Lofr9ph6qEsGEmLc8&s=Q_2FGcnCUXd2d1diXnIZ4vFcrzw91qJ3JQOptGd_lWQ&e=
> ]
>
>
> ___
> fipy mailing list
> fipy@nist.gov
>
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=j21xGXOI7t5RRCHtFFIovyygu3Lofr9ph6qEsGEmLc8&s=uWPzbXrhm_AJViYOrZmt-fspuUiMLR9EhlUF-7pFpsU&e=
>   [ NIST internal ONLY:
> https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=j21xGXOI7t5RRCHtFFIovyygu3Lofr9ph6qEsGEmLc8&s=Q_2FGcnCUXd2d1diXnIZ4vFcrzw91qJ3JQOptGd_lWQ&e=
> ]
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

determine cell neighbors.

[James Pringle]

I have to find a group of contagious cells which match a certain criteria,
which involves values on the cell center. I can figure out which cell
centers meet the criteria.

However, from the mesh2D variable, how do I figure out which cells are
contiguous to a given cell? Must I build a data structure from
variable.mesh.cellFaceIDs?
Is not that connectivity stored somewhere?

I am sorry for what must be a very basic question, but I have been
searching for a while.

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

Re: bug in linearBicgstabSolver.py

[James Pringle]

no worries, thanks

On Wed, Jun 22, 2016 at 10:58 AM, Daniel Wheeler 
wrote:

> Jamie, the new code is on the develop branch now
>
>
> https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_usnistgov_fipy_blob_develop_fipy_solvers_scipy_linearBicgstabSolver.py&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=FoPBLFKX5Z7TW3Sk-72AFrPn3EJESyjcGiyOfgsy7_s&s=Biiu3V9QVjqpBgmvVqdcLIoTV88lYbiG1-anl-4mUSI&e=
>
> It seemed to work for a small test case. Unfortunately, that solver
> isn't part of the test suite.
>
> On Wed, Jun 22, 2016 at 10:42 AM, James Pringle  wrote:
> > thanks
> > Jamie
>
> --
> Daniel Wheeler
> ___
> fipy mailing list
> fipy@nist.gov
>
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=FoPBLFKX5Z7TW3Sk-72AFrPn3EJESyjcGiyOfgsy7_s&s=J6QRk5zhR5NwcHCD_cwJ_TtPPK1BnT-89TSdrw28XPo&e=
>   [ NIST internal ONLY:
> https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=FoPBLFKX5Z7TW3Sk-72AFrPn3EJESyjcGiyOfgsy7_s&s=7RgqyeMB4ocaYTTrqqckSEauTlHtvJOG8aO9T9DBBHI&e=
> ]
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: bug in linearBicgstabSolver.py

[James Pringle]

thanks
Jamie

On Wed, Jun 22, 2016 at 10:04 AM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> That's terrible. Not to name any names, but Daniel did that about five
> years ago.
>
> I've submitted a fix at
> https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_usnistgov_fipy_pull_497&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=rqSLJ5XWBtgVE-7mbJie2VpGSqiUU-ARkbMDtlhb-Qk&s=ekHWBVJYr8pw7evCYx2IbG68LgM2t0bzf9V4FO08Psw&e=
>
> > On Jun 21, 2016, at 9:22 PM, James Pringle  wrote:
> >
> > The document string for fipy.solvers.scipy.linearBicgstabSolver.py
> claims it takes the arguments "tolerance," "iterations" and "precon", but
> it does not.
> >
> > To reproduce this bug, make sure you are using the scipy solvers, and
> compare the results for, e.g.
>  eq.solve(var=psi,solver=fipy.LinearPCGSolver(tolerance=1.0e-3)) to
> eq.solve(var=psi,solver=fipy.LinearBicgstabSolver(tolerance=1.0e-3)). The
> later throws the error:
> >
> > In [43]:
> eq.solve(var=psi,solver=fipy.LinearBicgstabSolver(tolerance=1.0e-3))
> >
> ---
> > TypeError Traceback (most recent call
> last)
> >  in ()
> > > 1
> eq.solve(var=psi,solver=fipy.LinearBicgstabSolver(tolerance=1.0e-3))
> >
> > TypeError: __init__() got an unexpected keyword argument 'tolerance'
> >
> >
> > In fipy version 3.1; the error appears to be in the latest git version
> of the code as well; it appears the problem is in the
> fipy/fipy/solvers/scipy/linearBicgstabSolver.py, where the arguments are
> simply not based to bicgstab.
> >
> > This is too bad, for LinearBicgstabSolver gives by far the best results
> on my problem, and the underlying scipy.sparse.linalg.bicgstab does take
> the argument "tol".
> >
> > I think I could fix the problem staring at the code, but the
> class/inheritance that is going on is above what I really understand, and I
> am afraid I would bolix things further.
> >
> > Thanks,
> > Jamie
> > ___
> > fipy mailing list
> > fipy@nist.gov
> >
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=rqSLJ5XWBtgVE-7mbJie2VpGSqiUU-ARkbMDtlhb-Qk&s=53Zt7rb_fGwvi-ZERkInTL4oEOiILcAkz4LrCqREm_o&e=
> >  [ NIST internal ONLY:
> https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=rqSLJ5XWBtgVE-7mbJie2VpGSqiUU-ARkbMDtlhb-Qk&s=xYOeFF5MIcAEsgelMexvlh_smIT-Bb36TV2ypQPDllg&e=
> ]
>
>
> ___
> fipy mailing list
> fipy@nist.gov
>
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=rqSLJ5XWBtgVE-7mbJie2VpGSqiUU-ARkbMDtlhb-Qk&s=53Zt7rb_fGwvi-ZERkInTL4oEOiILcAkz4LrCqREm_o&e=
>   [ NIST internal ONLY:
> https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=7HJI3EH6RqKWf16fbYIYKw&m=rqSLJ5XWBtgVE-7mbJie2VpGSqiUU-ARkbMDtlhb-Qk&s=xYOeFF5MIcAEsgelMexvlh_smIT-Bb36TV2ypQPDllg&e=
> ]
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

bug in linearBicgstabSolver.py

[James Pringle]

The document string for fipy.solvers.scipy.linearBicgstabSolver.py claims
it takes the arguments "tolerance," "iterations" and "precon", but it does
not.

To reproduce this bug, make sure you are using the scipy solvers, and
compare the results for, e.g.
  eq.solve(var=psi,solver=fipy.LinearPCGSolver(tolerance=1.0e-3)) to
eq.solve(var=psi,solver=fipy.LinearBicgstabSolver(tolerance=1.0e-3)). The
later throws the error:

In [43]:
eq.solve(var=psi,solver=fipy.LinearBicgstabSolver(tolerance=1.0e-3))
---
TypeError Traceback (most recent call last)
 in ()
> 1 eq.solve(var=psi,solver=fipy.LinearBicgstabSolver(tolerance=1.0e-3))

TypeError: __init__() got an unexpected keyword argument 'tolerance'



In fipy version 3.1; the error appears to be in the latest git version of
the code as well; it appears the problem is in the
fipy/fipy/solvers/scipy/linearBicgstabSolver.py, where the arguments are
simply not based to bicgstab.

This is too bad, for LinearBicgstabSolver gives by far the best results on
my problem, and the underlying scipy.sparse.linalg.bicgstab does take the
argument "tol".

I think I could fix the problem staring at the code, but the
class/inheritance that is going on is above what I really understand, and I
am afraid I would bolix things further.

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

appropriate solver/Convection term selection

[James Pringle]

Dear FiPy folks --

   I am solving a 2D elliptic equation (from a steady-state diffusion
   system) with strong advection/convection terms
   (e.g. hyperbolic-ish). The convection terms are often locally
   dominant, however the over all solution is determined by the
   diffusion/elliptic terms.

   One way to think about it is that the streamlines of the convection
   terms form closed loops in some parts of the domain, and so the
   solution in those loops is set by the diffusion across the closed
   stream lines. For those who care, this problem is similar to
   Stommel's solution for gyre scale circulation in the ocean, but
   with bathymetry.

   The diffusivity is slightly anisotropic (factor of no more than
   0.5). The velocity-like terms are incompressible -- they can be
   expressed as a stream function. I suspect the solution to this
   steady problem will not be overly sensitive to
   upwind-spurious-diffusion. The fipy formulation is:

   eq=(fipy.DiffusionTerm(var=psi,coeff=DiffCoeff)+
   fipy.ExponentialConvectionTerm(var=psi,coeff=convCoeff))

   The default scipy solver is going in the right direction, but the
   solution has not fully converged in the interior of the closed
   streamlines (I test this by setting all boundaries to psi=1, so the
   solution must be psi=1 everywhere). If I increase diffusion enough,
   the scipy solver does solve it correctly.

   The default pySparse solver dies after running out of memmory on a
   machine with lots of it. I am configuring trilinos now. In the
   distant past I had success with similar problems with the MUDPACK
   class of multigrid methods (using mud2).

   What solvers would you suggest to work with a problem like this?
   What guides you in your choice of solvers? What guides your choice
   of the type of the Convection term?

Thanks,
Jamie Pringle
University of New Hampshire
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: scipy's Delaunay output to fipy mesh2D object

[James Pringle]

Thanks --

   Once I have pounded on this a bit more, I will write the function. It
should be straightforward. I concentrated on making it simple and easy to
understand; since it is run only to make the grid, efficiency was not a
great issue for me.

Jamie

On Fri, Jun 17, 2016 at 10:03 AM, Daniel Wheeler 
wrote:

> James, this is awesome, thanks for posting. It would be a very good
> idea to have a DelaunayMesh in FiPy that used Scipy for the
> triangulation as it would give another way to make and test
> unstructured meshes without the Gmsh dependency. Something like,
>
> mesh = DelaunayMesh(points_in_2D)
>
> In the long run, it would be good to change your code so that it
> doesn'y loop in Python. I suspect that your code is quite inefficient
> due to these loops. However, you have something that works, which is
> great and it is only a one time cost at the beginning of the
> simulation. Is efficiency an issue for you with this?
>
> On Thu, Jun 16, 2016 at 4:29 PM, James Pringle  wrote:
> > For the use of others --
> >
> >  Sometimes when dealing with arbitrary grids, it is useful to use the
> > output of routines other than Gmsh to create the fipy mesh. In my case, I
> > need to create a complex grid from ocean bathymetry, and I ended up using
> > scipy.spatial.Delaunay to do so. With help from Daniel and Jonathan, I
> > created a code to create a mesh2D object from the output of Delaunay. It
> is
> > written to emphasize clarity over speed. I hope others find it useful.
> >
> > James Pringle
> >
> > #==
> > #This code creates a fipy grid from the output of a Delaunay
> > #triangulation. The code is written to prioratize clarity over speed.
> >
> > from pylab import *
> > from numpy import *
> > import fipy
> > from scipy.spatial import Delaunay
> >
> >
> > #make a simple example set of points to triangulate with delaunay
> > points=array([[ 0.,  0.],
> >[ 0.,  1.],
> >[ 1.,  0.],
> >[ 1.,  1.]])
> > tri=Delaunay(points)
> >
> > #when this code is run as is, it should produce
> > #  faceVertexIDs=
> > #   [[3, 1, 0, 2, 0
> > #[1, 0, 3, 3, 2]]
> > #and
> > #  cellFaceIds=
> > #[[0, 3],
> > # [1, 2],
> > # [2, 4]]
> >
> >
> > #now create the arrays that fipy.mesh.mesh2D needs.  vertexCoords is
> > #of shape (2,N), where N is the number of points, and contains the
> > #coordinates of the vertices. It is equivalent to tri.points in
> > #content, but reshaped.
> > print 'Making vertexCoords'
> > vertexCoords=tri.points.T
> >
> > #faceVertexID is of shape (2,M) where M is the number of faces/edges
> > #(faces is the fipy terminology, edges is the scipy.spatial.Delaunay
> > #terminology). faceVertexID contains the points (as indices into
> > #vertexCoords) which make up each face. This data can be extracted
> > #from simplices from Delaunay, which contains the faces of each
> > #triangle. HOWEVER, simplices contains many repeated faces, since each
> > #triangle will border on others. Only one copy of each face should be
> > #inserted into faceVertexID.
> >
> > print 'Making faceVertexIDs'
> > faceIn={} #this is a dictionary of all faces that have already been
> >   #inserted into the faceVertexID, with the key a tuple of
> >   #indices, sorted
> > faceVertexIDs_list=[] #structure to build
> >
> > for ntri in xrange(tri.simplices.shape[0]):
> > for nface in xrange(3):
> > if nface==0:
> > face=tri.simplices[ntri,[0,1]]
> > elif nface==1:
> > face=tri.simplices[ntri,[1,2]]
> > else:
> > face=tri.simplices[ntri,[2,0]]
> > faceKey=tuple(sort(face))
> > if not (faceKey in faceIn):
> > faceIn[faceKey]=True
> > faceVertexIDs_list.append(face)
> >
> > faceVertexIDs=array(faceVertexIDs_list).T
> >
> > #now create cellFaceIDs of shape (3,P) where P is the number of cells
> > #in the domain. It contains the faces (as indices into faceVertexIDs)
> > #that make up each cell.
> >
> > #first create dictionary with a key made up of a sorted tuple of vertex
> > #indices which maps from a face to its location in faceVertexIDs
> > print 'Making cellFaceIDs'
> > faceMap={}
> > for nface in xrange(faceVertexIDs.shape[1]):
> > faceKey=tuple(sort(faceVertexIDs[:,nface]))
> > faceMap[faceKey]=nface
> >
> > #

scipy's Delaunay output to fipy mesh2D object

[James Pringle]

For the use of others --

 Sometimes when dealing with arbitrary grids, it is useful to use the
output of routines other than Gmsh to create the fipy mesh. In my case, I
need to create a complex grid from ocean bathymetry, and I ended up using
scipy.spatial.Delaunay to do so. With help from Daniel and Jonathan, I
created a code to create a mesh2D object from the output of Delaunay. It is
written to emphasize clarity over speed. I hope others find it useful.

James Pringle

#==
#This code creates a fipy grid from the output of a Delaunay
#triangulation. The code is written to prioratize clarity over speed.

from pylab import *
from numpy import *
import fipy
from scipy.spatial import Delaunay


#make a simple example set of points to triangulate with delaunay
points=array([[ 0.,  0.],
   [ 0.,  1.],
   [ 1.,  0.],
   [ 1.,  1.]])
tri=Delaunay(points)

#when this code is run as is, it should produce
#  faceVertexIDs=
#   [[3, 1, 0, 2, 0
#[1, 0, 3, 3, 2]]
#and
#  cellFaceIds=
#[[0, 3],
# [1, 2],
# [2, 4]]


#now create the arrays that fipy.mesh.mesh2D needs.  vertexCoords is
#of shape (2,N), where N is the number of points, and contains the
#coordinates of the vertices. It is equivalent to tri.points in
#content, but reshaped.
print 'Making vertexCoords'
vertexCoords=tri.points.T

#faceVertexID is of shape (2,M) where M is the number of faces/edges
#(faces is the fipy terminology, edges is the scipy.spatial.Delaunay
#terminology). faceVertexID contains the points (as indices into
#vertexCoords) which make up each face. This data can be extracted
#from simplices from Delaunay, which contains the faces of each
#triangle. HOWEVER, simplices contains many repeated faces, since each
#triangle will border on others. Only one copy of each face should be
#inserted into faceVertexID.

print 'Making faceVertexIDs'
faceIn={} #this is a dictionary of all faces that have already been
  #inserted into the faceVertexID, with the key a tuple of
  #indices, sorted
faceVertexIDs_list=[] #structure to build

for ntri in xrange(tri.simplices.shape[0]):
for nface in xrange(3):
if nface==0:
face=tri.simplices[ntri,[0,1]]
elif nface==1:
face=tri.simplices[ntri,[1,2]]
else:
face=tri.simplices[ntri,[2,0]]
faceKey=tuple(sort(face))
if not (faceKey in faceIn):
faceIn[faceKey]=True
faceVertexIDs_list.append(face)

faceVertexIDs=array(faceVertexIDs_list).T

#now create cellFaceIDs of shape (3,P) where P is the number of cells
#in the domain. It contains the faces (as indices into faceVertexIDs)
#that make up each cell.

#first create dictionary with a key made up of a sorted tuple of vertex
#indices which maps from a face to its location in faceVertexIDs
print 'Making cellFaceIDs'
faceMap={}
for nface in xrange(faceVertexIDs.shape[1]):
faceKey=tuple(sort(faceVertexIDs[:,nface]))
faceMap[faceKey]=nface

#now loop over simplices, find each edge/face, and map from points to
#faces
cellFaceIDs=tri.simplices.T*0
for ntri in xrange(tri.simplices.shape[0]):
for nface in xrange(3):
if nface==0:
face=tri.simplices[ntri,[0,1]]
elif nface==1:
face=tri.simplices[ntri,[1,2]]
else:
face=tri.simplices[ntri,[2,0]]
faceKey=tuple(sort(face))
cellFaceIDs[nface,ntri]=faceMap[faceKey]

#ok, now instantiate a mesh2D object with this data
print 'Making mesh'
mesh=fipy.meshes.mesh2D.Mesh2D(vertexCoords,faceVertexIDs,cellFaceIDs)
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: how to create large grid with holes and odd geometry

[James Pringle]

I have successfully written and tested code to convert the output of
spicy.spatial.Delaunay's triangulation into a fipy mesh. I have posted it
to this mailing list under a separate message with the title "scipy's
Delaunay output to fipy mesh2D object"

Jamie Pringle
University of New Hampshire
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: how to create large grid with holes and odd geometry

[James Pringle]

Jonathan and Daniel -- Thanks. In a few days or so, I shall post what I
come up with for others to use.
Jamie

On Wed, Jun 15, 2016 at 3:58 PM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> You can generate a simple mesh in order to probe what's stored by doing
> something like:
>
>   mesh1 = fp.Grid2D(nx=2, ny=2)
>   mesh2 = fp.Grid2D(nx=2, ny=2) + [[2], [1]]
>   mesh = mesh1 + mesh2
>
>
> > On Jun 15, 2016, at 2:06 PM, James Pringle  wrote:
> >
> > Jonathan --
> >
> > Thank you, this is exactly what I need. Two more questions.
> >   • Is the order of the vertices in faceVertexIDs important?
> >   • Is the order of faces in cellFaceIDs important? Must they wind
> clockwise or counterclockwise?
> > Thanks,
> > Jamie
> >
> > On Wed, Jun 15, 2016 at 1:56 PM, Guyer, Jonathan E. Dr. (Fed) <
> jonathan.gu...@nist.gov> wrote:
> > The first three arguments to the Mesh2D constructor are (all that is)
> required:
> >
> >   class Mesh2D(Mesh):
> >   def __init__(self, vertexCoords, faceVertexIDs, cellFaceIDs, ...):
> >
> > All other arguments have default values assigned.
> >
> > For your case:
> >
> > vertexCoords is of shape (2, N) where N is the number of vertices
> > faceVertexIDs is of shape (2, M) where M is the number of faces
> > cellFaceIDs is of shape (3, P) where P is the number of cells
> >
> > faceVertexIDs and cellFaceIDs are 0-based, as they are indices into the
> preceding array
> >
> > > On Jun 15, 2016, at 1:29 PM, James Pringle  wrote:
> > >
> > > Well, I am motivated to give it a go, since I only have the summer to
> make progress on project and it is blocking my research progress right now.
> Can you give me a pointer to where the appropriate quantities are defined?
> I can certainly write code to make the transformations, but it is a bit
> hard without understanding precisely what is defined in the Mesh2D object.
> I have made a simple Mesh2D object, but I am not sure which of the
> attributes, etc, are absolutely required, or how they are precisely defined.
> > >
> > > Perhaps most useful for me would be
> > >   • a definition of which parts of the Mesh2D object must exist
> > >   • and the format of those parts, in particular the a face to
> vertex array and a
> > > cell to face array
> > > Or you could point me to the appropriate part of the Gmsh code so I
> can go from there. I presume poking around in fipy.Gmsh2D would be a good
> place to start? Is there a better place to start?
> > >
> > > I would love any documentation on the details of the Mesh2D object.
> > >
> > > Thanks,
> > > Jamie Pringle
> > > University of New Hampshire
> > >
> > > On Wed, Jun 15, 2016 at 12:21 PM, Daniel Wheeler <
> daniel.wheel...@gmail.com> wrote:
> > > Hi Jamie,
> > >
> > > There is no automated way to make a FiPy mesh from Scipy's Delaunay
> > > object. The problem is that FiPy needs a face to vertex array and a
> > > cell to face array while the Delaunay object just has the cell to
> > > vertex array. The functionality to extract the face to vertex array
> > > and cell to face array is in FiPy because it must be done for Gmsh,
> > > however, it is not abstracted in so that it can be reused.
> > >
> > > It is certainly possible to make the correct arrays with Numpy and
> > > pass them to the Mesh2D object, but it's a bit of work to write the
> > > code. If I find some time I might give it a go, but I don't know when
> > > I will get to it.
> > >
> > > Cheers,
> > >
> > > Daniel
> > >
> > > On Tue, Jun 14, 2016 at 4:15 PM, James Pringle 
> wrote:
> > > > Daniel et al. --
> > > >
> > > > As referenced earlier in this thread, I have a complex domain
> with holes
> > > > in it; I now have it broken up into triangles, based on the Delaunay
> package
> > > > in SciPy. I have
> > > >
> > > > Locations of all vertex coordinates,
> > > > list of vertices that make up faces
> > > > list of faces that make cells
> > > > list of faces that make up all the internal and external boundaries.
> > > >
> > > > Can you point me to any code or documentation that would help me
> understand
> > > > how to make this into a Mesh2D object? I am having a devil of a time
> > > > figuring out from the manual online. The best would be s

Re: how to create large grid with holes and odd geometry

[James Pringle]

Jonathan --

Thank you, this is exactly what I need. Two more questions.

   1. Is the order of the vertices in faceVertexIDs important?
   2. Is the order of faces in cellFaceIDs important? Must they wind
   clockwise or counterclockwise?

Thanks,
Jamie

On Wed, Jun 15, 2016 at 1:56 PM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> The first three arguments to the Mesh2D constructor are (all that is)
> required:
>
>   class Mesh2D(Mesh):
>   def __init__(self, vertexCoords, faceVertexIDs, cellFaceIDs, ...):
>
> All other arguments have default values assigned.
>
> For your case:
>
> vertexCoords is of shape (2, N) where N is the number of vertices
> faceVertexIDs is of shape (2, M) where M is the number of faces
> cellFaceIDs is of shape (3, P) where P is the number of cells
>
> faceVertexIDs and cellFaceIDs are 0-based, as they are indices into the
> preceding array
>
> > On Jun 15, 2016, at 1:29 PM, James Pringle  wrote:
> >
> > Well, I am motivated to give it a go, since I only have the summer to
> make progress on project and it is blocking my research progress right now.
> Can you give me a pointer to where the appropriate quantities are defined?
> I can certainly write code to make the transformations, but it is a bit
> hard without understanding precisely what is defined in the Mesh2D object.
> I have made a simple Mesh2D object, but I am not sure which of the
> attributes, etc, are absolutely required, or how they are precisely defined.
> >
> > Perhaps most useful for me would be
> >   • a definition of which parts of the Mesh2D object must exist
> >   • and the format of those parts, in particular the a face to
> vertex array and a
> > cell to face array
> > Or you could point me to the appropriate part of the Gmsh code so I can
> go from there. I presume poking around in fipy.Gmsh2D would be a good place
> to start? Is there a better place to start?
> >
> > I would love any documentation on the details of the Mesh2D object.
> >
> > Thanks,
> > Jamie Pringle
> > University of New Hampshire
> >
> > On Wed, Jun 15, 2016 at 12:21 PM, Daniel Wheeler <
> daniel.wheel...@gmail.com> wrote:
> > Hi Jamie,
> >
> > There is no automated way to make a FiPy mesh from Scipy's Delaunay
> > object. The problem is that FiPy needs a face to vertex array and a
> > cell to face array while the Delaunay object just has the cell to
> > vertex array. The functionality to extract the face to vertex array
> > and cell to face array is in FiPy because it must be done for Gmsh,
> > however, it is not abstracted in so that it can be reused.
> >
> > It is certainly possible to make the correct arrays with Numpy and
> > pass them to the Mesh2D object, but it's a bit of work to write the
> > code. If I find some time I might give it a go, but I don't know when
> > I will get to it.
> >
> > Cheers,
> >
> > Daniel
> >
> > On Tue, Jun 14, 2016 at 4:15 PM, James Pringle  wrote:
> > > Daniel et al. --
> > >
> > > As referenced earlier in this thread, I have a complex domain with
> holes
> > > in it; I now have it broken up into triangles, based on the Delaunay
> package
> > > in SciPy. I have
> > >
> > > Locations of all vertex coordinates,
> > > list of vertices that make up faces
> > > list of faces that make cells
> > > list of faces that make up all the internal and external boundaries.
> > >
> > > Can you point me to any code or documentation that would help me
> understand
> > > how to make this into a Mesh2D object? I am having a devil of a time
> > > figuring out from the manual online. The best would be something that
> used
> > > the output of either the triangles or scipy.spatial.Delaunay()
> packaged.  My
> > > equation is of the form
> > >
> > > 0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi)
> > >
> > >
> > > and I can get the coefficients A(x,y) and B(x,y) on either the faces
> or in
> > > the cell centers are needed.
> > >
> > > Thank you.
> > > Jamie Pringle
> > > University of New Hampshire
> >
> > --
> > Daniel Wheeler
> > ___
> > fipy mailing list
> > fipy@nist.gov
> >
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=oMf1-WHpUHeD7kN3S7612CDdF2TDPqDF9R-n-71Ks1Y&m=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yiO9QAq_3EUy9lVY5F76Tc9mnaB5Ubclil

Re: how to create large grid with holes and odd geometry

[James Pringle]

Well, I am motivated to give it a go, since I only have the summer to make
progress on project and it is blocking my research progress right now. Can
you give me a pointer to where the appropriate quantities are defined? I
can certainly write code to make the transformations, but it is a bit hard
without understanding precisely what is defined in the Mesh2D object. I
have made a simple Mesh2D object, but I am not sure which of the
attributes, etc, are absolutely required, or how they are precisely
defined.

Perhaps most useful for me would be

   1. a definition of which parts of the Mesh2D object must exist
   2. and the format of those parts, in particular the a face to vertex
   array and a
   cell to face array

Or you could point me to the appropriate part of the Gmsh code so I can go
from there. I presume poking around in fipy.Gmsh2D would be a good place to
start? Is there a better place to start?

I would love any documentation on the details of the Mesh2D object.

Thanks,
Jamie Pringle
University of New Hampshire

On Wed, Jun 15, 2016 at 12:21 PM, Daniel Wheeler 
wrote:

> Hi Jamie,
>
> There is no automated way to make a FiPy mesh from Scipy's Delaunay
> object. The problem is that FiPy needs a face to vertex array and a
> cell to face array while the Delaunay object just has the cell to
> vertex array. The functionality to extract the face to vertex array
> and cell to face array is in FiPy because it must be done for Gmsh,
> however, it is not abstracted in so that it can be reused.
>
> It is certainly possible to make the correct arrays with Numpy and
> pass them to the Mesh2D object, but it's a bit of work to write the
> code. If I find some time I might give it a go, but I don't know when
> I will get to it.
>
> Cheers,
>
> Daniel
>
> On Tue, Jun 14, 2016 at 4:15 PM, James Pringle  wrote:
> > Daniel et al. --
> >
> > As referenced earlier in this thread, I have a complex domain with
> holes
> > in it; I now have it broken up into triangles, based on the Delaunay
> package
> > in SciPy. I have
> >
> > Locations of all vertex coordinates,
> > list of vertices that make up faces
> > list of faces that make cells
> > list of faces that make up all the internal and external boundaries.
> >
> > Can you point me to any code or documentation that would help me
> understand
> > how to make this into a Mesh2D object? I am having a devil of a time
> > figuring out from the manual online. The best would be something that
> used
> > the output of either the triangles or scipy.spatial.Delaunay()
> packaged.  My
> > equation is of the form
> >
> > 0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi)
> >
> >
> > and I can get the coefficients A(x,y) and B(x,y) on either the faces or
> in
> > the cell centers are needed.
> >
> > Thank you.
> > Jamie Pringle
> > University of New Hampshire
>
> --
> Daniel Wheeler
> ___
> fipy mailing list
> fipy@nist.gov
>
> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=oMf1-WHpUHeD7kN3S7612CDdF2TDPqDF9R-n-71Ks1Y&m=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yiO9QAq_3EUy9lVY5F76Tc9mnaB5Ubclilpum_QCUOE&e=
>   [ NIST internal ONLY:
> https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=oMf1-WHpUHeD7kN3S7612CDdF2TDPqDF9R-n-71Ks1Y&m=Q6M6XDjzGDlxMDgrY-riZ7Q70Abl5Aq7SRrTwZNYsws&s=yIJJw-iQknimZfYCRrZmk7V2eHSJrCtcAehgVbGeiDk&e=
> ]
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: how to create large grid with holes and odd geometry

[James Pringle]

Daniel et al. --

As referenced earlier in this thread, I have a complex domain with
holes in it; I now have it broken up into triangles, based on the Delaunay
package in SciPy. I have

   1. Locations of all vertex coordinates,
   2. list of vertices that make up faces
   3. list of faces that make cells
   4. list of faces that make up all the internal and external boundaries.

Can you point me to any code or documentation that would help me understand
how to make this into a Mesh2D object? I am having a devil of a time
figuring out from the manual online. The best would be something that used
the output of either the triangles or scipy.spatial.Delaunay() packaged.
My equation is of the form

0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi)


and I can get the coefficients A(x,y) and B(x,y) on either the faces or in
the cell centers are needed.

Thank you.
Jamie Pringle
University of New Hampshire



On Wed, Jun 8, 2016 at 2:13 PM, Pringle, James 
wrote:

> Thank you. Because of the regular nature of the original data, it is easy
> to make it all triangles.
>
> Cheers,
> Jamie
>
>
> On Wed, Jun 8, 2016 at 2:10 PM, Guyer, Jonathan E. Dr. (Fed) <
> jonathan.gu...@nist.gov> wrote:
>
>> Meshes with holes are not a problem for FiPy. Daniel will be happy to
>> help you create a Mesh2D from the output of the triangle package.
>> Basically, you need a list of vertex coordinates, a list of vertex IDs that
>> make up faces, and a list of faces that make up cells. Having all triangles
>> should be pretty easy; it gets messy when you've got mixes of different
>> cell topologies.
>>
>> > On Jun 8, 2016, at 10:57 AM, James Pringle  wrote:
>> >
>> > Thank you; the axes were indeed in grid spacing, and your back of the
>> envelope calculation of sparsity was exactly correct.
>> >
>> > I am now playing around with the python triangle package to create a
>> triangular mesh with the appropriate holes and boundary defined. (e.g.
>> https://urldefense.proofpoint.com/v2/url?u=http-3A__dzhelil.info_triangle__data-2D4.png&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=oMf1-WHpUHeD7kN3S7612CDdF2TDPqDF9R-n-71Ks1Y&m=k71WczOQwQx8_RK9cU8wRK_X2spChhXjIrj5C0zWgw4&s=Qk6jAmccEUMVpP5TFa_gtI6T1b4bSBcRfmUOVIkzC-s&e=
>> ) .
>> >
>> > What is the best way to covert the triangle data (similar in form to
>> scipy.spatial.Delaunay output) into the form fipy likes -- or what is the
>> best documentation of the fipy mesh data structure?
>> >
>> > Does FiPy have the capability to deal with holes in a triangle mesh,
>> and have BC's on them? I did see the trick that Danial Wheeler mentioned.
>> Thanks!
>> >
>> > After I get the mesh, I was planing to reverse engineer the output of
>> fipy.meshes.gmshMesh.Gmsh2D to figure out how to get my mesh into fipy if
>> you don't mention something better.
>> >
>> > Thanks a bunch,
>> > Jamie
>> >
>> > On Wed, Jun 8, 2016 at 1:46 PM, Guyer, Jonathan E. Dr. (Fed) <
>> jonathan.gu...@nist.gov> wrote:
>> > If the domain were not so large and so sparse, I'd be inclined to
>> create a simple, rectilinear Grid2D of the full extent and then use known
>> coefficients to mask out (set B to zero?) the solution where you don't
>> know/care.
>> >
>> > Assuming the axes are labeled in grid spacings (?), then your mesh
>> would have around 5 million elements in it, with fewer than 1 million
>> actually being solved (although it looks like even less than 20% of the
>> domain is active). I don't think that would perform very well.
>> >
>> > I'm not thinking of anything clever with Gmsh off the top of my head.
>> >
>> > You could break the total domain into sub-grids, only instantiate the
>> corresponding Grid2D's if they're not empty, and then concatenate them
>> together. Sketching:
>> >
>> > A B C D
>> >  +-+-+-+-+
>> > 1|     | **  | | |
>> >  | |  ** | | |
>> >  +-+-+-+-+
>> > 2| |*|**   | |
>> >  | | |* *  | |
>> >  +-+-+-+-+
>> > 3| | | *** | |
>> >  | | |  ***| |
>> >  +-+-+-+-+
>> > 4| | |   **|*|
>> >  | | |*| *   |
>> >  +-+-+-+-+
>> >
>> >
>> > mesh = gridB1 + gridB2 + gridC2 + gridC3 + gridC4 + gridD4
>> >
>> >
>> >
>> >
>> > > On Jun 7, 2016, at 10:46

Re: how to create large grid with holes and odd geometry

[James Pringle]

Thank you. Because of the regular nature of the original data, it is easy
to make it all triangles.

Cheers,
Jamie

On Wed, Jun 8, 2016 at 2:10 PM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> Meshes with holes are not a problem for FiPy. Daniel will be happy to help
> you create a Mesh2D from the output of the triangle package. Basically, you
> need a list of vertex coordinates, a list of vertex IDs that make up faces,
> and a list of faces that make up cells. Having all triangles should be
> pretty easy; it gets messy when you've got mixes of different cell
> topologies.
>
> > On Jun 8, 2016, at 10:57 AM, James Pringle  wrote:
> >
> > Thank you; the axes were indeed in grid spacing, and your back of the
> envelope calculation of sparsity was exactly correct.
> >
> > I am now playing around with the python triangle package to create a
> triangular mesh with the appropriate holes and boundary defined. (e.g.
> https://urldefense.proofpoint.com/v2/url?u=http-3A__dzhelil.info_triangle__data-2D4.png&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=oMf1-WHpUHeD7kN3S7612CDdF2TDPqDF9R-n-71Ks1Y&m=k71WczOQwQx8_RK9cU8wRK_X2spChhXjIrj5C0zWgw4&s=Qk6jAmccEUMVpP5TFa_gtI6T1b4bSBcRfmUOVIkzC-s&e=
> ) .
> >
> > What is the best way to covert the triangle data (similar in form to
> scipy.spatial.Delaunay output) into the form fipy likes -- or what is the
> best documentation of the fipy mesh data structure?
> >
> > Does FiPy have the capability to deal with holes in a triangle mesh, and
> have BC's on them? I did see the trick that Danial Wheeler mentioned.
> Thanks!
> >
> > After I get the mesh, I was planing to reverse engineer the output of
> fipy.meshes.gmshMesh.Gmsh2D to figure out how to get my mesh into fipy if
> you don't mention something better.
> >
> > Thanks a bunch,
> > Jamie
> >
> > On Wed, Jun 8, 2016 at 1:46 PM, Guyer, Jonathan E. Dr. (Fed) <
> jonathan.gu...@nist.gov> wrote:
> > If the domain were not so large and so sparse, I'd be inclined to create
> a simple, rectilinear Grid2D of the full extent and then use known
> coefficients to mask out (set B to zero?) the solution where you don't
> know/care.
> >
> > Assuming the axes are labeled in grid spacings (?), then your mesh would
> have around 5 million elements in it, with fewer than 1 million actually
> being solved (although it looks like even less than 20% of the domain is
> active). I don't think that would perform very well.
> >
> > I'm not thinking of anything clever with Gmsh off the top of my head.
> >
> > You could break the total domain into sub-grids, only instantiate the
> corresponding Grid2D's if they're not empty, and then concatenate them
> together. Sketching:
> >
> > A B C D
> >  +-+-+-+-+
> > 1| | **  | | |
> >  | |  ** | | |
> >  +-+-+-+-+
> > 2| |*|**   | |
> >  | | |* *  | |
> >  +-+-+-+-+
> > 3| | | *** | |
> >  | | |  ***| |
> >  +-+-+-+-+
> > 4| | |   **|*|
> >  | | |*| *   |
> >  +-+-+-+-+
> >
> >
> > mesh = gridB1 + gridB2 + gridC2 + gridC3 + gridC4 + gridD4
> >
> >
> >
> >
> > > On Jun 7, 2016, at 10:46 AM, James Pringle  wrote:
> > >
> > > Dear mailing list & developers --
> > >
> > > I am looking for hints on the best way to proceed in creating a
> grid/mesh for a rather complex geometry. I am just looking for which method
> (Gmsh or something else?) to start with, so I can most efficiently start
> coding without exploring blind alleys.
> > >
> > > I am solving an elliptic/advective problem of the form
> > >
> > > 0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi)
> > >
> > > where Psi is the variable to solve for, and A(x,y) and B(x,y) are
> coefficients known on a set of discrete points shown as black in
> https://urldefense.proofpoint.com/v2/url?u=https-3A__dl.dropboxusercontent.com_u_382250_Grid01.png&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=oMf1-WHpUHeD7kN3S7612CDdF2TDPqDF9R-n-71Ks1Y&m=pa2VS0gonWbqlcwLDzR4QyP-iFHDpvnSHIqg-fk0QD4&s=ctz2HUgelmtveB6M5tJgv57BV8OqeCl3-MMJoIQOtmk&e=
> . The black appears solid because the grid is dense.
> > >
> > > The locations of the points where the coefficients are known define
> the grid. The number of points is large (911130 points) and they are evenly
> spaced where they exist. Note that there are holes 

Re: how to create large grid with holes and odd geometry

[James Pringle]

Thank you; the axes were indeed in grid spacing, and your back of the
envelope calculation of sparsity was exactly correct.

I am now playing around with the python triangle package to create a
triangular mesh with the appropriate holes and boundary defined. (e.g.
http://dzhelil.info/triangle//data-4.png ) .

What is the best way to covert the triangle data (similar in form to
scipy.spatial.Delaunay output) into the form fipy likes -- or what is the
best documentation of the fipy mesh data structure?

Does FiPy have the capability to deal with holes in a triangle mesh, and
have BC's on them? I did see the trick that Danial Wheeler mentioned.
Thanks!

After I get the mesh, I was planing to reverse engineer the output of
fipy.meshes.gmshMesh.Gmsh2D to figure out how to get my mesh into fipy if
you don't mention something better.

Thanks a bunch,
Jamie

On Wed, Jun 8, 2016 at 1:46 PM, Guyer, Jonathan E. Dr. (Fed) <
jonathan.gu...@nist.gov> wrote:

> If the domain were not so large and so sparse, I'd be inclined to create a
> simple, rectilinear Grid2D of the full extent and then use known
> coefficients to mask out (set B to zero?) the solution where you don't
> know/care.
>
> Assuming the axes are labeled in grid spacings (?), then your mesh would
> have around 5 million elements in it, with fewer than 1 million actually
> being solved (although it looks like even less than 20% of the domain is
> active). I don't think that would perform very well.
>
> I'm not thinking of anything clever with Gmsh off the top of my head.
>
> You could break the total domain into sub-grids, only instantiate the
> corresponding Grid2D's if they're not empty, and then concatenate them
> together. Sketching:
>
> A B C D
>  +-+-+-+-+
> 1| | **  | | |
>  | |  ** | | |
>  +-+-+-+-+
> 2| |*|**   | |
>  | | |* *  | |
>  +-+-+-+-+
> 3| | | *** | |
>  | | |  ***| |
>  +-+-+-+-+
> 4| | |   **|*|
>  | | |*| *   |
>  +-+-----+-+-+
>
>
> mesh = gridB1 + gridB2 + gridC2 + gridC3 + gridC4 + gridD4
>
>
>
>
> > On Jun 7, 2016, at 10:46 AM, James Pringle  wrote:
> >
> > Dear mailing list & developers --
> >
> > I am looking for hints on the best way to proceed in creating a
> grid/mesh for a rather complex geometry. I am just looking for which method
> (Gmsh or something else?) to start with, so I can most efficiently start
> coding without exploring blind alleys.
> >
> > I am solving an elliptic/advective problem of the form
> >
> > 0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi)
> >
> > where Psi is the variable to solve for, and A(x,y) and B(x,y) are
> coefficients known on a set of discrete points shown as black in
> https://urldefense.proofpoint.com/v2/url?u=https-3A__dl.dropboxusercontent.com_u_382250_Grid01.png&d=CwICAg&c=c6MrceVCY5m5A_KAUkrdoA&r=oMf1-WHpUHeD7kN3S7612CDdF2TDPqDF9R-n-71Ks1Y&m=pa2VS0gonWbqlcwLDzR4QyP-iFHDpvnSHIqg-fk0QD4&s=ctz2HUgelmtveB6M5tJgv57BV8OqeCl3-MMJoIQOtmk&e=
> . The black appears solid because the grid is dense.
> >
> > The locations of the points where the coefficients are known define the
> grid. The number of points is large (911130 points) and they are evenly
> spaced where they exist. Note that there are holes in the domain that
> represent actual islands in the ocean.
> >
> > I am happy to keep the resolution of the grid/mesh equal to the spacing
> of the points where the coefficients are known.
> >
> > What is the best way to approach creating a grid for this problem? I
> would love code, of course, but would be very happy with suggestions of the
> best way to start.
> >
> > Thanks
> > Jamie Pringle
> > University of New Hampshire
> > ___
> > fipy mailing list
>
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

how to create large grid with holes and odd geometry

[James Pringle]

Dear mailing list & developers --

I am looking for hints on the best way to proceed in creating a
grid/mesh for a rather complex geometry. I am just looking for which method
(Gmsh or something else?) to start with, so I can most efficiently start
coding without exploring blind alleys.

I am solving an elliptic/advective problem of the form

0=J(Psi,A(x,y)) + \Del(B(x,y)*\Del Psi)

where Psi is the variable to solve for, and A(x,y) and B(x,y) are
coefficients known on a set of discrete points shown as black in
https://dl.dropboxusercontent.com/u/382250/Grid01.png . The black appears
solid because the grid is dense.

The locations of the points where the coefficients are known define the
grid. The number of points is large (911130 points) and they are evenly
spaced where they exist. Note that there are holes in the domain that
represent actual islands in the ocean.

I am happy to keep the resolution of the grid/mesh equal to the spacing of
the points where the coefficients are known.

What is the best way to approach creating a grid for this problem? I would
love code, of course, but would be very happy with suggestions of the best
way to start.

Thanks
Jamie Pringle
University of New Hampshire
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: how to evaluate individual terms of equation

[James Pringle]

Thank you.
Jamie

On Thu, Feb 18, 2016 at 11:15 AM, Guyer, Jonathan E. Dr. <
jonathan.gu...@nist.gov> wrote:

> That's exactly what FiPy does if you pass in a CellVariable diffusion
> coefficient. There can be times you'll want to choose a different average,
> but if that's the case, you should have been doing it when you declared
> your DiffusionTerm.
>
> The ConvectionTerm is more complicated, as discussed in
> https://github.com/usnistgov/fipy/issues/461. The exact weighting between
> values at cell centers and values at faces depends on both the convection
> scheme you choose and on the Peclet number (which depends on the diffusion
> coefficient and how you interpolate that). What you've written is about as
> good as you're going to get without doing a *lot* more work.
>
> You should not try to cast convCoeff or Psi.faceGrad to a cell-centered
> value, as we provide no divergence operator on cell vectors.
>
> On Feb 18, 2016, at 9:35 AM, James Pringle  wrote:
>
> > p.s. the following code works fairly well --
> >
> > crossIsoAdvec_modGrid=(Psi.arithmeticFaceValue *
> convCoeff).divergence
> > botFricCurl_modGrid=(DiffCoeff.arithmeticFaceValue *
> Psi.faceGrad).divergence
> >
> > I just want to make sure it is somewhat sensible...
> >
> > On Thu, Feb 18, 2016 at 9:28 AM, James Pringle  wrote:
> > Dear Jon --
> >
> > I am having a bit of a problem with your technique. All of my
> coefficients vary in space, e.g. they are define as follows:
> >
> > mesh = Grid2D(Lx=Lx,Ly=Ly, nx=nx, ny=ny)
> >
> > convCoeff=FaceVariable(mesh=mesh,rank=1)
> > convCoeff.value[0,:]=
> > convCoeff.value[1,:]= 
> >
> > DiffCoeff=CellVariable(mesh=mesh, rank=0)
> > DiffCoeff.value = ...
> >
> > eq =
> (DiffusionTerm(var=Psi,coeff=DiffCoeff)+ExponentialConvectionTerm(var=Psi,coeff=convCoeff))
> > eq.solve(var=Psi)
> >
> > So when I try to implement your suggestions of (DiffCoeff *
> Psi.faceGrad).divergence or (Psi * convCoeff).divergence I get errors
> because DiffCoeff is a Cell variable, and Psi.faceGrad is a face variable.
> Likewise for Psi*convCoeff
> >
> > So what is the most numerically constant way to solve this issue -- map
> DiffCoeff to face, or Psi.faceGrad to cell? Likewise should I map convCoeff
> to cell or Psi to face? I want to do it in the way that is most consistent
> with the underlying numerics.
> >
> > Thanks
> > Jamie
> >
> > On Wed, Feb 17, 2016 at 9:43 AM, James Pringle  wrote:
> > Thank you, Jon, that is helpful. I would love to have the ability to
> find the individual terms of the equation as calculated in the matrix that
> is solved by fiPy. This would make closing budgets, etc, easier. But your
> solution is a good start.
> >
> > Thanks again
> > Jamie Pringle
> > University of New Hampshire
> >
> > On Tue, Feb 16, 2016 at 4:53 PM, Guyer, Jonathan E. Dr. <
> jonathan.gu...@nist.gov> wrote:
> > James -
> >
> > I don't think there's a straightforward way to get at this.
> >
> > You can certainly write the explicit forms, e.g.,
> >
> >   (DiffCoeff * Psi.faceGrad).divergence
> >
> > or
> >
> >   (Psi * convCoeff).divergence
> >
> > but this isn't exactly the same as the matrix FiPy builds, as discussed
> at https://github.com/usnistgov/fipy/issues/461.
> >
> > I can see the value for diagnosing and simply understanding mechanisms,
> so I'll think about ways we might provide access to this.
> >
> > - Jon
> >
> > On Feb 12, 2016, at 10:18 AM, James Pringle  wrote:
> >
> > > I feel this should be an elementary question, but I can't seem to
> figure out how to answer it. I am solving a simple linear elleptic-ish
> equation with
> > >
> > > eq =
> (DiffusionTerm(var=Psi,coeff=DiffCoeff)+ExponentialConvectionTerm(var=Psi,coeff=convCoeff))
> > > eq.solve(var=Psi)
> > >
> > > This works fine; the solution matches what I would expect.  I have two
> questions:
> > >
> > > First, how can I obtain the value of the individual terms of the
> equation, evaluated with the solution in Psi?
> > >
> > > Second, is there any way to define my own new Psi (which is not an
> exact solution to the equation), and easily evaluate the DiffusionTerm and
> ExponentialConvectionTerm for that Psi?
> > >
> > > I am trying to illustrate how various parts of the solution evolve
> over the solution space, and how various approximations to the full
> solution differ from the ful

Re: how to evaluate individual terms of equation

[James Pringle]

p.s. the following code works fairly well --

crossIsoAdvec_modGrid=(Psi.arithmeticFaceValue * convCoeff).divergence
botFricCurl_modGrid=(DiffCoeff.arithmeticFaceValue *
Psi.faceGrad).divergence

I just want to make sure it is somewhat sensible...

On Thu, Feb 18, 2016 at 9:28 AM, James Pringle  wrote:

> Dear Jon --
>
> I am having a bit of a problem with your technique. All of my
> coefficients vary in space, e.g. they are define as follows:
>
> mesh = Grid2D(Lx=Lx,Ly=Ly, nx=nx, ny=ny)
>
> convCoeff=FaceVariable(mesh=mesh,rank=1)
> convCoeff.value[0,:]=
> convCoeff.value[1,:]= 
>
> DiffCoeff=CellVariable(mesh=mesh, rank=0)
> DiffCoeff.value = ...
>
> eq =
> (DiffusionTerm(var=Psi,coeff=DiffCoeff)+ExponentialConvectionTerm(var=Psi,coeff=convCoeff))
> eq.solve(var=Psi)
>
> So when I try to implement your suggestions of (DiffCoeff *
> Psi.faceGrad).divergence or (Psi * convCoeff).divergence I get errors
> because DiffCoeff is a Cell variable, and Psi.faceGrad is a face variable.
> Likewise for Psi*convCoeff
>
> So what is the most numerically constant way to solve this issue -- map
> DiffCoeff to face, or Psi.faceGrad to cell? Likewise should I map convCoeff
> to cell or Psi to face? I want to do it in the way that is most consistent
> with the underlying numerics.
>
> Thanks
> Jamie
>
> On Wed, Feb 17, 2016 at 9:43 AM, James Pringle  wrote:
>
>> Thank you, Jon, that is helpful. I would love to have the ability to find
>> the individual terms of the equation as calculated in the matrix that is
>> solved by fiPy. This would make closing budgets, etc, easier. But your
>> solution is a good start.
>>
>> Thanks again
>> Jamie Pringle
>> University of New Hampshire
>>
>> On Tue, Feb 16, 2016 at 4:53 PM, Guyer, Jonathan E. Dr. <
>> jonathan.gu...@nist.gov> wrote:
>>
>>> James -
>>>
>>> I don't think there's a straightforward way to get at this.
>>>
>>> You can certainly write the explicit forms, e.g.,
>>>
>>>   (DiffCoeff * Psi.faceGrad).divergence
>>>
>>> or
>>>
>>>   (Psi * convCoeff).divergence
>>>
>>> but this isn't exactly the same as the matrix FiPy builds, as discussed
>>> at https://github.com/usnistgov/fipy/issues/461.
>>>
>>> I can see the value for diagnosing and simply understanding mechanisms,
>>> so I'll think about ways we might provide access to this.
>>>
>>> - Jon
>>>
>>> On Feb 12, 2016, at 10:18 AM, James Pringle  wrote:
>>>
>>> > I feel this should be an elementary question, but I can't seem to
>>> figure out how to answer it. I am solving a simple linear elleptic-ish
>>> equation with
>>> >
>>> > eq =
>>> (DiffusionTerm(var=Psi,coeff=DiffCoeff)+ExponentialConvectionTerm(var=Psi,coeff=convCoeff))
>>> > eq.solve(var=Psi)
>>> >
>>> > This works fine; the solution matches what I would expect.  I have two
>>> questions:
>>> >
>>> > First, how can I obtain the value of the individual terms of the
>>> equation, evaluated with the solution in Psi?
>>> >
>>> > Second, is there any way to define my own new Psi (which is not an
>>> exact solution to the equation), and easily evaluate the DiffusionTerm and
>>> ExponentialConvectionTerm for that Psi?
>>> >
>>> > I am trying to illustrate how various parts of the solution evolve
>>> over the solution space, and how various approximations to the full
>>> solution differ from the full solution.
>>> >
>>> > Thanks to all who glance at this
>>> > James Pringle
>>> > University of New Hampshire
>>> >
>>> > ___
>>> > 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: how to evaluate individual terms of equation

[James Pringle]

Dear Jon --

I am having a bit of a problem with your technique. All of my
coefficients vary in space, e.g. they are define as follows:

mesh = Grid2D(Lx=Lx,Ly=Ly, nx=nx, ny=ny)

convCoeff=FaceVariable(mesh=mesh,rank=1)
convCoeff.value[0,:]=
convCoeff.value[1,:]= 

DiffCoeff=CellVariable(mesh=mesh, rank=0)
DiffCoeff.value = ...

eq =
(DiffusionTerm(var=Psi,coeff=DiffCoeff)+ExponentialConvectionTerm(var=Psi,coeff=convCoeff))
eq.solve(var=Psi)

So when I try to implement your suggestions of (DiffCoeff *
Psi.faceGrad).divergence or (Psi * convCoeff).divergence I get errors
because DiffCoeff is a Cell variable, and Psi.faceGrad is a face variable.
Likewise for Psi*convCoeff

So what is the most numerically constant way to solve this issue -- map
DiffCoeff to face, or Psi.faceGrad to cell? Likewise should I map convCoeff
to cell or Psi to face? I want to do it in the way that is most consistent
with the underlying numerics.

Thanks
Jamie

On Wed, Feb 17, 2016 at 9:43 AM, James Pringle  wrote:

> Thank you, Jon, that is helpful. I would love to have the ability to find
> the individual terms of the equation as calculated in the matrix that is
> solved by fiPy. This would make closing budgets, etc, easier. But your
> solution is a good start.
>
> Thanks again
> Jamie Pringle
> University of New Hampshire
>
> On Tue, Feb 16, 2016 at 4:53 PM, Guyer, Jonathan E. Dr. <
> jonathan.gu...@nist.gov> wrote:
>
>> James -
>>
>> I don't think there's a straightforward way to get at this.
>>
>> You can certainly write the explicit forms, e.g.,
>>
>>   (DiffCoeff * Psi.faceGrad).divergence
>>
>> or
>>
>>   (Psi * convCoeff).divergence
>>
>> but this isn't exactly the same as the matrix FiPy builds, as discussed
>> at https://github.com/usnistgov/fipy/issues/461.
>>
>> I can see the value for diagnosing and simply understanding mechanisms,
>> so I'll think about ways we might provide access to this.
>>
>> - Jon
>>
>> On Feb 12, 2016, at 10:18 AM, James Pringle  wrote:
>>
>> > I feel this should be an elementary question, but I can't seem to
>> figure out how to answer it. I am solving a simple linear elleptic-ish
>> equation with
>> >
>> > eq =
>> (DiffusionTerm(var=Psi,coeff=DiffCoeff)+ExponentialConvectionTerm(var=Psi,coeff=convCoeff))
>> > eq.solve(var=Psi)
>> >
>> > This works fine; the solution matches what I would expect.  I have two
>> questions:
>> >
>> > First, how can I obtain the value of the individual terms of the
>> equation, evaluated with the solution in Psi?
>> >
>> > Second, is there any way to define my own new Psi (which is not an
>> exact solution to the equation), and easily evaluate the DiffusionTerm and
>> ExponentialConvectionTerm for that Psi?
>> >
>> > I am trying to illustrate how various parts of the solution evolve over
>> the solution space, and how various approximations to the full solution
>> differ from the full solution.
>> >
>> > Thanks to all who glance at this
>> > James Pringle
>> > University of New Hampshire
>> >
>> > ___
>> > 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: how to evaluate individual terms of equation

[James Pringle]

Thank you, Jon, that is helpful. I would love to have the ability to find
the individual terms of the equation as calculated in the matrix that is
solved by fiPy. This would make closing budgets, etc, easier. But your
solution is a good start.

Thanks again
Jamie Pringle
University of New Hampshire

On Tue, Feb 16, 2016 at 4:53 PM, Guyer, Jonathan E. Dr. <
jonathan.gu...@nist.gov> wrote:

> James -
>
> I don't think there's a straightforward way to get at this.
>
> You can certainly write the explicit forms, e.g.,
>
>   (DiffCoeff * Psi.faceGrad).divergence
>
> or
>
>   (Psi * convCoeff).divergence
>
> but this isn't exactly the same as the matrix FiPy builds, as discussed at
> https://github.com/usnistgov/fipy/issues/461.
>
> I can see the value for diagnosing and simply understanding mechanisms, so
> I'll think about ways we might provide access to this.
>
> - Jon
>
> On Feb 12, 2016, at 10:18 AM, James Pringle  wrote:
>
> > I feel this should be an elementary question, but I can't seem to figure
> out how to answer it. I am solving a simple linear elleptic-ish equation
> with
> >
> > eq =
> (DiffusionTerm(var=Psi,coeff=DiffCoeff)+ExponentialConvectionTerm(var=Psi,coeff=convCoeff))
> > eq.solve(var=Psi)
> >
> > This works fine; the solution matches what I would expect.  I have two
> questions:
> >
> > First, how can I obtain the value of the individual terms of the
> equation, evaluated with the solution in Psi?
> >
> > Second, is there any way to define my own new Psi (which is not an exact
> solution to the equation), and easily evaluate the DiffusionTerm and
> ExponentialConvectionTerm for that Psi?
> >
> > I am trying to illustrate how various parts of the solution evolve over
> the solution space, and how various approximations to the full solution
> differ from the full solution.
> >
> > Thanks to all who glance at this
> > James Pringle
> > University of New Hampshire
> >
> > ___
> > 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 ]

how to evaluate individual terms of equation

[James Pringle]

I feel this should be an elementary question, but I can't seem to figure
out how to answer it. I am solving a simple linear elleptic-ish equation
with

eq =
(DiffusionTerm(var=Psi,coeff=DiffCoeff)+ExponentialConvectionTerm(var=Psi,coeff=convCoeff))
eq.solve(var=Psi)


This works fine; the solution matches what I would expect.  I have two
questions:

First, how can I obtain the value of the individual terms of the equation,
evaluated with the solution in Psi?

Second, is there any way to define my own new Psi (which is not an exact
solution to the equation), and easily evaluate the DiffusionTerm and
ExponentialConvectionTerm for that Psi?

I am trying to illustrate how various parts of the solution evolve over the
solution space, and how various approximations to the full solution differ
from the full solution.

Thanks to all who glance at this
James Pringle
University of New Hampshire
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]