/[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 1414 - (hide annotations)
Thu Feb 14 10:01:43 2008 UTC (11 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 4318 byte(s)
a first verion of a Stokes solver
1 gross 1414 # $Id:$
2     #
3     #######################################################
4     #
5     # Copyright 2008 by University of Queensland
6     #
7     # http://esscc.uq.edu.au
8     # 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     #
12     #######################################################
13     #
14    
15     """
16     Some models for flow
17    
18     @var __author__: name of author
19     @var __copyright__: copyrights
20     @var __license__: licence agreement
21     @var __url__: url entry point on documentation
22     @var __version__: version
23     @var __date__: date of the version
24     """
25    
26     __author__="Lutz Gross, l.gross@uq.edu.au"
27     __copyright__=""" Copyright (c) 2008 by ACcESS MNRF
28     http://www.access.edu.au
29     Primary Business: Queensland, Australia"""
30     __license__="""Licensed under the Open Software License version 3.0
31     http://www.opensource.org/licenses/osl-3.0.php"""
32     __url__="http://www.iservo.edu.au/esys"
33     __version__="$Revision:$"
34     __date__="$Date:$"
35    
36     from escript import *
37     import util
38     from linearPDEs import LinearPDE
39     from pdetools import HomogeneousSaddlePointProblem
40    
41     class StokesProblemCartesian(HomogeneousSaddlePointProblem):
42     """
43     solves
44    
45     -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i
46     u_{i,i}=0
47    
48     u=0 where fixed_u_mask>0
49     eta*(u_{i,j}+u_{j,i})*n_j=surface_stress
50    
51     if surface_stress is not give 0 is assumed.
52    
53     typical usage:
54    
55     sp=StokesProblemCartesian(domain)
56     sp.setTolerance()
57     sp.initialize(...)
58     v,p=sp.solve(v0,p0)
59     """
60     def __init__(self,domain,**kwargs):
61     HomogeneousSaddlePointProblem.__init__(self,**kwargs)
62     self.domain=domain
63     self.vol=util.integrate(1.,Function(self.domain))
64     self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
65     self.__pde_u.setSymmetryOn()
66     self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
67    
68     self.__pde_prec=LinearPDE(domain)
69     self.__pde_prec.setReducedOrderOn()
70     self.__pde_prec.setSymmetryOn()
71    
72     self.__pde_proj=LinearPDE(domain)
73     self.__pde_proj.setReducedOrderOn()
74     self.__pde_proj.setSymmetryOn()
75     self.__pde_proj.setValue(D=1.)
76    
77     def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):
78     self.eta=eta
79     A =self.__pde_u.createCoefficientOfGeneralPDE("A")
80     self.__pde_u.setValue(A=Data())
81     for i in range(self.domain.getDim()):
82     for j in range(self.domain.getDim()):
83     A[i,j,j,i] += 1.
84     A[i,j,i,j] += 1.
85     self.__pde_prec.setValue(D=1./self.eta)
86     self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)
87    
88     def B(self,arg):
89     d=util.div(arg)
90     self.__pde_proj.setValue(Y=d)
91     self.__pde_proj.setTolerance(self.getSubProblemTolerance())
92     return self.__pde_proj.getSolution(verbose=self.show_details)
93    
94     def inner(self,p0,p1):
95     s0=util.interpolate(p0,Function(self.domain))
96     s1=util.interpolate(p1,Function(self.domain))
97     return util.integrate(s0*s1)
98    
99     def getStress(self,u):
100     mg=util.grad(u)
101     return 2.*self.eta*util.symmetric(mg)
102    
103     def solve_A(self,u,p):
104     """
105     solves Av=f-Au-B^*p (v=0 on fixed_u_mask)
106     """
107     self.__pde_u.setTolerance(self.getSubProblemTolerance())
108     self.__pde_u.setValue(X=-self.getStress(u)+p*util.kronecker(self.domain))
109     return self.__pde_u.getSolution(verbose=self.show_details)
110    
111     def solve_prec(self,p):
112     self.__pde_prec.setTolerance(self.getSubProblemTolerance())
113     self.__pde_prec.setValue(Y=p)
114     q=self.__pde_prec.getSolution(verbose=self.show_details)
115     return q
116     def stoppingcriterium(self,Bv,v,p):
117     n_r=util.sqrt(self.inner(Bv,Bv))
118     n_v=util.Lsup(v)
119     if self.verbose: print "PCG step %s: L2(div(v)) = %s, Lsup(v)=%s"%(self.iter,n_r,n_v)
120     self.iter+=1
121     if n_r <= self.vol**(1./2.-1./self.domain.getDim())*n_v*self.getTolerance():
122     if self.verbose: print "PCG terminated after %s steps."%self.iter
123     return True
124     else:
125     return False

  ViewVC Help
Powered by ViewVC 1.1.26