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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26