/[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 2625 by jfenwick, Fri Aug 21 06:30:25 2009 UTC revision 3074 by gross, Tue Jul 27 01:47:45 2010 UTC
# Line 1  Line 1 
1    # -*- coding: utf-8 -*-
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 10  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2009 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# Line 37  from linearPDEs import LinearPDE, Linear Line 38  from linearPDEs import LinearPDE, Linear
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39    
40  class DarcyFlow(object):  class DarcyFlow(object):
41       """
42       solves the problem
43      
44       *u_i+k_{ij}*p_{,j} = g_i*
45       *u_{i,i} = f*
46      
47       where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
48      
49       :note: The problem is solved in a least squares formulation.
50       """
51      
52       def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False):
53          """
54          initializes the Darcy flux problem
55          :param domain: domain of the problem
56          :type domain: `Domain`
57          :param useReduced: uses reduced oreder on flux and pressure
58          :type useReduced: ``bool``
59          :param adaptSubTolerance: switches on automatic subtolerance selection
60          :type adaptSubTolerance: ``bool``
61          :param solveForFlux: if True the solver solves for the flux (do not use!)
62          :type solveForFlux: ``bool``  
63          """
64          self.domain=domain
65          self.solveForFlux=solveForFlux
66          self.useReduced=useReduced
67          self.__adaptSubTolerance=adaptSubTolerance
68          self.verbose=False
69          
70          self.__pde_k=LinearPDESystem(domain)
71          self.__pde_k.setSymmetryOn()
72          if self.useReduced: self.__pde_k.setReducedOrderOn()
73    
74          self.__pde_p=LinearSinglePDE(domain)
75          self.__pde_p.setSymmetryOn()
76          if self.useReduced: self.__pde_p.setReducedOrderOn()
77    
78          self.__pde_l=LinearSinglePDE(domain)   # this is here for getSolverOptionsWeighting
79          # self.__pde_l.setSymmetryOn()
80          # if self.useReduced: self.__pde_l.setReducedOrderOn()
81          self.setTolerance()
82          self.setAbsoluteTolerance()
83          self.__f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
84          self.__g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
85          
86       def getSolverOptionsFlux(self):
87          """
88          Returns the solver options used to solve the flux problems
89          
90          *K^{-1} u=F*
91          
92          :return: `SolverOptions`
93          """
94          return self.__pde_k.getSolverOptions()
95          
96       def setSolverOptionsFlux(self, options=None):
97          """
98          Sets the solver options used to solve the flux problems
99          
100          *K^{-1}u=F*
101          
102          If ``options`` is not present, the options are reset to default
103          
104          :param options: `SolverOptions`
105          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
106          """
107          return self.__pde_v.setSolverOptions(options)
108        
109       def getSolverOptionsPressure(self):
110          """
111          Returns the solver options used to solve the pressure problems
112          
113          *(Q^* K Q)p=-Q^*G*
114          
115          :return: `SolverOptions`
116          """
117          return self.__pde_p.getSolverOptions()
118          
119       def setSolverOptionsPressure(self, options=None):
120          """
121          Sets the solver options used to solve the pressure problems
122          
123          *(Q^* K Q)p=-Q^*G*
124          
125          If ``options`` is not present, the options are reset to default
126          
127          :param options: `SolverOptions`
128          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
129          """
130          return self.__pde_p.setSolverOptions(options)
131    
132       def getSolverOptionsWeighting(self):
133          """
134          Returns the solver options used to solve the pressure problems
135    
136          *(D K D^*)p=-D F*
137    
138          :return: `SolverOptions`
139          """
140          return self.__pde_l.getSolverOptions()
141    
142       def setSolverOptionsWeighting(self, options=None):
143          """
144          Sets the solver options used to solve the pressure problems
145    
146          *(D K D^*)p=-D F*
147    
148          If ``options`` is not present, the options are reset to default
149    
150          :param options: `SolverOptions`
151          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
152          """
153          return self.__pde_l.setSolverOptions(options)
154    
155    
156       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
157          """
158          assigns values to model parameters
159          
160          :param f: volumetic sources/sinks
161          :type f: scalar value on the domain (e.g. `Data`)
162          :param g: flux sources/sinks
163          :type g: vector values on the domain (e.g. `Data`)
164          :param location_of_fixed_pressure: mask for locations where pressure is fixed
165          :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)
166          :param location_of_fixed_flux:  mask for locations where flux is fixed.
167          :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)
168          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
169          ``s`` on the main diagonal is used.
170          :type permeability: scalar or tensor values on the domain (e.g. `Data`)
171          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
172          :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
173          or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
174          is along the *x_i* axis.
175          """
176          if f !=None:
177         f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("X"))
178         if f.isEmpty():
179            f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
180         else:
181            if f.getRank()>0: raise ValueError,"illegal rank of f."
182            self.__f=f
183          if g !=None:
184         g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
185         if g.isEmpty():
186            g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187         else:
188            if not g.getShape()==(self.domain.getDim(),):
189               raise ValueError,"illegal shape of g"
190            self.__g=g
191          if location_of_fixed_pressure!=None:
192               self.__pde_p.setValue(q=location_of_fixed_pressure)
193               #self.__pde_l.setValue(q=location_of_fixed_pressure)
194          if location_of_fixed_flux!=None:
195               self.__pde_k.setValue(q=location_of_fixed_flux)
196                
197          if permeability!=None:
198         perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
199         if perm.getRank()==0:
200            perm_inv=(1./perm)*util.kronecker(self.domain.getDim())
201            perm=perm*util.kronecker(self.domain.getDim())
202         elif perm.getRank()==2:
203            perm_inv=util.inverse(perm)
204         else:
205            raise ValueError,"illegal rank of permeability."
206    
207         self.__permeability=perm
208         self.__permeability_inv=perm_inv
209         self.__l =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))/10
210         #self.__l =(self.domain.getSize()**2*util.length(self.__permeability_inv))/10
211         if  self.solveForFlux:
212            self.__pde_k.setValue(D=self.__permeability_inv)
213         else:
214            self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))
215         self.__pde_p.setValue(A=self.__permeability)
216         #self.__pde_l.setValue(D=1/self.__l)
217             #self.__pde_l.setValue(A=self.__permeability)
218    
219       def __applWeight(self, v, f=None):
220          # solves L p = f-Dv with p = 0
221          if self.getSolverOptionsWeighting().isVerbose() or self.verbose: print "DarcyFlux: Applying weighting operator"
222          if f == None:
223         return -util.div(v)*self.__l
224          else:
225         return (f-util.div(v))*self.__l
226          # if f == None:
227          #      self.__pde_l.setValue(Y=-util.div(v))  
228          # else:
229          #      return (f-util.div(v))/self.__l
230          # return self.__pde_l.getSolution()
231          
232       def __getPressure(self, v, p0, g=None):
233          # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0
234          if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"
235          if g == None:
236         self.__pde_p.setValue(X=-v, r=p0)
237          else:
238         self.__pde_p.setValue(X=g-v, r=p0)      
239          p=self.__pde_p.getSolution()
240          return p
241    
242       def __Aprod_v(self,dv):
243          # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv  
244          dp=self.__getPressure(dv, p0=Data()) # dp = (G*KG)^{-1} (0-G^*dv)
245          a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp
246          b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)
247          return ArithmeticTuple(a,b,-dp)
248    
249       def __Msolve_PCG_v(self,r):
250          # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv
251          if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"
252          self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=Data())
253          # self.__pde_p.getOperator().saveMM("prec.mm")
254          return self.__pde_k.getSolution()
255    
256       def __inner_PCG_v(self,v,r):
257          return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])
258          
259       def __Aprod_p(self,dp):
260          if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
261          Gdp=util.grad(dp)
262          self.__pde_k.setValue(Y=-Gdp,X=Data(), r=Data())
263          du=self.__pde_k.getSolution()
264          # self.__pde_v.getOperator().saveMM("proj.mm")
265          return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)
266    
267       def __getFlux(self,p, v0, f=None, g=None):
268          # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp
269          if f!=None:
270         self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))
271          self.__pde_k.setValue(r=v0)
272          g2=util.tensor_mult(self.__permeability_inv,g)
273          if p == None:
274         self.__pde_k.setValue(Y=g2)
275          else:
276         self.__pde_k.setValue(Y=g2-util.grad(p))
277          return self.__pde_k.getSolution()  
278          
279          #v=self.__getFlux(p, u0, f=self.__f, g=g2)      
280       def __Msolve_PCG_p(self,r):
281          if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
282          self.__pde_p.setValue(X=r[0]-r[1], Y=Data(), r=Data(), y=Data())
283          # self.__pde_p.getOperator().saveMM("prec.mm")
284          return self.__pde_p.getSolution()
285            
286       def __inner_PCG_p(self,p,r):
287           return util.integrate(util.inner(util.grad(p), r[0]-r[1]))
288    
289       def __L2(self,v):
290          return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
291    
292       def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
293          """
294          solves the problem.
295          
296          The iteration is terminated if the residual norm is less then self.getTolerance().
297    
298          :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
299          :type u0: vector value on the domain (e.g. `Data`).
300          :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
301          :type p0: scalar value on the domain (e.g. `Data`).
302          :param verbose: if set some information on iteration progress are printed
303          :type verbose: ``bool``
304          :return: flux and pressure
305          :rtype: ``tuple`` of `Data`.
306    
307          :note: The problem is solved as a least squares form
308          *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*
309          *G^*u+G^* K Gp=G^*g*
310    
311          where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.
312          """
313          self.verbose=verbose
314          rtol=self.getTolerance()
315          atol=self.getAbsoluteTolerance()
316          self.setSubProblemTolerance()
317          num_corrections=0
318          converged=False
319          norm_r=None
320          
321          # Eliminate the hydrostatic pressure:
322          if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."
323          self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))        
324          p0=self.__pde_p.getSolution()
325          g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0))
326          norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5    
327    
328          p=p0*0
329          if self.solveForFlux:
330         v=u0.copy()
331          else:
332         v=self.__getFlux(p, u0, f=self.__f, g=g2)
333    
334          while not converged and norm_g2 > 0:
335         Gp=util.grad(p)
336         KGp=util.tensor_mult(self.__permeability,Gp)
337         if self.verbose:
338            def_p=g2-(v+KGp)
339            def_v=self.__f-util.div(v)
340            print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))
341            print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v)))
342            print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5
343            print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2
344            print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5
345         ATOL=atol+rtol*norm_g2
346         if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)
347         if norm_r == None or norm_r>ATOL:
348            if num_corrections>max_num_corrections:
349               raise ValueError,"maximum number of correction steps reached."
350          
351            if self.solveForFlux:
352               # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)
353               v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),
354                   self.__Aprod_v,
355                   v,
356                   self.__Msolve_PCG_v,
357                   self.__inner_PCG_v,
358                   atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
359               p=r[2]
360            else:
361               # initial residual is r=G^*(g2-KGp - v)
362               p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),
363                     self.__Aprod_p,
364                     p,
365                     self.__Msolve_PCG_p,
366                     self.__inner_PCG_p,
367                     atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
368               v=r[1]
369            if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r
370            num_corrections+=1
371         else:
372            if self.verbose: print "DarcyFlux: stopping criterium reached."
373            converged=True
374          return v,p+p0
375       def setTolerance(self,rtol=1e-4):
376          """
377          sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
378    
379          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
380          
381          where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).
382          
383          :param rtol: relative tolerance for the pressure
384          :type rtol: non-negative ``float``
385          """
386          if rtol<0:
387         raise ValueError,"Relative tolerance needs to be non-negative."
388          self.__rtol=rtol
389       def getTolerance(self):
390          """
391          returns the relative tolerance
392          :return: current relative tolerance
393          :rtype: ``float``
394          """
395          return self.__rtol
396    
397       def setAbsoluteTolerance(self,atol=0.):
398          """
399          sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
400          
401          *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*
402    
403    
404          where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
405    
406          :param atol: absolute tolerance for the pressure
407          :type atol: non-negative ``float``
408          """
409          if atol<0:
410         raise ValueError,"Absolute tolerance needs to be non-negative."
411          self.__atol=atol
412       def getAbsoluteTolerance(self):
413          """
414          returns the absolute tolerance
415          :return: current absolute tolerance
416          :rtype: ``float``
417          """
418          return self.__atol
419       def getSubProblemTolerance(self):
420          """
421          Returns a suitable subtolerance
422          :type: ``float``
423          """
424          return max(util.EPSILON**(0.5),self.getTolerance()**2)
425    
426       def setSubProblemTolerance(self):
427          """
428          Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
429          """
430          if self.__adaptSubTolerance:
431         sub_tol=self.getSubProblemTolerance()
432         self.getSolverOptionsFlux().setTolerance(sub_tol)
433         self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
434         self.getSolverOptionsPressure().setTolerance(sub_tol)
435         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
436         self.getSolverOptionsWeighting().setTolerance(sub_tol)
437         self.getSolverOptionsWeighting().setAbsoluteTolerance(0.)
438         if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
439    
440    
441    class DarcyFlowOld(object):
442      """      """
443      solves the problem      solves the problem
444    
# Line 63  class DarcyFlow(object): Line 465  class DarcyFlow(object):
465             s=self.domain.getSize()             s=self.domain.getSize()
466             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
467             # self.__l=(3.*util.longestEdge(self.domain))**2             # self.__l=(3.*util.longestEdge(self.domain))**2
468             # self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2             #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
469          else:          else:
470             self.__l=weight             self.__l=weight
471          self.__pde_v=LinearPDESystem(domain)          self.__pde_v=LinearPDESystem(domain)
# Line 274  class DarcyFlow(object): Line 676  class DarcyFlow(object):
676           rtol=self.getTolerance()           rtol=self.getTolerance()
677           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
678       self.setSubProblemTolerance()       self.setSubProblemTolerance()
       
