/[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 1673 by gross, Thu Jul 24 22:28:50 2008 UTC revision 2208 by gross, Mon Jan 12 06:37:07 2009 UTC
# Line 1  Line 1 
1  # $Id:$  ########################################################
2  #  #
3  #######################################################  # Copyright (c) 2003-2008 by University of Queensland
4    # Earth Systems Science Computational Center (ESSCC)
5    # http://www.uq.edu.au/esscc
6  #  #
7  #       Copyright 2008 by University of Queensland  # Primary Business: Queensland, Australia
8  #  # Licensed under the Open Software License version 3.0
9  #                http://esscc.uq.edu.au  # http://www.opensource.org/licenses/osl-3.0.php
 #        Primary Business: Queensland, Australia  
 #  Licensed under the Open Software License version 3.0  
 #     http://www.opensource.org/licenses/osl-3.0.php  
 #  
 #######################################################  
10  #  #
11    ########################################################
12    
13    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
14    Earth Systems Science Computational Center (ESSCC)
15    http://www.uq.edu.au/esscc
16    Primary Business: Queensland, Australia"""
17    __license__="""Licensed under the Open Software License version 3.0
18    http://www.opensource.org/licenses/osl-3.0.php"""
19    __url__="http://www.uq.edu.au/esscc/escript-finley"
20    
21  """  """
22  Some models for flow  Some models for flow
# Line 24  Some models for flow Line 30  Some models for flow
30  """  """
31    
32  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2008 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys"  
 __version__="$Revision:$"  
 __date__="$Date:$"  
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.createCoefficientOfGeneralPDE("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)  
   
         # self.__pde_proj.setValue(D=1/eta)  
         # self.__pde_proj.setValue(Y=1.)  
         # self.__inv_eta=util.interpolate(self.__pde_proj.getSolution(),ReducedFunction(self.domain))  
         self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))  
   
       def B(self,arg):  
          a=util.div(arg, ReducedFunction(self.domain))  
          return a-util.integrate(a)/self.vol  
   
       def inner(self,p0,p1):  
          return util.integrate(p0*p1)  
          # return util.integrate(1/self.__inv_eta**2*p0*p1)  
   
       def getStress(self,u):  
          mg=util.grad(u)  
          return 2.*self.eta*util.symmetric(mg)  
       def getEtaEffective(self):  
          return self.eta  
