/[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 3515 - (hide annotations)
Thu May 19 08:20:57 2011 UTC (8 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 26387 byte(s)
some first work for DiracFunctions
1 gross 3071 # -*- coding: utf-8 -*-
2 ksteube 1809 ########################################################
3 gross 1414 #
4 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 gross 1414 #
8 ksteube 1809 # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11 gross 1414 #
12 ksteube 1809 ########################################################
13 gross 1414
14 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 ksteube 1809 Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18     __license__="""Licensed under the Open Software License version 3.0
19     http://www.opensource.org/licenses/osl-3.0.php"""
20 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1809
22 gross 1414 """
23     Some models for flow
24    
25 jfenwick 2625 :var __author__: name of author
26     :var __copyright__: copyrights
27     :var __license__: licence agreement
28     :var __url__: url entry point on documentation
29     :var __version__: version
30     :var __date__: date of the version
31 gross 1414 """
32    
33     __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35 jfenwick 3432 import escript
36 gross 1414 import util
37 gross 2474 from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38 gross 2264 from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39 gross 1414
40 gross 3502 class DarcyFlow(object):
41 gross 3499 """
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 gross 3502 :cvar SIMPLE: simple solver
50     :cvar POST: solver using global postprocessing of flux
51     :cvar STAB: solver uses (non-symmetric) stabilization
52     :cvar SYMSTAB: solver uses symmetric stabilization
53 gross 3499 """
54 gross 3502 SIMPLE="SIMPLE"
55     POST="POST"
56     STAB="STAB"
57     SYMSTAB="SYMSTAB"
58     def __init__(self, domain, useReduced=False, solver="SYMSTAB", verbose=False, w=1.):
59 gross 3499 """
60     initializes the Darcy flux problem
61     :param domain: domain of the problem
62     :type domain: `Domain`
63     :param useReduced: uses reduced oreder on flux and pressure
64     :type useReduced: ``bool``
65 gross 3502 :param solver: solver method
66     :type solver: in [`DarcyFlow.SIMPLE`, `DarcyFlow.POST', `DarcyFlow.STAB`, `DarcyFlow.SYMSTAB` ]
67     :param verbose: if ``True`` some information on the iteration progress are printed.
68     :type verbose: ``bool``
69     :param w: weighting factor for `DarcyFlow.POST` solver
70     :type w: ``float``
71    
72 gross 3499 """
73     self.domain=domain
74 gross 3502 self.solver=solver
75 gross 3499 self.useReduced=useReduced
76 gross 3502 self.verbose=verbose
77 gross 3499 self.scale=1.
78 gross 3500
79 gross 3499
80 gross 3502 self.__pde_v=LinearPDESystem(domain)
81     self.__pde_v.setSymmetryOn()
82     if self.useReduced: self.__pde_v.setReducedOrderOn()
83     self.__pde_p=LinearSinglePDE(domain)
84     self.__pde_p.setSymmetryOn()
85     if self.useReduced: self.__pde_p.setReducedOrderOn()
86 gross 3499
87 gross 3502 if self.solver == self.SIMPLE:
88     if self.verbose: print "DarcyFlow: simple solver is used."
89     self.__pde_v.setValue(D=util.kronecker(self.domain.getDim()))
90     elif self.solver == self.POST:
91     self.w=w
92     if util.inf(w)<0.:
93     raise ValueError,"Weighting factor must be non-negative."
94     if self.verbose: print "DarcyFlow: global postprocessing of flux is used."
95     elif self.solver == self.STAB:
96     if self.verbose: print "DarcyFlow: (non-symmetric) stabilization is used."
97     elif self.solver == self.SYMSTAB:
98     if self.verbose: print "DarcyFlow: symmetric stabilization is used."
99     else:
100     raise ValueError,"unknown solver %s"%self.solver
101     self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
102     self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
103     self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
104     self.location_of_fixed_flux = escript.Vector(0, self.__pde_v.getFunctionSpaceForCoefficient("q"))
105     self.setTolerance()
106    
107    
108 gross 3499 def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
109     """
110     assigns values to model parameters
111    
112     :param f: volumetic sources/sinks
113     :type f: scalar value on the domain (e.g. `escript.Data`)
114     :param g: flux sources/sinks
115     :type g: vector values on the domain (e.g. `escript.Data`)
116     :param location_of_fixed_pressure: mask for locations where pressure is fixed
117     :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
118     :param location_of_fixed_flux: mask for locations where flux is fixed.
119     :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
120     :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
121 gross 3502 :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
122 gross 3499
123     :note: the values of parameters which are not set by calling ``setValue`` are not altered.
124     :note: at any point on the boundary of the domain the pressure
125     (``location_of_fixed_pressure`` >0) or the normal component of the
126     flux (``location_of_fixed_flux[i]>0``) if direction of the normal
127     is along the *x_i* axis.
128    
129     """
130 gross 3502 if location_of_fixed_pressure!=None:
131     self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
132     self.__pde_p.setValue(q=self.location_of_fixed_pressure)
133     if location_of_fixed_flux!=None:
134     self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
135     self.__pde_v.setValue(q=self.location_of_fixed_flux)
136 gross 3499
137    
138 gross 3502 # pressure is rescaled by the factor 1/self.scale
139 gross 3499 if permeability!=None:
140 gross 3502
141     perm=util.interpolate(permeability,self.__pde_v.getFunctionSpaceForCoefficient("A"))
142 gross 3499 V=util.vol(self.domain)
143 gross 3502 l=V**(1./self.domain.getDim())
144    
145 gross 3499 if perm.getRank()==0:
146     perm_inv=(1./perm)
147 gross 3502 self.scale=util.integrate(perm_inv)/V*l
148 gross 3499 perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
149     perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
150 gross 3502
151    
152 gross 3499 elif perm.getRank()==2:
153     perm_inv=util.inverse(perm)
154 gross 3502 self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)/V)*l
155 gross 3499 perm_inv*=(1./self.scale)
156     perm=perm*self.scale
157     else:
158     raise ValueError,"illegal rank of permeability."
159    
160     self.__permeability=perm
161     self.__permeability_inv=perm_inv
162 gross 3502 if self.verbose: print "DarcyFlow: scaling factor for pressure is %e."%self.scale
163    
164     if self.solver == self.SIMPLE:
165     self.__pde_p.setValue(A=self.__permeability)
166     elif self.solver == self.POST:
167     self.__pde_p.setValue(A=self.__permeability)
168     k=util.kronecker(self.domain.getDim())
169     self.lamb = self.w*util.length(perm_inv)*l
170     self.__pde_v.setValue(D=self.__permeability_inv, A=self.lamb*self.domain.getSize()*util.outer(k,k))
171     elif self.solver == self.STAB:
172 gross 3499 self.__pde_p.setValue(A=0.5*self.__permeability)
173 gross 3502 self.__pde_v.setValue(D=0.5*self.__permeability_inv)
174     elif self.solver == self.SYMSTAB:
175     self.__pde_p.setValue(A=0.5*self.__permeability)
176     self.__pde_v.setValue(D=0.5*self.__permeability_inv)
177    
178 gross 3499 if g != None:
179 gross 3502 g=util.interpolate(g, self.__pde_v.getFunctionSpaceForCoefficient("Y"))
180 gross 3499 if g.isEmpty():
181 gross 3502 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
182 gross 3499 else:
183     if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
184 gross 3502 self.__g=g
185 gross 3499 if f !=None:
186 gross 3502 f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187     if f.isEmpty():
188     f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189 gross 3499 else:
190     if f.getRank()>0: raise ValueError,"illegal rank of f."
191 gross 3502 self.__f=f
192 gross 3071 def getSolverOptionsFlux(self):
193     """
194     Returns the solver options used to solve the flux problems
195     :return: `SolverOptions`
196     """
197 gross 3502 return self.__pde_v.getSolverOptions()
198 gross 3071
199     def setSolverOptionsFlux(self, options=None):
200     """
201     Sets the solver options used to solve the flux problems
202     If ``options`` is not present, the options are reset to default
203     :param options: `SolverOptions`
204     """
205     return self.__pde_v.setSolverOptions(options)
206    
207     def getSolverOptionsPressure(self):
208     """
209     Returns the solver options used to solve the pressure problems
210     :return: `SolverOptions`
211     """
212     return self.__pde_p.getSolverOptions()
213    
214     def setSolverOptionsPressure(self, options=None):
215     """
216     Sets the solver options used to solve the pressure problems
217     If ``options`` is not present, the options are reset to default
218    
219     :param options: `SolverOptions`
220     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
221     """
222     return self.__pde_p.setSolverOptions(options)
223    
224     def setTolerance(self,rtol=1e-4):
225     """
226 gross 3502 sets the relative tolerance ``rtol`` for the pressure for the stabelized solvers.
227 gross 3071
228     :param rtol: relative tolerance for the pressure
229     :type rtol: non-negative ``float``
230     """
231     if rtol<0:
232     raise ValueError,"Relative tolerance needs to be non-negative."
233     self.__rtol=rtol
234 gross 3502
235 gross 3071 def getTolerance(self):
236     """
237     returns the relative tolerance
238     :return: current relative tolerance
239     :rtype: ``float``
240     """
241     return self.__rtol
242 gross 3502
243     def solve(self,u0,p0, max_iter=100, iter_restart=20):
244 gross 3071 """
245 gross 3502 solves the problem.
246 gross 3071
247 gross 3502 The iteration is terminated if the residual norm is less then self.getTolerance().
248 gross 2960
249 gross 3502 :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.
250     :type u0: vector value on the domain (e.g. `escript.Data`).
251     :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.
252     :type p0: scalar value on the domain (e.g. `escript.Data`).
253     :param max_iter: maximum number of (outer) iteration steps for the stabilization solvers,
254     :type max_iter: ``int``
255     :param iter_restart: number of steps after which the iteration is restarted. The larger ``iter_restart`` the larger the required memory.
256     A small value for ``iter_restart`` may require a large number of iteration steps or may even lead to a failure
257     of the iteration. ``iter_restart`` is relevant for the stabilization solvers only.
258     :type iter_restart: ``int``
259     :return: flux and pressure
260     :rtype: ``tuple`` of `escript.Data`.
261 gross 2960
262 gross 3071 """
263 gross 3502 # rescale initial guess:
264     p0=p0/self.scale
265     if self.solver == self.SIMPLE or self.solver == self.POST :
266     self.__pde_p.setValue(X=self.__g ,
267     Y=self.__f,
268 gross 3510 y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),
269 gross 3502 r=p0)
270     p=self.__pde_p.getSolution()
271     u = self.getFlux(p, u0)
272     elif self.solver == self.STAB:
273     u,p = self.__solve_STAB(u0,p0, max_iter, iter_restart)
274     elif self.solver == self.SYMSTAB:
275     u,p = self.__solve_SYMSTAB(u0,p0, max_iter, iter_restart)
276 gross 2474
277 gross 3502 if self.verbose:
278     KGp=util.tensor_mult(self.__permeability,util.grad(p))
279     def_p=self.__g-(u+KGp)
280     def_v=self.__f-util.div(u, self.__pde_v.getFunctionSpaceForCoefficient("X"))
281     print "DarcyFlux: |g-u-K*grad(p)|_2 = %e (|u|_2 = %e)."%(self.__L2(def_p),self.__L2(u))
282     print "DarcyFlux: |f-div(u)|_2 = %e (|grad(u)|_2 = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
283     #rescale result
284     p=p*self.scale
285     return u,p
286    
287     def getFlux(self,p, u0=None):
288 gross 2100 """
289 gross 3502 returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
290 jfenwick 2625 on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
291 gross 3502 Notice that ``g`` and ``f`` are used, see `setValue`.
292 gross 2264
293 jfenwick 2625 :param p: pressure.
294 gross 3440 :type p: scalar value on the domain (e.g. `escript.Data`).
295 gross 3502 :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
296     :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
297 jfenwick 2625 :return: flux
298 gross 3440 :rtype: `escript.Data`
299 gross 2264 """
300 gross 3502 if self.solver == self.SIMPLE or self.solver == self.POST :
301     KGp=util.tensor_mult(self.__permeability,util.grad(p))
302     self.__pde_v.setValue(Y=self.__g-KGp, X=escript.Data())
303     if u0 == None:
304     self.__pde_v.setValue(r=escript.Data())
305     else:
306     self.__pde_v.setValue(r=u0)
307     u= self.__pde_v.getSolution()
308     elif self.solver == self.POST:
309     self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g)-util.grad(p),
310     X=self.lamb * self.__f * util.kronecker(self.domain.getDim()))
311     if u0 == None:
312     self.__pde_v.setValue(r=escript.Data())
313     else:
314     self.__pde_v.setValue(r=u0)
315     u= self.__pde_v.getSolution()
316     elif self.solver == self.STAB:
317     gp=util.grad(p)
318     self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)+gp),
319     X= p * util.kronecker(self.domain.getDim()),
320     y= - p * self.domain.getNormal())
321     if u0 == None:
322     self.__pde_v.setValue(r=escript.Data())
323     else:
324     self.__pde_v.setValue(r=u0)
325     u= self.__pde_v.getSolution()
326     elif self.solver == self.SYMSTAB:
327     gp=util.grad(p)
328     self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)-gp),
329     X= escript.Data() ,
330     y= escript.Data() )
331     if u0 == None:
332     self.__pde_v.setValue(r=escript.Data())
333     else:
334     self.__pde_v.setValue(r=u0)
335     u= self.__pde_v.getSolution()
336     return u
337    
338    
339     def __solve_STAB(self, u0, p0, max_iter, iter_restart):
340     # p0 is used as an initial guess
341     u=self.getFlux(p0, u0)
342     self.__pde_p.setValue( Y=self.__f-util.div(u),
343     X=0.5*(self.__g - u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
344     y= escript.Data(),
345     r=escript.Data())
346 gross 3510
347 gross 3502 dp=self.__pde_p.getSolution()
348     p=GMRES(dp,
349     self.__STAB_Aprod,
350     p0,
351     self.__inner,
352     atol=self.__norm(p0+dp)*self.getTolerance() ,
353     rtol=0.,
354     iter_max=max_iter,
355     iter_restart=iter_restart,
356     verbose=self.verbose,P_R=None)
357    
358     u=self.getFlux(p, u0)
359     return u,p
360 gross 2264
361 gross 3502 def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):
362     # p0 is used as an initial guess
363     u=self.getFlux(p0, u0)
364     self.__pde_p.setValue( Y= self.__f,
365     X= 0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
366 gross 3510 y= - util.inner(self.domain.getNormal(), u),
367 gross 3502 r=escript.Data())
368     dp=self.__pde_p.getSolution()
369 gross 3510
370 gross 3502 p=GMRES(dp,
371     self.__SYMSTAB_Aprod,
372     p0,
373     self.__inner,
374     atol=self.__norm(p0+dp)*self.getTolerance() ,
375     rtol=0.,
376     iter_max=max_iter,
377     iter_restart=iter_restart,
378     verbose=self.verbose,P_R=None)
379    
380     u=self.getFlux(p, u0)
381     return u,p
382    
383     def __L2(self,v):
384     return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))
385    
386     def __norm(self,r):
387     return util.sqrt(self.__inner(r,r))
388    
389     def __inner(self,r,s):
390     return util.integrate(util.inner(r,s), escript.Function(self.domain))
391    
392     def __STAB_Aprod(self,p):
393     gp=util.grad(p)
394     self.__pde_v.setValue(Y=-0.5*gp,
395     X=-p*util.kronecker(self.__pde_v.getDomain()),
396     y= p * self.domain.getNormal(),
397     r=escript.Data())
398     u = -self.__pde_v.getSolution()
399     self.__pde_p.setValue(Y=util.div(u),
400     X=0.5*(u+util.tensor_mult(self.__permeability,gp)),
401     y=escript.Data(),
402     r=escript.Data())
403    
404     return self.__pde_p.getSolution()
405    
406     def __SYMSTAB_Aprod(self,p):
407     gp=util.grad(p)
408     self.__pde_v.setValue(Y=0.5*gp ,
409     X=escript.Data(),
410     y=escript.Data(),
411     r=escript.Data())
412     u = -self.__pde_v.getSolution()
413     self.__pde_p.setValue(Y=escript.Data(),
414     X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),
415 gross 3510 y=escript.Data(),
416 gross 3502 r=escript.Data())
417    
418     return self.__pde_p.getSolution()
419    
420    
421 gross 1414 class StokesProblemCartesian(HomogeneousSaddlePointProblem):
422 gross 2251 """
423 gross 2264 solves
424 gross 1414
425 gross 2208 -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
426     u_{i,i}=0
427 gross 1414
428 gross 2208 u=0 where fixed_u_mask>0
429     eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
430 gross 1414
431 gross 2264 if surface_stress is not given 0 is assumed.
432 gross 1414
433 gross 2251 typical usage:
434 gross 1414
435 gross 2208 sp=StokesProblemCartesian(domain)
436     sp.setTolerance()
437     sp.initialize(...)
438     v,p=sp.solve(v0,p0)
439 gross 2251 """
440 gross 2719 def __init__(self,domain,**kwargs):
441 gross 2100 """
442 gross 2208 initialize the Stokes Problem
443 gross 2100
444 gross 2793 The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
445     LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
446     with macro elements for the pressure.
447    
448     :param domain: domain of the problem.
449 jfenwick 2625 :type domain: `Domain`
450 gross 2100 """
451 gross 2719 HomogeneousSaddlePointProblem.__init__(self,**kwargs)
452 gross 1414 self.domain=domain
453 gross 3502 self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
454     self.__pde_v.setSymmetryOn()
455 gross 2474
456 gross 1414 self.__pde_prec=LinearPDE(domain)
457     self.__pde_prec.setReducedOrderOn()
458     self.__pde_prec.setSymmetryOn()
459    
460 gross 2415 self.__pde_proj=LinearPDE(domain)
461     self.__pde_proj.setReducedOrderOn()
462     self.__pde_proj.setValue(D=1)
463     self.__pde_proj.setSymmetryOn()
464    
465 gross 2474 def getSolverOptionsVelocity(self):
466     """
467     returns the solver options used solve the equation for velocity.
468    
469 jfenwick 2625 :rtype: `SolverOptions`
470 gross 2474 """
471 gross 3502 return self.__pde_v.getSolverOptions()
472 gross 2474 def setSolverOptionsVelocity(self, options=None):
473     """
474     set the solver options for solving the equation for velocity.
475    
476 jfenwick 2625 :param options: new solver options
477     :type options: `SolverOptions`
478 gross 2474 """
479 gross 3502 self.__pde_v.setSolverOptions(options)
480 gross 2474 def getSolverOptionsPressure(self):
481     """
482     returns the solver options used solve the equation for pressure.
483 jfenwick 2625 :rtype: `SolverOptions`
484 gross 2474 """
485     return self.__pde_prec.getSolverOptions()
486     def setSolverOptionsPressure(self, options=None):
487     """
488     set the solver options for solving the equation for pressure.
489 jfenwick 2625 :param options: new solver options
490     :type options: `SolverOptions`
491 gross 2474 """
492     self.__pde_prec.setSolverOptions(options)
493 gross 2415
494 gross 2474 def setSolverOptionsDiv(self, options=None):
495     """
496     set the solver options for solving the equation to project the divergence of
497     the velocity onto the function space of presure.
498    
499 jfenwick 2625 :param options: new solver options
500     :type options: `SolverOptions`
501 gross 2474 """
502 artak 2689 self.__pde_proj.setSolverOptions(options)
503 gross 2474 def getSolverOptionsDiv(self):
504     """
505     returns the solver options for solving the equation to project the divergence of
506     the velocity onto the function space of presure.
507    
508 jfenwick 2625 :rtype: `SolverOptions`
509 gross 2474 """
510 artak 2689 return self.__pde_proj.getSolverOptions()
511 gross 2474
512 gross 2793 def updateStokesEquation(self, v, p):
513     """
514     updates the Stokes equation to consider dependencies from ``v`` and ``p``
515     :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
516     """
517     pass
518     def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
519     """
520     assigns new values to the model parameters.
521    
522     :param f: external force
523     :type f: `Vector` object in `FunctionSpace` `Function` or similar
524     :param fixed_u_mask: mask of locations with fixed velocity.
525     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
526     :param eta: viscosity
527     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
528     :param surface_stress: normal surface stress
529     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
530     :param stress: initial stress
531     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
532     """
533     if eta !=None:
534     k=util.kronecker(self.domain.getDim())
535     kk=util.outer(k,k)
536 jfenwick 3432 self.eta=util.interpolate(eta, escript.Function(self.domain))
537 gross 2793 self.__pde_prec.setValue(D=1/self.eta)
538 gross 3502 self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
539 gross 2793 if restoration_factor!=None:
540     n=self.domain.getNormal()
541 gross 3502 self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
542 gross 2793 if fixed_u_mask!=None:
543 gross 3502 self.__pde_v.setValue(q=fixed_u_mask)
544 gross 2793 if f!=None: self.__f=f
545     if surface_stress!=None: self.__surface_stress=surface_stress
546     if stress!=None: self.__stress=stress
547    
548 jfenwick 3432 def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
549 gross 2208 """
550     assigns values to the model parameters
551 gross 2100
552 jfenwick 2625 :param f: external force
553     :type f: `Vector` object in `FunctionSpace` `Function` or similar
554     :param fixed_u_mask: mask of locations with fixed velocity.
555     :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
556     :param eta: viscosity
557     :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
558     :param surface_stress: normal surface stress
559     :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
560     :param stress: initial stress
561     :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
562 gross 2208 """
563 gross 2793 self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
564 gross 1414
565 gross 2719 def Bv(self,v,tol):
566 gross 2251 """
567     returns inner product of element p and div(v)
568 gross 1414
569 jfenwick 2625 :param v: a residual
570     :return: inner product of element p and div(v)
571     :rtype: ``float``
572 gross 2100 """
573 gross 2793 self.__pde_proj.setValue(Y=-util.div(v))
574 gross 2719 self.getSolverOptionsDiv().setTolerance(tol)
575     self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
576     out=self.__pde_proj.getSolution()
577     return out
578 gross 2208
579 gross 2445 def inner_pBv(self,p,Bv):
580     """
581     returns inner product of element p and Bv=-div(v)
582    
583 jfenwick 2625 :param p: a pressure increment
584     :param Bv: a residual
585     :return: inner product of element p and Bv=-div(v)
586     :rtype: ``float``
587 gross 2445 """
588 jfenwick 3432 return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
589 gross 2445
590 gross 2251 def inner_p(self,p0,p1):
591 gross 2100 """
592 gross 2251 Returns inner product of p0 and p1
593 gross 1414
594 jfenwick 2625 :param p0: a pressure
595     :param p1: a pressure
596     :return: inner product of p0 and p1
597     :rtype: ``float``
598 gross 2100 """
599 jfenwick 3432 s0=util.interpolate(p0, escript.Function(self.domain))
600     s1=util.interpolate(p1, escript.Function(self.domain))
601 gross 2100 return util.integrate(s0*s1)
602 artak 1550
603 gross 2251 def norm_v(self,v):
604 gross 2100 """
605 gross 2251 returns the norm of v
606 gross 2208
607 jfenwick 2625 :param v: a velovity
608     :return: norm of v
609     :rtype: non-negative ``float``
610 gross 2100 """
611 gross 2719 return util.sqrt(util.integrate(util.length(util.grad(v))**2))
612 gross 2100
613 gross 2793
614 gross 2719 def getDV(self, p, v, tol):
615 gross 1414 """
616 gross 2251 return the value for v for a given p (overwrite)
617    
618 jfenwick 2625 :param p: a pressure
619 gross 2719 :param v: a initial guess for the value v to return.
620     :return: dv given as *Adv=(f-Av-B^*p)*
621 gross 1414 """
622 gross 2793 self.updateStokesEquation(v,p)
623 gross 3502 self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
624 gross 2719 self.getSolverOptionsVelocity().setTolerance(tol)
625     self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
626 gross 2100 if self.__stress.isEmpty():
627 gross 3502 self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
628 gross 2100 else:
629 gross 3502 self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
630     out=self.__pde_v.getSolution()
631 gross 2208 return out
632 gross 1414
633 gross 2445 def norm_Bv(self,Bv):
634 gross 2251 """
635     Returns Bv (overwrite).
636    
637 jfenwick 2625 :rtype: equal to the type of p
638     :note: boundary conditions on p should be zero!
639 gross 2251 """
640 jfenwick 3432 return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
641 gross 2251
642 gross 2719 def solve_AinvBt(self,p, tol):
643 gross 2251 """
644 gross 2719 Solves *Av=B^*p* with accuracy `tol`
645 gross 2251
646 jfenwick 2625 :param p: a pressure increment
647     :return: the solution of *Av=B^*p*
648     :note: boundary conditions on v should be zero!
649 gross 2251 """
650 gross 3502 self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
651     out=self.__pde_v.getSolution()
652 gross 2251 return out
653    
654 gross 2719 def solve_prec(self,Bv, tol):
655 gross 2251 """
656 jfenwick 2625 applies preconditioner for for *BA^{-1}B^** to *Bv*
657     with accuracy `self.getSubProblemTolerance()`
658 gross 2251
659 jfenwick 2625 :param Bv: velocity increment
660     :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
661     :note: boundary conditions on p are zero.
662 gross 2251 """
663 gross 2445 self.__pde_prec.setValue(Y=Bv)
664 gross 2719 self.getSolverOptionsPressure().setTolerance(tol)
665     self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
666     out=self.__pde_prec.getSolution()
667     return out

  ViewVC Help
Powered by ViewVC 1.1.26