FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
fstr_solve_NLGEOM.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
6 
8  use m_fstr
9  use m_static_lib
10  use m_static_output
12  use m_fstr_restart
14  use m_fstr_timeinc
15  use m_fstr_cutback
17 
18  implicit none
19 
20 contains
21 
22  !======================================================================!
28  subroutine fstr_solve_nlgeom(hecMESH,hecMAT,fstrSOLID,fstrMAT,fstrPARAM,conMAT)
29  type (hecmwST_local_mesh) :: hecMESH
30  type (hecmwST_matrix ) :: hecMAT
31  type (fstr_param ) :: fstrPARAM
32  type (fstr_solid ) :: fstrSOLID
33  type (fstrST_matrix_contact_lagrange) :: fstrMAT
34  type (fstr_info_contactChange) :: infoCTChange, infoCTChange_bak
35  type (hecmwST_matrix ),optional :: conMAT
36 
37  integer(kind=kint) :: ndof, nn
38  integer(kind=kint) :: j, i, tot_step, step_count, tot_step_print, CBbound
39  integer(kind=kint) :: sub_step
40  integer(kind=kint) :: restart_step_num, restart_substep_num
41  real(kind=kreal) :: ctime, dtime, endtime, factor
42  real(kind=kreal) :: time_1, time_2
43  logical :: ctchanged, is_OutPoint
44 
45  if(hecmesh%my_rank==0) call fstr_timeinc_printstatus_init
46 
47  hecmat%NDOF = hecmesh%n_dof
48 
49  ndof = hecmat%NDOF
50  nn = ndof*ndof
51 
52  if( fstrsolid%TEMP_ngrp_tot>0 .and. hecmesh%hecmw_flag_initcon==1 ) then
53  fstrsolid%last_temp = 0.0d0
54  fstrsolid%temperature = 0.0d0
55  do j=1, size(hecmesh%node_init_val_item)
56  i = hecmesh%node_init_val_index(j)
57  fstrsolid%last_temp(j) = hecmesh%node_init_val_item(i)
58  fstrsolid%temperature(j) = hecmesh%node_init_val_item(i)
59  end do
60  endif
61  if( fstrsolid%TEMP_ngrp_tot>0 .and. associated(g_initialcnd) ) then
62  fstrsolid%last_temp = 0.0d0
63  fstrsolid%temperature = 0.0d0
64  do j=1,size(g_initialcnd)
65  if( g_initialcnd(j)%cond_name=="temperature" ) then
66  if( .not. associated(fstrsolid%temperature) ) then
67  allocate( fstrsolid%temperature( hecmesh%n_node ) )
68  allocate( fstrsolid%temp_bak( hecmesh%n_node ) )
69  allocate( fstrsolid%last_temp( hecmesh%n_node ) )
70  endif
71  do i= 1, hecmesh%n_node
72  fstrsolid%last_temp(i) = g_initialcnd(j)%realval(i)
73  fstrsolid%temperature(i) = fstrsolid%last_temp(i)
74  enddo
75  endif
76  end do
77  endif
78 
79  if( associated( fstrsolid%contacts ) ) then
80  call initialize_contact_output_vectors(fstrsolid,hecmat)
81  call setup_contact_elesurf_for_area( 1, hecmesh, fstrsolid )
82  endif
83 
84  restart_step_num = 1
85  restart_substep_num = 1
86  fstrsolid%unode = 0.0d0
87  step_count = 0 !**
88  infoctchange%contactNode_previous = 0
89  infoctchange%contactNode_current = 0
90  if( fstrsolid%restart_nout < 0 ) then
91  call fstr_read_restart(restart_step_num,restart_substep_num,step_count,ctime,dtime,hecmesh,fstrsolid, &
92  fstrparam,infoctchange%contactNode_previous)
93  hecmat%Iarray(98) = 1
94  call fstr_set_time( ctime )
95  call fstr_set_timeinc_base( dtime )
96  fstrsolid%restart_nout = - fstrsolid%restart_nout
97  else
98  call fstr_static_output( 1, 0, 0.d0, hecmesh, fstrsolid, fstrparam, fstrpr%solution_type, .true. )
99  endif
100 
101  fstrsolid%FACTOR = 0.0d0
102  call fstr_cutback_init( hecmesh, fstrsolid, fstrparam )
103  call fstr_cutback_save( fstrsolid, infoctchange, infoctchange_bak )
104 
105  do tot_step=1, fstrsolid%nstep_tot
106  tot_step_print = tot_step+restart_step_num-1
107  if(hecmesh%my_rank==0) write(*,*) ''
108  if(hecmesh%my_rank==0) write(*,'(a,i5)') ' loading step=',tot_step_print
109 
110  if( fstrsolid%TEMP_ngrp_tot>0 ) then
111  do j=1, hecmesh%n_node
112  fstrsolid%temp_bak(j) = fstrsolid%temperature(j)
113  end do
114  endif
115  call fstr_updatestate( hecmesh, fstrsolid, 0.0d0 )
116 
117  ! -------------------------------------------------------------------------
118  ! STEP LOOP
119  ! -------------------------------------------------------------------------
120  sub_step = restart_substep_num
121  do while(.true.)
122 
123  ! ----- time history of factor
124  call fstr_timeinc_settimeincrement( fstrsolid%step_ctrl(tot_step), fstrparam, sub_step, &
125  & fstrsolid%NRstat_i, fstrsolid%NRstat_r, fstrsolid%AutoINC_stat, fstrsolid%CutBack_stat )
126  if( fstrsolid%TEMP_irres > 0 ) then
127  fstrsolid%FACTOR(1) = 0.d0
128  fstrsolid%FACTOR(2) = 1.d0
129  call table_nlsta(hecmesh,fstrsolid,tot_step,fstr_get_time()+fstr_get_timeinc(), factor)
130  fstrsolid%TEMP_FACTOR = factor
131  else
132  call table_nlsta(hecmesh,fstrsolid,tot_step,fstr_get_time(),factor)
133  fstrsolid%FACTOR(1) = factor
134  call table_nlsta(hecmesh,fstrsolid,tot_step,fstr_get_time()+fstr_get_timeinc(), factor)
135  fstrsolid%FACTOR(2) = factor
136  endif
137 
138  if(hecmesh%my_rank==0) then
139  write(*,'(A,I0,2(A,E12.4))') ' sub_step= ',sub_step,', &
140  & current_time=',fstr_get_time(), ', time_inc=',fstr_get_timeinc()
141  write(*,'(A,2f12.7)') ' loading_factor= ', fstrsolid%FACTOR
142  if( fstrsolid%TEMP_irres > 0 ) write(*,'(A,2f12.7)') ' readtemp_factor= ', fstrsolid%TEMP_FACTOR
143  endif
144 
145  time_1 = hecmw_wtime()
146 
147  ! analysis algorithm ( Newton-Rapshon Method )
148  if( .not. associated( fstrsolid%contacts ) ) then
149  call fstr_newton( tot_step, hecmesh, hecmat, fstrsolid, fstrparam, &
150  restart_step_num, sub_step, fstr_get_time(), fstr_get_timeinc() )
151  else
152  if( fstrparam%contact_algo == kcaslagrange ) then
153  ! For Parallel Contact with Multi-Partition Domains
154  if(paracontactflag.and.present(conmat)) then
155  call fstr_newton_contactslag( tot_step, hecmesh, hecmat, fstrsolid, fstrparam, fstrmat, &
156  restart_step_num, restart_substep_num, sub_step, fstr_get_time(), fstr_get_timeinc(), infoctchange, conmat )
157  else
158  call fstr_newton_contactslag( tot_step, hecmesh, hecmat, fstrsolid, fstrparam, fstrmat, &
159  restart_step_num, restart_substep_num, sub_step, fstr_get_time(), fstr_get_timeinc(), infoctchange )
160  endif
161  else if( fstrparam%contact_algo == kcaalagrange ) then
162  call fstr_newton_contactalag( tot_step, hecmesh, hecmat, fstrsolid, fstrparam, &
163  restart_step_num, restart_substep_num, sub_step, fstr_get_time(), fstr_get_timeinc(), infoctchange )
164  endif
165  endif
166 
167  ! Time Increment
168  if( hecmesh%my_rank == 0 ) call fstr_timeinc_printstatus( fstrsolid%step_ctrl(tot_step), fstrparam, &
169  & tot_step_print, sub_step, fstrsolid%NRstat_i, fstrsolid%NRstat_r, &
170  & fstrsolid%AutoINC_stat, fstrsolid%CutBack_stat )
171  if( fstr_cutback_active() ) then
172 
173  if( fstrsolid%CutBack_stat == 0 ) then ! converged
174  call fstr_cutback_save( fstrsolid, infoctchange, infoctchange_bak ) ! save analysis state
175  call fstr_proceed_time() ! current time += time increment
176 
177  else ! not converged
178  cbbound = fstrparam%ainc(fstrsolid%step_ctrl(tot_step)%AincParam_id)%CBbound
179  if( fstrsolid%CutBack_stat == cbbound ) then
180  if( hecmesh%my_rank == 0 ) then
181  write(*,*) 'Number of successive cutback reached max number: ',cbbound
182  call fstr_timeinc_printstatus_final(.false.)
183  endif
184  call hecmw_abort( hecmw_comm_get_comm() )
185  endif
186  call fstr_cutback_load( fstrsolid, infoctchange, infoctchange_bak ) ! load analysis state
187  call fstr_set_contact_active( infoctchange%contactNode_current > 0 )
188 
189  ! restore matrix structure for slagrange contact analysis
190  if( associated( fstrsolid%contacts ) .and. fstrparam%contact_algo == kcaslagrange ) then
191  if(paracontactflag.and.present(conmat)) then
192  call fstr_mat_con_contact( tot_step, hecmat, fstrsolid, fstrmat, infoctchange, conmat)
193  conmat%B(:) = 0.0d0
194  else
195  call fstr_mat_con_contact( tot_step, hecmat, fstrsolid, fstrmat, infoctchange)
196  endif
197  call solve_lineq_contact_init(hecmesh, hecmat, fstrmat, fstr_is_matrixstruct_symmetric(fstrsolid, hecmesh))
198  endif
199  if( hecmesh%my_rank == 0 ) write(*,*) '### State has been restored at time =',fstr_get_time()
200 
201  !stop if # of substeps reached upper bound.
202  if( sub_step == fstrsolid%step_ctrl(tot_step)%num_substep ) then
203  if( hecmesh%my_rank == 0 ) then
204  write(*,'(a,i5,a,f6.3)') '### Number of substeps reached max number: at total_step=', &
205  & tot_step_print, ' time=', fstr_get_time()
206  endif
207  call hecmw_abort( hecmw_comm_get_comm())
208  endif
209 
210  ! output time
211  time_2 = hecmw_wtime()
212  if( hecmesh%my_rank==0) write(imsg,'(a,",",2(I8,","),f10.2)') &
213  & 'step, substep, solve (sec) :', tot_step_print, sub_step, time_2 - time_1
214  cycle
215  endif
216  else
217  if( fstrsolid%CutBack_stat > 0 ) stop
218  call fstr_proceed_time() ! current time += time increment
219  endif
220 
221  step_count = step_count + 1
222 
223  ! ----- Restart
224  if( fstrsolid%restart_nout > 0 .and. mod(step_count,fstrsolid%restart_nout) == 0 ) then
225  call fstr_write_restart(tot_step,tot_step_print,sub_step,step_count,fstr_get_time(), &
226  & fstr_get_timeinc_base(), hecmesh,fstrsolid,fstrparam,.false.,infoctchange%contactNode_current)
227  endif
228 
229  ! ----- Result output (include visualize output)
230  is_outpoint = fstr_timeinc_istimepoint( fstrsolid%step_ctrl(tot_step), fstrparam ) &
231  & .or. fstr_timeinc_isstepfinished( fstrsolid%step_ctrl(tot_step) )
232  call fstr_static_output( tot_step, step_count, fstr_get_time(), hecmesh, fstrsolid, fstrparam, &
233  & fstrpr%solution_type, is_outpoint )
234 
235  time_2 = hecmw_wtime()
236  if( hecmesh%my_rank==0 ) then
237  write(imsg,'(A,",",2(I8,","),f10.2)') 'step, substep, solve (sec) :', tot_step_print, sub_step, time_2 - time_1
238  write(imsg,'(A,I0,",",1pE15.8)') '### stepcount (for output), time :', step_count, fstr_get_time()
239  endif
240 
241  !if time reached the end time of step, exit loop.
242  if( fstr_timeinc_isstepfinished( fstrsolid%step_ctrl(tot_step) ) ) exit
243 
244  if( sub_step == fstrsolid%step_ctrl(tot_step)%num_substep ) then
245  if( hecmesh%my_rank == 0 ) then
246  write(*,'(a,i5,a,f6.3)') '### Number of substeps reached max number: at total_step=', &
247  & tot_step_print, ' time=', fstr_get_time()
248  endif
249  if( hecmesh%my_rank == 0 ) call fstr_timeinc_printstatus_final(.false.)
250  stop !stop if # of substeps reached upper bound.
251  endif
252 
253  sub_step = sub_step + 1
254  enddo !--- end of substep loop
255 
256  ! ----- Restart at the end of step
257  if( fstrsolid%restart_nout > 0 ) then
258  call fstr_write_restart(tot_step,tot_step_print,sub_step,step_count,fstr_get_time(),fstr_get_timeinc_base(), &
259  & hecmesh,fstrsolid,fstrparam,.true.,infoctchange%contactNode_current)
260  endif
261  restart_substep_num = 1
262  if( fstrsolid%TEMP_irres > 0 ) exit
263  enddo !--- end of tot_step loop
264 
265  call fstr_cutback_finalize( fstrsolid )
266 
267  ! message
268  if(myrank == 0)then
269  call fstr_timeinc_printstatus_final(.true.)
270  write(imsg,'("### FSTR_SOLVE_NLGEOM FINISHED!")')
271  write(*,'("### FSTR_SOLVE_NLGEOM FINISHED!")')
272  endif
273 
274  end subroutine fstr_solve_nlgeom
275 
276  !C================================================================C
279  !C================================================================C
280  subroutine table_nlsta(hecMESH, fstrSOLID, cstep, time, f_t)
281  type ( hecmwST_local_mesh ), intent(in) :: hecMESH
282  type ( fstr_solid ), intent(in) :: fstrSOLID
283  integer(kind=kint), intent(in) :: cstep
284  real(kind=kreal), intent(in) :: time
285  real(kind=kreal), intent(out) :: f_t
286 
287  integer(kind=kint) :: i
288  integer(kind=kint) :: jj_n_amp, jj1, jj2
289  integer(kind=kint) :: s1, s2, flag
290  real(kind=kreal) :: t_1, t_2, t_t, f_1, f_2, tincre
291 
292  s1 = 0; s2 = 0
293  jj_n_amp = fstrsolid%step_ctrl( cstep )%amp_id
294 
295  if( jj_n_amp <= 0 ) then ! Amplitude not defined
296  f_t = (time-fstrsolid%step_ctrl(cstep)%starttime)/fstrsolid%step_ctrl(cstep)%elapsetime
297  if( f_t>1.d0 ) f_t=1.d0
298  else
299  tincre = fstrsolid%step_ctrl( cstep )%initdt
300  jj1 = hecmesh%amp%amp_index(jj_n_amp - 1)
301  jj2 = hecmesh%amp%amp_index(jj_n_amp)
302 
303  jj1 = jj1 + 2
304  t_t = time-fstrsolid%step_ctrl(cstep)%starttime
305 
306  ! if(jj2 .eq. 0) then
307  ! f_t = 1.0
308  if(t_t .gt. hecmesh%amp%amp_table(jj2)) then
309  f_t = hecmesh%amp%amp_val(jj2)
310  else if(t_t .le. hecmesh%amp%amp_table(jj2)) then
311  flag=0
312  do i = jj1, jj2
313  if(t_t .le. hecmesh%amp%amp_table(i)) then
314  s2 = i
315  s1 = i - 1
316  flag = 1
317  endif
318  if( flag == 1 ) exit
319  end do
320 
321  t_2 = hecmesh%amp%amp_table(s2)
322  t_1 = hecmesh%amp%amp_table(s1)
323  f_2 = hecmesh%amp%amp_val(s2)
324  f_1 = hecmesh%amp%amp_val(s1)
325  if( t_2-t_1 .lt. 1.0e-20) then
326  if(myrank == 0) then
327  write(imsg,*) 'stop due to t_2-t_1 <= 0'
328  endif
329  call hecmw_abort( hecmw_comm_get_comm())
330  endif
331  f_t = ((t_2*f_1 - t_1*f_2) + (f_2 - f_1)*t_t) / (t_2 - t_1)
332  endif
333 
334  endif
335 
336  end subroutine table_nlsta
337 
338 end module m_fstr_solve_nlgeom
This module provides functions of reconstructing.
logical function fstr_is_matrixstruct_symmetric(fstrSOLID, hecMESH)
this function judges whether sitiffness matrix is symmetric or not
subroutine fstr_mat_con_contact(cstep, hecMAT, fstrSOLID, fstrMAT, infoCTChange, conMAT)
this subroutine reconstructs node-based (stiffness) matrix structure \corresponding to contact state
This module provides functions to deal with cutback.
Definition: fstr_Cutback.f90:7
subroutine fstr_cutback_save(fstrSOLID, infoCTChange, infoCTChange_bak)
Save analysis status.
subroutine fstr_cutback_load(fstrSOLID, infoCTChange, infoCTChange_bak)
Load analysis status.
subroutine fstr_cutback_init(hecMESH, fstrSOLID, fstrPARAM)
Initializer of cutback variables.
subroutine fstr_cutback_finalize(fstrSOLID)
Finalizer of cutback variables.
logical function fstr_cutback_active()
This module provides functions on nonlinear analysis.
subroutine fstr_newton_contactalag(cstep, hecMESH, hecMAT, fstrSOLID, fstrPARAM, restart_step_num, restart_substep_num, sub_step, ctime, dtime, infoCTChange)
This subroutine solve nonlinear solid mechanics problems by Newton-Raphson method combined with Neste...
subroutine fstr_newton_contactslag(cstep, hecMESH, hecMAT, fstrSOLID, fstrPARAM, fstrMAT, restart_step_num, restart_substep_num, sub_step, ctime, dtime, infoCTChange, conMAT)
This subroutine solve nonlinear solid mechanics problems by Newton-Raphson method....
subroutine fstr_newton(cstep, hecMESH, hecMAT, fstrSOLID, fstrPARAM, restrt_step_num, sub_step, ctime, dtime)
This subroutine solve nonlinear solid mechanics problems by Newton-Raphson method.
This module provides functions to read in and write out restart fiels.
Definition: fstr_Restart.f90:8
subroutine fstr_read_restart(cstep, substep, step_count, ctime, dtime, hecMESH, fstrSOLID, fstrPARAM, contactNode)
Read in restart file.
subroutine fstr_write_restart(cstep, cstep_ext, substep, step_count, ctime, dtime, hecMESH, fstrSOLID, fstrPARAM, is_StepFinished, contactNode)
write out restart file
This module provides main suboruitne for nonliear calculation.
subroutine fstr_solve_nlgeom(hecMESH, hecMAT, fstrSOLID, fstrMAT, fstrPARAM, conMAT)
This module provides main suborutine for nonlinear calculation.
subroutine table_nlsta(hecMESH, fstrSOLID, cstep, time, f_t)
This subroutine decide the loading increment considering the amplitude definition.
This module provides functions to deal with time and increment of stress analysis.
real(kind=kreal) function fstr_get_timeinc()
logical function fstr_timeinc_istimepoint(stepinfo, fstrPARAM)
real(kind=kreal) function fstr_get_timeinc_base()
subroutine fstr_set_timeinc_base(dtime_base)
subroutine fstr_timeinc_settimeincrement(stepinfo, fstrPARAM, substep, NRstatI, NRstatR, AutoINC_stat, Cutback_stat)
real(kind=kreal) function fstr_get_time()
subroutine fstr_proceed_time()
subroutine fstr_timeinc_printstatus_final(success_flag)
subroutine fstr_timeinc_printstatus_init
subroutine fstr_set_time(time)
logical function fstr_timeinc_isstepfinished(stepinfo)
subroutine fstr_timeinc_printstatus(stepinfo, fstrPARAM, totstep, substep, NRstatI, NRstatR, AutoINC_stat, Cutback_stat)
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
integer(kind=kint) myrank
PARALLEL EXECUTION.
Definition: m_fstr.f90:80
integer(kind=kint), parameter imsg
Definition: m_fstr.f90:94
integer(kind=kint), parameter kcaslagrange
contact analysis algorithm
Definition: m_fstr.f90:53
integer(kind=kint), parameter kcaalagrange
Definition: m_fstr.f90:54
type(fstr_param), target fstrpr
GLOBAL VARIABLE INITIALIZED IN FSTR_SETUP.
Definition: m_fstr.f90:190
type(tinitialcondition), dimension(:), pointer, save g_initialcnd
Definition: m_fstr.f90:135
logical paracontactflag
PARALLEL CONTACT FLAG.
Definition: m_fstr.f90:84
This module provides functions to solve sparse system of \linear equitions in the case of contact ana...
subroutine, public solve_lineq_contact_init(hecMESH, hecMAT, fstrMAT, is_sym)
This subroutine.
This modules just summarizes all modules used in static analysis.
Definition: static_LIB.f90:6
This module provides functions to output result.
subroutine fstr_static_output(cstep, istep, time, hecMESH, fstrSOLID, fstrPARAM, flag, outflag)
Output result.