/[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 1151 - (hide annotations)
Thu May 17 04:11:18 2007 UTC (13 years, 8 months ago) by jongui
File MIME type: text/x-python
File size: 7638 byte(s)
More update.
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 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