doc/examples/wave_with_pyvisi.py

from esys.escript import *
from esys.escript.pdetools import Locator
from esys.escript.linearPDEs import LinearPDE
from esys.finley import Brick
from numarray import identity,zeros,ones
from esys.pyvisi import Scene, DataCollector, Ellipsoid, Camera
from esys.pyvisi.constant import *

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)
print "time step size = ",h

U0=0.01 # amplitude of point source

def wavePropagation(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())
   print "src_radius = ",src_radius

   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)
   print "u at point charge=",u_pc
  
   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))
 
   s = Scene(renderer = Renderer.OFFLINE_JPG, x_size = 500, y_size = 500)
   dc = DataCollector(source = Source.ESCRIPT)

   while t<tend:
     # ... 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
     print n,"-th time step t ",t
     u_pc=L.getValue(u)
     print "u at point charge=",u_pc
     
     u_pc_x=u_pc[0]
     u_pc_y=u_pc[1]
     u_pc_z=u_pc[2]
      
     # save displacements at point source to file for t > 0
     u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
 
     # ... save current acceleration in units of gravity and displacements 
     if n==1 or n%10==0: 
     #saveVTK("./data/usoln.%i.vtu"%(n/10), acceleration=length(a)/9.81, 
         #     displacement = u, tensor = stress, Ux = u[0])

         dc.setData(acceleration = length(a)/9.81, displacement = u, 
                 tensor = stress, Ux = u[0])

         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)

         c = Camera(scene = s, data_collector = dc, 
                 viewport = Viewport.SOUTH_WEST)
         c.isometricView()

         s.render(image_name = "wave_%02d.jpg" % (n/10))

   u_pc_data.close()
  
mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
wavePropagation(mydomain,h,tend,lam,mu,rho,U0)

1	from esys.escript import *
2	from esys.escript.pdetools import Locator
3	from esys.escript.linearPDEs import LinearPDE
4	from esys.finley import Brick
5	from numarray import identity,zeros,ones
6	from esys.pyvisi import Scene, DataCollector, Ellipsoid, Camera
7	from esys.pyvisi.constant import *
8
9	ne=32 # number of cells in x_0 and x_1 directions
10	width=10000. # length in x_0 and x_1 directions
11	lam=3.462e9
12	mu=3.462e9
13	rho=1154.
14	tend=60
15	h=(1./5.)sqrt(rho/(lam+2mu))*(width/ne)
16	print "time step size = ",h
17
18	U0=0.01 # amplitude of point source
19
20	def wavePropagation(domain,h,tend,lam,mu,rho,U0):
21	x=domain.getX()
22	# ... open new PDE ...
23	mypde=LinearPDE(domain)
24	mypde.setSolverMethod(LinearPDE.LUMPING)
25	kronecker=identity(mypde.getDim())
26
27	# spherical source at middle of bottom face
28
29	xc=[width/2.,width/2.,0.]
30	# define small radius around point xc
31	# Lsup(x) returns the maximum value of the argument x
32	src_radius = 0.1*Lsup(domain.getSize())
33	print "src_radius = ",src_radius
34
35	dunit=numarray.array([1.,0.,0.]) # defines direction of point source
36
37	mypde.setValue(D=kronecker*rho)
38	# ... set initial values ....
39	n=0
40	# initial value of displacement at point source is constant (U0=0.01)
41	# for first two time steps
42	u=U0whereNegative(length(x-xc)-src_radius)dunit
43	u_last=U0whereNegative(length(x-xc)-src_radius)dunit
44	t=0
45
46	# define the location of the point source
47	L=Locator(domain,numarray.array(xc))
48	# find potential at point source
49	u_pc=L.getValue(u)
50	print "u at point charge=",u_pc
51
52	u_pc_x = u_pc[0]
53	u_pc_y = u_pc[1]
54	u_pc_z = u_pc[2]
55
56	# open file to save displacement at point source
57	u_pc_data=open('./data/U_pc.out','w')
58	u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
59
60	s = Scene(renderer = Renderer.OFFLINE_JPG, x_size = 500, y_size = 500)
61	dc = DataCollector(source = Source.ESCRIPT)
62
63	while t<tend:
64	# ... get current stress ....
65	g=grad(u)
66	stress=lamtrace(g)kronecker+mu*(g+transpose(g))
67	# ... get new acceleration ....
68	mypde.setValue(X=-stress)
69	a=mypde.getSolution()
70	# ... get new displacement ...
71	u_new=2u-u_last+h2a
72	# ... shift displacements ....
73	u_last=u
74	u=u_new
75	t+=h
76	n+=1
77	print n,"-th time step t ",t
78	u_pc=L.getValue(u)
79	print "u at point charge=",u_pc
80
81	u_pc_x=u_pc[0]
82	u_pc_y=u_pc[1]
83	u_pc_z=u_pc[2]
84
85	# save displacements at point source to file for t > 0
86	u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
87
88	# ... save current acceleration in units of gravity and displacements
89	if n==1 or n%10==0:
90	#saveVTK("./data/usoln.%i.vtu"%(n/10), acceleration=length(a)/9.81,
91	# displacement = u, tensor = stress, Ux = u[0])
92
93	dc.setData(acceleration = length(a)/9.81, displacement = u,
94	tensor = stress, Ux = u[0])
95
96	e = Ellipsoid(scene = s, data_collector = dc,
97	viewport = Viewport.SOUTH_WEST,
98	lut = Lut.COLOR, cell_to_point = True, outline = True)
99	e.setScaleFactor(scale_factor = 0.7)
100	e.setMaxScaleFactor(max_scale_factor = 1000)
101	e.setRatio(ratio = 10)
102
103	c = Camera(scene = s, data_collector = dc,
104	viewport = Viewport.SOUTH_WEST)
105	c.isometricView()
106
107	s.render(image_name = "wave_%02d.jpg" % (n/10))
108
109	u_pc_data.close()
110
111	mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
112	wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
113