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

Contents of /trunk-mpi-branch/finley/test/python/stokes_problems.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 968 - (show annotations)
Tue Feb 13 22:56:57 2007 UTC (12 years, 9 months ago) by ksteube
File MIME type: text/x-python
File size: 4079 byte(s)
Branch for MPI solution of implicit problems
1 # $Id$
2
3 """
4 solvers for the stokes problem
5
6 @var __author__: name of author
7 @var __copyright__: copyrights
8 @var __license__: licence agreement
9 @var __url__: url entry point on documentation
10 @var __version__: version
11 @var __date__: date of the version
12 """
13
14 __author__="Lutz Gross, l.gross@uq.edu.au"
15 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
16 http://www.access.edu.au
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__="http://www.iservo.edu.au/esys"
21 __version__="$Revision$"
22 __date__="$Date$"
23
24 from esys.escript import *
25 from esys.escript.pdetools import SaddlePointProblem
26 from esys.escript.linearPDEs import LinearPDE
27 from esys.finley import Rectangle
28
29 class SimpleStokesProblem(SaddlePointProblem):
30 """
31 simple example of saddle point problem
32 """
33 def __init__(self,domain):
34 super(SimpleStokesProblem, self).__init__(self)
35
36 self.__pde_u=LinearPDE(domain)
37 self.__pde_u.setSymmetryOn()
38 self.__pde_u.setValue(A=identityTensor4(dom))
39
40 self.__pde_p=LinearPDE(domain)
41 self.__pde_p.setReducedOrderOn()
42 self.__pde_p.setSymmetryOn()
43 self.__pde_p.setValue(D=1.)
44
45 def initialize(self,f=Data(),fixed_u_mask=Data()):
46 self.__pde_u.setValue(q=fixed_u_mask,Y=f)
47 def inner(self,p0,p1):
48 return integrate(p0*p1,Function(self.__pde_p.getDomain()))
49
50 def solve_f(self,u,p,tol=1.e-8):
51 self.__pde_u.setTolerance(tol)
52 self.__pde_u.setValue(X=grad(u)+p*kronecker(self.__pde_u.getDomain()))
53 return self.__pde_u.getSolution()
54 def solve_g(self,u,tol=1.e-8):
55 self.__pde_p.setTolerance(tol)
56 self.__pde_p.setValue(X=-u)
57 dp=self.__pde_p.getSolution()
58 return dp
59
60 class StokesProblem(SaddlePointProblem):
61 """
62 simple example of saddle point problem
63 """
64 def __init__(self,domain):
65 super(StokesProblem, self).__init__(self)
66 self.domain=domain
67 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
68 self.__pde_u.setSymmetryOn()
69
70 self.__pde_p=LinearPDE(domain)
71 self.__pde_p.setReducedOrderOn()
72 self.__pde_p.setSymmetryOn()
73
74 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1):
75 self.eta=eta
76 A =self.__pde_u.createCoefficientOfGeneralPDE("A")
77 for i in range(self.domain.getDim()):
78 for j in range(self.domain.getDim()):
79 A[i,j,j,i] += self.eta
80 A[i,j,i,j] += self.eta
81 self.__pde_p.setValue(D=1./self.eta)
82 self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=f)
83
84 def inner(self,p0,p1):
85 return integrate(p0*p1,Function(self.__pde_p.getDomain()))
86
87 def solve_f(self,u,p,tol=1.e-8):
88 self.__pde_u.setTolerance(tol)
89 g=grad(u)
90 self.__pde_u.setValue(X=self.eta*symmetric(g)+p*kronecker(self.__pde_u.getDomain()))
91 return self.__pde_u.getSolution()
92
93 def solve_g(self,u,tol=1.e-8):
94 self.__pde_p.setTolerance(tol)
95 self.__pde_p.setValue(X=-u)
96 dp=self.__pde_p.getSolution()
97 return dp
98
99 NE=50
100 dom=Rectangle(NE,NE,order=2)
101 # prop=SimpleStokesProblem(dom)
102 prop=StokesProblem(dom)
103 x=dom.getX()
104 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)
105 u0=Vector(0.,Solution(dom))
106 u0[0]=x[1]*whereZero(x[1]-1.)
107 p0=Scalar(0,ReducedSolution(dom))
108 # prop.initialize(fixed_u_mask=mask)
109 prop.initialize(fixed_u_mask=mask,eta=10.)
110 u,p=prop.solve(u0,p0,tolerance=0.01)
111 # saveVTK("stokes.xml",u=u,p=p,m=mask,u0=u0)
112
113 eta=whereNegative(x[1]-0.5)*1.e6+whereNonNegative(x[1]-0.5)
114 prop.initialize(fixed_u_mask=mask,eta=eta)
115 u,p=prop.solve(u0,p0,tolerance=0.01,tolerance_u=0.1,accepted_reduction=0.8)
116 saveVTK("stokes.xml",u=u,p=p,m=mask,u0=u0)
117
118 # vim: expandtab shiftwidth=4:

  ViewVC Help
Powered by ViewVC 1.1.26