test/python/run_escript_with_lazy_evaluation.py

# Import the necessary modules.
from esys.escript import *
from esys.escript.pdetools import Locator
from esys.escript.linearPDEs import LinearPDE
from esys.finley import Rectangle, Brick
from numarray import identity,zeros,ones
from esys.pyvisi import Scene, DataCollector, Map, Velocity, Ellipsoid, Camera
from esys.pyvisi.constant import *
import unittest, os
from stat import ST_SIZE

try:
    PYVISI_WORKDIR=os.environ['PYVISI_WORKDIR']
except KeyError:
    PYVISI_WORKDIR='.'
try:
    PYVISI_TEST_DATA_ROOT=os.environ['PYVISI_TEST_DATA_ROOT']
except KeyError:
    PYVISI_TEST_DATA_ROOT='.'

PYVISI_TEST_ESCRIPT_REFERENCE_IMAGES_PATH = os.path.join(PYVISI_TEST_DATA_ROOT,\
        "data_reference_images", "escript")
PYVISI_TEST_ESCRIPT_IMAGES_PATH = os.path.join(PYVISI_WORKDIR, \
        "data_sample_images", "escript")

MIN_IMAGE_SIZE = 100
X_SIZE = 400
Y_SIZE = 400
JPG_RENDERER = Renderer.OFFLINE_JPG

class TestEscript:
    def tearDown(self):
        self.scene

    def render(self, file):
        self.scene.render(image_name = \
                os.path.join(PYVISI_TEST_ESCRIPT_IMAGES_PATH, file))

        self.failUnless(os.stat(os.path.join(PYVISI_TEST_ESCRIPT_IMAGES_PATH, \
                file))[ST_SIZE] > MIN_IMAGE_SIZE)

class TestEscriptMap(unittest.TestCase, TestEscript):
    def testEscriptMap(self):
        #... set some parameters ...
        xc=[0.02,0.002]
        r=0.001
        qc=50.e6
        Tref=0.
        rhocp=2.6e6
        eta=75.
        kappa=240.
        tend=5.
        # ... time, time step size and counter ...
        t=0
        h=0.1
        i=0
        #... generate domain ...
        mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)
        #... open PDE ...
        mypde=LinearPDE(mydomain)
        mypde.setSymmetryOn()
        mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref)
        # ... set heat source: ....
        x=mydomain.getX()
        qH=qc*whereNegative(length(x-xc)-r)
        # ... set initial temperature ....
        T=Tref

        # Create a Scene.
        s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)
        self.scene = s

        # Create a DataCollector reading directly from escript objects.
        dc = DataCollector(source = Source.ESCRIPT)

        # Create a Map.
        m = Map(scene = s, data_collector = dc,
                viewport = Viewport.SOUTH_WEST, lut = Lut.COLOR,
                cell_to_point = False, outline = True)

        # Create a Camera.
        c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)

        # ... start iteration:
        while t<0.4:
            i+=1
            t+=h
            mypde.setValue(Y=qH+rhocp/h*T)
            T=mypde.getSolution()

            dc.setData(temp = T)

            # Render the object.
            self.render("TestEscriptMap%02d.jpg" % i)

class TestEscriptVelocity(unittest.TestCase, TestEscript):
    def testEscriptVelocity(self):
        #... set some parameters ...
        lam=1.
        mu=0.1
        alpha=1.e-6
        xc=[0.3,0.3,1.]
        beta=8.
        T_ref=0.
        T_0=1.
        #... generate domain ...
        mydomain = Brick(l0=1.,l1=1., l2=1.,n0=10, n1=10, n2=10)
        x=mydomain.getX()
        #... set temperature ...
        T=T_0*exp(-beta*length(x-xc))
        #... open symmetric PDE ...
        mypde=LinearPDE(mydomain)
        mypde.setSymmetryOn()
        #... set coefficients ...
        C=Tensor4(0.,Function(mydomain))
        for i in range(mydomain.getDim()):
          for j in range(mydomain.getDim()):
             C[i,i,j,j]+=lam
             C[j,i,j,i]+=mu
             C[j,i,i,j]+=mu
        msk=whereZero(x[0])*[1.,0.,0.] \
           +whereZero(x[1])*[0.,1.,0.] \
           +whereZero(x[2])*[0.,0.,1.]
        sigma0=(lam+2./3.*mu)*alpha*(T-T_ref)*kronecker(mydomain)
        mypde.setValue(A=C,X=sigma0,q=msk)
        #... solve pde ...
        u=mypde.getSolution()
        #... calculate von-Misses
        g=grad(u)
        sigma=mu*(g+transpose(g))+lam*trace(g)*kronecker(mydomain)-sigma0
        sigma_mises=\
                sqrt(((sigma[0,0]-sigma[1,1])**2+(sigma[1,1]-sigma[2,2])**2+ \
                (sigma[2,2]-sigma[0,0])**2)/6. \
                +sigma[0,1]**2 + sigma[1,2]**2 + sigma[2,0]**2)
         
        # Create a Scene.
        s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)
        self.scene = s

        # Create a DataCollector reading directly from an escript object.
        dc = DataCollector(source = Source.ESCRIPT)
        dc.setData(disp = u, stress = sigma_mises)

        # Create a Velocity.
        v = Velocity(scene = s, data_collector = dc, 
                viewport = Viewport.SOUTH_WEST,
                arrow = Arrow.THREE_D, color_mode = ColorMode.SCALAR,
                lut = Lut.COLOR, cell_to_point = True, outline = True)
        v.setScaleFactor(scale_factor = 0.3)

        # Create a Camera.
        c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
        c.isometricView()

        # Render the object.
        self.render("TestEscriptVelocity.jpg")

