/[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 878 - (show annotations)
Wed Oct 25 03:42:56 2006 UTC (12 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 4422 byte(s)
bigger mesh
1 # $Id: pdetools.py 867 2006-10-09 06:50:09Z gross $
2
3 """
4 Provides some tools related to PDEs.
5
6 Currently includes:
7 - Projector - to project a discontinuous
8 - Locator - to trace values in data objects at a certain location
9 - TimeIntegrationManager - to handel extraplotion in time
10 - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
11
12 @var __author__: name of author
13 @var __copyright__: copyrights
14 @var __license__: licence agreement
15 @var __url__: url entry point on documentation
16 @var __version__: version
17 @var __date__: date of the version
18 """
19
20 __author__="Lutz Gross, l.gross@uq.edu.au"
21 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
22 http://www.access.edu.au
23 Primary Business: Queensland, Australia"""
24 __license__="""Licensed under the Open Software License version 3.0
25 http://www.opensource.org/licenses/osl-3.0.php"""
26 __url__="http://www.iservo.edu.au/esys"
27 __version__="$Revision$"
28 __date__="$Date:$"
29
30 from esys.escript import *
31 from esys.escript.pdetools import SaddlePointProblem
32 from esys.escript.linearPDEs import LinearPDE
33 from esys.finley import Rectangle
34
35 class SimpleStokesProblem(SaddlePointProblem):
36 """
37 simple example of saddle point problem
38 """
39 def __init__(self,domain):
40 super(SimpleStokesProblem, self).__init__(self)
41
42 self.__pde_u=LinearPDE(domain)
43 self.__pde_u.setSymmetryOn()
44 self.__pde_u.setValue(A=identityTensor4(dom))
45
46 self.__pde_p=LinearPDE(domain)
47 self.__pde_p.setReducedOrderOn()
48 self.__pde_p.setSymmetryOn()
49 self.__pde_p.setValue(D=1.)
50
51 def initialize(self,f=Data(),fixed_u_mask=Data()):
52 self.__pde_u.setValue(q=fixed_u_mask,Y=f)
53 def inner(self,p0,p1):
54 return integrate(p0*p1,Function(self.__pde_p.getDomain()))
55
56 def solve_f(self,u,p,tol=1.e-8):
57 self.__pde_u.setTolerance(tol)
58 self.__pde_u.setValue(X=grad(u)+p*kronecker(self.__pde_u.getDomain()))
59 return self.__pde_u.getSolution()
60 def solve_g(self,u,tol=1.e-8):
61 self.__pde_p.setTolerance(tol)
62 self.__pde_p.setValue(X=-u)
63 dp=self.__pde_p.getSolution()
64 return dp
65
66 class StokesProblem(SaddlePointProblem):
67 """
68 simple example of saddle point problem
69 """
70 def __init__(self,domain):
71 super(StokesProblem, self).__init__(self)
72 self.domain=domain
73 self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
74 self.__pde_u.setSymmetryOn()
75
76 self.__pde_p=LinearPDE(domain)
77 self.__pde_p.setReducedOrderOn()
78 self.__pde_p.setSymmetryOn()
79
80 def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1):
81 self.eta=eta
82 A =self.__pde_u.createCoefficientOfGeneralPDE("A")
83 for i in range(self.domain.getDim()):
84 for j in range(self.domain.getDim()):
85 A[i,j,j,i] += self.eta
86 A[i,j,i,j] += self.eta
87 self.__pde_p.setValue(D=1./self.eta)
88 self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=f)
89
90 def inner(self,p0,p1):
91 return integrate(p0*p1,Function(self.__pde_p.getDomain()))
92
93 def solve_f(self,u,p,tol=1.e-8):
94 self.__pde_u.setTolerance(tol)
95 g=grad(u)
96 self.__pde_u.setValue(X=self.eta*symmetric(g)+p*kronecker(self.__pde_u.getDomain()))
97 return self.__pde_u.getSolution()
98
99 def solve_g(self,u,tol=1.e-8):
100 self.__pde_p.setTolerance(tol)
101 self.__pde_p.setValue(X=-u)
102 dp=self.__pde_p.getSolution()
103 return dp
104
105 NE=50
106 dom=Rectangle(NE,NE,order=2)
107 # prop=SimpleStokesProblem(dom)
108 prop=StokesProblem(dom)
109 x=dom.getX()
110 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)
111 u0=Vector(0.,Solution(dom))
112 u0[0]=x[1]*whereZero(x[1]-1.)
113 p0=Scalar(0,ReducedSolution(dom))
114 # prop.initialize(fixed_u_mask=mask)
115 prop.initialize(fixed_u_mask=mask,eta=10.)
116 u,p=prop.solve(u0,p0,tolerance=0.01)
117 # saveVTK("stokes.xml",u=u,p=p,m=mask,u0=u0)
118
119 eta=whereNegative(x[1]-0.5)*1.e6+whereNonNegative(x[1]-0.5)
120 prop.initialize(fixed_u_mask=mask,eta=eta)
121 u,p=prop.solve(u0,p0,tolerance=0.01,tolerance_u=0.1,accepted_reduction=0.8)
122 saveVTK("stokes.xml",u=u,p=p,m=mask,u0=u0)
123
124 # vim: expandtab shiftwidth=4:

  ViewVC Help
Powered by ViewVC 1.1.26