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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1414 - (hide annotations)
Thu Feb 14 10:01:43 2008 UTC (11 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 8656 byte(s)
a first verion of a Stokes solver
1 gross 1414 #######################################################
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