SPEED
MAKE_ANELASTIC_COEFFICIENTS.f90
Go to the documentation of this file.
1! Copyright (C) 2012 The SPEED FOUNDATION
2! Author: Ilario Mazzieri
3!
4! This file is part of SPEED.
5!
6! SPEED is free software; you can redistribute it and/or modify it
7! under the terms of the GNU Affero General Public License as
8! published by the Free Software Foundation, either version 3 of the
9! License, or (at your option) any later version.
10!
11! SPEED is distributed in the hope that it will be useful, but
12! WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14! Affero General Public License for more details.
15!
16! You should have received a copy of the GNU Affero General Public License
17! along with SPEED. If not, see <http://www.gnu.org/licenses/>.
18
35
36
37 subroutine make_anelastic_coefficients(nmat, N_SLS, prop_mat, QS, QP, f_val, &
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
209 end subroutine make_anelastic_coefficients
210
211
212! Copyright (C) 2012 The SPEED FOUNDATION
213! Author: Ilario Mazzieri
214!
215! This file is part of SPEED.
216!
217! SPEED is free software; you can redistribute it and/or modify it
218! under the terms of the GNU Affero General Public License as
219! published by the Free Software Foundation, either version 3 of the
220! License, or (at your option) any later version.
221!
222! SPEED is distributed in the hope that it will be useful, but
223! WITHOUT ANY WARRANTY; without even the implied warranty of
224! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
225! Affero General Public License for more details.
226!
227! You should have received a copy of the GNU Affero General Public License
228! along with SPEED. If not, see <http://www.gnu.org/licenses/>.
229
230 subroutine factorize_matrix(A,ndim,pivot)
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
313 end subroutine factorize_matrix
314
315! Copyright (C) 2012 The SPEED FOUNDATION
316! Author: Ilario Mazzieri
317!
318! This file is part of SPEED.
319!
320! SPEED is free software; you can redistribute it and/or modify it
321! under the terms of the GNU Affero General Public License as
322! published by the Free Software Foundation, either version 3 of the
323! License, or (at your option) any later version.
324!
325! SPEED is distributed in the hope that it will be useful, but
326! WITHOUT ANY WARRANTY; without even the implied warranty of
327! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
328! Affero General Public License for more details.
329!
330! You should have received a copy of the GNU Affero General Public License
331! along with SPEED. If not, see <http://www.gnu.org/licenses/>.
332
333
334 subroutine direct_lu_solver(A,tn,ndim,x_sol,pivot)
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
375 end subroutine direct_lu_solver
376
377
378
379
380
381
382
383
make_anelastic_coefficients
subroutine make_anelastic_coefficients(nmat, n_sls, prop_mat, qs,
Compute anelastic coefficients for Standard Linear Solid model USE THE LIBRARY QR_SOLVE!
Definition MAKE_ANELASTIC_COEFFICIENTS.f90:38
factorize_matrix
subroutine factorize_matrix(a, ndim, pivot)
Definition MAKE_ANELASTIC_COEFFICIENTS.f90:231
direct_lu_solver
subroutine direct_lu_solver(a, tn, ndim, x_sol, pivot)
Definition MAKE_ANELASTIC_COEFFICIENTS.f90:335
qr_solve
subroutine qr_solve(m, n, a, b, x)
Definition QR_SOLVE.f90:2384
speed_exit_codes
SPEED exit codes.
Definition MODULES.f90:25
speed_exit_codes::exit_singularmtx
integer, parameter exit_singularmtx
Definition MODULES.f90:34