/[escript]/trunk/finley/test/python/stokes_problems.py
ViewVC logotype

Contents of /trunk/finley/test/python/stokes_problems.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (show annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 7 months ago) by jfenwick
File MIME type: text/x-python
File size: 4460 byte(s)
Don't panic.
Updating copyright stamps

1
2 ########################################################
3 #
4 # Copyright (c) 2003-2010 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
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 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 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 __url__="https://launchpad.net/escript-finley"
21
22 """
23 solvers for the stokes problem
24
25 :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 """
32
33 __author__="Lutz Gross, l.gross@uq.edu.au"
34
35 from esys.escript import *
36 from esys.escript.pdetools import SaddlePointProblem
37 from esys.escript.linearPDEs import LinearPDE
38 from esys.finley import Rectangle
39
40 class SimpleStokesProblem(SaddlePointProblem):
41 """
42 simple example of saddle point problem
43 """
44 def __init__(self,domain):
45 super(SimpleStokesProblem, self).__init__(self)
46
47 self.__pde_u=LinearPDE(domain)
48 self.__pde_u.setSymmetryOn()
49 self.__pde_u.setValue(A=identityTensor4(dom))
50
51 self.__pde_p=LinearPDE(domain)
52 self.__pde_p.setReducedOrderOn()
53 self.__pde_p.setSymmetryOn()
54 self.__pde_p.setValue(D=1.)
55
56 def initialize(self,f=Data(),fixed_u_mask=Data()):
57 self.__pde_u.setValue(q=fixed_u_mask,Y=f)
58 def inner(self,p0,p1):
59 return integrate(p0*p1,Function(self.__pde_p.getDomain()))
60
61 def solve_f(self,u,p,tol=1.e-8):
62 self.__pde_u.setTolerance(tol)
63 self.__pde_u.setValue(X=grad(u)+p*kronecker(self.__pde_u.getDomain()))
64 return self.__pde_u.getSolution()
65 def solve_g(self,u,tol=1.e-8):
66 self.__pde_p.setTolerance(tol)
67 self.__pde_p.setValue(X=-u)
68 dp=self.__pde_p.getSolution()
69 return dp
70
71 class StokesProblem(SaddlePointProblem):
72 """
73 simple example of saddle point problem
74 """
75 def __init__(self,domain):
76 super(StokesProblem, self).__init__(self)
77 self.domain=domain
78 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
79 self.__pde_u.setSymmetryOn()
80
81 self.__pde_p=LinearPDE(domain)
82 self.__pde_p.setReducedOrderOn()
83 self.__pde_p.setSymmetryOn()
84
85 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1):
86 self.eta=eta
87 A =self.__pde_u.createCoefficientOfGeneralPDE("A")
88 for i in range(self.domain.getDim()):
89 for j in range(self.domain.getDim()):
90 A[i,j,j,i] += self.eta
91 A[i,j,i,j] += self.eta
92 self.__pde_p.setValue(D=1./self.eta)
93 self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=f)
94
95 def inner(self,p0,p1):
96 return integrate(p0*p1,Function(self.__pde_p.getDomain()))
97
98 def solve_f(self,u,p,tol=1.e-8):
99 self.__pde_u.setTolerance(tol)
100 g=grad(u)
101 self.__pde_u.setValue(X=self.eta*symmetric(g)+p*kronecker(self.__pde_u.getDomain()))
102 return self.__pde_u.getSolution()
103
104 def solve_g(self,u,tol=1.e-8):
105 self.__pde_p.setTolerance(tol)
106 self.__pde_p.setValue(X=-u)
107 dp=self.__pde_p.getSolution()
108 return dp
109
110 NE=50
111 dom=Rectangle(NE,NE,order=2)
112 # prop=SimpleStokesProblem(dom)
113 prop=StokesProblem(dom)
114 x=dom.getX()
115 mask=(whereZero(x[0])+whereZero(x[0]-1.)+whereZero(x[1]-1.))*unitVector(0,dom)+(whereZero(x[1]-1.)+whereZero(x[1]))*unitVector(1,dom)
116 u0=Vector(0.,Solution(dom))
117 u0[0]=x[1]*whereZero(x[1]-1.)
118 p0=Scalar(0,ReducedSolution(dom))
119 # prop.initialize(fixed_u_mask=mask)
120 prop.initialize(fixed_u_mask=mask,eta=10.)
121 u,p=prop.solve(u0,p0,tolerance=0.01)
122 # saveVTK("stokes.vtu",u=u,p=p,m=mask,u0=u0)
123
124 eta=whereNegative(x[1]-0.5)*1.e6+whereNonNegative(x[1]-0.5)
125 prop.initialize(fixed_u_mask=mask,eta=eta)
126 u,p=prop.solve(u0,p0,tolerance=0.01,tolerance_u=0.1,accepted_reduction=0.8)
127 saveVTK("stokes.vtu",u=u,p=p,m=mask,u0=u0)
128
129 # vim: expandtab shiftwidth=4:

  ViewVC Help
Powered by ViewVC 1.1.26