679           num_corrections=0           num_corrections=0
680           converged=False           converged=False
681           p=p0           p=p0
# Line 297  class DarcyFlow(object): Line 698  class DarcyFlow(object):
698                 ATOL=(atol+rtol*fac)                 ATOL=(atol+rtol*fac)
699                 if self.verbose:                 if self.verbose:
700                      print "DarcyFlux: L2 norm of v = %e."%norm_v                      print "DarcyFlux: L2 norm of v = %e."%norm_v
701                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp                      print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp
702                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)
703                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
704                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
# Line 378  class StokesProblemCartesian(Homogeneous Line 779  class StokesProblemCartesian(Homogeneous
779              sp.initialize(...)              sp.initialize(...)
780              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
781       """       """
782       def __init__(self,domain,adaptSubTolerance=True, **kwargs):       def __init__(self,domain,**kwargs):
783           """           """
784           initialize the Stokes Problem           initialize the Stokes Problem
785    
786           :param domain: domain of the problem. The approximation order needs to be two.           The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
787             LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
788             with macro elements for the pressure.
789    
790             :param domain: domain of the problem.
791           :type domain: `Domain`           :type domain: `Domain`
      :param adaptSubTolerance: If True the tolerance for subproblem is set automatically.  
      :type adaptSubTolerance: ``bool``  
          :warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.  
