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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2620 - (hide annotations)
Thu Aug 20 06:24:00 2009 UTC (10 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 23389 byte(s)
some small additions to pycad to make life a bit easier.
1 ksteube 1809 ########################################################
2 gross 1414 #
3 jfenwick 2548 # Copyright (c) 2003-2009 by University of Queensland
4 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
5     # http://www.uq.edu.au/esscc
6 gross 1414 #
7 ksteube 1809 # Primary Business: Queensland, Australia
8     # Licensed under the Open Software License version 3.0
9     # http://www.opensource.org/licenses/osl-3.0.php
10 gross 1414 #
11 ksteube 1809 ########################################################
12 gross 1414
13 jfenwick 2549 __copyright__="""Copyright (c) 2003-2009 by University of Queensland
14 ksteube 1809 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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
20 ksteube 1809
21 gross 1414 """
22     Some models for flow
23    
24     @var __author__: name of author
25     @var __copyright__: copyrights
26     @var __license__: licence agreement
27     @var __url__: url entry point on documentation
28     @var __version__: version
29     @var __date__: date of the version
30     """
31    
32     __author__="Lutz Gross, l.gross@uq.edu.au"
33    
34     from escript import *
35     import util
36 gross 2474 from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
37 gross 2264 from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
38 gross 1414
39 gross 2100 class DarcyFlow(object):
40     """
41 gross 2264 solves the problem
42 gross 1659
43 gross 2100 M{u_i+k_{ij}*p_{,j} = g_i}
44     M{u_{i,i} = f}
45 gross 1659
46 gross 2264 where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,
47 gross 1659
48 gross 2100 @note: The problem is solved in a least squares formulation.
49     """
50 gross 1659
51 gross 2474 def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):
52 gross 2100 """
53 gross 2208 initializes the Darcy flux problem
54 gross 2100 @param domain: domain of the problem
55     @type domain: L{Domain}
56 gross 2474 @param useReduced: uses reduced oreder on flux and pressure
57     @type useReduced: C{bool}
58     @param adaptSubTolerance: switches on automatic subtolerance selection
59     @type adaptSubTolerance: C{bool}
60 gross 2100 """
61     self.domain=domain
62 gross 2351 if weight == None:
63 gross 2354 s=self.domain.getSize()
64     self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
65 gross 2486 # self.__l=(3.*util.longestEdge(self.domain))**2
66     # self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2
67 gross 2351 else:
68     self.__l=weight
69 gross 2100 self.__pde_v=LinearPDESystem(domain)
70 gross 2208 if useReduced: self.__pde_v.setReducedOrderOn()
71     self.__pde_v.setSymmetryOn()
72 gross 2267 self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
73 gross 2100 self.__pde_p=LinearSinglePDE(domain)
74     self.__pde_p.setSymmetryOn()
75 gross 2208 if useReduced: self.__pde_p.setReducedOrderOn()
76 gross 2100 self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
77     self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
78 gross 2264 self.setTolerance()
79     self.setAbsoluteTolerance()
80 gross 2474 self.__adaptSubTolerance=adaptSubTolerance
81     self.verbose=False
82     def getSolverOptionsFlux(self):
83     """
84     Returns the solver options used to solve the flux problems
85    
86     M{(I+D^*D)u=F}
87    
88     @return: L{SolverOptions}
89     """
90     return self.__pde_v.getSolverOptions()
91     def setSolverOptionsFlux(self, options=None):
92     """
93     Sets the solver options used to solve the flux problems
94    
95     M{(I+D^*D)u=F}
96    
97     If C{options} is not present, the options are reset to default
98     @param options: L{SolverOptions}
99     @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.
100     """
101     return self.__pde_v.setSolverOptions(options)
102     def getSolverOptionsPressure(self):
103     """
104     Returns the solver options used to solve the pressure problems
105    
106     M{(Q^*Q)p=Q^*G}
107    
108     @return: L{SolverOptions}
109     """
110     return self.__pde_p.getSolverOptions()
111     def setSolverOptionsPressure(self, options=None):
112     """
113     Sets the solver options used to solve the pressure problems
114    
115     M{(Q^*Q)p=Q^*G}
116    
117     If C{options} is not present, the options are reset to default
118     @param options: L{SolverOptions}
119     @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.
120     """
121     return self.__pde_p.setSolverOptions(options)
122 gross 1659
123 gross 2100 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
124     """
125 gross 2208 assigns values to model parameters
126 gross 1659
127 gross 2208 @param f: volumetic sources/sinks
128     @type f: scalar value on the domain (e.g. L{Data})
129 gross 2100 @param g: flux sources/sinks
130 gross 2208 @type g: vector values on the domain (e.g. L{Data})
131 gross 2100 @param location_of_fixed_pressure: mask for locations where pressure is fixed
132 gross 2208 @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
133     @param location_of_fixed_flux: mask for locations where flux is fixed.
134     @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
135 gross 2264 @param permeability: permeability tensor. If scalar C{s} is given the tensor with
136     C{s} on the main diagonal is used. If vector C{v} is given the tensor with
137 gross 2208 C{v} on the main diagonal is used.
138     @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
139 gross 1659
140 gross 2208 @note: the values of parameters which are not set by calling C{setValue} are not altered.
141     @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
142     or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
143     is along the M{x_i} axis.
144 gross 2100 """
145 gross 2264 if f !=None:
146 gross 2100 f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
147     if f.isEmpty():
148     f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
149     else:
150     if f.getRank()>0: raise ValueError,"illegal rank of f."
151 gross 2267 self.__f=f
152 gross 2264 if g !=None:
153 gross 2100 g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
154     if g.isEmpty():
155     g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
156     else:
157     if not g.getShape()==(self.domain.getDim(),):
158     raise ValueError,"illegal shape of g"
159     self.__g=g
160 gross 1659
161 gross 2100 if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
162     if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
163 gross 1659
164 gross 2100 if permeability!=None:
165     perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
166     if perm.getRank()==0:
167     perm=perm*util.kronecker(self.domain.getDim())
168     elif perm.getRank()==1:
169     perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
170     for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
171     elif perm.getRank()==2:
172     pass
173     else:
174     raise ValueError,"illegal rank of permeability."
175     self.__permeability=perm
176     self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
177 gross 1659
178 gross 2264 def setTolerance(self,rtol=1e-4):
179     """
180     sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
181 gross 1659
182 gross 2264 M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
183    
184     where C{atol} is an absolut tolerance (see L{setAbsoluteTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
185    
186     @param rtol: relative tolerance for the pressure
187     @type rtol: non-negative C{float}
188 gross 2100 """
189 gross 2264 if rtol<0:
190     raise ValueError,"Relative tolerance needs to be non-negative."
191     self.__rtol=rtol
192     def getTolerance(self):
193 gross 2100 """
194 gross 2264 returns the relative tolerance
195 gross 1659
196 gross 2264 @return: current relative tolerance
197     @rtype: C{float}
198 caltinay 2169 """
199 gross 2264 return self.__rtol
200 gross 1659
201 gross 2264 def setAbsoluteTolerance(self,atol=0.):
202 gross 2208 """
203 gross 2264 sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
204 gross 1659
205 gross 2264 M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
206 gross 2208
207 gross 2264 where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
208 gross 2208
209 gross 2264 @param atol: absolute tolerance for the pressure
210     @type atol: non-negative C{float}
211     """
212     if atol<0:
213     raise ValueError,"Absolute tolerance needs to be non-negative."
214     self.__atol=atol
215     def getAbsoluteTolerance(self):
216     """
217     returns the absolute tolerance
218    
219     @return: current absolute tolerance
220     @rtype: C{float}
221     """
222     return self.__atol
223     def getSubProblemTolerance(self):
224 gross 2474 """
225     Returns a suitable subtolerance
226     @type: C{float}
227     """
228     return max(util.EPSILON**(0.75),self.getTolerance()**2)
229     def setSubProblemTolerance(self):
230 gross 2264 """
231 gross 2474 Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.
232 gross 2264 """
233 gross 2474 if self.__adaptSubTolerance:
234     sub_tol=self.getSubProblemTolerance()
235     self.getSolverOptionsFlux().setTolerance(sub_tol)
236     self.getSolverOptionsFlux().setAbsoluteTolerance(0.)
237     self.getSolverOptionsPressure().setTolerance(sub_tol)
238     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
239     if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol
240 gross 2208
241 gross 2474 def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):
242 gross 2264 """
243     solves the problem.
244 gross 1659
245 gross 2208 The iteration is terminated if the residual norm is less then self.getTolerance().
246 gross 1659
247 gross 2208 @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.
248     @type u0: vector value on the domain (e.g. L{Data}).
249     @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.
250     @type p0: scalar value on the domain (e.g. L{Data}).
251     @param verbose: if set some information on iteration progress are printed
252     @type verbose: C{bool}
253     @return: flux and pressure
254     @rtype: C{tuple} of L{Data}.
255 gross 2264
256 gross 2208 @note: The problem is solved as a least squares form
257 gross 2100
258 gross 2208 M{(I+D^*D)u+Qp=D^*f+g}
259     M{Q^*u+Q^*Qp=Q^*g}
260 gross 2100
261 gross 2264 where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
262     We eliminate the flux form the problem by setting
263 caltinay 2169
264 gross 2208 M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
265 caltinay 2169
266 gross 2208 form the first equation. Inserted into the second equation we get
267 caltinay 2169
268 gross 2208 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
269 gross 2264
270     which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
271 gross 2208 PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
272     """
273 gross 2386 self.verbose=verbose
274 gross 2264 rtol=self.getTolerance()
275     atol=self.getAbsoluteTolerance()
276 gross 2474 self.setSubProblemTolerance()
277    
278 gross 2264 num_corrections=0
279     converged=False
280     p=p0
281     norm_r=None
282     while not converged:
283 gross 2474 v=self.getFlux(p, fixed_flux=u0)
284 gross 2264 Qp=self.__Q(p)
285     norm_v=self.__L2(v)
286     norm_Qp=self.__L2(Qp)
287     if norm_v == 0.:
288     if norm_Qp == 0.:
289     return v,p
290     else:
291     fac=norm_Qp
292     else:
293     if norm_Qp == 0.:
294     fac=norm_v
295     else:
296     fac=2./(1./norm_v+1./norm_Qp)
297     ATOL=(atol+rtol*fac)
298     if self.verbose:
299     print "DarcyFlux: L2 norm of v = %e."%norm_v
300     print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
301 gross 2486 print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)
302     print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
303 gross 2264 print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
304     if norm_r == None or norm_r>ATOL:
305     if num_corrections>max_num_corrections:
306     raise ValueError,"maximum number of correction steps reached."
307 gross 2354 p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
308 gross 2264 num_corrections+=1
309     else:
310     converged=True
311     return v,p
312     def __L2(self,v):
313     return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
314    
315     def __Q(self,p):
316     return util.tensor_mult(self.__permeability,util.grad(p))
317    
318     def __Aprod(self,dp):
319 gross 2474 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
320 gross 2264 Qdp=self.__Q(dp)
321     self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
322 gross 2474 du=self.__pde_v.getSolution()
323 gross 2370 # self.__pde_v.getOperator().saveMM("proj.mm")
324 gross 2264 return Qdp+du
325     def __inner_GMRES(self,r,s):
326     return util.integrate(util.inner(r,s))
327    
328 gross 2100 def __inner_PCG(self,p,r):
329 gross 2264 return util.integrate(util.inner(self.__Q(p), r))
330 gross 2100
331     def __Msolve_PCG(self,r):
332 gross 2474 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
333 gross 2264 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
334 gross 2370 # self.__pde_p.getOperator().saveMM("prec.mm")
335 gross 2474 return self.__pde_p.getSolution()
336 gross 2100
337 gross 2474 def getFlux(self,p=None, fixed_flux=Data()):
338 gross 2264 """
339     returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
340     on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
341     Note that C{g} and C{f} are used, see L{setValue}.
342    
343     @param p: pressure.
344     @type p: scalar value on the domain (e.g. L{Data}).
345     @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
346     @type fixed_flux: vector values on the domain (e.g. L{Data}).
347     @param tol: relative tolerance to be used.
348     @type tol: positive C{float}.
349     @return: flux
350     @rtype: L{Data}
351     @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}}
352     for the permeability M{k_{ij}}
353     """
354 gross 2474 self.setSubProblemTolerance()
355 gross 2264 g=self.__g
356     f=self.__f
357 gross 2267 self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
358 gross 2264 if p == None:
359     self.__pde_v.setValue(Y=g)
360     else:
361     self.__pde_v.setValue(Y=g-self.__Q(p))
362 gross 2474 return self.__pde_v.getSolution()
363 gross 2264
364 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
365 gross 2251 """
366 gross 2264 solves
367 gross 1414
368 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
369     u_{i,i}=0
370 gross 1414
371 gross 2208 u=0 where fixed_u_mask>0
372     eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
373 gross 1414
374 gross 2264 if surface_stress is not given 0 is assumed.
375 gross 1414
376 gross 2251 typical usage:
377 gross 1414
378 gross 2208 sp=StokesProblemCartesian(domain)
379     sp.setTolerance()
380     sp.initialize(...)
381     v,p=sp.solve(v0,p0)
382 gross 2251 """
383 gross 2474 def __init__(self,domain,adaptSubTolerance=True, **kwargs):
384 gross 2100 """
385 gross 2208 initialize the Stokes Problem
386 gross 2100
387 gross 2208 @param domain: domain of the problem. The approximation order needs to be two.
388 gross 2100 @type domain: L{Domain}
389 gross 2474 @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
390     @type adaptSubTolerance: C{bool}
391 gross 2208 @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
392 gross 2100 """
393 gross 2474 HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**kwargs)
394 gross 1414 self.domain=domain
395     self.vol=util.integrate(1.,Function(self.domain))
396     self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
397     self.__pde_u.setSymmetryOn()
398 gross 2474
399 gross 1414 self.__pde_prec=LinearPDE(domain)
400     self.__pde_prec.setReducedOrderOn()
401     self.__pde_prec.setSymmetryOn()
402    
403 gross 2415 self.__pde_proj=LinearPDE(domain)
404     self.__pde_proj.setReducedOrderOn()
405     self.__pde_proj.setValue(D=1)
406     self.__pde_proj.setSymmetryOn()
407    
408 gross 2474 def getSolverOptionsVelocity(self):
409     """
410     returns the solver options used solve the equation for velocity.
411    
412     @rtype: L{SolverOptions}
413     """
414     return self.__pde_u.getSolverOptions()
415     def setSolverOptionsVelocity(self, options=None):
416     """
417     set the solver options for solving the equation for velocity.
418    
419     @param options: new solver options
420     @type options: L{SolverOptions}
421     """
422     self.__pde_u.setSolverOptions(options)
423     def getSolverOptionsPressure(self):
424     """
425     returns the solver options used solve the equation for pressure.
426     @rtype: L{SolverOptions}
427     """
428     return self.__pde_prec.getSolverOptions()
429     def setSolverOptionsPressure(self, options=None):
430     """
431     set the solver options for solving the equation for pressure.
432     @param options: new solver options
433     @type options: L{SolverOptions}
434     """
435     self.__pde_prec.setSolverOptions(options)
436 gross 2415
437 gross 2474 def setSolverOptionsDiv(self, options=None):
438     """
439     set the solver options for solving the equation to project the divergence of
440     the velocity onto the function space of presure.
441    
442     @param options: new solver options
443     @type options: L{SolverOptions}
444     """
445     self.__pde_prec.setSolverOptions(options)
446     def getSolverOptionsDiv(self):
447     """
448     returns the solver options for solving the equation to project the divergence of
449     the velocity onto the function space of presure.
450    
451     @rtype: L{SolverOptions}
452     """
453     return self.__pde_prec.getSolverOptions()
454     def setSubProblemTolerance(self):
455     """
456     Updates the tolerance for subproblems
457     """
458     if self.adaptSubTolerance():
459     sub_tol=self.getSubProblemTolerance()
460     self.getSolverOptionsDiv().setTolerance(sub_tol)
461     self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
462     self.getSolverOptionsPressure().setTolerance(sub_tol)
463     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
464     self.getSolverOptionsVelocity().setTolerance(sub_tol)
465     self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
466    
467    
468 gross 2620 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
469 gross 2208 """
470     assigns values to the model parameters
471 gross 2100
472 gross 2208 @param f: external force
473     @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
474     @param fixed_u_mask: mask of locations with fixed velocity.
475     @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
476     @param eta: viscosity
477     @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
478     @param surface_stress: normal surface stress
479     @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
480     @param stress: initial stress
481     @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
482     @note: All values needs to be set.
483     """
484     self.eta=eta
485     A =self.__pde_u.createCoefficient("A")
486     self.__pde_u.setValue(A=Data())
487     for i in range(self.domain.getDim()):
488     for j in range(self.domain.getDim()):
489 gross 2264 A[i,j,j,i] += 1.
490 gross 2208 A[i,j,i,j] += 1.
491 gross 2620 n=self.domain.getNormal()
492 gross 2264 self.__pde_prec.setValue(D=1/self.eta)
493 gross 2620 self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask, d=restoration_factor*util.outer(n,n))
494 gross 2251 self.__f=f
495     self.__surface_stress=surface_stress
496 gross 2208 self.__stress=stress
497 gross 1414
498 gross 2445 def Bv(self,v):
499 gross 2251 """
500     returns inner product of element p and div(v)
501 gross 1414
502 gross 2251 @param p: a pressure increment
503     @param v: a residual
504     @return: inner product of element p and div(v)
505     @rtype: C{float}
506 gross 2100 """
507 gross 2445 self.__pde_proj.setValue(Y=-util.div(v))
508     return self.__pde_proj.getSolution()
509 gross 2208
510 gross 2445 def inner_pBv(self,p,Bv):
511     """
512     returns inner product of element p and Bv=-div(v)
513    
514     @param p: a pressure increment
515     @param v: a residual
516     @return: inner product of element p and Bv=-div(v)
517     @rtype: C{float}
518     """
519     return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
520    
521 gross 2251 def inner_p(self,p0,p1):
522 gross 2100 """
523 gross 2251 Returns inner product of p0 and p1
524 gross 1414
525 gross 2251 @param p0: a pressure
526     @param p1: a pressure
527     @return: inner product of p0 and p1
528 gross 2208 @rtype: C{float}
529 gross 2100 """
530     s0=util.interpolate(p0/self.eta,Function(self.domain))
531     s1=util.interpolate(p1/self.eta,Function(self.domain))
532     return util.integrate(s0*s1)
533 artak 1550
534 gross 2251 def norm_v(self,v):
535 gross 2100 """
536 gross 2251 returns the norm of v
537 gross 2208
538 gross 2251 @param v: a velovity
539     @return: norm of v
540     @rtype: non-negative C{float}
541 gross 2100 """
542 gross 2251 return util.sqrt(util.integrate(util.length(util.grad(v))))
543 gross 2100
544 gross 2251 def getV(self, p, v0):
545 gross 1414 """
546 gross 2251 return the value for v for a given p (overwrite)
547    
548     @param p: a pressure
549 gross 2264 @param v0: a initial guess for the value v to return.
550 gross 2251 @return: v given as M{v= A^{-1} (f-B^*p)}
551 gross 1414 """
552 gross 2251 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
553 gross 2100 if self.__stress.isEmpty():
554 gross 2251 self.__pde_u.setValue(X=p*util.kronecker(self.domain))
555 gross 2100 else:
556 gross 2251 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
557 gross 2474 out=self.__pde_u.getSolution()
558 gross 2208 return out
559 gross 1414
560 gross 2445 def norm_Bv(self,Bv):
561 gross 2251 """
562     Returns Bv (overwrite).
563    
564     @rtype: equal to the type of p
565     @note: boundary conditions on p should be zero!
566     """
567 gross 2445 return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
568 gross 2251
569     def solve_AinvBt(self,p):
570     """
571     Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
572    
573     @param p: a pressure increment
574 gross 2264 @return: the solution of M{Av=B^*p}
575 gross 2251 @note: boundary conditions on v should be zero!
576     """
577     self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
578 gross 2474 out=self.__pde_u.getSolution()
579 gross 2251 return out
580    
581 gross 2445 def solve_prec(self,Bv):
582 gross 2251 """
583     applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
584 gross 2445 with accuracy L{self.getSubProblemTolerance()}
585 gross 2251
586     @param v: velocity increment
587     @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
588     @note: boundary conditions on p are zero.
589     """
590 gross 2445 self.__pde_prec.setValue(Y=Bv)
591 gross 2474 return self.__pde_prec.getSolution()

  ViewVC Help
Powered by ViewVC 1.1.26