/[escript]/trunk/doc/examples/cookbook/cblib/wavesolver2d.py
ViewVC logotype

Contents of /trunk/doc/examples/cookbook/cblib/wavesolver2d.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2667 - (show annotations)
Thu Sep 17 01:49:11 2009 UTC (10 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 4087 byte(s)
Renamed the main cookbook tex file to match our convention.
Replaced doc/cookbook/figures/heatrefraction002contqu.pdf with
a version which is actually pdf. However it needs to be regnerated since
it it sideways.

The examples have had their copyright notices fixed (dates were too early).
sb2.py has been removed since it uses pyvisi.

scons will now build the cookbook as parts of a docs build.
Also in reposnse to :
scons cookbook_pdf


1
2 ########################################################
3 #
4 # Copyright (c) 2009 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2009 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 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 # You can shorten the execution time by reducing variable tend from 60 to 0.5
23
24 # Importing all the necessary modules required.
25 from esys.escript import *
26 from esys.escript.pdetools import Locator
27 from esys.escript.linearPDEs import LinearPDE
28 from esys.finley import Rectangle
29 from numpy import identity,zeros,ones
30 import os
31 from phones import *
32
33 ########################################################
34 # subroutine: wavesolver2d
35 # Can solve a generic 2D wave propagation problem with a
36 # point source in a homogeneous medium.
37 # Arguments:
38 # domain : domain to solve over
39 # h : time step
40 # tend : end time
41 # lam, mu : lames constants for domain
42 # rho : density of domain
43 # U0 : magnitude of source
44 # xc : source location in domain (Vector)
45 # savepath: where to output the data files
46 ########################################################
47 def wavesolver2d(domain,h,tend,lam,mu,rho,U0,xc,savepath):
48 x=domain.getX()
49 # ... open new PDE ...
50 mypde=LinearPDE(domain)
51 #mypde.setSolverMethod(LinearPDE.LUMPING)
52 mypde.setSymmetryOn()
53 kmat = kronecker(domain)
54 mypde.setValue(D=kmat*rho)
55
56 # define small radius around point xc
57 # Lsup(x) returns the maximum value of the argument x
58 src_radius = 50#2*Lsup(domain.getSize())
59 print "src_radius = ",src_radius
60
61 dunit=numpy.array([0.,1.]) # defines direction of point source
62
63
64 # ... set initial values ....
65 n=0
66 # initial value of displacement at point source is constant (U0=0.01)
67 # for first two time steps
68 u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)*dunit
69 u_m1=u
70 t=0
71
72 u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2)
73 u_pc_x1 = u_pot[0,0]
74 u_pc_y1 = u_pot[0,1]
75 u_pc_x2 = u_pot[1,0]
76 u_pc_y2 = u_pot[1,1]
77 u_pc_x3 = u_pot[2,0]
78 u_pc_y3 = u_pot[2,1]
79
80 # open file to save displacement at point source
81 u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
82 u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
83
84 while t<tend:
85 # ... get current stress ....
86
87 ##OLD WAY
88 g=grad(u)
89 stress=lam*trace(g)*kmat+mu*(g+transpose(g))
90 ### ... get new acceleration ....
91 #mypde.setValue(X=-stress)
92 #a=mypde.getSolution()
93 ### ... get new displacement ...
94 #u_p1=2*u-u_m1+h*h*a
95 ###NEW WAY
96 mypde.setValue(X=-stress*(h*h),Y=(rho*2*u-rho*u_m1))
97 u_p1 = mypde.getSolution()
98 # ... shift displacements ....
99 u_m1=u
100 u=u_p1
101 #stress =
102 t+=h
103 n+=1
104 print n,"-th time step t ",t
105 u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)
106
107 # print "u at point charge=",u_pc
108 u_pc_x1 = u_pot[0,0]
109 u_pc_y1 = u_pot[0,1]
110 u_pc_x2 = u_pot[1,0]
111 u_pc_y2 = u_pot[1,1]
112 u_pc_x3 = u_pot[2,0]
113 u_pc_y3 = u_pot[2,1]
114
115 # save displacements at point source to file for t > 0
116 u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
117
118 # ... save current acceleration in units of gravity and displacements
119 #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
120 #displacement = length(u), tensor = stress, Ux = u[0] )
121 saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
122
123 u_pc_data.close()

  ViewVC Help
Powered by ViewVC 1.1.26