| 1 |
|
| 2 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 3 |
% |
| 4 |
% Copyright (c) 2003-2010 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 |
\section{Seismic Wave Propagation in Two Dimensions} |
| 15 |
|
| 16 |
\sslist{example08a.py} |
| 17 |
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 not only a scalar |
| 19 |
magnitude as for the pressure wave solution, but also a direction. This type of |
| 20 |
scenario is apparent in wave types that exhibit compressional and transverse |
| 21 |
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 |
| 24 |
\begin{equation} \label{eqn:wav} \index{wave equation} |
| 25 |
\rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2} - \frac{\partial |
| 26 |
\sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0 |
| 27 |
\end{equation} |
| 28 |
where $\sigma$ is the stress given by |
| 29 |
\begin{equation} \label{eqn:sigw} |
| 30 |
\sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( |
| 31 |
u\hackscore{i,j} + u\hackscore{j,i}) |
| 32 |
\end{equation} |
| 33 |
and $\lambda$ and $\mu$ represent Lame's parameters. Specifically for seismic |
| 34 |
waves, $\mu$ is the propagation materials shear modulus. |
| 35 |
In a similar process to the previous chapter, we will use the acceleration |
| 36 |
solution to solve this PDE. By substituting $a$ directly for |
| 37 |
$\frac{\partial^{2}u\hackscore{i}}{\partial t^2}$ we can derive the |
| 38 |
displacement solution. Using $a$ we can see that \autoref{eqn:wav} becomes |
| 39 |
\begin{equation} \label{eqn:wava} |
| 40 |
\rho a\hackscore{i} - \frac{\partial |
| 41 |
\sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0 |
| 42 |
\end{equation} |
| 43 |
Thus the problem will be solved for acceleration and then converted to |
| 44 |
displacement using the backwards difference approximation. |
| 45 |
|
| 46 |
Consider now the stress $\sigma$. One can see that the stress consists of two |
| 47 |
distinct terms: |
| 48 |
\begin{subequations} |
| 49 |
\begin{equation} \label{eqn:sigtrace} |
| 50 |
\lambda u\hackscore{k,k} \delta\hackscore{ij} |
| 51 |
\end{equation} |
| 52 |
\begin{equation} \label{eqn:sigtrans} |
| 53 |
\mu (u\hackscore{i,j} + u\hackscore{j,i}) |
| 54 |
\end{equation} |
| 55 |
\end{subequations} |
| 56 |
One simply recognizes in \autoref{eqn:sigtrace} that $u\hackscore{k,k}$ is the |
| 57 |
trace of the displacement solution and that $\delta\hackscore{ij}$ is the |
| 58 |
kronecker delta function with dimensions equivalent to $u$. The second term |
| 59 |
\autoref{eqn:sigtrans} is the sum of $u$ with its own transpose. Putting these |
| 60 |
facts together we see that the spatial differential of the stress is given by the |
| 61 |
gradient of $u$ and the aforementioned opperations. This value is then submitted |
| 62 |
to the \esc PDE as $X$. |
| 63 |
\begin{python} |
| 64 |
g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g)) |
| 65 |
mypde.setValue(X=-stress) # set PDE values |
| 66 |
\end{python} |
| 67 |
The solution is then obtained via the usual method and the displacement is |
| 68 |
calculated so that the memory variables can be updated for the next time |
| 69 |
iteration. |
| 70 |
\begin{python} |
| 71 |
accel = mypde.getSolution() #get PDE solution for accelleration |
| 72 |
u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement |
| 73 |
u_m1=u; u=u_p1 # shift values by 1 |
| 74 |
\end{python} |
| 75 |
|
| 76 |
Saving the data has been handled slightly differently in this example. The VTK |
| 77 |
files generated can be quite large and take a significant amount of time to save |
| 78 |
to the hard disk. To avoid doing this at every iteration a test is devised which |
| 79 |
saves only at specific time intervals. |
| 80 |
|
| 81 |
To do this there are two new parameters in our script. |
| 82 |
\begin{python} |
| 83 |
# data recording times |
| 84 |
rtime=0.0 # first time to record |
| 85 |
rtime_inc=tend/20.0 # time increment to record |
| 86 |
\end{python} |
| 87 |
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 |
| 89 |
statement is introduced to the time loop |
| 90 |
\begin{python} |
| 91 |
if (t >= rtime): |
| 92 |
saveVTK(os.path.join(savepath,"ex08a.%05d.vtu"%n),displacement=length(u),\ |
| 93 |
acceleration=length(accel),tensor=stress) |
| 94 |
rtime=rtime+rtime_inc #increment data save time |
| 95 |
\end{python} |
| 96 |
\verb!t! is the time counter. Whenever the recording time \verb!rtime! is less |
| 97 |
then \verb!t! the solution is saved and \verb!rtime! is incremented. This |
| 98 |
limits the number of outputs and increases the speed of the solver. |
| 99 |
|
| 100 |
\section{Multi-threading} |
| 101 |
The wave equation solution can be quite demanding on cpu time. Enhancements can |
| 102 |
be made by accessing multiple threads or cores on your computer. This does not |
| 103 |
require any modification to the solution script and only comes into play when |
| 104 |
escript is called from the shell. To use multiple threads \esc is called using |
| 105 |
the \verb!-t! option with an interger argument for the number of threads |
| 106 |
required. For example |
| 107 |
\begin{verbatim} |
| 108 |
$escript -t 4 example08a.py |
| 109 |
\end{verbatim} |
| 110 |
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 |
| 113 |
\autoref{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} |
| 136 |
\sslist{example08b.py} |
| 137 |
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 |
| 139 |
be |
| 140 |
a radially propagating wave but due to the vector nature of the PDE used, a |
| 141 |
direction will need to be applied to the source. |
| 142 |
|
| 143 |
The first step is to choose an amplitude and create the source as in the |
| 144 |
previous chapter. |
| 145 |
\begin{python} |
| 146 |
U0=0.01 # amplitude of point source |
| 147 |
# will introduce a spherical source at middle left of bottom face |
| 148 |
xc=[ndx/2,0] |
| 149 |
|
| 150 |
############################################FIRST TIME STEPS AND SOURCE |
| 151 |
# define small radius around point xc |
| 152 |
src_length = 40; print "src_length = ",src_length |
| 153 |
# set initial values for first two time steps with source terms |
| 154 |
xb=FunctionOnBoundary(domain).getX() |
| 155 |
y=source[0]*(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} |
| 160 |
where \verb xc is the source point on the boundary of the model. Note that |
| 161 |
because the source is specifically located on the boundary, we have used the |
| 162 |
\verb!FunctionOnBoundary! call to ensure the nodes are located upon the |
| 163 |
boundary only. These boundary nodes are passed to source as \verb!xb!. The |
| 164 |
source direction is then defined as an $(x,y)$ array and multiplied by the |
| 165 |
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} |
| 169 |
src_dir=numpy.array([0.,-1.]) # defines direction of point source as down |
| 170 |
y=y*src_dir |
| 171 |
\end{python} |
| 172 |
The function can then be applied as a boundary condition by setting it equal to |
| 173 |
$y$ in the general form. |
| 174 |
\begin{python} |
| 175 |
mypde.setValue(y=y) #set the source as a function on the boundary |
| 176 |
\end{python} |
| 177 |
The final step is to qualify the initial conditions. Due to the fact that we are |
| 178 |
no longer using the source to define our initial condition to the model, we |
| 179 |
must set the model state to zero for the first two time steps. |
| 180 |
\begin{python} |
| 181 |
# initial value of displacement at point source is constant (U0=0.01) |
| 182 |
# for first two time steps |
| 183 |
u=[0.0,0.0]*wherePositive(x) |
| 184 |
u_m1=u |
| 185 |
\end{python} |
| 186 |
|
| 187 |
If the source is time progressive, $y$ can be updated during the |
| 188 |
iteration stage. This is covered in the following section. |
| 189 |
|
| 190 |
\begin{figure}[htp] |
| 191 |
\centering |
| 192 |
\subfigure[Example 08a at 0.025s ]{ |
| 193 |
\includegraphics[width=3in]{figures/ex08pw50.png} |
| 194 |
\label{fig:ex08pw50} |
| 195 |
} |
| 196 |
\subfigure[Example 08a at 0.175s ]{ |
| 197 |
\includegraphics[width=3in]{figures/ex08pw350.png} |
| 198 |
\label{fig:ex08pw350} |
| 199 |
} \\ |
| 200 |
\subfigure[Example 08a at 0.325s ]{ |
| 201 |
\includegraphics[width=3in]{figures/ex08pw650.png} |
| 202 |
\label{fig:ex08pw650} |
| 203 |
} |
| 204 |
\subfigure[Example 08a at 0.475s ]{ |
| 205 |
\includegraphics[width=3in]{figures/ex08pw950.png} |
| 206 |
\label{fig:ex08pw950} |
| 207 |
} |
| 208 |
\label{fig:ex08pw} |
| 209 |
\caption{Results of Example 08 at various times.} |
| 210 |
\end{figure} |
| 211 |
\clearpage |
| 212 |
|
| 213 |
\section{Time variant source} |
| 214 |
|
| 215 |
\sslist{example08b.py} |
| 216 |
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, |
| 218 |
advantageous to have a time smoothed source as it can reduce the temporal |
| 219 |
frequency range and thus mitigate aliasing in the solution. |
| 220 |
|
| 221 |
It is quite |
| 222 |
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 |
| 224 |
this example the time variant source will be the derivative of a gausian curve |
| 225 |
defined by the required dominant frequency (\autoref{fig:tvsource}). |
| 226 |
\begin{python} |
| 227 |
#Creating the time function of the source. |
| 228 |
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 |
| 235 |
ampmax=0 |
| 236 |
for it in range(0,ls): |
| 237 |
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} |
| 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; |
| 253 |
\begin{python} |
| 254 |
# set initial values for first two time steps with source terms |
| 255 |
y=source[0] |
| 256 |
*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length) |
| 257 |
src_dir=numpy.array([0.,-1.]) # defines direction of point source as down |
| 258 |
y=y*src_dir |
| 259 |
mypde.setValue(y=y) #set the source as a function on the boundary |
| 260 |
# initial value of displacement at point source is constant (U0=0.01) |
| 261 |
# for first two time steps |
| 262 |
u=[0.0,0.0]*whereNegative(x) |
| 263 |
u_m1=u |
| 264 |
\end{python} |
| 265 |
|
| 266 |
Finally, for the length of the source, we are required to update each new |
| 267 |
solution in the itterative section of the solver. This is done via; |
| 268 |
\begin{python} |
| 269 |
# increment loop values |
| 270 |
t=t+h; n=n+1 |
| 271 |
if (n < ls): |
| 272 |
y=source[n]**(cos(length(x-xc)*3.1415/src_length)+1)*\ |
| 273 |
whereNegative(length(x-xc)-src_length) |
| 274 |
y=y*src_dir; mypde.setValue(y=y) #set the source as a function on the |
| 275 |
boundary |
| 276 |
\end{python} |
| 277 |
|
| 278 |
\section{Absorbing Boundary Conditions} |
| 279 |
To mitigate the effect of the boundary on the model, absorbing boundary |
| 280 |
conditions can be introduced. These conditions effectively dampen the wave |
| 281 |
energy as they approach the bounday and thus prevent that energy from being |
| 282 |
reflected. This type of approach is used typically when a model only represents |
| 283 |
a small portion of the entire model, which in reality may have infinite bounds. |
| 284 |
It is inpractical to calculate the solution for an infinite model and thus ABCs |
| 285 |
allow us the create an approximate solution with small to zero boundary effects |
| 286 |
on a model with a solvable size. |
| 287 |
|
| 288 |
To dampen the waves, the method of \citet{Cerjan1985} |
| 289 |
where the solution and the stress are multiplied by a damping function defined |
| 290 |
on $n$ nodes of the domain adjacent to the boundary, given by; |
| 291 |
\begin{equation} |
| 292 |
\gamma =\sqrt{\frac{| -log( \gamma \hackscore{b} ) |}{n^2}} |
| 293 |
\end{equation} |
| 294 |
\begin{equation} |
| 295 |
y=e^{-(\gamma x)^2} |
| 296 |
\end{equation} |
| 297 |
This is applied to the bounding 20-50 pts of the model using the location |
| 298 |
specifiers of \esc; |
| 299 |
\begin{python} |
| 300 |
# Define where the boundary decay will be applied. |
| 301 |
bn=30. |
| 302 |
bleft=xstep*bn; bright=mx-(xstep*bn); bbot=my-(ystep*bn) |
| 303 |
# btop=ystep*bn # don't apply to force boundary!!! |
| 304 |
|
| 305 |
# locate these points in the domain |
| 306 |
left=x[0]-bleft; right=x[0]-bright; bottom=x[1]-bbot |
| 307 |
|
| 308 |
tgamma=0.98 # decay value for exponential function |
| 309 |
def calc_gamma(G,npts): |
| 310 |
func=np.sqrt(abs(-1.*np.log(G)/(npts**2.))) |
| 311 |
return func |
| 312 |
|
| 313 |
gleft = calc_gamma(tgamma,bleft) |
| 314 |
gright = calc_gamma(tgamma,bleft) |
| 315 |
gbottom= calc_gamma(tgamma,ystep*bn) |
| 316 |
|
| 317 |
print 'gamma', gleft,gright,gbottom |
| 318 |
|
| 319 |
# calculate decay functions |
| 320 |
def abc_bfunc(gamma,loc,x,G): |
| 321 |
func=exp(-1.*(gamma*abs(loc-x))**2.) |
| 322 |
return func |
| 323 |
|
| 324 |
fleft=abc_bfunc(gleft,bleft,x[0],tgamma) |
| 325 |
fright=abc_bfunc(gright,bright,x[0],tgamma) |
| 326 |
fbottom=abc_bfunc(gbottom,bbot,x[1],tgamma) |
| 327 |
# apply these functions only where relevant |
| 328 |
abcleft=fleft*whereNegative(left) |
| 329 |
abcright=fright*wherePositive(right) |
| 330 |
abcbottom=fbottom*wherePositive(bottom) |
| 331 |
# make sure the inside of the abc is value 1 |
| 332 |
abcleft=abcleft+whereZero(abcleft) |
| 333 |
abcright=abcright+whereZero(abcright) |
| 334 |
abcbottom=abcbottom+whereZero(abcbottom) |
| 335 |
# multiply the conditions together to get a smooth result |
| 336 |
abc=abcleft*abcright*abcbottom |
| 337 |
\end{python} |
| 338 |
Note that the boundary conditions are not applied to the surface, as this is |
| 339 |
effectively a free surface where normal reflections would be experienced. |
| 340 |
Special conditions can be introduced at this surface if they are known. The |
| 341 |
resulting boundary damping function can be viewed in |
| 342 |
\autoref{fig:abconds}. |
| 343 |
|
| 344 |
\section{Second order Meshing} |
| 345 |
For stiff problems like the wave equation it is often prudent to implement |
| 346 |
second order meshing. This creates a more accurate mesh approximation with some |
| 347 |
increased processing cost. To turn second order meshing on, the \verb!rectangle! |
| 348 |
function accpets an \verb!order! keyword argument. |
| 349 |
\begin{python} |
| 350 |
domain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy,order=2) # create the domain |
| 351 |
\end{python} |
| 352 |
Other pycad functions and objects have similar keyword arguments for higher |
| 353 |
order meshing. |
| 354 |
|
| 355 |
Note that when implementing second order meshing, a smaller timestep is required |
| 356 |
then for first order meshes as the second order essentially reduces the size of |
| 357 |
the mesh by half. |
| 358 |
|
| 359 |
\begin{figure}[ht] |
| 360 |
\centering |
| 361 |
\includegraphics[width=5in]{figures/ex08babc.png} |
| 362 |
\label{fig:abconds} |
| 363 |
\caption{Absorbing boundary conditions for example08b.py} |
| 364 |
\end{figure} |
| 365 |
|
| 366 |
\begin{figure}[htp] |
| 367 |
\centering |
| 368 |
\subfigure[Example 08b at 0.03s ]{ |
| 369 |
\includegraphics[width=3in]{figures/ex08sw060.png} |
| 370 |
\label{fig:ex08pw060} |
| 371 |
} |
| 372 |
\subfigure[Example 08b at 0.16s ]{ |
| 373 |
\includegraphics[width=3in]{figures/ex08sw320.png} |
| 374 |
\label{fig:ex08pw320} |
| 375 |
} \\ |
| 376 |
\subfigure[Example 08b at 0.33s ]{ |
| 377 |
\includegraphics[width=3in]{figures/ex08sw660.png} |
| 378 |
\label{fig:ex08pw660} |
| 379 |
} |
| 380 |
\subfigure[Example 08b at 0.44s ]{ |
| 381 |
\includegraphics[width=3in]{figures/ex08sw880.png} |
| 382 |
\label{fig:ex08pw880} |
| 383 |
} |
| 384 |
\label{fig:ex08pw} |
| 385 |
\caption{Results of Example 08b at various times.} |
| 386 |
\end{figure} |
| 387 |
\clearpage |
| 388 |
|
| 389 |
\section{Pycad example} |
| 390 |
\sslist{example08c.py} |
| 391 |
To make the problem more interesting we will now introduce an interface to the |
| 392 |
middle of the domain. Infact we will use the same domain as we did for heat flux |
| 393 |
in \autoref{CHAP HEAT 2}. The domain contains a syncline with two set of |
| 394 |
material properties on either side of the interface. |
| 395 |
|
| 396 |
\begin{figure}[ht] |
| 397 |
\begin{center} |
| 398 |
\includegraphics[width=5in]{figures/gmsh-example08c.png} |
| 399 |
\caption{Domain geometry for example08c.py showing line tangents.} |
| 400 |
\label{fig:ex08cgeo} |
| 401 |
\end{center} |
| 402 |
\end{figure} |
| 403 |
|
| 404 |
It is simple enough to slightly modify the scripts of the previous sections to |
| 405 |
accept this domain. Multiple material parameters must now be deined and assigned |
| 406 |
to specific tagged areas. Again this is done via |
| 407 |
\begin{python} |
| 408 |
lam=Scalar(0,Function(domain)) |
| 409 |
lam.setTaggedValue("top",lam1) |
| 410 |
lam.setTaggedValue("bottom",lam2) |
| 411 |
mu=Scalar(0,Function(domain)) |
| 412 |
mu.setTaggedValue("top",mu1) |
| 413 |
mu.setTaggedValue("bottom",mu2) |
| 414 |
rho=Scalar(0,Function(domain)) |
| 415 |
rho.setTaggedValue("top",rho1) |
| 416 |
rho.setTaggedValue("bottom",rho2) |
| 417 |
\end{python} |
| 418 |
Don't forget that teh source boudnary must also be tagged and added so it can be |
| 419 |
referenced |
| 420 |
\begin{python} |
| 421 |
# Add the subdomains and flux boundaries. |
| 422 |
d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\ |
| 423 |
PropertySet("linetop",l30)) |
| 424 |
\end{python} |
| 425 |
It is now possible to solve the script as in the previous examples. |
| 426 |
|
| 427 |
\begin{figure}[ht] |
| 428 |
\centering |
| 429 |
\includegraphics[width=4in]{figures/ex08c2601.png} |
| 430 |
\caption{Modelling results of example08c.py at 0.2601 seconds. Notice the |
| 431 |
refraction of the wave front about the boundary between the two layers.} |
| 432 |
\label{fig:ex08cres} |
| 433 |
\end{figure} |
| 434 |
|
| 435 |
|
Parent Directory
| 
Revision Log