Dear Fipy users and developers,
I was trying to solve a system of partial differential equations using fipy
but I am keep getting error of "" . I would appreciate if you could help me
on this error. Here is my code and 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
Kn=sqrt((3.14*R*T/(2*M))*muf/(Pf*hf))
mesh = Grid1D(nx=nx, dx=dx)
X = mesh.faceCenters[0]
#meshd = Grid1D(nx=nx, dx=dx)
############################### VARIABLE DEFINITION
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)
############################# 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-55-c706252968c2> in <module>() 73
Pd.updateOld() 74 Pu.updateOld()---> 75 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 10 is out of bounds for axis 1 with size 10
---------------------------end
--
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 ]
