SPEED
Functions/Subroutines
MAKE_ABC_FORCE.f90 File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine make_abc_force (nn, xq, wq, dd, rho, lambda, mu, aa11, aa12, aa13, aa21, aa22, aa23, aa31, aa32, aa33, bb11, bb12, bb13, bb21, bb22, bb23, bb31, bb32, bb33, gg1, gg2, gg3, dd1, dd2, dd3, ia, ib, ja, jb, ka, kb, ux, uy, uz, vx, vy, vz, fx, f
 Computes ABC's.
 

Function/Subroutine Documentation

◆ make_abc_force()

subroutine make_abc_force ( integer*4  nn,
real*8, dimension(nn)  xq,
real*8, dimension(nn)  wq,
real*8, dimension(nn,nn)  dd,
real*8, dimension(nn,nn,nn)  rho,
real*8, dimension(nn,nn,nn)  lambda,
real*8, dimension(nn,nn,nn)  mu,
real*8  aa11,
real*8  aa12,
real*8  aa13,
real*8  aa21,
real*8  aa22,
real*8  aa23,
real*8  aa31,
real*8  aa32,
real*8  aa33,
real*8  bb11,
real*8  bb12,
real*8  bb13,
real*8  bb21,
real*8  bb22,
real*8  bb23,
real*8  bb31,
real*8  bb32,
real*8  bb33,
real*8  gg1,
real*8  gg2,
real*8  gg3,
real*8  dd1,
real*8  dd2,
real*8  dd3,
integer*4  ia,
integer*4  ib,
integer*4  ja,
integer*4  jb,
integer*4  ka,
integer*4  kb,
real*8, dimension(nn,nn,nn)  ux,
real*8, dimension(nn,nn,nn)  uy,
real*8, dimension(nn,nn,nn)  uz,
real*8, dimension(nn,nn,nn)  vx,
real*8, dimension(nn,nn,nn)  vy,
real*8, dimension(nn,nn,nn)  vz,
real*8, dimension(nn,nn,nn)  fx,
  f 
)

Computes ABC's.

Author
Ilario Mazzieri
Date
September, 2013
Version
1.0
Parameters
[in]nnnumber of 1D GLL nodes
[in]xqGLL nodes
[in]wqGLL weights
[in]ddspectral derivative matrix
[in]rho
[in]lambda
[in]mumaterial properties given node by node
[in]AA11costant value for the bilinear map
[in]AA12costant value for the bilinear map
[in]AA13costant value for the bilinear map
[in]AA21costant value for the bilinear map
[in]AA22costant value for the bilinear map
[in]AA23costant value for the bilinear map
[in]AA31costant value for the bilinear map
[in]AA32costant value for the bilinear map
[in]AA33costant value for the bilinear map
[in]BB11costant value for the bilinear map
[in]BB12costant value for the bilinear map
[in]BB13costant value for the bilinear map
[in]BB21costant value for the bilinear map
[in]BB22costant value for the bilinear map
[in]BB23costant value for the bilinear map
[in]BB31costant value for the bilinear map
[in]BB32costant value for the bilinear map
[in]BB33costant value for the bilinear map
[in]GG1costant value for the bilinear map
[in]GG2costant value for the bilinear map
[in]GG3costant value for the bilinear map
[in]DD1costant value for the bilinear map
[in]DD2costant value for the bilinear map
[in]DD3costant value for the bilinear map
[in]iaindex for selecting the element absorbing face
[in]ibindex for selecting the element absorbing face
[in]jaindex for selecting the element absorbing face
[in]jbindex for selecting the element absorbing face
[in]kaindex for selecting the element absorbing face
[in]kbindex for selecting the element absorbing face
[in]uxx-displacement
[in]uyy-displacement
[in]uzz-displacement
[in]vxx-velocity
[in]vyy-velocity
[in]vzz-velocity
[out]fxx-absorbing forces
[out]fyy-absorbing forces
[out]fzz-absorbing forces

Definition at line 71 of file MAKE_ABC_FORCE.f90.

77
78!
79!**************************************************************************************************
80
81 implicit none
82
83 integer*4 :: nn
84 integer*4 :: ia,ib,ja,jb,ka,kb
85 integer*4 :: i,j,k,p,q,r
86
87 real*8 :: aa11,aa12,aa13,aa21,aa22,aa23,aa31,aa32,aa33
88 real*8 :: bb11,bb12,bb13,bb21,bb22,bb23,bb31,bb32,bb33
89 real*8 :: gg1,gg2,gg3,dd1,dd2,dd3
90 real*8 :: dxdx,dxdy,dxdz,dydx,dydy,dydz,dzdx,dzdy,dzdz,det_j
91 real*8 :: duxdx,duxdy,duxdz,duydx,duydy,duydz,duzdx,duzdy,duzdz
92 real*8 :: dut1dt1,dut2dt2,dut3dt1,dut3dt2
93 real*8 :: dxdu,dxdv,dydu,dydv,dzdu,dzdv
94 real*8 :: s1,s2,s3,vt1,vt2,vt3,r8t
95 real*8 :: t1x,t1y,t1z,t2x,t2y,t2z,t3x,t3y,t3z
96 real*8 :: t1ux,t2ux,t3ux,t1uy,t2uy,t3uy,t1uz,t2uz,t3uz
97
98 real*8, dimension(nn) :: xq,wq
99
100 real*8, dimension(nn,nn) :: dd
101
102 real*8, dimension(nn,nn,nn) :: rho,lambda,mu
103 real*8, dimension(nn,nn,nn) :: ux,uy,uz
104 real*8, dimension(nn,nn,nn) :: vx,vy,vz
105 real*8, dimension(nn,nn,nn) :: fx,fy,fz
106 real*8, dimension(nn,nn,nn) :: alfa,beta
107
108 do r = 1,nn
109 do q = 1,nn
110 do p = 1,nn
111 alfa(p,q,r) = dsqrt((lambda(p,q,r) + 2.0d0*mu(p,q,r))/rho(p,q,r))
112 beta(p,q,r) = dsqrt(mu(p,q,r)/rho(p,q,r))
113 enddo
114 enddo
115 enddo
116
117 do r = 1,nn
118 do q = 1,nn
119 do p = 1,nn
120 fx(p,q,r) = 0.0d0; fy(p,q,r) = 0.0d0; fz(p,q,r) = 0.0d0
121 enddo
122 enddo
123 enddo
124
125 do r = ka,kb
126 do q = ja,jb
127 do p = ia,ib
128 t1ux = 0.d0; t1uy = 0.d0; t1uz = 0.d0
129 t2ux = 0.d0; t2uy = 0.d0; t2uz = 0.d0
130 t3ux = 0.d0; t3uy = 0.d0; t3uz = 0.d0
131
132 dxdx = aa11 + bb12*xq(r) + bb13*xq(q) + gg1*xq(q)*xq(r)
133 dydx = aa21 + bb22*xq(r) + bb23*xq(q) + gg2*xq(q)*xq(r)
134 dzdx = aa31 + bb32*xq(r) + bb33*xq(q) + gg3*xq(q)*xq(r)
135
136 dxdy = aa12 + bb11*xq(r) + bb13*xq(p) + gg1*xq(r)*xq(p)
137 dydy = aa22 + bb21*xq(r) + bb23*xq(p) + gg2*xq(r)*xq(p)
138 dzdy = aa32 + bb31*xq(r) + bb33*xq(p) + gg3*xq(r)*xq(p)
139
140 dxdz = aa13 + bb11*xq(q) + bb12*xq(p) + gg1*xq(p)*xq(q)
141 dydz = aa23 + bb21*xq(q) + bb22*xq(p) + gg2*xq(p)*xq(q)
142 dzdz = aa33 + bb31*xq(q) + bb32*xq(p) + gg3*xq(p)*xq(q)
143
144 det_j = dxdz * (dydx*dzdy - dzdx*dydy) &
145 - dydz * (dxdx*dzdy - dzdx*dxdy) &
146 + dzdz * (dxdx*dydy - dydx*dxdy)
147
148
149 do i = 1,nn
150 t1ux = t1ux + ux(i,q,r) * dd(p,i)
151 t1uy = t1uy + uy(i,q,r) * dd(p,i)
152 t1uz = t1uz + uz(i,q,r) * dd(p,i)
153 enddo
154
155 do j = 1,nn
156 t2ux = t2ux + ux(p,j,r) * dd(q,j)
157 t2uy = t2uy + uy(p,j,r) * dd(q,j)
158 t2uz = t2uz + uz(p,j,r) * dd(q,j)
159 enddo
160
161 do k = 1,nn
162 t3ux = t3ux + ux(p,q,k) * dd(r,k)
163 t3uy = t3uy + uy(p,q,k) * dd(r,k)
164 t3uz = t3uz + uz(p,q,k) * dd(r,k)
165 enddo
166
167
168 duxdx = 1.0d0 / det_j *(&
169 ((dydy*dzdz - dydz*dzdy) * t1ux) + &
170 ((dydz*dzdx - dydx*dzdz) * t2ux) + &
171 ((dydx*dzdy - dydy*dzdx) * t3ux))
172
173 duydx = 1.0d0 / det_j *(&
174 ((dydy*dzdz - dydz*dzdy) * t1uy) + &
175 ((dydz*dzdx - dydx*dzdz) * t2uy) + &
176 ((dydx*dzdy - dydy*dzdx) * t3uy))
177
178 duzdx = 1.0d0 / det_j *(&
179 ((dydy*dzdz - dydz*dzdy) * t1uz) + &
180 ((dydz*dzdx - dydx*dzdz) * t2uz) + &
181 ((dydx*dzdy - dydy*dzdx) * t3uz))
182
183 duxdy = 1.0d0 / det_j *(&
184 ((dzdy*dxdz - dzdz*dxdy) * t1ux) + &
185 ((dzdz*dxdx - dzdx*dxdz) * t2ux) + &
186 ((dzdx*dxdy - dzdy*dxdx) * t3ux))
187
188 duydy = 1.0d0 / det_j *(&
189 ((dzdy*dxdz - dzdz*dxdy) * t1uy) + &
190 ((dzdz*dxdx - dzdx*dxdz) * t2uy) + &
191 ((dzdx*dxdy - dzdy*dxdx) * t3uy))
192
193 duzdy = 1.0d0 / det_j *(&
194 ((dzdy*dxdz - dzdz*dxdy) * t1uz) + &
195 ((dzdz*dxdx - dzdx*dxdz) * t2uz) + &
196 ((dzdx*dxdy - dzdy*dxdx) * t3uz))
197
198 duxdz = 1.0d0 / det_j *(&
199 ((dxdy*dydz - dxdz*dydy) * t1ux) + &
200 ((dxdz*dydx - dxdx*dydz) * t2ux) + &
201 ((dxdx*dydy - dxdy*dydx) * t3ux))
202
203 duydz = 1.0d0 / det_j *(&
204 ((dxdy*dydz - dxdz*dydy) * t1uy) + &
205 ((dxdz*dydx - dxdx*dydz) * t2uy) + &
206 ((dxdx*dydy - dxdy*dydx) * t3uy))
207
208 duzdz = 1.0d0 / det_j *(&
209 ((dxdy*dydz - dxdz*dydy) * t1uz) + &
210 ((dxdz*dydx - dxdx*dydz) * t2uz) + &
211 ((dxdx*dydy - dxdy*dydx) * t3uz))
212
213 if ((ia.eq.ib).and.(ia.eq.1)) then
214 dxdu = dxdy
215 dydu = dydy
216 dzdu = dzdy
217
218 dxdv = -dxdz
219 dydv = -dydz
220 dzdv = -dzdz
221 endif
222
223 if ((ja.eq.jb).and.(ja.eq.1)) then
224 dxdu = dxdz
225 dydu = dydz
226 dzdu = dzdz
227
228 dxdv = -dxdx
229 dydv = -dydx
230 dzdv = -dzdx
231 endif
232
233 if ((ka.eq.kb).and.(ka.eq.1)) then
234 dxdu = dxdx
235 dydu = dydx
236 dzdu = dzdx
237
238 dxdv = -dxdy
239 dydv = -dydy
240 dzdv = -dzdy
241 endif
242
243 if ((ia.eq.ib).and.(ia.eq.nn)) then
244 dxdu = dxdy
245 dydu = dydy
246 dzdu = dzdy
247
248 dxdv = dxdz
249 dydv = dydz
250 dzdv = dzdz
251 endif
252
253 if ((ja.eq.jb).and.(ja.eq.nn)) then
254 dxdu = dxdz
255 dydu = dydz
256 dzdu = dzdz
257
258 dxdv = dxdx
259 dydv = dydx
260 dzdv = dzdx
261 endif
262
263 if ((ka.eq.kb).and.(ka.eq.nn)) then
264 dxdu = dxdx
265 dydu = dydx
266 dzdu = dzdx
267
268 dxdv = dxdy
269 dydv = dydy
270 dzdv = dzdy
271 endif
272
273 t1x = dxdu
274 t1y = dydu
275 t1z = dzdu
276
277 t2x = dxdv
278 t2y = dydv
279 t2z = dzdv
280
281 t3x = t1y*t2z - t1z*t2y
282 t3y = t1z*t2x - t1x*t2z
283 t3z = t1x*t2y - t1y*t2x
284
285 r8t = dsqrt(t1x*t1x +t1y*t1y +t1z*t1z)
286 t1x = t1x/r8t
287 t1y = t1y/r8t
288 t1z = t1z/r8t
289
290 r8t = dsqrt(t3x*t3x +t3y*t3y +t3z*t3z)
291 t3x = t3x/r8t
292 t3y = t3y/r8t
293 t3z = t3z/r8t
294
295 t2x = t3y*t1z - t3z*t1y
296 t2y = t3z*t1x - t3x*t1z
297 t2z = t3x*t1y - t3y*t1x
298
299
300 dut1dt1 = t1x*t1x*duxdx + t1x*t1y*duydx + t1x*t1z*duzdx &
301 + t1y*t1x*duxdy + t1y*t1y*duydy + t1y*t1z*duzdy &
302 + t1z*t1x*duxdz + t1z*t1y*duydz + t1z*t1z*duzdz
303
304 dut2dt2 = t2x*t2x*duxdx + t2x*t2y*duydx + t2x*t2z*duzdx &
305 + t2y*t2x*duxdy + t2y*t2y*duydy + t2y*t2z*duzdy &
306 + t2z*t2x*duxdz + t2z*t2y*duydz + t2z*t2z*duzdz
307
308 dut3dt1 = t1x*t3x*duxdx + t1x*t3y*duydx + t1x*t3z*duzdx &
309 + t1y*t3x*duxdy + t1y*t3y*duydy + t1y*t3z*duzdy &
310 + t1z*t3x*duxdz + t1z*t3y*duydz + t1z*t3z*duzdz
311
312 dut3dt2 = t2x*t3x*duxdx + t2x*t3y*duydx + t2x*t3z*duzdx &
313 + t2y*t3x*duxdy + t2y*t3y*duydy + t2y*t3z*duzdy &
314 + t2z*t3x*duxdz + t2z*t3y*duydz + t2z*t3z*duzdz
315
316 vt1 = t1x*vx(p,q,r) + t1y*vy(p,q,r) + t1z*vz(p,q,r)
317 vt2 = t2x*vx(p,q,r) + t2y*vy(p,q,r) + t2z*vz(p,q,r)
318 vt3 = t3x*vx(p,q,r) + t3y*vy(p,q,r) + t3z*vz(p,q,r)
319
320 s1 = (mu(p,q,r)*(2.0d0*beta(p,q,r) -alfa(p,q,r))) &
321 /beta(p,q,r) * dut3dt1 - mu(p,q,r)/beta(p,q,r) * vt1
322 s2 = (mu(p,q,r)*(2.0d0*beta(p,q,r) -alfa(p,q,r))) &
323 /beta(p,q,r) * dut3dt2 - mu(p,q,r)/beta(p,q,r) * vt2
324 s3 = (lambda(p,q,r)*beta(p,q,r) +2.0d0*mu(p,q,r)* &
325 (beta(p,q,r) -alfa(p,q,r)))/alfa(p,q,r) &
326 * (dut1dt1 +dut2dt2) - (lambda(p,q,r) +2.0d0*mu(p,q,r))/alfa(p,q,r) * vt3
327
328 det_j = dsqrt(((dxdu*dxdu +dydu*dydu +dzdu*dzdu) &
329 * (dxdv*dxdv +dydv*dydv +dzdv*dzdv)) &
330 -((dxdu*dxdv +dydu*dydv +dzdu*dzdv) &
331 * (dxdu*dxdv +dydu*dydv +dzdu*dzdv)))
332
333 fx(p,q,r) = fx(p,q,r) - wq(p)*wq(q)*wq(r)*det_j/wq(1) &
334 * (s1*t1x + s2*t2x + s3*t3x)
335 fy(p,q,r) = fy(p,q,r) - wq(p)*wq(q)*wq(r)*det_j/wq(1) &
336 * (s1*t1y + s2*t2y + s3*t3y)
337 fz(p,q,r) = fz(p,q,r) - wq(p)*wq(q)*wq(r)*det_j/wq(1) &
338 * (s1*t1z + s2*t2z + s3*t3z)
339
340 enddo
341 enddo
342 enddo
343
344 return
345