/[escript]/branches/scons-dev/escript/py_src/levelset.py
ViewVC logotype

Diff of /branches/scons-dev/escript/py_src/levelset.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1680 by ksteube, Thu Jul 24 04:46:46 2008 UTC revision 1681 by ksteube, Thu Jul 31 06:20:18 2008 UTC
# Line 27  class LevelSet(object): Line 27  class LevelSet(object):
27           self.velocity = None           self.velocity = None
28    
29           #==================================================           #==================================================
30           self.__fc=TransportPDE(self.__domain,num_equations=1,theta=0.5)           self.__FC=True
31           # self.__fc.setReducedOrderOn()           if self.__FC:
32           self.__fc.setValue(M=Scalar(1.,Function(self.__domain)))              self.__fc=TransportPDE(self.__domain,num_equations=1,theta=0.5)
33           self.__fc.setInitialSolution(phi+self.__diam)              # self.__fc.setReducedOrderOn()
34                self.__fc.setValue(M=Scalar(1.,Function(self.__domain)))
35                self.__fc.setInitialSolution(phi+self.__diam)
36             else:
37                self.__fcpde = LinearPDE(self.__domain,numEquations=1, numSolutions=1)
38                self.__fcpde.setSolverMethod(solver=LinearPDE.LUMPING)      
39                self.__fcpde.setReducedOrderOn()
40                self.__fcpde.setSymmetryOn()
41                self.__fcpde.setValue(D=1.)
42    
43           #=======================================           #=======================================
44           self.__reinitPde = LinearPDE(self.__domain,numEquations=1, numSolutions=1)           self.__reinitFC=False
45           self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING)                 if self.__reinitFC:
46           self.__reinitPde.setReducedOrderOn()              self.__reinitfc=TransportPDE(self.__domain,num_equations=1,theta=1.)
47           self.__reinitPde.setSymmetryOn()              self.__reinitfc.setValue(M=Scalar(1.,Function(self.__domain)))
48           self.__reinitPde.setValue(D=1.)              self.__reinitfc.setTolerance(1.e-5)
49           # self.__reinitPde.setTolerance(1.e-8)           else:
50           # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)              self.__reinitPde = LinearPDE(self.__domain,numEquations=1, numSolutions=1)
51                # self.__reinitPde.setSolverMethod(solver=LinearPDE.LUMPING)      
52                self.__reinitPde.setReducedOrderOn()
53                self.__reinitPde.setSymmetryOn()
54                self.__reinitPde.setValue(D=1.)
55                # self.__reinitPde.setTolerance(1.e-8)
56                # self.__reinitPde.setSolverMethod(preconditioner=LinearPDE.ILU0)
57           #=======================================           #=======================================
58           self.__updateInterface()           self.__updateInterface()
59           print "phi range:",inf(phi), sup(phi)           print "phi range:",inf(phi), sup(phi)
# Line 49  class LevelSet(object): Line 63  class LevelSet(object):
63           # self.__fc.setValue(q=q)           # self.__fc.setValue(q=q)
64           self.__reinitPde.setValue(q=q)           self.__reinitPde.setValue(q=q)
65    
66         def getH(self):
67             return self.__h
68       def getDomain(self):       def getDomain(self):
69           return self.__domain           return self.__domain
70    
# Line 57  class LevelSet(object): Line 73  class LevelSet(object):
73           returns a new dt for a given velocity using the courant coundition           returns a new dt for a given velocity using the courant coundition
74           """           """
75           self.velocity=velocity           self.velocity=velocity
76           self.__fc.setValue(C=-interpolate(velocity,Function(self.__domain)))           if self.__FC:
77           dt=self.__fc.getSafeTimeStepSize()              self.__fc.setValue(C=-interpolate(velocity,Function(self.__domain)))
78                dt=self.__fc.getSafeTimeStepSize()
79             else:
80                dt=0.5*self.__h/sup(length(velocity))
81           return dt           return dt
82    
83       def getLevelSetFunction(self):       def getLevelSetFunction(self):
# Line 101  class LevelSet(object): Line 120  class LevelSet(object):
120    
121           @param dt: time step forward           @param dt: time step forward
122           """           """
123           self.__phi=self.__fc.solve(dt, verbose=False)-self.__diam           if self.__FC:
124                self.__phi=self.__fc.solve(dt, verbose=False)-self.__diam
125             else:
126                self.__fcpde.setValue(Y = self.__phi-(dt/2.)*inner(self.velocity,grad(self.__phi)))
127                phi_half = self.__fcpde.getSolution()
128                self.__fcpde.setValue(Y = self.__phi-dt*inner(self.velocity,grad(phi_half)))
129                self.__phi= self.__fcpde.getSolution()
130           self.__update_count += 1           self.__update_count += 1
131           if self.__update_count%self.__reinitialization_each == 0:           if self.__update_count%self.__reinitialization_each == 0:
132              self.__phi=self.__reinitialise(self.__phi)              self.__phi=self.__reinitialise(self.__phi)
133              self.__fc.setInitialSolution(self.__phi+self.__diam)              if self.__FC:
134                    self.__fc.setInitialSolution(self.__phi+self.__diam)
135           self.__updateInterface()           self.__updateInterface()
136    
137       def setTolerance(self,tolerance=1e-3):       def setTolerance(self,tolerance=1e-3):
138           self.__fc.setTolerance(tolerance)           if self.__FC:
139                self.__fc.setTolerance(tolerance)
140    
141       def __reinitialise(self,phi):       def __reinitialise(self,phi):
142           print "reintialization started:"           print "reintialization started:"
143           s=self.__makeInterface(phi,1.)           s=self.__makeInterface(phi,1.5)
144           g=grad(phi)           g=grad(phi)
145           w = s*g/(length(g)+EPSILON)           w = s*g/(length(g)+EPSILON)
146           dtau = 0.5*inf(Function(self.__domain).getSize())           if self.__reinitFC:
147           print "step size: dt = ",dtau               self.__reinitfc.setValue(C=-w,Y=s)
148                 self.__reinitfc.setInitialSolution(phi+self.__diam)
149                 # dtau=self.__reinitfc.getSafeTimeStepSize()
150                 dtau = 0.5*inf(Function(self.__domain).getSize())
151             else:
152                 dtau = 0.5*inf(Function(self.__domain).getSize())
153             print "step size: dt = ",dtau,inf(abs(phi.interpolate(Function(self.__domain)))/abs(s-inner(w,grad(phi))))
154           iter =0           iter =0
155           # self.__reinitPde.setValue(q=whereNegative(abs(phi)-self.__h), r=phi)           # self.__reinitPde.setValue(q=whereNegative(abs(phi)-2*self.__h), r=phi)
156           # self.__reinitPde.setValue(r=phi)           # self.__reinitPde.setValue(r=phi)
157           while (iter<=self.__reinitialization_steps_max):           while (iter<=self.__reinitialization_steps_max):
158                   phi_old=phi                   phi_old=phi
159                   self.__reinitPde.setValue(Y = phi+(dtau/2.)*(s-inner(w,grad(phi))))                   if self.__reinitFC:
160                   phi_half = self.__reinitPde.getSolution()                     phi = self.__reinitfc.solve(dtau)-self.__diam
161                   self.__reinitPde.setValue(Y = phi+dtau*(s-inner(w,grad(phi_half))))                   else:
162                   phi = self.__reinitPde.getSolution()                     self.__reinitPde.setValue(Y = phi+(dtau/2.)*(s-inner(w,grad(phi))))
163                       phi_half = self.__reinitPde.getSolution()
164                       self.__reinitPde.setValue(Y = phi+dtau*(s-inner(w,grad(phi_half))))
165                       # g=grad(phi)
166                       # S=inner(w,grad(phi))
167                       # self.__reinitPde.setValue(Y = phi+dtau*(s-S),X=dtau**2/2*w*(s-S))
168                       phi = self.__reinitPde.getSolution()
169                   change = Lsup(phi-phi_old)/self.__diam                   change = Lsup(phi-phi_old)/self.__diam
170                   print "phi range:",inf(phi), sup(phi)                   print "phi range:",inf(phi), sup(phi)
171                   print "iteration :", iter, " error:", change                   print "iteration :", iter, " error:", change
Legend:
Removed from v.1680  
changed lines
  Added in v.1681
  ViewVC Help
Powered by ViewVC 1.1.26