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