SPEED
GET_MAX_VALUES.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
19
45
46 subroutine get_max_values(nmonitpgm, elem_mpgm, local_el_num, ne_loc, &
47 sdeg_mat, nm, cs_loc, cs_nnz_loc, u1, u0, u_1,&
48 omega, nnod_loc, &
49 xr_mpgm,yr_mpgm,zr_mpgm,&
50 max_u, max_v, max_a, max_o)
51
52
53 implicit none
54
55 integer*4 :: nnod_loc, ne_loc, cs_nnz_loc, nm
56 integer*4 :: nn, im, imon, iaz
57 integer*4 :: i, j, k, is, in, id
58 integer*4 :: ie, ielem
59 integer*4, dimension(0:cs_nnz_loc) :: cs_loc
60 integer*4, dimension(nm) :: sdeg_mat
61 integer*4, dimension(ne_loc) :: local_el_num
62
63 real*8 :: uxm, uym, uzm
64 real*8 :: vxm, vym, vzm
65 real*8 :: axm, aym, azm
66
67
68 real*8, dimension(:), allocatable :: ct,ww
69 real*8, dimension(:,:), allocatable :: dd
70 real*8, dimension(:,:,:), allocatable :: ux_el, uy_el, uz_el
71
72 real*8, dimension(3*nnod_loc) :: u_1,u0,u1
73
74 integer*4 :: imonpgm, nmonitpgm,ndt_monitpgm
75 integer*4, dimension(nmonitpgm) :: node_mpgm
76 integer*4, dimension(nmonitpgm) :: elem_mpgm
77 real*8 :: tmp, dt
78 real*8 :: rotang_monitpgm
79 real*8 :: upm,unm,vpm,vnm,apm,anm
80 real*8, dimension(nmonitpgm) :: xr_mpgm,yr_mpgm,zr_mpgm
81 real*8, dimension(nmonitpgm,9) :: max_u,max_v,max_a
82 real*8, dimension(nmonitpgm,3) :: max_o
83
84
85 real*8 :: variable1m,variable2m,variable3m
86 real*8 :: variable4m,variable5m,variable6m
87
88 real*8, dimension(3*nnod_loc) :: omega
89
90 real*8, dimension(:,:,:), allocatable :: variable1_el,variable2_el,variable3_el
91 real*8, dimension(:,:,:), allocatable :: variable4_el,variable5_el,variable6_el
92
93
94
95 upm = 0.0d0; unm = 0.0d0; vpm = 0.0d0; vnm = 0.0d0; apm = 0.0d0; anm = 0.0d0
96
97
98 do imonpgm = 1,nmonitpgm
99
100 ielem = elem_mpgm(imonpgm)
101 call get_indloc_from_indglo(local_el_num, ne_loc, ielem, ie)
102
103 if (ie .ne. 0) then
104
105 im = cs_loc(cs_loc(ie -1) +0)
106 nn = sdeg_mat(im) +1
107 allocate(ct(nn),ww(nn),dd(nn,nn))
108 allocate(ux_el(nn,nn,nn),uy_el(nn,nn,nn),uz_el(nn,nn,nn))
109 allocate(variable1_el(nn,nn,nn),variable2_el(nn,nn,nn),variable3_el(nn,nn,nn))
110 allocate(variable4_el(nn,nn,nn),variable5_el(nn,nn,nn),variable6_el(nn,nn,nn))
111 call make_lgl_nw(nn,ct,ww,dd)
112
113 do k = 1,nn
114 do j = 1,nn
115 do i = 1,nn
116 is = nn*nn*(k -1) +nn*(j -1) +i
117 in = cs_loc(cs_loc(ie -1) + is)
118
119 iaz = 3*(in -1) +1
120 ux_el(i,j,k) = u1(iaz)
121 iaz = 3*(in -1) +2
122 uy_el(i,j,k) = u1(iaz)
123 iaz = 3*(in -1) +3
124 uz_el(i,j,k) = u1(iaz)
125 enddo
126 enddo
127 enddo
128
129 call get_monitor_value(nn,ct,ux_el,&
130 xr_mpgm(imonpgm),yr_mpgm(imonpgm),zr_mpgm(imonpgm),uxm)
131 call get_monitor_value(nn,ct,uy_el,&
132 xr_mpgm(imonpgm),yr_mpgm(imonpgm),zr_mpgm(imonpgm),uym)
133 call get_monitor_value(nn,ct,uz_el,&
134 xr_mpgm(imonpgm),yr_mpgm(imonpgm),zr_mpgm(imonpgm),uzm)
135
136 do k = 1,nn
137 do j = 1,nn
138 do i = 1,nn
139 is = nn*nn*(k -1) +nn*(j -1) +i
140 in = cs_loc(cs_loc(ie -1) + is)
141
142 iaz = 3*(in -1) +1
143 ux_el(i,j,k) = (3.0d0*u1(iaz) - 4.0d0*u0(iaz) + 1.0d0*u_1(iaz)) / (2.0d0*dt)
144 iaz = 3*(in -1) +2
145 uy_el(i,j,k) = (3.0d0*u1(iaz) - 4.0d0*u0(iaz) + 1.0d0*u_1(iaz)) / (2.0d0*dt)
146 iaz = 3*(in -1) +3
147 uz_el(i,j,k) = (3.0d0*u1(iaz) - 4.0d0*u0(iaz) + 1.0d0*u_1(iaz)) / (2.0d0*dt)
148 enddo
149 enddo
150 enddo
151
152 call get_monitor_value(nn,ct,ux_el,&
153 xr_mpgm(imonpgm),yr_mpgm(imonpgm),zr_mpgm(imonpgm),vxm)
154 call get_monitor_value(nn,ct,uy_el,&
155 xr_mpgm(imonpgm),yr_mpgm(imonpgm),zr_mpgm(imonpgm),vym)
156 call get_monitor_value(nn,ct,uz_el,&
157 xr_mpgm(imonpgm),yr_mpgm(imonpgm),zr_mpgm(imonpgm),vzm)
158
159 do k = 1,nn
160 do j = 1,nn
161 do i = 1,nn
162 is = nn*nn*(k -1) +nn*(j -1) +i
163 in = cs_loc(cs_loc(ie -1) + is)
164
165 iaz = 3*(in -1) +1
166 ux_el(i,j,k) = (u1(iaz) -2.0*u0(iaz) +u_1(iaz)) / dt**2
167 iaz = 3*(in -1) +2
168 uy_el(i,j,k) = (u1(iaz) -2.0*u0(iaz) +u_1(iaz)) / dt**2
169 iaz = 3*(in -1) +3
170 uz_el(i,j,k) = (u1(iaz) -2.0*u0(iaz) +u_1(iaz)) / dt**2
171 enddo
172 enddo
173 enddo
174
175 call get_monitor_value(nn,ct,ux_el,&
176 xr_mpgm(imonpgm),yr_mpgm(imonpgm),zr_mpgm(imonpgm),axm)
177 call get_monitor_value(nn,ct,uy_el,&
178 xr_mpgm(imonpgm),yr_mpgm(imonpgm),zr_mpgm(imonpgm),aym)
179 call get_monitor_value(nn,ct,uz_el,&
180 xr_mpgm(imonpgm),yr_mpgm(imonpgm),zr_mpgm(imonpgm),azm)
181
182 if (dabs(uxm).lt.1.0e-30) uxm = 0.0e+00
183 if (dabs(uym).lt.1.0e-30) uym = 0.0e+00
184 if (dabs(uzm).lt.1.0e-30) uzm = 0.0e+00
185 if (dabs(vxm).lt.1.0e-30) vxm = 0.0e+00
186 if (dabs(vym).lt.1.0e-30) vym = 0.0e+00
187 if (dabs(vzm).lt.1.0e-30) vzm = 0.0e+00
188 if (dabs(axm).lt.1.0e-30) axm = 0.0e+00
189 if (dabs(aym).lt.1.0e-30) aym = 0.0e+00
190 if (dabs(azm).lt.1.0e-30) azm = 0.0e+00
191
192 if (rotang_monitpgm.ne.0.0d0) then
193
194 ! upm = Displacement fault Parallel Max
195 ! unm = Displacement fault Normal Max
196
197 upm = uxm * cos(rotang_monitpgm) + uym * sin(rotang_monitpgm)
198 unm = -uxm * sin(rotang_monitpgm) + uym * cos(rotang_monitpgm)
199 vpm = vxm * cos(rotang_monitpgm) + vym * sin(rotang_monitpgm)
200 vnm = -vxm * sin(rotang_monitpgm) + vym * cos(rotang_monitpgm)
201 apm = axm * cos(rotang_monitpgm) + aym * sin(rotang_monitpgm)
202 anm = -axm * sin(rotang_monitpgm) + aym * cos(rotang_monitpgm)
203
204 uxm = upm
205 uym = unm
206 vxm = vpm
207 vym = vnm
208 axm = apm
209 aym = anm
210 endif
211
212 if (dabs(uxm).gt.max_u(imonpgm,1)) max_u(imonpgm,1) = dabs(uxm)
213 if (dabs(uym).gt.max_u(imonpgm,2)) max_u(imonpgm,2) = dabs(uym)
214 if (dabs(uzm).gt.max_u(imonpgm,3)) max_u(imonpgm,3) = dabs(uzm)
215
216 tmp = dsqrt(dabs(uxm*uym))
217 if (tmp.gt.max_u(imonpgm,4)) max_u(imonpgm,4) = tmp
218
219 tmp = ((dabs(uxm)+dabs(uym))/2)
220 if (tmp.gt.max_u(imonpgm,5)) max_u(imonpgm,5) = tmp
221
222 tmp = (dabs(uxm*uym*uzm))**(0.3333333)
223 if (tmp.gt.max_u(imonpgm,6)) max_u(imonpgm,6) = tmp
224
225 tmp = ((dabs(uxm)+dabs(uym)+dabs(uzm))/3)
226 if (tmp.gt.max_u(imonpgm,7)) max_u(imonpgm,7) = tmp
227
228 tmp = dsqrt((uxm)**2 + (uym)**2)
229 if (tmp.gt.max_u(imonpgm,8)) max_u(imonpgm,8) = tmp
230
231 tmp = dsqrt((uxm)**2 + (uym)**2 + (uzm)**2)
232 if (tmp.gt.max_u(imonpgm,9)) max_u(imonpgm,9) = tmp
233
234
235 if (dabs(vxm).gt.max_v(imonpgm,1)) max_v(imonpgm,1) = dabs(vxm)
236 if (dabs(vym).gt.max_v(imonpgm,2)) max_v(imonpgm,2) = dabs(vym)
237 if (dabs(vzm).gt.max_v(imonpgm,3)) max_v(imonpgm,3) = dabs(vzm)
238
239 tmp = dsqrt(dabs(vxm*vym))
240 if (tmp.gt.max_v(imonpgm,4)) max_v(imonpgm,4) = tmp
241
242 tmp = ((dabs(vxm)+dabs(vym))/2)
243 if (tmp.gt.max_v(imonpgm,5)) max_v(imonpgm,5) = tmp
244
245 tmp = (dabs(vxm*vym*vzm))**(0.3333333)
246 if (tmp.gt.max_v(imonpgm,6)) max_v(imonpgm,6) = tmp
247
248 tmp = ((dabs(vxm)+dabs(vym)+dabs(vzm))/3)
249 if (tmp.gt.max_v(imonpgm,7)) max_v(imonpgm,7) = tmp
250
251 tmp = dsqrt((vxm)**2 + (vym)**2)
252 if (tmp.gt.max_v(imonpgm,8)) max_v(imonpgm,8) = tmp
253
254 tmp = dsqrt((vxm)**2 + (vym)**2 + (vzm)**2)
255 if (tmp.gt.max_v(imonpgm,9)) max_v(imonpgm,9) = tmp
256
257 if (dabs(axm).gt.max_a(imonpgm,1)) max_a(imonpgm,1) = dabs(axm)
258 if (dabs(aym).gt.max_a(imonpgm,2)) max_a(imonpgm,2) = dabs(aym)
259 if (dabs(azm).gt.max_a(imonpgm,3)) max_a(imonpgm,3) = dabs(azm)
260
261 tmp = dsqrt(dabs(axm*aym))
262 if (tmp.gt.max_a(imonpgm,4)) max_a(imonpgm,4) = tmp
263
264 tmp = ((dabs(axm)+dabs(aym))/2)
265 if (tmp.gt.max_a(imonpgm,5)) max_a(imonpgm,5) = tmp
266
267 tmp = (dabs(axm*aym*azm))**(0.3333333)
268 if (tmp.gt.max_a(imonpgm,6)) max_a(imonpgm,6) = tmp
269
270 tmp = ((dabs(axm)+dabs(aym)+dabs(azm))/3)
271 if (tmp.gt.max_a(imonpgm,7)) max_a(imonpgm,7) = tmp
272
273 tmp = dsqrt((axm)**2 + (aym)**2)
274 if (tmp.gt.max_a(imonpgm,8)) max_a(imonpgm,8) = tmp
275
276 tmp = dsqrt((axm)**2 + (aym)**2 + (azm)**2)
277 if (tmp.gt.max_a(imonpgm,9)) max_a(imonpgm,9) = tmp
278
279!********************************************************
280! OMEGA
281!********************************************************
282
283 do k = 1,nn
284 do j = 1,nn
285 do i = 1,nn
286 is = nn*nn*(k -1) +nn*(j -1) +i
287 in = cs_loc(cs_loc(ie -1) + is)
288
289 iaz = 3*(in -1) +1
290 variable1_el(i,j,k) = omega(iaz)
291 iaz = 3*(in -1) +2
292 variable2_el(i,j,k) = omega(iaz)
293 iaz = 3*(in -1) +3
294 variable3_el(i,j,k) = omega(iaz)
295 enddo
296 enddo
297 enddo
298
299
300 call get_monitor_value(nn,ct,variable1_el,&
301 xr_mpgm(imonpgm),yr_mpgm(imonpgm),zr_mpgm(imonpgm),variable1m)
302
303 call get_monitor_value(nn,ct,variable2_el,&
304 xr_mpgm(imonpgm),yr_mpgm(imonpgm),zr_mpgm(imonpgm),variable2m)
305
306 call get_monitor_value(nn,ct,variable3_el,&
307 xr_mpgm(imonpgm),yr_mpgm(imonpgm),zr_mpgm(imonpgm),variable3m)
308
309 if (dabs(variable1m).lt.1.0e-30) variable1m = 0.0e+00
310 if (dabs(variable2m).lt.1.0e-30) variable2m = 0.0e+00
311 if (dabs(variable3m).lt.1.0e-30) variable3m = 0.0e+00
312 if (dabs(variable1m).gt.max_o(imonpgm,1)) max_o(imonpgm,1) = dabs(variable1m)
313 if (dabs(variable2m).gt.max_o(imonpgm,2)) max_o(imonpgm,2) = dabs(variable2m)
314 if (dabs(variable3m).gt.max_o(imonpgm,3)) max_o(imonpgm,3) = dabs(variable3m)
315 !-------------------------------------------------------------
316
317 deallocate(ct,ww,dd)
318 deallocate(ux_el,uy_el,uz_el)
319 deallocate(variable1_el,variable2_el,variable3_el)
320 deallocate(variable4_el,variable5_el,variable6_el)
321
322
323 endif
324 enddo
325
326 end subroutine get_max_values
get_indloc_from_indglo
subroutine get_indloc_from_indglo(local_el, nel_loc, ie, ic)
Returns local id from global id.
Definition GET_INDLOC_FROM_INDGLO.f90:30
get_max_values
subroutine get_max_values(nmonitpgm, elem_mpgm, local_el_num, ne_loc, sdeg_mat, nm, cs_loc, cs_nnz_loc, u1, u0, u_1, omega, nnod_loc, xr_mpgm, yr_mpgm, zr_mpgm, max_u, max_v, max_a, max_o)
Computes max values for Peak Ground Map.
Definition GET_MAX_VALUES.f90:51
get_monitor_value
subroutine get_monitor_value(nb_nod, xq, val, xref, yref, zref, re
Computes the mean value.
Definition GET_MONITOR_VALUE.f90:33
make_lgl_nw
subroutine make_lgl_nw(nb_pnt, xq, wq, dd)
Makes Gauss-Legendre-Lobatto nodes, weigths and spectral derivatives.
Definition MAKE_LGL_NW.f90:29