gmane.comp.python.fipy

Coupled PDE and Non-Liner Algebraic Equations

2012-05-22T08:46:11

Dear FiPy developpers and users,

I'm trying to solve the following system of PDE and Non-Liner Algebraic 
Equations (NLAE):

the PDEs read:
\frac{\partial{\phi}}{\partial t} = \nabla ( D(y) \nabla y)
where both phi and y are arrays of cellVariable and the number of PDE 
equals the length of the \phi array

and the set of NLAE

F(\phi, y)= 0 at every node in the domain
where y is the local root of F, i.e. \frac{\partial F}{\partial y} is a 
square full matrix.
This local problem is non-linear and I use SciPy.optimize.fsolve as 
solver.


My problem:

I don't know an explict way of writing y = G{\phi} (and its derivatives) 
as suggested for example here:
http://www.ctcms.nist.gov/fipy/examples/phase/generated/examples.phase.simple.html?highlight=expansion%20source

Therefore I wrote a local solver for F(\phi, y) (+ numerical derivatives) 
and then use a series of setValue() getValue() to interact with FiPy 
variables and solve the PDE sequentially with sweeps....and the whole 
process becomes very slow as you

please unsubscribe me from the mailing list

Yeung Au — 2012-05-22T01:15:49

Thanks.
_______________________________________________
fipy mailing list
fipy-R3+/ord2DXQ< at >public.gmane.org
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: RE-visiting advection-diffusion equations in 2D

Yun Tao — 2012-05-22T01:11:06

Very interesting stuff! For others' reference, in fipy 2.2-dev5124, there
are the following syntax changes:

mesh._getOrientedAreaProjections()    ->    mesh.orientedAreaProjections
mesh._getOrientedFaceNormals()    ->    mesh.orientedFaceNormals
conv._getGeomCoeff(mesh)    ->    conv._calcGeomCoeff(var)
eq._getDiffusiveGeomCoeff(mesh)    ->    diff._calcGeomCoeff(var)

Two follow-up questions:

1. When I called mesh._getCellDistances() of a 2D mesh (nx=ny=3, dx=dy=1),
the values in the array are non-uniform:




This seems odd given how I imaged the grid to be, as a square divided
equally into 9 parts since if the latter's true, wouldn't the distances all
be equal to 1?

2. My version of the code is attached here. Though modifying either the
advection or the diffusion coefficients doesn't seem to affect the plot. Is
that because I'm not using a specific vector viewer?

Thanks

1. syntax change under fipy 2.2dev
2. mesh._getCellDistances() not uniform
3. changing advection or diffusion coefficients don't alt

Re: PCG Solver

Daniel Wheeler — 2012-05-21T21:26:30

It's just a deprecation warning from pysparse. Nothing to worry about
though. I think this is because we haven't updated the pysparse calls
so that they use the new pysparse interface.

Cheers

On Mon, May 21, 2012 at 4:36 PM, wang yunbo <wang1320-olO2ZdjDehc3uPMLIKxrzw< at >public.gmane.org> wrote:

PCG Solver

wang yunbo — 2012-05-21T20:36:39

Hi there,

I need to use PCG solver in my code for fipy trunk. 

I built the solver like this:
solver = LinearPCGSolver(tolerance=1.0e-15, iterations=1000)

and for each step:
eqch.solve(concentration, boundaryConditions=BCs, solver=solver, dt=dt)

Here I got a deprecation warning:

/home/wang1320/build/lib/python2.6/site-packages/FiPy-3.0_dev5154-py2.6.egg/fipy/solvers/pysparse/pysparseSolver.py:76: DeprecationWarning: PyArray_FromDimsAndDataAndDescr: use PyArray_NewFromDescr.
  self.iterations, P)

What should I change to make the PCG solver work in this case?  Thank you.


Best,

Yunbo

 
_______________________________________________
fipy mailing list
fipy-R3+/ord2DXQ< at >public.gmane.org
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: Boundary condition in fipy trunk

Jonathan Guyer — 2012-05-21T19:13:25

On May 21, 2012, at 2:17 PM, wang yunbo wrote:



m.exteriorFaces or m.getExteriorFaces(), but not m.exteriorFaces()


_______________________________________________
fipy mailing list
fipy-R3+/ord2DXQ< at >public.gmane.org
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Boundary condition in fipy trunk

