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

  ViewVC Help
Powered by ViewVC 1.1.26