# Annotation of /trunk/doc/examples/usersguide/wave.py

Revision 2474 - (hide annotations)
Tue Jun 16 06:32:15 2009 UTC (9 years, 8 months ago) by gross
Original Path: trunk/doc/examples/wave.py
File MIME type: text/x-python
File size: 3663 byte(s)
```linearPDEs has is now using the SolverOptions class to talk to PASO
```
 1 ksteube 1215 2 ksteube 1811 ######################################################## 3 # 4 # Copyright (c) 2003-2008 by University of Queensland 5 # Earth Systems Science Computational Center (ESSCC) 6 7 # 8 # Primary Business: Queensland, Australia 9 # Licensed under the Open Software License version 3.0 10 11 # 12 ######################################################## 13 14 __copyright__="""Copyright (c) 2003-2008 by University of Queensland 15 Earth Systems Science Computational Center (ESSCC) 16 http://www.uq.edu.au/esscc 17 Primary Business: Queensland, Australia""" 18 __license__="""Licensed under the Open Software License version 3.0 19 20 jfenwick 2344 __url__= 21 ksteube 1811 22 ksteube 1215 # You can shorten the execution time by reducing variable tend from 60 to 0.5 23 24 jgs 108 from esys.escript import * 25 lkettle 582 from esys.escript.pdetools import Locator 26 gross 327 from esys.escript.linearPDEs import LinearPDE 27 jgs 110 from esys.finley import Brick 28 jfenwick 2455 from numpy import identity,zeros,ones 29 lkettle 582 30 ksteube 1215 if not os.path.isdir("data"): 31 print "\nCreating subdirectory 'data'\n" 32 os.mkdir("data") 33 34 lkettle 582 ne=32 # number of cells in x_0 and x_1 directions 35 width=10000. # length in x_0 and x_1 directions 36 jgs 110 lam=3.462e9 37 mu=3.462e9 38 rho=1154. 39 tend=60 40 lkettle 582 h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne) 41 jgs 110 print "time step size = ",h 42 jgs 108 43 lkettle 582 U0=0.01 # amplitude of point source 44 jgs 110 45 lkettle 582 def wavePropagation(domain,h,tend,lam,mu,rho,U0): 46 jgs 108 x=domain.getX() 47 # ... open new PDE ... 48 jgs 110 mypde=LinearPDE(domain) 49 gross 2474 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING) 50 lkettle 582 kronecker=identity(mypde.getDim()) 51 52 # spherical source at middle of bottom face 53 54 xc=[width/2.,width/2.,0.] 55 # define small radius around point xc 56 # Lsup(x) returns the maximum value of the argument x 57 src_radius = 0.1*Lsup(domain.getSize()) 58 print "src_radius = ",src_radius 59 60 jfenwick 2455 dunit=numpy.array([1.,0.,0.]) # defines direction of point source 61 lkettle 582 62 mypde.setValue(D=kronecker*rho) 63 jgs 110 # ... set initial values .... 64 jgs 108 n=0 65 lkettle 582 # initial value of displacement at point source is constant (U0=0.01) 66 # for first two time steps 67 u=U0*whereNegative(length(x-xc)-src_radius)*dunit 68 u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit 69 jgs 110 t=0 70 lkettle 582 71 # define the location of the point source 72 jfenwick 2455 L=Locator(domain,numpy.array(xc)) 73 lkettle 582 # find potential at point source 74 u_pc=L.getValue(u) 75 print "u at point charge=",u_pc 76 77 u_pc_x = u_pc[0] 78 u_pc_y = u_pc[1] 79 u_pc_z = u_pc[2] 80 81 # open file to save displacement at point source 82 gross 2420 u_pc_data=FileWriter('./data/U_pc.out') 83 lkettle 582 u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z)) 84 85 jgs 108 while t

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision