FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
fstr_solve_NonLinear.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 
9  use m_fstr
10  use m_static_lib
11  use m_static_output
12 
13  use m_fstr_spring
15  use m_fstr_update
16  use m_fstr_ass_load
17  use m_fstr_addbc
18  use m_fstr_residual
19  use m_fstr_restart
21 
22  implicit none
23 
24 contains
25 
26 
29  subroutine fstr_newton( cstep, hecMESH, hecMAT, fstrSOLID, fstrPARAM, &
30  restrt_step_num, sub_step, ctime, dtime )
31 
32  integer, intent(in) :: cstep
33  type (hecmwST_local_mesh) :: hecMESH
34  type (hecmwST_matrix) :: hecMAT
35  type (fstr_solid) :: fstrSOLID
36  integer, intent(in) :: sub_step
37  real(kind=kreal), intent(in) :: ctime
38  real(kind=kreal), intent(in) :: dtime
39  type (fstr_param) :: fstrPARAM
40  type (fstrST_matrix_contact_lagrange) :: fstrMAT
41 
42  type (hecmwST_local_mesh), pointer :: hecMESHmpc
43  type (hecmwST_matrix), pointer :: hecMATmpc
44  integer(kind=kint) :: ndof
45  integer(kind=kint) :: i, iter
46  integer(kind=kint) :: stepcnt
47  integer(kind=kint) :: restrt_step_num
48  real(kind=kreal) :: tt0, tt, res, qnrm, rres, tincr, xnrm, dunrm, rxnrm
49  real(kind=kreal), allocatable :: coord(:), p(:)
50  logical :: isLinear = .false.
51 
52  call hecmw_mpc_mat_init(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
53 
54  if(.not. fstrpr%nlgeom)then
55  islinear = .true.
56  endif
57 
58  hecmat%NDOF = hecmesh%n_dof
59  ndof = hecmat%NDOF
60 
61  allocate(p(hecmesh%n_node*ndof))
62  allocate(coord(hecmesh%n_node*ndof))
63 
64  tincr = dtime
65  if( fstrsolid%step_ctrl(cstep)%solution == stepstatic ) tincr = 0.d0
66 
67  p = 0.0d0
68  stepcnt = 0
69  fstrsolid%dunode(:) = 0.0d0
70  fstrsolid%NRstat_i(:) = 0 ! logging newton iteration(init)
71 
72  call fstr_ass_load(cstep, ctime+dtime, hecmesh, hecmat, fstrsolid, fstrparam)
73 
74  ! ----- Inner Iteration, lagrange multiplier constant
75  do iter=1,fstrsolid%step_ctrl(cstep)%max_iter
76  stepcnt = stepcnt+1
77 
78  call fstr_stiffmatrix( hecmesh, hecmat, fstrsolid, ctime, tincr )
79  call fstr_addspring(cstep, hecmesh, hecmat, fstrsolid, fstrparam)
80 
81  ! ----- Set Boundary condition
82  call hecmw_mpc_mat_ass(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
83  call hecmw_mpc_trans_rhs(hecmesh, hecmat, hecmatmpc)
84  call fstr_addbc(cstep, hecmesh, hecmatmpc, fstrsolid, fstrparam, fstrmat, stepcnt)
85 
86  !----- SOLVE [Kt]{du}={R}
87  if( sub_step == restrt_step_num .and. iter == 1 ) hecmatmpc%Iarray(98) = 1
88  if( iter == 1 ) then
89  hecmatmpc%Iarray(97) = 2 !Force numerical factorization
90  else
91  hecmatmpc%Iarray(97) = 1 !Need numerical factorization
92  endif
93  hecmatmpc%X = 0.0d0
94  call fstr_set_current_config_to_mesh(hecmeshmpc,fstrsolid,coord)
95  call solve_lineq(hecmeshmpc,hecmatmpc)
96  call fstr_recover_initial_config_to_mesh(hecmeshmpc,fstrsolid,coord)
97  call hecmw_mpc_tback_sol(hecmesh, hecmat, hecmatmpc)
98 
99  ! ----- update the small displacement and the displacement for 1step
100  ! \delta u^k => solver's solution
101  ! \Delta u_{n+1}^{k} = \Delta u_{n+1}^{k-1} + \delta u^k
102  do i = 1, hecmesh%n_node*ndof
103  fstrsolid%dunode(i) = fstrsolid%dunode(i) + hecmat%X(i)
104  enddo
105 
106  ! ----- update the strain, stress, and internal force
107  call fstr_updatenewton(hecmesh, hecmat, fstrsolid, ctime, tincr, iter)
108 
109  ! ----- Set residual
110  if( fstrsolid%DLOAD_follow /= 0 .or. fstrsolid%CLOAD_ngrp_rot /= 0 ) &
111  & call fstr_ass_load(cstep, ctime+dtime, hecmesh, hecmat, fstrsolid, fstrparam )
112 
113  call fstr_update_ndforce(cstep, hecmesh, hecmat, fstrsolid)
114 
115  if( islinear ) exit
116 
117  ! ----- check convergence
118  call hecmw_innerproduct_r(hecmesh, ndof, hecmat%B, hecmat%B, res)
119  res = sqrt(res)
120  call hecmw_innerproduct_r(hecmesh, ndof, hecmat%X, hecmat%X, xnrm)
121  xnrm = sqrt(xnrm)
122  call hecmw_innerproduct_r(hecmesh, ndof, fstrsolid%QFORCE, fstrsolid%QFORCE, qnrm)
123  qnrm = sqrt(qnrm)
124  if (qnrm < 1.0d-8) qnrm = 1.0d0
125  if( iter == 1 ) then
126  dunrm = xnrm
127  else
128  call hecmw_innerproduct_r(hecmesh, ndof, fstrsolid%dunode, fstrsolid%dunode, dunrm)
129  dunrm = sqrt(dunrm)
130  endif
131  rres = res/qnrm
132  rxnrm = xnrm/dunrm
133  if( hecmesh%my_rank == 0 ) then
134  if (qnrm == 1.0d0) then
135  write(*,"(a,i8,a,1pe11.4,a,1pe11.4)")" iter:", iter, ", residual(abs):", rres, ", disp.corr.:", rxnrm
136  else
137  write(*,"(a,i8,a,1pe11.4,a,1pe11.4)")" iter:", iter, ", residual:", rres, ", disp.corr.:", rxnrm
138  endif
139  endif
140  if( hecmw_mat_get_flag_diverged(hecmat) == kno ) then
141  if( rres < fstrsolid%step_ctrl(cstep)%converg ) exit
142  if( rxnrm < fstrsolid%step_ctrl(cstep)%converg ) exit
143  endif
144 
145  ! ----- check divergence and NaN
146  if( iter == fstrsolid%step_ctrl(cstep)%max_iter .or. rres > fstrsolid%step_ctrl(cstep)%maxres .or. rres /= rres ) then
147  if( hecmesh%my_rank == 0) then
148  write(ilog,'(a,i5,a,i5)') '### Fail to Converge : at total_step=', cstep, ' sub_step=', sub_step
149  write( *,'(a,i5,a,i5)') ' ### Fail to Converge : at total_step=', cstep, ' sub_step=', sub_step
150  end if
151  fstrsolid%NRstat_i(knstmaxit) = max(fstrsolid%NRstat_i(knstmaxit),iter) ! logging newton iteration(maxtier)
152  fstrsolid%NRstat_i(knstsumit) = fstrsolid%NRstat_i(knstsumit) + iter ! logging newton iteration(sumofiter)
153  fstrsolid%CutBack_stat = fstrsolid%CutBack_stat + 1
154  if( iter == fstrsolid%step_ctrl(cstep)%max_iter ) fstrsolid%NRstat_i(knstdresn) = 1
155  if( rres > fstrsolid%step_ctrl(cstep)%maxres .or. rres /= rres ) fstrsolid%NRstat_i(knstdresn) = 2
156  return
157  end if
158  enddo
159  ! ----- end of inner loop
160 
161  fstrsolid%NRstat_i(knstmaxit) = max(fstrsolid%NRstat_i(knstmaxit),iter) ! logging newton iteration(maxtier)
162  fstrsolid%NRstat_i(knstsumit) = fstrsolid%NRstat_i(knstsumit) + iter ! logging newton iteration(sum of iter)
163 
164  ! ----- update the total displacement
165  ! u_{n+1} = u_{n} + \Delta u_{n+1}
166  do i=1,hecmesh%n_node*ndof
167  fstrsolid%unode(i) = fstrsolid%unode(i) + fstrsolid%dunode(i)
168  enddo
169 
170  call fstr_updatestate( hecmesh, fstrsolid, tincr )
171 
172  fstrsolid%CutBack_stat = 0
173  deallocate(coord)
174  deallocate(p)
175  call hecmw_mpc_mat_finalize(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
176  end subroutine fstr_newton
177 
178 
182  subroutine fstr_newton_contactalag( cstep, hecMESH, hecMAT, fstrSOLID, fstrPARAM, &
183  restart_step_num, restart_substep_num, sub_step, ctime, dtime, infoCTChange )
184  use mcontact
185 
186  integer, intent(in) :: cstep
187  type (hecmwST_local_mesh) :: hecMESH
188  type (hecmwST_matrix) :: hecMAT
189  type (fstr_solid) :: fstrSOLID
190  integer, intent(in) :: sub_step
191  real(kind=kreal), intent(in) :: ctime
192  real(kind=kreal), intent(in) :: dtime
193  type (fstr_param) :: fstrPARAM
194  type (fstr_info_contactChange) :: infoCTChange
195  type (fstrST_matrix_contact_lagrange) :: fstrMAT
196 
197  type (hecmwST_local_mesh), pointer :: hecMESHmpc
198  type (hecmwST_matrix), pointer :: hecMATmpc
199  integer(kind=kint) :: ndof
200  integer(kind=kint) :: ctAlgo
201  integer(kind=kint) :: i, iter
202  integer(kind=kint) :: al_step, n_al_step, stepcnt
203  real(kind=kreal) :: tt0, tt, res, res0, res1, maxv, relres, tincr
204  integer(kind=kint) :: restart_step_num, restart_substep_num
205  logical :: convg, ctchange
206  integer(kind=kint) :: n_node_global
207  real(kind=kreal), allocatable :: coord(:)
208 
209  call hecmw_mpc_mat_init(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
210 
211  ! sum of n_node among all subdomains (to be used to calc res)
212  n_node_global = hecmesh%nn_internal
213  call hecmw_allreduce_i1(hecmesh,n_node_global,hecmw_sum)
214 
215  ctalgo = fstrparam%contact_algo
216 
217  hecmat%NDOF = hecmesh%n_dof
218  ndof = hecmat%NDOF
219 
220  fstrsolid%NRstat_i(:) = 0 ! logging newton iteration(init)
221 
222  allocate(coord(hecmesh%n_node*ndof))
223 
224  tincr = dtime
225  if( fstrsolid%step_ctrl(cstep)%solution == stepstatic ) tincr = 0.0d0
226 
227  if( cstep == 1 .and. sub_step == restart_substep_num ) then
228  if(hecmesh%my_rank==0) write(*,*) "---Scanning initial contact state---"
229  call fstr_scan_contact_state( cstep, sub_step, 0, dtime, ctalgo, hecmesh, fstrsolid, infoctchange )
230  if(hecmesh%my_rank==0) write(*,*)
231  endif
232 
233  hecmat%X = 0.0d0
234 
235  stepcnt = 0
236 
237  call fstr_ass_load(cstep, ctime+dtime, hecmesh, hecmat, fstrsolid, fstrparam)
238 
239  ! ----- Augmentation loop. In case of no contact, it is inactive
240  n_al_step = fstrsolid%step_ctrl(cstep)%max_contiter
241  if( .not. fstr_is_contact_active() ) n_al_step = 1
242 
243  do al_step = 1, n_al_step
244 
245  if( hecmesh%my_rank == 0) then
246  if( n_al_step > 1 ) then
247  print *, "Augmentation step:", al_step
248  write(imsg, *) "Augmentation step:", al_step
249  endif
250  end if
251 
252  fstrsolid%dunode(:) = 0.0d0
253 
254  ! ----- Inner Iteration, lagrange multiplier constant
255  res0 = 0.0d0
256  res1 = 0.0d0
257  relres = 1.0d0
258 
259  do iter = 1,fstrsolid%step_ctrl(cstep)%max_iter
260  stepcnt = stepcnt+1
261 
262  call fstr_stiffmatrix( hecmesh, hecmat, fstrsolid, ctime, tincr )
263  call fstr_addspring(cstep, hecmesh, hecmat, fstrsolid, fstrparam)
264 
265  ! ----- Contact
266  if( fstr_is_contact_active() .and. stepcnt==1 ) then
267  maxv = hecmw_mat_diag_max( hecmat, hecmesh )
268  call fstr_set_contact_penalty( maxv )
269  endif
270  if( fstr_is_contact_active() ) then
271  call fstr_contactbc( iter, hecmesh, hecmat, fstrsolid )
272  endif
273 
274  ! ----- Set Boundary condition
275  call hecmw_mpc_mat_ass(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
276  call hecmw_mpc_trans_rhs(hecmesh, hecmat, hecmatmpc)
277  call fstr_addbc(cstep, hecmesh,hecmatmpc,fstrsolid,fstrparam,fstrmat,stepcnt)
278 
279  !----- SOLVE [Kt]{du}={R}
280  if( sub_step == restart_step_num .and. iter == 1 ) hecmatmpc%Iarray(98) = 1
281  if( iter == 1 ) then
282  hecmatmpc%Iarray(97) = 2 !Force numerical factorization
283  else
284  hecmatmpc%Iarray(97) = 1 !Need numerical factorization
285  endif
286  hecmatmpc%X = 0.0d0
287  call fstr_set_current_config_to_mesh(hecmeshmpc,fstrsolid,coord)
288  call solve_lineq(hecmeshmpc,hecmatmpc)
289  call fstr_recover_initial_config_to_mesh(hecmeshmpc,fstrsolid,coord)
290  call hecmw_mpc_tback_sol(hecmesh, hecmat, hecmatmpc)
291 
292  if( hecmesh%n_dof == 3 ) then
293  call hecmw_update_3_r (hecmesh, hecmat%X, hecmat%NP)
294  if( hecmesh%my_rank == 0 ) then
295  write(imsg, *) 'hecmw_update_3_R: OK'
296  endif
297  else if( hecmesh%n_dof == 2 ) then
298  call hecmw_update_2_r (hecmesh, hecmat%X, hecmat%NP)
299  if( hecmesh%my_rank == 0 ) then
300  write(imsg, *) 'hecmw_update_2_R: OK'
301  endif
302  endif
303 
304  ! ----- update the small displacement and the displacement for 1step
305  ! \delta u^k => solver's solution
306  ! \Delta u_{n+1}^{k} = \Delta u_{n+1}^{k-1} + \delta u^k
307  do i = 1, hecmesh%n_node*ndof
308  fstrsolid%dunode(i) = fstrsolid%dunode(i)+hecmat%X(i)
309  enddo
310 
311  ! ----- update the strain, stress, and internal force
312  call fstr_updatenewton(hecmesh, hecmat, fstrsolid, ctime, tincr, iter)
313 
314  ! ----- Set residual
315  if( fstrsolid%DLOAD_follow /= 0 .or. fstrsolid%CLOAD_ngrp_rot /= 0 ) &
316  call fstr_ass_load(cstep, ctime+dtime, hecmesh, hecmat, fstrsolid, fstrparam)
317 
318  call fstr_update_ndforce(cstep, hecmesh, hecmat, fstrsolid)
319 
320  if( fstr_is_contact_active() ) then
321  call fstr_update_contact0(hecmesh, fstrsolid, hecmat%B)
322  endif
323 
324  res = fstr_get_residual(hecmat%B, hecmesh)
325 
326  ! ----- Gather global residual
327  res = sqrt(res)/n_node_global
328  if( iter == 1 ) res0 = res
329  if( res0 == 0.0d0 ) then
330  res0 = 1.0d0
331  else
332  relres = dabs( res1-res )/res0
333  endif
334 
335  if( hecmesh%my_rank == 0 ) then
336  write(*, '(a,i3,a,2e15.7)') ' - Residual(',iter,') =', res, relres
337  endif
338 
339  ! ----- check convergence
340  if( res < fstrsolid%step_ctrl(cstep)%converg .or. &
341  relres < fstrsolid%step_ctrl(cstep)%converg ) exit
342  res1 = res
343 
344  ! ----- check divergence and NaN
345  if( iter == fstrsolid%step_ctrl(cstep)%max_iter .or. res > fstrsolid%step_ctrl(cstep)%maxres .or. res /= res ) then
346  if( hecmesh%my_rank == 0) then
347  write( *,'(a,i5,a,i5)') ' ### Fail to Converge : at total_step=', cstep, ' sub_step=', sub_step
348  end if
349  fstrsolid%NRstat_i(knstmaxit) = max(fstrsolid%NRstat_i(knstmaxit),iter) ! logging newton iteration(maxtier)
350  fstrsolid%NRstat_i(knstsumit) = fstrsolid%NRstat_i(knstsumit) + iter ! logging newton iteration(sumofiter)
351  fstrsolid%NRstat_i(knstciter) = al_step ! logging contact iteration
352  fstrsolid%CutBack_stat = fstrsolid%CutBack_stat + 1
353  if( iter == fstrsolid%step_ctrl(cstep)%max_iter ) fstrsolid%NRstat_i(knstdresn) = 1
354  if( res > fstrsolid%step_ctrl(cstep)%maxres .or. res /= res ) fstrsolid%NRstat_i(knstdresn) = 2
355  return
356  end if
357 
358  enddo
359  ! ----- end of inner loop
360 
361  fstrsolid%NRstat_i(knstmaxit) = max(fstrsolid%NRstat_i(knstmaxit),iter) ! logging newton iteration(maxtier)
362  fstrsolid%NRstat_i(knstsumit) = fstrsolid%NRstat_i(knstsumit) + iter ! logging newton iteration(sum of iter)
363 
364  ! ----- deal with contact boundary
365  convg = .true.
366  ctchange = .false.
367  if( associated(fstrsolid%contacts) ) then
368  call fstr_update_contact_multiplier( hecmesh, fstrsolid, ctchange )
369  call fstr_scan_contact_state( cstep, sub_step, al_step, dtime, ctalgo, hecmesh, fstrsolid, infoctchange, hecmat%B )
370  if( infoctchange%contact2free+infoctchange%contact2neighbor+infoctchange%free2contact > 0 ) &
371  ctchange = .true.
372  endif
373  if( fstr_is_contact_active() ) then
374  gnt(2)=gnt(2)/iter
375  convg = fstr_is_contact_conv(ctalgo,infoctchange,hecmesh)
376  endif
377 
378  ! ----- update the total displacement
379  ! u_{n+1} = u_{n} + \Delta u_{n+1}
380  do i=1,hecmesh%n_node*ndof
381  fstrsolid%unode(i) = fstrsolid%unode(i)+fstrsolid%dunode(i)
382  enddo
383 
384  if( convg .and. (.not.ctchange) ) exit
385 
386  ! ----- check divergence
387  if( al_step >= fstrsolid%step_ctrl(cstep)%max_contiter ) then
388  if( hecmesh%my_rank == 0) then
389  write( *,'(a,i5,a,i5)') ' ### Contact failed to Converge : at total_step=', cstep, ' sub_step=', sub_step
390  end if
391  fstrsolid%NRstat_i(knstciter) = al_step ! logging contact iteration
392  fstrsolid%CutBack_stat = fstrsolid%CutBack_stat + 1
393  fstrsolid%NRstat_i(knstdresn) = 3
394  return
395  end if
396 
397  enddo
398  ! ----- end of augmentation loop
399 
400  fstrsolid%NRstat_i(knstciter) = al_step-1 ! logging contact iteration
401 
402  call fstr_updatestate( hecmesh, fstrsolid, tincr )
403 
404  deallocate(coord)
405  fstrsolid%CutBack_stat = 0
406  call hecmw_mpc_mat_finalize(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
407  end subroutine fstr_newton_contactalag
408 
409 
412  subroutine fstr_newton_contactslag( cstep, hecMESH, hecMAT, fstrSOLID, fstrPARAM, fstrMAT, &
413  restart_step_num, restart_substep_num, sub_step, ctime, dtime, infoCTChange, conMAT )
414 
415  use mcontact
418 
419  integer, intent(in) :: cstep
420  type (hecmwST_local_mesh) :: hecMESH
421  type (hecmwST_matrix) :: hecMAT
422  type (fstr_solid) :: fstrSOLID
423  integer, intent(in) :: sub_step
424  real(kind=kreal), intent(in) :: ctime
425  real(kind=kreal), intent(in) :: dtime
426  type (fstr_param) :: fstrPARAM
427  type (fstr_info_contactChange) :: infoCTChange
428  type (fstrST_matrix_contact_lagrange) :: fstrMAT
429  type (hecmwST_matrix), optional :: conMAT
430 
431  type (hecmwST_local_mesh), pointer :: hecMESHmpc
432  type (hecmwST_matrix), pointer :: hecMATmpc
433  integer(kind=kint) :: ndof
434  integer(kind=kint) :: ctAlgo
435  integer(kind=kint) :: i, iter, max_iter_contact
436  integer(kind=kint) :: stepcnt, count_step
437  real(kind=kreal) :: tt0, tt, res, res0, res1, relres, tincr, resx
438  integer(kind=kint) :: restart_step_num, restart_substep_num
439  logical :: is_mat_symmetric
440  integer(kind=kint) :: n_node_global
441  integer(kind=kint) :: contact_changed_global
442  integer(kint) :: nndof
443  real(kreal) :: q_residual,x_residual
444  real(kind=kreal), allocatable :: coord(:)
445  integer(kind=kint) :: istat
446 
447  call hecmw_mpc_mat_init(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
448 
449  ! sum of n_node among all subdomains (to be used to calc res)
450  n_node_global = hecmesh%nn_internal
451  call hecmw_allreduce_i1(hecmesh,n_node_global,hecmw_sum)
452 
453  if( hecmat%Iarray(99) == 4 .and. .not. fstr_is_matrixstruct_symmetric(fstrsolid, hecmesh) ) then
454  write(*, *) ' This type of direct solver is not yet available in such case ! '
455  write(*, *) ' Please use intel MKL direct solver !'
456  call hecmw_abort( hecmw_comm_get_comm() )
457  endif
458 
459  ctalgo = fstrparam%contact_algo
460 
461  hecmat%NDOF = hecmesh%n_dof
462  ndof = hecmat%NDOF
463 
464  fstrsolid%NRstat_i(:) = 0 ! logging newton iteration(init)
465 
466  allocate(coord(hecmesh%n_node*ndof))
467 
468  tincr = dtime
469  if( fstrsolid%step_ctrl(cstep)%solution == stepstatic ) tincr = 0.0d0
470 
471  fstrsolid%dunode(:) = 0.0d0
472 
473  if( cstep==1 .and. sub_step==restart_substep_num ) then
475  call fstr_scan_contact_state( cstep, sub_step, 0, dtime, ctalgo, hecmesh, fstrsolid, infoctchange, hecmat%B )
476  if(paracontactflag.and.present(conmat)) then
477  call hecmw_mat_copy_profile( hecmat, conmat )
478  endif
479  if ( fstr_is_contact_active() ) then
480  ! ---- For Parallel Contact with Multi-Partition Domains
481  if(paracontactflag.and.present(conmat)) then
482  call fstr_mat_con_contact(cstep, hecmat, fstrsolid, fstrmat, infoctchange, conmat)
483  else
484  call fstr_mat_con_contact(cstep, hecmat, fstrsolid, fstrmat, infoctchange)
485  endif
486  elseif( hecmat%Iarray(99)==4 ) then
487  write(*, *) ' This type of direct solver is not yet available in such case ! '
488  write(*, *) ' Please change the solver type to intel MKL direct solver !'
489  call hecmw_abort(hecmw_comm_get_comm())
490  endif
491  is_mat_symmetric = fstr_is_matrixstruct_symmetric(fstrsolid, hecmesh)
492  call solve_lineq_contact_init(hecmesh, hecmat, fstrmat, is_mat_symmetric)
493  endif
494 
495  stepcnt = 0
496 
497  call fstr_ass_load(cstep, ctime+dtime, hecmesh, hecmat, fstrsolid, fstrparam)
498 
499  if( paracontactflag.and.present(conmat) ) then
500  call hecmw_mat_clear_b(conmat)
501  endif
502 
503  if( fstr_is_contact_active() ) then
504  ! ---- For Parallel Contact with Multi-Partition Domains
505  if(paracontactflag.and.present(conmat)) then
506  call fstr_ass_load_contact(cstep, hecmesh, conmat, fstrsolid, fstrmat)
507  else
508  call fstr_ass_load_contact(cstep, hecmesh, hecmat, fstrsolid, fstrmat)
509  endif
510  endif
511 
512  fstrsolid%dunode(:) = 0.0d0
513 
514  count_step = 0
515 
516  loopforcontactanalysis: do while( .true. )
517  count_step = count_step+1
518 
519  ! ----- Inner Iteration
520  res0 = 0.d0
521  res1 = 0.d0
522  relres = 1.d0
523 
524  do iter = 1, fstrsolid%step_ctrl(cstep)%max_iter
525  call hecmw_barrier(hecmesh)
526  if( myrank == 0 ) print *,'-------------------------------------------------'
527  call hecmw_barrier(hecmesh)
528  stepcnt = stepcnt+1
529 
530  call fstr_stiffmatrix(hecmesh, hecmat, fstrsolid, ctime, tincr)
531  call fstr_addspring(cstep, hecmesh, hecmat, fstrsolid, fstrparam)
532 
533  if( paracontactflag .and. present(conmat) ) then
534  call hecmw_mat_clear( conmat )
535  conmat%X = 0.0d0
536  endif
537 
538  if( fstr_is_contact_active() ) then
539  ! ---- For Parallel Contact with Multi-Partition Domains
540  if( paracontactflag .and. present(conmat) ) then
541  call fstr_addcontactstiffness(cstep,iter,conmat,fstrmat,fstrsolid)
542  else
543  call fstr_addcontactstiffness(cstep,iter,hecmat,fstrmat,fstrsolid)
544  endif
545  endif
546 
547  ! ----- Set Boundary condition
548  call hecmw_mpc_mat_ass(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
549  call hecmw_mpc_trans_rhs(hecmesh, hecmat, hecmatmpc)
550  if(paracontactflag.and.present(conmat)) then
551  call fstr_addbc(cstep, hecmesh, hecmatmpc, fstrsolid, fstrparam, fstrmat, stepcnt, conmat)
552  else
553  call fstr_addbc(cstep, hecmesh, hecmatmpc, fstrsolid, fstrparam, fstrmat, stepcnt)
554  endif
555 
556  nndof = hecmat%N*hecmat%ndof
557 
558  !----- SOLVE [Kt]{du}={R}
559  ! ---- For Parallel Contact with Multi-Partition Domains
560  hecmatmpc%X = 0.0d0
561  call fstr_set_current_config_to_mesh(hecmeshmpc,fstrsolid,coord)
562  if(paracontactflag.and.present(conmat)) then
563  q_residual = fstr_get_norm_para_contact(hecmatmpc,fstrmat,conmat,hecmesh)
564  call solve_lineq_contact(hecmeshmpc, hecmatmpc, fstrmat, istat, 1.0d0, conmat)
565  else
566  q_residual = fstr_get_norm_contact('residualForce',hecmesh,hecmatmpc,fstrsolid,fstrmat)
567  call solve_lineq_contact(hecmeshmpc, hecmatmpc, fstrmat, istat)
568  endif
569  call fstr_recover_initial_config_to_mesh(hecmeshmpc,fstrsolid,coord)
570  call hecmw_mpc_tback_sol(hecmesh, hecmat, hecmatmpc)
571  ! ----- check matrix solver error
572  if( istat /= 0 ) then
573  if( hecmesh%my_rank == 0) then
574  write( *,'(a,i5,a,i5)') ' ### Fail to Converge : at total_step=', cstep, ' sub_step=', sub_step
575  end if
576  fstrsolid%NRstat_i(knstdresn) = 4
577  fstrsolid%CutBack_stat = fstrsolid%CutBack_stat + 1
578  return
579  end if
580 
581  x_residual = fstr_get_x_norm_contact(hecmat,fstrmat,hecmesh)
582 
583  call hecmw_innerproduct_r(hecmesh,ndof,hecmat%X,hecmat%X,resx)
584  resx = sqrt(resx)/n_node_global
585 
586  if( hecmesh%my_rank==0 ) then
587  write(*,'(a,i3,a,e15.7)') ' - ResidualX (',iter,') =',resx
588  write(*,'(a,i3,a,e15.7)') ' - ResidualX+LAG(',iter,') =',sqrt(x_residual)/n_node_global
589  write(*,'(a,i3,a,e15.7)') ' - ResidualQ (',iter,') =',sqrt(q_residual)/n_node_global
590  endif
591 
592  ! ----- update the small displacement and the displacement for 1step
593  do i = 1, hecmesh%n_node*ndof
594  fstrsolid%dunode(i) = fstrsolid%dunode(i) + hecmat%X(i)
595  enddo
596 
597  ! ----- update the Lagrange multipliers
598  if( fstr_is_contact_active() ) then
599  do i = 1, fstrmat%num_lagrange
600  fstrmat%lagrange(i) = fstrmat%lagrange(i)+hecmat%X(hecmesh%n_node*ndof+i)
601  enddo
602  endif
603 
604  ! ----- update the strain, stress, and internal force (only QFORCE)
605  call fstr_updatenewton(hecmesh, hecmat, fstrsolid, ctime, tincr, iter)
606 
607  ! ----- Set residual
608  if( fstrsolid%DLOAD_follow /= 0 .or. fstrsolid%CLOAD_ngrp_rot /= 0 ) &
609  call fstr_ass_load(cstep, ctime+dtime, hecmesh, hecmat, fstrsolid, fstrparam)
610 
611  if(paracontactflag.and.present(conmat)) then
612  call fstr_update_ndforce(cstep,hecmesh,hecmat,fstrsolid,conmat )
613  else
614  call fstr_update_ndforce(cstep,hecmesh,hecmat,fstrsolid)
615  endif
616 
617  if( fstr_is_contact_active() ) then
618  if(paracontactflag.and.present(conmat)) then
619  call hecmw_mat_clear_b( conmat )
620  call fstr_update_ndforce_contact(cstep,hecmesh,hecmat,fstrmat,fstrsolid,conmat)
621  else
622  call fstr_update_ndforce_contact(cstep,hecmesh,hecmat,fstrmat,fstrsolid)
623  endif
624  endif
625 
626  if( paracontactflag .and. present(conmat) ) then
627  res = fstr_get_norm_para_contact(hecmat,fstrmat,conmat,hecmesh)
628  else
629  res = fstr_get_norm_contact('residualForce',hecmesh,hecmat,fstrsolid,fstrmat)
630  endif
631 
632  res = sqrt(res)/n_node_global
633  if( iter == 1 ) res0 = res
634  if( res0 == 0.0d0 ) then
635  res0 =1.0d0
636  else
637  relres = dabs( res1-res )/res0
638  endif
639  if( hecmesh%my_rank == 0 ) then
640  write(*, '(a,i3,a,2e15.7)') ' - Residual(',iter,') =',res,relres
641  endif
642 
643  ! ----- check convergence
644  if( res < fstrsolid%step_ctrl(cstep)%converg .or. &
645  relres < fstrsolid%step_ctrl(cstep)%converg ) exit
646  res1 = res
647 
648  ! ----- check divergence and NaN
649  if( iter == fstrsolid%step_ctrl(cstep)%max_iter .or. res > fstrsolid%step_ctrl(cstep)%maxres .or. res /= res ) then
650  if( hecmesh%my_rank == 0) then
651  write( *,'(a,i5,a,i5)') ' ### Fail to Converge : at total_step=', cstep, ' sub_step=', sub_step
652  end if
653  fstrsolid%NRstat_i(knstmaxit) = max(fstrsolid%NRstat_i(knstmaxit),iter) ! logging newton iteration(maxtier)
654  fstrsolid%NRstat_i(knstsumit) = fstrsolid%NRstat_i(knstsumit) + iter ! logging newton iteration(sumofiter)
655  fstrsolid%NRstat_i(knstciter) = count_step ! logging contact iteration
656  fstrsolid%CutBack_stat = fstrsolid%CutBack_stat + 1
657  if( iter == fstrsolid%step_ctrl(cstep)%max_iter ) fstrsolid%NRstat_i(knstdresn) = 1
658  if( res > fstrsolid%step_ctrl(cstep)%maxres .or. res /= res ) fstrsolid%NRstat_i(knstdresn) = 2
659  return
660  end if
661 
662  enddo
663  ! ----- end of inner loop
664 
665  fstrsolid%NRstat_i(knstmaxit) = max(fstrsolid%NRstat_i(knstmaxit),iter) ! logging newton iteration(maxtier)
666  fstrsolid%NRstat_i(knstsumit) = fstrsolid%NRstat_i(knstsumit) + iter ! logging newton iteration(sum of iter)
667 
668  call fstr_scan_contact_state( cstep, sub_step, count_step, dtime, ctalgo, hecmesh, fstrsolid, infoctchange, hecmat%B )
669 
670  if( hecmat%Iarray(99) == 4 .and. .not. fstr_is_contact_active() ) then
671  write(*, *) ' This type of direct solver is not yet available in such case ! '
672  write(*, *) ' Please use intel MKL direct solver !'
673  call hecmw_abort( hecmw_comm_get_comm() )
674  endif
675 
676  is_mat_symmetric = fstr_is_matrixstruct_symmetric(fstrsolid, hecmesh)
677  contact_changed_global = 0
678  if( fstr_is_matrixstructure_changed(infoctchange) ) then
679  ! ---- For Parallel Contact with Multi-Partition Domains
680  if(paracontactflag.and.present(conmat)) then
681  call fstr_mat_con_contact( cstep, hecmat, fstrsolid, fstrmat, infoctchange, conmat)
682  else
683  call fstr_mat_con_contact( cstep, hecmat, fstrsolid, fstrmat, infoctchange)
684  endif
685  contact_changed_global = 1
686  endif
687 
688  if( fstr_is_contact_conv(ctalgo,infoctchange,hecmesh) ) exit loopforcontactanalysis
689 
690  call hecmw_allreduce_i1(hecmesh, contact_changed_global, hecmw_max)
691 
692  if (contact_changed_global > 0) then
693  call hecmw_mat_clear_b( hecmat )
694  if( paracontactflag .and. present(conmat) ) then
695  call hecmw_mat_clear_b( conmat )
696  endif
697  call solve_lineq_contact_init(hecmesh, hecmat, fstrmat, is_mat_symmetric)
698  endif
699 
700  ! ----- check divergence
701  if( count_step >= fstrsolid%step_ctrl(cstep)%max_contiter ) then
702  if( hecmesh%my_rank == 0) then
703  write( *,'(a,i5,a,i5)') ' ### Contact failed to Converge : at total_step=', cstep, ' sub_step=', sub_step
704  end if
705  fstrsolid%NRstat_i(knstciter) = count_step ! logging contact iteration
706  fstrsolid%CutBack_stat = fstrsolid%CutBack_stat + 1
707  fstrsolid%NRstat_i(knstdresn) = 3
708  return
709  end if
710 
711  ! ----- Set residual for next newton iteration
712  if(paracontactflag.and.present(conmat)) then
713  call fstr_update_ndforce(cstep,hecmesh,hecmat,fstrsolid,conmat )
714  else
715  call fstr_update_ndforce(cstep,hecmesh,hecmat,fstrsolid)
716  endif
717 
718  if( fstr_is_contact_active() ) then
719  if(paracontactflag.and.present(conmat)) then
720  call hecmw_mat_clear_b( conmat )
721  call fstr_update_ndforce_contact(cstep,hecmesh,hecmat,fstrmat,fstrsolid,conmat)
722  else
723  call fstr_update_ndforce_contact(cstep,hecmesh,hecmat,fstrmat,fstrsolid)
724  endif
725  endif
726 
727  enddo loopforcontactanalysis
728 
729  fstrsolid%NRstat_i(knstciter) = count_step ! logging contact iteration
730 
731  ! ----- update the total displacement
732  ! u_{n+1} = u_{n} + \Delta u_{n+1}
733  do i = 1, hecmesh%n_node*ndof
734  fstrsolid%unode(i) = fstrsolid%unode(i)+fstrsolid%dunode(i)
735  enddo
736 
737  call fstr_updatestate(hecmesh, fstrsolid, tincr)
738  call fstr_update_contact_tangentforce( fstrsolid )
739 
740  deallocate(coord)
741  fstrsolid%CutBack_stat = 0
742  call hecmw_mpc_mat_finalize(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
743  end subroutine fstr_newton_contactslag
744 
745 
746 end module m_fstr_nonlinearmethod
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_save_originalmatrixstructure(hecMAT)
This subroutine saves original matrix structure constructed originally by hecMW_matrix.
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: 1) obtain contact stiffness matrix of each contact pair and assemble ...
subroutine, public fstr_update_ndforce_contact(cstep, hecMESH, hecMAT, fstrMAT, fstrSOLID, conMAT)
This subroutine obtains contact nodal force vector of each contact pair and assembles it into right-h...
subroutine, public fstr_ass_load_contact(cstep, hecMESH, hecMAT, fstrSOLID, fstrMAT)
This subroutine adds initial contact force vector to the right-hand side vector \at the beginning of ...
subroutine, public fstr_addcontactstiffness(cstep, iter, hecMAT, fstrMAT, fstrSOLID)
This subroutine obtains contact stiffness matrix of each contact pair and assembles it into global st...
This module provides a function to deal with prescribed displacement.
Definition: fstr_AddBC.f90:7
subroutine fstr_addbc(cstep, hecMESH, hecMAT, fstrSOLID, fstrPARAM, fstrMAT, iter, conMAT)
Add Essential Boundary Conditions.
Definition: fstr_AddBC.f90:14
This module provides functions to take into acount external load.
subroutine fstr_ass_load(cstep, ctime, hecMESH, hecMAT, fstrSOLID, fstrPARAM)
This subroutine assmble following external force into fstrSOLIDGL and hecMATB afterwards.
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 function to calcualte residual of nodal force.
subroutine, public fstr_update_ndforce(cstep, hecMESH, hecMAT, fstrSOLID, conMAT)
real(kind=kreal) function, public fstr_get_x_norm_contact(hecMAT, fstrMAT, hecMESH)
real(kind=kreal) function, public fstr_get_residual(force, hecMESH)
Calculate magnitude of a real vector.
real(kind=kreal) function, public fstr_get_norm_para_contact(hecMAT, fstrMAT, conMAT, hecMESH)
real(kind=kreal) function, public fstr_get_norm_contact(flag, hecMESH, hecMAT, fstrSOLID, fstrMAT)
Calculate square norm.
This module provides functions to read in and write out restart fiels.
Definition: fstr_Restart.f90:8
This module provides functions to deal with spring force.
Definition: fstr_Spring.f90:7
subroutine fstr_addspring(cstep, hecMESH, hecMAT, fstrSOLID, fstrPARAM)
Definition: fstr_Spring.f90:12
This module provides function to calcualte tangent stiffness matrix.
subroutine, public fstr_stiffmatrix(hecMESH, hecMAT, fstrSOLID, time, tincr)
接線剛性マトリックスを作成するサブルーチン
This module provides function to calcualte to do updates.
Definition: fstr_Update.f90:6
subroutine fstr_updatestate(hecMESH, fstrSOLID, tincr)
Update elastiplastic status.
subroutine fstr_updatenewton(hecMESH, hecMAT, fstrSOLID, time, tincr, iter, strainEnergy)
変位/応力・ひずみ/内力のアップデート
Definition: fstr_Update.f90:26
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
subroutine fstr_recover_initial_config_to_mesh(hecMESH, fstrSOLID, coord)
Definition: m_fstr.f90:1070
integer(kind=kint), parameter imsg
Definition: m_fstr.f90:94
integer(kind=kint), parameter ilog
FILE HANDLER.
Definition: m_fstr.f90:91
subroutine fstr_set_current_config_to_mesh(hecMESH, fstrSOLID, coord)
Definition: m_fstr.f90:1057
type(fstr_param), target fstrpr
GLOBAL VARIABLE INITIALIZED IN FSTR_SETUP.
Definition: m_fstr.f90:190
integer(kind=kint), parameter kno
Definition: m_fstr.f90:31
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.
subroutine, public solve_lineq_contact(hecMESH, hecMAT, fstrMAT, istat, rf, conMAT)
This subroutine.
This modules just summarizes all modules used in static analysis.
Definition: static_LIB.f90:6
This module provides functions to output result.
This module provides functions to calcualte contact stiff matrix.
Definition: fstr_contact.f90:6
subroutine fstr_update_contact_multiplier(hecMESH, fstrSOLID, ctchanged)
Update lagrangian multiplier.
logical function fstr_is_matrixstructure_changed(infoCTChange)
logical function fstr_is_contact_conv(ctAlgo, infoCTChange, hecMESH)
subroutine fstr_contactbc(iter, hecMESH, hecMAT, fstrSOLID)
Introduce contact stiff into global stiff matrix or mpc conditions into hecMESH.
subroutine fstr_scan_contact_state(cstep, sub_step, cont_step, dt, ctAlgo, hecMESH, fstrSOLID, infoCTChange, B)
Scanning contact state.
subroutine fstr_set_contact_penalty(maxv)
real(kind=kreal), dimension(2), save gnt
1:current avarage penetration; 2:current relative tangent displacement
subroutine fstr_update_contact_tangentforce(fstrSOLID)
Update tangent force.
subroutine fstr_update_contact0(hecMESH, fstrSOLID, B)
Update lagrangian multiplier.
logical function fstr_is_contact_active()