/[escript]/trunk/escript/py_src/flows.py
ViewVC logotype

Diff of /trunk/escript/py_src/flows.py

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

revision 1467 by gross, Wed Apr 2 08:10:37 2008 UTC revision 1673 by gross, Thu Jul 24 22:28:50 2008 UTC
# Line 36  __date__="$Date:$" Line 36  __date__="$Date:$"
36  from escript import *  from escript import *
37  import util  import util
38  from linearPDEs import LinearPDE  from linearPDEs import LinearPDE
39  from pdetools import HomogeneousSaddlePointProblem  from pdetools import HomogeneousSaddlePointProblem,Projector
40    
41    class StokesProblemCartesian_DC(HomogeneousSaddlePointProblem):
42          """
43          solves
44    
45              -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i
46                    u_{i,i}=0
47    
48              u=0 where  fixed_u_mask>0
49              eta*(u_{i,j}+u_{j,i})*n_j=surface_stress
50    
51          if surface_stress is not give 0 is assumed.
52    
53          typical usage:
54    
55                sp=StokesProblemCartesian(domain)
56                sp.setTolerance()
57                sp.initialize(...)
58                v,p=sp.solve(v0,p0)
59          """
60          def __init__(self,domain,**kwargs):
61             HomogeneousSaddlePointProblem.__init__(self,**kwargs)
62             self.domain=domain
63             self.vol=util.integrate(1.,Function(self.domain))
64             self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
65             self.__pde_u.setSymmetryOn()
66             # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
67                
68             # self.__pde_proj=LinearPDE(domain,numEquations=1,numSolutions=1)
69             # self.__pde_proj.setReducedOrderOn()
70             # self.__pde_proj.setSymmetryOn()
71             # self.__pde_proj.setSolverMethod(LinearPDE.LUMPING)
72    
73          def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):
74            self.eta=eta
75            A =self.__pde_u.createCoefficientOfGeneralPDE("A")
76        self.__pde_u.setValue(A=Data())
77            for i in range(self.domain.getDim()):
78            for j in range(self.domain.getDim()):
79                A[i,j,j,i] += 1.
80                A[i,j,i,j] += 1.
81            # self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))
82            self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
83    
84            # self.__pde_proj.setValue(D=1/eta)
85            # self.__pde_proj.setValue(Y=1.)
86            # self.__inv_eta=util.interpolate(self.__pde_proj.getSolution(),ReducedFunction(self.domain))
87            self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))
88    
89          def B(self,arg):
90             a=util.div(arg, ReducedFunction(self.domain))
91             return a-util.integrate(a)/self.vol
92    
93          def inner(self,p0,p1):
94             return util.integrate(p0*p1)
95             # return util.integrate(1/self.__inv_eta**2*p0*p1)
96    
97          def getStress(self,u):
98             mg=util.grad(u)
99             return 2.*self.eta*util.symmetric(mg)
100          def getEtaEffective(self):
101             return self.eta
102    
103          def solve_A(self,u,p):
104             """
105             solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
106             """
107             self.__pde_u.setTolerance(self.getSubProblemTolerance())
108             self.__pde_u.setValue(X=-self.getStress(u),X_reduced=-p*util.kronecker(self.domain))
109             return  self.__pde_u.getSolution(verbose=self.show_details)
110    
111    
112          def solve_prec(self,p):
113            a=self.__inv_eta*p
114            return a-util.integrate(a)/self.vol
115    
116          def stoppingcriterium(self,Bv,v,p):
117              n_r=util.sqrt(self.inner(Bv,Bv))
118              n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))
119              if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v) , util.Lsup(v)
120              if self.iter == 0: self.__n_v=n_v;
121              self.__n_v, n_v_old =n_v, self.__n_v
122              self.iter+=1
123              if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():
124                  if self.verbose: print "PCG terminated after %s steps."%self.iter
125                  return True
126              else:
127                  return False
128          def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
129          if TOL==None:
130                 TOL=self.getTolerance()
131              if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)
132              self.iter+=1
133              
134              if norm_r <= norm_b*TOL:
135                  if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)
136                  return True
137              else:
138                  return False
139    
140    
141  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
142        """        """
# Line 63  class StokesProblemCartesian(Homogeneous Line 163  class StokesProblemCartesian(Homogeneous
163           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
164           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
165           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
166           self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
167                            
168           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
169           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 82  class StokesProblemCartesian(Homogeneous Line 182  class StokesProblemCartesian(Homogeneous
182          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
183              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
184              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
185      self.__pde_prec.setValue(D=1./self.eta)      self.__pde_prec.setValue(D=1/self.eta)
186          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
187    
188        def B(self,arg):        def B(self,arg):
# Line 96  class StokesProblemCartesian(Homogeneous Line 196  class StokesProblemCartesian(Homogeneous
196           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1,Function(self.domain))
197           return util.integrate(s0*s1)           return util.integrate(s0*s1)
198    
199          def inner_a(self,a0,a1):
200             p0=util.interpolate(a0[1],Function(self.domain))
201             p1=util.interpolate(a1[1],Function(self.domain))
202             alfa=(1/self.vol)*util.integrate(p0)
203             beta=(1/self.vol)*util.integrate(p1)
204         v0=util.grad(a0[0])
205         v1=util.grad(a1[0])
206             return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))
207    
208    
209        def getStress(self,u):        def getStress(self,u):
210           mg=util.grad(u)           mg=util.grad(u)
211           return 2.*self.eta*util.symmetric(mg)           return 2.*self.eta*util.symmetric(mg)
212          def getEtaEffective(self):
213             return self.eta
214    
215        def solve_A(self,u,p):        def solve_A(self,u,p):
216           """           """
217           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
218           """           """
219           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
220           self.__pde_u.setValue(X=-self.getStress(u)+p*util.kronecker(self.domain))           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))
221           return  self.__pde_u.getSolution(verbose=self.show_details)           return  self.__pde_u.getSolution(verbose=self.show_details)
222    
223    
224        def solve_prec(self,p):        def solve_prec(self,p):
225         #proj=Projector(domain=self.domain, reduce = True, fast=False)
226           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
227           self.__pde_prec.setValue(Y=p)           self.__pde_prec.setValue(Y=p)
228           q=self.__pde_prec.getSolution(verbose=self.show_details)           q=self.__pde_prec.getSolution(verbose=self.show_details)
229           return q           return q
230    
231          def solve_prec1(self,p):
232         #proj=Projector(domain=self.domain, reduce = True, fast=False)
233             self.__pde_prec.setTolerance(self.getSubProblemTolerance())
234             self.__pde_prec.setValue(Y=p)
235             q=self.__pde_prec.getSolution(verbose=self.show_details)
236         q0=util.interpolate(q,Function(self.domain))
237             print util.inf(q*q0),util.sup(q*q0)
238             q-=(1/self.vol)*util.integrate(q0)
239             print util.inf(q*q0),util.sup(q*q0)
240             return q
241    
242        def stoppingcriterium(self,Bv,v,p):        def stoppingcriterium(self,Bv,v,p):
243            n_r=util.sqrt(self.inner(Bv,Bv))            n_r=util.sqrt(self.inner(Bv,Bv))
244            n_v=util.Lsup(v)            n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))
245            if self.verbose: print "PCG step %s: L2(div(v)) = %s, Lsup(v)=%s"%(self.iter,n_r,n_v)            if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v)
246              if self.iter == 0: self.__n_v=n_v;
247              self.__n_v, n_v_old =n_v, self.__n_v
248            self.iter+=1            self.iter+=1
249            if n_r <= self.vol**(1./2.-1./self.domain.getDim())*n_v*self.getTolerance():            if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():
250                if self.verbose: print "PCG terminated after %s steps."%self.iter                if self.verbose: print "PCG terminated after %s steps."%self.iter
251                return True                return True
252            else:            else:
253                return False                return False
254        def stoppingcriterium_GMRES(self,norm_r,norm_b):        def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
255            if self.verbose: print "GMRES step %s: L2(r) = %s, L2(b)*TOL=%s"%(self.iter,norm_r,norm_b*self.getTolerance())        if TOL==None:
256                 TOL=self.getTolerance()
257              if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)
258            self.iter+=1            self.iter+=1
259            if norm_r <= norm_b*self.getTolerance():            
260                if self.verbose: print "GMRES terminated after %s steps."%self.iter            if norm_r <= norm_b*TOL:
261                  if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)
262                return True                return True
263            else:            else:
264                return False                return False
265    
266    

Legend:
Removed from v.1467  
changed lines
  Added in v.1673

  ViewVC Help
Powered by ViewVC 1.1.26