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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1158 by jongui, Tue May 22 04:24:01 2007 UTC revision 1159 by jongui, Tue May 22 05:09:31 2007 UTC
# Line 42  class TestEscriptMap(unittest.TestCase, Line 42  class TestEscriptMap(unittest.TestCase,
42    
43      def testEscriptMap(self):      def testEscriptMap(self):
44          #... set some parameters ...          #... set some parameters ...
45          #xc=[0.02,0.002]          xc=[0.02,0.002]
46          #r=0.001          r=0.001
47          #qc=50.e6          qc=50.e6
48          #Tref=0.          Tref=0.
49          #rhocp=2.6e6          rhocp=2.6e6
50          #eta=75.          eta=75.
51          #kappa=240.          kappa=240.
52          #tend=5.          tend=5.
53          # ... time, time step size and counter ...          # ... time, time step size and counter ...
54          t=0          t=0
55          #h=0.1          h=0.1
56          i=0          i=0
57    
58          #... generate domain ...          #... generate domain ...
59          mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)          mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)
60          #... open PDE ...          #... open PDE ...
61          #mypde=LinearPDE(mydomain)          mypde=LinearPDE(mydomain)
62          #mypde.setSymmetryOn()          mypde.setSymmetryOn()
63          #mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref)          mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref)
64          # ... set heat source: ....          # ... set heat source: ....
65          #x=mydomain.getX()          x=mydomain.getX()
66          #qH=qc*whereNegative(length(x-xc)-r)          qH=qc*whereNegative(length(x-xc)-r)
67          # ... set initial temperature ....          # ... set initial temperature ....
68          #T=Tref          T=Tref
69    
70          # Create a Scene.          # Create a Scene.
71          s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)          s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)
# Line 84  class TestEscriptMap(unittest.TestCase, Line 84  class TestEscriptMap(unittest.TestCase,
84    
85          # ... start iteration:          # ... start iteration:
86          while t<0.5:          while t<0.5:
             t+=0.1  
87              i+=1              i+=1
88              #t+=h              t+=h
89              #mypde.setValue(Y=qH+rhocp/h*T)              mypde.setValue(Y=qH+rhocp/h*T)
90              #T=mypde.getSolution()              T=mypde.getSolution()
91    
             #======= syntetic data =========  
             T=length(mydomain.getX())  
92              dc.setData(temp = T)              dc.setData(temp = T)
93    
94              # Render the object.              # Render the object.
# Line 103  class TestEscriptVelocity(unittest.TestC Line 100  class TestEscriptVelocity(unittest.TestC
100    
101      def testEscriptVelocity(self):      def testEscriptVelocity(self):
102          #... set some parameters ...          #... set some parameters ...
103          #lam=1.          lam=1.
104          #mu=0.1          mu=0.1
105          #alpha=1.e-6          alpha=1.e-6
106          #xc=[0.3,0.3,1.]          xc=[0.3,0.3,1.]
107          #beta=8.          beta=8.
108          #T_ref=0.          T_ref=0.
109          #T_0=1.          T_0=1.
110          #... generate domain ...          #... generate domain ...
111          mydomain = Brick(l0=1.,l1=1., l2=1.,n0=10, n1=10, n2=10)          mydomain = Brick(l0=1.,l1=1., l2=1.,n0=10, n1=10, n2=10)
112          #x=mydomain.getX()          x=mydomain.getX()
113          #... set temperature ...          #... set temperature ...
114          #T=T_0*exp(-beta*length(x-xc))          T=T_0*exp(-beta*length(x-xc))
115          #... open symmetric PDE ...          #... open symmetric PDE ...
116          #mypde=LinearPDE(mydomain)          mypde=LinearPDE(mydomain)
117          #mypde.setSymmetryOn()          mypde.setSymmetryOn()
118          #... set coefficients ...          #... set coefficients ...
119          #C=Tensor4(0.,Function(mydomain))          C=Tensor4(0.,Function(mydomain))
120          #for i in range(mydomain.getDim()):          for i in range(mydomain.getDim()):
121          #  for j in range(mydomain.getDim()):            for j in range(mydomain.getDim()):
122          #    C[i,i,j,j]+=lam               C[i,i,j,j]+=lam
123          #    C[j,i,j,i]+=mu               C[j,i,j,i]+=mu
124          #    C[j,i,i,j]+=mu               C[j,i,i,j]+=mu
125          #msk=whereZero(x[0])*[1.,0.,0.] \          msk=whereZero(x[0])*[1.,0.,0.] \
126          #   +whereZero(x[1])*[0.,1.,0.] \             +whereZero(x[1])*[0.,1.,0.] \
127          #   +whereZero(x[2])*[0.,0.,1.]             +whereZero(x[2])*[0.,0.,1.]
128          #sigma0=(lam+2./3.*mu)*alpha*(T-T_ref)*kronecker(mydomain)          sigma0=(lam+2./3.*mu)*alpha*(T-T_ref)*kronecker(mydomain)
129          #mypde.setValue(A=C,X=sigma0,q=msk)          mypde.setValue(A=C,X=sigma0,q=msk)
130          #... solve pde ...          #... solve pde ...
131          #u=mypde.getSolution()          u=mypde.getSolution()
132          #==============          #==============
133          #... calculate von-Misses          # ... calculate von-Misses
134          #u=mydomain.getX()          u=mydomain.getX()
135          #g=grad(u)          g=grad(u)
136          #sigma=mu*(g+transpose(g))+lam*trace(g)*kronecker(mydomain)-sigma0          sigma=mu*(g+transpose(g))+lam*trace(g)*kronecker(mydomain)-sigma0
137          #sigma_mises=\          sigma_mises=\
138          #       sqrt(((sigma[0,0]-sigma[1,1])**2+(sigma[1,1]-sigma[2,2])**2+ \                  sqrt(((sigma[0,0]-sigma[1,1])**2+(sigma[1,1]-sigma[2,2])**2+ \
139          #       (sigma[2,2]-sigma[0,0])**2)/6. \                  (sigma[2,2]-sigma[0,0])**2)/6. \
140          #       +sigma[0,1]**2 + sigma[1,2]**2 + sigma[2,0]**2)                  +sigma[0,1]**2 + sigma[1,2]**2 + sigma[2,0]**2)
141          #          
142          # Create a Scene.          # Create a Scene.
143          s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)          s = Scene(renderer = JPG_RENDERER, x_size = X_SIZE, y_size = Y_SIZE)
144          self.scene = s          self.scene = s
# Line 149  class TestEscriptVelocity(unittest.TestC Line 146  class TestEscriptVelocity(unittest.TestC
146          # Create a DataCollector reading directly from an escript object.          # Create a DataCollector reading directly from an escript object.
147          dc = DataCollector(source = Source.ESCRIPT)          dc = DataCollector(source = Source.ESCRIPT)
148          T=length(mydomain.getX())          T=length(mydomain.getX())
149          #dc.setData(disp = u, stress = sigma_mises)          dc.setData(disp = u, stress = sigma_mises)
         dc.setData(disp = T, stress = T)  
