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

Diff of /trunk/doc/examples/cookbook/example09a.py

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

revision 5287 by jfenwick, Wed Apr 9 05:41:57 2014 UTC revision 5288 by sshaw, Tue Dec 2 23:18:40 2014 UTC
# Line 1  Line 1 
1  from __future__ import division  from __future__ import division, print_function
 from __future__ import print_function  
2  ##############################################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2009-2014 by University of Queensland  # Copyright (c) 2009-2014 by University of Queensland
# Line 32  __url__="https://launchpad.net/escript-f Line 31  __url__="https://launchpad.net/escript-f
31  import matplotlib  import matplotlib
32  matplotlib.use('agg') #It's just here for automated testing  matplotlib.use('agg') #It's just here for automated testing
33  from esys.escript import *  from esys.escript import *
 from esys.finley import Rectangle  
34  import os  import os
35  # smoothing operator  # smoothing operator
36  from esys.escript.pdetools import Projector, Locator  from esys.escript.pdetools import Projector, Locator
# Line 42  import numpy as np Line 40  import numpy as np
40  import pylab as pl  import pylab as pl
41  import matplotlib.cm as cm  import matplotlib.cm as cm
42  from esys.escript.linearPDEs import LinearPDE, SolverOptions  from esys.escript.linearPDEs import LinearPDE, SolverOptions
 from esys.finley import ReadMesh  