wang yunbo — 2012-05-21T18:17:52

Hi there,

I want to add a fixed-flux boundary condition like this in fipy trunk, 

BCs = (FixedFlux(faces=m.exteriorFaces(),value=0),) 

An error is encountered. 

Traceback (most recent call last):
  File "phasefieldtest.py", line 67, in <module>
    BCs = (FixedFlux(faces=m.exteriorFaces(),value=0),)
  File "/home/wang1320/build/lib/python2.6/site-packages/FiPy-3.0_dev5154-py2.6.egg/fipy/boundaryConditions/fixedFlux.py", line 65, in __init__
    BoundaryCondition.__init__(self,faces,value)
  File "/home/wang1320/build/lib/python2.6/site-packages/FiPy-3.0_dev5154-py2.6.egg/fipy/boundaryConditions/boundaryCondition.py", line 64, in __init__
    if not (self.faces | self.faces.mesh.exteriorFaces 
AttributeError: 'numpy.ndarray' object has no attribute 'mesh'

What do I need to change to make the boundary conditions work?


Thanks.

Best,

Yunbo
_______________________________________________
fipy mailing list
fipy-R3+/ord2DXQ< at >public.gmane.org
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://ema

Re: Solving a convection diffusion reaction equation

Jonathan Guyer — 2012-05-21T15:08:12

On May 19, 2012, at 8:09 AM, Matej Svejda wrote:


Definitely slower, although for 1D you may not notice it much.

Basically, you should never (and should never have to) iterate over cells directly. You want to write vectorized expressions and let NumPy do the iterations for you. 

At a minimum, your H update inside your step loop should look like

  alpha = mesh.getCellCenters()[0]
  H.setValue(h(alpha, t))

but even better is to make the following definition *before* your step loop and then let FiPy take care of updating H as t changes.

  alpha = CellVariable(mesh=mesh, value=mesh.getCellCenters()[0])
  t = Variable()
  H = h(alpha, t)

For your renormalization, it's much more efficient to write

  phi.setValue(value=phi.getValue() / phi.getCellVolumeAverage())



You need to be sure that phi is defined with hasOld=True and you need to perform

  phi.updateOld() 

at the beginning of your step loop but outside your sweep loop. 

See http://www.ctcms.nist.gov/fipy/documentation/FAQ.html#iterations-timeste

Re: Solving a convection diffusion reaction equation

Matej Svejda — 2012-05-19T12:09:00

First of all: Thank you so much for your help!

Thanks for the tip, but it didn't seem to change much...


I am updating the function H in the following way:

newHValues = [0.0] * len(centers)
for i in range(len(centers)):
alpha = centers[i]
newHValues[i] = h(alpha, t)
H.setValue(newHValues)

The normalization is performed the same way:

for value in values:
normalization += value * dAlpha
for i in range(length):
newValues[i] = values[i] / normalization
phi.setValue(value=newValues)

Is this slover than using variables?


How would I combine the renormalization and the sweeping? I tried
something along the lines of:

for step in range(timeSteps):
#update the H function variable
for i in range(sweepSteps):
eq.sweep(var=phi, boundaryConditions=BCs, dt=timeStepDuration)
normalize(phi, dAlpha)
phi.updateOld()

Where normalize directly sets the value of phi (as can be seen in the
code I posted earlier). Is this the right way? It doesn't seem to make
any difference...

When compared to a solution that I

Re: Solving a convection diffusion reaction equation

Jonathan Guyer — 2012-05-17T12:38:44

On May 16, 2012, at 5:05 PM, Matej Svejda wrote:


I would use .getFaceGrad() here, not .getFaceGradAverage(). I didn't realize we even had such a function (for five years!), and contrary to its documentation, I think it has to be 1st order accurate, not 2nd. For what it's worth, it doesn't appear that we use it anywhere.



Manually updating phi as a result of your renormalization is probably necessary, although it's possible that you're not doing what you intend. If you show me the code for your update, I can check that it has the effect you want.

I don't think it should be necessary to manually update H if you define it as a Variable expression and if elapsed time is defined as a Variable. There are a number of ways that you might not be getting that, though, so it's again probably best to just show us the code you use to define and to update H.

That said, it's not the worst thing in the world to manually update H. It just clutters up your iteration loop.



Switching from .getFaceGradAverage() to .get

Solving a convection diffusion reaction equation

Matej Svejda — 2012-05-16T21:05:09

Hi,

I'm having some problems with solving a PDE using FiPy. The equation
is of the form:

C * \frac{\partial \varphi(\alpha, t)}{\partial t} = \frac{\partial^2
\varphi(\alpha, t)}{\partial \alpha^2} + \frac{\partial}{\partial
\alpha}  (\frac{\partial H(\alpha, t)}{\partial \alpha} *
\varphi(\alpha, t))

Where C is a constant and H(\alpha, t) is a (rather complicated)
function, defined piecewise on each interval of \alpha. Right know I'm
solving it with FiPy in the following way:

eq = TransientTerm(coeff = C) == ImplicitDiffusionTerm() +
ExponentialConvectionTerm(coeff=H.getFaceGradAverage())

And then I simply iterate over the time-interval using
eq.solve(var=phi, boundaryConditions=bcs, dt=timeStepDuration). For
each time-step I update the CellVariable H manually (since I couldn't
manage for it to be updated automatically by updating the
time-variable). I also perform a "normalization" after each time step
of the solution phi, since it represents a probability distribution
and therefore the integral over

Re: RE-visiting advection-diffusion equations in 2D

Jonathan Guyer — 2012-05-16T13:39:11

On May 15, 2012, at 5:15 PM, Yun Tao wrote:


The actual calculation of the Peclet number is performed at

http://matforge.org/fipy/browser/branches/version-2_1/fipy/terms/convectionTerm.py#L143

The calculation of the advection strength is at 

http://matforge.org/fipy/browser/branches/version-2_1/fipy/terms/convectionTerm.py#L123

That code is really old and so now I'd write

  numerix.dot(self.coeff, mesh._getOrientedAreaProjections())



As to visualization, that's going to be tricky. The Peclet number you calculate this way is a scalar FaceVariable and we don't have any tools to visualize such a thing. I'd be inclined to just visualize the vector field of u * dAP /D (you can get dAP with mesh._getCellDistances()). Hmmm... that will tend to over-represent advection vectors that lie parallel to a face. Maybe better to calculate the scalar Peclet number and then multiply it by the face normals to get a vector field that FiPy can plot, e.g.,


_______________________________________________
fipy mailing li

Re: Second derivative of function in Source

Matej Svejda — 2012-05-16T08:34:57

Thank you very much for your help :)

