/[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 1841 by gross, Fri Oct 3 03:57:52 2008 UTC revision 2156 by gross, Mon Dec 15 05:09:02 2008 UTC
# Line 1  Line 1 
   
1  ########################################################  ########################################################
2  #  #
3  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2008 by University of Queensland
# Line 34  __author__="Lutz Gross, l.gross@uq.edu.a Line 33  __author__="Lutz Gross, l.gross@uq.edu.a
33    
34  from escript import *  from escript import *
35  import util  import util
36  from linearPDEs import LinearPDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE
37  from pdetools import HomogeneousSaddlePointProblem,Projector  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm
   
 class StokesProblemCartesian_DC(HomogeneousSaddlePointProblem):  
       """  
       solves  
   
           -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i  
                 u_{i,i}=0  
   
           u=0 where  fixed_u_mask>0  
           eta*(u_{i,j}+u_{j,i})*n_j=surface_stress  
   
       if surface_stress is not give 0 is assumed.  
   
       typical usage:  
   
             sp=StokesProblemCartesian(domain)  
             sp.setTolerance()  
             sp.initialize(...)  
             v,p=sp.solve(v0,p0)  
       """  
       def __init__(self,domain,**kwargs):  
          HomogeneousSaddlePointProblem.__init__(self,**kwargs)  
          self.domain=domain  
          self.vol=util.integrate(1.,Function(self.domain))  
          self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())  
          self.__pde_u.setSymmetryOn()  
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)  
               
          # self.__pde_proj=LinearPDE(domain,numEquations=1,numSolutions=1)  
          # self.__pde_proj.setReducedOrderOn()  
          # self.__pde_proj.setSymmetryOn()  
          # self.__pde_proj.setSolverMethod(LinearPDE.LUMPING)  
   
       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):  
         self.eta=eta  
         A =self.__pde_u.createCoefficient("A")  
     self.__pde_u.setValue(A=Data())  
         for i in range(self.domain.getDim()):  
         for j in range(self.domain.getDim()):  
             A[i,j,j,i] += 1.  
             A[i,j,i,j] += 1.  
         # self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)  