150    
151          # Create a Velocity.          # Create a Velocity.
152          v = Velocity(scene = s, data_collector = dc,          v = Velocity(scene = s, data_collector = dc,
# Line 168  class TestEscriptVelocity(unittest.TestC Line 164  class TestEscriptVelocity(unittest.TestC
164    
165  ne=32 # number of cells in x_0 and x_1 directions  ne=32 # number of cells in x_0 and x_1 directions
166  width=10000. # length in x_0 and x_1 directions  width=10000. # length in x_0 and x_1 directions
167  #lam=3.462e9  lam=3.462e9
168  #mu=3.462e9  mu=3.462e9
169  #rho=1154.  rho=1154.
170  #tend=60  tend=60
171  #h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)  h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
172    
173  #U0=0.01 # amplitude of point source  U0=0.01 # amplitude of point source
174  class TestEscriptEllipsoid(unittest.TestCase, TestEscript):  class TestEscriptEllipsoid(unittest.TestCase, TestEscript):
175      def tearDown(self):      def tearDown(self):
176          del self.scene          del self.scene
177    
178      def testEscriptEllipsoid(self):      def testEscriptEllipsoid(self):
179          mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)          mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
180          #self.wavePropagation(mydomain,h,tend,lam,mu,rho,U0)          self.wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
         self.wavePropagation(mydomain)  
181    
182      def wavePropagation(self, mydomain):      def wavePropagation(self,domain,h,tend,lam,mu,rho,U0):
183          #x=domain.getX()          x=domain.getX()
184          # ... open new PDE ...          # ... open new PDE ...
185          #mypde=LinearPDE(domain)          mypde=LinearPDE(domain)
186          #mypde.setSolverMethod(LinearPDE.LUMPING)          mypde.setSolverMethod(LinearPDE.LUMPING)
187          #kronecker=identity(mypde.getDim())          kronecker=identity(mypde.getDim())
188    
189          #  spherical source at middle of bottom face          #  spherical source at middle of bottom face
190    
191          #xc=[width/2.,width/2.,0.]          xc=[width/2.,width/2.,0.]
192          # define small radius around point xc          # define small radius around point xc
193          # Lsup(x) returns the maximum value of the argument x          # Lsup(x) returns the maximum value of the argument x
194          #src_radius = 0.1*Lsup(domain.getSize())          src_radius = 0.1*Lsup(domain.getSize())
195          #dunit=numarray.array([1.,0.,0.]) # defines direction of point source          dunit=numarray.array([1.,0.,0.]) # defines direction of point source
196    
197          #mypde.setValue(D=kronecker*rho)          mypde.setValue(D=kronecker*rho)
198          # ... set initial values ....          # ... set initial values ....
199          #n=0          n=0
200          # initial value of displacement at point source is constant          # initial value of displacement at point source is constant
201          # (U0=0.01) for first two time steps          # (U0=0.01) for first two time steps
202          #u=U0*whereNegative(length(x-xc)-src_radius)*dunit          u=U0*whereNegative(length(x-xc)-src_radius)*dunit
203          #u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit          u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit
204          t=0          t=0
205    
206          # define the location of the point source          # define the location of the point source
207          #L=Locator(domain,numarray.array(xc))          L=Locator(domain,numarray.array(xc))
208          # find potential at point source          # find potential at point source
209          #u_pc=L.getValue(u)          u_pc=L.getValue(u)
210    
211          #u_pc_x = u_pc[0]          u_pc_x = u_pc[0]
212          #u_pc_y = u_pc[1]          u_pc_y = u_pc[1]
213          #u_pc_z = u_pc[2]          u_pc_z = u_pc[2]
214    
215          # open file to save displacement at point source          # open file to save displacement at point source
216          #u_pc_data=open('./data/U_pc.out','w')          #u_pc_data=open('./data/U_pc.out','w')
# Line 229  class TestEscriptEllipsoid(unittest.Test Line 224  class TestEscriptEllipsoid(unittest.Test
224          dc = DataCollector(source = Source.ESCRIPT)          dc = DataCollector(source = Source.ESCRIPT)
225    
226          # Create an Ellipsoid.          # Create an Ellipsoid.
227          e = Map(scene = s, data_collector = dc,          e = Ellipsoid(scene = s, data_collector = dc,
228                 viewport = Viewport.SOUTH_WEST,                 viewport = Viewport.SOUTH_WEST,
229                 lut = Lut.COLOR, cell_to_point = True, outline = True)                 lut = Lut.COLOR, cell_to_point = True, outline = True)
230          #e.setScaleFactor(scale_factor = 0.7)          e.setScaleFactor(scale_factor = 0.7)
231          #e.setMaxScaleFactor(max_scale_factor = 1000)          e.setMaxScaleFactor(max_scale_factor = 1000)
232          #e.setRatio(ratio = 10)          e.setRatio(ratio = 10)
233    
234          # Create a Camera.          # Create a Camera.
235          c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)          c = Camera(scene = s, viewport = Viewport.SOUTH_WEST)
236          c.isometricView()          c.isometricView()
237    
238          while t<0.4:          while t<0.4:
             t+=0.1  
