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

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

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

revision 1639 by gross, Mon Jul 14 08:55:25 2008 UTC revision 1659 by gross, Fri Jul 18 02:28:13 2008 UTC
# Line 38  import util Line 38  import util
38  from linearPDEs import LinearPDE  from linearPDEs import LinearPDE
39  from pdetools import HomogeneousSaddlePointProblem,Projector  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    
88          def B(self,arg):
89             a=util.div(arg, ReducedFunction(self.domain))
90             return a-util.integrate(a)/self.vol
91    
92          def inner(self,p0,p1):
93             return util.integrate(p0*p1)
94    
95          def getStress(self,u):
96             mg=util.grad(u)
97             return 2.*self.eta*util.symmetric(mg)
98          def getEtaEffective(self):
99             return self.eta
100    
101          def solve_A(self,u,p):
102             """
103             solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
104             """
105             self.__pde_u.setTolerance(self.getSubProblemTolerance())
106             self.__pde_u.setValue(X=-self.getStress(u),X_reduced=-p*util.kronecker(self.domain))
107             return  self.__pde_u.getSolution(verbose=self.show_details)
108    
109    
110          def solve_prec(self,p):
111            a=self.__inv_eta*p
112            return a-util.integrate(a)/self.vol
113    
114          def stoppingcriterium(self,Bv,v,p):
115              n_r=util.sqrt(self.inner(Bv,Bv))
116              n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))
117              if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v) , util.Lsup(v)
118              if self.iter == 0: self.__n_v=n_v;
119              self.__n_v, n_v_old =n_v, self.__n_v
120              self.iter+=1
121              if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():
122                  if self.verbose: print "PCG terminated after %s steps."%self.iter
123                  return True
124              else:
125                  return False
126    
127    
128  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
129        """        """
130        solves        solves
# Line 103  class StokesProblemCartesian(Homogeneous Line 190  class StokesProblemCartesian(Homogeneous
190           beta=(1/self.vol)*util.integrate(p1)           beta=(1/self.vol)*util.integrate(p1)
191       v0=util.grad(a0[0])       v0=util.grad(a0[0])
192       v1=util.grad(a1[0])       v1=util.grad(a1[0])
      #print "NORM",alfa,beta,util.integrate((p0-alfa)*(p1-beta))+util.integrate(util.inner(v0,v1))  
193           return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))           return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))
194    
195    
# Line 147  class StokesProblemCartesian(Homogeneous Line 233  class StokesProblemCartesian(Homogeneous
233            if self.iter == 0: self.__n_v=n_v;            if self.iter == 0: self.__n_v=n_v;
234            self.__n_v, n_v_old =n_v, self.__n_v            self.__n_v, n_v_old =n_v, self.__n_v
235            self.iter+=1            self.iter+=1
           print abs(n_v_old-self.__n_v), n_v, self.getTolerance()  
236            if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= 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():
237                if self.verbose: print "PCG terminated after %s steps."%self.iter                if self.verbose: print "PCG terminated after %s steps."%self.iter
238                return True                return True

Legend:
Removed from v.1639  
changed lines
  Added in v.1659

  ViewVC Help
Powered by ViewVC 1.1.26