FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
12
17 use m_fstr_addbc
21
22 implicit none
23
24contains
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
746end module m_fstr_nonlinearmethod
This module provides functions of reconstructing.
subroutine fstr_save_originalmatrixstructure(hecmat)
This subroutine saves original matrix structure constructed originally by hecMW_matrix.
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: 1) obtain contact stiffness matrix of each contact pair and assemble ...
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_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_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_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.
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...
This module provides function to calcualte residual of nodal force.
real(kind=kreal) function, public fstr_get_x_norm_contact(hecmat, fstrmat, hecmesh)
real(kind=kreal) function, public fstr_get_norm_contact(flag, hecmesh, hecmat, fstrsolid, fstrmat)
Calculate square norm.
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)
subroutine, public fstr_update_ndforce(cstep, hecmesh, hecmat, fstrsolid, conmat)
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_updatenewton(hecmesh, hecmat, fstrsolid, time, tincr, iter, strainenergy)
変位/応力・ひずみ/内力のアップデート
Definition: fstr_Update.f90:26
subroutine fstr_updatestate(hecmesh, fstrsolid, tincr)
Update elastiplastic status.
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_set_current_config_to_mesh(hecmesh, fstrsolid, coord)
Definition: m_fstr.f90:1057
integer(kind=kint), parameter imsg
Definition: m_fstr.f90:94
integer(kind=kint), parameter ilog
FILE HANDLER.
Definition: m_fstr.f90:91
type(fstr_param), target fstrpr
GLOBAL VARIABLE INITIALIZED IN FSTR_SETUP.
Definition: m_fstr.f90:190
subroutine fstr_recover_initial_config_to_mesh(hecmesh, fstrsolid, coord)
Definition: m_fstr.f90:1070
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
logical function fstr_is_matrixstructure_changed(infoctchange)
subroutine fstr_update_contact_tangentforce(fstrsolid)
Update tangent force.
subroutine fstr_scan_contact_state(cstep, sub_step, cont_step, dt, ctalgo, hecmesh, fstrsolid, infoctchange, b)
Scanning contact state.
subroutine fstr_contactbc(iter, hecmesh, hecmat, fstrsolid)
Introduce contact stiff into global stiff matrix or mpc conditions into hecMESH.
logical function fstr_is_contact_conv(ctalgo, infoctchange, hecmesh)
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_multiplier(hecmesh, fstrsolid, ctchanged)
Update lagrangian multiplier.
subroutine fstr_update_contact0(hecmesh, fstrsolid, b)
Update lagrangian multiplier.
logical function fstr_is_contact_active()