SPEED
MAKE_NH_Enhanced_ASSIGN_PROP.f90
Go to the documentation of this file.
1
2! Copyright (C) 2012 The SPEED FOUNDATION
3! Author: Ilario Mazzieri
4!
5! This file is part of SPEED.
6!
7! SPEED is free software; you can redistribute it and/or modify it
8! under the terms of the GNU Affero General Public License as
9! published by the Free Software Foundation, either version 3 of the
10! License, or (at your option) any later version.
11!
12! SPEED is distributed in the hope that it will be useful, but
13! WITHOUT ANY WARRANTY; without even the implied warranty of
14! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15! Affero General Public License for more details.
16!
17! You should have received a copy of the GNU Affero General Public License
18! along with SPEED. If not, see <http://www.gnu.org/licenses/>.
19
20
25
36
44
51
52
53 subroutine make_nh_enhanced_assign_prop(nn_loc, nmat, prop_mat, sdeg_mat, &
54 ne_loc, cs_nnz_loc, cs_loc, &
55 node_nhe_flag, count, NN_src_ind_loc, QS, QP, &
56 lambda_nhe, mu_nhe, rho_nhe, &
57 Qs_nhe_el, Qp_nhe_el, mpi_id, mpi_comm) !gamma_dmp_nhe_el
58
59 !use speed_par
60
61 implicit none
62
63 include 'SPEED.MPI'
64
65 integer*4 :: nn_loc, mpi_id, nmat, count, mpi_comm, mpi_ierr
66 integer*4, dimension(nmat) :: sdeg_mat
67 integer*4, dimension(nn_loc) :: node_nhe_flag
68 integer*4, dimension(count) :: NN_src_ind_loc
69
70 real*8, dimension(nmat) :: qs, qp
71 real*8, dimension(nn_loc), intent(inout) :: lambda_nhe, mu_nhe, rho_nhe
72 real*4, dimension(nn_loc) :: qs_nhe, qp_nhe !gamma_dmp_nhe
73
74 integer*4 :: ne_loc, cs_nnz_loc
75 integer*4, dimension(0:cs_nnz_loc) :: cs_loc
76 real*4, dimension(ne_loc) :: qs_nhe_el, qp_nhe_el !gamma_dmp_nhe_el
77
78 character*70 :: file_tomo
79 integer*4 :: npts_tomo, stat, error
80 real*4, dimension(:), allocatable :: tomo_rho, tomo_vs, tomo_vp, tomo_qs, tomo_qp
81 real*8, dimension(nmat,4) :: prop_mat
82
83 real*8 :: t0, t1, time_elapsed, dummy, vs_dum, vp_dum
84 integer*4 :: i, j, ipt, inode, ie
85 integer*4 :: im, istart, iend, nn
86
87 rho_nhe = 0.d0; mu_nhe = 0.d0; lambda_nhe = 0.d0;
88 qs_nhe_el = 0.d0; qp_nhe_el = 0.d0;
89 qs_nhe = 0.d0; qp_nhe = 0.d0;
90
91 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
92 ! Reading Tomography Grid Points data
93 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
94 ! tomo_xyz_mech.in = text file with tomography data
95 ! 1st Line = No. of Tomo Points (npts_tomo)
96 ! 2nd to (npts_tomo+1)th line contains 8 columns:
97 ! x_cord y_cord z_cord Rho Vs Vp Qs Qp (Units must be same as given in *.mate file)
98 file_tomo = 'tomo_xyz_mech.in'
99
100 if(mpi_id .eq.0) then
101 open(124,file=file_tomo)
102 read(124,*)
103 read(124,*) npts_tomo
104 endif
105 call mpi_bcast(npts_tomo,1,speed_integer,0,mpi_comm,mpi_ierr)
106
107 if(mpi_id.eq.0) then
108 allocate(tomo_rho(npts_tomo),tomo_vs(npts_tomo),tomo_vp(npts_tomo),tomo_qs(npts_tomo),tomo_qp(npts_tomo))
109 do ipt = 1, npts_tomo
110 read(124,*)dummy, dummy, dummy, tomo_rho(ipt), tomo_vs(ipt), tomo_vp(ipt), tomo_qs(ipt), tomo_qp(ipt)
111 enddo
112 close(124)
113 endif
114
115
116 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
117 ! Assigning Properties to Each node based on Their Nearest Neighbor
118 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
119
120 !!!!!!!!!!!!!!! Rho !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
121 if (mpi_id.ne.0) then
122 allocate(tomo_rho(npts_tomo))
123 endif
124
125 call mpi_barrier(mpi_comm, mpi_ierr)
126 call mpi_bcast(tomo_rho, npts_tomo, speed_integer, 0, mpi_comm, mpi_ierr)
127
128 i = 0
129 do inode=1,nn_loc
130 if ((node_nhe_flag(inode).eq.999)) then
131 i = i + 1
132 rho_nhe(inode) = tomo_rho(nn_src_ind_loc(i))
133 else
134 rho_nhe(inode) = prop_mat(node_nhe_flag(inode),1)
135 endif
136 enddo
137 deallocate(tomo_rho)
138
139
140 !!!!!!!!!!!!!!!!!!!!!!!!!! Mu !!!!!!!!!!!!!!!!!!!!!!!
141 if (mpi_id.ne.0) then
142 allocate(tomo_vs(npts_tomo))
143 endif
144
145 call mpi_barrier(mpi_comm, mpi_ierr)
146 call mpi_bcast(tomo_vs, npts_tomo, speed_integer, 0, mpi_comm, mpi_ierr)
147
148 i = 0
149 do inode=1,nn_loc
150 if ((node_nhe_flag(inode).eq.999)) then
151 i = i + 1
152 vs_dum = tomo_vs(nn_src_ind_loc(i))
153 mu_nhe(inode) = rho_nhe(inode) * vs_dum**2
154 else
155 mu_nhe(inode) = prop_mat(node_nhe_flag(inode),3)
156 endif
157 enddo
158 deallocate(tomo_vs)
159
160 !!!!!!!!!!!!!!!!!!!!!!!!!! Lambda !!!!!!!!!!!!!!!!!!!!!!!
161 if (mpi_id.ne.0) then
162 allocate(tomo_vp(npts_tomo))
163 endif
164
165 call mpi_barrier(mpi_comm, mpi_ierr)
166 call mpi_bcast(tomo_vp, npts_tomo, speed_integer, 0, mpi_comm, mpi_ierr)
167
168 i = 0
169 do inode=1,nn_loc
170 if ((node_nhe_flag(inode).eq.999)) then
171 i = i + 1
172 vp_dum = tomo_vp(nn_src_ind_loc(i))
173 vs_dum = mu_nhe(inode)/rho_nhe(inode)
174 lambda_nhe(inode) = rho_nhe(inode) * (vp_dum**2 - 2*vs_dum)
175 else
176 lambda_nhe(inode) = prop_mat(node_nhe_flag(inode),2)
177 endif
178 enddo
179 deallocate(tomo_vp)
180
181
182 !!!!!!!!!!!!!!!!!!!!!!!!!! Qs !!!!!!!!!!!!!!!!!!!!!!!
183 if (mpi_id.ne.0) then
184 allocate(tomo_qs(npts_tomo))
185 endif
186
187 call mpi_barrier(mpi_comm, mpi_ierr)
188 call mpi_bcast(tomo_qs, npts_tomo, speed_integer, 0, mpi_comm, mpi_ierr)
189
190 i = 0
191 do inode=1,nn_loc
192 if ((node_nhe_flag(inode).eq.999)) then
193 i = i + 1
194 qs_nhe(inode) = tomo_qs(nn_src_ind_loc(i))
195 else
196 qs_nhe(inode) = qs(node_nhe_flag(inode))
197 endif
198 enddo
199 deallocate(tomo_qs)
200
201 !!!!!!!!!!!!!!!!!!!!!!!!!! Qp !!!!!!!!!!!!!!!!!!!!!!!
202 if (mpi_id.ne.0) then
203 allocate(tomo_qp(npts_tomo))
204 endif
205
206 call mpi_barrier(mpi_comm, mpi_ierr)
207 call mpi_bcast(tomo_qp, npts_tomo, speed_integer, 0, mpi_comm, mpi_ierr)
208
209 i = 0
210 do inode=1,nn_loc
211 if ((node_nhe_flag(inode).eq.999)) then
212 i = i + 1
213 qp_nhe(inode) = tomo_qp(nn_src_ind_loc(i))
214 else
215 qp_nhe(inode) = qp(node_nhe_flag(inode))
216 endif
217 enddo
218 deallocate(tomo_qp)
219
220 call mpi_barrier(mpi_comm, mpi_ierr)
221
222 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
223 !!!!! Making Qs, Qp for Element instead of Each node (To save Memory)
224
225 qs_nhe_el = 0;
226 qp_nhe_el = 0;
227
228 do ie = 1,ne_loc
229 im = cs_loc(cs_loc(ie -1))
230 nn = sdeg_mat(im) + 1
231 istart = cs_loc(ie-1) + 1
232 iend = cs_loc(ie) - 1
233
234 do j = istart,iend
235 qs_nhe_el(ie) = qs_nhe_el(ie) + qs_nhe(cs_loc(j))
236 qp_nhe_el(ie) = qp_nhe_el(ie) + qp_nhe(cs_loc(j))
237 enddo
238 qs_nhe_el(ie) = qs_nhe_el(ie)/(nn**3)
239 qp_nhe_el(ie) = qp_nhe_el(ie)/(nn**3)
240 enddo !ie=1,ne_loc
241
242 end subroutine MAKE_NH_Enhanced_ASSIGN_PROP
make_nh_enhanced_assign_prop
subroutine make_nh_enhanced_assign_prop(nn_loc, nmat, prop_mat, sdeg_mat, ne_loc, cs_nnz_loc, cs_loc, node_nhe_flag, count, nn_src_ind_loc, qs, qp, lambda_nhe, mu_nhe, rho_nhe, qs_nhe_el, qp_nhe_el, mpi_id, mpi_comm)
...Not-Honoring Enhanced (NHE) Implementation
Definition MAKE_NH_Enhanced_ASSIGN_PROP.f90:58