38    
39          # self.__pde_proj.setValue(D=1/eta)  class DarcyFlow(object):
40          # self.__pde_proj.setValue(Y=1.)      """
41          # self.__inv_eta=util.interpolate(self.__pde_proj.getSolution(),ReducedFunction(self.domain))      solves the problem
42          self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))  
43        M{u_i+k_{ij}*p_{,j} = g_i}
44        def B(self,arg):      M{u_{i,i} = f}
45           a=util.div(arg, ReducedFunction(self.domain))  
46           return a-util.integrate(a)/self.vol      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,
47    
48        def inner(self,p0,p1):      @note: The problem is solved in a least squares formulation.
49           return util.integrate(p0*p1)      """
50           # return util.integrate(1/self.__inv_eta**2*p0*p1)  
51        def __init__(self, domain):
52        def getStress(self,u):          """
53           mg=util.grad(u)          initializes the Darcy flux problem
54           return 2.*self.eta*util.symmetric(mg)          @param domain: domain of the problem
55        def getEtaEffective(self):          @type domain: L{Domain}
56           return self.eta          """
57            self.domain=domain
58        def solve_A(self,u,p):          self.__pde_v=LinearPDESystem(domain)
59           """          self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))
60           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)          self.__pde_v.setSymmetryOn()
61            self.__pde_p=LinearSinglePDE(domain)
62            self.__pde_p.setSymmetryOn()
63            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
64            self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
65    
66        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
67            """
68            assigns values to model parameters
69    
70            @param f: volumetic sources/sinks
71            @type f: scalar value on the domain (e.g. L{Data})
72            @param g: flux sources/sinks
73            @type g: vector values on the domain (e.g. L{Data})
74            @param location_of_fixed_pressure: mask for locations where pressure is fixed
75            @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
76            @param location_of_fixed_flux:  mask for locations where flux is fixed.
77            @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
78            @param permeability: permeability tensor. If scalar C{s} is given the tensor with
79                                 C{s} on the main diagonal is used. If vector C{v} is given the tensor with
80                                 C{v} on the main diagonal is used.
81            @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
82    
83            @note: the values of parameters which are not set by calling C{setValue} are not altered.
84            @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
85                   or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
86                   is along the M{x_i} axis.
87            """
88            if f !=None:
89               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
90               if f.isEmpty():
91                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
92               else:
93                   if f.getRank()>0: raise ValueError,"illegal rank of f."
94               self.f=f
95            if g !=None:  
96               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
97               if g.isEmpty():
98                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
99               else:
100                 if not g.getShape()==(self.domain.getDim(),):
101                   raise ValueError,"illegal shape of g"
102               self.__g=g
103    
104            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
105            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
106    
107            if permeability!=None:
108               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
109               if perm.getRank()==0:
110                   perm=perm*util.kronecker(self.domain.getDim())
111               elif perm.getRank()==1:
112                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
113                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
114               elif perm.getRank()==2:
115                  pass
116               else:
117                  raise ValueError,"illegal rank of permeability."
118               self.__permeability=perm
119               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
120    
121    
122        def getFlux(self,p, fixed_flux=Data(),tol=1.e-8, show_details=False):
123            """
124            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
125            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
126            Note that C{g} and C{f} are used, L{setValue}.
127            
128            @param p: pressure.
129            @type p: scalar value on the domain (e.g. L{Data}).
130            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
131            @type fixed_flux: vector values on the domain (e.g. L{Data}).
132            @param tol: relative tolerance to be used.
133            @type tol: positive float.
134            @return: flux
135            @rtype: L{Data}
136            @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}
137                   for the permeability M{k_{ij}}
138            """
139            self.__pde_v.setTolerance(tol)
140            self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux)
141            return self.__pde_v.getSolution(verbose=show_details)
142    
143        def solve(self,u0,p0,atol=0,rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
144             """
145             solves the problem.
146    
147             The iteration is terminated if the error in the pressure is less then C{rtol * |q| + atol} where
148             C{|q|} denotes the norm of the right hand side (see escript user's guide for details).
149    
150             @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged.
151             @type u0: vector value on the domain (e.g. L{Data}).
152             @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.
153             @type p0: scalar value on the domain (e.g. L{Data}).
154             @param atol: absolute tolerance for the pressure
155             @type atol: non-negative C{float}
156             @param rtol: relative tolerance for the pressure
157             @type rtol: non-negative C{float}
158             @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}
159             @type sub_rtol: positive-negative C{float}
160             @param verbose: if set some information on iteration progress are printed
161             @type verbose: C{bool}
162             @param show_details:  if set information on the subiteration process are printed.
163             @type show_details: C{bool}
164             @return: flux and pressure
165             @rtype: C{tuple} of L{Data}.
166    
167             @note: The problem is solved as a least squares form
168    
169             M{(I+D^*D)u+Qp=D^*f+g}
170             M{Q^*u+Q^*Qp=Q^*g}
171    
172             where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
173             We eliminate the flux form the problem by setting
174    
175             M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
176    
177             form the first equation. Inserted into the second equation we get
178    
179             M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with p=p0  on location_of_fixed_pressure
180            
181             which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
182             PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
183           """           """
184           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.verbose=verbose
185           self.__pde_u.setValue(X=-self.getStress(u),X_reduced=-p*util.kronecker(self.domain))           self.show_details= show_details and self.verbose
186           return  self.__pde_u.getSolution(verbose=self.show_details)           self.__pde_v.setTolerance(sub_rtol)
187             self.__pde_p.setTolerance(sub_rtol)
188             u2=u0*self.__pde_v.getCoefficient("q")
189        def solve_prec(self,p):           #
190          a=self.__inv_eta*p           # first the reference velocity is calculated from
191          return a-util.integrate(a)/self.vol           #
192             #   (I+D^*D)u_ref=D^*f+g (including bundray conditions for u)
193        def stoppingcriterium(self,Bv,v,p):           #
194            n_r=util.sqrt(self.inner(Bv,Bv))           self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=u0)
195            n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))           u_ref=self.__pde_v.getSolution(verbose=show_details)
196            if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v) , util.Lsup(v)           if self.verbose: print "DarcyFlux: maximum reference flux = ",util.Lsup(u_ref)
197            if self.iter == 0: self.__n_v=n_v;           self.__pde_v.setValue(r=Data())
198            self.__n_v, n_v_old =n_v, self.__n_v           #
199            self.iter+=1           #   and then we calculate a reference pressure
200            if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():           #
201                if self.verbose: print "PCG terminated after %s steps."%self.iter           #       Q^*Qp_ref=Q^*g-Q^*u_ref ((including bundray conditions for p)
202                return True           #
203            else:           self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,(self.__g-u_ref)), r=p0)
204                return False           p_ref=self.__pde_p.getSolution(verbose=self.show_details)
205        def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):           if self.verbose: print "DarcyFlux: maximum reference pressure = ",util.Lsup(p_ref)
206        if TOL==None:           self.__pde_p.setValue(r=Data())
207               TOL=self.getTolerance()           #
208            if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)           #   (I+D^*D)du + Qdp = - Qp_ref                       u=du+u_ref
209            self.iter+=1           #   Q^*du + Q^*Qdp = Q^*g-Q^*u_ref-Q^*Qp_ref=0        p=dp+pref
210                       #
211            if norm_r <= norm_b*TOL:           #      du= -(I+D^*D)^(-1} Q(p_ref+dp)  u = u_ref+du
212                if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)           #
213                return True           #  => Q^*(I-(I+D^*D)^(-1})Q dp = Q^*(I+D^*D)^(-1} Qp_ref
214            else:           #  or Q^*(I-(I+D^*D)^(-1})Q p = Q^*Qp_ref
215                return False           #
216             #   r= Q^*( (I+D^*D)^(-1} Qp_ref - Q dp + (I+D^*D)^(-1})Q dp) = Q^*(-du-Q dp)
217             #            with du=-(I+D^*D)^(-1} Q(p_ref+dp)
218             #
219             #  we use the (du,Qdp) to represent the resudual
220             #  Q^*Q is a preconditioner
221             #  
222             #  <(Q^*Q)^{-1}r,r> -> right hand side norm is <Qp_ref,Qp_ref>
223             #
224             Qp_ref=util.tensor_mult(self.__permeability,util.grad(p_ref))
225             norm_rhs=util.sqrt(util.integrate(util.inner(Qp_ref,Qp_ref)))
226             ATOL=max(norm_rhs*rtol +atol, 200. * util.EPSILON * norm_rhs)
227             if not ATOL>0:
228                 raise ValueError,"Negative absolute tolerance (rtol = %e, norm right hand side =%, atol =%e)."%(rtol, norm_rhs, atol)
229             if self.verbose: print "DarcyFlux: norm of right hand side = %e (absolute tolerance = %e)"%(norm_rhs,ATOL)
230             #
231             #   caclulate the initial residual
232             #
233             self.__pde_v.setValue(X=Data(), Y=-util.tensor_mult(self.__permeability,util.grad(p0)), r=Data())
234             du=self.__pde_v.getSolution(verbose=show_details)
235             r=ArithmeticTuple(util.tensor_mult(self.__permeability,util.grad(p0-p_ref)), du)
236             dp,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
237             util.saveVTK("d.vtu",p=dp,p_ref=p_ref)
238             return u_ref+r[1],dp
239            
240        def __Aprod_PCG(self,p):
241              if self.show_details: print "DarcyFlux: Applying operator"
242              Qp=util.tensor_mult(self.__permeability,util.grad(p))
243              self.__pde_v.setValue(Y=Qp,X=Data())
244              w=self.__pde_v.getSolution(verbose=self.show_details)
245              return ArithmeticTuple(-Qp,w)
246    
247        def __inner_PCG(self,p,r):
248             a=util.tensor_mult(self.__permeability,util.grad(p))
249             out=-util.integrate(util.inner(a,r[0]+r[1]))
250             return out
251    
252        def __Msolve_PCG(self,r):
253              if self.show_details: print "DarcyFlux: Applying preconditioner"
254              self.__pde_p.setValue(X=-util.transposed_tensor_mult(self.__permeability,r[0]+r[1]))
255              return self.__pde_p.getSolution(verbose=self.show_details)
256    
257  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
258        """        """
259        solves        solves
260    
261            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
262                  u_{i,i}=0                  u_{i,i}=0
263    
264            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
265            eta*(u_{i,j}+u_{j,i})*n_j=surface_stress            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
266    
267        if surface_stress is not give 0 is assumed.        if surface_stress is not given 0 is assumed.
268    
269        typical usage:        typical usage:
270    
# Line 157  class StokesProblemCartesian(Homogeneous Line 274  class StokesProblemCartesian(Homogeneous
274              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
275        """        """
276        def __init__(self,domain,**kwargs):        def __init__(self,domain,**kwargs):
277             """
278             initialize the Stokes Problem
279    
280             @param domain: domain of the problem. The approximation order needs to be two.
281             @type domain: L{Domain}
282             @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
283             """
284           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
285           self.domain=domain           self.domain=domain
286           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
287           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())
288           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
289           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
290             # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
291                            
292           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
293           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
294             # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
295           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
296    
297           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
# Line 173  class StokesProblemCartesian(Homogeneous Line 299  class StokesProblemCartesian(Homogeneous
299           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
300           self.__pde_proj.setValue(D=1.)           self.__pde_proj.setValue(D=1.)
301    
302        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):
303            """
304            assigns values to the model parameters
305    
306            @param f: external force
307            @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
308            @param fixed_u_mask: mask of locations with fixed velocity.
309            @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
310            @param eta: viscosity
311            @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
312            @param surface_stress: normal surface stress
313            @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
314            @param stress: initial stress
315        @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
316            @note: All values needs to be set.
317    
318            """
319          self.eta=eta          self.eta=eta
320          A =self.__pde_u.createCoefficient("A")          A =self.__pde_u.createCoefficient("A")
321      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
# Line 183  class StokesProblemCartesian(Homogeneous Line 325  class StokesProblemCartesian(Homogeneous
325              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
326      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
327          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)
328            self.__stress=stress
329    
330          def B(self,v):
331            """
332            returns div(v)
333            @rtype: equal to the type of p
334    
335            @note: boundary conditions on p should be zero!
336            """
337            if self.show_details: print "apply divergence:"
338            self.__pde_proj.setValue(Y=-util.div(v))
339            self.__pde_proj.setTolerance(self.getSubProblemTolerance())
340            return self.__pde_proj.getSolution(verbose=self.show_details)
341    
342          def inner_pBv(self,p,Bv):
343             """
344             returns inner product of element p and Bv  (overwrite)
345            
346             @type p: equal to the type of p
347             @type Bv: equal to the type of result of operator B
348             @rtype: C{float}
349    
350        def B(self,arg):           @rtype: equal to the type of p
351           d=util.div(arg)           """
352           self.__pde_proj.setValue(Y=d)           s0=util.interpolate(p,Function(self.domain))
353           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           s1=util.interpolate(Bv,Function(self.domain))
          return self.__pde_proj.getSolution(verbose=self.show_details)  
   
       def inner(self,p0,p1):  
          s0=util.interpolate(p0,Function(self.domain))  
          s1=util.interpolate(p1,Function(self.domain))  
354           return util.integrate(s0*s1)           return util.integrate(s0*s1)
355    
356        def inner_a(self,a0,a1):        def inner_p(self,p0,p1):
357           p0=util.interpolate(a0[1],Function(self.domain))           """
358           p1=util.interpolate(a1[1],Function(self.domain))           returns inner product of element p0 and p1  (overwrite)
359           alfa=(1/self.vol)*util.integrate(p0)          
360           beta=(1/self.vol)*util.integrate(p1)           @type p0: equal to the type of p
361       v0=util.grad(a0[0])           @type p1: equal to the type of p
362       v1=util.grad(a1[0])           @rtype: C{float}
363           return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))  
364             @rtype: equal to the type of p
365             """
366             s0=util.interpolate(p0/self.eta,Function(self.domain))
367             s1=util.interpolate(p1/self.eta,Function(self.domain))
368             return util.integrate(s0*s1)
369    
370          def inner_v(self,v0,v1):
371             """
372             returns inner product of two element v0 and v1  (overwrite)
373            
374             @type v0: equal to the type of v
375             @type v1: equal to the type of v
376             @rtype: C{float}
377    
378        def getStress(self,u):           @rtype: equal to the type of v
379           mg=util.grad(u)           """
380           return 2.*self.eta*util.symmetric(mg)       gv0=util.grad(v0)
381        def getEtaEffective(self):       gv1=util.grad(v1)
382           return self.eta           return util.integrate(util.inner(gv0,gv1))
383    
384        def solve_A(self,u,p):        def solve_A(self,u,p):
385           """           """
386           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
387           """           """
388             if self.show_details: print "solve for velocity:"
389           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
390           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))           if self.__stress.isEmpty():
391           return  self.__pde_u.getSolution(verbose=self.show_details)              self.__pde_u.setValue(X=-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
392             else:
393                self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
394             out=self.__pde_u.getSolution(verbose=self.show_details)
395             return  out
396    
397        def solve_prec(self,p):        def solve_prec(self,p):
398       #proj=Projector(domain=self.domain, reduce = True, fast=False)           if self.show_details: print "apply preconditioner:"
399           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
400           self.__pde_prec.setValue(Y=p)           self.__pde_prec.setValue(Y=p)
401           q=self.__pde_prec.getSolution(verbose=self.show_details)           q=self.__pde_prec.getSolution(verbose=self.show_details)
402           return q           return q
   
       def solve_prec1(self,p):  
      #proj=Projector(domain=self.domain, reduce = True, fast=False)  
          self.__pde_prec.setTolerance(self.getSubProblemTolerance())  
          self.__pde_prec.setValue(Y=p)  
          q=self.__pde_prec.getSolution(verbose=self.show_details)  
      q0=util.interpolate(q,Function(self.domain))  
          print util.inf(q*q0),util.sup(q*q0)  
          q-=(1/self.vol)*util.integrate(q0)  
          print util.inf(q*q0),util.sup(q*q0)  
          return q  
   
       def stoppingcriterium(self,Bv,v,p):  
           n_r=util.sqrt(self.inner(Bv,Bv))  
           n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))  
           if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v)  
           if self.iter == 0: self.__n_v=n_v;  
           self.__n_v, n_v_old =n_v, self.__n_v  
           self.iter+=1  
           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.verbose: print "PCG terminated after %s steps."%self.iter  
               return True  
           else:  
               return False  
       def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):  
       if TOL==None:  
              TOL=self.getTolerance()  
           if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)  
           self.iter+=1  
             
           if norm_r <= norm_b*TOL:  
               if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)  
               return True  
           else:  
               return False  
   
   

Legend:
Removed from v.1841  
changed lines
  Added in v.2156

  ViewVC Help
Powered by ViewVC 1.1.26