38    
39        def solve_A(self,u,p):  class DarcyFlow(object):
40        """
41        solves the problem
42    
43        M{u_i+k_{ij}*p_{,j} = g_i}
44        M{u_{i,i} = f}
45    
46        where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,
47    
48        @note: The problem is solved in a least squares formulation.
49        """
50    
51        def __init__(self, domain,useReduced=False):
52            """
53            initializes the Darcy flux problem
54            @param domain: domain of the problem
55            @type domain: L{Domain}
56            """
57            self.domain=domain
58            self.__pde_v=LinearPDESystem(domain)
59            if useReduced: self.__pde_v.setReducedOrderOn()
60            self.__pde_v.setSymmetryOn()
61            self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))
62            self.__pde_p=LinearSinglePDE(domain)
63            self.__pde_p.setSymmetryOn()
64            if useReduced: self.__pde_p.setReducedOrderOn()
65            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
66            self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
67            self.__ATOL= None
68    
69        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
70            """
71            assigns values to model parameters
72    
73            @param f: volumetic sources/sinks
74            @type f: scalar value on the domain (e.g. L{Data})
75            @param g: flux sources/sinks
76            @type g: vector values on the domain (e.g. L{Data})
77            @param location_of_fixed_pressure: mask for locations where pressure is fixed
78            @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
79            @param location_of_fixed_flux:  mask for locations where flux is fixed.
80            @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
81            @param permeability: permeability tensor. If scalar C{s} is given the tensor with
82                                 C{s} on the main diagonal is used. If vector C{v} is given the tensor with
83                                 C{v} on the main diagonal is used.
84            @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
85    
86            @note: the values of parameters which are not set by calling C{setValue} are not altered.
87            @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
88                   or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
89                   is along the M{x_i} axis.
90            """
91            if f !=None:
92               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
93               if f.isEmpty():
94                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
95               else:
96                   if f.getRank()>0: raise ValueError,"illegal rank of f."
97               self.f=f
98            if g !=None:  
99               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
100               if g.isEmpty():
101                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
102               else:
103                 if not g.getShape()==(self.domain.getDim(),):
104                   raise ValueError,"illegal shape of g"
105               self.__g=g
106    
107            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
108            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
109    
110            if permeability!=None:
111               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
112               if perm.getRank()==0:
113                   perm=perm*util.kronecker(self.domain.getDim())
114               elif perm.getRank()==1:
115                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
116                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
117               elif perm.getRank()==2:
118                  pass
119               else:
120                  raise ValueError,"illegal rank of permeability."
121               self.__permeability=perm
122               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
123    
124    
125        def getFlux(self,p=None, fixed_flux=Data(),tol=1.e-8, show_details=False):
126            """
127            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
128            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
129            Note that C{g} and C{f} are used, see L{setValue}.
130            
131            @param p: pressure.
132            @type p: scalar value on the domain (e.g. L{Data}).
133            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
134            @type fixed_flux: vector values on the domain (e.g. L{Data}).
135            @param tol: relative tolerance to be used.
136            @type tol: positive C{float}.
137            @return: flux
138            @rtype: L{Data}
139            @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}}
140                   for the permeability M{k_{ij}}
141            """
142            self.__pde_v.setTolerance(tol)
143            g=self.__g
144            f=self.__f
145            self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux)
146            if p == None:
147               self.__pde_v.setValue(Y=g)
148            else:
149               self.__pde_v.setValue(Y=g-util.tensor_mult(self.__permeability,util.grad(p)))
150            return self.__pde_v.getSolution(verbose=show_details)
151    
152        def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):
153            """
154            returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}
155            on locations where C{location_of_fixed_pressure} is positive (see L{setValue}).
156            Note that C{g} is used, see L{setValue}.
157            
158            @param v: flux.
159            @type v: vector-valued on the domain (e.g. L{Data}).
160            @param fixed_pressure: pressure on the locations of the domain marked be C{location_of_fixed_pressure}.
161            @type fixed_pressure: vector values on the domain (e.g. L{Data}).
162            @param tol: relative tolerance to be used.
163            @type tol: positive C{float}.
164            @return: pressure
165            @rtype: L{Data}
166            @note: the method uses the least squares solution M{p=(Q^*Q)^{-1}Q^*(g-u)} where and M{(Qp)_i=k_{ij}p_{,j}}
167                   for the permeability M{k_{ij}}
168            """
169            self.__pde_v.setTolerance(tol)
170            g=self.__g
171            self.__pde_p.setValue(r=fixed_pressure)
172            if v == None:
173               self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-v))
174            else:
175               self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g))
176            return self.__pde_p.getSolution(verbose=show_details)
177    
178        def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None):
179            """
180            set the tolerance C{ATOL} used to terminate the solution process. It is used
181    
182            M{ATOL = atol + rtol * max( |g-v_ref|, |Qp_ref| )}
183    
184            where M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}. If C{v_ref} or C{p_ref} is not present zero is assumed.
185    
186            The iteration is terminated if for the current approximation C{p}, flux C{v=(I+D^*D)^{-1}(D^*f-g-Qp)} and their residual
187    
188            M{r=Q^*(g-Qp-v)}
189    
190            the condition
191    
192            M{<(Q^*Q)^{-1} r,r> <= ATOL}
193    
194            holds. M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}
195    
196            @param atol: absolute tolerance for the pressure
197            @type atol: non-negative C{float}
198            @param rtol: relative tolerance for the pressure
199            @type rtol: non-negative C{float}
200            @param p_ref: reference pressure. If not present zero is used. You may use physical arguments to set a resonable value for C{p_ref}, use the
201            L{getPressure} method or use  the value from a previous time step.
202            @type p_ref: scalar value on the domain (e.g. L{Data}).
203            @param v_ref: reference velocity.  If not present zero is used. You may use physical arguments to set a resonable value for C{v_ref}, use the
204            L{getFlux} method or use  the value from a previous time step.
205            @type v_ref: vector-valued on the domain (e.g. L{Data}).
206            @return: used absolute tolerance.
207            @rtype: positive C{float}
208            """
209            g=self.__g
210            if not v_ref == None:
211               f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)
212            else:
213               f1=util.integrate(util.length(util.interpolate(g))**2)
214            if not p_ref == None:
215               f2=util.integrate(util.length(util.tensor_mult(self.__permeability,util.grad(p_ref)))**2)
216            else:
217               f2=0
218            self.__ATOL= atol + rtol * util.sqrt(max(f1,f2))
219            if self.__ATOL<=0:
220               raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL
221            return self.__ATOL
222            
223        def getTolerance(self):
224            """
225            returns the current tolerance.
226      
227            @return: used absolute tolerance.
228            @rtype: positive C{float}
229            """
230            if self.__ATOL==None:
231               raise ValueError,"no tolerance is defined."
232            return self.__ATOL
233    
234        def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):
235             """
236             solves the problem.
237    
238             The iteration is terminated if the residual norm is less then self.getTolerance().
239    
240             @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.
241             @type u0: vector value on the domain (e.g. L{Data}).
242             @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.
243             @type p0: scalar value on the domain (e.g. L{Data}).
244             @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}
245             @type sub_rtol: positive-negative C{float}
246             @param verbose: if set some information on iteration progress are printed
247             @type verbose: C{bool}
248             @param show_details:  if set information on the subiteration process are printed.
249             @type show_details: C{bool}
250             @return: flux and pressure
251             @rtype: C{tuple} of L{Data}.
252    
253             @note: The problem is solved as a least squares form
254    
255             M{(I+D^*D)u+Qp=D^*f+g}
256             M{Q^*u+Q^*Qp=Q^*g}
257    
258             where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
259             We eliminate the flux form the problem by setting
260    
261             M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
262    
263             form the first equation. Inserted into the second equation we get
264    
265             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
266            
267             which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
268             PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
269           """           """
270           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           self.verbose=verbose
271           """           self.show_details= show_details and self.verbose
272           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_v.setTolerance(sub_rtol)
273           self.__pde_u.setValue(X=-self.getStress(u),X_reduced=-p*util.kronecker(self.domain))           self.__pde_p.setTolerance(sub_rtol)
274           return  self.__pde_u.getSolution(verbose=self.show_details)           ATOL=self.getTolerance()
275             if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL
276             #########################################################################################################################
277        def solve_prec(self,p):           #
278          a=self.__inv_eta*p           #   we solve:
279          return a-util.integrate(a)/self.vol           #  
280             #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )
281        def stoppingcriterium(self,Bv,v,p):           #
282            n_r=util.sqrt(self.inner(Bv,Bv))           #   residual is
283            n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))           #
284            if self.verbose: print "PCG step %s: L2(div(v)) = %s, L2(grad(v))=%s"%(self.iter,n_r,n_v) , util.Lsup(v)           #    r=  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) - Q dp +(I+D^*D)^{-1})Q dp ) = Q^* (g - Qp - v)
285            if self.iter == 0: self.__n_v=n_v;           #
286            self.__n_v, n_v_old =n_v, self.__n_v           #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC
287            self.iter+=1           #
288            if self.iter>1 and n_r <= n_v*self.getTolerance() and abs(n_v_old-self.__n_v) <= n_v * self.getTolerance():           #    we use (g - Qp, v) to represent the residual. not that
289                if self.verbose: print "PCG terminated after %s steps."%self.iter           #
290                return True           #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)
291            else:           #
292                return False           #   while the initial residual is
293        def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):           #
294        if TOL==None:           #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC
295               TOL=self.getTolerance()           #  
296            if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)           d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))
297            self.iter+=1           self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)
298                       v00=self.__pde_v.getSolution(verbose=show_details)
299            if norm_r <= norm_b*TOL:           if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)
300                if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)           self.__pde_v.setValue(r=Data())
301                return True           # start CG
302            else:           r=ArithmeticTuple(d0, v00)
303                return False           p,r=PCG(r,self.__Aprod_PCG,p0,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
304             return r[1],p
305    
306        def __Aprod_PCG(self,dp):
307              if self.show_details: print "DarcyFlux: Applying operator"
308              #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)
309              mQdp=util.tensor_mult(self.__permeability,util.grad(dp))
310              self.__pde_v.setValue(Y=mQdp,X=Data(), r=Data())
311              du=self.__pde_v.getSolution(verbose=self.show_details)
312              return ArithmeticTuple(mQdp,du)
313    
314        def __inner_PCG(self,p,r):
315             a=util.tensor_mult(self.__permeability,util.grad(p))
316             f0=util.integrate(util.inner(a,r[0]))
317             f1=util.integrate(util.inner(a,r[1]))
318             # print "__inner_PCG:",f0,f1,"->",f0-f1
319             return f0-f1
320    
321        def __Msolve_PCG(self,r):
322              if self.show_details: print "DarcyFlux: Applying preconditioner"
323              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1]), r=Data())
324              return self.__pde_p.getSolution(verbose=self.show_details)
325    
326  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
327        """        """
328        solves        solves
329    
330            -(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}
331                  u_{i,i}=0                  u_{i,i}=0
332    
333            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
334            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
335    
336        if surface_stress is not give 0 is assumed.        if surface_stress is not given 0 is assumed.
337    
338        typical usage:        typical usage:
339    
# Line 158  class StokesProblemCartesian(Homogeneous Line 343  class StokesProblemCartesian(Homogeneous
343              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
344        """        """
345        def __init__(self,domain,**kwargs):        def __init__(self,domain,**kwargs):
346             """
347             initialize the Stokes Problem
348    
349             @param domain: domain of the problem. The approximation order needs to be two.
350             @type domain: L{Domain}
351             @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
352             """
353           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
354           self.domain=domain           self.domain=domain
355           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
356           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())
357           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
358           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
359             # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)
360                            
361           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
362           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
363             # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
364           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
365    
366           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
# Line 174  class StokesProblemCartesian(Homogeneous Line 368  class StokesProblemCartesian(Homogeneous
368           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
369           self.__pde_proj.setValue(D=1.)           self.__pde_proj.setValue(D=1.)
370    
371        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()):
372            """
373            assigns values to the model parameters
374    
375            @param f: external force
376            @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
377            @param fixed_u_mask: mask of locations with fixed velocity.
378            @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
379            @param eta: viscosity
380            @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
381            @param surface_stress: normal surface stress
382            @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
383            @param stress: initial stress
384        @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
385            @note: All values needs to be set.
386    
387            """
388          self.eta=eta          self.eta=eta
389          A =self.__pde_u.createCoefficientOfGeneralPDE("A")          A =self.__pde_u.createCoefficient("A")
390      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
391          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
392          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
# Line 184  class StokesProblemCartesian(Homogeneous Line 394  class StokesProblemCartesian(Homogeneous
394              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
395      self.__pde_prec.setValue(D=1/self.eta)      self.__pde_prec.setValue(D=1/self.eta)
396          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)
397            self.__stress=stress
398    
399          def B(self,v):
400            """
401            returns div(v)
402            @rtype: equal to the type of p
403    
404            @note: boundary conditions on p should be zero!
405            """
406            if self.show_details: print "apply divergence:"
407            self.__pde_proj.setValue(Y=-util.div(v))
408            self.__pde_proj.setTolerance(self.getSubProblemTolerance())
409            return self.__pde_proj.getSolution(verbose=self.show_details)
410    
411          def inner_pBv(self,p,Bv):
412             """
413             returns inner product of element p and Bv  (overwrite)
414            
415             @type p: equal to the type of p
416             @type Bv: equal to the type of result of operator B
417             @rtype: C{float}
418    
419        def B(self,arg):           @rtype: equal to the type of p
420           d=util.div(arg)           """
421           self.__pde_proj.setValue(Y=d)           s0=util.interpolate(p,Function(self.domain))
422           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))  
423           return util.integrate(s0*s1)           return util.integrate(s0*s1)
424    
425        def inner_a(self,a0,a1):        def inner_p(self,p0,p1):
426           p0=util.interpolate(a0[1],Function(self.domain))           """
427           p1=util.interpolate(a1[1],Function(self.domain))           returns inner product of element p0 and p1  (overwrite)
428           alfa=(1/self.vol)*util.integrate(p0)          
429           beta=(1/self.vol)*util.integrate(p1)           @type p0: equal to the type of p
430       v0=util.grad(a0[0])           @type p1: equal to the type of p
431       v1=util.grad(a1[0])           @rtype: C{float}
          return util.integrate((p0-alfa)*(p1-beta)+((1/self.eta)**2)*util.inner(v0,v1))  
432    
433             @rtype: equal to the type of p
434             """
435             s0=util.interpolate(p0/self.eta,Function(self.domain))
436             s1=util.interpolate(p1/self.eta,Function(self.domain))
437             return util.integrate(s0*s1)
438    
439        def getStress(self,u):        def inner_v(self,v0,v1):
440           mg=util.grad(u)           """
441           return 2.*self.eta*util.symmetric(mg)           returns inner product of two element v0 and v1  (overwrite)
442        def getEtaEffective(self):          
443           return self.eta           @type v0: equal to the type of v
444             @type v1: equal to the type of v
445             @rtype: C{float}
446    
447             @rtype: equal to the type of v
448             """
449         gv0=util.grad(v0)
450         gv1=util.grad(v1)
451             return util.integrate(util.inner(gv0,gv1))
452    
453        def solve_A(self,u,p):        def solve_A(self,u,p):
454           """           """
455           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
456           """           """
457             if self.show_details: print "solve for velocity:"
458           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
459           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))           if self.__stress.isEmpty():
460           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))
461             else:
462                self.__pde_u.setValue(X=self.__stress-2*self.eta*util.symmetric(util.grad(u))+p*util.kronecker(self.domain))
463             out=self.__pde_u.getSolution(verbose=self.show_details)
464             return  out
465    
466        def solve_prec(self,p):        def solve_prec(self,p):
467       #proj=Projector(domain=self.domain, reduce = True, fast=False)           if self.show_details: print "apply preconditioner:"
468           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.__pde_prec.setTolerance(self.getSubProblemTolerance())
469           self.__pde_prec.setValue(Y=p)           self.__pde_prec.setValue(Y=p)
470           q=self.__pde_prec.getSolution(verbose=self.show_details)           q=self.__pde_prec.getSolution(verbose=self.show_details)
471           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.1673  
changed lines
  Added in v.2208

  ViewVC Help
Powered by ViewVC 1.1.26