ne=32 # number of cells in x_0 and x_1 directions
width=10000. # length in x_0 and x_1 directions
lam=3.462e9
mu=3.462e9
rho=1154.
tend=60
h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)

U0=0.01 # amplitude of point source
class TestEscriptEllipsoid(unittest.TestCase, TestEscript):
    def testEscriptEllipsoid(self):
        mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
        self.wavePropagation(mydomain,h,tend,lam,mu,rho,U0)

    def tearDown(self):
        ne
        width
        lam
        mu
        rho
        tend
        h

        U0

    def wavePropagation(self, domain,h,tend,lam,mu,rho,U0):
        x=domain.getX()
        # ... open new PDE ...
        mypde=LinearPDE(domain)
        mypde.setSolverMethod(LinearPDE.LUMPING)
        kronecker=identity(mypde.getDim())

        #  spherical source at middle of bottom face

        xc=[width/2.,width/2.,0.]
        # define small radius around point xc
        # Lsup(x) returns the maximum value of the argument x
        src_radius = 0.1*Lsup(domain.getSize())
        dunit=numarray.array([1.,0.,0.]) # defines direction of point source

        mypde.setValue(D=kronecker*rho)
        # ... set initial values ....
        n=0
        # initial value of displacement at point source is constant 
        # (U0=0.01) for first two time steps
        u=U0*whereNegative(length(x-xc)-src_radius)*dunit
        u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit
        t=0

        # define the location of the point source 
        L=Locator(domain,numarray.array(xc))
        # find potential at point source
        u_pc=L.getValue(u)

        u_pc_x = u_pc[0]
        u_pc_y = u_pc[1]
        u_pc_z = u_pc[2]

        # open file to save displacement at point source
        #u_pc_data=open('./data/U_pc.out','w')
        #u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))

        # Create a Scene.
        s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)
        self.scene = s

        # Create a DataCollector reading directly from escript objects.
        dc = DataCollector(source = Source.ESCRIPT)

        # Create an Ellipsoid.
        e = Ellipsoid(scene = s, data_collector = dc, 
               viewport = Viewport.SOUTH_WEST,
               lut = Lut.COLOR, cell_to_point = True, outline = True)
        e.setScaleFactor(scale_factor = 0.7)
        e.setMaxScaleFactor(max_scale_factor = 1000)
        e.setRatio(ratio = 10)

        # Create a Camera.
        c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
        c.isometricView()

        while t<0.4:
            # ... get current stress ....
            g=grad(u)
            stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
            # ... get new acceleration ....
            mypde.setValue(X=-stress)          
            a=mypde.getSolution()
            # ... get new displacement ...
            u_new=2*u-u_last+h**2*a
            # ... shift displacements ....
            u_last=u
            u=u_new
            t+=h
            n+=1
            u_pc=L.getValue(u)

            u_pc_x=u_pc[0]
            u_pc_y=u_pc[1]
            u_pc_z=u_pc[2]

            # ... save current acceleration in units of gravity and displacements 
            if n==1 or n%10==0: 
             dc.setData(acceleration = length(a)/9.81, displacement = u, 
                     tensor = stress, Ux = u[0])
              # Render the object.
             self.render("TestEscriptEllipsoid%02d.jpg" % (n/10))


if __name__ == '__main__':
    suite = unittest.TestSuite()
    suite.addTest(unittest.TestLoader().loadTestsFromTestCase(TestEscriptMap))
    suite.addTest(unittest.TestLoader().loadTestsFromTestCase(TestEscriptVelocity))
    suite.addTest(unittest.TestLoader().loadTestsFromTestCase(TestEscriptEllipsoid))
    unittest.TextTestRunner(verbosity=2).run(suite)

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=kappakronecker(mydomain),D=rhocp/h,d=eta,y=etaTref)
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_0exp(-betalength(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))+lamtrace(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+2mu))*(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=U0whereNegative(length(x-xc)-src_radius)dunit
204	u_last=U0whereNegative(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=lamtrace(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=2u-u_last+h2a
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