# 21 septembre : testing explicit source term (without c2 commented) #example from : http://www.ctcms.nist.gov/fipy/examples/levelSet/generated/examples.levelSet.electroChem.howToWriteAScript.html#module-examples.levelSet.electroChem.howToWriteAScript #Importing Libraries 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 # File 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 print "the radius is CX=",CX 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]) #Hand parameters circular=0 NOISE=0 # if 0 no noise, or else this is the value of the variance dtmin=30 #Nnarrowband=100 Nnarrowband=1 __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) #this is the typical absorbing boundary (email from Daniel dated july 12,2011) implicit term return self.distanceVar._cellInterfaceFlag * 1e+20 class _CUSTOMMetalIonSourceVariableEXPLICIT(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) value=0.25 #this is the absorbing boundary used before + a source term explicit term, email from Daniel dated july 12,2011 #debug section debug=0 if debug==1: print "DEBUGGAGE ON !!" print self.distanceVar._cellInterfaceFlag print len(self.distanceVar._cellInterfaceFlag) print type(self.distanceVar._cellInterfaceFlag) figure() X, Y = mesh.getCellCenters() for i in range(len(self.distanceVar._cellInterfaceFlag)): if self.distanceVar._cellInterfaceFlag[i]!=0: print "i=",i," ", print self.distanceVar._cellInterfaceFlag[i] print "X=",X[i] print "Y=",Y[i] plot(X[i],Y[i],'+') show() #raw_input("Press ENTER to exit") return self.distanceVar._cellInterfaceFlag * 1e+20 * value #* (x**2+y**2<(CX)**2) 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) coeffEXPLICIT=_CUSTOMMetalIonSourceVariableEXPLICIT(ionVar = ionVar, distanceVar = distanceVar, depositionRate = depositionRate, metalIonMolarVolume = metalIonMolarVolume) # set a given value at the interface return eq + ImplicitSourceTerm(coeff) - coeffEXPLICIT #or just absorbing boundary # return eq + ImplicitSourceTerm(coeff) #Mesh definition if circular==0: mesh = Grid2D(nx=nx, ny=ny, dx=dx,dy=dy) else: #circular mesh : radius=CX mesh = GmshImporter2D(''' cellSize = %(cellSize)g; radius = %(radius)g; Point(1) = {0, 0, 0, cellSize}; Point(2) = {-radius, 0, 0, cellSize}; Point(3) = {0, radius, 0, cellSize}; Point(4) = {radius, 0, 0, cellSize}; Point(5) = {0, -radius, 0, cellSize}; Circle(6) = {2, 1, 3}; Circle(7) = {3, 1, 4}; Circle(8) = {4, 1, 5}; Circle(9) = {5, 1, 2}; Line Loop(10) = {6, 7, 8, 9}; Plane Surface(11) = {10}; ''' % locals()) #Variables definitions c = CellVariable(name="c", mesh=mesh,value=0.,hasOld=1) #definition of intermediate variable (will be min of c and another concentration c2 later ...) cPRODUCT=c #phi the distance variable narrowBandWidth = Nnarrowband* cellSize phi = DistanceVariable(name='phi',mesh= mesh,value=1.,narrowBandWidth=narrowBandWidth,hasOld=1) x, y = mesh.getCellCenters() if circular==0: phi.setValue(-1., where=(((x-CX)/Rstart)**2+(0.1*(y-CY)/Rstart)**2 < 1 )) else: phi.setValue(-1., where=((x/Rstart)**2+(y/Rstart)**2 < 1 )) print "Rstart is ",Rstart phi.calcDistanceFunction() #depositionrate if NOISE==0: #no noise: depositionRateVariable = (abs(v0*cPRODUCT.getGrad().dot(phi.getGrad())) + v0*cPRODUCT.getGrad().dot(phi.getGrad()) ) /2.0 + 1E-18 else: #with noise depositionRateVariable = abs(GaussianNoiseVariable(mesh=mesh, mean=1., variance=NOISE))*( (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 div 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) # BC BCs = FixedValue(faces=mesh.getExteriorFaces(), value=1) X, Y = mesh.getFaceCenters() XX = FaceVariable(mesh=mesh, value=X) YY = FaceVariable(mesh=mesh, value=Y) if circular==0: sourcec2 = ((XX-CX)**2+(YY-CY)**2 <= Rstart*Rstart ) else: sourcec2 = ((XX)**2+(YY)**2 <= Rstart*Rstart ) vc=Viewer(c, datamin=0, datamax=1) vp=Viewer(phi, datamin=-1e-9, datamax=1e-9) ve=Viewer(extensionVelocityVariable, datamin=0, datamax=extensionVelocityVariable.max()) ##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() phi.calcDistanceFunction() extensionVelocityVariable.setValue(depositionRateVariable()) phi.updateOld() c.updateOld() phi.extendVariable(extensionVelocityVariable) dt = min(dt0 * cellSize / extensionVelocityVariable.max().getValue(), dtmin) absolutetime=absolutetime+dt advectionEquation.solve(phi, dt=dt) ceq.solve(var=c, dt=dt,boundaryConditions=BCs) title='snapt'+str(step+10000)+'.png' 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)