239              # ... get current stress ....              # ... get current stress ....
240              #g=grad(u)              g=grad(u)
241              #stress=lam*trace(g)*kronecker+mu*(g+transpose(g))              stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
242              # ... get new acceleration ....              # ... get new acceleration ....
243              #mypde.setValue(X=-stress)                        mypde.setValue(X=-stress)          
244              #a=mypde.getSolution()              a=mypde.getSolution()
245              # ... get new displacement ...              # ... get new displacement ...
246              #u_new=2*u-u_last+h**2*a              u_new=2*u-u_last+h**2*a
247              # ... shift displacements ....              # ... shift displacements ....
248              #u_last=u              u_last=u
249              #u=u_new              u=u_new
250              #t+=h              t+=h
251              #n+=1              n+=1
252              #u_pc=L.getValue(u)              u_pc=L.getValue(u)
253    
254              #u_pc_x=u_pc[0]              u_pc_x=u_pc[0]
255              #u_pc_y=u_pc[1]              u_pc_y=u_pc[1]
256              #u_pc_z=u_pc[2]              u_pc_z=u_pc[2]
257    
258              # ... save current acceleration in units of gravity and displacements              # ... save current acceleration in units of gravity and displacements
259              # === syntetic data              # === syntetic data
260              u=domain.getX()              u=domain.getX()
261              #T=length(mydomain.getX())              if n==1 or n%10==0:
262              #if n==1 or n%10==0:                  dc.setData(acceleration = length(a)/9.81, displacement = u,
263               #dc.setData(acceleration = length(a)/9.81, displacement = u,                       tensor = stress, Ux = u[0])
264              #        tensor = stress, Ux = u[0])  
265              dc.setData(acceleration = length(u), displacement = u, tensor = grad(u), Ux = u[0])               # Render the object.
               # Render the object.  
266              self.render("TestEscriptEllipsoid%02d.jpg" % (t/10))              self.render("TestEscriptEllipsoid%02d.jpg" % (t/10))
267    
268    

Legend:
Removed from v.1158  
changed lines
  Added in v.1159

  ViewVC Help
Powered by ViewVC 1.1.26