Sorry for the confusion. For some reasons, it did not show me that error,
but I fixed the Pf definition. Here is my code which still give me the
index error.
CODE
--------------------------------------------------------
from __future__ import division
from numpy import *
from pylab import *
from matplotlib import rc
from fipy import *
#################### INITIALLIZATION
###########################################
L = 0.0508 #+0.0021+0.0014
nx = 10.0
dx = L / nx
Pdi=6894757.29
Pui=7170547.58
D=1e-6
tmax = 1000
sigma = 400
k0 = 10e-16
T = 300
M = 0.0399
R = 8.314
muf=1.96e-5
hf=30e-6
phif=0.25
A1=1.14
A2=0.28
wkf=4
wkd=4
wku=4
mesh = Grid1D(nx=nx, dx=dx)
############################### VARIABLE DEFINITIN
###################################
Pk = CellVariable(mesh=mesh, value=Pdi,hasOld=True)
Pf = CellVariable(mesh=mesh, value=Pdi,hasOld=True)
Pd = CellVariable(mesh=mesh, value=Pdi,hasOld=True)
Pu = CellVariable(mesh=mesh, value=Pui,hasOld=True)
Kn=sqrt((3.14*R*T/(2*M))*muf/(Pf*hf))
############################# Kerogen mass balance
Pk.constrain([Pd], mesh.facesLeft)
Pk.constrain([Pu], mesh.facesRight)
eqpk = TransientTerm(var=Pk) == DiffusionTerm(D, var=Pk)
############################# Natural Fracture mass balance
Df = k0 * (1 + 6 * A1 * Kn + 12 * A2 * Kn**2)*Pf/muf
wkfD=wkf*D
Pf.constrain([Pd], mesh.facesLeft)
Pf.constrain([Pu], mesh.facesRight)
eqpf = TransientTerm(phif,var=Pf) == DiffusionTerm(Df,var=Pf) -
ImplicitSourceTerm(coeff=wkfD,var=Pf) +
ImplicitSourceTerm(coeff=wkfD,var=Pk)
################################ Downstream tank
sigmaD = (k0*(1 + 6 * A1 * Kn + 12 * A2 * Kn**2)*sigma)/(muf)
eqpd = TransientTerm(var=Pd) == ImplicitSourceTerm(coeff=sigmaD,var=Pf) -
ImplicitSourceTerm(coeff=sigmaD,var=Pd)-
ImplicitSourceTerm(coeff=wkfD,var=Pd) +
ImplicitSourceTerm(coeff=wkfD,var=Pk)
############################## Upstream tank
eqpu = - TransientTerm(var=Pu) == ImplicitSourceTerm(coeff=sigmaD,var=Pf) -
ImplicitSourceTerm(coeff=sigmaD,var=Pu)-
ImplicitSourceTerm(coeff=wkfD,var=Pu) +
ImplicitSourceTerm(coeff=wkfD,var=Pk)
########################## TIMESTEP
dt = 0.9 * dx**2 / (2 * D)
steps = 100
######################### Coupling
eqn = eqpk & eqpf & eqpd & eqpu
######################## Solver
for t in range(1):
Pk.updateOld()
Pf.updateOld()
Pd.updateOld()
Pu.updateOld()
eqn.sweep(dt=dt)
----------------------------------------------------------------------
Error:
----------------------
--------------------------------------------------------------------------IndexError
Traceback (most recent call
last)<ipython-input-2-02f95fecdd85> in <module>() 72
Pd.updateOld() 73 Pu.updateOld()---> 74 eqn.sweep(dt=dt)
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/term.pyc
in sweep(self, var, solver, boundaryConditions, dt, underRelaxation,
residualFn, cacheResidual, cacheError) 234 235 """-->
236 solver = self._prepareLinearSystem(var=var, solver=solver,
boundaryConditions=boundaryConditions, dt=dt) 237
solver._applyUnderRelaxation(underRelaxation=underRelaxation) 238
residual = solver._calcResidual(residualFn=residualFn)
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/term.pyc
in _prepareLinearSystem(self, var, solver, boundaryConditions, dt)
168
transientGeomCoeff=self._getTransientGeomCoeff(var), 169
diffusionGeomCoeff=self._getDiffusionGeomCoeff(var),--> 170
buildExplicitIfOther=self._buildExplcitIfOther) 171 172
self._buildCache(matrix, RHSvector)
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/coupledBinaryTerm.pyc
in _buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions,
dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
120
transientGeomCoeff=uncoupledTerm._getTransientGeomCoeff(tmpVar),
121
diffusionGeomCoeff=uncoupledTerm._getDiffusionGeomCoeff(tmpVar),-->
122
buildExplicitIfOther=buildExplicitIfOther) 123
124 termMatrix += tmpMatrix
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/binaryTerm.pyc
in _buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions,
dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
66
transientGeomCoeff=transientGeomCoeff, 67
diffusionGeomCoeff=diffusionGeomCoeff,---> 68
buildExplicitIfOther=buildExplicitIfOther) 69 70
matrix += tmpMatrix
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/unaryTerm.pyc
in _buildAndAddMatrices(self, var, SparseMatrix, boundaryConditions,
dt, transientGeomCoeff, diffusionGeomCoeff, buildExplicitIfOther)
97 dt=dt,
98
transientGeomCoeff=transientGeomCoeff,---> 99
diffusionGeomCoeff=diffusionGeomCoeff)
100 elif buildExplicitIfOther: 101 _, matrix,
RHSvector = self._buildMatrix(self.var,
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/terms/abstractDiffusionTerm.pyc
in _buildMatrix(self, var, SparseMatrix, boundaryConditions, dt,
transientGeomCoeff, diffusionGeomCoeff) 335
coeff = self.nthCoeff 336 --> 337
nthCoeffFaceGrad = coeff[numerix.newaxis] *
var.faceGrad[:,numerix.newaxis] 338 s =
(slice(0,None,None),) + (numerix.newaxis,) * (len(coeff.shape) - 1) +
(slice(0,None,None),) 339 normalsNthCoeff =
coeff[numerix.newaxis] * normals[s]
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/variables/variable.pyc
in __getitem__(self, index) 1573
operatorClass=self._getitemClass(index=index), 1574
opShape=numerix._indexShape(index=index, arrayShape=self.shape),->
1575 unit=self.unit, 1576
canInline=False) 1577
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/variables/variable.pyc
in _getUnit(self) 253 <PhysicalUnit m> 254
"""--> 255 return self._extractUnit(self.value) 256 257
@getsetDeprecated
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/variables/variable.pyc
in _getValue(self) 536 537 if self.stale or not
self._isCached() or self._value is None:--> 538 value =
self._calcValue() 539 if self._isCached(): 540
self._setValueInternal(value=value)
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/variables/faceGradVariable.pyc
in _calcValue(self) 80 return self._calcValueInline()
81 else:---> 82 return
self._calcValueNoInline() 83 84 def
_calcValueInline(self):
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/variables/faceGradVariable.pyc
in _calcValueNoInline(self) 154 s = (faceMask,) 155
--> 156 N2[s] = self.var.faceValue[s] 157 158 N
= (N2 - numerix.take(self.var,id1, axis=-1)) / dAP
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/variables/variable.pyc
in __getitem__(self, index) 1573
operatorClass=self._getitemClass(index=index), 1574
opShape=numerix._indexShape(index=index, arrayShape=self.shape),->
1575 unit=self.unit, 1576
canInline=False) 1577
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/variables/variable.pyc
in _getUnit(self) 253 <PhysicalUnit m> 254
"""--> 255 return self._extractUnit(self.value) 256 257
@getsetDeprecated
/home/mkhome/anaconda/lib/python2.7/site-packages/fipy/variables/variable.pyc
in _getValue(self) 559 value[..., mask]
= constraint.value 560 except:--> 561
value[..., mask] =
numerix.array(constraint.value)[..., mask] 562 563
return value
IndexError: index 1000 is out of bounds for axis 1 with size 1000
On Fri, May 6, 2016 at 1:17 PM, Raymond Smith <[email protected]> wrote:
> I think Daniel's point is that Pf is not defined yet when it is first
> called in the script. You could fix this by moving the line using Pf that
> Daniel pointed to after the definition that you refer to. I think the
> general point, however, is that it makes it much easier to help if you try
> to make sure that the script you send gives the errors that you see.
>
> Best,
> Ray
>
>
> <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=webmail&utm_term=icon>
> Virus-free.
> www.avast.com
> <https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=webmail&utm_term=link>
> <#m_8816831225450020302_m_-68387843283325897_DDB4FAA8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
>
> On Fri, May 6, 2016 at 1:00 PM, Mohammad Kazemi <[email protected]>
> wrote:
>
>> Thanks Daniel for your response. I thought defining variables as
>> cellvariable should be enough.
>>
>> Pf = CellVariable(mesh=mesh, value=Pdi,hasOld=True)
>>
>>
>> Should I define it in some other way?
>>
>> Thanks
>>
>> On Fri, May 6, 2016 at 12:10 PM, Daniel Wheeler <
>> [email protected]> wrote:
>>
>>> On Fri, May 6, 2016 at 12:09 PM, Daniel Wheeler
>>> <[email protected]> wrote:
>>> > Hi Mo,
>>>
>>> Apologies for that.
>>>
>>> --
>>> Daniel Wheeler
>>> _______________________________________________
>>> fipy mailing list
>>> [email protected]
>>> http://www.ctcms.nist.gov/fipy
>>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>>>
>>
>>
>>
>> --
>> Mohammad Kazemi
>> West Virginia University
>> Department of Petroleum and Natural Gas Engineering
>>
>> _______________________________________________
>> fipy mailing list
>> [email protected]
>> http://www.ctcms.nist.gov/fipy
>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>>
>>
>
> _______________________________________________
> fipy mailing list
> [email protected]
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>
--
Mohammad Kazemi
West Virginia University
Department of Petroleum and Natural Gas Engineering
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
