1!------------------------------------------------------------------------------- 2! Copyright (c) 2019 FrontISTR Commons 3! This software is released under the MIT License, see LICENSE.txt 4!------------------------------------------------------------------------------- 5!> \brief This module contains subroutines controlling dynamic calculation 6 7module fstr_solver_dynamic 8 9 use m_fstr 10 use fstr_dynamic_nlexplicit 11 use fstr_dynamic_nlimplicit 12 use fstr_matrix_con_contact 13 use fstr_frequency_analysis !Frequency analysis module 14 15contains 16 17 !C================================================================C 18 !> Master subroutine for dynamic analysis 19 !C================================================================C 20 subroutine fstr_solve_DYNAMIC(hecMESH,hecMAT,fstrSOLID,fstrEIG & 21 ,fstrDYNAMIC,fstrRESULT,fstrPARAM & 22 ,fstrCPL,fstrFREQ, fstrMAT & 23 ,conMAT ) 24 use m_fstr_setup 25 implicit none 26 type(hecmwST_local_mesh) :: hecMESH 27 type(hecmwST_matrix) :: hecMAT 28 type(fstr_eigen) :: fstrEIG 29 type(fstr_solid) :: fstrSOLID 30 type(hecmwST_result_data) :: fstrRESULT 31 type(fstr_param) :: fstrPARAM 32 type(fstr_dynamic) :: fstrDYNAMIC 33 type(fstr_couple) :: fstrCPL !for COUPLE 34 type(fstr_freqanalysis) :: fstrFREQ 35 type(fstrST_matrix_contact_lagrange) :: fstrMAT !< type fstrST_matrix_contact_lagrange 36 type(fstr_info_contactChange) :: infoCTChange !< type fstr_info_contactChange 37 type(hecmwST_matrix), optional :: conMAT 38 integer(kind=kint) :: i, j, num_monit, ig, is, iE, ik, in, ing, iunitS, iunit, ierror, flag, limit 39 character(len=HECMW_FILENAME_LEN) :: fname, header 40 integer(kind=kint) :: restrt_step_num, ndof 41 integer(kind=kint) :: restrt_step(1) 42 43 num_monit = 0 44 45 if(dabs(fstrDYNAMIC%t_delta) < 1.0e-20) then 46 if( hecMESH%my_rank == 0 ) then 47 write(imsg,*) 'stop due to fstrDYNAMIC%t_delta = 0' 48 end if 49 call hecmw_abort( hecmw_comm_get_comm()) 50 end if 51 call fstr_dynamic_alloc( hecMESH, fstrDYNAMIC ) 52 53 !C-- file open for local use 54 !C 55 if(fstrDYNAMIC%idx_resp == 1) then ! time history analysis 56 call hecmw_ctrl_is_subdir( flag, limit ) 57 if( flag == 0)then 58 header = "" 59 else 60 header = adjustl("MONITOR/") 61 call hecmw_ctrl_make_subdir( adjustl("MONITOR/test.txt"), limit ) 62 endif 63 ig = fstrDYNAMIC%ngrp_monit 64 is = hecMESH%node_group%grp_index(ig-1)+1 65 iE = hecMESH%node_group%grp_index(ig) 66 do ik=is,iE 67 in = hecMESH%node_group%grp_item(ik) 68 if (in > hecMESH%nn_internal) cycle 69 num_monit = num_monit+1 70 ing = hecMESH%global_node_id(in) 71 iunitS = 10*(num_monit-1) 72 73 iunit = iunitS + fstrDYNAMIC%dynamic_IW4 74 write(fname,'(a,i0,a)') trim(header)//'dyna_disp_',ing,'.txt' 75 if(fstrDYNAMIC%restart_nout < 0 ) then 76 open(iunit,file=fname, position = 'append', iostat=ierror) 77 else 78 open(iunit,file=fname, status = 'replace', iostat=ierror) 79 endif 80 if( ierror /= 0 ) then 81 write(*,*) 'stop due to file opening error:',trim(fname) 82 call hecmw_abort( hecmw_comm_get_comm()) 83 end if 84 85 iunit = iunitS + fstrDYNAMIC%dynamic_IW5 86 write(fname,'(a,i0,a)') trim(header)//'dyna_velo_',ing,'.txt' 87 if(fstrDYNAMIC%restart_nout < 0 ) then 88 open(iunit,file=fname, position = 'append', iostat=ierror) 89 else 90 open(iunit,file=fname, status = 'replace', iostat=ierror) 91 endif 92 if( ierror /= 0 ) then 93 write(*,*) 'stop due to file opening error',trim(fname) 94 call hecmw_abort( hecmw_comm_get_comm()) 95 end if 96 97 iunit = iunitS + fstrDYNAMIC%dynamic_IW6 98 write(fname,'(a,i0,a)') trim(header)//'dyna_acce_',ing,'.txt' 99 if(fstrDYNAMIC%restart_nout < 0 ) then 100 open(iunit,file=fname, position = 'append', iostat=ierror) 101 else 102 open(iunit,file=fname, status = 'replace', iostat=ierror) 103 endif 104 if( ierror /= 0 ) then 105 write(*,*) 'stop due to file opening error',trim(fname) 106 call hecmw_abort( hecmw_comm_get_comm()) 107 end if 108 109 iunit = iunitS + fstrDYNAMIC%dynamic_IW7 110 write(fname,'(a,i0,a)') trim(header)//'dyna_force_',ing,'.txt' 111 if(fstrDYNAMIC%restart_nout < 0 ) then 112 open(iunit,file=fname, position = 'append', iostat=ierror) 113 else 114 open(iunit,file=fname, status = 'replace', iostat=ierror) 115 endif 116 if( ierror /= 0 ) then 117 write(*,*) 'stop due to file opening error',trim(fname) 118 call hecmw_abort( hecmw_comm_get_comm()) 119 end if 120 iunit = iunitS + fstrDYNAMIC%dynamic_IW8 121 write(fname,'(a,i0,a)') trim(header)//'dyna_strain_',ing,'.txt' 122 if(fstrDYNAMIC%restart_nout < 0 ) then 123 open(iunit,file=fname, position = 'append', iostat=ierror) 124 else 125 open(iunit,file=fname, status = 'replace', iostat=ierror) 126 endif 127 if( ierror /= 0 ) then 128 write(*,*) 'stop due to file opening error',trim(fname) 129 call hecmw_abort( hecmw_comm_get_comm()) 130 end if 131 132 iunit = iunitS + fstrDYNAMIC%dynamic_IW9 133 write(fname,'(a,i0,a)') trim(header)//'dyna_stress_',ing,'.txt' 134 if(fstrDYNAMIC%restart_nout < 0 ) then 135 open(iunit,file=fname, position = 'append', iostat=ierror) 136 else 137 open(iunit,file=fname, status = 'replace', iostat=ierror) 138 endif 139 if( ierror /= 0 ) then 140 write(*,*) 'stop due to file opening error',trim(fname) 141 call hecmw_abort( hecmw_comm_get_comm()) 142 end if 143 enddo 144 endif 145 146 !C 147 !C-- initial condition 148 !C 149 fstrDYNAMIC%DISP = 0.d0 150 fstrDYNAMIC%VEL = 0.d0 151 fstrDYNAMIC%ACC = 0.d0 152 fstrSOLID%unode(:) =0.d0 153 fstrSOLID%QFORCE(:) =0.d0 154 fstrDYNAMIC%kineticEnergy=0.d0 155 fstrDYNAMIC%strainEnergy=0.d0 156 fstrDYNAMIC%totalEnergy=0.d0 157 call fstr_UpdateState( hecMESH, fstrSOLID, 0.d0 ) 158 159 ! ---- restart 160 161 restrt_step_num = 1 162 fstrDYNAMIC%i_step = 0 163 infoCTChange%contactNode_previous = 0 164 165 if(associated(g_InitialCnd))then 166 ndof = HECMAT%NDOF 167 do j = 1, size(g_InitialCnd) 168 if(g_InitialCnd(j)%cond_name == "velocity")then 169 do i= 1, hecMESH%n_node 170 ing = g_InitialCnd(j)%intval(i) 171 if(ing <= 0) cycle 172 fstrDYNAMIC%VEL(ndof*i-(ndof-ing),1) = g_InitialCnd(j)%realval(i) 173 end do 174 elseif(g_InitialCnd(j)%cond_name == "acceleration")then 175 do i = 1, hecMESH%n_node 176 ing = g_InitialCnd(j)%intval(i) 177 if(ing <= 0) cycle 178 fstrDYNAMIC%ACC(ndof*i-(ndof-ing),1) = g_InitialCnd(j)%realval(i) 179 enddo 180 endif 181 enddo 182 endif 183 184 if(fstrDYNAMIC%restart_nout >= 0 ) then 185 call dynamic_bc_init (hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC) 186 call dynamic_bc_init_vl(hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC) 187 call dynamic_bc_init_ac(hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC) 188 endif 189 190 !restart 191 if(fstrDYNAMIC%restart_nout < 0 ) then 192 if( .not. associated( fstrSOLID%contacts ) ) then 193 call fstr_read_restart_dyna_nl(restrt_step_num,hecMESH,fstrSOLID,fstrDYNAMIC,fstrPARAM) 194 else 195 call fstr_read_restart_dyna_nl(restrt_step_num,hecMESH,fstrSOLID,fstrDYNAMIC,fstrPARAM,& 196 infoCTChange%contactNode_previous) 197 endif 198 restrt_step_num = restrt_step_num + 1 199 fstrDYNAMIC%restart_nout = - fstrDYNAMIC%restart_nout 200 hecMAT%Iarray(98) = 1 201 end if 202 203 if(fstrDYNAMIC%idx_resp == 1) then ! time history analysis 204 205 if(fstrDYNAMIC%idx_eqa == 1) then ! implicit dynamic analysis 206 if(.not. associated( fstrSOLID%contacts ) ) then 207 call fstr_solve_dynamic_nlimplicit(1, hecMESH,hecMAT,fstrSOLID,fstrEIG & 208 ,fstrDYNAMIC,fstrRESULT,fstrPARAM & 209 ,fstrCPL, restrt_step_num ) 210 elseif( fstrPARAM%contact_algo == kcaSLagrange ) then 211 ! ---- For Parallel Contact with Multi-Partition Domains 212 if(paraContactFlag.and.present(conMAT)) then 213 call fstr_solve_dynamic_nlimplicit_contactSLag(1, hecMESH,hecMAT,fstrSOLID,fstrEIG & 214 ,fstrDYNAMIC,fstrRESULT,fstrPARAM & 215 ,fstrCPL,fstrMAT,restrt_step_num,infoCTChange & 216 ,conMAT ) 217 else 218 call fstr_solve_dynamic_nlimplicit_contactSLag(1, hecMESH,hecMAT,fstrSOLID,fstrEIG & 219 ,fstrDYNAMIC,fstrRESULT,fstrPARAM & 220 ,fstrCPL,fstrMAT,restrt_step_num,infoCTChange ) 221 endif 222 endif 223 224 else if(fstrDYNAMIC%idx_eqa == 11) then ! explicit dynamic analysis 225 call fstr_solve_dynamic_nlexplicit(hecMESH,hecMAT,fstrSOLID,fstrEIG & 226 ,fstrDYNAMIC,fstrRESULT,fstrPARAM,infoCTChange & 227 ,fstrCPL, restrt_step_num ) 228 endif 229 230 else if(fstrDYNAMIC%idx_resp == 2) then ! frequency response analysis 231 232 if( fstrPARAM%nlgeom ) then 233 if( hecMESH%my_rank .eq. 0 ) then 234 write(imsg,*) 'stop: steady-state harmonic response analysis is not yet available !' 235 end if 236 call hecmw_abort( hecmw_comm_get_comm()) 237 end if 238 239 if( hecMESH%my_rank .eq. 0 ) then 240 call fstr_solve_frequency_analysis(hecMESH, hecMAT, fstrSOLID, fstrEIG, fstrDYNAMIC, & 241 fstrRESULT, fstrPARAM, fstrCPL, fstrFREQ, fstrMAT, & 242 restrt_step_num) 243 end if 244 end if 245 246 !C-- file close for local use 247 if(fstrDYNAMIC%idx_resp == 1) then ! time history analysis 248 do i=1,num_monit 249 iunitS = 10*(i-1) 250 close(iunitS + fstrDYNAMIC%dynamic_IW4) 251 close(iunitS + fstrDYNAMIC%dynamic_IW5) 252 close(iunitS + fstrDYNAMIC%dynamic_IW6) 253 close(iunitS + fstrDYNAMIC%dynamic_IW7) 254 close(iunitS + fstrDYNAMIC%dynamic_IW8) 255 close(iunitS + fstrDYNAMIC%dynamic_IW9) 256 enddo 257 endif 258 !C-- end of finalization 259 260 call fstr_dynamic_finalize( fstrDYNAMIC ) 261 call hecMAT_finalize( hecMAT ) 262 263 deallocate( fstrEIG%mass ) 264 265 end subroutine fstr_solve_DYNAMIC 266 267end module fstr_solver_dynamic 268