On Mon, May 14, 2012 at 3:47 PM, Jonathan Guyer <guyer-R3+/ord2DXQ< at >public.gmane.org> wrote:

_______________________________________________
fipy mailing list
fipy-R3+/ord2DXQ< at >public.gmane.org
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: RE-visiting advection-diffusion equations in 2D

Yun Tao — 2012-05-15T21:15:31

Thank you both, Jonathan and Daniel, for the most remarkably detailed
assistance and insights!! Your replies have clarified numerous confusions
someone who's a novice in numerical analysis would probably have. I ran
through the routine again with the fixes made and the video is attached. As
expected, only the CentralDifference method fails under a high Peclet
value. By the way, is there a way for FiPy to calculate and display that?
According to the FiPy doc
http://www.ctcms.nist.gov/fipy/documentation/numerical/discret.html#linear-equations,
the advection strength requires the dot product of the advection component
and the normal vector at each face, yet the latter value is not obvious to
me as an end user not knowing the topology of the solution domain.

Cheers,

Yun

On Wed, May 2, 2012 at 8:57 AM, Jonathan Guyer <guyer-R3+/ord2DXQ< at >public.gmane.org> wrote:

Re: bridging 2D mesh with bivariate_normal from matplotlib

Yun Tao — 2012-05-14T19:11:37

Oops, don't know how I missed that. Guess bivariate_normal()'s more
flexible than I thought. Thanks!

On Mon, May 14, 2012 at 6:37 AM, Jonathan Guyer <guyer-R3+/ord2DXQ< at >public.gmane.org> wrote:

Re: Second derivative of function in Source

Jonathan Guyer — 2012-05-14T13:47:34

On May 12, 2012, at 9:44 AM, Matej Svejda wrote:


I think this is covered in http://www.ctcms.nist.gov/fipy/documentation/FAQ.html#how-do-i-represent-a-term-that-doesn-t-involve-the-dependent-variable



