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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1414 - (show annotations)
Thu Feb 14 10:01:43 2008 UTC (12 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 8656 byte(s)
a first verion of a Stokes solver
1 #######################################################
2 #
3 # Copyright 2003-2007 by ACceSS MNRF
4 # Copyright 2007 by University of Queensland
5 #
6 # http://esscc.uq.edu.au
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 #######################################################
12 #
13
14 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
15 http://www.access.edu.au
16 Primary Business: Queensland, Australia"""
17 __license__="""Licensed under the Open Software License version 3.0
18 http://www.opensource.org/licenses/osl-3.0.php"""
19 import unittest
20 import tempfile
21
22 from esys.escript import *
23 from esys.finley import Rectangle
24 import sys
25 import os
26 try:
27 FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
28 except KeyError:
29 FINLEY_WORKDIR='.'
30
31
32 NE=6
33 TOL=1.e-5
34 VERBOSE=False
35
36 from esys.escript import *
37 from esys.escript.models import StokesProblemCartesian
38 from esys.finley import Rectangle, Brick
39 class Test_Simple2DModels(unittest.TestCase):
40 def setUp(self):
41 self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)
42 def tearDown(self):
43 del self.domain
44 def test_StokesProblemCartesian_P_0(self):
45 ETA=1.
46 P1=0.
47
48 x=self.domain.getX()
49 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
50 mask=whereZero(x[0]) * [1.,1.] \
51 +whereZero(x[0]-1) * [1.,1.] \
52 +whereZero(x[1]) * [1.,0.] \
53 +whereZero(x[1]-1) * [1.,1.]
54
55 sp=StokesProblemCartesian(self.domain)
56
57 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
58 u0=(1-x[0])*x[0]*[0.,1.]
59 p0=Scalar(P1,ReducedSolution(self.domain))
60 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100)
61
62 error_v0=Lsup(u[0]-u0[0])
63 error_v1=Lsup(u[1]-u0[1])/0.25
64 zz=P1*x[0]*x[1]-p
65 error_p=Lsup(zz-integrate(zz))
66 # print error_p, error_v0,error_v1
67 self.failUnless(error_p<TOL,"pressure error too large.")
68 self.failUnless(error_v0<TOL,"0-velocity error too large.")
69 self.failUnless(error_v1<TOL,"1-velocity error too large.")
70
71 def test_StokesProblemCartesian_P_small(self):
72 ETA=1.
73 P1=1.
74
75 x=self.domain.getX()
76 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
77 mask=whereZero(x[0]) * [1.,1.] \
78 +whereZero(x[0]-1) * [1.,1.] \
79 +whereZero(x[1]) * [1.,0.] \
80 +whereZero(x[1]-1) * [1.,1.]
81
82 sp=StokesProblemCartesian(self.domain)
83
84 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
85 u0=(1-x[0])*x[0]*[0.,1.]
86 p0=Scalar(P1,ReducedSolution(self.domain))
87 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100)
88
89 error_v0=Lsup(u[0]-u0[0])
90 error_v1=Lsup(u[1]-u0[1])/0.25
91 zz=P1*x[0]*x[1]-p
92 error_p=Lsup(zz-integrate(zz))
93 # print error_p, error_v0,error_v1
94 self.failUnless(error_p<TOL,"pressure error too large.")
95 self.failUnless(error_v0<TOL,"0-velocity error too large.")
96 self.failUnless(error_v1<TOL,"1-velocity error too large.")
97
98 def test_StokesProblemCartesian_P_large(self):
99 ETA=1.
100 P1=1000.
101
102 x=self.domain.getX()
103 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
104 mask=whereZero(x[0]) * [1.,1.] \
105 +whereZero(x[0]-1) * [1.,1.] \
106 +whereZero(x[1]) * [1.,0.] \
107 +whereZero(x[1]-1) * [1.,1.]
108
109 sp=StokesProblemCartesian(self.domain)
110
111 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
112 u0=(1-x[0])*x[0]*[0.,1.]
113 p0=Scalar(P1,ReducedSolution(self.domain))
114 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100)
115
116 error_v0=Lsup(u[0]-u0[0])
117 error_v1=Lsup(u[1]-u0[1])/0.25
118 zz=P1*x[0]*x[1]-p
119 error_p=Lsup(zz-integrate(zz))/P1
120 # print error_p, error_v0,error_v1
121 self.failUnless(error_p<TOL,"pressure error too large.")
122 self.failUnless(error_v0<TOL,"0-velocity error too large.")
123 self.failUnless(error_v1<TOL,"1-velocity error too large.")
124
125
126 class Test_Simple3DModels(unittest.TestCase):
127 def setUp(self):
128 self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)
129 def tearDown(self):
130 del self.domain
131 def test_StokesProblemCartesian_P_0(self):
132 ETA=1.
133 P1=0.
134
135 x=self.domain.getX()
136 F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
137 x=self.domain.getX()
138 mask=whereZero(x[0]) * [1.,1.,1.] \
139 +whereZero(x[0]-1) * [1.,1.,1.] \
140 +whereZero(x[1]) * [1.,1.,1.] \
141 +whereZero(x[1]-1) * [1.,1.,1.] \
142 +whereZero(x[2]) * [1.,1.,0.] \
143 +whereZero(x[2]-1) * [1.,1.,1.]
144
145
146 sp=StokesProblemCartesian(self.domain)
147
148 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
149 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
150 p0=Scalar(P1,ReducedSolution(self.domain))
151 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100)
152
153 error_v0=Lsup(u[0]-u0[0])
154 error_v1=Lsup(u[1]-u0[1])
155 error_v2=Lsup(u[2]-u0[2])/0.25**2
156 zz=P1*x[0]*x[1]*x[2]-p
157 error_p=Lsup(zz-integrate(zz))
158 # print error_p, error_v0,error_v1,error_v2
159 self.failUnless(error_p<TOL,"pressure error too large.")
160 self.failUnless(error_v0<TOL,"0-velocity error too large.")
161 self.failUnless(error_v1<TOL,"1-velocity error too large.")
162 self.failUnless(error_v2<TOL,"2-velocity error too large.")
163 def test_StokesProblemCartesian_P_small(self):
164 ETA=1.
165 P1=1.
166
167 x=self.domain.getX()
168 F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
169 mask=whereZero(x[0]) * [1.,1.,1.] \
170 +whereZero(x[0]-1) * [1.,1.,1.] \
171 +whereZero(x[1]) * [1.,1.,1.] \
172 +whereZero(x[1]-1) * [1.,1.,1.] \
173 +whereZero(x[2]) * [1.,1.,0.] \
174 +whereZero(x[2]-1) * [1.,1.,1.]
175
176
177 sp=StokesProblemCartesian(self.domain)
178
179 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
180 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
181 p0=Scalar(P1,ReducedSolution(self.domain))
182 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100)
183
184 error_v0=Lsup(u[0]-u0[0])
185 error_v1=Lsup(u[1]-u0[1])
186 error_v2=Lsup(u[2]-u0[2])/0.25**2
187 zz=P1*x[0]*x[1]*x[2]-p
188 error_p=Lsup(zz-integrate(zz))/P1
189 # print error_p, error_v0,error_v1,error_v2
190 self.failUnless(error_p<TOL,"pressure error too large.")
191 self.failUnless(error_v0<TOL,"0-velocity error too large.")
192 self.failUnless(error_v1<TOL,"1-velocity error too large.")
193 self.failUnless(error_v2<TOL,"2-velocity error too large.")
194 def test_StokesProblemCartesian_P_large(self):
195 ETA=1.
196 P1=1000.
197
198 x=self.domain.getX()
199 F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
200 mask=whereZero(x[0]) * [1.,1.,1.] \
201 +whereZero(x[0]-1) * [1.,1.,1.] \
202 +whereZero(x[1]) * [1.,1.,1.] \
203 +whereZero(x[1]-1) * [1.,1.,1.] \
204 +whereZero(x[2]) * [1.,1.,0.] \
205 +whereZero(x[2]-1) * [1.,1.,1.]
206
207
208 sp=StokesProblemCartesian(self.domain)
209
210 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
211 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
212 p0=Scalar(P1,ReducedSolution(self.domain))
213 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100)
214
215 error_v0=Lsup(u[0]-u0[0])
216 error_v1=Lsup(u[1]-u0[1])
217 error_v2=Lsup(u[2]-u0[2])/0.25**2
218 zz=P1*x[0]*x[1]*x[2]-p
219 error_p=Lsup(zz-integrate(zz))/P1
220 # print error_p, error_v0,error_v1,error_v2
221 self.failUnless(error_p<TOL,"pressure error too large.")
222 self.failUnless(error_v0<TOL,"0-velocity error too large.")
223 self.failUnless(error_v1<TOL,"1-velocity error too large.")
224 self.failUnless(error_v2<TOL,"2-velocity error too large.")
225
226 if __name__ == '__main__':
227 suite = unittest.TestSuite()
228 suite.addTest(unittest.makeSuite(Test_Simple2DModels))
229 suite.addTest(unittest.makeSuite(Test_Simple3DModels))
230 s=unittest.TextTestRunner(verbosity=2).run(suite)
231 if not s.wasSuccessful(): sys.exit(1)
232

  ViewVC Help
Powered by ViewVC 1.1.26