Re: fipy diffusion setup
On Tue, Apr 19, 2016 at 4:10 PM, Phil Battle wrote: > Hi Daniel > Thank you for following up- > I did try your suggestion, ( C is my concentration dependent Cell variable) > > Dx=Dpe_x*(a_x+((1.-a_x)/(b_x*C+g_x))) > Dz=Dpe_z*(a_z+((1.-a_z)/(b_z*C+g_z))) > diff = Dx.dot([[1,0],[0,0]])+Dz.dot([[0,0],[0,1]]) > > and program runs and "looks correct" so good to go!! Great news. > But before your suggestion I was trying to create the suggested rank 2 > diffusion coefficient a different way and I did not get far before hitting a > road block- If you have time to answer this specific question it may help > with understanding the structure of variables in this program > > I created a 2x2 Cell Variable on a 2x2 2Dim mesh. I can show that code if > necessary, but the result for the initialized CellVariable is shown below > > print(D.value) > > [[[ 0. 2. 4. 6.] > [ 1. 1. 1. 1.]] > > [[ 1. 1. 1. 1.] > [ 1. 1. 1. 1.]]] > > > Now the part that I am having trouble with. How can I change the 6 to a 0.? I think you just need to do D[0, 0, 3] = 6. However, it isn't always a good practice to access the spatial index since if the grid changes then the index has a different spatial location > The closet I can come is changing all the values at that mesh position to 0. > that is using > > D.setValue(0,where=[[0,0,0,1]]) I'm not sure why the [[0, 0, 0, 1]] does what it does. Numpy does the following: ~~~ In [9]: a = np.arange(16).reshape([2, 2, 4]) In [10]: a Out[10]: array([[[ 0, 1, 2, 3], [ 4, 5, 6, 7]], [[ 8, 9, 10, 11], [12, 13, 14, 15]]]) In [12]: a[[[0, 0, 0, 1]]] = -1 In [13]: print a [[[-1 -1 -1 -1] [-1 -1 -1 -1]] [[-1 -1 -1 -1] [-1 -1 -1 -1]]] ~~~ If you're trying to do fancy indexing with FiPy, you might want check with Numpy first so that you can trust the FiPy result. -- 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: fipy diffusion setup
Hi Daniel Thank you for following up- I did try your suggestion, ( C is my concentration dependent Cell variable) Dx=Dpe_x*(a_x+((1.-a_x)/(b_x*C+g_x))) Dz=Dpe_z*(a_z+((1.-a_z)/(b_z*C+g_z))) diff = Dx.dot([[1,0],[0,0]])+Dz.dot([[0,0],[0,1]]) and program runs and "looks correct" so good to go!! But before your suggestion I was trying to create the suggested rank 2 diffusion coefficient a different way and I did not get far before hitting a road block- If you have time to answer this specific question it may help with understanding the structure of variables in this program I created a 2x2 Cell Variable on a 2x2 2Dim mesh. I can show that code if necessary, but the result for the initialized CellVariable is shown below print(D.value) [[[ 0. 2. 4. 6.] [ 1. 1. 1. 1.]] [[ 1. 1. 1. 1.] [ 1. 1. 1. 1.]]] Now the part that I am having trouble with. How can I change the 6 to a 0.? The closet I can come is changing all the values at that mesh position to 0. that is using D.setValue(0,where=[[0,0,0,1]]) print(D.value) [[[ 0. 2. 4. 0.] [ 1. 1. 1. 0.]] [[ 1. 1. 1. 0.] [ 1. 1. 1. 0.]]] again what I was angling for was an output that looked like [[[ 0. 2. 4. 0.] [ 1. 1. 1. 1.]] [[ 1. 1. 1. 1.] [ 1. 1. 1. 1.]]] thanks for your help Phil On 4/19/2016 6:51 AM, Daniel Wheeler wrote: Phil, just to follow up and complete the answer. I neglected to explain how to create a rank 2 CellVariable from two rank 0 CellVariables to preserve the dependency. The dot operator helps with that. Given two rank 0 CellVariables, v00 and v11, then to get the rank 2 CellVariable use, vv = v00.dot([[1, 0], [0, 0]]) + v11.dot([[0, 0], [0, 1]]) On Mon, Apr 18, 2016 at 3:33 PM, Daniel Wheeler wrote: On Mon, Apr 18, 2016 at 12:33 PM, Phil Battle wrote: Ok, below works as I would expect ( note Dx and Dz are constants) Thanks for that. I can confirm that it's broken in Python 2.7 as well when the diffusion coefficient is [[Dx, Dy]]. -- 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: fipy diffusion setup
Phil, just to follow up and complete the answer. I neglected to explain how to create a rank 2 CellVariable from two rank 0 CellVariables to preserve the dependency. The dot operator helps with that. Given two rank 0 CellVariables, v00 and v11, then to get the rank 2 CellVariable use, vv = v00.dot([[1, 0], [0, 0]]) + v11.dot([[0, 0], [0, 1]]) On Mon, Apr 18, 2016 at 3:33 PM, Daniel Wheeler wrote: > On Mon, Apr 18, 2016 at 12:33 PM, Phil Battle wrote: >> Ok, below works as I would expect ( note Dx and Dz are constants) > > Thanks for that. I can confirm that it's broken in Python 2.7 as well > when the diffusion coefficient is [[Dx, Dy]]. > > > -- > Daniel Wheeler -- 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: fipy diffusion setup
On Mon, Apr 18, 2016 at 12:33 PM, Phil Battle wrote: > Ok, below works as I would expect ( note Dx and Dz are constants) Thanks for that. I can confirm that it's broken in Python 2.7 as well when the diffusion coefficient is [[Dx, Dy]]. -- 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: fipy diffusion setup
Ok, below works as I would expect ( note Dx and Dz are constants)
"
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from fipy import *
nx = 50.
ny = 50.
dx = .4
dy = .4
L = dx*nx
mesh = Grid2D(ny=ny,dy=dy, dx=dx,nx=nx)
w=10 #waveguide width um
dpe = 3 #initial Proton ex depth
dxpe = .098 # pe diffusion
dzpe=.056 #um^2/h
umax=0.2*dpe*(dxpe/dzpe)**.5
C = CellVariable(name = "solution variable",mesh = mesh, value =
0.,hasOld=1)
X,Z=mesh.cellCenters
wdpe = (Z<=dpe) & ((L/2-w/2)< X) & (X < (w/2 +L/2 ))
C.setValue(1,where = wdpe)
viewer = Viewer(vars=C, datamin=0., datamax=1.1)
Dx=1.
Dz=5.
eqX = TransientTerm() == DiffusionTerm([[[Dx,0],[0,Dz]]])
timeStepDuration = 2
steps = 5
time=2
while time<=10:
C.updateOld()
input('')
viewer.plot()
for step in range(steps):
print(eqX.sweep(var=C,dt=timeStepDuration))
#input('d')
time+=2
"
But when
[[[Dx,0],[0,Dz]]] replaced by [[Dx,Dz]]
is has the following error
---
AssertionError Traceback (most recent call last)
in() 42 viewer.plot()43 for step in range(steps):--->
44print(eqX.sweep(var=C,dt=timeStepDuration))45 #input('d')46
C:\Users\philb\AppData\Local\Continuum\Anaconda3\lib\site-packages\fipy-unknown-py3.5.egg\fipy\terms\term.py
in sweep(self, var, solver, boundaryConditions, dt, underRelaxation,
residualFn, cacheResidual, cacheError) 221 222 """ --> 223solver =
self._prepareLinearSystem(var=var, solver=solver,
boundaryConditions=boundaryConditions, dt=dt)224
solver._applyUnderRelaxation(underRelaxation=underRelaxation)225
residual =
solver._calcResidual(residualFn=residualFn)C:\Users\philb\AppData\Local\Continuum\Anaconda3\lib\site-packages\fipy-unknown-py3.5.egg\fipy\terms\term.py
in _prepareLinearSystem(self, var, solver, boundaryConditions, dt) 155
transientGeomCoeff=self._getTransientGeomCoeff(var),156
diffusionGeomCoeff=self._getDiffusionGeomCoeff(var),-->
157buildExplicitIfOther=self._buildExplcitIfOther) 158 159
self._buildCache(matrix,
RHSvector)C:\Users\philb\AppData\Local\Continuum\Anaconda3\lib\site-packages\fipy-unknown-py3.5.egg\fipy\terms\binaryTerm.py
in _buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions, dt,
transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther) 66
transientGeomCoeff=transientGeomCoeff,67
diffusionGeomCoeff=diffusionGeomCoeff,--->
68buildExplicitIfOther=buildExplicitIfOther) 69 70 matrix +=
tmpMatrixC:\Users\philb\AppData\Local\Continuum\Anaconda3\lib\site-packages\fipy-unknown-py3.5.egg\fipy\terms\unaryTerm.py
in _buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions, dt,
transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther) 97
dt=dt,98 transientGeomCoeff=transientGeomCoeff,--->
99diffusionGeomCoeff=diffusionGeomCoeff) 100 elif
buildExplicitIfOther:101 _, matrix, RHSvector =
self._buildMatrix(self.var,
C:\Users\philb\AppData\Local\Continuum\Anaconda3\lib\site-packages\fipy-unknown-py3.5.egg\fipy\terms\abstractDiffusionTerm.py
in _buildMatrix(self, var, SparseMatrix, boundaryConditions, dt,
transientGeomCoeff, diffusionGeomCoeff) 351 352 ids =
self._reshapeIDs(var, numerix.arange(mesh.numberOfCells))-->
353L.addAt(self.constraintL.ravel(), ids.ravel(),
ids.swapaxes(0,1).ravel())354 b +=
numerix.reshape(self.constraintB.ravel(), ids.shape).sum(-2).ravel()355
C:\Users\philb\AppData\Local\Continuum\Anaconda3\lib\site-packages\fipy-unknown-py3.5.egg\fipy\matrices\scipyMatrix.py
in addAt(self, vector, id1, id2) 251 2.50 --- 2.20252 """ -->
253assert(len(id1) == len(id2) == len(vector))254 255 temp =
sp.csr_matrix((vector, (id1, id2)), self.matrix.shape)AssertionError:
On 4/18/2016 10:04 AM, Daniel Wheeler wrote:
On Mon, Apr 18, 2016 at 11:49 AM, Phil Battle wrote:
is it possible the conversion to python 3 is the issue with the
equivalency in the case of a constant anisotropic diffusion coeff. ?
I was wondering about that. Can you post a snippet (enough to demo the
issue) and I'll try it in Python 2.7 if you don't have FiPy working in
2.7.
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: fipy diffusion setup
On Mon, Apr 18, 2016 at 11:49 AM, Phil Battle wrote: > > is it possible the conversion to python 3 is the issue with the > equivalency in the case of a constant anisotropic diffusion coeff. ? I was wondering about that. Can you post a snippet (enough to demo the issue) and I'll try it in Python 2.7 if you don't have FiPy working in 2.7. -- 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: fipy diffusion setup
thank you for the response - I will try your suggestion on the anisotropic diffusion coeff. is it possible the conversion to python 3 is the issue with the equivalency in the case of a constant anisotropic diffusion coeff. ? My feeling is that [[1., 5.]] shoud be equivalent to [[[1., 0.], [0., 5.]]]. I'm not sure why that isn't working correctly. On 4/18/2016 8:54 AM, Daniel Wheeler wrote: > On Fri, Apr 15, 2016 at 9:29 AM, Phil wrote: >> But when written like >> >> TransientTerm() == DiffusionTerm([[1.,5.]]), which I think the >> manual said was equivalent ( under freq asked questions) >> I get an assertion error, which I can post but before posting should >> these two diffusion terms be equivalent? > My feeling is that [[1., 5.]] shoud be equivalent to [[[1., 0.], [0., > 5.]]]. I'm not sure why that isn't working correctly. > >> Sticking with the complete matrix form of the diffusion Coefficient >> for Dx and Dz ( both similar functions of the cellVariable C) >> >> when I use this form >> >> TransientTerm() == DiffusionTerm([[[Dx,0],[0,Dz]]]) >> >> I get the following traceback error > That's not surprising. Use a single FaceVariable or CellVariable of > rank 2 for the diffusion coefficient. It's safer to do that for > anisotropy. > >> I see questions posted on this list, Gist and stack overflow - is there a >> preference for posting ? > I prefer Stackoverflow, but depends on the nature of the question. > Stackoverflow isn't good for very open ended questions or > conversations. > ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: fipy diffusion setup
> On Apr 18, 2016, at 10:54 AM, Daniel Wheeler > wrote: > > On Fri, Apr 15, 2016 at 9:29 AM, Phil wrote: >> I see questions posted on this list, Gist and stack overflow - is there a >> preference for posting ? > > I prefer Stackoverflow, but depends on the nature of the question. > Stackoverflow isn't good for very open ended questions or > conversations. We monitor both our own list (obviously!) and StackOverflow. There's really no way to monitor gist, but it's a good place to put code snippets referenced from the list or SO. As Daniel says, SO isn't very good for open-ended conversations. On the other hand, it's easier to search and has a lower barrier to entry (you don't need to subscribe) than the mailing list. Really, it's whatever you feel most comfortable with. ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
Re: fipy diffusion setup
On Fri, Apr 15, 2016 at 9:29 AM, Phil wrote: > > But when written like > > TransientTerm() == DiffusionTerm([[1.,5.]]), which I think the > manual said was equivalent ( under freq asked questions) > I get an assertion error, which I can post but before posting should > these two diffusion terms be equivalent? My feeling is that [[1., 5.]] shoud be equivalent to [[[1., 0.], [0., 5.]]]. I'm not sure why that isn't working correctly. > Sticking with the complete matrix form of the diffusion Coefficient > for Dx and Dz ( both similar functions of the cellVariable C) > > when I use this form > > TransientTerm() == DiffusionTerm([[[Dx,0],[0,Dz]]]) > > I get the following traceback error That's not surprising. Use a single FaceVariable or CellVariable of rank 2 for the diffusion coefficient. It's safer to do that for anisotropy. > I see questions posted on this list, Gist and stack overflow - is there a > preference for posting ? I prefer Stackoverflow, but depends on the nature of the question. Stackoverflow isn't good for very open ended questions or conversations. -- Daniel Wheeler ___ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
fipy diffusion setup
Hello
I am am using FiPy with Python 3.4 - I think the conversion went
well enough but I do have some questions.
I have successfully implemented a 1-dim concentration dependent
diffusion problem in fipy, where success was defined by reproducing
what other researchers had produced in a publication.
My diffusion coefficient looks like
Dx=Dpe_x*(a_x+((1.-a_x)/(b_x*C+g_x)))
where C is a CellVariable on a 1 dim grid and the other terms are
constants.
As mentioned this all worked well. However when I went to extend to
2 Dimension and include second diffusion coeff, Dz (C) , for
anisotropic diffusion, I am running into some issues.
first for constant diffusion coeff of Dx =1, Dz = 5 on 2 dim
rectangular grid (Grid2D)
When written as
TransientTerm() ==DiffusionTerm([[[1.,0],[0,5.]]])
the results behave like I would expect ( on a 2dim rectangular mesh)
But when written like
TransientTerm() == DiffusionTerm([[1.,5.]]), which I think the
manual said was equivalent ( under freq asked questions)
I get an assertion error, which I can post but before posting should
these two diffusion terms be equivalent?
Sticking with the complete matrix form of the diffusion Coefficient
for Dx and Dz ( both similar functions of the cellVariable C)
when I use this form
TransientTerm() == DiffusionTerm([[[Dx,0],[0,Dz]]])
I get the following traceback error
---
ValueError Traceback (most recent call last)
in() 24
Dz=Dpe_z*(a_z+((1.-a_z)/(b_z*C+g_z)))
25
---> 26eqX = TransientTerm() ==
DiffusionTerm([[[Dx,0],[0,Dz]]])C:\Users\PhilB\Anaconda3\lib\site-packages\fipy-unknown-py3.4.egg\fipy\terms\abstractDiffusionTerm.py
in __init__(self, coeff, var) 62 from fipy.variables.variable import
Variable63 if not isinstance(self.nthCoeff, Variable):--->
64self.nthCoeff = Variable(value = self.nthCoeff)65 66 from
fipy.variables.cellVariable import
CellVariableC:\Users\PhilB\Anaconda3\lib\site-packages\fipy-unknown-py3.4.egg\fipy\variables\variable.py
in __init__(self, value, unit, array, name, cached) 109 unit = None110
--> 111self._setValueInternal(value=value, unit=unit, array=array)112
113 self._name =
nameC:\Users\PhilB\Anaconda3\lib\site-packages\fipy-unknown-py3.4.egg\fipy\variables\variable.py
in _setValueInternal(self, value, unit, array) 638 639 def
_setValueInternal(self, value, unit=None, array=None):--> 640self._value
= self._makeValue(value=value, unit=unit, array=array)641 642 def
_makeValue(self, value, unit=None,
array=None):C:\Users\PhilB\Anaconda3\lib\site-packages\fipy-unknown-py3.4.egg\fipy\variables\variable.py
in _makeValue(self, value, unit, array) 668 669 if unit is not None or
type(value) in [type(''), type(()), type([])]:--> 670value =
PF(value=value, unit=unit, array=array)671 elif array is not None:672
array[:] =
valueC:\Users\PhilB\Anaconda3\lib\site-packages\fipy-unknown-py3.4.egg\fipy\tools\dimensions\physicalField.py
in __init__(self, value, unit, array) 203 204 if type(value) in
[type([]),type< span class="ansiyellow" style="box-sizing: border-box;
color: rgb(196, 160, 0);">(())]:--> 205value = [PhysicalField(item,unit)
for item in value]206 if unit is None:207 unit =
value[0].unitC:\Users\PhilB\Anaconda3\lib\site-packages\fipy-unknown-py3.4.egg\fipy\tools\dimensions\physicalField.py
in (.0) 203 204 if type(value) in [type([]),type< span
class="ansiyellow" style="box-sizing: border-box; color: rgb(196, 160,
0);">(())]:--> 205value = [PhysicalField(item,unit) for item in
value]206 if unit is None:207 unit =
value[0].unitC:\Users\PhilB\Anaconda3\lib\site-packages\fipy-unknown-py3.4.egg\fipy\tools\dimensions\physicalField.py
in __init__(self, value, unit, array) 212 else:213 normalized +=
[item.inUnitsOf(unit).value]--> 214value = numerix.array(normalized)215
216 if unit is None:ValueError: setting an array element with a sequence.
Is there a more appropriate way of writing the 2Dim anisotropic
concentration dependent diff coeff?
Note if I use isotropic 2 dim concentration dependent diffusion the
results look correct ( in 2D)
TransientTerm() == DiffusionTerm(Dx)
Let me know if this challenge makes sense and if there is a solution
thanks
Phil
ps
I see questions posted on this list, Gist and stack overflow - is there
a preference for posting ?
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