43  from esys.weipa import saveVTK  from esys.weipa import saveVTK
44    try:
45        # This imports the rectangle domain function
46        from esys.finley import Rectangle, ReadMesh
47        HAVE_FINLEY = True
48    except ImportError:
49        print("Finley module not available")
50        HAVE_FINLEY = False
51  ########################################################MPI WORLD CHECK  ########################################################MPI WORLD CHECK
52  if getMPISizeWorld() > 1:  if getMPISizeWorld() > 1:
53          import sys          import sys
54          print("This example will not run in an MPI world.")          print("This example will not run in an MPI world.")
55          sys.exit(0)          sys.exit(0)
56    
57  #################################################ESTABLISHING VARIABLES  if HAVE_FINLEY:
58  # where to save output data      #################################################ESTABLISHING VARIABLES
59  savepath = "data/example09a"      # where to save output data
60  meshpath = "data/example09m"      savepath = "data/example09a"
61  mkDir(savepath)      meshpath = "data/example09m"
62  #Geometric and material property related variables.      mkDir(savepath)
63  step=10.0 # the element size      #Geometric and material property related variables.
64        step=10.0 # the element size
65  vel2=1800.;   vel1=3000.  
66  rho2=2300.;   rho1=3100. #density      vel2=1800.;   vel1=3000.
67  mu2=vel2**2.*rho2/4.;  mu1=vel1**2.*rho1/4.  #bulk modulus      rho2=2300.;   rho1=3100. #density
68  lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2.  #lames constant      mu2=vel2**2.*rho2/4.;  mu1=vel1**2.*rho1/4.  #bulk modulus
69        lam2=vel2**2.*rho2/2.; lam1=vel1**2.*rho1/2.  #lames constant
70  ####################################################TESTING SWITCH  
71  testing=True      ####################################################TESTING SWITCH
72  if testing:      testing=True
73      print('The testing end time is currently selected. This severely limits the number of time iterations.')      if testing:
74      print("Try changing testing to False for more iterations.")          print('The testing end time is currently selected. This severely limits the number of time iterations.')
75      tend=0.001          print("Try changing testing to False for more iterations.")
76      #Model Parameters              tend=0.001
77      mx=40.          #Model Parameters    
78      my=40.          mx=40.
79      mz=20.          my=40.
80      outputs=5          mz=20.
81  else:          outputs=5
82      tend=0.1    # end time      else:
83      #Model Parameters          tend=0.1    # end time
84      mx=100.0   #x width of model          #Model Parameters
85      my=100.0   #y width of model          mx=100.0   #x width of model
86      mz=50.0   #depth of model          my=100.0   #y width of model
87      outputs=200          mz=50.0   #depth of model
88            outputs=200
89  ####################################################TIME RELATED VARIABLES  
90  h=0.00005    # time step      ####################################################TIME RELATED VARIABLES
91  # data recording times      h=0.00005    # time step
92  rtime=0.0 # first time to record      # data recording times
93  rtime_inc=tend/outputs # time increment to record      rtime=0.0 # first time to record
94  #Check to make sure number of time steps is not too large.      rtime_inc=tend/outputs # time increment to record
95  print("Time step size= ",h, "Expected number of outputs= ",tend/h)      #Check to make sure number of time steps is not too large.
96        print("Time step size= ",h, "Expected number of outputs= ",tend/h)
97  ####################################################CREATING THE SOURCE FUNCTION  
98  U0=0.1 # amplitude of point source      ####################################################CREATING THE SOURCE FUNCTION
99  ls=500   # length of the source      U0=0.1 # amplitude of point source
100  source=np.zeros(ls,'float') # source array      ls=500   # length of the source
101  decay1=np.zeros(ls,'float') # decay curve one      source=np.zeros(ls,'float') # source array
102  decay2=np.zeros(ls,'float') # decay curve two      decay1=np.zeros(ls,'float') # decay curve one
103  time=np.zeros(ls,'float')   # time values      decay2=np.zeros(ls,'float') # decay curve two
104  g=np.log(0.01)/ls      time=np.zeros(ls,'float')   # time values
105        g=np.log(0.01)/ls
106  dfeq=50 #Dominant Frequency  
107  a = 2.0 * (np.pi * dfeq)**2.0      dfeq=50 #Dominant Frequency
108  t0 = 5.0 / (2.0 * np.pi * dfeq)      a = 2.0 * (np.pi * dfeq)**2.0
109  srclength = 5. * t0      t0 = 5.0 / (2.0 * np.pi * dfeq)
110  ls = int(srclength/h)      srclength = 5. * t0
111  print('source length',ls)      ls = int(srclength/h)
112  source=np.zeros(ls,'float') # source array      print('source length',ls)
113  ampmax=0      source=np.zeros(ls,'float') # source array
114  for it in range(0,ls):      ampmax=0
115      t = it*h      for it in range(0,ls):
116      tt = t-t0          t = it*h
117      dum1 = np.exp(-a * tt * tt)          tt = t-t0
118      source[it] = -2. * a * tt * dum1          dum1 = np.exp(-a * tt * tt)
119      if (abs(source[it]) > ampmax):          source[it] = -2. * a * tt * dum1
120          ampmax = abs(source[it])          if (abs(source[it]) > ampmax):
121      time[t]=t*h              ampmax = abs(source[it])
122            time[t]=t*h
123  # will introduce a spherical source at middle left of bottom face  
124  xc=[mx/2,my/2,0]      # will introduce a spherical source at middle left of bottom face
125        xc=[mx/2,my/2,0]
126  ####################################################DOMAIN CONSTRUCTION  
127  domain=ReadMesh(os.path.join(meshpath,'example09m.fly')) # create the domain      ####################################################DOMAIN CONSTRUCTION
128  x=domain.getX() # get the locations of the nodes in the domain      domain=ReadMesh(os.path.join(meshpath,'example09m.fly')) # create the domain
129        x=domain.getX() # get the locations of the nodes in the domain
130  lam=Scalar(0,Function(domain))  
131  lam.setTaggedValue("vintfa",lam1)      lam=Scalar(0,Function(domain))
132  lam.setTaggedValue("vintfb",lam2)      lam.setTaggedValue("vintfa",lam1)
133  mu=Scalar(0,Function(domain))      lam.setTaggedValue("vintfb",lam2)
134  mu.setTaggedValue("vintfa",mu1)      mu=Scalar(0,Function(domain))
135  mu.setTaggedValue("vintfb",mu2)      mu.setTaggedValue("vintfa",mu1)
136  rho=Scalar(0,Function(domain))      mu.setTaggedValue("vintfb",mu2)
137  rho.setTaggedValue("vintfa",rho1)      rho=Scalar(0,Function(domain))
138  rho.setTaggedValue("vintfb",rho2)      rho.setTaggedValue("vintfa",rho1)
139        rho.setTaggedValue("vintfb",rho2)
140  ##########################################################ESTABLISH PDE  
141  mypde=LinearPDE(domain) # create pde      ##########################################################ESTABLISH PDE
142  mypde.setSymmetryOn() # turn symmetry on      mypde=LinearPDE(domain) # create pde
143  # turn lumping on for more efficient solving      mypde.setSymmetryOn() # turn symmetry on
144  #mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)      # turn lumping on for more efficient solving
145  kmat = kronecker(domain) # create the kronecker delta function of the domain      #mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
146  mypde.setValue(D=rho*kmat) #set the general form value D      kmat = kronecker(domain) # create the kronecker delta function of the domain
147        mypde.setValue(D=rho*kmat) #set the general form value D
148  ############################################FIRST TIME STEPS AND SOURCE  
149  # define small radius around point xc      ############################################FIRST TIME STEPS AND SOURCE
150  src_rad = 20; print("sourc radius= ",src_rad)      # define small radius around point xc
151  # set initial values for first two time steps with source terms      src_rad = 20; print("sourc radius= ",src_rad)
152  xb=FunctionOnBoundary(domain).getX()      # set initial values for first two time steps with source terms
153  yx=(cos(length(xb-xc)*3.1415/src_rad)+1)*whereNegative(length(xb-xc)-src_rad)      xb=FunctionOnBoundary(domain).getX()
154  stop=Scalar(0.0,FunctionOnBoundary(domain))      yx=(cos(length(xb-xc)*3.1415/src_rad)+1)*whereNegative(length(xb-xc)-src_rad)
155  stop.setTaggedValue("stop",1.0)      stop=Scalar(0.0,FunctionOnBoundary(domain))
156  src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down      stop.setTaggedValue("stop",1.0)
157        src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
158  mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary  
159  # initial value of displacement at point source is constant (U0=0.01)      mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
160  # for first two time steps      # initial value of displacement at point source is constant (U0=0.01)
161  u=[0.0,0.0,0.0]*wherePositive(x)      # for first two time steps
162  u_m1=u      u=[0.0,0.0,0.0]*wherePositive(x)
163        u_m1=u
164  ####################################################ITERATION VARIABLES  
165  n=0 # iteration counter      ####################################################ITERATION VARIABLES
166  t=0 # time counter      n=0 # iteration counter
167  ##############################################################ITERATION      t=0 # time counter
168  while t<tend:      ##############################################################ITERATION
169      # get current stress      while t<tend:
170      g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc          # get current stress
171      mypde.setValue(X=-stress) # set PDE values          g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))#*abc
172      accel = mypde.getSolution() #get PDE solution for accelleration          mypde.setValue(X=-stress) # set PDE values
173      u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement          accel = mypde.getSolution() #get PDE solution for accelleration
174      u_p1=u_p1#*abc          # apply boundary conditions          u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
175      u_m1=u; u=u_p1 # shift values by 1          u_p1=u_p1#*abc          # apply boundary conditions
176      # save current displacement, acceleration and pressure          u_m1=u; u=u_p1 # shift values by 1
177      if (t >= rtime):          # save current displacement, acceleration and pressure
178          saveVTK(os.path.join(savepath,"ex09a.%05d.vtu"%n),displacement=length(u),\          if (t >= rtime):
179                                      acceleration=length(accel),tensor=stress)              saveVTK(os.path.join(savepath,"ex09a.%05d.vtu"%n),displacement=length(u),\
180          rtime=rtime+rtime_inc #increment data save time                                          acceleration=length(accel),tensor=stress)
181      # increment loop values              rtime=rtime+rtime_inc #increment data save time
182      t=t+h; n=n+1          # increment loop values
183      if (n < ls):          t=t+h; n=n+1
184          mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary          if (n < ls):
185      print("time step %d, t=%s"%(n,t))              mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
186            print("time step %d, t=%s"%(n,t))

Legend:
Removed from v.5287  
changed lines
  Added in v.5288

  ViewVC Help
Powered by ViewVC 1.1.26