 # Diff of /trunk/doc/cookbook/example08.tex

revision 3152 by ahallam, Thu Jul 15 02:57:46 2010 UTC revision 3153 by ahallam, Sun Sep 5 01:31:49 2010 UTC
# Line 14  Line 14
14  \section{Seismic Wave Propagation in Two Dimensions}  \section{Seismic Wave Propagation in Two Dimensions}
15
16  \sslist{example08a.py}  \sslist{example08a.py}

17  We will now expand upon the previous chapter by introducing a vector form of  We will now expand upon the previous chapter by introducing a vector form of
18  the wave equation. This means that the waves will have both a scalar magnitude,  the wave equation. This means that the waves will have not only a scalar
19  but also a direction. This type of scenario is apparent in wave forms that  magnitude as for the pressure wave solution, but also a direction. This type of
20  exhibit compressional and transverse particle motion. A common type of wave  scenario is apparent in wave types that exhibit compressional and transverse
21  that obeys this principle are seismic waves.  particle motion. An example of this would be seismic waves.
22
23  Wave propagation in the earth can be described by the elastic wave equation  Wave propagation in the earth can be described by the elastic wave equation
24  \begin{equation} \label{eqn:wav} \index{wave equation}  \begin{equation} \label{eqn:wav} \index{wave equation}
# Line 31  where $\sigma$ is the stress given by Line 30  where $\sigma$ is the stress given by
30   \sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu (   \sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu (
31  u\hackscore{i,j} + u\hackscore{j,i})  u\hackscore{i,j} + u\hackscore{j,i})
32  \end{equation}  \end{equation}
33  with $\lambda$ and $\mu$ representing the Lame Coefficients. Specifically  and $\lambda$ and $\mu$ represent Lame's parameters. Specifically for seismic
34  $\mu$ is the bulk modulus.  waves, $\mu$ is the propagation materials shear modulus.
35  In a similar process to the previous chapter, we will use the acceleration  In a similar process to the previous chapter, we will use the acceleration
36  solution to solve this PDE. By substituting $a$ directly for  solution to solve this PDE. By substituting $a$ directly for
37  $\frac{\partial^{2}u\hackscore{i}}{\partial t^2}$ we can derive the  $\frac{\partial^{2}u\hackscore{i}}{\partial t^2}$ we can derive the
# Line 41  displacement solution. Using $a$ we can Line 40  displacement solution. Using $a$ we can
40  \rho a\hackscore{i} - \frac{\partial  \rho a\hackscore{i} - \frac{\partial
41  \sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0  \sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0
42  \end{equation}  \end{equation}
43  This is very similar to our acceleration solution in the previous chapter. Thus  Thus the problem will be solved for acceleration and then converted to
44  the problem will be solved for acceleration and for displacement using the  displacement using the backwards difference approximation.
backwards difference approximation.
45
46  Consider now the stress $\sigma$. One can see that the stress consists of two  Consider now the stress $\sigma$. One can see that the stress consists of two
47  distinct terms:  distinct terms:
# Line 87  rtime=0.0 # first time to record Line 85  rtime=0.0 # first time to record
85  rtime_inc=tend/20.0 # time increment to record  rtime_inc=tend/20.0 # time increment to record
86  \end{python}  \end{python}
87  Currently the PDE solution will be saved to file $20$ times between the start of  Currently the PDE solution will be saved to file $20$ times between the start of
88  the modelling and the final time step. With these parameters set an if catch is  the modelling and the final time step. With these parameters set, an if
89  introduced to the time loop  statement is introduced to the time loop
90  \begin{python}  \begin{python}
91  if (t >= rtime):  if (t >= rtime):
92      saveVTK(os.path.join(savepath,"ex08a.%05d.vtu"%n),displacement=length(u),\      saveVTK(os.path.join(savepath,"ex08a.%05d.vtu"%n),displacement=length(u),\
# Line 111  $escript -t 4 example08a.py Line 109$escript -t 4 example08a.py
109  \end{verbatim}  \end{verbatim}
110  would call the script in this section and solve it using 4 threads.  would call the script in this section and solve it using 4 threads.
111
112    The computation times on an increasing number of cores is outlines in Table
113    \ref{tab:wpcores}
114
115    \begin{table}[ht]
116    \begin{center}
117    \caption{Computation times for an increasing number of cores.}
118    \label{tab:wpcores}
119    \begin{tabular}{| c | c |}
120    \hline
121    Number of Cores & Time (s) \\
122    \hline
123    1 & 691.0 \\
124    2 & 400.0 \\
125    3 & 305.0 \\
126    4 & 328.0 \\
127    5 & 323.0 \\
128    6 & 292.0 \\
129    7 & 282.0 \\
130    8 & 445.0 \\ \hline
131    \end{tabular}
132    \end{center}
133    \end{table}
134
135  \section{Vector source on the boundary}  \section{Vector source on the boundary}
136    \sslist{example08b.py}
137  For this particular example, we will introduce the source by applying a  For this particular example, we will introduce the source by applying a
138  displacment to the boundary during the initial time steps. The source will again  displacment to the boundary during the initial time steps. The source will again
139  be  be
# Line 127  xc=[ndx/2,0] Line 149  xc=[ndx/2,0]
149
150  ############################################FIRST TIME STEPS AND SOURCE  ############################################FIRST TIME STEPS AND SOURCE
151  # define small radius around point xc  # define small radius around point xc
152  src_length = 20; print "src_length = ",src_length  src_length = 40; print "src_length = ",src_length
153  # set initial values for first two time steps with source terms  # set initial values for first two time steps with source terms
154  y=U0*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)\  xb=FunctionOnBoundary(domain).getX()
155                                      -src_length)  y=source*(cos(length(x-xc)*3.1415/src_length)+1)*\
156                    whereNegative(length(xb-src_length))
157    src_dir=numpy.array([0.,1.]) # defines direction of point source as down
158    y=y*src_dir
159  \end{python}  \end{python}
160  where \verb xc  is the source point on the boundary of the model. The source  where \verb xc  is the source point on the boundary of the model. Note that
161  direction is then defined as an $(x,y)$ array and multiplied by the source  because the source is specifically located on the boundary, we have used the
162  function. The directional array must have a magnitude of $\left| 1 \right|$  \verb!FunctionOnBoundary! call to ensure the nodes are located upon the
163  otherwise the  boundary only. These boundary nodes are passed to source as \verb!xb!. The
164  amplitude of the source will become modified. For this example, the source is  source direction is then defined as an $(x,y)$ array and multiplied by the
165  directed in the $-y$ direction.  source function. The directional array must have a magnitude of $\left| 1 166 \right|$ otherwise the amplitude of the source will become modified. For this
167    example, the source is directed in the $-y$ direction.
168  \begin{python}  \begin{python}
169  src_dir=numpy.array([0.,-1.]) # defines direction of point source as down  src_dir=numpy.array([0.,-1.]) # defines direction of point source as down
170  y=y*src_dir  y=y*src_dir
# Line 153  must set the model state to zero for the Line 180  must set the model state to zero for the
180  \begin{python}  \begin{python}
181  # initial value of displacement at point source is constant (U0=0.01)  # initial value of displacement at point source is constant (U0=0.01)
182  # for first two time steps  # for first two time steps
183  u=[0.0,0.0]*whereNegative(x)  u=[0.0,0.0]*wherePositive(x)
184  u_m1=u  u_m1=u
185  \end{python}  \end{python}
186
# Line 181  iteration stage. This is covered in the Line 208  iteration stage. This is covered in the
208  \label{fig:ex08pw}  \label{fig:ex08pw}
209  \caption{Results of Example 08 at various times.}  \caption{Results of Example 08 at various times.}
210  \end{figure}  \end{figure}
211    \clearpage
212
213  \section{Time variant source}  \section{Time variant source}
214
215  \sslist{example08a.py}  \sslist{example08b.py}

216  Until this point, all of the wave propagation examples in this cookbook have  Until this point, all of the wave propagation examples in this cookbook have
217  used impulsive sources which are smooth in space but not time. It is however,  used impulsive sources which are smooth in space but not time. It is however,
218  advantageous to have a time smoothed source as it can reduce the temporal  advantageous to have a time smoothed source as it can reduce the temporal
219  frequency range and thus mitigate aliasing in the solution.  frequency range and thus mitigate aliasing in the solution.
220

221  It is quite  It is quite
222  simple to implement a source which is smooth in time. In addition to the  simple to implement a source which is smooth in time. In addition to the
223  original source function the only extra requirement is a time function. For  original source function the only extra requirement is a time function. For
224  this example it will be a decaying sinusoidal curve defined by;  this example the time variant source will be the derivative of a gausian curve
225    defined by the required dominant frequency (\refig{fig:tvsource}).
226  \begin{python}  \begin{python}
227  U0=0.1 # amplitude of point source  #Creating the time function of the source.
228  ls=100   # length of the source  dfeq=50 #Dominant Frequency
229    a = 2.0 * (np.pi * dfeq)**2.0
230    t0 = 5.0 / (2.0 * np.pi * dfeq)
231    srclength = 5. * t0
232    ls = int(srclength/h)
233    print 'source length',ls
234  source=np.zeros(ls,'float') # source array  source=np.zeros(ls,'float') # source array
235  g=np.log(0.01)/ls  ampmax=0
236  for t in range(0,ls):  for it in range(0,ls):
237      source[t]=np.exp(g*t)*U0*np.sin(2.*np.pi*t/(0.75*ls))      t = it*h
238        tt = t-t0
239        dum1 = np.exp(-a * tt * tt)
240        source[it] = -2. * a * tt * dum1
241        if (abs(source[it]) > ampmax):
242            ampmax = abs(source[it])
243        time[t]=t*h
244  \end{python}  \end{python}
245    \begin{figure}[ht]
246    \centering
247    \includegraphics[width=3in]{figures/source.png}
248    \caption{Time variant source with a dominant frequency of 50Hz.}
249    \label{fig:tvsource}
250    \end{figure}
251
252  We then build the source and the first two time steps via;  We then build the source and the first two time steps via;
253  \begin{python}  \begin{python}
# Line 247  elastic wave equations}, 1985, Cerjan C, Line 291  elastic wave equations}, 1985, Cerjan C,
291  where the solution and the stress are multiplied by a damping function defined  where the solution and the stress are multiplied by a damping function defined
292  on $n$ nodes of the domain adjacent to the boundary, given by;  on $n$ nodes of the domain adjacent to the boundary, given by;
293  \begin{equation}  \begin{equation}
294   y=  \gamma =\sqrt\left(\frac{|-log(\gamma\hackscore{b})|}{n^2}\right)
295    \end{equation}
296    \begin{equation}
297    y=e^{-(\gamma x)^2}
298  \end{equation}  \end{equation}
299  This is applied to the bounding 20-50 pts of the model using the location  This is applied to the bounding 20-50 pts of the model using the location
300  specifiers of \esc;  specifiers of \esc;
# Line 291  abcbottom=abcbottom+whereZero(abcbottom) Line 338  abcbottom=abcbottom+whereZero(abcbottom)
338  abc=abcleft*abcright*abcbottom  abc=abcleft*abcright*abcbottom
339  \end{python}  \end{python}
340  Note that the boundary conditions are not applied to the surface, as this is  Note that the boundary conditions are not applied to the surface, as this is
341  effectively a free surface where normal reflections would be experienced. The  effectively a free surface where normal reflections would be experienced.
342  resulting boundary damping function can be viewed in Figure \ref{fig:abconds}  Special conditions can be introduced at this surface if they are known. The
343    resulting boundary damping function can be viewed in Figure
344    \ref{fig:abconds}.
345
346    \section{Second order Meshing}
347    For stiff problems like the wave equation it is often prudent to implement
348    second order meshing. This creates a more accurate mesh approximation with some
349    increased processing cost. To turn second order meshing on, the \verb!rectangle!
350    function accpets an \verb!order! keyword argument.
351    \begin{python}
352    domain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy,order=2) # create the domain
353    \end{python}
354    Other pycad functions and objects have similar keyword arguments for higher
355    order meshing.
356
357    Note that when implementing second order meshing, a smaller timestep is required
358    then for first order meshes as the second order essentially reduces the size of
359    the mesh by half.
360
361  \begin{figure}[ht]  \begin{figure}[ht]
362   \centering   \centering
# Line 322  resulting boundary damping function can Line 386  resulting boundary damping function can
386  \label{fig:ex08pw}  \label{fig:ex08pw}
387  \caption{Results of Example 08b at various times.}  \caption{Results of Example 08b at various times.}
388  \end{figure}  \end{figure}
389    \clearpage
390
392    \sslist{example08c.py}
393    To make the problem more interesting we will now introduce an interface to the
394    middle of the domain. Infact we will use the same domain as we did for heat flux
395    in Chapter \ref{CHAP HEAT 2}. The domain contains a syncline with two set of
396    material properties on either side of the interface.
397
398    \begin{figure}[ht]
399    \begin{center}
400    \includegraphics[width=5in]{figures/example08cgeo.png}
401    \caption{Domain geometry for example08c.py showing line tangents.}
402    \label{fig:ex08cgeo}
403    \end{center}
404    \end{figure}
405
406    It is simple enough to slightly modify the scripts of the previous sections to
407    accept this domain. Multiple material parameters must now be deined and assigned
408    to specific tagged areas. Again this is done via
409    \begin{python}
410    lam=Scalar(0,Function(domain))
411    lam.setTaggedValue("top",lam1)
412    lam.setTaggedValue("bottom",lam2)
413    mu=Scalar(0,Function(domain))
414    mu.setTaggedValue("top",mu1)
415    mu.setTaggedValue("bottom",mu2)
416    rho=Scalar(0,Function(domain))
417    rho.setTaggedValue("top",rho1)
418    rho.setTaggedValue("bottom",rho2)
419    \end{python}
420    Don't forget that teh source boudnary must also be tagged and added so it can be
421    referenced
422    \begin{python}
423    # Add the subdomains and flux boundaries.
425                                         PropertySet("linetop",l30))
426    \end{python}
427    It is now possible to solve the script as in the previous examples.
428
429    \begin{figure}[ht]
430    \centering
431    \includegraphics[width=4in]{figures/ex08c2601.png}
432    \caption{Modelling results of example08c.py at 0.2601 seconds. Notice the
433    refraction of the wave front about the boundary between the two layers.}
434    \label{fig:ex08cres}
435    \end{figure}
436
437

Legend:
 Removed from v.3152 changed lines Added in v.3153