Hi, thanks for this!
I've updated but it seems that the program can not pass line 267-268 for phase 1 (phase 0 was ok). Is this because of the convergence problem or possibly some other reasons?
Thanks and regards,
hzhatlboro
I've updated but it seems that the program can not pass line 267-268 for phase 1 (phase 0 was ok). Is this because of the convergence problem or possibly some other reasons?
Thanks and regards,
hzhatlboro
----- Original Message ----
From: Daniel Wheeler <[EMAIL PROTECTED]>
To: Multiple recipients of list <[email protected]>
Sent: Friday, 13 October, 2006 2:59:30 PM
Subject: Re: solvers in Fipy
From: Daniel Wheeler <[EMAIL PROTECTED]>
To: Multiple recipients of list <[email protected]>
Sent: Friday, 13 October, 2006 2:59:30 PM
Subject: Re: solvers in Fipy
On Oct 12, 2006, at 5:45 PM, Daniel Wheeler wrote:
You have "2 / p" in a number of your coefficients. Replace that with "2 / (p==0 * 1e-20 + p)". That will stop the divide by 0 error.
Apologies. The above should be "2 / ((p == 0) * 1e-20 + p)".
--
Daniel Wheeler
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
cViewer=viewers.make(vars=c, limits={'datamin':0., 'datamax':1.})
phaseViewer=viewers.make(vars=phases, limits={'datamin':0., 'datamax':1.})
cViewer.plot()
phaseViewer.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 i in range(steps):
print "Step ", i+1
c.updateOld()
cRes=1e+10
print "Solving equation of c..."
while cRes > 1e-3:
cRes=c.equation.sweep(var=c, dt=dt)
print "Equation of c is solved..."
for phii in phases:
print "Solving equation of phase ", counti, "..."
for j in range(numPhases):
s[j]=(phases[j]>0)
#print "s ",j,"=",s[j]
p=0
for pj in s:
p+=pj
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==0) * 1e-20 + p))*s[counti]*s[j]*M[counti,j]*epsilon[counti,j]**2/2
S0+=-(2 / ((p==0) * 1e-20 + p))*s[counti]*s[j]*M[counti,j]*.5*epsilon[counti,j]**2* \
(phases[j].getFaceGrad()).getDivergence()
S1+=-(2 / ((p==0) * 1e-20 + p))*s[counti]*s[j]*M[counti,j]*omega[counti,j]*(phases[j])
imScoeff+=(2 / ((p==0) * 1e-20 + 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==0) * 1e-20 + 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
raw_input("Formulation of phase equation done. Press <Return> to continue...")
phii.updateOld()
phaseRes=1e+10
while phaseRes > 1e-3:
phaseRes=phii.equation.sweep(var=phii, dt=dt)
print "phase ", counti, "equation is solved..."
counti=counti+1
if __name__=='__main__':
from fipy import viewers
cViewer=viewers.make(vars=c,limits={'datamin':0.,'datamax':1.})
phaseViewer=viewers.make(vars=phases,limits={'datamin':0.,'datamax':1.})
cViewer.plot()
phaseViewer.plot()
raw_input("Press <Return> to continue...")
