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

Diff of /trunk/doc/examples/cookbook/example09c.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  from esys.weipa import saveVTK  from esys.weipa import saveVTK
35  import os  import os
36  # smoothing operator  # smoothing operator
# Line 43  import numpy as np Line 41  import numpy as np
41  import pylab as pl  import pylab as pl
42  import matplotlib.cm as cm  import matplotlib.cm as cm
43  from esys.escript.linearPDEs import LinearPDE, SolverOptions  from esys.escript.linearPDEs import LinearPDE, SolverOptions
44  from esys.finley import ReadMesh  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/example09c"      # where to save output data
60  meshpath = "data/example09n"      savepath = "data/example09c"
61  mkDir(savepath)      meshpath = "data/example09n"
62  #Geometric and material property related variables.      mkDir(savepath)
63  domain=ReadMesh(os.path.join(savepath,'example09n.fly')) # create the domain      #Geometric and material property related variables.
64  x=Solution(domain).getX()      domain=ReadMesh(os.path.join(savepath,'example09n.fly')) # create the domain
65  #parameters layers 1,2,3,4 and fault      x=Solution(domain).getX()
66  prho=np.array([2200.,2500.,3200.,4500.,5500.]) #density      #parameters layers 1,2,3,4 and fault
67  pvel=np.array([1500.,2200.,3000.,3200.,5000.]) #velocity      prho=np.array([2200.,2500.,3200.,4500.,5500.]) #density
68  pmu=pvel**2.*prho/4.                              #bulk modulus      pvel=np.array([1500.,2200.,3000.,3200.,5000.]) #velocity
69  plam=pvel**2.*prho/2.                             #lames constant      pmu=pvel**2.*prho/4.                              #bulk modulus
70  nlayers=4      plam=pvel**2.*prho/2.                             #lames constant
71  width=300.0      nlayers=4
72  rho=Scalar(0,Function(domain))      width=300.0
73  vel=Scalar(0,Function(domain))      rho=Scalar(0,Function(domain))
74  mu=Scalar(0,Function(domain))      vel=Scalar(0,Function(domain))
75  lam=Scalar(0,Function(domain))      mu=Scalar(0,Function(domain))
76        lam=Scalar(0,Function(domain))
77  print(0.5*np.sqrt(prho/(plam+2*pmu))*0.5)  
78        print(0.5*np.sqrt(prho/(plam+2*pmu))*0.5)
79  for i in range(0,nlayers):  
80      rho.setTaggedValue('lblock%d'%i,prho[i])      for i in range(0,nlayers):
81      rho.setTaggedValue('rblock%d'%i,prho[i])          rho.setTaggedValue('lblock%d'%i,prho[i])
82      vel.setTaggedValue('lblock%d'%i,pvel[i])          rho.setTaggedValue('rblock%d'%i,prho[i])
83      vel.setTaggedValue('rblock%d'%i,pvel[i])          vel.setTaggedValue('lblock%d'%i,pvel[i])
84      mu.setTaggedValue('lblock%d'%i,pmu[i])          vel.setTaggedValue('rblock%d'%i,pvel[i])
85      mu.setTaggedValue('rblock%d'%i,pmu[i])          mu.setTaggedValue('lblock%d'%i,pmu[i])
86      lam.setTaggedValue('lblock%d'%i,plam[i])          mu.setTaggedValue('rblock%d'%i,pmu[i])
87      lam.setTaggedValue('rblock%d'%i,plam[i])          lam.setTaggedValue('lblock%d'%i,plam[i])
88  i=nlayers          lam.setTaggedValue('rblock%d'%i,plam[i])
89  rho.setTaggedValue('fault',prho[i])      i=nlayers
90  vel.setTaggedValue('fault',pvel[i])      rho.setTaggedValue('fault',prho[i])
91  mu.setTaggedValue('fault',pmu[i])      vel.setTaggedValue('fault',pvel[i])
92  lam.setTaggedValue('fault',plam[i])      mu.setTaggedValue('fault',pmu[i])
93        lam.setTaggedValue('fault',plam[i])
94    
95  # Time related variables.  
96  testing=False      # Time related variables.
97  if testing:      testing=False
98      print('The testing end time is currently selected. This severely limits the number of time iterations.')      if testing:
99      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.')
100      tend=0.1          print("Try changing testing to False for more iterations.")
101  else:          tend=0.1
102      tend=0.1    # end time      else:
103            tend=0.1    # end time
104  h=0.00001    # time step  
105  # data recording times      h=0.00001    # time step
106  rtime=0.0 # first time to record      # data recording times
107  rtime_inc=tend/750.0 # time increment to record      rtime=0.0 # first time to record
108  #Check to make sure number of time steps is not too large.      rtime_inc=tend/750.0 # time increment to record
109  print("Time step size= ",h, "Expected number of outputs= ",tend/h)      #Check to make sure number of time steps is not too large.
110        print("Time step size= ",h, "Expected number of outputs= ",tend/h)
111  U0=0.1 # amplitude of point source  
112  ls=500   # length of the source      U0=0.1 # amplitude of point source
113  source=np.zeros(ls,'float') # source array      ls=500   # length of the source
114  decay1=np.zeros(ls,'float') # decay curve one      source=np.zeros(ls,'float') # source array
115  decay2=np.zeros(ls,'float') # decay curve two      decay1=np.zeros(ls,'float') # decay curve one
116  time=np.zeros(ls,'float')   # time values      decay2=np.zeros(ls,'float') # decay curve two
117  g=np.log(0.01)/ls      time=np.zeros(ls,'float')   # time values
118        g=np.log(0.01)/ls
119  dfeq=50 #Dominant Frequency  
120  a = 2.0 * (np.pi * dfeq)**2.0      dfeq=50 #Dominant Frequency
121  t0 = 5.0 / (2.0 * np.pi * dfeq)      a = 2.0 * (np.pi * dfeq)**2.0
122  srclength = 5. * t0      t0 = 5.0 / (2.0 * np.pi * dfeq)
123  ls = int(srclength/h)      srclength = 5. * t0
124  print('source length',ls)      ls = int(srclength/h)
125  source=np.zeros(ls,'float') # source array      print('source length',ls)
126  ampmax=0      source=np.zeros(ls,'float') # source array
127  for it in range(0,ls):      ampmax=0
128      t = it*h      for it in range(0,ls):
129      tt = t-t0          t = it*h
130      dum1 = np.exp(-a * tt * tt)          tt = t-t0
131      source[it] = -2. * a * tt * dum1          dum1 = np.exp(-a * tt * tt)
132      if (abs(source[it]) > ampmax):          source[it] = -2. * a * tt * dum1
133          ampmax = abs(source[it])          if (abs(source[it]) > ampmax):
134      time[t]=t*h              ampmax = abs(source[it])
135            time[t]=t*h
136  # will introduce a spherical source at middle left of bottom face  
137  xc=[150,0]      # will introduce a spherical source at middle left of bottom face
138        xc=[150,0]
139  ##########################################################ESTABLISH PDE  
140  mypde=LinearPDE(domain) # create pde      ##########################################################ESTABLISH PDE
141  mypde.setSymmetryOn() # turn symmetry on      mypde=LinearPDE(domain) # create pde
142  # turn lumping on for more efficient solving      mypde.setSymmetryOn() # turn symmetry on
143  mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)      # turn lumping on for more efficient solving
144  kmat = kronecker(domain) # create the kronecker delta function of the domain      mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
145  mypde.setValue(D=rho*kmat) #set the general form value D      kmat = kronecker(domain) # create the kronecker delta function of the domain
146        mypde.setValue(D=rho*kmat) #set the general form value D
147  ############################################FIRST TIME STEPS AND SOURCE  
148  # define small radius around point xc      ############################################FIRST TIME STEPS AND SOURCE
149  src_length = 10; print("src_length = ",src_length)      # define small radius around point xc
150  # set initial values for first two time steps with source terms      src_length = 10; print("src_length = ",src_length)
151  xb=FunctionOnBoundary(domain).getX()      # set initial values for first two time steps with source terms
152  yx=(cos(length(xb-xc)*3.1415/src_length)+1)*whereNegative(length(xb-xc)-src_length)      xb=FunctionOnBoundary(domain).getX()
153  stop=Scalar(0.0,FunctionOnBoundary(domain))      yx=(cos(length(xb-xc)*3.1415/src_length)+1)*whereNegative(length(xb-xc)-src_length)
154  stop.setTaggedValue("top",1.0)      stop=Scalar(0.0,FunctionOnBoundary(domain))
155  src_dir=numpy.array([0.,-1.]) # defines direction of point source as down      stop.setTaggedValue("top",1.0)
156        src_dir=numpy.array([0.,-1.]) # defines direction of point source as down
157  mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary  
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)  
160  # for first two time steps      # initial value of displacement at point source is constant (U0=0.01)
161  u=[0.0,0.0]*x      # for first two time steps
162  u_m1=u      u=[0.0,0.0]*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,"ex09c.%05d.vtu"%n),displacement=length(u),\          if (t >= rtime):
179                                      acceleration=length(accel),tensor=stress)              saveVTK(os.path.join(savepath,"ex09c.%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[n]*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[n]*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