/[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 2793 - (hide annotations)
Tue Dec 1 06:10:10 2009 UTC (9 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 24362 byte(s)
some new util functions
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 jfenwick 2625 :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 gross 1414 """
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 jfenwick 2625 *u_i+k_{ij}*p_{,j} = g_i*
44     *u_{i,i} = f*
45 gross 1659
46 jfenwick 2625 where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
47 gross 1659
48 jfenwick 2625 :note: The problem is solved in a least squares formulation.
49 gross 2100 """
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 jfenwick 2625 :param domain: domain of the problem
55     :type domain: `Domain`
56     :param useReduced: uses reduced oreder on flux and pressure
57     :type useReduced: ``bool``
58     :param adaptSubTolerance: switches on automatic subtolerance selection
59     :type adaptSubTolerance: ``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 artak 2689 #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 jfenwick 2625 *(I+D^*D)u=F*
87 gross 2474
88 jfenwick 2625 :return: `SolverOptions`
89 gross 2474 """
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 jfenwick 2625 *(I+D^*D)u=F*
96 gross 2474
97 jfenwick 2625 If ``options`` is not present, the options are reset to default
98     :param options: `SolverOptions`
99     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
100 gross 2474 """
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 jfenwick 2625 *(Q^*Q)p=Q^*G*
107 gross 2474
108 jfenwick 2625 :return: `SolverOptions`
109 gross 2474 """
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 jfenwick 2625 *(Q^*Q)p=Q^*G*
116 gross 2474
117 jfenwick 2625 If ``options`` is not present, the options are reset to default
118     :param options: `SolverOptions`
119     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
120 gross 2474 """
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 jfenwick 2625 :param f: volumetic sources/sinks
128     :type f: scalar value on the domain (e.g. `Data`)
129     :param g: flux sources/sinks
130     :type g: vector values on the domain (e.g. `Data`)
131     :param location_of_fixed_pressure: mask for locations where pressure is fixed
132     :type location_of_fixed_pressure: scalar value on the domain (e.g. `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. `Data`)
135     :param permeability: permeability tensor. If scalar ``s`` is given the tensor with
136     ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with
137     ``v`` on the main diagonal is used.
138     :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)
139 gross 1659
140 jfenwick 2625 :note: the values of parameters which are not set by calling ``setValue`` are not altered.
141     :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)
142     or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal
143     is along the *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 jfenwick 2625 sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if
181 gross 1659
182 jfenwick 2625 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
183 gross 2264
184 jfenwick 2625 where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
185 gross 2264
186 jfenwick 2625 :param rtol: relative tolerance for the pressure
187     :type rtol: non-negative ``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 jfenwick 2625 :return: current relative tolerance
197     :rtype: ``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 jfenwick 2625 sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if
204 gross 1659
205 jfenwick 2625 *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*
206 gross 2208
207 jfenwick 2625 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}*.
208 gross 2208
209 jfenwick 2625 :param atol: absolute tolerance for the pressure
210     :type atol: non-negative ``float``
211 gross 2264 """
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 jfenwick 2625 :return: current absolute tolerance
220     :rtype: ``float``
221 gross 2264 """
222     return self.__atol
223     def getSubProblemTolerance(self):
224 gross 2474 """
225     Returns a suitable subtolerance
226 jfenwick 2625 @type: ``float``
227 gross 2474 """
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 jfenwick 2625 :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.
248     :type u0: vector value on the domain (e.g. `Data`).
249     :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.
250     :type p0: scalar value on the domain (e.g. `Data`).
251     :param verbose: if set some information on iteration progress are printed
252     :type verbose: ``bool``
253     :return: flux and pressure
254     :rtype: ``tuple`` of `Data`.
255 gross 2264
256 jfenwick 2625 :note: The problem is solved as a least squares form
257 gross 2100
258 jfenwick 2625 *(I+D^*D)u+Qp=D^*f+g*
259     *Q^*u+Q^*Qp=Q^*g*
260 gross 2100
261 jfenwick 2625 where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.
262 gross 2264 We eliminate the flux form the problem by setting
263 caltinay 2169
264 jfenwick 2625 *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 jfenwick 2625 *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 jfenwick 2625 which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step
271     PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.
272 gross 2208 """
273 gross 2386 self.verbose=verbose
274 gross 2264 rtol=self.getTolerance()
275     atol=self.getAbsoluteTolerance()
276 gross 2474 self.setSubProblemTolerance()
277 gross 2264 num_corrections=0
278     converged=False
279     p=p0
280     norm_r=None
281     while not converged:
282 gross 2474 v=self.getFlux(p, fixed_flux=u0)
283 gross 2264 Qp=self.__Q(p)
284     norm_v=self.__L2(v)
285     norm_Qp=self.__L2(Qp)
286     if norm_v == 0.:
287     if norm_Qp == 0.:
288     return v,p
289     else:
290     fac=norm_Qp
291     else:
292     if norm_Qp == 0.:
293     fac=norm_v
294     else:
295     fac=2./(1./norm_v+1./norm_Qp)
296     ATOL=(atol+rtol*fac)
297     if self.verbose:
298     print "DarcyFlux: L2 norm of v = %e."%norm_v
299     print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
300 gross 2486 print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)
301     print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)
302 gross 2264 print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
303     if norm_r == None or norm_r>ATOL:
304     if num_corrections>max_num_corrections:
305     raise ValueError,"maximum number of correction steps reached."
306 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)
307 gross 2264 num_corrections+=1
308     else:
309     converged=True
310     return v,p
311     def __L2(self,v):
312     return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
313    
314     def __Q(self,p):
315     return util.tensor_mult(self.__permeability,util.grad(p))
316    
317     def __Aprod(self,dp):
318 gross 2474 if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"
319 gross 2264 Qdp=self.__Q(dp)
320     self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
321 gross 2474 du=self.__pde_v.getSolution()
322 gross 2370 # self.__pde_v.getOperator().saveMM("proj.mm")
323 gross 2264 return Qdp+du
324     def __inner_GMRES(self,r,s):
325     return util.integrate(util.inner(r,s))
326    
327 gross 2100 def __inner_PCG(self,p,r):
328 gross 2264 return util.integrate(util.inner(self.__Q(p), r))
329 gross 2100
330     def __Msolve_PCG(self,r):
331 gross 2474 if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"
332 gross 2264 self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
333 gross 2370 # self.__pde_p.getOperator().saveMM("prec.mm")
334 gross 2474 return self.__pde_p.getSolution()
335 gross 2100
336 gross 2474 def getFlux(self,p=None, fixed_flux=Data()):
337 gross 2264 """
338 jfenwick 2625 returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``
339     on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
340     Note that ``g`` and ``f`` are used, see `setValue`.
341 gross 2264
342 jfenwick 2625 :param p: pressure.
343     :type p: scalar value on the domain (e.g. `Data`).
344     :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.
345     :type fixed_flux: vector values on the domain (e.g. `Data`).
346     :return: flux
347     :rtype: `Data`
348     :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*
349     for the permeability *k_{ij}*
350 gross 2264 """
351 gross 2474 self.setSubProblemTolerance()
352 gross 2264 g=self.__g
353     f=self.__f
354 gross 2267 self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
355 gross 2264 if p == None:
356     self.__pde_v.setValue(Y=g)
357     else:
358     self.__pde_v.setValue(Y=g-self.__Q(p))
359 gross 2474 return self.__pde_v.getSolution()
360 gross 2264
361 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
362 gross 2251 """
363 gross 2264 solves
364 gross 1414
365 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
366     u_{i,i}=0
367 gross 1414
368 gross 2208 u=0 where fixed_u_mask>0
369     eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
370 gross 1414
371 gross 2264 if surface_stress is not given 0 is assumed.
372 gross 1414
373 gross 2251 typical usage:
374 gross 1414
375 gross 2208 sp=StokesProblemCartesian(domain)
376     sp.setTolerance()
377     sp.initialize(...)
378     v,p=sp.solve(v0,p0)
379 gross 2251 """
380 gross 2719 def __init__(self,domain,**kwargs):
381 gross 2100 """
382 gross 2208 initialize the Stokes Problem
383 gross 2100
384 gross 2793 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
385     LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
386     with macro elements for the pressure.
387    
388     :param domain: domain of the problem.
389 jfenwick 2625 :type domain: `Domain`
390 gross 2100 """
391 gross 2719 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
392 gross 1414 self.domain=domain
393     self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
394     self.__pde_u.setSymmetryOn()
395 gross 2474
396 gross 1414 self.__pde_prec=LinearPDE(domain)
397     self.__pde_prec.setReducedOrderOn()
398     self.__pde_prec.setSymmetryOn()
399    
400 gross 2415 self.__pde_proj=LinearPDE(domain)
401     self.__pde_proj.setReducedOrderOn()
402     self.__pde_proj.setValue(D=1)
403     self.__pde_proj.setSymmetryOn()
404    
405 gross 2474 def getSolverOptionsVelocity(self):
406     """
407     returns the solver options used solve the equation for velocity.
408    
409 jfenwick 2625 :rtype: `SolverOptions`
410 gross 2474 """
411     return self.__pde_u.getSolverOptions()
412     def setSolverOptionsVelocity(self, options=None):
413     """
414     set the solver options for solving the equation for velocity.
415    
416 jfenwick 2625 :param options: new solver options
417     :type options: `SolverOptions`
418 gross 2474 """
419     self.__pde_u.setSolverOptions(options)
420     def getSolverOptionsPressure(self):
421     """
422     returns the solver options used solve the equation for pressure.
423 jfenwick 2625 :rtype: `SolverOptions`
424 gross 2474 """
425     return self.__pde_prec.getSolverOptions()
426     def setSolverOptionsPressure(self, options=None):
427     """
428     set the solver options for solving the equation for pressure.
429 jfenwick 2625 :param options: new solver options
430     :type options: `SolverOptions`
431 gross 2474 """
432     self.__pde_prec.setSolverOptions(options)
433 gross 2415
434 gross 2474 def setSolverOptionsDiv(self, options=None):
435     """
436     set the solver options for solving the equation to project the divergence of
437     the velocity onto the function space of presure.
438    
439 jfenwick 2625 :param options: new solver options
440     :type options: `SolverOptions`
441 gross 2474 """
442 artak 2689 self.__pde_proj.setSolverOptions(options)
443 gross 2474 def getSolverOptionsDiv(self):
444     """
445     returns the solver options for solving the equation to project the divergence of
446     the velocity onto the function space of presure.
447    
448 jfenwick 2625 :rtype: `SolverOptions`
449 gross 2474 """
450 artak 2689 return self.__pde_proj.getSolverOptions()
451 gross 2474
452 gross 2793 def updateStokesEquation(self, v, p):
453     """
454     updates the Stokes equation to consider dependencies from ``v`` and ``p``
455     :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
456     """
457     pass
458     def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
459     """
460     assigns new values to the model parameters.
461    
462     :param f: external force
463     :type f: `Vector` object in `FunctionSpace` `Function` or similar
464     :param fixed_u_mask: mask of locations with fixed velocity.
465     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
466     :param eta: viscosity
467     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
468     :param surface_stress: normal surface stress
469     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
470     :param stress: initial stress
471     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
472     """
473     if eta !=None:
474     k=util.kronecker(self.domain.getDim())
475     kk=util.outer(k,k)
476     self.eta=util.interpolate(eta, Function(self.domain))
477     self.__pde_prec.setValue(D=1/self.eta)
478     self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
479     if restoration_factor!=None:
480     n=self.domain.getNormal()
481     self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))
482     if fixed_u_mask!=None:
483     self.__pde_u.setValue(q=fixed_u_mask)
484     if f!=None: self.__f=f
485     if surface_stress!=None: self.__surface_stress=surface_stress
486     if stress!=None: self.__stress=stress
487    
488 gross 2620 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):
489 gross 2208 """
490     assigns values to the model parameters
491 gross 2100
492 jfenwick 2625 :param f: external force
493     :type f: `Vector` object in `FunctionSpace` `Function` or similar
494     :param fixed_u_mask: mask of locations with fixed velocity.
495     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
496     :param eta: viscosity
497     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
498     :param surface_stress: normal surface stress
499     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
500     :param stress: initial stress
501     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
502 gross 2208 """
503 gross 2793 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
504 gross 1414
505 gross 2719 def Bv(self,v,tol):
506 gross 2251 """
507     returns inner product of element p and div(v)
508 gross 1414
509 jfenwick 2625 :param v: a residual
510     :return: inner product of element p and div(v)
511     :rtype: ``float``
512 gross 2100 """
513 gross 2793 self.__pde_proj.setValue(Y=-util.div(v))
514 gross 2719 self.getSolverOptionsDiv().setTolerance(tol)
515     self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
516     out=self.__pde_proj.getSolution()
517     return out
518 gross 2208
519 gross 2445 def inner_pBv(self,p,Bv):
520     """
521     returns inner product of element p and Bv=-div(v)
522    
523 jfenwick 2625 :param p: a pressure increment
524     :param Bv: a residual
525     :return: inner product of element p and Bv=-div(v)
526     :rtype: ``float``
527 gross 2445 """
528     return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
529    
530 gross 2251 def inner_p(self,p0,p1):
531 gross 2100 """
532 gross 2251 Returns inner product of p0 and p1
533 gross 1414
534 jfenwick 2625 :param p0: a pressure
535     :param p1: a pressure
536     :return: inner product of p0 and p1
537     :rtype: ``float``
538 gross 2100 """
539 gross 2719 s0=util.interpolate(p0,Function(self.domain))
540     s1=util.interpolate(p1,Function(self.domain))
541 gross 2100 return util.integrate(s0*s1)
542 artak 1550
543 gross 2251 def norm_v(self,v):
544 gross 2100 """
545 gross 2251 returns the norm of v
546 gross 2208
547 jfenwick 2625 :param v: a velovity
548     :return: norm of v
549     :rtype: non-negative ``float``
550 gross 2100 """
551 gross 2719 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
552 gross 2100
553 gross 2793
554 gross 2719 def getDV(self, p, v, tol):
555 gross 1414 """
556 gross 2251 return the value for v for a given p (overwrite)
557    
558 jfenwick 2625 :param p: a pressure
559 gross 2719 :param v: a initial guess for the value v to return.
560     :return: dv given as *Adv=(f-Av-B^*p)*
561 gross 1414 """
562 gross 2793 self.updateStokesEquation(v,p)
563 gross 2719 self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)
564     self.getSolverOptionsVelocity().setTolerance(tol)
565     self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
566 gross 2100 if self.__stress.isEmpty():
567 gross 2719 self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
568 gross 2100 else:
569 gross 2719 self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
570 gross 2474 out=self.__pde_u.getSolution()
571 gross 2208 return out
572 gross 1414
573 gross 2445 def norm_Bv(self,Bv):
574 gross 2251 """
575     Returns Bv (overwrite).
576    
577 jfenwick 2625 :rtype: equal to the type of p
578     :note: boundary conditions on p should be zero!
579 gross 2251 """
580 gross 2445 return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
581 gross 2251
582 gross 2719 def solve_AinvBt(self,p, tol):
583 gross 2251 """
584 gross 2719 Solves *Av=B^*p* with accuracy `tol`
585 gross 2251
586 jfenwick 2625 :param p: a pressure increment
587     :return: the solution of *Av=B^*p*
588     :note: boundary conditions on v should be zero!
589 gross 2251 """
590 gross 2719 self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))
591 gross 2474 out=self.__pde_u.getSolution()
592 gross 2251 return out
593    
594 gross 2719 def solve_prec(self,Bv, tol):
595 gross 2251 """
596 jfenwick 2625 applies preconditioner for for *BA^{-1}B^** to *Bv*
597     with accuracy `self.getSubProblemTolerance()`
598 gross 2251
599 jfenwick 2625 :param Bv: velocity increment
600     :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
601     :note: boundary conditions on p are zero.
602 gross 2251 """
603 gross 2445 self.__pde_prec.setValue(Y=Bv)
604 gross 2719 self.getSolverOptionsPressure().setTolerance(tol)
605     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
606     out=self.__pde_prec.getSolution()
607     return out

  ViewVC Help
Powered by ViewVC 1.1.26