792           """           """
793           HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
794           self.domain=domain           self.domain=domain
          self.vol=util.integrate(1.,Function(self.domain))  
795           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())
796           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
797            
# Line 440  class StokesProblemCartesian(Homogeneous Line 841  class StokesProblemCartesian(Homogeneous
841       :param options: new solver options       :param options: new solver options
842       :type options: `SolverOptions`       :type options: `SolverOptions`
843       """       """
844       self.__pde_prec.setSolverOptions(options)       self.__pde_proj.setSolverOptions(options)
845       def getSolverOptionsDiv(self):       def getSolverOptionsDiv(self):
846           """           """
847       returns the solver options for solving the equation to project the divergence of       returns the solver options for solving the equation to project the divergence of
# Line 448  class StokesProblemCartesian(Homogeneous Line 849  class StokesProblemCartesian(Homogeneous
849            
850       :rtype: `SolverOptions`       :rtype: `SolverOptions`
851       """       """
852       return self.__pde_prec.getSolverOptions()       return self.__pde_proj.getSolverOptions()
853       def setSubProblemTolerance(self):  
854         def updateStokesEquation(self, v, p):
855           """           """
856       Updates the tolerance for subproblems           updates the Stokes equation to consider dependencies from ``v`` and ``p``
857             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
858           """           """
859       if self.adaptSubTolerance():           pass
860               sub_tol=self.getSubProblemTolerance()       def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
861           self.getSolverOptionsDiv().setTolerance(sub_tol)          """
862           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)          assigns new values to the model parameters.
863           self.getSolverOptionsPressure().setTolerance(sub_tol)  
864           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)          :param f: external force
865           self.getSolverOptionsVelocity().setTolerance(sub_tol)          :type f: `Vector` object in `FunctionSpace` `Function` or similar
866           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)          :param fixed_u_mask: mask of locations with fixed velocity.
867                    :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
868            :param eta: viscosity
869            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
870            :param surface_stress: normal surface stress
871            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
872            :param stress: initial stress
873        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
874            """
875            if eta !=None:
876                k=util.kronecker(self.domain.getDim())
877                kk=util.outer(k,k)
878                self.eta=util.interpolate(eta, Function(self.domain))
879            self.__pde_prec.setValue(D=1/self.eta)
880                self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
881            if restoration_factor!=None:
882                n=self.domain.getNormal()
883                self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))
884            if fixed_u_mask!=None:
885                self.__pde_u.setValue(q=fixed_u_mask)
886            if f!=None: self.__f=f
887            if surface_stress!=None: self.__surface_stress=surface_stress
888            if stress!=None: self.__stress=stress
889    
890       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
891          """          """
# Line 477  class StokesProblemCartesian(Homogeneous Line 901  class StokesProblemCartesian(Homogeneous
901          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
902          :param stress: initial stress          :param stress: initial stress
903      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
         :note: All values needs to be set.  
