Hi, Daniel and Jon,
Thanks for your useful support and I've now reached a stage to choose a proper
solver for my problem (there might be some hidden errors, though). I did find
that it is much easier to formulate my equations by using FiPy.
Attached please find the program and I have tried to use "sweep", which is a
concept not quite familiar to me, to solve the equations. However, I got the
following error. I would be grateful if you could advise what kind of solvers
are available in FiPy and which one should be used for which circumstance. In
particular, which solver I should choose for my problem. I apologize if this
info was in the manual already.
Thanks and regards,
hzhatlboro
Traceback (most recent call last):
File "C:\Python24\hzh\pf2dtest.py", line 251, in ?
phaseRes=phii.equation.sweep(var=phii, dt=dt)
File "C:\Python24\Lib\site-packages\fipy\terms\term.py", line 149, in sweep
matrix, RHSvector = self.__buildMatrix(var, boundaryConditions, dt)
File "C:\Python24\Lib\site-packages\fipy\terms\term.py", line 95, in
__buildMatrix
return self._buildMatrix(var, boundaryConditions, dt)
File "C:\Python24\Lib\site-packages\fipy\terms\binaryTerm.py", line 66, in
_buildMatrix
termMatrix, termRHSvector = self.term2._buildMatrix(var,
boundaryConditions, dt = dt)
File "C:\Python24\Lib\site-packages\fipy\terms\binaryTerm.py", line 64, in
_buildMatrix
matrix, RHSvector = self.term1._buildMatrix(var, boundaryConditions, dt =
dt)
File "C:\Python24\Lib\site-packages\fipy\terms\binaryTerm.py", line 64, in
_buildMatrix
matrix, RHSvector = self.term1._buildMatrix(var, boundaryConditions, dt =
dt)
File "C:\Python24\Lib\site-packages\fipy\terms\binaryTerm.py", line 64, in
_buildMatrix
matrix, RHSvector = self.term1._buildMatrix(var, boundaryConditions, dt =
dt)
File "C:\Python24\Lib\site-packages\fipy\terms\binaryTerm.py", line 64, in
_buildMatrix
matrix, RHSvector = self.term1._buildMatrix(var, boundaryConditions, dt =
dt)
File "C:\Python24\Lib\site-packages\fipy\terms\binaryTerm.py", line 66, in
_buildMatrix
termMatrix, termRHSvector = self.term2._buildMatrix(var,
boundaryConditions, dt = dt)
File "C:\Python24\Lib\site-packages\fipy\terms\cellTerm.py", line 121, in
_buildMatrix
inline._optionalInline(self._buildMatrixIn, self._buildMatrixPy, L,
var.getOld(), b, dt, coeffVectors)
File "C:\Python24\Lib\site-packages\fipy\tools\inline\inline.py", line 10, in
_optionalInline
return pythonFn(*args)
File "C:\Python24\Lib\site-packages\fipy\terms\cellTerm.py", line 88, in
_buildMatrixPy
b += Numeric.array(oldArray) * Numeric.array(coeffVectors['old value']) / dt
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 162, in
__array__
return Numeric.array(self.getValue(), t)
File "C:\Python24\Lib\site-packages\fipy\variables\cellVariable.py", line
108, in getValue
return Variable.getValue(self)
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 342, in
getValue
value = self._calcValue()
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 1708,
in _calcValue
return self.op(self.var[0].getValue(), val1)
File "C:\Python24\Lib\site-packages\fipy\variables\cellVariable.py", line
108, in getValue
return Variable.getValue(self)
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 342, in
getValue
value = self._calcValue()
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 1708,
in _calcValue
return self.op(self.var[0].getValue(), val1)
File "C:\Python24\Lib\site-packages\fipy\variables\cellVariable.py", line
108, in getValue
return Variable.getValue(self)
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 342, in
getValue
value = self._calcValue()
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 1702,
in _calcValue
val1 = self.var[1].getValue()
File "C:\Python24\Lib\site-packages\fipy\variables\cellVariable.py", line
108, in getValue
return Variable.getValue(self)
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 342, in
getValue
value = self._calcValue()
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 1708,
in _calcValue
return self.op(self.var[0].getValue(), val1)
File "C:\Python24\Lib\site-packages\fipy\variables\cellVariable.py", line
108, in getValue
return Variable.getValue(self)
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 342, in
getValue
value = self._calcValue()
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 1708,
in _calcValue
return self.op(self.var[0].getValue(), val1)
File "C:\Python24\Lib\site-packages\fipy\variables\cellVariable.py", line
108, in getValue
return Variable.getValue(self)
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 342, in
getValue
value = self._calcValue()
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 1708,
in _calcValue
return self.op(self.var[0].getValue(), val1)
File "C:\Python24\Lib\site-packages\fipy\variables\cellVariable.py", line
108, in getValue
return Variable.getValue(self)
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 342, in
getValue
value = self._calcValue()
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 1708,
in _calcValue
return self.op(self.var[0].getValue(), val1)
File "C:\Python24\Lib\site-packages\fipy\variables\cellVariable.py", line
108, in getValue
return Variable.getValue(self)
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 342, in
getValue
value = self._calcValue()
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 1708,
in _calcValue
return self.op(self.var[0].getValue(), val1)
File "C:\Python24\Lib\site-packages\fipy\variables\variable.py", line 1753,
in <lambda>
return self._getBinaryOperatorVariable(lambda a,b: b/a, other)
ZeroDivisionError: divide by zero
Send instant messages to your online friends http://uk.messenger.yahoo.com numPhases=4
H=0.25
Calpha0=0.002
Cimc0=0.4105
Cliq0=0.9769
Wij=4*H
sigmaLiqEta=0.1
sigmaEtaAlpha=0.3
sigmaGb=sigmaEtaAlpha
Mgb=7e4
MetaAlpha=Mgb
MliqEta=1e6
Dalpha=1e-6
Deta=2e-5
DliqEta=2e-1
Dliq=1.
DetaAlpha=2e-3
Dgb=2e-2
dt=0.0125
steps=10
import scipy
from fipy.tools import numerix
# sigma & M
sigma=sigmaGb*numerix.ones((numPhases,numPhases))
M=Mgb*numerix.ones((numPhases,numPhases))
sigma-=sigmaGb*numerix.array(scipy.eye(numPhases,numPhases))
M-=Mgb*numerix.array(scipy.eye(numPhases,numPhases))
sigma[:,numPhases-1]=sigmaLiqEta
M[:,numPhases-1]=MliqEta
sigma[numPhases-1,:]=sigmaLiqEta
M[numPhases-1,:]=MliqEta
sigma[0,numPhases-1]=0
M[0,numPhases-1]=0
sigma[numPhases-1,0]=0
M[numPhases-1,0]=0
sigma[numPhases-1,numPhases-1]=0
M[numPhases-1,numPhases-1]=0
#epsilon
from fipy.tools.numerix import sqrt
from fipy.tools import numerix
epsilon=sqrt(8*sigma*4*H)/numerix.pi
#wij
omega=sigma/H
print "sigma = "
print sigma, '\n'
print "M = "
print M, '\n'
print "epsilon = "
print epsilon, '\n'
print "omega = "
print omega, '\n'
nx=16
ny=90
dx=H
dy=dx
from fipy.meshes.grid2D import Grid2D
mesh=Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
from fipy.variables.cellVariable import CellVariable
c=CellVariable(name="Concentration",
mesh=mesh,
value=Calpha0,
hasOld=1)
phases=[
CellVariable(mesh=mesh, name='phi0', value=0, hasOld=1),
CellVariable(mesh=mesh, name='phi1', value=0, hasOld=1),
CellVariable(mesh=mesh, name='phi2', value=0, hasOld=1),
CellVariable(mesh=mesh, name='phi3', value=0, hasOld=1)
]
x, y=mesh.getCellCenters()[...,0], mesh.getCellCenters()[...,1]
c.setValue(Cimc0, where=(41*H<y) & (y<45*H))
c.setValue(Cliq0, where=(y>45*H))
phases[0].setValue(1, where=(y<41*H))
phases[1].setValue(1, where=(41*H<y) & (y<45*H) & (x<8*H))
phases[2].setValue(1, where=(41*H<y) & (y<45*H) & (x>8*H))
phases[3].setValue(1, where=(y>45*H))
# save intial conditions to data files
#from fipy.viewers.tsvViewer import TSVViewer
#TSVViewer(vars=phases[3]).plot(filename="pf2dtestPhi3init.tsv")
from fipy import viewers
viewer=viewers.make(vars=c, limits={'datamin':0., 'datamax':1.})
viewer.plot()
from fipy.terms.transientTerm import TransientTerm
from fipy.terms.implicitDiffusionTerm import ImplicitDiffusionTerm
from fipy.terms.implicitSourceTerm import ImplicitSourceTerm
def Dphi(phases):
g1=(phases[0]>0.8)
g2=(phases[3]>0.8)
g3=((phases[1]>0.8) | (phases[2]>0.8))
g4=((phases[3]>=0.2) & (phases[3]<=0.8))
g5=(phases[3]<0.2)
g6=((((1-phases[1])/2)<phases[0]) & (phases[0]<=0.8))
g7=((((1-phases[2])/2)<phases[0]) & (phases[0]<=0.8))
DphiValue=Dgb*(1-g1-g2-g3-g4-g5*g6*g7)+g1*Dalpha+g2*Dliq+g3*Deta+g4*DliqEta+DetaAlpha*(g5*g6*g7)
return DphiValue
def sumPhases(phases):
sumPhasesValue=0.
for phii in phases:
sumPhasesValue+=phii
return sumPhasesValue
c.equation=TransientTerm()==ImplicitDiffusionTerm(coeff=(Dphi(phases)*sumPhases(phases)))
G0alphaCu=-1.9073e4
G0alphaSn=-2.7281e4
Lalpha0=-1.1448e4
Lalpha1=-1.1694e4
T=523.
R=8.314
Vm=16.29e-6
GserSn=3.4616e5
G0LCu=-1.1083e4
G0LSn=-2.8963e4
LL0=-1.0487e4
LL1=-1.8198e4
LL2=-1.05284e4
from fipy.tools.numerix import log
def falpha(c):
falphaValue=((1-c)*G0alphaCu+c*G0alphaSn+R*T*((1-c)*log(1-c)+c*log(c))+c*(1-c)*(Lalpha0+Lalpha1*(1-2*c)))/Vm
return falphaValue
def falpha1(c):
falpha1Value=(-G0alphaCu+G0alphaSn+R*T*(-log(1-c)+log(c))+(1-c)*(Lalpha0+Lalpha1*(1-2*c))-c*(Lalpha0+Lalpha1* \
(1-2*c))-2*c*(1-c)*Lalpha1)/Vm
return falpha1Value
def feta(c):
fetaValue=(2.e5*(c-0.435)**2+0.545*G0alphaCu+0.455*GserSn-6869.5-0.1589*T)/Vm
return fetaValue
def feta1(c):
feta1Value=(4e5*c-1.74e5)/Vm
return feta1Value
def fliq(c):
fliqValue=((1-c)*G0LCu+c*G0LSn+R*T*((1-c)*log(1-c)+c*log(c))+c*(1-c)*(LL0+LL1*(1-2*c)+LL2*(1-4*c+4*c**2)))/Vm
return fliqValue
def fliq1(c):
fliq1Value=(-G0LCu+G0LSn+R*T*(-log(1-c)+log(c))+(1-c)*(LL0+LL1*(1-2*c)+LL2*(1-4*c+4*c**2))-c*(LL0+LL1*(1-2*c)+ \
LL2*(1-4*c+4*c**2))+c*(1-c)*(-2*LL1+LL2*(-4+8*c)))/Vm
return fliq1Value
counti=0;
s=[
CellVariable(mesh=mesh, name='s0', value=(phases[0]>0)),
CellVariable(mesh=mesh, name='s1', value=(phases[1]>0)),
CellVariable(mesh=mesh, name='s2', value=(phases[2]>0)),
CellVariable(mesh=mesh, name='s3', value=(phases[3]>0))
]
for phii in phases:
for j in range(numPhases):
s[j]=(phases[j]>0)
print "s ",j,"=",s[j]
p=0
for pj in s:
p+=pj
print "p = ",p
diffCoe=0
S0=0
S1=0
S20=0
S2=0
S3=0
imScoeff=0
for j in range(numPhases):
if j!=counti:
diffCoe+=(2/p)*s[counti]*s[j]*M[counti,j]*epsilon[counti,j]**2/2
S0+=-(2/p)*s[counti]*s[j]*M[counti,j]*.5*epsilon[counti,j]**2*(phases[j].getFaceGrad()).getDivergence()
S1+=-(2/p)*s[counti]*s[j]*M[counti,j]*omega[counti,j]*(phases[j])
imScoeff+=(2/p)*s[counti]*s[j]*M[counti,j]*omega[counti,j]
for k in range(numPhases):
if k!=counti & k!=j:
S20+=s[k]*((.5*epsilon[counti,k]**2-.5*epsilon[j,k]**2)*(phases[k].getFaceGrad()).getDivergence() \
+(omega[counti,k]-omega[j,k])*phases[k])
S2+=-(2/p)*s[counti]*s[j]*M[counti,j]*S20
if counti==0:
ChemEnergyI=falpha(c)
ChemEnergyI1=falpha1(c)
elif counti==numPhases-1:
ChemEnergyI=fliq(c)
ChemEnergyI1=fliq1(c)
else:
ChemEnergyI=feta(c)
ChemEnergyI1=feta1(c)
if j==0:
ChemEnergyJ=falpha(c)
ChemEnergyJ1=falpha1(c)
elif j==numPhases-1:
ChemEnergyJ=fliq(c)
ChemEnergyJ1=fliq1(c)
else:
ChemEnergyJ=feta(c)
ChemEnergyJ1=feta1(c)
S3+=(ChemEnergyI-c*ChemEnergyI1)-(ChemEnergyJ-c*ChemEnergyJ1)
phii.equation=TransientTerm()==ImplicitDiffusionTerm(coeff=(diffCoe))+ImplicitSourceTerm(coeff=imScoeff)+S0+S1+S2+S3
for i in range(steps):
c.updateOld()
phii.updateOld()
cRes=1e+10
phaseRes=1e+10
while cRes > 1e-3 or phaseRes > 1e-3:
cRes=c.equation.sweep(var=c, dt=dt)
phaseRes=phii.equation.sweep(var=phii, dt=dt)
if __name__=='__main__':
viewer.plot()
counti+=1
