/[escript]/trunk/pyvisi/test/python/run_escript_with_lazy_evaluation.py
ViewVC logotype

Annotation of /trunk/pyvisi/test/python/run_escript_with_lazy_evaluation.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1166 - (hide annotations)
Thu May 24 00:56:00 2007 UTC (12 years, 8 months ago) by jongui
File MIME type: text/x-python
File size: 7655 byte(s)
Updated some test.
1 jongui 1151 # Import the necessary modules.
2     from esys.escript import *
3     from esys.escript.pdetools import Locator
4     from esys.escript.linearPDEs import LinearPDE
5     from esys.finley import Rectangle, Brick
6     from numarray import identity,zeros,ones
7     from esys.pyvisi import Scene, DataCollector, Map, Velocity, Ellipsoid, Camera
8     from esys.pyvisi.constant import *
9     import unittest, os
10     from stat import ST_SIZE
11    
12     try:
13     PYVISI_WORKDIR=os.environ['PYVISI_WORKDIR']
14     except KeyError:
15     PYVISI_WORKDIR='.'
16     try:
17     PYVISI_TEST_DATA_ROOT=os.environ['PYVISI_TEST_DATA_ROOT']
18     except KeyError:
19     PYVISI_TEST_DATA_ROOT='.'
20    
21     PYVISI_TEST_ESCRIPT_REFERENCE_IMAGES_PATH = os.path.join(PYVISI_TEST_DATA_ROOT,\
22     "data_reference_images", "escript")
23     PYVISI_TEST_ESCRIPT_IMAGES_PATH = os.path.join(PYVISI_WORKDIR, \
24     "data_sample_images", "escript")
25    
26     MIN_IMAGE_SIZE = 100
27     X_SIZE = 400
28     Y_SIZE = 400
29     JPG_RENDERER = Renderer.OFFLINE_JPG
30    
31     class TestEscript:
32     def render(self, file):
33     self.scene.render(image_name = \
34     os.path.join(PYVISI_TEST_ESCRIPT_IMAGES_PATH, file))
35    
36     self.failUnless(os.stat(os.path.join(PYVISI_TEST_ESCRIPT_IMAGES_PATH, \
37     file))[ST_SIZE] > MIN_IMAGE_SIZE)
38    
39     class TestEscriptMap(unittest.TestCase, TestEscript):
40 jongui 1158 def tearDown(self):
41     del self.scene
42    
43 jongui 1151 def testEscriptMap(self):
44     #... set some parameters ...
45 jongui 1159 xc=[0.02,0.002]
46     r=0.001
47     qc=50.e6
48     Tref=0.
49     rhocp=2.6e6
50     eta=75.
51     kappa=240.
52     tend=5.
53 jongui 1151 # ... time, time step size and counter ...
54     t=0
55 jongui 1159 h=0.1
56 jongui 1151 i=0
57 jongui 1158
58 jongui 1151 #... generate domain ...
59     mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)
60     #... open PDE ...
61 jongui 1159 mypde=LinearPDE(mydomain)
62     mypde.setSymmetryOn()
63     mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref)
64 jongui 1151 # ... set heat source: ....
65 jongui 1159 x=mydomain.getX()
66     qH=qc*whereNegative(length(x-xc)-r)
67 jongui 1151 # ... set initial temperature ....
68 jongui 1159 T=Tref
69 jongui 1151
70     # Create a Scene.
71     s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)
72     self.scene = s
73    
74     # Create a DataCollector reading directly from escript objects.
75     dc = DataCollector(source = Source.ESCRIPT)
76    
77     # Create a Map.
78 jongui 1166 m = Map(scene = s, data_collector = dc, \
79     viewport = Viewport.SOUTH_WEST, lut = Lut.COLOR, \
80 jongui 1151 cell_to_point = False, outline = True)
81    
82     # Create a Camera.
83     c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
84    
85     # ... start iteration:
86 jongui 1166 while t<0.4:
87 jongui 1151 i+=1
88 jongui 1159 t+=h
89     mypde.setValue(Y=qH+rhocp/h*T)
90     T=mypde.getSolution()
91 jongui 1151
92     dc.setData(temp = T)
93    
94     # Render the object.
95     self.render("TestEscriptMap%02d.jpg" % i)
96    
97     class TestEscriptVelocity(unittest.TestCase, TestEscript):
98 jongui 1158 def tearDown(self):
99     del self.scene
100    
101 jongui 1151 def testEscriptVelocity(self):
102     #... set some parameters ...
103 jongui 1159 lam=1.
104     mu=0.1
105     alpha=1.e-6
106     xc=[0.3,0.3,1.]
107     beta=8.
108     T_ref=0.
109     T_0=1.
110 jongui 1166
111 jongui 1151 #... generate domain ...
112     mydomain = Brick(l0=1.,l1=1., l2=1.,n0=10, n1=10, n2=10)
113 jongui 1159 x=mydomain.getX()
114 jongui 1151 #... set temperature ...
115 jongui 1159 T=T_0*exp(-beta*length(x-xc))
116 jongui 1151 #... open symmetric PDE ...
117 jongui 1159 mypde=LinearPDE(mydomain)
118     mypde.setSymmetryOn()
119 jongui 1166
120 jongui 1151 #... set coefficients ...
121 jongui 1159 C=Tensor4(0.,Function(mydomain))
122     for i in range(mydomain.getDim()):
123 jongui 1166 for j in range(mydomain.getDim()):
124     C[i,i,j,j]+=lam
125     C[j,i,j,i]+=mu
126     C[j,i,i,j]+=mu
127    
128 jongui 1159 msk=whereZero(x[0])*[1.,0.,0.] \
129 jongui 1166 +whereZero(x[1])*[0.,1.,0.] \
130     +whereZero(x[2])*[0.,0.,1.]
131 jongui 1159 sigma0=(lam+2./3.*mu)*alpha*(T-T_ref)*kronecker(mydomain)
132     mypde.setValue(A=C,X=sigma0,q=msk)
133 jongui 1166
134 jongui 1151 #... solve pde ...
135 jongui 1159 u=mypde.getSolution()
136 jongui 1166 #... calculate von-Misses
137 jongui 1159 g=grad(u)
138     sigma=mu*(g+transpose(g))+lam*trace(g)*kronecker(mydomain)-sigma0
139 jongui 1166 sigma_mises=sqrt(((sigma[0,0]-sigma[1,1])**2+\
140     (sigma[1,1]-sigma[2,2])**2+ \
141 jongui 1159 (sigma[2,2]-sigma[0,0])**2)/6. \
142     +sigma[0,1]**2 + sigma[1,2]**2 + sigma[2,0]**2)
143 jongui 1166
144 jongui 1151 # Create a Scene.
145     s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)
146     self.scene = s
147    
148     # Create a DataCollector reading directly from an escript object.
149     dc = DataCollector(source = Source.ESCRIPT)
150 jongui 1159 dc.setData(disp = u, stress = sigma_mises)
151 jongui 1151
152     # Create a Velocity.
153 jongui 1166 v = Velocity(scene = s, data_collector = dc, \
154     viewport = Viewport.SOUTH_WEST, \
155     arrow = Arrow.THREE_D, color_mode = ColorMode.SCALAR, \
156 jongui 1151 lut = Lut.COLOR, cell_to_point = True, outline = True)
157     v.setScaleFactor(scale_factor = 0.3)
158    
159     # Create a Camera.
160     c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
161     c.isometricView()
162    
163     # Render the object.
164     self.render("TestEscriptVelocity.jpg")
165    
166     ne=32 # number of cells in x_0 and x_1 directions
167     width=10000. # length in x_0 and x_1 directions
168 jongui 1159 lam=3.462e9
169     mu=3.462e9
170     rho=1154.
171     tend=60
172     h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
173 jongui 1151
174 jongui 1159 U0=0.01 # amplitude of point source
175 jongui 1151 class TestEscriptEllipsoid(unittest.TestCase, TestEscript):
176 jongui 1158 def tearDown(self):
177     del self.scene
178    
179 jongui 1151 def testEscriptEllipsoid(self):
180     mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
181 jongui 1159 self.wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
182 jongui 1151
183 jongui 1159 def wavePropagation(self,domain,h,tend,lam,mu,rho,U0):
184     x=domain.getX()
185 jongui 1151 # ... open new PDE ...
186 jongui 1159 mypde=LinearPDE(domain)
187     mypde.setSolverMethod(LinearPDE.LUMPING)
188     kronecker=identity(mypde.getDim())
189 jongui 1151
190     # spherical source at middle of bottom face
191    
192 jongui 1159 xc=[width/2.,width/2.,0.]
193 jongui 1151 # define small radius around point xc
194     # Lsup(x) returns the maximum value of the argument x
195 jongui 1159 src_radius = 0.1*Lsup(domain.getSize())
196     dunit=numarray.array([1.,0.,0.]) # defines direction of point source
197 jongui 1151
198 jongui 1159 mypde.setValue(D=kronecker*rho)
199 jongui 1151 # ... set initial values ....
200 jongui 1159 n=0
201 jongui 1166 # initial value of displacement at point source is constant (U0=0.01)
202     # for first two time steps
203 jongui 1159 u=U0*whereNegative(length(x-xc)-src_radius)*dunit
204     u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit
205 jongui 1151 t=0
206    
207     # define the location of the point source
208 jongui 1159 L=Locator(domain,numarray.array(xc))
209 jongui 1151 # find potential at point source
210 jongui 1159 u_pc=L.getValue(u)
211 jongui 1151
212 jongui 1159 u_pc_x = u_pc[0]
213     u_pc_y = u_pc[1]
214     u_pc_z = u_pc[2]
215 jongui 1151
216     # open file to save displacement at point source
217     #u_pc_data=open('./data/U_pc.out','w')
218     #u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
219    
220     # Create a Scene.
221     s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)
222     self.scene = s
223    
224     # Create a DataCollector reading directly from escript objects.
225     dc = DataCollector(source = Source.ESCRIPT)
226    
227     # Create an Ellipsoid.
228 jongui 1159 e = Ellipsoid(scene = s, data_collector = dc,
229 jongui 1166 viewport = Viewport.SOUTH_WEST,
230     lut = Lut.COLOR, cell_to_point = True, outline = True)
231 jongui 1159 e.setScaleFactor(scale_factor = 0.7)
232 jongui 1166 e.setMaxScaleFactor(max_scale_factor = 1000)
233 jongui 1159 e.setRatio(ratio = 10)
234 jongui 1151
235     # Create a Camera.
236     c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
237     c.isometricView()
238    
239 jongui 1166 while t<0.5:
240 jongui 1151 # ... get current stress ....
241 jongui 1159 g=grad(u)
242     stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
243 jongui 1151 # ... get new acceleration ....
244 jongui 1159 mypde.setValue(X=-stress)
245     a=mypde.getSolution()
246 jongui 1151 # ... get new displacement ...
247 jongui 1159 u_new=2*u-u_last+h**2*a
248 jongui 1151 # ... shift displacements ....
249 jongui 1159 u_last=u
250     u=u_new
251     t+=h
252     n+=1
253     u_pc=L.getValue(u)
254 jongui 1151
255 jongui 1159 u_pc_x=u_pc[0]
256     u_pc_y=u_pc[1]
257     u_pc_z=u_pc[2]
258 jongui 1151
259 jongui 1166 # ... save current acceleration in units of gravity and
260     # displacements
261 jongui 1159 if n==1 or n%10==0:
262 jongui 1166 dc.setData(acceleration = length(a)/9.81, displacement = u,
263     tensor = stress, Ux = u[0])
264 jongui 1159
265 jongui 1166 # Render the object.
266     self.render("TestEscriptEllipsoid%02d.jpg" % (n/10))
267 jongui 1151
268    
269     if __name__ == '__main__':
270     suite = unittest.TestSuite()
271     suite.addTest(unittest.TestLoader().loadTestsFromTestCase(TestEscriptMap))
272     suite.addTest(unittest.TestLoader().loadTestsFromTestCase(TestEscriptVelocity))
273     suite.addTest(unittest.TestLoader().loadTestsFromTestCase(TestEscriptEllipsoid))
274     unittest.TextTestRunner(verbosity=2).run(suite)
275    

  ViewVC Help
Powered by ViewVC 1.1.26