SPEED
Functions/Subroutines
MAKE_ANELASTIC_COEFFICIENTS.f90 File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine make_anelastic_coefficients (nmat, n_sls, prop_mat, qs,
 Compute anelastic coefficients for Standard Linear Solid model USE THE LIBRARY QR_SOLVE!
 
subroutine factorize_matrix (a, ndim, pivot)
 
subroutine direct_lu_solver (a, tn, ndim, x_sol, pivot)
 

Function/Subroutine Documentation

◆ direct_lu_solver()

subroutine direct_lu_solver ( real*8, dimension(ndim,ndim), intent(inout)  a,
real*8, dimension(ndim), intent(inout)  tn,
integer*4, intent(in)  ndim,
real*8, dimension(ndim), intent(inout)  x_sol,
integer*4, dimension(ndim), intent(in)  pivot 
)

Definition at line 334 of file MAKE_ANELASTIC_COEFFICIENTS.f90.

335
336
337 integer*4, intent(in) :: ndim
338 real*8, dimension(ndim), intent(inout) :: x_sol
339 real*8, dimension(ndim,ndim), intent(inout) :: a
340 real*8, dimension(ndim), intent(inout) :: tn
341 integer*4, dimension(ndim), intent(in) :: pivot
342
343 integer*4 :: ii,ij,ik,iflag,ip
344 real*8, dimension(ndim) ::y_sol
345 real*8, dimension(ndim) :: dd
346 real*8 :: aw,cmax,rmax,temp,sun
347
348 if (ndim .le. 1) then
349 x_sol(1) = tn(1) / a(1, 1)
350 return
351 endif
352
353 ip = pivot(1)
354 y_sol(1) = tn(ip)
355 do ii=2,ndim
356 sun = 0.d0
357 do ij=1,ii-1
358 sun = a(ii, ij) * y_sol(ij) + sun
359 enddo
360 ip = pivot(ii)
361 y_sol(ii) = tn(ip) - sun
362 enddo
363
364 x_sol(ndim) = y_sol(ndim) / a(ndim, ndim)
365
366 do ii=ndim-1,1,-1
367 sun = 0.d0
368 do ij=ii+1,ndim
369 sun = a(ii, ij) * x_sol(ij) + sun
370 enddo
371 x_sol(ii) = (y_sol(ii) - sun) / a(ii, ii)
372 enddo
373
374 return

◆ factorize_matrix()

subroutine factorize_matrix ( real*8, dimension(ndim,ndim), intent(inout)  a,
integer*4, intent(in)  ndim,
integer*4, dimension(ndim), intent(out)  pivot 
)

Definition at line 230 of file MAKE_ANELASTIC_COEFFICIENTS.f90.

231
232 use speed_exit_codes
233
234 implicit none
235
236 real*8, dimension(ndim,ndim), intent(inout) :: a
237 integer*4, intent(in) :: ndim
238 integer*4, dimension(ndim), intent(out) :: pivot
239
240
241 real*8, dimension(ndim) :: dd
242 integer*4 :: ii,ij,ik,iflag,ip,is
243 real*8 :: aw,cmax,rmax,temp,sun
244
245 pivot(:) = 0
246 iflag = 1
247 do ii=1,ndim
248 pivot(ii) = ii
249 rmax = 0.d0
250 do ij=1,ndim
251 rmax = max(rmax, abs(a(ii, ij)))
252 enddo
253
254 if (rmax.eq.0) then
255 iflag = 0
256 rmax = 1.d0
257 endif
258 dd(ii) = rmax
259 enddo
260
261 if (ndim.le.1) return
262
263 do ik=1,ndim-1
264
265 cmax = abs(a(ik, ik)) / dd(ik)
266 is=ik
267 do ii=ik+1,ndim
268 aw = abs(a(ii,ik)) / dd(ii)
269 if (aw.gt.cmax) then
270 cmax = aw
271 is=ii
272 endif
273 enddo
274
275 if (cmax.eq.0) then
276 iflag = 0
277 else
278 if (is.gt.ik) then
279 iflag = -iflag
280 ii=pivot(is)
281 pivot(is)=pivot(ik)
282 pivot(ik) = ii
283 temp=dd(is)
284 dd(is)=dd(ik)
285 dd(ik) = temp
286 do ij=1,ndim
287 temp =a(is,ij)
288 a(is,ij)=a(ik,ij)
289 a(ik, ij) = temp
290 enddo
291
292 endif
293
294 do ii=ik+1,ndim
295 a(ii, ik) = a(ii, ik) / a(ik, ik)
296 do ij=ik+1,ndim
297 a(ii, ij) = a(ii, ij) - a(ii, ik) * a(ik, ij)
298 enddo
299 enddo
300
301 endif
302 enddo
303 if (a(ndim,ndim) .eq. 0.d0) then
304 iflag = 0
305 endif
306
307 if (iflag .eq. 0 ) then
308 write(*,*) 'SINGULAR MATRIX'
309 call exit(exit_singularmtx)
310 endif
311
312
speed_exit_codes
SPEED exit codes.
Definition MODULES.f90:25
speed_exit_codes::exit_singularmtx
integer, parameter exit_singularmtx
Definition MODULES.f90:34

References speed_exit_codes::exit_singularmtx.

◆ make_anelastic_coefficients()

subroutine make_anelastic_coefficients ( integer*4  nmat,
real*8  n_sls,
real*8, dimension(nmat,4)  prop_mat,
real*8, dimension(nmat)  qs 
)

Compute anelastic coefficients for Standard Linear Solid model USE THE LIBRARY QR_SOLVE!

Author
Ilario Mazzieri
Date
November, 2014
Version
1.0
Parameters
[in]nmatnumber of materials
[in]N_SLSnumber of subdivision in the standard linear solid (dafult N_SLS=3) possible values are N_SLS=3,4,5 (see MODULES.f90)
[in]prop_matmaterial properties
[in]QSquality factors for s-waves
[in]QPquality factors for p-waves
[in]f_valvalue for the reference frequency
[in]mpi_idid of the mpi-process
[out]Y_lambdaanelastic coefficient (Lame 1st modulus)
[out]Y_muanelastic coefficient (Lame 2nd modulus)
[out]frequency_rangefrequency range where quality factors are assumed to be constant

Definition at line 37 of file MAKE_ANELASTIC_COEFFICIENTS.f90.

38 y_lambda, y_mu, frequency_range, mpi_id)
39
40
41 use speed_exit_codes
42
43 implicit none
44
45 integer*4 :: nmat,mpi_id, im, N_SLS, i, j, k, &
46 pivot_qp(n_sls),pivot_qs(n_sls)
47
48 real*8 :: f_val, fmin, fmax, esp1, esp2, deltax, &
49 rho, lambda, mu, vp2, vs2, pi
50
51 real*8 :: prop_mat(nmat,4), qs(nmat), qp(nmat), &
52 y_lambda(nmat,n_sls),&
53 y_mu(nmat,n_sls), &
54 frequency_range(n_sls),&
55 frequency_range_sampling(2*n_sls-1),&
56 ys(n_sls),yp(n_sls), &
57 a_qp(2*n_sls-1,n_sls),rhs_qp(2*n_sls-1),a_qs(2*n_sls-1,n_sls),rhs_qs(2*n_sls-1)
58
59 pi = 4.d0*datan(1.d0);
60
61 if (f_val .le. 10) then
62 fmin = 0.05d0; !fmax = 5.d0;
63 elseif(f_val .le. 100) then
64 fmin = 5.d0; !fmax = 100.d0;
65 elseif(f_val .le. 1000) then
66 fmin = 50.d0; !fmax = 1000.d0;
67 elseif(f_val .le. 10000) then
68 fmin = 500.d0; !fmax = 10000.d0;
69 elseif(f_val .le. 100000) then
70 fmin = 5000.d0; !fmax = 100000.d0;
71 endif
72
73 fmax=100*fmin;
74
75 !fmin and fmax not used - fixed value for fmin = f_val*0.1 and fmax = f_val*10
76
77 f_val = 1;
78
79 if(mpi_id .eq. 0) then
80 write(*,*) 'Frequency range where Quality Factor is assumed to be constant is:'
81 write(*,*) 'FMIN =' ,f_val*0.1, ' FMAX =', f_val*10;
82 endif
83
84
85 esp1 = log10(fmin); esp2 = log10(fmax);
86 deltax = (esp2-esp1)/(n_sls-1);
87
88
89
90 if (n_sls .eq. 3) then
91
92 frequency_range(1) = f_val*0.1 !2.*pi*fmin;
93 frequency_range(2) = f_val*1 !*fmin;
94 frequency_range(3) = f_val*10 !*fmax;
95 frequency_range_sampling(1) = frequency_range(1);
96 frequency_range_sampling(2) = 0.5*(frequency_range(1)+frequency_range(2))
97 frequency_range_sampling(3) = frequency_range(2);
98 frequency_range_sampling(4) = 0.5*(frequency_range(2)+frequency_range(3))
99 frequency_range_sampling(5) = frequency_range(3);
100
101
102 elseif (n_sls .eq. 4) then
103
104 frequency_range(1) = 2.*pi*fmin;
105 frequency_range(2) = 2.*pi*10**(esp1+deltax)
106 frequency_range(3) = 2.*pi*10**(esp1+2*deltax)
107 frequency_range(4) = 2.*pi*fmax;
108 frequency_range_sampling(1) = frequency_range(1);
109 frequency_range_sampling(2) = 0.5*(frequency_range(1)+frequency_range(2))
110 frequency_range_sampling(3) = frequency_range(2);
111 frequency_range_sampling(4) = 0.5*(frequency_range(2)+frequency_range(3))
112 frequency_range_sampling(5) = frequency_range(3);
113 frequency_range_sampling(6) = 0.5*(frequency_range(3)+frequency_range(4))
114 frequency_range_sampling(7) = frequency_range(4);
115
116
117 elseif (n_sls .eq. 5) then
118
119 frequency_range(1) = f_val*0.1 ! 2.*pi*fmin;
120 frequency_range(2) = f_val*0.5 !2.*pi*10**(esp1+deltax)
121 frequency_range(3) = f_val*1 !2.*pi*10**(esp1+2*deltax)
122 frequency_range(4) = f_val*5 !2.*pi*10**(esp1+3*deltax)
123 frequency_range(5) = f_val*10 !2.*pi*fmax;
124 frequency_range_sampling(1) = frequency_range(1);
125 frequency_range_sampling(2) = 0.5*(frequency_range(1)+frequency_range(2))
126 frequency_range_sampling(3) = frequency_range(2);
127 frequency_range_sampling(4) = 0.5*(frequency_range(2)+frequency_range(3))
128 frequency_range_sampling(5) = frequency_range(3);
129 frequency_range_sampling(6) = 0.5*(frequency_range(3)+frequency_range(4))
130 frequency_range_sampling(7) = frequency_range(4);
131 frequency_range_sampling(8) = 0.5*(frequency_range(4)+frequency_range(5))
132 frequency_range_sampling(9) = frequency_range(5);
133
134
135 endif
136
137 if(mpi_id .eq. 0) then
138 write(*,*) 'Sampled frequencies:'
139 write(*,*) 'FREQ =' ,frequency_range_sampling
140 endif
141
142
143
144 do im = 1, nmat
145
146 rho = prop_mat(im,1); lambda=prop_mat(im,2); mu=prop_mat(im,3);
147 vp2 = (lambda + 2.d0*mu)/rho; vs2 = mu/rho;
148
149 y_mu(im,:) = 0.d0; y_lambda(im,:) = 0.d0;
150 if (qp(im) .ne. 0.d0 .and. qs(im) .ne. 0.d0) then
151
152
153 do i = 1, 2*n_sls-1
154 do j = 1, n_sls
155
156 a_qp(i,j) = (frequency_range(j)*frequency_range_sampling(i) + frequency_range(j)**2*(1.d0/qp(im))) &
157 / (frequency_range(j)**2 + frequency_range_sampling(i)**2);
158
159 a_qs(i,j) = (frequency_range(j)*frequency_range_sampling(i) + frequency_range(j)**2*(1.d0/qs(im))) &
160 / (frequency_range(j)**2 + frequency_range_sampling(i)**2);
161
162 enddo
163
164 rhs_qp(i) = 1.d0/qp(im);
165 rhs_qs(i) = 1.d0/qs(im);
166
167 enddo
168
169
170 !call FACTORIZE_MATRIX(A_QP,N_SLS,PIVOT_QP)
171 !call FACTORIZE_MATRIX(A_QS,N_SLS,PIVOT_QS)
172
173 !call DIRECT_LU_SOLVER(A_QP,RHS_QP,N_SLS,YP,PIVOT_QP)
174 !call DIRECT_LU_SOLVER(A_QS,RHS_QS,N_SLS,YS,PIVOT_QS)
175
176
177
178 call qr_solve(2*n_sls-1,n_sls, a_qp, rhs_qp, yp)
179 call qr_solve(2*n_sls-1,n_sls, a_qs, rhs_qs, ys)
180
181 y_mu(im,:) = ys;
182 y_lambda(im,:) = (vp2*yp - 2.d0*vs2*ys)/(vp2-2.d0*vs2);
183
184 endif
185
186! do i = 1, N_SLS
187! if (mpi_id .eq. 0 .and.Y_lambda(im,i) .lt. 0.d0 .or. Y_mu(im,i) .lt. 0) then
188! write(*,*) 'ERROR! ANELASTIC COEFFICIENTS HAVE TO BE POSITIVE!'
189! write(*,*) 'ANELASTIC COEFFICIENTS FOR MATERIAL ', im, ' ARE '
190! write(*,*) 'Y_MU = ', Y_mu(im,:)
191! write(*,*) 'Y_LAMBDA = ', Y_lambda(im,:)
192! write(*,*) 'Y_P = ', YP
193! write(*,*) '----------------------------------------------------'
194! endif
195! enddo
196
197 if(mpi_id .eq. 0) then
198 write(*,*) 'ANELASTIC COEFFICIENTS FOR MATERIAL ', im, ' ARE '
199 write(*,*) 'Y_MU = ', y_mu(im,:)
200 write(*,*) 'Y_LAMBDA = ', y_lambda(im,:)
201 write(*,*) 'Y_P = ', yp
202 write(*,*) '----------------------------------------------------'
203 endif
204
205 enddo
206
207 return
208
qr_solve
subroutine qr_solve(m, n, a, b, x)
Definition QR_SOLVE.f90:2384

References qr_solve().

Here is the call graph for this function: