SPEED
Functions/Subroutines
MAKE_ANELASTIC_COEFFICIENTS_NH.f90 File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine make_anelastic_coefficients_nh (nn, n_sls, lambda_el, mu
 Compute anelastic coefficients for Standard Linear Solid model, when the not-honoring strategy is applied USE THE LIBRARY QR_SOLVE!
 

Function/Subroutine Documentation

◆ make_anelastic_coefficients_nh()

subroutine make_anelastic_coefficients_nh ( integer*4  nn,
real*8  n_sls,
  lambda_el,
  mu 
)

Compute anelastic coefficients for Standard Linear Solid model, when the not-honoring strategy is applied 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]lambda_elLame 1st modulus given node by node
[in]mu_elLame 2nd modulus given node by node
[in]QS_nhquality factors for s-waves
[in]QP_nhquality factors for p-waves
[in]f_valvalue for the reference frequency
[in]mpi_idid of the mpi-process
[out]Y_lambda_nhanelastic coefficient (Lame 1st modulus)
[out]Y_mu_nhanelastic 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_NH.f90.

38 rho_el, qs_nh, qp_nh, f_val, &
39 y_lambda_nh, y_mu_nh, frequency_range, mpi_id)
40
41 use speed_exit_codes
42
43 implicit none
44
45 integer*4 :: nn, mpi_id, 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 vp2, vs2, pi
50
51 real*8 :: qs_nh, qp_nh, &
52 y_lambda_nh(1,n_sls),&
53 y_mu_nh(1,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 rho, rho_el(nn,nn,nn), &
59 lambda, lambda_el(nn,nn,nn), &
60 mu, mu_el(nn,nn,nn), &
61 rho_all, mu_all, lambda_all
62
63 rho_all = 0.d0; lambda_all = 0.d0; mu_all = 0.d0
64
65 do i =1, nn
66 do j =1, nn
67 do k =1, nn
68 mu_all = mu_all + mu_el(i,j,k)
69 lambda_all = lambda_all + lambda_el(i,j,k)
70 rho_all = rho_all + rho_el(i,j,k)
71 enddo
72 enddo
73 enddo
74
75 mu = mu_all/(nn**3)
76 lambda = lambda_all/(nn**3)
77 rho = rho_all/(nn**3)
78
79 pi = 4.d0*datan(1.d0);
80
81 if (f_val .le. 10) then
82 fmin = 0.05d0; !fmax = 5.d0;
83 elseif(f_val .le. 100) then
84 fmin = 5.d0; !fmax = 100.d0;
85 elseif(f_val .le. 1000) then
86 fmin = 50.d0; !fmax = 1000.d0;
87 elseif(f_val .le. 10000) then
88 fmin = 500.d0; !fmax = 10000.d0;
89 elseif(f_val .le. 100000) then
90 fmin = 5000.d0; !fmax = 100000.d0;
91 endif
92
93 f_val = 1;
94
95 fmax=100*fmin;
96
97 !if(mpi_id .eq. 0) then
98 ! write(*,*) 'Frequency range where Quality Factor is assumed to be constant is:'
99 ! write(*,*) 'FMIN =' ,fmin, ' FMAX =', fmax;
100 !endif
101
102 esp1 = log10(fmin); esp2 = log10(fmax);
103 deltax = (esp2-esp1)/(n_sls-1);
104
105
106 if (n_sls .eq. 3) then
107
108 frequency_range(1) = f_val*0.1 !2.*pi*fmin;
109 frequency_range(2) = f_val*1 !*fmin;
110 frequency_range(3) = f_val*10 !*fmax;
111 frequency_range_sampling(1) = frequency_range(1);
112 frequency_range_sampling(2) = 0.5*(frequency_range(1)+frequency_range(2))
113 frequency_range_sampling(3) = frequency_range(2);
114 frequency_range_sampling(4) = 0.5*(frequency_range(2)+frequency_range(3))
115 frequency_range_sampling(5) = frequency_range(3);
116
117
118 elseif (n_sls .eq. 4) then
119
120 frequency_range(1) = 2.*pi*fmin;
121 frequency_range(2) = 2.*pi*10**(esp1+deltax)
122 frequency_range(3) = 2.*pi*10**(esp1+2*deltax)
123 frequency_range(4) = 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
132
133 elseif (n_sls .eq. 5) then
134
135 frequency_range(1) = f_val*0.1 ! 2.*pi*fmin;
136 frequency_range(2) = f_val*0.5 !2.*pi*10**(esp1+deltax)
137 frequency_range(3) = f_val*1 !2.*pi*10**(esp1+2*deltax)
138 frequency_range(4) = f_val*5 !2.*pi*10**(esp1+3*deltax)
139 frequency_range(5) = f_val*10 !2.*pi*fmax;
140 frequency_range_sampling(1) = frequency_range(1);
141 frequency_range_sampling(2) = 0.5*(frequency_range(1)+frequency_range(2))
142 frequency_range_sampling(3) = frequency_range(2);
143 frequency_range_sampling(4) = 0.5*(frequency_range(2)+frequency_range(3))
144 frequency_range_sampling(5) = frequency_range(3);
145 frequency_range_sampling(6) = 0.5*(frequency_range(3)+frequency_range(4))
146 frequency_range_sampling(7) = frequency_range(4);
147 frequency_range_sampling(6) = 0.5*(frequency_range(4)+frequency_range(5))
148 frequency_range_sampling(7) = frequency_range(5);
149
150
151 endif
152
153
154
155
156
157 vp2 = (lambda + 2.d0*mu)/rho; vs2 = mu/rho;
158
159 y_mu_nh(1,:) = 0.d0; y_lambda_nh(1,:) = 0.d0;
160 if (qp_nh .ne. 0.d0 .and. qs_nh .ne. 0.d0) then
161
162 do i = 1, 2*n_sls-1
163 do j = 1, n_sls
164
165 a_qp(i,j) = (frequency_range(j)*frequency_range_sampling(i) + frequency_range(j)**2*(1.d0/qp_nh)) &
166 / (frequency_range(j)**2 + frequency_range_sampling(i)**2);
167
168 a_qs(i,j) = (frequency_range(j)*frequency_range_sampling(i) + frequency_range(j)**2*(1.d0/qs_nh)) &
169 / (frequency_range(j)**2 + frequency_range_sampling(i)**2);
170
171 enddo
172
173 rhs_qp(i) = 1.d0/qp_nh;
174 rhs_qs(i) = 1.d0/qs_nh;
175
176 enddo
177
178 !call FACTORIZE_MATRIX(A_QP,N_SLS,PIVOT_QP)
179 !call FACTORIZE_MATRIX(A_QS,N_SLS,PIVOT_QS)
180
181 !call DIRECT_LU_SOLVER(A_QP,RHS_QP,N_SLS,YP,PIVOT_QP)
182 !call DIRECT_LU_SOLVER(A_QS,RHS_QS,N_SLS,YS,PIVOT_QS)
183
184 call qr_solve(2*n_sls-1,n_sls, a_qp, rhs_qp, yp)
185 call qr_solve(2*n_sls-1,n_sls, a_qs, rhs_qs, ys)
186
187
188 y_mu_nh(1,:) = ys;
189 y_lambda_nh(1,:) = (vp2*yp - 2.d0*vs2*ys)/(vp2-2.d0*vs2);
190
191 endif
192
193
194! do i = 1, N_SLS
195! if (Y_lambda_nh(1,i) .lt. 0.d0 ) then
196! write(*,*) 'ERROR! ANELASTIC COEFFICIENTS HAVE TO BE POSITIVE!'
197! write(*,*) 'Y_LAMBDA = ', Y_lambda_nh(1,i)
198! write(*,*) 'SETTING Y_LAMBDA = 0'
199! Y_lambda_nh(1,i) = 0.d0
200! endif
201! if (Y_mu_nh(1,i) .lt. 0.d0) then
202! write(*,*) 'ERROR! ANELASTIC COEFFICIENTS HAVE TO BE POSITIVE!'
203! write(*,*) 'Y_MU = ', Y_mu_nh(1,i)
204! write(*,*) 'SETTING Y_MU = 0'
205! Y_mu_nh(1,i) = 0.d0
206! write(*,*) '----------------------------------------------------'
207! endif
208! enddo
209
210
211
212 return
213
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

References qr_solve().

Here is the call graph for this function: