Hi all, thanks to your continuous help, I have now a running code for simulating typical diffusion limited aggregation with fipy. to remove anisotropy, I replaced the 2D carthesian grid, with an irregular mesh. It works great(anisotropy is removed) except I got know an additional warning (not really an error as the code is running with consistent behaviour) but it seems to slow down a lot the code :
here is anexample : step #106 time=155.476602908 (dt=0.351844720871) circularlevelset.py:226: MaximumIterationWarning: Iterations: 1001. Relative error: 6.38479e-08 ceq.solve(var=c, dt=dt,boundaryConditions=BCs) step #107 time=155.873131134 (dt=0.396528225831) circularlevelset.py:226: MaximumIterationWarning: Iterations: 1001. Relative error: 4.1797e-09 ceq.solve(var=c, dt=dt,boundaryConditions=BCs) step #108 time=156.473576822 (dt=0.600445687855) circularlevelset.py:226: MaximumIterationWarning: Iterations: 1001. Relative error: 2.54615e-08 ceq.solve(var=c, dt=dt,boundaryConditions=BCs) I attached the code, if someone get the courage and the time to look at it ..or do you have ideas just like that ? thanks again so much, Julien
# 16 mai 2011 : tentative de mettre une mesh circular #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]) NOISE=0 # if 0 no noise, or else this is the value of the variance dtmin=100 Nnarrowband=100 __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) #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) #phi the distance variable narrowBandWidth = Nnarrowband* cellSize phi = DistanceVariable(name='phi',mesh= mesh,value=1.,narrowBandWidth=narrowBandWidth,hasOld=1) x, y = mesh.getCellCenters() #phi.setValue(-1., where=((x-CX)**2+(y-CY)**2 < Rstart**2 )) phi.setValue(-1., where=((x)**2+(y)**2 < Rstart**2 )) phi.calcDistanceFunction() #depositionrate if NOISE==0: #no noise: depositionRateVariable = (abs(v0*c.getGrad().dot(phi.getGrad())) + v0*c.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 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) #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) ###OLDERROR ??dt = min(dt0 * cellSize / extensionVelocityVariable.max(),0.01) 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' #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/")
parameters.input
Description: parameters.input