H.getGrad().getGrad() should be rank-2 (and that's what I get with both FiPy 2.1 and trunk/)


The gradient of a scalar is a vector. The gradient of a vector is a rank-2 tensor. The dimensionality doesn't matter. Working in 1D can often mislead like this.

\nabla^2 H is not equivalent to \nabla(\nabla H). It is equivalent to \nabla\cdot(\nabla H), and in FiPy you do that with H.getFaceGrad().getDivergence()


_______________________________________________
fipy mailing list
fipy-R3+/ord2DXQ< at >public.gmane.org
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: bridging 2D mesh with bivariate_normal from matplotlib

Jonathan Guyer — 2012-05-14T13:37:05



I don't believe you need to reshape and flatten. bivariate_normal() returns the same shape that you pass in, so all you need to do is

x, y = mesh.getCellCenters()
 
z = bivariate_normal(x, y, 1., 1., 0., 0.)
 
phi = CellVariable(mesh=mesh, value=z)


_______________________________________________
fipy mailing list
fipy-R3+/ord2DXQ< at >public.gmane.org
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Second derivative of function in Source

Matej Svejda — 2012-05-12T13:44:41

Hello,

I'm using FiPy to solve a Convection Diffusion Reaction Equation (one
spacial dimension, plus time). I want my source to:

da^2(H(a, t)) * P(a, t)

Where da^2(H(a,t)) is the second special derivative of a function H (a
known function and NOT  the function I'm solving for) and P the
function I'm solving for. Both are represented as CellVariable.

I have tried to write: H.getGrad().getGrad() * P but the ranks don't
match (H has rank-0 and H.getGrad().getGrad() has rank-3). I don't get
why the rank is increased when I take the gradient of a function that
is defined on a 1-d mesh. Any ideas?

Thanks in advance,
Matej
_______________________________________________
fipy mailing list
fipy-R3+/ord2DXQ< at >public.gmane.org
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: bridging 2D mesh with bivariate_normal from matplotlib

Yun Tao — 2012-05-11T22:01:19

attachment), thanks to your help. For those who're curious, the crucial
thing to do is to first reshape the axial variables, feed it through
bivariate normal, and finally flatten the result for CellVariable. To
illustrate:

x, y = mesh.getCellCenters()








Cheers,

Yun




_______________________________________________

Re: Request: Scharfetter-Gummel discretization forsolvingdrift-diffusion equations of pn diode

Jonathan Guyer — 2012-05-10T18:01:25

On May 10, 2012, at 6:36 AM, jimmy.123 wrote:


I've posted a script that I wrote quite awhile ago to http://matforge.org/fipy/wiki/ScharfetterGummel. It's not the way my current PV codes work, but they require a much more complicated infrastructure that I'm not quite ready to release yet.

Although that script still runs with FiPy 2.1, it diverges with the PySparse solvers and the current is not stable when using the Trilinos solvers. I don't have time to diagnose why, but the code should at least give you an idea of how to go about it.

In this case, I implemented Scharfetter-Gummel as an interpolation scheme. It calculates the concentration at a face based on the Fermi levels in the neighboring cells. This is a little bit different from the more typical presentation of Scharfetter-Gummel as a specialized discretization of the fluxes, but the end effect is essentially the same.
_______________________________________________
fipy mailing list
fipy-R3+/ord2DXQ< at >public.gmane.org
http://www.ctcms.nist.gov/fip

Request: Scharfetter-Gummel discretization for solvingdrift-diffusion equations of pn diode

jimmy.123 — 2012-05-10T10:36:05

 Dear guys,

I am trying to solve the drift-diffusion equations of electrons and holes in a semiconductor, in which Scharfetter-Gummel discretization is needed. Hence, I would like to kindly ask whether it is possible for you (or one related) to provide me one script to solve the drift-diffusion equations or a script to calculate the current density from carrier density and electric field by using Scharfetter-Gummel discretization? (I supposed it is possible because I saw in one previous reply it was mentioned request of the script was possible).

If there is no available script at present, could you please provide an example on the scharfetterGummelFaceVariable module so that I can understand it and make further changes based on it?

Thank you for the help.

Best regards,
Jiaguo Zhang
_______________________________________________
fipy mailing list
fipy-R3+/ord2DXQ< at >public.gmane.org
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]