SPEED
TIME_LOOP.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
26
27 subroutine time_loop(el_new)
28
29 use max_var
30 use dgjump
31 use speed_timeloop
32 use speed_sci
33
34 implicit none
35
36 include 'SPEED.MPI'
37
38 type(el4loop), dimension(nelem_dg), intent(in) :: el_new
39
40! GENERAL
41 character*70 :: file_monitor, file_outU1, file_outU0, file_outV1, &
42 file_outsv, filepg, file_trav_load, file_outxyz
43 character*70 :: file_fe0, file_fe1
44 integer*4 :: unit_fe0, unit_fe1
45
46 integer*4 :: mm, iface, ifacene, iene, imne, ie_curr, iene_curr, trofa, posit, jshift, &
47 ic, ih, ik, il, it, ipos, ineg, &
48 isnap, its, fn, nn, ip, im, imon, iaz, &
49 i, j, k, is, in, id, istage, itime, &
50 ie, ielem, count_monitor, find_tag, icase, irand, isdof
51
52! Debug SSI
53 character*70 :: file_getval_name
54 integer*4 :: file_getval_id
55 ! real*8 ::dum_ind
56
57! Elapsed time print-out iteration divisor
58 integer*4 :: it_divisor = 1000
59
60 integer*4, dimension(:), allocatable :: node_counter, vtk_numbering
61
62
63 real*8 :: start_time, final_time, time_in_seconds, &
64 start, finish, &
65 tt0, tt1, tt2, deltat2, tt_int, &
66 get_func_value
67
68 real*8 :: tol = 1.d-1
69
70
71 real*8, dimension(:), allocatable :: jumpx, jumpy, jumpz
72 real*8, dimension(:), allocatable :: func_value
73
74! ELEMENT BY ELEMENT
75 real*8 :: rho, lambda, mu, tref
76
77 real*8, dimension(:), allocatable :: ct,ww
78 real*8, dimension(:,:), allocatable :: dd
79 real*8, dimension(:,:,:), allocatable :: ux_el, uy_el, uz_el, &
80 duxdx_el, duydx_el, duzdx_el, &
81 duxdy_el, duydy_el, duzdy_el, &
82 duxdz_el, duydz_el, duzdz_el, &
83 sxx_el, syy_el, szz_el, &
84 syz_el, szx_el, sxy_el, &
85 fx_el, fy_el, fz_el, &
86 mv_el, &
87 vx_el, vy_el, vz_el, &
88 div_el, rotx_el, roty_el, rotz_el, &
89 strainxx_el, strainyy_el, strainzz_el, &
90 strainxy_el, strainyz_el, strainzx_el, &
91 stressxx_el, stressyy_el, stresszz_el, &
92 stressxy_el, stressyz_el, stresszx_el, &
93 rho_el, lambda_el, mu_el, gamma_el
94
95! MPI
96 integer*4, dimension(mpi_np) :: send_length, recv_length, small_send_length, small_recv_length
97 integer*4, dimension(mpi_np) :: send_length_jump, recv_length_jump, double_send_length, double_recv_length
98
99 integer*4, dimension(:), allocatable :: small_send_buffer,small_recv_buffer
100
101 real*8, dimension(:), allocatable :: u_1, u0, u1, u2, v1, v2, jump, jump_minus, &
102 u_p, u_m, jump_all, &
103 fk, fe, mv, &
104 send_buffer, recv_buffer, &
105 double_send_buffer, double_recv_buffer, &
106 send_buffer_jump, recv_buffer_jump, sdofrecv_temp
107
108! DAMPING
109 real*8 :: gamma, qs_nh, qp_nh
110 real*8, dimension(1,N_SLS) :: y_lambda_nh, y_mu_nh
111 real*8, dimension(:), allocatable :: mc, mck, damp_term
112 real*8, dimension(:,:,:), allocatable :: mc_el, mck_el
113 real*8, dimension(:,:), allocatable :: strain_visc
114 real*8, dimension(:,:,:,:), allocatable :: strain_visc_xx,strain_visc_xy,&
115 strain_visc_xz,strain_visc_yy,&
116 strain_visc_yz,strain_visc_zz
117
118
119! PEAK GROUND MAP
120 real*8, dimension (:,:), allocatable :: max_u,max_v,max_a, max_o
121
122! NON LINEAR ELASTIC MATERIAL
123 integer*4 :: im_nle, which_mat_nle, y_or_n_node
124 integer*4, dimension(:), allocatable :: con_spx_loc_nle
125 integer*4, dimension(:), allocatable :: node_nle_4_mpi
126 integer*4, dimension(:,:,:), allocatable :: yes_or_not_nle_el
127 real*8 :: depth_nle
128 real*8, dimension(:,:,:), allocatable :: r_el
129 real*8, dimension(:), allocatable :: depth_nle_el
130
131! OUTPUT
132 real*8, dimension(:), allocatable :: strain, omega, stress, strainold !Added strainold - SSI AH
133
134! RUNGE-KUTTA
135 real*8, dimension(:), allocatable :: a, b,c, u_int, v_int
136 real*8, dimension(:,:), allocatable :: kappau, kappav
137
138! DEBUG
139 real*8, dimension(:,:,:), allocatable :: gamma_new
140 real*8, dimension(:), allocatable :: mu_nle, gamma_nle
141 real*8 :: val_debug
142
143! TRAVELLING LOAD
144 integer*4 :: nb_tra_load
145 integer*4, dimension(:), allocatable :: node_tra
146 real*8, dimension(:), allocatable :: dist_tra
147
148! TEST NOT-HONORING
149 integer*4 :: check_not_hon
150
151! INSTABILITY CONTROL
152 ! Set to TRUE when instability control is enabled and instability is detected.
153 logical :: b_instability_abort
154 logical, dimension(mpi_np) :: b_instability_abort_all
155
156 b_instability_abort = .false.
157
158!---------------------------------------------------------------------------
159
160 if(mpi_id .eq. 0) start = mpi_wtime()
161
162
163 allocate(func_value(nfunc))
164
165!---------------------------------------------------------------------------
166! Monitor LiST
167!---------------------------------------------------------------------------
168
169 if (nmonitors_lst.ge.1) then
170
171 call make_monitor_files(nmonitors_lst,el_monitor_lst,local_el_num, nelem_loc,&
172 count_monitor, file_monitor,monitor_file, mpi_id, &
173 opt_out_var)
174
175 endif
176
177 if(nmonitors_pgm .ge. 1) then
178 allocate(max_u(nmonitors_pgm,9), max_v(nmonitors_pgm,9), max_a(nmonitors_pgm,9), max_o(nmonitors_pgm,3))
179 max_u = 0.d0; max_v = 0.d0; max_a = 0.d0; max_o = 0.d0
180 endif
181
182!---------------------------------------------------------------------------
183 ! INITIALIZE OSCILLATOR AH
184 !--------------------------------------------------------------------------
185 if (sdofnum.gt.0) then
186 !n_bld is greater than 1 only for mpi_id =0
187 ! The force-calculations od SDOF are done only using mpi_id = 0 process
188 call read_sdof_input_files
189
190 allocate(sdofinputd(3*sdofnum)) ! base displacement
191 allocate(sdofinput(3*sdofnum)) !base acceleration (x,y,z)
192 allocate(sdofinputbuffer(3*mpi_np*sdofnum))
193 allocate(sdofrecv_temp(3*sdofnum))
194 allocate(sdofforceinput(3*sdofnum)) !(Fx,Fy,Fz)
195 sdofforceinput=0
196 allocate(sdofforceinputbuffer(3*sdofnum*mpi_np))
197
198 call make_sdof_output_files
199
200 !Debug - Writing files with values of Interaction forces at each time step - SS
201 ! if (mpi_id.eq.0) then
202 ! do isdof=1,SDOFnum
203 ! file_getval_id = 2433 + isdof
204 ! write(file_getval_name,'(A12,I4.4,A4)') 'GETVAL_SDOF_', isdof, '.txt'
205 ! open(file_getval_id, file=file_getval_name, status='replace')
206 ! close(file_getval_id)
207 ! enddo
208 ! endif
209
210 endif
211!---------------------------------------------------------------------------
212! Travelling Load
213!---------------------------------------------------------------------------
214 if(nload_traz_el .gt. 0) then
215 file_trav_load = 'TRAVPOINTS.input'
216 call read_dime_filepg(file_trav_load,nb_tra_load)
217 !write(*,*) nb_tra_load
218 !read(*,*)
219 allocate(node_tra(nb_tra_load), dist_tra(nb_tra_load))
220
221 call read_file_trav_point(file_trav_load, nb_tra_load, node_tra, dist_tra)
222
223
224 ! if(mpi_id .eq. 0) then
225 ! write(*,*) nb_tra_load
226 ! write(*,*) node_tra
227 ! write(*,*) dist_tra
228 ! endif
229 ! call MPI_BARRIER(mpi_comm,mpi_ierr)
230 endif
231
232
233!---------------------------------------------------------------------------
234! MPI INITIALIZATION --- NEW NUMERATION recv & send
235!---------------------------------------------------------------------------
236
237 allocate(send_buffer(3*nsend)); send_buffer= 0.0d0
238 allocate(recv_buffer(3*nrecv)); recv_buffer = 0.0d0
239
240 do i = 1,mpi_np
241 send_length(i) = 3*proc_send(i); recv_length(i) = 3*proc_recv(i)
242 if(opt_out_var(4) .eq. 1 .or. opt_out_var(5) .eq. 1 .or. damping_type .eq. 2) then
243 double_send_length(i) = 6*proc_send(i); double_recv_length(i) = 6*proc_recv(i)
244 endif
245 enddo
246
247 do i = 1,nrecv
248 in = node_recv(i)
249 call get_indloc_from_indglo(local_node_num, nnod_loc, in, id)
250 node_recv(i) = id
251 enddo
252
253 do i = 1,nsend
254 in = node_send(i)
255 call get_indloc_from_indglo(local_node_num, nnod_loc, in, id)
256 node_send(i) = id
257 enddo
258
259!---------------------------------------------------------------------------
260! MPI INITIALIZATION --- NON CONFORMING DG
261!---------------------------------------------------------------------------
262
263 if (nelem_dg_glo .gt. 0) then
264
265 allocate(send_buffer_jump(3*nsend_jump), recv_buffer_jump(3*nrecv_jump))
266 send_buffer_jump = 0.0d0; recv_buffer_jump = 0.0d0
267
268 do i = 1,mpi_np
269 send_length_jump(i) = 3*proc_send_jump(i)
270 recv_length_jump(i) = 3*proc_recv_jump(i)
271 enddo
272
273 allocate(jump_minus(3*nrecv_jump))
274
275 do i = 1,nsend_jump
276 in = node_send_jump(i)
277 call get_indloc_from_indglo(local_node_num, nnod_loc, in, id)
278 node_send_jump(i) = id
279 enddo
280
281 endif
282
283!---------------------------------------------------------------------------
284! UNKNOWNS INITIALIZATION
285!---------------------------------------------------------------------------
286
287 allocate(fe(3*nnod_loc),fk(3*nnod_loc),mv(3*nnod_loc)); fe = 0.0d0; fk = 0.0d0; mv = 0.0d0
288
289 if(rk_scheme .eq. 'RUNGEKUTTA') then
290 allocate(u1(3*nnod_loc),u2(3*nnod_loc),v1(3*nnod_loc),v2(3*nnod_loc),jump(3*nnod_loc))
291 u1 = 0.0d0; u2 = 0.0d0; v1 = 0.0d0; jump = 0.0d0;
292 if(testmode .eq. 1) call get_initial_conditions(nnod_loc, u2, u1, v1, nelem_loc, con_spx_loc, con_nnz_loc,&
293 nmat,prop_mat, sdeg_mat, local_node_num, &
294 xx_spx_loc,yy_spx_loc,zz_spx_loc, mpi_id,deltat)
295 u2 = 0.d0
296 else
297 allocate(u0(3*nnod_loc),u1(3*nnod_loc),u2(3*nnod_loc),v1(3*nnod_loc),jump(3*nnod_loc))
298 u0 = 0.0d0; u1 = 0.0d0; u2 = 0.0d0; v1 = 0.0d0; jump = 0.0d0;
299 if(testmode .eq. 1) call get_initial_conditions(nnod_loc, u0, u1, v1, nelem_loc, con_spx_loc, con_nnz_loc,&
300 nmat, prop_mat, sdeg_mat, local_node_num, &
301 xx_spx_loc,yy_spx_loc,zz_spx_loc, mpi_id,deltat)
302
303
304 endif
305
306!---------------------------------------------------------------------------
307! UNKNOWNS INITIALIZATION FOR SSI - AH
308!---------------------------------------------------------------------------
309
310 if (sdofnum.gt.0) then
311 allocate(ug1(3*sdofnum), ug2(3*sdofnum), ug3(3*sdofnum)) !!! base displacement at time n+1,n,n-1
312 ug1=0; ug2=0; ug3=0
313 endif
314
315
316!---------------------------------------------------------------------------
317! UNKNOWNS INITIALIZATION FOR DAMPING
318!---------------------------------------------------------------------------
319
320 if (make_damping_yes_or_not.eq. 1 .and. damping_type .eq. 1) then
321 allocate(mc(3*nnod_loc), mck(3*nnod_loc)); mc = 0.d0; mck = 0.d0
322 elseif(damping_type .eq. 2) then
323 allocate(strain_visc(6*nnod_loc,n_sls)); strain_visc = 0.d0
324 elseif(damping_type .eq. 3) then
325 allocate(mc(3*nnod_loc),damp_term((3*nnod_loc))); mc = 0.d0; damp_term=0.d0;
326 endif
327
328!---------------------------------------------------------------------------
329! UNKNOWNS INITIALIZATION FOR OUTOUT
330!---------------------------------------------------------------------------
331
332 if(opt_out_var(4) .eq. 1) then
333 allocate(stress(6*nnod_loc)); stress = 0.d0
334 endif
335 if(opt_out_var(5) .eq. 1 .or. damping_type .eq. 2) then !Modified in SSI - AH
336 allocate(strain(6*nnod_loc)); strain = 0.d0
337 endif
338 if(opt_out_var(6) .eq. 1) then
339 allocate(omega(3*nnod_loc)); omega = 0.d0
340 endif
341
342!---------------------------------------------------------------------------
343! UNKNOWNS INITIALIZATION FOR DEBUG
344!---------------------------------------------------------------------------
345
346 if(debug .eq. 1) then
347 allocate(mu_nle(nnod_loc), gamma_nle(nnod_loc))
348 mu_nle = 0.d0; gamma_nle = 0.d0
349 endif
350
351!---------------------------------------------------------------------------
352! NODE COUNTER INITIALIZATION FOR OUTPUT
353!---------------------------------------------------------------------------
354
355
356 if(opt_out_var(4) .eq. 1 .or. opt_out_var(5) .eq. 1 .or. opt_out_var(6) .eq. 1 &
357 .or. damping_type .eq. 2) then
358
359 allocate(node_counter(nnod_loc)) ; node_counter = 0
360
361 do ie = 1,nelem_loc
362 im = con_spx_loc(con_spx_loc(ie -1)); nn = sdeg_mat(im) +1
363 do k = 1,nn
364 do j = 1,nn
365 do i = 1,nn
366 is = nn*nn*(k -1) +nn*(j -1) +i; in = con_spx_loc(con_spx_loc(ie -1) + is)
367 node_counter(in) = node_counter(in) + 1
368 enddo
369 enddo
370 enddo
371 enddo
372
373 allocate(small_recv_buffer(nrecv), small_send_buffer(nsend))
374 small_recv_buffer = 0; small_send_buffer = 0
375
376 do i = 1,mpi_np
377 small_send_length(i) = proc_send(i); small_recv_length(i) = proc_recv(i)
378 enddo
379
380 do i = 1,nsend
381 small_send_buffer(i) = node_counter(node_send(i))
382 enddo
383
384 call exchange_integer(nsend,small_send_buffer,nrecv,small_recv_buffer,&
385 mpi_np,small_send_length,small_recv_length,&
386 mpi_comm,mpi_stat,mpi_ierr,mpi_id)
387
388 do i = 1,nrecv
389 node_counter(node_recv(i)) = node_counter(node_recv(i)) + small_recv_buffer(i)
390 enddo
391
392 deallocate(small_recv_buffer, small_send_buffer)
393
394 do i = 1, nnod_loc
395 if(node_counter(i) .eq. 0) write(*,'(A)') 'Error in node counter!!!'
396 enddo
397 endif
398
399 ! if (mpi_id .eq. 0) write(*,*) node_counter
400!---------------------------------------------------------------------------
401! STOP & GO <---> READING BACKUP FILES (If present)
402!---------------------------------------------------------------------------
403
404 if(len_trim(bkp_file) .ne. 70) then
405 file_outu0 = bkp_file(1:len_trim(bkp_file)) // '/U0_'
406 file_outu1 = bkp_file(1:len_trim(bkp_file)) // '/U1_'
407 file_outv1 = bkp_file(1:len_trim(bkp_file)) // '/V1_'
408 file_outxyz = bkp_file(1:len_trim(bkp_file)) // '/XYZ_'
409 if (damping_type .eq. 2 ) then
410 file_outsv = bkp_file(1:len_trim(bkp_file)) // '/SV_'
411 endif
412 else
413 file_outu0 = 'U0_'; file_outu1 = 'U1_'; file_outv1 = 'V1_'
414 file_outxyz = 'XYZ_';
415
416 if (damping_type .eq. 2 ) file_outsv = 'SV_'
417
418 endif
419
420 if (tstart .gt. 0.0d0) then
421 isnap = nint(tstart/(deltat*trestart));
422! tstart = tstart + deltat;
423
424 call read_fileout(file_outu0,isnap,mpi_id,3*nnod_loc,u0)
425 call read_fileout(file_outu1,isnap,mpi_id,3*nnod_loc,u1)
426 call read_fileout(file_outv1,isnap,mpi_id,3*nnod_loc,v1)
427
428 if (damping_type .eq. 2) then
429 call read_fileout_vs(file_outsv,isnap,mpi_id,6*nnod_loc,n_sls,strain_visc)
430 endif
431
432 isnap = isnap + 1;
433
434 else
435 isnap = 1;
436 endif
437
438!---------------------------------------------------------------------------
439! MAKING MASS MATRIX
440!---------------------------------------------------------------------------
441
442 if (nelem_loc.gt.0) then
443
444!$OMP PARALLEL &
445!$OMP PRIVATE(ie,im,nn,ct,ww,dd,mv_el,rho_el,lambda_el,mu_el,gamma_el,k,j,i,is,in,iaz)
446
447!$OMP DO
448 do ie = 1, nelem_loc
449
450 im = con_spx_loc(con_spx_loc(ie - 1)); nn = sdeg_mat(im) + 1
451
452 allocate(ct(nn), ww(nn), dd(nn,nn), mv_el(nn,nn,nn))
453 allocate(rho_el(nn,nn,nn),lambda_el(nn,nn,nn),mu_el(nn,nn,nn),gamma_el(nn,nn,nn))
454 call make_lgl_nw(nn,ct,ww,dd)
455
456 rho_el = prop_mat(im,1); lambda_el = prop_mat(im,2)
457 mu_el = prop_mat(im,3); gamma_el = prop_mat(im,4)
458
459 ! < STANDARD CASE: >
460 ! < if you DON'T go THROUGH the following loop the >
461 ! < mechanical properties are given as usual >
462 ! < BLOCK BY BLOCK >
463
464 if (n_case.gt.0) then
465
466 do icase = 1, n_case
467
468 if (val_case(icase) .eq. tag_mat(im)) then
469 if (num_testcase .ne. 0) check_not_hon = 1;
470
471 call make_eltensor_for_cases(tag_case(icase), val_case(icase),&
472 nn, rho_el, lambda_el, mu_el, gamma_el,&
473 nnod_loc, zs_elev, zs_all, vs_tria, thick, &
474 con_nnz_loc, con_spx_loc, ie,&
475 sub_tag_all, zz_spx_loc, mpi_id, local_node_num, &
476 damping_type, qs_nh, qp_nh, &
477 xx_spx_loc, yy_spx_loc, check_not_hon,label_testcase)
478
479
480 endif
481 enddo
482
483 endif
484
485 !!!!!!!!!!! NOT-HONORING ENHANCED !!!!!!!!!!!!!!!!!!!!!!
486 if (nmat_nhe.gt.0) then
487 qs_nh = qs_nhe_el(ie)
488 qp_nh = qp_nhe_el(ie)
489 call get_mech_prop_nh_enhanced(ie, nn, nnod_loc, con_nnz_loc, con_spx_loc, &
490 rho_nhe, lambda_nhe, mu_nhe, qs_nh, fmax, &
491 rho_el, lambda_el, mu_el, gamma_el)
492 endif
493
495!Random materials begin
497
498 if (nmat_rnd.gt.0) then
499 do irand = 1, nmat_rnd
500 if (rand_mat(irand) .eq. tag_mat(im)) then
501 do k = 1,nn
502 do j = 1,nn
503 do i = 1,nn
504 is = nn*nn*(k -1) +nn*(j -1) +i
505 in = con_spx_loc(con_spx_loc(ie -1) + is)
506 rho_el(i,j,k) = rho_rnd(in);
507 lambda_el(i,j,k) = lambda_rnd(in);
508 mu_el(i,j,k) = mu_rnd(in);
509 enddo
510 enddo
511 enddo
512 endif
513 enddo
514 endif
515
517!Random materials end
519
520
521
522 call make_mass_matrix(nn,ct,ww,dd,rho_el,&
523 alfa11(ie),alfa12(ie),alfa13(ie),&
524 alfa21(ie),alfa22(ie),alfa23(ie),&
525 alfa31(ie),alfa32(ie),alfa33(ie),&
526 beta11(ie),beta12(ie),beta13(ie),&
527 beta21(ie),beta22(ie),beta23(ie),&
528 beta31(ie),beta32(ie),beta33(ie),&
529 gamma1(ie),gamma2(ie),gamma3(ie),&
530 delta1(ie),delta2(ie),delta3(ie),&
531 mv_el)
532
533
534
535 do k = 1,nn
536 do j = 1,nn
537 do i = 1,nn
538 is = nn*nn*(k -1) +nn*(j -1) +i
539 in = con_spx_loc(con_spx_loc(ie -1) + is)
540
541 iaz = 3*(in -1) +1; mv(iaz) = mv(iaz) + mv_el(i,j,k)
542 iaz = 3*(in -1) +2; mv(iaz) = mv(iaz) + mv_el(i,j,k)
543 iaz = 3*(in -1) +3; mv(iaz) = mv(iaz) + mv_el(i,j,k)
544
545 enddo
546 enddo
547 enddo
548
549 deallocate(ct, ww, dd, mv_el, rho_el, lambda_el, mu_el, gamma_el)
550
551 enddo ! loop on elements
552!$OMP END DO
553!$OMP END PARALLEL
554 endif
555
556
557 do i = 1,nrecv
558 in = node_recv(i)
559 recv_buffer(3*(i -1) +1) = mv(3*(in -1) +1)
560 recv_buffer(3*(i -1) +2) = mv(3*(in -1) +2)
561 recv_buffer(3*(i -1) +3) = mv(3*(in -1) +3)
562 enddo
563
564 call exchange_double(3*nrecv,recv_buffer,3*nsend,send_buffer,&
565 mpi_np,recv_length,send_length,&
566 mpi_comm,mpi_stat,mpi_ierr,mpi_id)
567
568 do i = 1,nsend
569 in = node_send(i)
570 mv(3*(in -1) +1) = mv(3*(in -1) +1) +send_buffer(3*(i -1) +1)
571 mv(3*(in -1) +2) = mv(3*(in -1) +2) +send_buffer(3*(i -1) +2)
572 mv(3*(in -1) +3) = mv(3*(in -1) +3) +send_buffer(3*(i -1) +3)
573 enddo
574
575!---------------------------------------------------------------------------
576! MAKING DAMPING MATRIX
577!---------------------------------------------------------------------------
578
579 if (make_damping_yes_or_not .eq. 1 .and. damping_type .ne. 2) then
580 if (nelem_loc .gt. 0) then
581
582!$OMP PARALLEL &
583!$OMP PRIVATE(ie,im,nn,ct,ww,dd,rho_el,lambda_el,mu_el,gamma_el,mc_el,mck_el,i,j,k,is,in,iaz) &
584!$OMP PRIVATE(irand)
585
586!$OMP DO
587 do ie = 1,nelem_loc
588
589 im = con_spx_loc(con_spx_loc(ie -1)); nn = sdeg_mat(im) +1
590
591 allocate(ct(nn),ww(nn),dd(nn,nn))
592 allocate(rho_el(nn,nn,nn),lambda_el(nn,nn,nn),mu_el(nn,nn,nn),gamma_el(nn,nn,nn))
593 call make_lgl_nw(nn,ct,ww,dd)
594
595 rho_el = prop_mat(im,1); lambda_el = prop_mat(im,2)
596 mu_el = prop_mat(im,3); gamma_el = prop_mat(im,4)
597
598 if (damping_type .eq. 1) allocate(mc_el(nn,nn,nn),mck_el(nn,nn,nn))
599 if (damping_type .eq. 3) allocate(mc_el(nn,nn,nn))
600
601
602 ! < STANDARD CASE: >
603 ! < if you DON'T go THROUGH the following loop the >
604 ! < mechanical properties are given as usual >
605 ! < BLOCK BY BLOCK >
606
607 if (n_case.gt.0) then
608
609 do icase = 1, n_case
610
611 if (val_case(icase) .eq. tag_mat(im)) then
612
613
614 call make_eltensor_for_cases(tag_case(icase), val_case(icase),&
615 nn, rho_el, lambda_el, mu_el, gamma_el,&
616 nnod_loc, zs_elev, zs_all, vs_tria, thick, &
617 con_nnz_loc, con_spx_loc, ie,&
618 sub_tag_all, zz_spx_loc, mpi_id, local_node_num, &
619 damping_type, qs_nh, qp_nh, &
620 xx_spx_loc, yy_spx_loc, 0,label_testcase)
621 endif
622 enddo
623
624 endif
625
626 !!!!!!!!!!! NOT-HONORING ENHANCED !!!!!!!!!!!!!!!!!!!!!!
627 if (nmat_nhe.gt.0) then
628 qs_nh = qs_nhe_el(ie)
629 qp_nh = qp_nhe_el(ie)
630 call get_mech_prop_nh_enhanced(ie, nn, nnod_loc, con_nnz_loc, con_spx_loc, &
631 rho_nhe, lambda_nhe, mu_nhe, qs_nh, fmax, &
632 rho_el, lambda_el, mu_el, gamma_el)
633 endif
634
636!Random materials begin
638
639 if (nmat_rnd.gt.0) then
640 do irand = 1, nmat_rnd
641 if (rand_mat(irand) .eq. tag_mat(im)) then
642 do k = 1,nn
643 do j = 1,nn
644 do i = 1,nn
645 is = nn*nn*(k -1) +nn*(j -1) +i
646 in = con_spx_loc(con_spx_loc(ie -1) + is)
647 rho_el(i,j,k) = rho_rnd(in);
648 lambda_el(i,j,k) = lambda_rnd(in);
649 mu_el(i,j,k) = mu_rnd(in);
650 enddo
651 enddo
652 enddo
653 endif
654 enddo
655 endif
656
658!Random materials end
660
661
662 if (damping_type .eq. 1) then
663
664 call make_damping_matrix(nn,ct,ww,dd,&
665 rho_el,gamma_el,&
666 alfa11(ie),alfa12(ie),alfa13(ie),&
667 alfa21(ie),alfa22(ie),alfa23(ie),&
668 alfa31(ie),alfa32(ie),alfa33(ie),&
669 beta11(ie),beta12(ie),beta13(ie),&
670 beta21(ie),beta22(ie),beta23(ie),&
671 beta31(ie),beta32(ie),beta33(ie),&
672 gamma1(ie),gamma2(ie),gamma3(ie),&
673 delta1(ie),delta2(ie),delta3(ie),&
674 mc_el,mck_el)
675 do k = 1,nn
676 do j = 1,nn
677 do i = 1,nn
678 is = nn*nn*(k -1) +nn*(j -1) +i
679 in = con_spx_loc(con_spx_loc(ie -1) + is)
680
681 iaz = 3*(in -1) +1
682 mc(iaz) = mc(iaz) +mc_el(i,j,k)
683 mck(iaz) = mck(iaz) +mck_el(i,j,k)
684
685 iaz = 3*(in -1) +2
686 mc(iaz) = mc(iaz) +mc_el(i,j,k)
687 mck(iaz) = mck(iaz) +mck_el(i,j,k)
688
689 iaz = 3*(in -1) +3
690 mc(iaz) = mc(iaz) +mc_el(i,j,k)
691 mck(iaz) = mck(iaz) +mck_el(i,j,k)
692
693 enddo
694 enddo
695 enddo
696
697 elseif(damping_type .eq. 3) then
698
699 call make_rayleigh_damping_matrix(nn,ct,ww,dd,&
700 rho_el,&
701 alfa11(ie),alfa12(ie),alfa13(ie),&
702 alfa21(ie),alfa22(ie),alfa23(ie),&
703 alfa31(ie),alfa32(ie),alfa33(ie),&
704 beta11(ie),beta12(ie),beta13(ie),&
705 beta21(ie),beta22(ie),beta23(ie),&
706 beta31(ie),beta32(ie),beta33(ie),&
707 gamma1(ie),gamma2(ie),gamma3(ie),&
708 delta1(ie),delta2(ie),delta3(ie),&
709 mc_el)
710
711 do k = 1,nn
712 do j = 1,nn
713 do i = 1,nn
714 is = nn*nn*(k -1) +nn*(j -1) +i
715 in = con_spx_loc(con_spx_loc(ie -1) + is)
716
717 iaz = 3*(in -1) +1; mc(iaz) = mc(iaz) + a0_ray(im)*mc_el(i,j,k);
718 iaz = 3*(in -1) +2; mc(iaz) = mc(iaz) + a0_ray(im)*mc_el(i,j,k);
719 iaz = 3*(in -1) +3; mc(iaz) = mc(iaz) + a0_ray(im)*mc_el(i,j,k);
720
721 enddo
722 enddo
723 enddo
724
725 endif
726
727 deallocate(ct,ww,dd,rho_el,lambda_el,mu_el,gamma_el,mc_el)
728 if (damping_type .eq. 1 ) deallocate(mck_el)
729
730 enddo !loop on lements
731!$OMP END DO
732!$OMP END PARALLEL
733
734 endif
735
736 if (damping_type .eq. 1 .or. damping_type .eq. 3) then
737 do i = 1,nrecv
738 in = node_recv(i)
739 recv_buffer(3*(i -1) +1) = mc(3*(in -1) +1)
740 recv_buffer(3*(i -1) +2) = mc(3*(in -1) +2)
741 recv_buffer(3*(i -1) +3) = mc(3*(in -1) +3)
742 enddo
743
744 call exchange_double(3*nrecv,recv_buffer,3*nsend,send_buffer,&
745 mpi_np,recv_length,send_length,&
746 mpi_comm,mpi_stat,mpi_ierr,mpi_id)
747
748 do i = 1,nsend
749 in = node_send(i)
750 mc(3*(in -1) +1) = mc(3*(in -1) +1) +send_buffer(3*(i -1) +1)
751 mc(3*(in -1) +2) = mc(3*(in -1) +2) +send_buffer(3*(i -1) +2)
752 mc(3*(in -1) +3) = mc(3*(in -1) +3) +send_buffer(3*(i -1) +3)
753 enddo
754
755 endif
756
757 if (damping_type .eq. 1) then
758
759 do i = 1,nrecv
760 in = node_recv(i)
761 recv_buffer(3*(i -1) +1) = mck(3*(in -1) +1)
762 recv_buffer(3*(i -1) +2) = mck(3*(in -1) +2)
763 recv_buffer(3*(i -1) +3) = mck(3*(in -1) +3)
764 enddo
765
766 call exchange_double(3*nrecv,recv_buffer,3*nsend,send_buffer,&
767 mpi_np,recv_length,send_length,&
768 mpi_comm,mpi_stat,mpi_ierr,mpi_id)
769
770 do i = 1,nsend
771 in = node_send(i)
772 mck(3*(in -1) +1) = mck(3*(in -1) +1) +send_buffer(3*(i -1) +1)
773 mck(3*(in -1) +2) = mck(3*(in -1) +2) +send_buffer(3*(i -1) +2)
774 mck(3*(in -1) +3) = mck(3*(in -1) +3) +send_buffer(3*(i -1) +3)
775 enddo
776
777 endif
778!---------------------------------------------------------------------------
779 endif
780
781!---------------------------------------------------------------------------
782! SET VARAIBLES FOR LOW STORAGE RUNGE-KUTTA SCHEME
783
784 if (rk_scheme .eq. 'RUNGEKUTTA') then
785
786 allocate(a(rk_stages), b(rk_stages), c(rk_stages))
787 a=0.d0; b=0.d0; c=0.d0
788 call make_butcherarray(rk_order, rk_stages, a,b,c)
789
790 else !LEAP-FROG
791 allocate(a(1),b(1),c(1))
792 rk_stages=1; rk_order=0;
793 a=0.d0 ; b=0.d0; c=0.d0;
794 endif
795
796 allocate(u_int(3*nnod_loc), v_int(3*nnod_loc)); u_int = 0.d0; v_int = 0.d0;
797
798!---------------------------------------------------------------------------
799! SET VARAIBLES FOR NON-LINEAR ELASTIC FORCES
800! NEW NON LINEAR IMPLEMETATION
801 if (nmat_nle .gt. 0) then
802 allocate(con_spx_loc_nle(0:con_nnz_loc));con_spx_loc_nle = con_spx_loc
803 allocate(depth_nle_el(nelem_loc)); depth_nle_el = 0.d0;
804 allocate(node_nle_4_mpi(1:nnod_loc)); node_nle_4_mpi = 0;
805
806 do ie = 1,nelem_loc
807 im = con_spx_loc(con_spx_loc(ie -1)); nn = sdeg_mat(im) +1
808
809
810 which_mat_nle = 0
811 do im_nle = 1,nmat_nle
812 if (tag_mat(im).eq.tag_mat_nle(im_nle)) then
813 which_mat_nle = im_nle
814 depth_nle = val_mat_nle(im_nle,1)
815 endif
816 enddo
817 con_spx_loc_nle(con_spx_loc_nle(ie -1)) = which_mat_nle
818 depth_nle_el(ie) = depth_nle;
819 if (which_mat_nle .eq. 0) con_spx_loc_nle(con_spx_loc_nle(ie-1)+1:con_spx_loc_nle(ie)-1) = 0;
820 if (which_mat_nle .ne. 0) then
821 do k = 1,nn
822 do j = 1,nn
823 do i = 1,nn
824 is = nn*nn*(k -1) +nn*(j -1) +i
825 in = con_spx_loc(con_spx_loc(ie -1) + is)
826 con_spx_loc_nle(con_spx_loc_nle(ie -1) + is) = 1
827 node_nle_4_mpi(in) = 1;
828 if (n_case .gt. 0) then
829 if (tag_case(1) .ne. 16 .and. tag_case(1) .ne. 21) then
830 if (depth_nle_el(ie) .lt. zs_elev(in) .or. zs_elev(in) .lt. 0.d0) then
831 con_spx_loc_nle(con_spx_loc_nle(ie -1) + is) = 0; node_nle_4_mpi(in) = 0
832 endif
833 if ((tag_case(1) .gt. 0).and.(zs_all(in) .lt. 0.0d0)) then
834 con_spx_loc_nle(con_spx_loc_nle(ie -1) + is) = 0; node_nle_4_mpi(in) = 0
835 endif
836
837 elseif(tag_case(1) .eq. 16) then
838 if (depth_nle_el(ie) .lt. zs_elev(in) .or. zs_elev(in) .lt. 0.d0) then
839 con_spx_loc_nle(con_spx_loc_nle(ie -1) + is) = 0; node_nle_4_mpi(in) = 0
840 endif
841 if (vs_tria(in) .ge. 450.d0) then
842 con_spx_loc_nle(con_spx_loc_nle(ie -1) + is) = 0; node_nle_4_mpi(in) = 0
843 endif
844
845 elseif(tag_case(1) .eq. 21) then
846 if (depth_nle_el(ie) .lt. zs_elev(in) .or. zs_elev(in) .lt. 0.d0) then
847 con_spx_loc_nle(con_spx_loc_nle(ie -1) + is) = 0; node_nle_4_mpi(in) = 0
848 endif
849 if (vs_tria(in) .ge. 600.d0 .or. (zs_all(in) .lt. 0.0d0)) then
850 con_spx_loc_nle(con_spx_loc_nle(ie -1) + is) = 0; node_nle_4_mpi(in) = 0
851 endif
852 endif
853 endif
854 enddo
855 enddo
856 enddo
857 endif
858 enddo
859 endif
860
861! write(*,*) 'nle', con_spx_loc_nle
862
863!---------------------------------------------------------------------------
864
865!---------------------------------------------------------------------------
866! BEGINNING OF THE TIME LOOP
867!---------------------------------------------------------------------------
868
869 itime = 1 !CONVERGENCE TEST
870
871
872 if(mpi_id .eq. 0) write(*,'(A)') '--------------Starting the loop------------------------'
873
874 tt0 = tstart - dfloat(1)*deltat; tt1 = tstart; tt2 = tstart + dfloat(1) * deltat
875 deltat = tt1 - tt0; deltat2 = deltat*deltat;
876 its = 0;
877
878 do its = 0,nts
879 if (opt_out_var(4) .eq. 1) stress = 0.d0
880 if (opt_out_var(5) .eq. 1 .or. damping_type .eq. 2) strain = 0.d0 !Modified in SSI - AH
881 if (opt_out_var(6) .eq. 1) omega = 0.d0
882
883 if (mpi_id.eq.0) write(*,'(A,E14.6)')'TIME = ',tt1
884 if (its.le. 200) then
885 start_time = mpi_wtime()
886 elseif (mod(its, it_divisor) .eq. 0) then
887 ! Print elapsed time since 200-th iteration every it_divisor iterations
888 if (mpi_id .eq. 0) then
889 final_time = mpi_wtime()
890 write(*,'(A5,I9,A16,F12.3)') '[Iter ', its, '] Elapsed time: ', final_time - start_time
891 endif
892 endif
893
894 u_int = 0.d0; v_int = 0.d0;
895
896!---------------------------------------------------------------------------
897! LOOP FOR RK SCHEMES
898!---------------------------------------------------------------------------
899
900 do istage = 1, rk_stages
901
902 tt_int = tt1 + c(istage)*deltat
903
904 fk = 0.0d0; jump = 0.d0; fe = 0.0d0;
905 if (damping_type .eq. 3) damp_term = 0.d0;
906
907
908 do id = 1, nnod_loc
909 do fn = 1,nfunc
910 !if(nload_traX_el .gt. 0) then !TRAVELING POINT LOAD X
911 ! call FIND_POS(nload_traX_el,fun_traX_el,tag_func(fn),find_tag)
912 ! do i = 1, nb_tra_load
913 ! if(local_node_num(id) .eq. node_tra(i) .and. find_tag .ne. 0) then
914 ! iaz = 3*(id -1) +1; fe(iaz) = fe(iaz) + Fel(fn,(3*(id -1) +1)) * &
915 ! GET_FUNC_VALUE(nfunc,func_type,func_indx,func_data,&
916 ! nfunc_data,fn,tt_int,dist_tra(i),val_traX_el(find_tag,1))
917 ! endif
918 ! enddo
919 !elseif(nload_traY_el .gt. 0) then !TRAVELING POINT LOAD X
920 ! call FIND_POS(nload_traY_el,fun_traY_el,tag_func(fn),find_tag)
921 ! do i = 1, nb_tra_load
922 ! if(local_node_num(id) .eq. node_tra(i) .and. find_tag .ne. 0) then
923 ! iaz = 3*(id -1) +2; fe(iaz) = fe(iaz) + Fel(fn,(3*(id -1) +2)) * &
924 ! GET_FUNC_VALUE(nfunc,func_type,func_indx,func_data,&
925 ! nfunc_data,fn,tt_int,dist_tra(i),val_traY_el(find_tag,1))
926 ! endif
927 ! enddo
928 if(nload_traz_el .gt. 0) then !TRAVELING POINT LOAD Z
929 call find_pos(nload_traz_el,fun_traz_el,tag_func(fn),find_tag)
930
931 !write(*,*) 'Inside travelling load', mpi_id, nb_tra_load, find_tag
932
933
934 do i = 1, nb_tra_load
935 ! write(*,*) 'Inside travelling load'
936 ! write(*,*) find_tag, local_node_num(id) , node_tra(i)
937 ! read(*,*)
938
939 if(local_node_num(id) .eq. node_tra(i) .and. find_tag .ne. 0) then
940 iaz = 3*(id -1) +3; fe(iaz) = fe(iaz) + fel(fn,(3*(id -1) +3)) * &
941 get_func_value(nfunc,func_type,func_indx,func_data,&
942 nfunc_data,fn,tt_int,dist_tra(i),val_traz_el(find_tag,1),its)
943 ! val_debug = GET_FUNC_VALUE(nfunc,func_type,func_indx,func_data,&
944 ! nfunc_data,fn,tt_int,dist_tra(i),val_traZ_el(find_tag,1))
945
946 ! write(*,*) 'I am here', local_node_num(id)
947 ! write(*,*) 'input', find_tag
948 endif
949 enddo
950
951
952 else !EXTERNAL FORCES
953 ! Point Source - SS
954 if (fel(fn,(3*(id -1) +1)) .ne. 0.d0) then
955 iaz = 3*(id -1) +1; fe(iaz) = fe(iaz) + fel(fn,(3*(id -1) +1)) * &
956 get_func_value(nfunc,func_type,func_indx,func_data,nfunc_data,fn,tt_int,0,0,its)
957 endif
958 if(fel(fn,(3*(id -1) +2)) .ne. 0.d0) then
959 iaz = 3*(id -1) +2; fe(iaz) = fe(iaz) + fel(fn,(3*(id -1) +2)) * &
960 get_func_value(nfunc,func_type,func_indx,func_data,nfunc_data,fn,tt_int,0,0,its)
961 endif
962 if(fel(fn,(3*(id -1) +3)) .ne. 0.d0) then
963 iaz = 3*(id -1) +3; fe(iaz) = fe(iaz) + fel(fn,(3*(id -1) +3)) * &
964 get_func_value(nfunc,func_type,func_indx,func_data,nfunc_data,fn,tt_int,0,0,its)
965 endif
966 endif
967 enddo
968
969 !! SSI
970 if (sdofnum.gt.0) then
971 if (locnode_buildid_map(id,1).gt.0) then
972 iaz = 3*(id -1);
973 do i=1, locnode_buildid_map(id,1)
974 isdof = locnode_buildid_map(id,i+1)
975 fe(iaz+1) = fe(iaz+1) + (1 /(node_counter_sdof(isdof))) * &
976 sdofforceinput(3*(isdof - 1) + 1)
977
978 fe(iaz+2) = fe(iaz+2) + (1 /node_counter_sdof(isdof)) * &
979 sdofforceinput(3*(isdof - 1) + 2)
980
981 fe(iaz+3) = fe(iaz+3) + (1 /node_counter_sdof(isdof)) * &
982 sdofforceinput(3*(isdof - 1) + 3)
983 enddo
984 endif
985 endif
986
987 enddo
988
989 do i = 1,nrecv
990 in = node_recv(i)
991 recv_buffer(3*(i -1) +1) = fe(3*(in -1) +1)
992 recv_buffer(3*(i -1) +2) = fe(3*(in -1) +2)
993 recv_buffer(3*(i -1) +3) = fe(3*(in -1) +3)
994 enddo
995
996 call exchange_double(3*nrecv,recv_buffer,3*nsend,send_buffer,&
997 mpi_np,recv_length,send_length,&
998 mpi_comm,mpi_stat,mpi_ierr,mpi_id)
999
1000 do i = 1,nsend
1001 in = node_send(i)
1002 fe(3*(in -1) +1) = fe(3*(in -1) +1) +send_buffer(3*(i -1) +1)
1003 fe(3*(in -1) +2) = fe(3*(in -1) +2) +send_buffer(3*(i -1) +2)
1004 fe(3*(in -1) +3) = fe(3*(in -1) +3) +send_buffer(3*(i -1) +3)
1005 enddo
1006
1007! do id = 1, nnod_loc
1008! if (fe(3*(id -1) +1) .ne. 0) then
1009! print*, id, fe(3*(id -1) +1)
1010! elseif (fe(3*(id -1) +2) .ne. 0) then
1011! print*, id, fe(3*(id -1) +2)
1012! elseif (fe(3*(id -1) +3) .ne. 0) then
1013! print*, id, fe(3*(id -1) +3)
1014! endif
1015! enddo
1016! read(*,*)
1017
1018 ! call MPI_BARRIER(mpi_comm, mpi_ierr)
1019 ! stop
1020!---------------------------------------------------------------------------
1021! EXCHANGE U0, U1 & V1 VALUES
1022!---------------------------------------------------------------------------
1023! if (restart .eq. 1) then
1024! do i = 1,nsend
1025! in = node_send(i)
1026! send_buffer(3*(i -1) +1) = u0(3*(in -1) +1)
1027! send_buffer(3*(i -1) +2) = u0(3*(in -1) +2)
1028! send_buffer(3*(i -1) +3) = u0(3*(in -1) +3)
1029! enddo
1030
1031! call EXCHANGE_DOUBLE(3*nsend,send_buffer,3*nrecv,recv_buffer,&
1032! mpi_np,send_length,recv_length,&
1033! mpi_comm,mpi_stat,mpi_ierr,mpi_id)
1034
1035! do i = 1,nrecv
1036! in = node_recv(i)
1037! u0(3*(in -1) +1) = recv_buffer(3*(i -1) +1)
1038! u0(3*(in -1) +2) = recv_buffer(3*(i -1) +2)
1039! u0(3*(in -1) +3) = recv_buffer(3*(i -1) +3)
1040! enddo
1041! endif
1042
1043 do i = 1,nsend
1044 in = node_send(i)
1045 send_buffer(3*(i -1) +1) = u1(3*(in -1) +1)
1046 send_buffer(3*(i -1) +2) = u1(3*(in -1) +2)
1047 send_buffer(3*(i -1) +3) = u1(3*(in -1) +3)
1048 enddo
1049
1050 call exchange_double(3*nsend,send_buffer,3*nrecv,recv_buffer,&
1051 mpi_np,send_length,recv_length,&
1052 mpi_comm,mpi_stat,mpi_ierr,mpi_id)
1053
1054 do i = 1,nrecv
1055 in = node_recv(i)
1056 u1(3*(in -1) +1) = recv_buffer(3*(i -1) +1)
1057 u1(3*(in -1) +2) = recv_buffer(3*(i -1) +2)
1058 u1(3*(in -1) +3) = recv_buffer(3*(i -1) +3)
1059 enddo
1060
1061 do i = 1,nsend
1062 in = node_send(i)
1063 send_buffer(3*(i -1) +1) = v1(3*(in -1) +1)
1064 send_buffer(3*(i -1) +2) = v1(3*(in -1) +2)
1065 send_buffer(3*(i -1) +3) = v1(3*(in -1) +3)
1066 enddo
1067
1068 call exchange_double(3*nsend,send_buffer,3*nrecv,recv_buffer,&
1069 mpi_np,send_length,recv_length,&
1070 mpi_comm,mpi_stat,mpi_ierr,mpi_id)
1071
1072 do i = 1,nrecv
1073 in = node_recv(i)
1074 v1(3*(in -1) +1) = recv_buffer(3*(i -1) +1)
1075 v1(3*(in -1) +2) = recv_buffer(3*(i -1) +2)
1076 v1(3*(in -1) +3) = recv_buffer(3*(i -1) +3)
1077 enddo
1078
1079
1080
1081!---------------------------------------------------------------------------
1082! MAKE INTERNAL FORCES
1083!---------------------------------------------------------------------------
1084
1085
1086
1087 if (nelem_loc.gt.0) then
1088
1089 ! write(*,*) 'nmat_nle', nmat_nle
1090
1091!$OMP PARALLEL &
1092!$OMP PRIVATE(ie,im,nn,ct,ww,dd,ux_el,uy_el,uz_el) &
1093!$OMP PRIVATE(duxdx_el,duydx_el,duzdx_el,duxdy_el,duydy_el,duzdy_el,duxdz_el,duydz_el,duzdz_el) &
1094!$OMP PRIVATE(sxx_el,syy_el,szz_el,syz_el,szx_el,sxy_el,fx_el,fy_el,fz_el,rho_el,lambda_el) &
1095!$OMP PRIVATE(mu_el,gamma_el,mc_el,mck_el,R_el,yes_or_not_nle_el) &
1096!$OMP PRIVATE(im_nle,which_mat_nle,Depth_nle,gamma_new) &
1097!$OMP PRIVATE(k,j,i,is,in,iaz,strain_visc_xx,strain_visc_xy) &
1098!$OMP PRIVATE(strain_visc_xz,strain_visc_yy,strain_visc_yz,strain_visc_zz,QS_nh,QP_nh) &
1099!$OMP PRIVATE(Y_lambda_nh,Y_mu_nh)
1100
1101!$OMP DO
1102
1103 do ie = 1,nelem_loc
1104
1105 im = con_spx_loc(con_spx_loc(ie -1)); nn = sdeg_mat(im) +1
1106
1107 allocate(ct(nn),ww(nn),dd(nn,nn))
1108 allocate(ux_el(nn,nn,nn),uy_el(nn,nn,nn),uz_el(nn,nn,nn))
1109 allocate(duxdx_el(nn,nn,nn),duydx_el(nn,nn,nn),duzdx_el(nn,nn,nn))
1110 allocate(duxdy_el(nn,nn,nn),duydy_el(nn,nn,nn),duzdy_el(nn,nn,nn))
1111 allocate(duxdz_el(nn,nn,nn),duydz_el(nn,nn,nn),duzdz_el(nn,nn,nn))
1112 allocate(sxx_el(nn,nn,nn),syy_el(nn,nn,nn),szz_el(nn,nn,nn))
1113 allocate(syz_el(nn,nn,nn),szx_el(nn,nn,nn),sxy_el(nn,nn,nn))
1114 allocate(fx_el(nn,nn,nn),fy_el(nn,nn,nn),fz_el(nn,nn,nn))
1115 allocate(rho_el(nn,nn,nn),lambda_el(nn,nn,nn),mu_el(nn,nn,nn),gamma_el(nn,nn,nn))
1116 call make_lgl_nw(nn,ct,ww,dd)
1117
1118 if (make_damping_yes_or_not .eq. 1 .and. damping_type .eq. 1) allocate(mc_el(nn,nn,nn),mck_el(nn,nn,nn))
1119
1120 if (damping_type .eq. 2) &
1121 allocate(strain_visc_xx(nn,nn,nn,n_sls),strain_visc_xy(nn,nn,nn,n_sls),&
1122 strain_visc_xz(nn,nn,nn,n_sls),strain_visc_yy(nn,nn,nn,n_sls),&
1123 strain_visc_yz(nn,nn,nn,n_sls),strain_visc_zz(nn,nn,nn,n_sls))
1124
1125 if (nmat_nle .gt. 0) then
1126 allocate(r_el(nn,nn,nn), yes_or_not_nle_el(nn,nn,nn),gamma_new(nn,nn,nn))
1127 r_el = 0; yes_or_not_nle_el = 0; gamma_new = 0;
1128
1129 !write(*,*) 'con', 'ie', con_spx_loc_nle(con_spx_loc_nle(ie -1))
1130 !write(*,*) 'con', 'ie', con_spx_loc_nle(con_spx_loc_nle(ie))
1131 !read(*,*)
1132 if (con_spx_loc_nle(con_spx_loc_nle(ie -1)) .ne. 0) then
1133 do k = 1,nn
1134 do j = 1,nn
1135 do i = 1,nn
1136 is = nn*nn*(k -1) +nn*(j -1) +i
1137 !write(*,*) con_spx_loc_nle(con_spx_loc_nle(ie -1) + is)
1138 !read(*,*)
1139 yes_or_not_nle_el(i,j,k) = con_spx_loc_nle(con_spx_loc_nle(ie -1) + is)
1140 in = con_spx_loc(con_spx_loc(ie -1) + is)
1141
1142 if (yes_or_not_nle_el(i,j,k) .eq. 1) then
1143 iaz = 3*(in -1) +1; mc(iaz) = 0.0d0; mck(iaz) = 0.0d0
1144 iaz = 3*(in -1) +2; mc(iaz) = 0.0d0; mck(iaz) = 0.0d0
1145 iaz = 3*(in -1) +3; mc(iaz) = 0.0d0; mck(iaz) = 0.0d0
1146 endif
1147
1148 if (debug .eq. 1) then
1149 mu_nle(in) = 0.0d0; gamma_nle(in) = 0.0d0
1150 endif
1151 enddo
1152 enddo
1153 enddo
1154 endif
1155 endif
1156 ! write(*,*) nmat_nle,'ie',ie, yes_or_not_nle_el
1157 ! read(*,*)
1158
1159
1160 !|
1161 !+----------------------------------------------------------------------------------------
1162 rho_el = prop_mat(im,1); lambda_el = prop_mat(im,2)
1163 mu_el = prop_mat(im,3); gamma_el = prop_mat(im,4)
1164
1165 do k = 1,nn
1166 do j = 1,nn
1167 do i = 1,nn
1168 is = nn*nn*(k -1) +nn*(j -1) +i
1169 in = con_spx_loc(con_spx_loc(ie -1) + is)
1170
1171 if (damping_type .eq. 3) then
1172
1173 iaz = 3*(in -1) +1; ux_el(i,j,k) = u1(iaz) + a1_ray(im)*v1(iaz)
1174 iaz = 3*(in -1) +2; uy_el(i,j,k) = u1(iaz) + a1_ray(im)*v1(iaz)
1175 iaz = 3*(in -1) +3; uz_el(i,j,k) = u1(iaz) + a1_ray(im)*v1(iaz)
1176 else
1177 iaz = 3*(in -1) +1; ux_el(i,j,k) = u1(iaz)
1178 iaz = 3*(in -1) +2; uy_el(i,j,k) = u1(iaz)
1179 iaz = 3*(in -1) +3; uz_el(i,j,k) = u1(iaz)
1180 endif
1181
1182 if(damping_type .eq. 2) then
1183 strain_visc_xx(i,j,k,:) = strain_visc(6*(in -1) +1,:)
1184 strain_visc_yy(i,j,k,:) = strain_visc(6*(in -1) +2,:)
1185 strain_visc_zz(i,j,k,:) = strain_visc(6*(in -1) +3,:)
1186 strain_visc_xy(i,j,k,:) = strain_visc(6*(in -1) +4,:)
1187 strain_visc_yz(i,j,k,:) = strain_visc(6*(in -1) +5,:)
1188 strain_visc_xz(i,j,k,:) = strain_visc(6*(in -1) +6,:)
1189 endif
1190 enddo
1191 enddo
1192 enddo
1193
1194 call make_strain_tensor(nn,ct,ww,dd,&
1195 alfa11(ie),alfa12(ie),alfa13(ie),&
1196 alfa21(ie),alfa22(ie),alfa23(ie),&
1197 alfa31(ie),alfa32(ie),alfa33(ie),&
1198 beta11(ie),beta12(ie),beta13(ie),&
1199 beta21(ie),beta22(ie),beta23(ie),&
1200 beta31(ie),beta32(ie),beta33(ie),&
1201 gamma1(ie),gamma2(ie),gamma3(ie),&
1202 delta1(ie),delta2(ie),delta3(ie),&
1203 ux_el,uy_el,uz_el,&
1204 duxdx_el,duydx_el,duzdx_el,&
1205 duxdy_el,duydy_el,duzdy_el,&
1206 duxdz_el,duydz_el,duzdz_el)
1207
1208 ! < STANDARD CASE: >
1209 ! < if you DON'T go THROUGH the following loop the >
1210 ! < mechanical properties are given as usual >
1211 ! < BLOCK BY BLOCK >
1212
1213 if (n_case.gt.0) then
1214 do icase = 1, n_case
1215 if (val_case(icase) .eq. tag_mat(im)) then
1216 call make_eltensor_for_cases(tag_case(icase), val_case(icase),&
1217 nn, rho_el, lambda_el, mu_el, gamma_el,&
1218 nnod_loc, zs_elev, zs_all, vs_tria, thick, &
1219 con_nnz_loc, con_spx_loc, ie,&
1220 sub_tag_all, zz_spx_loc, mpi_id, local_node_num, &
1221 damping_type, qs_nh, qp_nh, &
1222 xx_spx_loc, yy_spx_loc, 0,label_testcase)
1223 endif
1224 enddo
1225 endif
1226
1227 !!!!!!!!!!! NOT-HONORING ENHANCED !!!!!!!!!!!!!!!!!!!!!!
1228 if (nmat_nhe.gt.0) then
1229 qs_nh = qs_nhe_el(ie)
1230 qp_nh = qp_nhe_el(ie)
1231 call get_mech_prop_nh_enhanced(ie, nn, nnod_loc, con_nnz_loc, con_spx_loc, &
1232 rho_nhe, lambda_nhe, mu_nhe, qs_nh, fmax, &
1233 rho_el, lambda_el, mu_el, gamma_el)
1234 endif
1235
1236
1238!Random materials begin
1240
1241 if (nmat_rnd.gt.0) then
1242 do irand = 1, nmat_rnd
1243 if (rand_mat(irand) .eq. tag_mat(im)) then
1244 do k = 1,nn
1245 do j = 1,nn
1246 do i = 1,nn
1247 is = nn*nn*(k -1) +nn*(j -1) +i
1248 in = con_spx_loc(con_spx_loc(ie -1) + is)
1249 rho_el(i,j,k) = rho_rnd(in);
1250 lambda_el(i,j,k) = lambda_rnd(in);
1251 mu_el(i,j,k) = mu_rnd(in);
1252 enddo
1253 enddo
1254 enddo
1255 endif
1256 enddo
1257 endif
1258
1260!Random materials end
1262
1263 !+----------------------------------------------------------------------------------------
1264 !|
1265 ! NON LINEAR ELASTIC IMPLEMENTATION:
1266 ! MODIFYING lambda_el, mu_el, gamma_el
1267
1268! if (which_mat_nle.ne.0 .and. damping_type .eq. 1) then
1269 if (nmat_nle .gt. 0 .and. damping_type .eq. 1) then
1270
1271 which_mat_nle = con_spx_loc_nle(con_spx_loc_nle(ie-1))
1272 ! write(*,*) which_mat_nle,ie
1273 ! read(*,*)
1274
1275 if (which_mat_nle .ne. 0) then
1276
1277 ! Invariants of the strain tensor sigma_ij
1278 ! & main strain computations
1279
1280 call make_invariants_and_main_strain(nn,&
1281 duxdx_el,duydx_el,duzdx_el,&
1282 duxdy_el,duydy_el,duzdy_el,&
1283 duxdz_el,duydz_el,duzdz_el,&
1284 r_el,yes_or_not_nle_el)
1285
1286 call make_eltensor_for_cases_nle(prop_mat_nle(which_mat_nle,1),r_el, &
1287 nn,rho_el,lambda_el,mu_el,gamma_el,&
1288 con_nnz_loc,con_spx_loc,ie,&
1289 func_type,func_indx,func_data,nfunc_data,&
1290 nfunc,tt_int,tag_func,yes_or_not_nle_el,tag_case(1),&
1291 nnod_loc,vs_tria)
1292
1293
1294 call make_damping_matrix_nle(nn,ct,ww,dd,&
1295 rho_el,gamma_el,&
1296 alfa11(ie),alfa12(ie),alfa13(ie),&
1297 alfa21(ie),alfa22(ie),alfa23(ie),&
1298 alfa31(ie),alfa32(ie),alfa33(ie),&
1299 beta11(ie),beta12(ie),beta13(ie),&
1300 beta21(ie),beta22(ie),beta23(ie),&
1301 beta31(ie),beta32(ie),beta33(ie),&
1302 gamma1(ie),gamma2(ie),gamma3(ie),&
1303 delta1(ie),delta2(ie),delta3(ie),&
1304 mc_el,mck_el,&
1305 con_nnz_loc,con_spx_loc,ie,&
1306 prop_mat_nle(which_mat_nle,1),r_el,fpeak, &
1307 func_type,func_indx,func_data,nfunc_data,&
1308 nfunc,tt_int,tag_func,yes_or_not_nle_el,gamma_new, &
1309 tag_case(1),nnod_loc,vs_tria)
1310
1311
1312 do k = 1,nn
1313 do j = 1,nn
1314 do i = 1,nn
1315
1316 if (yes_or_not_nle_el(i,j,k) .eq. 1) then
1317
1318 is = nn*nn*(k -1) +nn*(j -1) +i
1319 in = con_spx_loc(con_spx_loc(ie -1) + is)
1320
1321 iaz = 3*(in -1) +1
1322 mc(iaz) = mc(iaz) +mc_el(i,j,k)
1323 mck(iaz) = mck(iaz) +mck_el(i,j,k)
1324 if (debug .eq. 1) mu_nle(in) = mu_el(i,j,k)
1325 if (debug .eq. 1) gamma_nle(in) = gamma_new(i,j,k)
1326
1327 iaz = 3*(in -1) +2
1328 mc(iaz) = mc(iaz) +mc_el(i,j,k)
1329 mck(iaz) = mck(iaz) +mck_el(i,j,k)
1330 iaz = 3*(in -1) +3
1331 mc(iaz) = mc(iaz) + mc_el(i,j,k)
1332 mck(iaz) = mck(iaz) + mck_el(i,j,k)
1333
1334 endif !if (yes_or_not_nle_el(i,j,k).eq.1) then
1335
1336 enddo
1337 enddo
1338 enddo
1339 endif !which_mat_nle .ne. 0
1340 endif !nmat_nle .ne. 0
1341 !|
1342 !+----------------------------------------------------------------------------------------
1343
1344 !+----------------------------------------------------------------------------------------
1345 !|
1346 ! HERE VISCOPLASTIC IMPLEMENTATION:
1347 ! MODIFYING lambda_el and mu_el
1348 !|
1349 !+----------------------------------------------------------------------------------------
1350
1351 if (damping_type .eq. 2) then
1352
1353
1354 if (n_case.gt.0) then
1355 do icase = 1, n_case
1356 if (val_case(icase) .eq. tag_mat(im)) then
1357 call make_anelastic_coefficients_nh(nn, n_sls, lambda_el, mu_el,&
1358 rho_el, qs_nh, qp_nh, fmax, &
1359 y_lambda_nh, y_mu_nh, frequency_range, mpi_id)
1360 endif
1361 enddo
1362 else
1363 y_lambda_nh(1,:) = y_lambda(im,:);
1364 y_mu_nh(1,:) = y_mu(im,:);
1365
1366 endif
1367
1368 if (nmat_nhe.gt.0) then
1369 call make_anelastic_coefficients_nh(nn, n_sls, lambda_el, mu_el,&
1370 rho_el, qs_nh, qp_nh, fmax, &
1371 y_lambda_nh, y_mu_nh, frequency_range, mpi_id)
1372 endif
1373
1374
1375 call make_stress_tensor_damped(nn,lambda_el,mu_el,&
1376 duxdx_el,duydx_el,duzdx_el,&
1377 duxdy_el,duydy_el,duzdy_el,&
1378 duxdz_el,duydz_el,duzdz_el,&
1379 strain_visc_xx, strain_visc_xy, strain_visc_xz, &
1380 strain_visc_yy, strain_visc_yz, strain_visc_zz, &
1381 y_lambda_nh(1,:),y_mu_nh(1,:), n_sls, &
1382 sxx_el,syy_el,szz_el,&
1383 syz_el,szx_el,sxy_el,ie)
1384
1385
1386
1387 do k = 1,nn
1388 do j = 1,nn
1389 do i = 1,nn
1390 is = nn*nn*(k -1) +nn*(j -1) +i
1391 in = con_spx_loc(con_spx_loc(ie -1) + is)
1392
1393 iaz = 6*(in -1) +1; strain(iaz) = strain(iaz) + duxdx_el(i,j,k)/node_counter(in) !xx
1394 !if (local_el_num(ie) .eq. 453) write(*,*) mpi_id, &
1395 ! local_node_num(in), duxdx_el(i,j,k)
1396 iaz = 6*(in -1) +2; strain(iaz) = strain(iaz) + duydy_el(i,j,k)/node_counter(in) !yy
1397 !if (local_el_num(ie) .eq. 453) write(*,*) mpi_id, &
1398 ! local_node_num(in), duydy_el(i,j,k)
1399 iaz = 6*(in -1) +3; strain(iaz) = strain(iaz) + duzdz_el(i,j,k)/node_counter(in) !zz
1400 !if (local_el_num(ie) .eq. 453) write(*,*) mpi_id, &
1401 ! local_node_num(in), duzdz_el(i,j,k)
1402 iaz = 6*(in -1) +4; strain(iaz) = strain(iaz) &
1403 + 0.5*(duxdy_el(i,j,k) + duydx_el(i,j,k))/node_counter(in) !xy
1404 !if (local_el_num(ie) .eq. 453) write(*,*) mpi_id, &
1405 ! local_node_num(in), duxdy_el(i,j,k)
1406 iaz = 6*(in -1) +5; strain(iaz) = strain(iaz) &
1407 + 0.5*(duydz_el(i,j,k) + duzdy_el(i,j,k))/node_counter(in) !yz
1408 !if (local_el_num(ie) .eq. 453) write(*,*) mpi_id, &
1409 ! local_node_num(in), duydz_el(i,j,k)
1410 iaz = 6*(in -1) +6; strain(iaz) = strain(iaz) &
1411 + 0.5*(duxdz_el(i,j,k) + duzdx_el(i,j,k))/node_counter(in) !zx
1412 !if (local_el_num(ie) .eq. 453) write(*,*) mpi_id, &
1413 ! local_node_num(in), duxdz_el(i,j,k)
1414 enddo
1415 enddo
1416 enddo
1417
1418
1419
1420 else
1421 call make_stress_tensor(nn,lambda_el,mu_el,&
1422 duxdx_el,duydx_el,duzdx_el,&
1423 duxdy_el,duydy_el,duzdy_el,&
1424 duxdz_el,duydz_el,duzdz_el,&
1425 sxx_el,syy_el,szz_el,&
1426 syz_el,szx_el,sxy_el)
1427 endif
1428
1429
1430
1431 if (nload_sism_el.gt.0) then
1432
1433 if(length_check_node_sism .gt. 0) then
1434
1435 call make_seismic_moment(nn,ct,ww,dd,&
1436 sxx_el,syy_el,szz_el,syz_el,szx_el,sxy_el,&
1437 nload_sism_el,&
1438 check_node_sism,check_dist_node_sism,&
1439 length_check_node_sism,ie,factor_seismic_moment,&
1440 tau_seismic_moment,&
1441 func_type,func_indx,func_data,nfunc_data,nfunc,tt_int,&
1442 con_nnz_loc,con_spx_loc,tag_func, &
1443 nelem_loc, local_el_num, &
1444 nnod_loc, local_node_num)
1445 endif
1446 endif
1447
1448 if (nload_expl_el.gt.0) then
1449
1450 if(length_check_node_expl .gt. 0) then
1451
1452 call make_expl_source(nn,ct,ww,dd,&
1453 sxx_el,syy_el,szz_el,syz_el,szx_el,sxy_el,&
1454 nload_expl_el,&
1455 check_node_expl,check_dist_node_expl,&
1456 length_check_node_expl,ie,factor_explosive_source,&
1457 func_type,func_indx,func_data,nfunc_data,nfunc,tt_int,&
1458 con_nnz_loc,con_spx_loc,tag_func, &
1459 nelem_loc, local_el_num, &
1460 nnod_loc, local_node_num)
1461 endif
1462 endif
1463
1464 if(opt_out_var(4) .eq. 1) then
1465 do k = 1,nn
1466 do j = 1,nn
1467 do i = 1,nn
1468 is = nn*nn*(k -1) +nn*(j -1) +i
1469 in = con_spx_loc(con_spx_loc(ie -1) + is)
1470 iaz = 6*(in -1) +1; stress(iaz) = stress(iaz) + sxx_el(i,j,k)/node_counter(in)
1471 iaz = 6*(in -1) +2; stress(iaz) = stress(iaz) + syy_el(i,j,k)/node_counter(in)
1472 iaz = 6*(in -1) +3; stress(iaz) = stress(iaz) + szz_el(i,j,k)/node_counter(in)
1473 iaz = 6*(in -1) +4; stress(iaz) = stress(iaz) + sxy_el(i,j,k)/node_counter(in)
1474 iaz = 6*(in -1) +5; stress(iaz) = stress(iaz) + syz_el(i,j,k)/node_counter(in)
1475 iaz = 6*(in -1) +6; stress(iaz) = stress(iaz) + szx_el(i,j,k)/node_counter(in)
1476 enddo
1477 enddo
1478 enddo
1479 endif
1480 if(opt_out_var(5) .eq. 1 .and. damping_type .eq. 1) then
1481 do k = 1,nn
1482 do j = 1,nn
1483 do i = 1,nn
1484 is = nn*nn*(k -1) +nn*(j -1) +i
1485 in = con_spx_loc(con_spx_loc(ie -1) + is)
1486 iaz = 6*(in -1) +1; strain(iaz) = strain(iaz) + duxdx_el(i,j,k) /node_counter(in)
1487 iaz = 6*(in -1) +2; strain(iaz) = strain(iaz) + duydy_el(i,j,k) /node_counter(in)
1488 iaz = 6*(in -1) +3; strain(iaz) = strain(iaz) + duzdz_el(i,j,k) /node_counter(in)
1489 iaz = 6*(in -1) +4; strain(iaz) = strain(iaz) + 0.5*(duxdy_el(i,j,k) + duydx_el(i,j,k)) &
1490 /node_counter(in)
1491 iaz = 6*(in -1) +5; strain(iaz) = strain(iaz) + 0.5*(duydz_el(i,j,k) + duzdy_el(i,j,k)) &
1492 /node_counter(in)
1493 iaz = 6*(in -1) +6; strain(iaz) = strain(iaz) + 0.5*(duxdz_el(i,j,k) + duzdx_el(i,j,k)) &
1494 /node_counter(in)
1495 enddo
1496 enddo
1497 enddo
1498 endif
1499
1500 if(opt_out_var(6) .eq. 1) then
1501 do k = 1,nn
1502 do j = 1,nn
1503 do i = 1,nn
1504 is = nn*nn*(k -1) +nn*(j -1) +i
1505 in = con_spx_loc(con_spx_loc(ie -1) + is)
1506 iaz = 3*(in -1) +1; omega(iaz) = omega(iaz) + 0.5*(duxdy_el(i,j,k) - duydx_el(i,j,k)) &
1507 /node_counter(in)
1508 iaz = 3*(in -1) +2; omega(iaz) = omega(iaz) + 0.5*(duydz_el(i,j,k) - duzdy_el(i,j,k)) &
1509 /node_counter(in)
1510 iaz = 3*(in -1) +3; omega(iaz) = omega(iaz) + 0.5*(duxdz_el(i,j,k) - duzdx_el(i,j,k)) &
1511 /node_counter(in)
1512 enddo
1513 enddo
1514 enddo
1515 endif
1516
1517
1518
1519 call make_internal_force(nn,ct,ww,dd,&
1520 alfa11(ie),alfa12(ie),alfa13(ie),&
1521 alfa21(ie),alfa22(ie),alfa23(ie),&
1522 alfa31(ie),alfa32(ie),alfa33(ie),&
1523 beta11(ie),beta12(ie),beta13(ie),&
1524 beta21(ie),beta22(ie),beta23(ie),&
1525 beta31(ie),beta32(ie),beta33(ie),&
1526 gamma1(ie),gamma2(ie),gamma3(ie),&
1527 delta1(ie),delta2(ie),delta3(ie),&
1528 sxx_el,syy_el,szz_el,syz_el,szx_el,sxy_el,&
1529 fx_el,fy_el,fz_el)
1530
1531
1532
1533 do k = 1,nn
1534 do j = 1,nn
1535 do i = 1,nn
1536 is = nn*nn*(k -1) +nn*(j -1) +i
1537 in = con_spx_loc(con_spx_loc(ie -1) + is)
1538
1539 iaz = 3*(in -1) +1; fk(iaz) = fk(iaz) +fx_el(i,j,k);
1540 iaz = 3*(in -1) +2; fk(iaz) = fk(iaz) +fy_el(i,j,k)
1541 iaz = 3*(in -1) +3; fk(iaz) = fk(iaz) +fz_el(i,j,k)
1542
1543 if (damping_type .eq. 3) then
1544 iaz = 3*(in -1) +1; damp_term(iaz) = damp_term(iaz) + ux_el(i,j,k)
1545 iaz = 3*(in -1) +2; damp_term(iaz) = damp_term(iaz) + uy_el(i,j,k)
1546 iaz = 3*(in -1) +3; damp_term(iaz) = damp_term(iaz) + uz_el(i,j,k)
1547 endif
1548
1549 enddo
1550 enddo
1551 enddo
1552
1553
1554 deallocate(ct,ww,dd,ux_el,uy_el,uz_el)
1555 deallocate(duxdx_el,duydx_el,duzdx_el,duxdy_el,duydy_el,duzdy_el,duxdz_el,duydz_el,duzdz_el)
1556 deallocate(sxx_el,syy_el,szz_el,syz_el,szx_el,sxy_el,fx_el,fy_el,fz_el)
1557 deallocate(rho_el,lambda_el,mu_el,gamma_el)
1558 if(damping_type .eq. 2) deallocate(strain_visc_xx,strain_visc_xy,strain_visc_xz,&
1559 strain_visc_yy,strain_visc_yz,strain_visc_zz)
1560
1561 if (make_damping_yes_or_not .eq. 1 .and. damping_type .eq. 1) deallocate(mc_el, mck_el)
1562 if (nmat_nle .gt. 0) deallocate(r_el,yes_or_not_nle_el,gamma_new)
1563
1564 enddo !loop on elements
1565!$OMP END DO
1566!$OMP END PARALLEL
1567
1568 endif
1569
1570
1571
1572
1573!---------------------------------------------------------------------------
1574! MAKE INTERNAL FORCES FOR ABC
1575!---------------------------------------------------------------------------
1576
1577 if (nelem_abc.gt.0) then
1578
1579!$OMP PARALLEL &
1580!$OMP PRIVATE(ie,ielem, ie_curr, im, nn, ct,ww,dd, ux_el, uy_el, uz_el, fx_el, fy_el, fz_el) &
1581!$OMP PRIVATE(vx_el, vy_el, vz_el, rho_el, lambda_el, mu_el, gamma_el) &
1582!$OMP PRIVATE( k, j, i, is, in, iaz, QS_nh, QP_nh)
1583
1584!$OMP DO
1585
1586 do ie = 1,nelem_abc
1587
1588 ielem = ielem_abc(ie,1)
1589 call get_indloc_from_indglo(local_el_num, nelem_loc, ielem, ie_curr)
1590
1591 im = con_spx_loc(con_spx_loc(ie_curr -1)); nn = sdeg_mat(im) +1
1592
1593 allocate(ct(nn),ww(nn),dd(nn,nn))
1594 allocate(ux_el(nn,nn,nn),uy_el(nn,nn,nn),uz_el(nn,nn,nn))
1595 allocate(fx_el(nn,nn,nn),fy_el(nn,nn,nn),fz_el(nn,nn,nn))
1596 allocate(vx_el(nn,nn,nn),vy_el(nn,nn,nn),vz_el(nn,nn,nn))
1597 allocate(rho_el(nn,nn,nn),lambda_el(nn,nn,nn),mu_el(nn,nn,nn),gamma_el(nn,nn,nn))
1598
1599 call make_lgl_nw(nn,ct,ww,dd)
1600
1601 rho_el = prop_mat(im,1); lambda_el = prop_mat(im,2)
1602 mu_el = prop_mat(im,3); gamma_el = prop_mat(im,4)
1603
1604 do k = 1,nn
1605 do j = 1,nn
1606 do i = 1,nn
1607 is = nn*nn*(k -1) +nn*(j -1) +i
1608 in = con_spx_loc(con_spx_loc(ie_curr -1) + is)
1609
1610 iaz = 3*(in -1) +1
1611 ux_el(i,j,k) = u1(iaz) !u1(iaz)
1612 vx_el(i,j,k) = v1(iaz) !(3.0d0*u1(iaz) - 4.0d0*u0(iaz) + 1.0d0*u_1(iaz)) / (2.0d0*deltat)
1613 iaz = 3*(in -1) +2
1614 uy_el(i,j,k) = u1(iaz) !u1(iaz)
1615 vy_el(i,j,k) = v1(iaz) !(3.0d0*u1(iaz) - 4.0d0*u0(iaz) + 1.0d0*u_1(iaz)) / (2.0d0*deltat)
1616 iaz = 3*(in -1) +3
1617 uz_el(i,j,k) = u1(iaz) !u1(iaz)
1618 vz_el(i,j,k) = v1(iaz) !(3.0d0*u1(iaz) - 4.0d0*u0(iaz) + 1.0d0*u_1(iaz)) / (2.0d0*deltat)
1619
1620 enddo
1621 enddo
1622 enddo
1623
1624 ! < STANDARD CASE: >
1625 ! < if you DON'T go THROUGH the following loop the >
1626 ! < mechanical properties are given as usual >
1627 ! < BLOCK BY BLOCK >
1628
1629 if (n_case.gt.0) then
1630 do icase = 1, n_case
1631 if (val_case(icase) .eq. tag_mat(im)) then
1632 call make_eltensor_for_cases(tag_case(icase), val_case(icase),&
1633 nn, rho_el, lambda_el, mu_el, gamma_el,&
1634 nnod_loc, zs_elev, zs_all, vs_tria, thick, &
1635 con_nnz_loc, con_spx_loc, ie,&
1636 sub_tag_all, zz_spx_loc, mpi_id, local_node_num, &
1637 damping_type, qs_nh, qp_nh, &
1638 xx_spx_loc, yy_spx_loc, 0,label_testcase)
1639 endif
1640 enddo
1641 endif
1642
1643
1644 !!!!!!!!!!! NOT-HONORING ENHANCED !!!!!!!!!!!!!!!!!!!!!!
1645 if (nmat_nhe.gt.0) then
1646 qs_nh = qs_nhe_el(ie)
1647 qp_nh = qp_nhe_el(ie)
1648 call get_mech_prop_nh_enhanced(ie, nn, nnod_loc, con_nnz_loc, con_spx_loc, &
1649 rho_nhe, lambda_nhe, mu_nhe, qs_nh, fmax, &
1650 rho_el, lambda_el, mu_el, gamma_el)
1651 endif
1652
1654!Random materials begin
1656
1657 if (nmat_rnd.gt.0) then
1658 do irand = 1, nmat_rnd
1659 if (rand_mat(irand) .eq. tag_mat(im)) then
1660 do k = 1,nn
1661 do j = 1,nn
1662 do i = 1,nn
1663 is = nn*nn*(k -1) +nn*(j -1) +i
1664 in = con_spx_loc(con_spx_loc(ie -1) + is)
1665 rho_el(i,j,k) = rho_rnd(in);
1666 lambda_el(i,j,k) = lambda_rnd(in);
1667 mu_el(i,j,k) = mu_rnd(in);
1668 enddo
1669 enddo
1670 enddo
1671 endif
1672 enddo
1673 endif
1674
1676!Random materials end
1678
1679
1680
1681 call make_abc_force(nn,ct,ww,dd,&
1682 rho_el,lambda_el,mu_el,&
1683 alfa11(ie_curr),alfa12(ie_curr),alfa13(ie_curr),&
1684 alfa21(ie_curr),alfa22(ie_curr),alfa23(ie_curr),&
1685 alfa31(ie_curr),alfa32(ie_curr),alfa33(ie_curr),&
1686 beta11(ie_curr),beta12(ie_curr),beta13(ie_curr),&
1687 beta21(ie_curr),beta22(ie_curr),beta23(ie_curr),&
1688 beta31(ie_curr),beta32(ie_curr),beta33(ie_curr),&
1689 gamma1(ie_curr),gamma2(ie_curr),gamma3(ie_curr),&
1690 delta1(ie_curr),delta2(ie_curr),delta3(ie_curr),&
1691 ielem_abc(ie,2),ielem_abc(ie,3),&
1692 ielem_abc(ie,4),ielem_abc(ie,5),&
1693 ielem_abc(ie,6),ielem_abc(ie,7),&
1694 ux_el,uy_el,uz_el,vx_el,vy_el,vz_el,&
1695 fx_el,fy_el,fz_el)
1696
1697 do k = 1,nn
1698 do j = 1,nn
1699 do i = 1,nn
1700 is = nn*nn*(k -1) +nn*(j -1) +i
1701 in = con_spx_loc(con_spx_loc(ie_curr -1) + is)
1702
1703 iaz = 3*(in -1) +1; fk(iaz) = fk(iaz) +fx_el(i,j,k)
1704 iaz = 3*(in -1) +2; fk(iaz) = fk(iaz) +fy_el(i,j,k)
1705 iaz = 3*(in -1) +3; fk(iaz) = fk(iaz) +fz_el(i,j,k)
1706 enddo
1707 enddo
1708 enddo
1709
1710 deallocate(ct,ww,dd,ux_el,uy_el,uz_el)
1711 deallocate(rho_el,lambda_el,mu_el,gamma_el)
1712 deallocate(vx_el,vy_el,vz_el)
1713 deallocate(fx_el,fy_el,fz_el)
1714
1715 enddo !loop on materials
1716!$OMP END DO
1717!$OMP END PARALLEL
1718
1719 endif
1720
1721!---------------------------------------------------------------------------
1722! MAKE DG JUMPS
1723!---------------------------------------------------------------------------
1724
1725 if(nelem_dg_glo .gt. 0) then
1726
1727 jump_minus = 0.d0
1728
1729 do i = 1,nsend_jump
1730 in = node_send_jump(i)
1731
1732 if (testmode .eq. 1) then
1733 send_buffer_jump(3*(i -1) +1) = u1(3*(in -1) +1) !+ v1(3*(in -1) +1)
1734 send_buffer_jump(3*(i -1) +2) = u1(3*(in -1) +2) !+ v1(3*(in -1) +2)
1735 send_buffer_jump(3*(i -1) +3) = u1(3*(in -1) +3) !+ v1(3*(in -1) +3)
1736 elseif (damping_type .eq. 3) then
1737 send_buffer_jump(3*(i -1) +1) = damp_term(3*(in -1) +1)
1738 send_buffer_jump(3*(i -1) +2) = damp_term(3*(in -1) +2)
1739 send_buffer_jump(3*(i -1) +3) = damp_term(3*(in -1) +3)
1740 else
1741 send_buffer_jump(3*(i -1) +1) = u1(3*(in -1) +1)
1742 send_buffer_jump(3*(i -1) +2) = u1(3*(in -1) +2)
1743 send_buffer_jump(3*(i -1) +3) = u1(3*(in -1) +3)
1744 endif
1745
1746 enddo
1747
1748 call exchange_double(3*nsend_jump, send_buffer_jump, 3*nrecv_jump, recv_buffer_jump,&
1749 mpi_np, send_length_jump, recv_length_jump,&
1750 mpi_comm, mpi_stat, mpi_ierr,mpi_id)
1751
1752 do i = 1,nrecv_jump
1753 in = node_recv_jump(i)
1754 jump_minus(3*(i -1) +1) = recv_buffer_jump(3*(in -1) +1)
1755 jump_minus(3*(i -1) +2) = recv_buffer_jump(3*(in -1) +2)
1756 jump_minus(3*(i -1) +3) = recv_buffer_jump(3*(in -1) +3)
1757 enddo
1758
1759
1760!$OMP PARALLEL &
1761!$OMP PRIVATE(ie,ielem, nn, ie_curr, u_p, jumpx, jumpy, jumpz, i, j, k, is, in, iaz) &
1762!$OMP PRIVATE(u_m, jshift, ic, iene, imne, mm, trofa, posit, id, iene_curr, im) &
1763!$OMP PRIVATE(jump_all)
1764
1765!$OMP DO
1766
1767 do ie = 1, nelem_dg
1768
1769 ielem = el_new(ie)%ind; nn = el_new(ie)%deg
1770 im = el_new(ie)%mate;
1771 call get_indloc_from_indglo(local_el_num, nelem_loc, ielem, ie_curr)
1772
1773 allocate(u_p(3*nn**3),jumpx(nn**3),jumpy(nn**3),jumpz(nn**3))
1774
1775 u_p = 0.d0; jumpx = 0.d0; jumpy = 0.d0; jumpz = 0.d0
1776
1777 do k = 1,nn
1778 do j = 1,nn
1779 do i = 1,nn
1780 is = nn*nn*(k -1) +nn*(j -1) +i
1781 in = con_spx_loc(con_spx_loc(ie_curr -1) + is)
1782
1783 if (testmode .eq. 1) then
1784 iaz = 3*(in -1) +1; u_p(is) = u1(iaz) !+ v1(iaz)
1785 iaz = 3*(in -1) +2; u_p(nn**3+is) = u1(iaz) !+ v1(iaz)
1786 iaz = 3*(in -1) +3; u_p(2*nn**3 + is) = u1(iaz) !+ v1(iaz)
1787 elseif(damping_type .eq. 3) then
1788 iaz = 3*(in -1) +1; u_p(is) = u1(iaz) + a1_ray(im)*v1(iaz)
1789 iaz = 3*(in -1) +2; u_p(nn**3+is) = u1(iaz) + a1_ray(im)*v1(iaz)
1790 iaz = 3*(in -1) +3; u_p(2*nn**3 + is) = u1(iaz) + a1_ray(im)*v1(iaz)
1791 else
1792 iaz = 3*(in -1) +1; u_p(is) = u1(iaz)
1793 iaz = 3*(in -1) +2; u_p(nn**3+is) = u1(iaz)
1794 iaz = 3*(in -1) +3; u_p(2*nn**3 + is) = u1(iaz)
1795 endif
1796
1797 enddo
1798 enddo
1799 enddo
1800
1801
1802 allocate(u_m(el_new(ie)%nnz_col)); u_m = 0.d0; jshift = 0
1803
1804 do ic = 1, el_new(ie)%num_of_ne
1805
1806 iene = el_new(ie)%el_conf(ic,1)
1807 imne = el_new(ie)%el_conf(ic,0)
1808
1809 mm = 2
1810 do i = 1, nmat
1811 if(tag_mat(i) .eq. imne ) then
1812 mm = sdeg_mat(i) +1
1813 endif
1814 enddo
1815
1816 trofa = 0
1817 if(mpi_np .gt. 1) call check_mpi(con_nnz_dg, con_spx_dg, iene, trofa, posit)
1818
1819 if(trofa .eq. 1) then
1820 do k = 1,mm
1821 do j = 1,mm
1822 do i = 1,mm
1823 is = mm*mm*(k -1) +mm*(j -1) +i
1824 id = con_spx_dg(con_spx_dg(posit-1)+is)
1825
1826 iaz = 3*(id -1) +1; u_m(is + jshift) = jump_minus(iaz)
1827 iaz = 3*(id -1) +2; u_m(mm**3 + is + jshift) = jump_minus(iaz)
1828 iaz = 3*(id -1) +3; u_m(2*(mm**3) + is + jshift) = jump_minus(iaz)
1829 enddo
1830 enddo
1831 enddo
1832 else
1833
1834 call get_indloc_from_indglo(local_el_num, nelem_loc, iene, iene_curr)
1835
1836 do k = 1,mm
1837 do j = 1,mm
1838 do i = 1,mm
1839 is = mm*mm*(k -1) +mm*(j -1) +i
1840 id = con_spx_loc(con_spx_loc(iene_curr-1)+is)
1841
1842 if (testmode .eq. 1) then
1843 iaz = 3*(id -1) +1; u_m(is + jshift) = u1(iaz) !+ v1(iaz)
1844 iaz = 3*(id -1) +2; u_m(mm**3 + is + jshift) = u1(iaz) !+ v1(iaz)
1845 iaz = 3*(id -1) +3; u_m(2*(mm**3) + is + jshift) = u1(iaz) !+ v1(iaz)
1846 elseif (damping_type .eq. 3) then
1847 iaz = 3*(id -1) +1; u_m(is + jshift) = u1(iaz) + a1_ray(imne)*v1(iaz)
1848 iaz = 3*(id -1) +2; u_m(mm**3 + is + jshift) = u1(iaz) + a1_ray(imne)*v1(iaz)
1849 iaz = 3*(id -1) +3; u_m(2*(mm**3) + is + jshift) = u1(iaz) + a1_ray(imne)*v1(iaz)
1850 else
1851 iaz = 3*(id -1) +1; u_m(is + jshift) = u1(iaz)
1852 iaz = 3*(id -1) +2; u_m(mm**3 + is + jshift) = u1(iaz)
1853 iaz = 3*(id -1) +3; u_m(2*(mm**3) + is + jshift) = u1(iaz)
1854 endif
1855 enddo
1856 enddo
1857 enddo
1858 endif
1859
1860 jshift = jshift + 3*mm**3
1861
1862 enddo
1863
1864 allocate(jump_all(3*nn**3))
1865 jump_all = 0.d0
1866
1867
1868 call matmul_sparse(el_new(ie)%matMin, el_new(ie)%nnz_minus, el_new(ie)%JMin, &
1869 el_new(ie)%IMin, jump_all, 3*nn**3, &
1870 u_m, el_new(ie)%nnz_col,1)
1871
1872 jumpx = jumpx + jump_all(1:nn**3)
1873 jumpy = jumpy + jump_all(nn**3+1 : 2*nn**3)
1874 jumpz = jumpz + jump_all(2*nn**3+1 : 3*nn**3)
1875
1876 deallocate(u_m, jump_all)
1877
1878 allocate(jump_all(3*nn**3)); jump_all = 0.d0
1879
1880 call matmul_sparse(el_new(ie)%matPlus, el_new(ie)%nnz_plus,el_new(ie)%JPlus, &
1881 el_new(ie)%IPlus, jump_all, 3*nn**3, &
1882 u_p, 3*nn**3,1)
1883
1884
1885 jumpx = jumpx + jump_all(1:nn**3)
1886 jumpy = jumpy + jump_all(nn**3+1 : 2*nn**3)
1887 jumpz = jumpz + jump_all(2*nn**3+1 : 3*nn**3)
1888
1889 do k = 1,nn
1890 do j = 1,nn
1891 do i = 1,nn
1892 is = nn*nn*(k -1) +nn*(j -1) +i
1893 in = con_spx_loc(con_spx_loc(ie_curr -1) + is)
1894
1895 iaz = 3*(in -1) + 1; jump(iaz) = jump(iaz) + jumpx(is)
1896 iaz = 3*(in -1) + 2; jump(iaz) = jump(iaz) + jumpy(is)
1897 iaz = 3*(in -1) + 3; jump(iaz) = jump(iaz) + jumpz(is)
1898
1899 enddo
1900 enddo
1901 enddo
1902
1903 deallocate(jumpx,jumpy,jumpz)
1904 deallocate(u_p, jump_all)
1905
1906 enddo
1907
1908!$OMP END DO
1909!$OMP END PARALLEL
1910
1911 endif !if (nelem_dg_glo .gt. 0)
1912
1913
1914
1915!---------------------------------------------------------------------------
1916! EXCHANGE DAMPING MATRIX FOR NONLINEAR ELASTICITY
1917!---------------------------------------------------------------------------
1918
1919
1920 if (nmat_nle.gt.0) then
1921
1922 do i = 1,nrecv
1923 in = node_recv(i)
1924 recv_buffer(3*(i -1) +1) = mc(3*(in -1) +1)
1925 recv_buffer(3*(i -1) +2) = mc(3*(in -1) +2)
1926 recv_buffer(3*(i -1) +3) = mc(3*(in -1) +3)
1927 enddo
1928
1929 call exchange_double(3*nrecv,recv_buffer,3*nsend,send_buffer,&
1930 mpi_np,recv_length,send_length,&
1931 mpi_comm,mpi_stat,mpi_ierr,mpi_id)
1932
1933 do i = 1,nsend
1934 in = node_send(i)
1935! y_or_n_node = 1
1936!
1937! if (n_case.gt.0 .and. tag_case .ne. 16) then
1938! if (Depth_nle .lt. zs_elev(in) .or. zs_elev(in) .lt. 0) y_or_n_node = 0
1939! if ((tag_case .gt. 0).and.(zs_all(in) .lt. 0.0d0)) y_or_n_node = 0
1940
1941! elseif(n_case .gt. 0 .and. tag_case .eq. 16) then
1942! if (Depth_nle .lt. zs_elev(in) .or. zs_elev(in) .lt. 0) y_or_n_node = 0
1943! if (vs_tria(in) .ge. 450.d0) y_or_n_node = 0
1944
1945! endif
1946
1947
1948! if(y_or_n_node .eq. 1) then
1949 if(node_nle_4_mpi(in) .eq. 1) then
1950 mc(3*(in -1) +1) = mc(3*(in -1) +1) +send_buffer(3*(i -1) +1)
1951 mc(3*(in -1) +2) = mc(3*(in -1) +2) +send_buffer(3*(i -1) +2)
1952 mc(3*(in -1) +3) = mc(3*(in -1) +3) +send_buffer(3*(i -1) +3)
1953 endif
1954 enddo
1955
1956
1957 do i = 1,nrecv
1958 in = node_recv(i)
1959 recv_buffer(3*(i -1) +1) = mck(3*(in -1) +1)
1960 recv_buffer(3*(i -1) +2) = mck(3*(in -1) +2)
1961 recv_buffer(3*(i -1) +3) = mck(3*(in -1) +3)
1962 enddo
1963
1964 call exchange_double(3*nrecv,recv_buffer,3*nsend,send_buffer,&
1965 mpi_np,recv_length,send_length,&
1966 mpi_comm,mpi_stat,mpi_ierr,mpi_id)
1967
1968 do i = 1,nsend
1969 in = node_send(i)
1970! y_or_n_node = 1
1971! if (n_case.gt.0 .and. tag_case .ne. 16) then
1972! if (Depth_nle .lt. zs_elev(in) .or. zs_elev(in) .lt. 0) y_or_n_node = 0
1973! if ((tag_case .gt. 0).and.(zs_all(in) .lt. 0.0d0)) y_or_n_node = 0
1974!
1975! elseif(n_case .gt. 0 .and. tag_case .eq. 16) then
1976! if (Depth_nle .lt. zs_elev(in) .or. zs_elev(in) .lt. 0) y_or_n_node = 0
1977! if (vs_tria(in) .ge. 450.d0) y_or_n_node = 0
1978!
1979! endif
1980!
1981! if(y_or_n_node .eq. 1) then
1982 if(node_nle_4_mpi(in) .eq. 1) then
1983 mck(3*(in -1) +1) = mck(3*(in -1) +1) +send_buffer(3*(i -1) +1)
1984 mck(3*(in -1) +2) = mck(3*(in -1) +2) +send_buffer(3*(i -1) +2)
1985 mck(3*(in -1) +3) = mck(3*(in -1) +3) +send_buffer(3*(i -1) +3)
1986 endif
1987 enddo
1988
1989 endif !if (nmat_nle.gt.0) then
1990
1991
1992 do i = 1,nrecv
1993 in = node_recv(i)
1994 recv_buffer(3*(i -1) +1) = fk(3*(in -1) +1)
1995 recv_buffer(3*(i -1) +2) = fk(3*(in -1) +2)
1996 recv_buffer(3*(i -1) +3) = fk(3*(in -1) +3)
1997 enddo
1998
1999 call exchange_double(3*nrecv,recv_buffer,3*nsend,send_buffer,&
2000 mpi_np,recv_length,send_length,&
2001 mpi_comm,mpi_stat,mpi_ierr,mpi_id)
2002
2003 do i = 1,nsend
2004 in = node_send(i)
2005 fk(3*(in -1) +1) = fk(3*(in -1) +1) +send_buffer(3*(i -1) +1)
2006 fk(3*(in -1) +2) = fk(3*(in -1) +2) +send_buffer(3*(i -1) +2)
2007 fk(3*(in -1) +3) = fk(3*(in -1) +3) +send_buffer(3*(i -1) +3)
2008 enddo
2009
2010
2011!---------------------------------------------------------------------------
2012! EXCHANGE JUMPS
2013!---------------------------------------------------------------------------
2014
2015 if (nelem_dg_glo.gt.0) then
2016
2017 do i = 1,nrecv
2018 in = node_recv(i)
2019 recv_buffer(3*(i -1) +1) = jump(3*(in -1) +1)
2020 recv_buffer(3*(i -1) +2) = jump(3*(in -1) +2)
2021 recv_buffer(3*(i -1) +3) = jump(3*(in -1) +3)
2022 enddo
2023
2024 call exchange_double(3*nrecv,recv_buffer,3*nsend,send_buffer,&
2025 mpi_np,recv_length,send_length,&
2026 mpi_comm,mpi_stat,mpi_ierr,mpi_id)
2027
2028 do i = 1,nsend
2029 in = node_send(i)
2030 jump(3*(in -1) +1) = jump(3*(in -1) +1) +send_buffer(3*(i -1) +1)
2031 jump(3*(in -1) +2) = jump(3*(in -1) +2) +send_buffer(3*(i -1) +2)
2032 jump(3*(in -1) +3) = jump(3*(in -1) +3) +send_buffer(3*(i -1) +3)
2033 enddo
2034
2035 endif
2036
2037
2038 if (rk_scheme .eq. 'RUNGEKUTTA') then
2039 !Low storage Kappau = dq for u
2040 if (istage .eq. 1) then
2041 u_int = deltat*v1
2042 if(make_damping_yes_or_not.eq. 1 .and. damping_type .eq. 1) &
2043 v_int = -deltat*(fk + jump + mck*u1 + mc*v1 - fe)/mv
2044 if(make_damping_yes_or_not.eq. 1 .and. damping_type .eq. 3) &
2045 v_int = -deltat*(fk + jump + mc*v1 - fe)/mv
2046 if(make_damping_yes_or_not.eq. 0 .or. damping_type .eq. 2) v_int = -deltat*(fk + jump - fe)/mv
2047
2048 else
2049 u_int = a(istage)*u_int + deltat*v1
2050 if(make_damping_yes_or_not.eq. 1 .and. damping_type .eq. 1) &
2051 v_int = a(istage)*v_int - deltat*(fk + jump + mck*u1 + mc*v1 - fe)/mv
2052 if(make_damping_yes_or_not.eq. 1 .and. damping_type .eq. 3) &
2053 v_int = a(istage)*v_int - deltat*(fk + jump + mc*v1 - fe)/mv
2054 if(make_damping_yes_or_not.eq. 0 .or. damping_type .eq. 2) v_int = a(istage)*v_int - deltat*(fk + jump - fe)/mv
2055
2056 endif
2057 u1 = u1 + b(istage)*u_int
2058 v1 = v1 + b(istage)*v_int
2059
2060 endif
2061
2062 enddo !istage
2063
2064
2065 if(rk_scheme .eq. 'RUNGEKUTTA') then
2066 u2 = u1; v2 = v1
2067
2068 else
2069
2070! file_fe0 = "fe0.txt"; unit_fe0 = 1000;
2071! file_fe1 = "fe1.txt"; unit_fe1 = 1001;
2072
2073! if (mpi_id .eq. 0) open(unit_fe0,file=file_fe0,position='APPEND')
2074! if (mpi_id .eq. 1) open(unit_fe1,file=file_fe1,position='APPEND')
2075
2076
2077!---------------------------------------------------------------------------
2078! TIME STEP FOR LEAP FROG
2079!---------------------------------------------------------------------------
2080! if(its .eq. 0 .and. testmode .ne. 1) then
2081
2082!!$OMP PARALLEL DO PRIVATE(iaz,id) !Check this with Ilario
2083
2084! do id = 1, nnod_loc
2085! iaz = 3*(id -1) +1
2086! u1(iaz) = u2(iaz) -v1(iaz)*deltat &
2087! +0.5d0*((fe(iaz) -fk(iaz) - jump(iaz)) /mv(iaz) *deltat2)
2088! u0(iaz) = u1(iaz) -v1(iaz)*deltat
2089
2090! iaz = 3*(id -1) +2
2091! u1(iaz) = u2(iaz) -v1(iaz)*deltat &
2092! +0.5d0*((fe(iaz) -fk(iaz) - jump(iaz)) /mv(iaz) *deltat2)
2093! u0(iaz) = u1(iaz) -v1(iaz)*deltat
2094
2095! iaz = 3*(id -1) +3
2096! u1(iaz) = u2(iaz) -v1(iaz)*deltat &
2097! +0.5d0*((fe(iaz) -fk(iaz) - jump(iaz)) /mv(iaz) *deltat2)
2098! u0(iaz) = u1(iaz) -v1(iaz)*deltat
2099! enddo
2100
2101! else
2102
2103 do id = 1, nnod_loc
2104
2105 iaz = 3*(id -1) +1
2106 if (make_damping_yes_or_not .eq. 0 .or. damping_type .eq. 2) then
2107 u2(iaz) = 2.0d0 * u1(iaz) - u0(iaz) &
2108 + (fe(iaz) - fk(iaz) - jump(iaz)) / mv(iaz) *deltat2
2109 else
2110 u2(iaz) = ( fe(iaz) - fk(iaz) - jump(iaz) - mck(iaz) * u1(iaz) &
2111 + (mc(iaz)/(2*deltat)) * u0(iaz) + (mv(iaz)/deltat2) &
2112 * ( 2.0d0 * u1(iaz) - u0(iaz) ) ) / (mv(iaz)/deltat2 + mc(iaz)/(2*deltat) )
2113
2114 endif
2115 if(abs(u2(iaz)).lt. 1.0e-30) u2(iaz) = 0.d0
2116
2117 iaz = 3*(id -1) +2
2118 if (make_damping_yes_or_not.eq. 0 .or. damping_type .eq. 2) then
2119 u2(iaz) = 2.0d0 * u1(iaz) - u0(iaz) &
2120 + (fe(iaz) - fk(iaz) - jump(iaz)) / mv(iaz) *deltat2
2121 else
2122 u2(iaz) = ( fe(iaz) - fk(iaz) - jump(iaz) - mck(iaz) * u1(iaz) &
2123 + (mc(iaz)/(2*deltat)) * u0(iaz) + (mv(iaz)/deltat2) &
2124 * ( 2.0d0 * u1(iaz) - u0(iaz) ) ) &
2125 / (mv(iaz)/deltat2 + mc(iaz)/(2*deltat) )
2126 endif
2127 if(abs(u2(iaz)).lt. 1.0e-30) u2(iaz) = 0.d0
2128
2129 iaz = 3*(id -1) +3
2130 if (make_damping_yes_or_not.eq. 0 .or. damping_type .eq. 2) then
2131 u2(iaz) = 2.0d0 * u1(iaz) - u0(iaz) &
2132 + (fe(iaz) - fk(iaz) - jump(iaz)) / mv(iaz) *deltat2
2133 else
2134 u2(iaz) = ( fe(iaz) - fk(iaz) - jump(iaz) - mck(iaz) * u1(iaz) &
2135 + (mc(iaz)/(2*deltat)) * u0(iaz) + (mv(iaz)/deltat2) &
2136 * ( 2.0d0 * u1(iaz) - u0(iaz) ) ) / (mv(iaz)/deltat2 + mc(iaz)/(2*deltat) )
2137 endif
2138 if(abs(u2(iaz)).lt. 1.0e-30) u2(iaz) = 0.d0
2139
2140
2141 enddo
2142
2143!$END OMP PARALLEL DO !Check this with Ilario
2144! endif
2145
2146 endif
2147
2148 do fn = 1,nfunc
2149 func_value(fn) = get_func_value(nfunc,func_type,func_indx,func_data,nfunc_data,fn,tt2,0,0,its)
2150! write(*,*) tt2, func_value(fn)
2151 enddo
2152! read(*,*)
2153
2154! if (mpi_id .eq. 0) close(unit_fe0)
2155! if (mpi_id .eq. 1) close(unit_fe1)
2156!---------------------------------------------------------------------------
2157! SET DIRICHLET BOUNDARY CONDITIONS
2158!---------------------------------------------------------------------------
2159
2160 if (nnode_dirx.gt.0) then
2161 do i = 1,nnode_dirx
2162 in = inode_dirx(i)
2163 !write(*,*) in
2164 iaz = 3*(in -1) +1;
2165 u2(iaz) = 0.0d0; u1(iaz) = 0.d0; v1(iaz) = 0.d0
2166 if(rk_scheme .eq. 'RUNGEKUTTA') v2(iaz) = 0.d0
2167 do fn = 1,nfunc
2168 ! Different in SSI-AH
2169 u2(iaz) = u2(iaz) + fel(fn,(3*(in -1) +1)) * func_value(fn)
2170 if(rk_scheme .eq. 'RUNGEKUTTA') v2(iaz) = v2(iaz) + fel(fn,(3*(in -1) +1)) * func_value(fn)
2171 enddo
2172 enddo
2173 endif
2174
2175 if (nnode_diry.gt.0) then
2176 do i = 1,nnode_diry
2177 in = inode_diry(i)
2178 iaz = 3*(in -1) +2;
2179 u2(iaz) = 0.0d0; u1(iaz) = 0.d0; v1(iaz) = 0.d0
2180 if(rk_scheme .eq. 'RUNGEKUTTA') v2(iaz) = 0.d0
2181 do fn = 1,nfunc
2182 u2(iaz) = u2(iaz) + fel(fn,(3*(in -1) +2)) * func_value(fn)
2183 if(rk_scheme .eq. 'RUNGEKUTTA') v2(iaz) = v2(iaz) + fel(fn,(3*(in -1) +2)) * func_value(fn)
2184 enddo
2185 enddo
2186 endif
2187
2188 if (nnode_dirz.gt.0) then
2189 do i = 1,nnode_dirz
2190 in = inode_dirz(i)
2191 iaz = 3*(in -1) +3;
2192 u2(iaz) = 0.0d0; u1(iaz) = 0.d0; v1(iaz) = 0.d0
2193 if(rk_scheme .eq. 'RUNGEKUTTA') v2(iaz) = 0.d0
2194 do fn = 1,nfunc
2195 u2(iaz) = u2(iaz) + fel(fn,(3*(in -1) +3)) * func_value(fn)
2196 if(rk_scheme .eq. 'RUNGEKUTTA') v2(iaz) = v2(iaz) + fel(fn,(3*(in -1) +3)) * func_value(fn)
2197 enddo
2198 enddo
2199 endif
2200
2201
2202!---------------------------------------------------------------------------
2203! EXCHANGE STRESS-STRAIN-ROTATION FOR OUTPUT
2204!---------------------------------------------------------------------------
2205
2206!! OMP DO COULD BE INSERTED HERE
2207
2208 if(opt_out_var(4) .eq. 1) then
2209 allocate(double_send_buffer(6*nsend)); double_send_buffer = 0.0d0
2210 allocate(double_recv_buffer(6*nrecv)); double_recv_buffer = 0.0d0
2211
2212 do i = 1,nsend
2213 in = node_send(i)
2214 double_send_buffer(6*(i -1) +1) = stress(6*(in -1) +1)
2215 double_send_buffer(6*(i -1) +2) = stress(6*(in -1) +2)
2216 double_send_buffer(6*(i -1) +3) = stress(6*(in -1) +3)
2217 double_send_buffer(6*(i -1) +4) = stress(6*(in -1) +4)
2218 double_send_buffer(6*(i -1) +5) = stress(6*(in -1) +5)
2219 double_send_buffer(6*(i -1) +6) = stress(6*(in -1) +6)
2220 enddo
2221
2222 call exchange_double(6*nsend,double_send_buffer,6*nrecv,double_recv_buffer,&
2223 mpi_np,double_send_length,double_recv_length,&
2224 mpi_comm,mpi_stat,mpi_ierr,mpi_id)
2225
2226 do i = 1,nrecv
2227 in = node_recv(i)
2228 stress(6*(in -1) +1) = double_recv_buffer(6*(i -1) +1)
2229 stress(6*(in -1) +2) = double_recv_buffer(6*(i -1) +2)
2230 stress(6*(in -1) +3) = double_recv_buffer(6*(i -1) +3)
2231 stress(6*(in -1) +4) = double_recv_buffer(6*(i -1) +4)
2232 stress(6*(in -1) +5) = double_recv_buffer(6*(i -1) +5)
2233 stress(6*(in -1) +6) = double_recv_buffer(6*(i -1) +6)
2234 enddo
2235
2236 deallocate(double_send_buffer,double_recv_buffer)
2237
2238 endif
2239
2240
2241 if(opt_out_var(5) .eq. 1 .or. damping_type .eq. 2) then
2242
2243 allocate(double_send_buffer(6*nsend)); double_send_buffer = 0.0d0
2244 allocate(double_recv_buffer(6*nrecv)); double_recv_buffer = 0.0d0
2245
2246 do i = 1,nsend
2247 in = node_send(i)
2248 double_send_buffer(6*(i -1) +1) = strain(6*(in -1) +1)
2249 double_send_buffer(6*(i -1) +2) = strain(6*(in -1) +2)
2250 double_send_buffer(6*(i -1) +3) = strain(6*(in -1) +3)
2251 double_send_buffer(6*(i -1) +4) = strain(6*(in -1) +4)
2252 double_send_buffer(6*(i -1) +5) = strain(6*(in -1) +5)
2253 double_send_buffer(6*(i -1) +6) = strain(6*(in -1) +6)
2254 enddo
2255
2256 call exchange_double(6*nsend,double_send_buffer,6*nrecv,double_recv_buffer,&
2257 mpi_np,double_send_length,double_recv_length,&
2258 mpi_comm,mpi_stat,mpi_ierr,mpi_id)
2259
2260 do i = 1,nrecv
2261 in = node_recv(i)
2262 strain(6*(in -1) +1) = double_recv_buffer(6*(i -1) +1)
2263 strain(6*(in -1) +2) = double_recv_buffer(6*(i -1) +2)
2264 strain(6*(in -1) +3) = double_recv_buffer(6*(i -1) +3)
2265 strain(6*(in -1) +4) = double_recv_buffer(6*(i -1) +4)
2266 strain(6*(in -1) +5) = double_recv_buffer(6*(i -1) +5)
2267 strain(6*(in -1) +6) = double_recv_buffer(6*(i -1) +6)
2268 enddo
2269
2270
2271 deallocate(double_send_buffer,double_recv_buffer)
2272
2273 endif
2274
2275 if(opt_out_var(6) .eq. 1) then
2276
2277 do i = 1,nsend
2278 in = node_send(i)
2279 send_buffer(3*(i -1) +1) = omega(3*(in -1) +1)
2280 send_buffer(3*(i -1) +2) = omega(3*(in -1) +2)
2281 send_buffer(3*(i -1) +3) = omega(3*(in -1) +3)
2282 enddo
2283
2284 call exchange_double(3*nsend,send_buffer,3*nrecv,recv_buffer,&
2285 mpi_np,send_length,recv_length,&
2286 mpi_comm,mpi_stat,mpi_ierr,mpi_id)
2287
2288 do i = 1,nrecv
2289 in = node_recv(i)
2290 omega(3*(in -1) +1) = recv_buffer(3*(i -1) +1)
2291 omega(3*(in -1) +2) = recv_buffer(3*(i -1) +2)
2292 omega(3*(in -1) +3) = recv_buffer(3*(i -1) +3)
2293
2294 enddo
2295 endif
2296
2297! endif
2298
2299!---------------------------------------------------------------------------
2300! UPDATE DAMPING TENSOR
2301!---------------------------------------------------------------------------
2302
2303 if (damping_type .eq. 2) then
2304 ! Different for SSI-AH
2305 call update_damping_tensor(nnod_loc, n_sls, frequency_range, &
2306 strain_visc, strain, deltat)
2307 endif
2308
2309!---------------------------------------------------------------------------
2310! WRITE MONITORS IN OUTPUT
2311!---------------------------------------------------------------------------
2312
2313 if (nmonitors_lst .ge. 1) then
2314
2315 if (mod(its,ndt_mon_lst) .eq. 0 .and. count_monitor .gt. 0 ) then
2316
2317
2318
2319 call write_output(nmonitors_lst, mpi_id, el_monitor_lst, local_el_num, nelem_loc, &
2320 con_spx_loc, con_nnz_loc, sdeg_mat, nmat, &
2321 u2, u1, u0, v1, nnod_loc, &
2322 stress, strain, omega, &
2323 xr_monitor_lst, yr_monitor_lst, zr_monitor_lst, &
2324 tt1, deltat, deltat2, &
2325 opt_out_var, count_monitor, monitor_file,debug, &
2326 mu_nle, gamma_nle, b_instabilitycontrol, instability_maxval, b_instability_abort)
2327
2328 endif
2329
2330 endif
2331
2332!---------------------------------------------------------------------------
2333! WRITE VTK files
2334!---------------------------------------------------------------------------
2335 ! if ((tt_int.ge.0.8).and.(tt_int.le.2.3)) then
2336 ! if (mod(its,ndt_mon_lst) .eq. 0) then
2337 ! call WRITE_VTK_VOLUME_TIMESERIES(its, mpi_id, nnod_loc, local_node_num, nmat, sdeg_mat, &
2338 ! xx_spx_loc, yy_spx_loc, zz_spx_loc, con_nnz_loc, con_spx_loc, u2, stress)
2339 ! endif
2340 ! endif
2341
2342 !---------------------------------------------------------------------------
2343 ! COMPUTE INPUT OF SDOF SYSTEM
2344 !---------------------------------------------------------------------------
2345
2346 if (sdofnum.gt.0) then
2347 sdofinput = 0
2348 sdofinputd = 0
2349 call compute_sdof_input(sdofnum, mpi_id, el_system_lst, local_el_num, nelem_loc, &
2350 con_spx_loc, con_nnz_loc, sdeg_mat, nmat, &
2351 u2, nnod_loc, &
2352 xr_system_lst, yr_system_lst, zr_system_lst, &
2353 deltat2, &
2354 sdofinput, sdofinputd, mpi_np, &
2355 ug1, ug2, ug3)
2356
2357 !!! ug1, ug2, ug3 : displacements at the base of the structure at times n+1(current), n and n-1 respectively;
2358 !!! SDOFinput : base acceleration
2359 endif
2360
2361 !---------------------------------------------------------------------------
2362 ! exchange ground motion of sdof system
2363 !---------------------------------------------------------------------------
2364
2365 if (sdofnum.gt.0) then
2366 ! base acceleration
2367 sdofrecv_temp = 0;
2368 call mpi_barrier(mpi_comm, mpi_ierr)
2369 call mpi_allreduce(sdofinput, sdofrecv_temp, 3*sdofnum, speed_double, mpi_sum, mpi_comm, mpi_ierr)
2370 sdofinput = sdofrecv_temp;
2371
2372 ! base displacement
2373 sdofrecv_temp = 0;
2374 call mpi_allreduce(sdofinputd, sdofrecv_temp, 3*sdofnum, speed_double, mpi_sum, mpi_comm, mpi_ierr)
2375 sdofinputd = sdofrecv_temp;
2376
2377 ! !Debug - SSI - SS
2378 ! if (mpi_id.eq.0) then
2379 ! do isdof=1,SDOFnum
2380 ! file_getval_id = 2433 + isdof
2381 ! write(file_getval_name,'(A12,I4.4,A4)') 'GETVAL_SDOF_', isdof, '.txt'
2382 ! open(file_getval_id, file=file_getval_name, position='append')
2383 ! !dum_ind = 3*(isdof - 1) + 1
2384 ! write(file_getval_id,*) tt_int, SDOFforceinput(1), SDOFinput(1), SDOFinputD(1)
2385 ! close(file_getval_id)
2386 ! enddo
2387 ! endif
2388 endif
2389
2390
2391 !write(*,*) 'MPI ID :', mpi_id, ', gd (after):', SDOFinputD(1:3), 'ga (after): ', SDOFinput(1:3)
2392
2393 !---------------------------------------------------------------------------
2394 ! compute output of sdof/mdof system(distribute this calculations among mpi_processess in next versions)
2395 !---------------------------------------------------------------------------
2396 if (sdofnum.gt.0) sdofforceinput=0
2397 ! This is only performed in process with mpi_id = 0
2398 if(n_bld .gt. 0) then
2399 do i=1,n_bld !!! number of oscillators
2400
2401 ! Rotating the input ground acceleration into principal coordinates of the building
2402 sdofag(i,1) = rot_cos(i)*sdofinput(3*(i-1)+1) + rot_sin(i)*sdofinput(3*(i-1)+2)
2403 sdofag(i,2) = -rot_sin(i)*sdofinput(3*(i-1)+1) + rot_cos(i)*sdofinput(3*(i-1)+2)
2404 sdofag(i,3) = sdofinput(3*(i-1)+3)
2405 gr_acc_rot = -1*sdofag(i,1:3)
2406
2407 sdofgd(i,1) = rot_cos(i)*sdofinputd(3*(i-1)+1) + rot_sin(i)*sdofinputd(3*(i-1)+2)
2408 sdofgd(i,2) = -rot_sin(i)*sdofinputd(3*(i-1)+1) + rot_cos(i)*sdofinputd(3*(i-1)+2)
2409 sdofgd(i,3) = sdofinputd(3*(i-1)+3)
2410
2411 do j=1,sys(i)%ndt
2412 if(sys(i)%SFS.eq.0) then
2413 call sdof_shear_model(i, sys(i)%NDOF, gr_acc_rot(1),1) !!! direction x
2414 call sdof_shear_model(i, sys(i)%NDOF, gr_acc_rot(2),2) !!! direction y
2415 !call SDOF_SHEAR_MODEL(I,gr_acc_rot(3),3) !!! direction z
2416 elseif (sys(i)%SFS.eq.1) then
2417 call sdof_sfs_model(i,gr_acc_rot(:),1) !!! direction x
2418 call sdof_sfs_model(i,gr_acc_rot(:),2) !!! direction y
2419 endif
2420 enddo
2421
2422 !Rotating the base reactions from building axes on to the global X,Y axes
2423 gr_intf(1) = rot_cos(i)*sys(i)%IntForce(1,1) - rot_sin(i)*sys(i)%IntForce(1,2)
2424 gr_intf(2) = rot_sin(i)*sys(i)%IntForce(1,1) + rot_cos(i)*sys(i)%IntForce(1,2)
2425 gr_intf(3) = sys(i)%IntForce(1,3)
2426 sdofforceinput((3*(i-1)+1):(3*(i-1)+3))=gr_intf(1:3) !!! Interaction force (Ku)
2427 enddo
2428
2429 if(mod(its,ndt_mon_lst) .eq. 0) then
2430 call write_sdof_output_files(tt1)
2431 endif
2432 endif !n_bld.gt.0
2433
2434 !---------------------------------------------------------------------------
2435 ! exchange force output of sdof/mdof system
2436 !---------------------------------------------------------------------------
2437
2438 if (sdofnum.gt.0) then
2439
2440 ! Shear Force in the Spring need to modify this MPI_ALLREDUCE, when Structural calculations are done in multiple MPI processors
2441 sdofrecv_temp = 0;
2442 call mpi_barrier(mpi_comm, mpi_ierr)
2443 call mpi_allreduce(sdofforceinput, sdofrecv_temp, 3*sdofnum, speed_double, mpi_sum, mpi_comm, mpi_ierr)
2444 sdofforceinput = sdofrecv_temp;
2445
2446 endif
2447
2448
2449!---------------------- EXIT Incase of Instability -----------------------------------------------------------
2450 if (b_instabilitycontrol) then
2451 if (b_instability_abort) then
2452 write(*,*) 'Worker ', mpi_id, ': instability detected. Signalling to other workers.'
2453 endif
2454
2455 ! Send & read all workers statuses
2456 call mpi_allgather(b_instability_abort, 1, mpi_logical, b_instability_abort_all, 1, mpi_logical, mpi_comm, mpi_ierr)
2457
2458 ! Abort if one worker fails
2459 if (any(b_instability_abort_all)) then
2460 if (mpi_id .eq. 0) then
2461 write(*,*) 'Instability detected, aborting.'
2462 endif
2463
2464 final_time = mpi_wtime()
2465 call mpi_finalize(mpi_ierr)
2466
2467 if (mpi_id .eq. 0) then
2468 write(*,*) 'Jobs aborted. '
2469 write(*,*) 'Elapsed wall clock time: ', final_time - start_time, ' seconds'
2470 write(*,*) 'Instability detected at iteration = ', its, ', TIME =', tt1
2471 endif
2472 call exit(exit_instab)
2473
2474 endif
2475
2476 endif
2477
2478
2479!---------------------------------------------------------------------------
2480! SETUP MAX VALUES FOR PEAK GROUND MAP
2481!---------------------------------------------------------------------------
2482
2483
2484! if (nmonitors_pgm.ge.1 .and. mod(its,ndt_mon_pgm).eq.0 .and. time_est .le. tol) then
2485!
2486! call GET_MAX_VALUES(nmonitors_pgm, el_monitor_pgm, local_el_num, nelem_loc, &
2487! sdeg_mat, nmat, con_spx_loc, con_nnz_loc, u1, u0, v1,&
2488! omega, nnod_loc, &
2489! xr_monitor_pgm,yr_monitor_pgm,zr_monitor_pgm, &
2490! max_u, max_v, max_a, max_o)
2491!
2492!
2493! endif
2494!---------------------------------------------------------------------------
2495! UPDATE VARIABLES
2496!---------------------------------------------------------------------------
2497
2498 if (rk_scheme .eq. 'RUNGEKUTTA') then
2499 v1 = v2
2500 u1 = u2
2501 else
2502 v1 = (u2-u0)/(2.0d0*deltat)
2503 u0 = u1
2504 u1 = u2
2505 endif
2506
2507 if (its .le. 200) then
2508 call mpi_barrier(mpi_comm, mpi_ierr)
2509 final_time = mpi_wtime()
2510 if (mpi_id .eq. 0) write(*,*) 'Time for a single step', final_time - start_time
2511 endif
2512
2513
2514
2515!---------------------------------------------------------------------------
2516! WRITE BACKUP FILES
2517!---------------------------------------------------------------------------
2518
2519 if(restart .eq. 1 .and. (dabs(tt2 - trestart*deltat*isnap) .le. 1.e-6)) then
2520
2521! write(*,*) dabs(tt2 - restart*deltat*isnap), tt2, restart, deltat, isnap
2522! read(*,*)
2523
2524 if (mpi_id.eq.0) write(*,*) 'Writing Backup files at TIME : ', tt2
2525
2526
2527 call write_fileout_grid(file_outu0,file_outxyz,isnap,mpi_id,3*nnod_loc,u0,&
2528 xx_spx_loc,yy_spx_loc,zz_spx_loc,&
2529 local_node_num,tstart)
2530
2531 call write_fileout_grid(file_outu1,file_outxyz,isnap,mpi_id,3*nnod_loc,u1,&
2532 xx_spx_loc,yy_spx_loc,zz_spx_loc,&
2533 local_node_num,tstart)
2534
2535 call write_fileout_grid(file_outv1,file_outxyz,isnap,mpi_id,3*nnod_loc,v1,&
2536 xx_spx_loc,yy_spx_loc,zz_spx_loc,&
2537 local_node_num,tstart)
2538
2539 if (damping_type .eq. 2) then
2540
2541 call write_fileout_grid_sv(file_outsv,file_outxyz,isnap,mpi_id,6*nnod_loc,n_sls,strain_visc,&
2542 xx_spx_loc,yy_spx_loc,zz_spx_loc,&
2543 local_node_num,tstart)
2544 endif
2545
2546
2547
2548 isnap = isnap + 1
2549
2550 endif
2551
2552
2553
2554!---------------------------------------------------------------------------
2555! VTK FILES
2556!---------------------------------------------------------------------------
2557 ! if (its.eq.0) then
2558 ! ! For Now, this only doesnot work for DG, or mesh with different Spectral Degree
2559 ! nn = sdeg_mat(con_spx_loc(con_spx_loc(0))) + 1
2560 ! allocate(vtk_numbering(nn*nn*nn))
2561
2562 ! call VTK_NODE_NUMBERING_MAP(nn, vtk_numbering)
2563
2564 ! if (mpi_id.eq.0) call WRITE_PVD_TIMEDATA(mpi_np, nts, 13*ndt_mon_lst,deltat)
2565
2566 ! call WRITE_VTK_MECH_PROP(nnod_loc, con_nnz_loc, con_spx_loc, nmat, sdeg_mat, &
2567 ! prop_mat,tag_mat, nmat_nlp, tag_mat_nlp, &
2568 ! xx_spx_loc, yy_spx_loc, zz_spx_loc, mpi_id , nn, vtk_numbering)
2569 ! call EXIT()
2570 ! endif
2571
2572 ! if (mod(its,13*ndt_mon_lst) .eq. 0) then
2573 ! call WRITE_VTU_TIMEDATA(nnod_loc, con_nnz_loc, con_spx_loc, &
2574 ! nmat, sdeg_mat, prop_mat,tag_mat, &
2575 ! xx_spx_loc, yy_spx_loc, zz_spx_loc, mpi_id, mpi_np, &
2576 ! nn, vtk_numbering, &
2577 ! its, strain, stress, u1)
2578 ! endif
2579
2580!---------------------------------------------------------------------------
2581! COMPUTE ERROR
2582!---------------------------------------------------------------------------
2583 if (itime .le. ntime_err) then
2584
2585 if (testmode .eq. 1 .and. abs(time_error(itime)-tt1) .le. deltat) then
2586
2587 call compute_energy_error(nnod_loc, u1, v1, tt1, nelem_loc, con_spx_loc, con_nnz_loc,&
2588 nmat, prop_mat, sdeg_mat, tag_mat, &
2589 local_node_num, local_el_num, nelem_dg,&
2590 alfa11,alfa12,alfa13,&
2591 alfa21,alfa22,alfa23,&
2592 alfa31,alfa32,alfa33,&
2593 beta11,beta12,beta13,&
2594 beta21,beta22,beta23,&
2595 beta31,beta32,beta33,&
2596 gamma1,gamma2,gamma3,&
2597 delta1,delta2,delta3,&
2598 xx_spx_loc, yy_spx_loc, zz_spx_loc, mpi_id, el_new, &
2599 nelem_dg_glo, &
2600 nsend_jump,node_send_jump, &
2601 nrecv_jump,node_recv_jump, &
2602 send_length_jump, recv_length_jump, &
2603 mpi_np, mpi_comm, &
2604 con_nnz_dg, con_spx_dg)
2605
2606 itime = itime + 1
2607
2608 endif
2609 endif
2610
2611
2612 tt0 = tt0 + deltat; tt1 = tt1 + deltat; tt2 = tt2 + deltat;
2613
2614 enddo !loop on its
2615
2616!---------------------------------------------------------------------------
2617! END TIME LOOP
2618!---------------------------------------------------------------------------
2619
2620 if(mpi_id .eq. 0) then
2621 finish = mpi_wtime(); time_in_seconds = finish - start
2622 endif
2623
2624 if (mpi_id.eq.0) then
2625 write(*,'(A)')
2626 write(*,'(A,F20.3,A)')'Time loop completed in ',time_in_seconds,' s'
2627 write(*,'(A)')'-------------------------------------------------------'
2628 write(*,'(A)')
2629 endif
2630
2631!---------------------------------------------------------------------------
2632! WRITING PEAK GROUD MAP
2633!---------------------------------------------------------------------------
2634
2635 if (nmonitors_pgm.ge.1) then
2636
2637 if(mpi_id.eq.0) then
2638 write(*,'(A)')
2639 write(*,'(A)')'--------------Writing Peak Ground Maps-----------------'
2640 write(*,'(A)')
2641 endif
2642
2643 filepg='PGD'
2644 call write_fileout_pg(monitor_file,filepg,mpi_id,nmonitors_pgm,max_u,9)
2645 filepg='PGV'
2646 call write_fileout_pg(monitor_file,filepg,mpi_id,nmonitors_pgm,max_v,9)
2647 filepg='PGA'
2648 call write_fileout_pg(monitor_file,filepg,mpi_id,nmonitors_pgm,max_a,9)
2649 filepg='PGO'
2650 call write_fileout_pg(monitor_file,filepg,mpi_id,nmonitors_pgm,max_o,3)
2651 endif
2652
2653
2654!---------------------------------------------------------------------------
2655! DEALLOCATION
2656!---------------------------------------------------------------------------
2657
2658 if (debug .eq. 1) deallocate(mu_nle,gamma_nle)
2659 if (damping_type .eq. 2) deallocate(strain_visc)
2660 if (damping_type .eq. 3) deallocate(damp_term)
2661
2662
2663 deallocate(func_value, u1, u2, v1, fe, fk, mv, jump,send_buffer,recv_buffer,v_int,u_int)
2664 if (rk_scheme .ne. 'RUNGEKUTTA') deallocate(u0)
2665
2666 if (nelem_dg_glo .gt. 0) deallocate(send_buffer_jump, recv_buffer_jump, jump_minus)
2667 if (opt_out_var(4) .eq. 1 .or. opt_out_var(5) .eq. 1 .or. opt_out_var(6) .eq. 1) then
2668 deallocate(node_counter)
2669 endif
2670 if (rk_scheme .eq. 'RUNGEKUTTA') deallocate(v2, a,b,c)
2671 if (nmonitors_pgm .ge. 1) deallocate(max_u, max_v, max_a, max_o)
2672 if (nmat_nle.gt.0) deallocate(con_spx_loc_nle,depth_nle_el,node_nle_4_mpi)
2673
2674 if(nload_traz_el .gt. 0) deallocate(node_tra, dist_tra)
2675
2676 !!! AH
2677 if(n_bld.gt.0) then
2678 deallocate(sys)
2679 deallocate(sdofag,sdofgd)
2680 endif
2681
2682 if(sdofnum.gt.0) then
2683 deallocate(sdofinput, sdofinputd, sdofforceinput)
2684 deallocate(sdofinputbuffer, sdofforceinputbuffer)
2685 deallocate(ug1, ug2, ug3, sdofrecv_temp)
2686 endif
2687 !!! AH
2688
2689 !ADAPTIVE TIME STEP
2690 !if (time_adpt .eq. 1) deallocate(fe_1,fe_2,fe_3,fext,fext_1,fext_2,fext_3,fk_1,fk_2,eta_1,eta,dg,d2g)
2691
2692
2693 end subroutine time_loop
check_mpi
subroutine check_mpi(n, v, ind, tt, pos)
Checks if a an element is present on a vector and give its position.
Definition CHECK_MPI.f90:30
compute_energy_error
subroutine compute_energy_error(nnod_loc, u1, v1, time, ne_loc, cs_loc, cs_nnz_loc, nm, prop_mat, sdeg_mat, tag_mat, loc_n_num, local_el_num, nelem_dg_local, alfa11, alfa12, alfa13, alfa21, alfa22, alfa23, alfa31, alfa32, alfa33, beta11, beta12, beta13, beta21, beta22, beta23, beta31, beta32, beta33, gamma1, gamma2, gamma3, delta1, delta2, delta3, xs_loc, ys_loc, zs_loc, mpi_id, el_new, nelem_dg_global, nsd_jump, node_sd_jump, nrv_jump, node_rv_jump, send_length_jump, recv_length_jump, mpi_np, mpi_comm, cs_nnz_dg, cs_dg)
Computes norms of error in test mode case.
Definition COMPUTE_ENERGY_ERROR.f90:100
compute_sdof_input
subroutine compute_sdof_input(sdof_num, mpi_id, elem_mlst, local_el_num, ne_loc, cs_loc, cs_nnz_loc, sdeg_mat, nmat, u2, nnod_loc, xr_mlst, yr_mlst, zr_mlst, dt2, sdofinputab, sdofinputdispl, mpi_np, ub1, ub2, ub3)
Computes acceleration of sdof systems.
Definition COMPUTE_SDOF_INPUT.f90:48
exchange_double
subroutine exchange_double(nsend, buff_send, nrecv, buff_recv, nproc, proc_send, proc_recv, comm, status, ierr, myid)
Exchanges double between MPI processes.
Definition EXCHANGE_DOUBLE.f90:39
exchange_integer
subroutine exchange_integer(nsend, buff_send, nrecv, buff_recv, nproc, proc_send, proc_recv, comm, status, ierr, myid)
Exchanges integers between MPI processes.
Definition EXCHANGE_INTEGER.f90:38
find_pos
subroutine find_pos(n, vec, tar, ind)
Find element position on a vector.
Definition FIND_POS.f90:29
get_func_value
real *8 function get_func_value(nb_fnc, type_fnc, ind_fnc, data_fnc, nb_data_fnc, id_fnc, time
Computes time evolution function.
Definition GET_FUNC_VALUE.f90:38
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_initial_conditions
subroutine get_initial_conditions(nnod_loc, u0, u1, v1, ne_loc, cs_loc, cs_nnz_loc, nm, prop_mat, sdeg_mat, loc_n_num, xs_loc, ys_loc, zs_loc, mpi_id, dt)
Set initial conditions.
Definition GET_INITIAL_CONDITIONS.f90:48
get_mech_prop_nh_enhanced
subroutine get_mech_prop_nh_enhanced(ie, nn, nn_loc, cs_nnz_loc, cs_loc, rho_nhe, lambda_nhe, mu_nhe, qs_nh, fmax, rho_el, lambda_el, mu_el, gamma_el)
...Not-Honoring Enhanced (NHE) Implementation
Definition GET_MECH_PROP_NH_ENHANCED.f90:28
make_abc_force
subroutine make_abc_force(nn, xq, wq, dd, rho, lambda, mu, aa11, aa12, aa13, aa21, aa22, aa23, aa31, aa32, aa33, bb11, bb12, bb13, bb21, bb22, bb23, bb31, bb32, bb33, gg1, gg2, gg3, dd1, dd2, dd3, ia, ib, ja, jb, ka, kb, ux, uy, uz, vx, vy, vz, fx, f
Computes ABC's.
Definition MAKE_ABC_FORCE.f90:77
make_anelastic_coefficients_nh
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 app...
Definition MAKE_ANELASTIC_COEFFICIENTS_NH.f90:38
make_butcherarray
subroutine make_butcherarray(order, stages, a_low, b_low, c)
Makes Butcher array for Runge-Kutta scheme.
Definition MAKE_BUTCHERARRAY.f90:36
make_damping_matrix
subroutine make_damping_matrix(nn, ct, ww, dd, rho, gamma, a11, a12, a13, a21, a22, a23, a31, a32, a33, b11, b12, b13, b21, b22, b23, b31, b32, b33, gg1, gg2, gg3, dd1, dd2, dd3, mc_el, mck_el)
Make damping matrices.
Definition MAKE_DAMPING_MATRIX.f90:63
make_damping_matrix_nle
subroutine make_damping_matrix_nle(nn, ct, ww, dd, rho, gamma0, alfa11, alfa12, alfa13, alfa21, alfa22, alfa23, alfa31, alfa32, alfa33, beta11, beta12, beta13, beta21, beta22, beta23, beta31, beta32, beta33, gamma1, gamma2, gamma3, delta1, delta2, delta3, mc_el, mck_el, cs_nnz, cs, ielem, vcase, r_el, fpeak, func_type, func_indx, func_data, nfdata, nf, t_str
Make damping matrices for non-linear elastic materials.
Definition MAKE_DAMPING_MATRIX_NLE.f90:82
make_eltensor_for_cases
subroutine make_eltensor_for_cases(tcase, vcase, nn, rho_el, lambda_el, mu_el, gamma_el, nn_loc, zs_elev, zs_all, vs_nodes, thic
Assignes material properties node by node.
Definition MAKE_ELTENSOR_FOR_CASES.f90:53
make_eltensor_for_cases_nle
subroutine make_eltensor_for_cases_nle(vcase, r_el, nn, rho_el, lambda_el, mu_el, gamma_el, cs_nnz, cs, ielem, func_type, func_indx, func_data, nfdata, nf, t_stress, tag_func, yon, tcase, nnod_loc, vs_tria)
Assignes material properties node by node for non linear elasticity.
Definition MAKE_ELTENSOR_FOR_CASES_NLE.f90:50
make_expl_source
subroutine make_expl_source(nn, ct, ww, dd, sxx, syy, szz, syz, szx, sxy, nl_expl, check_ne, check_dist_ne, length_cne, ielem, facsexpl, func_type, func_indx, func_data, nfdata, nf, t
Computes the explosive moment tensor.
Definition MAKE_EXPL_SOURCE.f90:59
make_internal_force
subroutine make_internal_force(nn, xq, wq, dd, aa11, aa12, aa13, aa21, aa22, aa23, aa31, aa32, aa33, bb11, bb12, bb13, bb21, bb22, bb23, bb31, bb32, bb33, gg1, gg2, gg3, dd1, dd2, dd3, sxx, syy, szz, syz, szx, sxy, fx, fy, fz)
Makes internal forces.
Definition MAKE_INTERNAL_FORCE.f90:67
make_invariants_and_main_strain
subroutine make_invariants_and_main_strain(nn, duxdx, duydx, duzdx, duxdy, duydy, duzdy, duxdz, duydz, duzdz, r, yon)
Makes invariants for the strain tensor and computes the main strain.
Definition MAKE_INVARIANTS_AND_MAIN_STRAIN.f90:43
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
make_mass_matrix
subroutine make_mass_matrix(nn, xq, wq, dd, rho, aa11, aa12, aa13, aa21, aa22, aa23, aa31, aa32, aa33, bb11, bb12, bb13, bb21, bb22, bb23, bb31, bb32, bb33, gg1, gg2, gg3, dd1, dd2, dd3, mass)
Makes mass matrix.
Definition MAKE_MASS_MATRIX.f90:60
make_rayleigh_damping_matrix
subroutine make_rayleigh_damping_matrix(nn, ct, ww, dd, rho, a11, a12, a13, a21, a22, a23, a31, a32, a33, b11, b12, b13, b21, b22, b23, b31, b32, b33, gg1, gg2, gg3, dd1, dd2, dd3, mc_el)
Compute the Rayleigh damping matrix C = A0 M + A1 K, M = mass matrix, K= stiffness matrix.
Definition MAKE_RAYLEIGH_DAMPING_MATRIX.f90:38
make_sdof_output_files
subroutine make_sdof_output_files
Creates files to store oscillator motion.
Definition MAKE_SDOF_OUTPUT_FILES.f90:25
make_seismic_moment
subroutine make_seismic_moment(nn, ct, ww, dd, sxx, syy, szz, syz, szx, sxy, nl_sism, check_ns, check_dist_ns, length_cns, ielem, facsmom, tausmom, func_type, func_indx, func_data, nfdata, nf, t
Creates the seismic moment tensor.
Definition MAKE_SEISMIC_MOMENT.f90:60
make_strain_tensor
subroutine make_strain_tensor(nn, ct, ww, dd, alfa11, alfa12, alfa13, alfa21, alfa22, alfa23
Computes spatial derivatives of displacement.
Definition MAKE_STRAIN_TENSOR.f90:67
make_stress_tensor
subroutine make_stress_tensor(nn, lambda, mu, duxdx, duydx, duzdx, duxdy, duydy, duzdy, duxdz, duydz, duzdz, sxx, syy, szz, syz, szx, sxy)
Computes the stress tensor.
Definition MAKE_STRESS_TENSOR.f90:49
make_stress_tensor_damped
subroutine make_stress_tensor_damped(nn, lambda, mu, duxdx, duydx, duzdx, duxdy, duydy, duzdy, duxdz, duydz, duzdz, strain_visc_xx, strain_visc_xy, strain_vi
Computes the stress tensor.
Definition MAKE_STRESS_TENSOR_DAMPED.f90:57
matmul_sparse
subroutine matmul_sparse(as, nnz, jsp, isp, vet_out, n, vet_in,
Matrix-vector multiplication for sparse matrices (RCS format).
Definition MATMUL_SPARSE.f90:34
read_dime_filepg
subroutine read_dime_filepg(filec, num_nodes)
Reads dimension of monitor files.
Definition READ_DIME_FILEPG.f90:27
read_fileout
subroutine read_fileout(filename, counter, procs, nb_vec, vec)
Reads files *.out for the restart of the simulation.
Definition READ_FILEOUT.f90:30
read_fileout_vs
subroutine read_fileout_vs(filename, counter, procs, nb_vec, ncol, vec)
Reads files *.out for the restart of the simulation.
Definition READ_FILEOUT_VS.f90:30
read_file_trav_point
subroutine read_file_trav_point(filec, nmonitors_pgm, n_monitor_pgm, dist_monitor_pgm)
Reads files such as MLST.input, MLST.position or MPGM.input.
Definition READ_FILE_TRAV_POINT.f90:32
read_sdof_input_files
subroutine read_sdof_input_files
Reads inputs files with oscillator properties.
Definition READ_SDOF_INPUT_FILES.f90:25
sdof_sfs_model
subroutine sdof_sfs_model(sid, gr_acc, direction)
Response computation for 4DOF model.
Definition SDOF_SFS_MODEL.f90:28
sdof_shear_model
subroutine sdof_shear_model(sid, ndof, gr_acc, direction)
Computation for SDOF shear model.
Definition SDOF_SHEAR_MODEL.f90:28
update_damping_tensor
subroutine update_damping_tensor(nnod_loc, n_sls, frequency_range,
Update the damping tensor with the Heun method.
Definition UPDATE_DAMPING_TENSOR.f90:32
write_fileout_grid
subroutine write_fileout_grid(file_name, file_xyz, count, proc, nv, vec
Writes output results for Restart.
Definition WRITE_FILEOUT_GRID.f90:32
write_fileout_grid_sv
subroutine write_fileout_grid_sv(file_name, file_xyz, count, proc, n
Writes output results for Restart.
Definition WRITE_FILEOUT_GRID_SV.f90:32
write_output
subroutine write_output(nmonitlst, mpi_id, elem_mlst, local_el_num, ne_loc, cs_loc, cs_nnz_loc, sdeg_mat, nmat, u2, u1, u0, v1, nnod_loc, stress, strain, omega, xr_mlst, yr_mlst, zr_mlst, tt1, dt, dt2, option_out_var, count_mon, monitor_file, dbg, mu, gamma, b_instabilitycontrol, instability_maxval, b_instability_abort)
Writes output results.
Definition WRITE_OUTPUT.f90:67
write_sdof_output_files
subroutine write_sdof_output_files(tt1tmp)
Write system output files.
Definition WRITE_SDOF_OUTPUT_FILES.f90:26
dgjump
Contains structure for jump matrices.
Definition MODULES.f90:155
max_var
Set maximal bounds.
Definition MODULES.f90:54
speed_sci
Contains parameters for MDOF.
Definition MODULES.f90:725
speed_sci::ug1
real *8, dimension(:), allocatable ug1
Definition MODULES.f90:781
speed_sci::sdofforceinputbuffer
real *8, dimension(:), allocatable sdofforceinputbuffer
Definition MODULES.f90:784
speed_sci::sdofinputd
real *8, dimension(:), allocatable sdofinputd
Definition MODULES.f90:783
speed_sci::n_bld
integer *4 n_bld
Definition MODULES.f90:773
speed_sci::gr_acc_rot
real *8, dimension(:), allocatable gr_acc_rot
Definition MODULES.f90:779
speed_sci::gr_intf
real *8, dimension(:), allocatable gr_intf
Definition MODULES.f90:779
speed_sci::ug2
real *8, dimension(:), allocatable ug2
Definition MODULES.f90:781
speed_sci::sdofforceinput
real *8, dimension(:), allocatable sdofforceinput
Definition MODULES.f90:783
speed_sci::ug3
real *8, dimension(:), allocatable ug3
Definition MODULES.f90:781
speed_sci::rot_sin
real *8, dimension(:), allocatable rot_sin
Definition MODULES.f90:780
speed_sci::sdofgd
real *8, dimension(:,:), allocatable sdofgd
Definition MODULES.f90:782
speed_sci::sdofinput
real *8, dimension(:), allocatable sdofinput
Definition MODULES.f90:783
speed_sci::sdofag
real *8, dimension(:,:), allocatable sdofag
Definition MODULES.f90:782
speed_sci::rot_cos
real *8, dimension(:), allocatable rot_cos
Definition MODULES.f90:780
speed_sci::sdofnum
integer *4 sdofnum
Definition MODULES.f90:772
speed_sci::sys
type(system), dimension(:), allocatable sys
SDOF system.
Definition MODULES.f90:771
speed_sci::sdofinputbuffer
real *8, dimension(:), allocatable sdofinputbuffer
Definition MODULES.f90:784
speed_timeloop
Contains a subset of SPEED paramters (used in TIME_LOOP)
Definition MODULES.f90:519
speed_sci::system
Definition MODULES.f90:729