904          """          """
905          self.eta=eta          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
         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.  
         n=self.domain.getNormal()  
     self.__pde_prec.setValue(D=1/self.eta)  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask, d=restoration_factor*util.outer(n,n))  
         self.__f=f  
         self.__surface_stress=surface_stress  
         self.__stress=stress  
906    
907       def Bv(self,v):       def Bv(self,v,tol):
908           """           """
909           returns inner product of element p and div(v)           returns inner product of element p and div(v)
910    
# Line 501  class StokesProblemCartesian(Homogeneous Line 912  class StokesProblemCartesian(Homogeneous
912           :return: inner product of element p and div(v)           :return: inner product of element p and div(v)
913           :rtype: ``float``           :rtype: ``float``
914           """           """
915           self.__pde_proj.setValue(Y=-util.div(v))           self.__pde_proj.setValue(Y=-util.div(v))
916           return self.__pde_proj.getSolution()       self.getSolverOptionsDiv().setTolerance(tol)
917         self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
918             out=self.__pde_proj.getSolution()
919             return out
920    
921       def inner_pBv(self,p,Bv):       def inner_pBv(self,p,Bv):
922           """           """
# Line 524  class StokesProblemCartesian(Homogeneous Line 938  class StokesProblemCartesian(Homogeneous
938           :return: inner product of p0 and p1           :return: inner product of p0 and p1
939           :rtype: ``float``           :rtype: ``float``
940           """           """
941           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0,Function(self.domain))
942           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1,Function(self.domain))
943           return util.integrate(s0*s1)           return util.integrate(s0*s1)
944    
945       def norm_v(self,v):       def norm_v(self,v):
# Line 536  class StokesProblemCartesian(Homogeneous Line 950  class StokesProblemCartesian(Homogeneous
950           :return: norm of v           :return: norm of v
951           :rtype: non-negative ``float``           :rtype: non-negative ``float``
952           """           """
953           return util.sqrt(util.integrate(util.length(util.grad(v))))           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
954    
955    
956       def getV(self, p, v0):       def getDV(self, p, v, tol):
957           """           """
958           return the value for v for a given p (overwrite)           return the value for v for a given p (overwrite)
959    
960           :param p: a pressure           :param p: a pressure
961           :param v0: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
962           :return: v given as *v= A^{-1} (f-B^*p)*           :return: dv given as *Adv=(f-Av-B^*p)*
963           """           """
964           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)           self.updateStokesEquation(v,p)
965             self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
966         self.getSolverOptionsVelocity().setTolerance(tol)
967         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
968           if self.__stress.isEmpty():           if self.__stress.isEmpty():
969              self.__pde_u.setValue(X=p*util.kronecker(self.domain))              self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
970           else:           else:
971              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
972           out=self.__pde_u.getSolution()           out=self.__pde_u.getSolution()
973           return  out           return  out
974    
# Line 563  class StokesProblemCartesian(Homogeneous Line 981  class StokesProblemCartesian(Homogeneous
981          """          """
982          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
983    
984       def solve_AinvBt(self,p):       def solve_AinvBt(self,p, tol):
985           """           """
986           Solves *Av=B^*p* with accuracy `self.getSubProblemTolerance()`           Solves *Av=B^*p* with accuracy `tol`
987    
988           :param p: a pressure increment           :param p: a pressure increment
989           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
990           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
991           """           """
992           self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))
993           out=self.__pde_u.getSolution()           out=self.__pde_u.getSolution()
994           return  out           return  out
995    
996       def solve_prec(self,Bv):       def solve_prec(self,Bv, tol):
997           """           """
998           applies preconditioner for for *BA^{-1}B^** to *Bv*           applies preconditioner for for *BA^{-1}B^** to *Bv*
999           with accuracy `self.getSubProblemTolerance()`           with accuracy `self.getSubProblemTolerance()`
# Line 585  class StokesProblemCartesian(Homogeneous Line 1003  class StokesProblemCartesian(Homogeneous
1003           :note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
1004           """           """
1005           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)
1006           return self.__pde_prec.getSolution()       self.getSolverOptionsPressure().setTolerance(tol)
1007         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
1008             out=self.__pde_prec.getSolution()
1009             return out

Legend:
Removed from v.2625  
changed lines
  Added in v.3074

  ViewVC Help
Powered by ViewVC 1.1.26