#example from : http://www.ctcms.nist.gov/fipy/examples/levelSet/generated/examples.levelSet.electroChem.howToWriteAScript.html#module-examples.levelSet.electroChem.howToWriteAScript from fipy import * from pylab import * from matplotlib.pyplot import * from mpl_toolkits.axes_grid.axislines import Subplot import sys import os from fipy.terms.implicitSourceTerm import ImplicitSourceTerm from fipy.models.levelSet.distanceFunction.levelSetDiffusionEquation import _buildLevelSetDiffusionEquation from fipy.tools import numerix from fipy.variables.cellVariable import CellVariable #parameters fparam=open('parameters.input','r') params=fparam.readline().split() cellSize=float(params[0]) nx = float(params[1]) ny = float(params[2]) dx = dy = cellSize CX=nx*cellSize/2.0 CY=ny*cellSize/2.0 Rstart=cellSize*float(params[3]) D=float(params[4]) v0=float(params[5]) numberOfsteps=int(params[6]) dt0 = float(params[7]) everynimages=int (params[8]) __docformat__ = 'restructuredtext' class _CUSTOMMetalIonSourceVariable(CellVariable): def __init__(self, ionVar = None, distanceVar = None, depositionRate = None, metalIonMolarVolume = None): CellVariable.__init__(self, distanceVar.mesh, hasOld = 0) self.ionVar = self._requires(ionVar) self.distanceVar = self._requires(distanceVar) self.depositionRate = self._requires(depositionRate) self.metalIonMolarVolume = metalIonMolarVolume def _calcValue(self): ionVar = numerix.array(self.ionVar) ionVar = numerix.where(ionVar > 1e-20, ionVar, 1e-20) #return numerix.array(self.depositionRate) * self.distanceVar.cellInterfaceAreas / (self.mesh.cellVolumes * self.metalIonMolarVolume) / ionVar #changed by the following line to get the boundary condition : c=0 at interface return self.distanceVar._cellInterfaceFlag * 1e+20 def CUSTOMbuildMetalIonDiffusionEquation(ionVar = None, distanceVar = None, depositionRate = 1, transientCoeff = 1, diffusionCoeff = 1, metalIonMolarVolume = 1): eq = _buildLevelSetDiffusionEquation(ionVar = ionVar, distanceVar = distanceVar, transientCoeff = transientCoeff, diffusionCoeff = diffusionCoeff) coeff = _CUSTOMMetalIonSourceVariable(ionVar = ionVar, distanceVar = distanceVar, depositionRate = depositionRate, metalIonMolarVolume = metalIonMolarVolume) return eq + ImplicitSourceTerm(coeff) #mesh definition mesh = Grid2D(nx=nx, ny=ny, dx=dx,dy=dy) #Variables definitions c = CellVariable(name="c", mesh=mesh,value=0.,hasOld=1) #phi the distance variable narrowBandWidth = 10* cellSize phi = DistanceVariable(name='phi',mesh= mesh,value=1.,narrowBandWidth=narrowBandWidth,hasOld=1) x, y = mesh.getCellCenters() #phi.setValue(-1., where=( ((x-CX)**2 < Rstart**2 ) & ((y-CY)**2 < Rstart**2 ) ) ) phi.setValue(-1., where=((x-CX)**2+(y-CY)**2 < Rstart**2 )) phi.calcDistanceFunction() #depositionrate #no noise: #depositionRateVariable = (abs(v0*c.getGrad().dot(phi.getGrad())) + v0*c.getGrad().dot(phi.getGrad()) ) /2.0 + 1E-18 #with noise depositionRateVariable = abs(GaussianNoiseVariable(mesh=mesh, mean=1., variance=0.02))*( (abs(v0*c.getGrad().dot(phi.getGrad())) + v0*c.getGrad().dot(phi.getGrad()) ) /2.0 ) + 1E-18 #and velocityvariable extensionVelocityVariable = CellVariable( name = 'extension velocity', mesh = mesh, value = depositionRateVariable) #Advection equation for phi : dphi dt+vextgradphi=0 ou divergence ? signe ? #HERE advectionEquation = buildHigherOrderAdvectionEquation(advectionCoeff=extensionVelocityVariable) #metal equation (the particules concentration diffuses only where phi=1) ceq= CUSTOMbuildMetalIonDiffusionEquation(ionVar=c,distanceVar=phi,depositionRate=depositionRateVariable,diffusionCoeff=D) BCs = FixedValue(faces=mesh.getExteriorFaces(), value=1) #X, Y = mesh.getFaceCenters() #XX = FaceVariable(mesh=mesh, value=X) #YY = FaceVariable(mesh=mesh, value=Y) #facesOutside = ((XX-CX)**2+(YY-CY)**2 >= CX*CY -0.1 ) #c.constrain(1., where=facesOutside) #visualization: #try: # viewer = MayaviSurfactantViewer(phi, # zoomFactor=1e6, # datamax=1.0, # datamin=0.0, # smooth=1) #except: print "MAYAVI not tried!!" #viewer = MultiViewer(viewers=( # Viewer(c, datamin=0, datamax=1),Viewer(phi, datamin=-1e-9, datamax=1e-9),Viewer(extensionVelocityVariable, datamin=0, datamax=extensionVelocityVariable.max()))) vc=Viewer(c, datamin=0, datamax=1) vp=Viewer(phi, datamin=-1e-9, datamax=1e-9) #vpzo=Viewer(phi, datamin=-1, datamax=1) #ve=Viewer(extensionVelocityVariable, datamin=-extensionVelocityVariable.max(), datamax=extensionVelocityVariable.max()) ve=Viewer(extensionVelocityVariable, datamin=0, datamax=extensionVelocityVariable.max()) vg=Viewer(c.getGrad().dot(c.getGrad())) #RECORD INTERFACE IN DATA FILE def record_interface(): title='snapt'+str(step+10000)+'.dat' #f=open(title,'w') TSVViewer(vars=(phi)).plot(filename=title) # for field observation #f.close() ##main loop absolutetime=0 dt=0 sauvegarde=open("steps.dat","w") sauvegarde.close() for step in range(numberOfsteps): stepline= "step #"+str(step)+" time="+str(absolutetime)+" (dt="+str(dt)+")" print stepline sauvegarde=open("steps.dat","a") sauvegarde.write(stepline+"\n") sauvegarde.close() #viewer.plot() phi.calcDistanceFunction() extensionVelocityVariable.setValue(depositionRateVariable()) phi.updateOld() c.updateOld() phi.extendVariable(extensionVelocityVariable) dt = min(dt0 * cellSize / extensionVelocityVariable.max(),0.01) absolutetime=absolutetime+dt advectionEquation.solve(phi, dt=dt) ceq.solve(var=c, dt=dt,boundaryConditions=BCs) # ceq.solve(var=c, dt=dt) title='snapt'+str(step+10000)+'.png' #savefig(title,transparent=True) vc.plot() if step%everynimages==0: savefig('c'+title) vp.plot() if step%everynimages==0: savefig('p'+title) ve.plot() if step%everynimages==0: savefig('e'+title) vg.plot() if step%everynimages==0: savefig('g'+title) os.system("mv *.png images/")