FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
fstr_solve_QuasiNewton.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 
7 module m_fstr_quasinewton
8 
9  use m_fstr_nonlinearmethod
10 
11  implicit none
12  ! parameters for line search
13  real(kind=kreal), parameter :: c_line_search=2.0, psi0_line_search=1.0d-2
14  real(kind=kreal), parameter :: delta_wolfe=0.2, sigma_wolfe=0.9, eps_wolfe=1.0d-3
15  real(kind=kreal), parameter :: omega_wolfe=0.001, delta_approx_wolfe=0.7
16  real(kind=kreal) :: c_wolfe, q_wolfe
17  integer, parameter :: n_mem_max=10
18  integer(kind=kint), parameter :: maxiter_ls = 10
19  integer(kind=kint), parameter :: pot_type=1
20 
21 contains
22 
25  subroutine fstr_quasi_newton( cstep, hecMESH, hecMAT, fstrSOLID, fstrPARAM, &
26  restrt_step_num, sub_step, ctime, dtime )
27  implicit none
28 
29  integer, intent(in) :: cstep
30  type (hecmwST_local_mesh) :: hecMESH
31  type (hecmwST_matrix) :: hecMAT
32  type (fstr_solid) :: fstrSOLID
33  integer, intent(in) :: sub_step
34  real(kind=kreal), intent(in) :: ctime
35  real(kind=kreal), intent(in) :: dtime
36  type (fstr_param) :: fstrPARAM
37  type (hecmwST_matrix_lagrange) :: hecLagMAT
38 
39  type (hecmwST_local_mesh), pointer :: hecMESHmpc
40  type (hecmwST_matrix), pointer :: hecMATmpc
41  integer(kind=kint) :: ndof
42  integer(kind=kint) :: i, iter
43  integer(kind=kint) :: stepcnt
44  integer(kind=kint) :: restrt_step_num
45  real(kind=kreal) :: tt0, tt, res, qnrm, rres, tincr, xnrm, dunrm, rxnrm
46  logical :: isLinear = .false.
47  integer(kind=kint) :: iterStatus
48 
49  real(kind=kreal), allocatable :: z_k(:), s_k(:,:), y_k(:,:), g_prev(:), rho_k(:)
50  real(kind=kreal) :: sdoty
51  integer :: n_mem
52  integer :: len_vector
53  integer(kind=kint) :: k
54 
55  integer :: u_debug
56  integer :: max_iter_bak
57  logical :: flag_approx_Wolfe
58 
59  max_iter_bak = fstrsolid%step_ctrl(cstep)%max_iter
60  fstrsolid%step_ctrl(cstep)%max_iter = 100*fstrsolid%step_ctrl(cstep)%max_iter
61 
62  call hecmw_mpc_mat_init(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
63 
64  if(.not. fstrpr%nlgeom)then
65  islinear = .true.
66  endif
67 
68  hecmat%NDOF = hecmesh%n_dof
69  ndof = hecmat%NDOF
70 
71  stepcnt = 0
72 
73  tincr = dtime
74  if( fstrsolid%step_ctrl(cstep)%solution == stepstatic ) tincr = 0.d0
75  call fstr_init_newton(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, heclagmat, ndof)
76  fstrsolid%GL0(:) = fstrsolid%GL(:) !store external load at du=0
77 
78  !! initialize du for non-zero Dirichlet condition
79  call fstr_addbc(cstep, hecmesh, hecmat, fstrsolid, fstrparam, heclagmat, 1, rhsvector=fstrsolid%dunode)
80  !! update residual vector
81  call fstr_calc_residual_vector(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam)
82 
83  len_vector = hecmesh%n_node*ndof
84  allocate(z_k(len_vector))
85  allocate(s_k(len_vector, n_mem_max))
86  allocate(y_k(len_vector, n_mem_max))
87  allocate(g_prev(len_vector))
88  allocate(rho_k(n_mem_max))
89  z_k(:) = 0.0d0
90  s_k(:,:) = 0.0d0
91  y_k(:,:) = 0.0d0
92  g_prev(:) = 0.0d0
93  rho_k(:) = 0.0d0
94  do i=1,len_vector
95  y_k(i,1) = -hecmat%B(i)
96  enddo
97 
98  ! parameter to judge Wolfe/approx Wolfe selection
99  c_wolfe = 0.0d0
100  q_wolfe = 0.0d0
101  flag_approx_wolfe = .false.
102  n_mem = 1
103  ! ----- Inner Iteration, lagrange multiplier constant
104  do iter=1,fstrsolid%step_ctrl(cstep)%max_iter
105  stepcnt = stepcnt+1
106 
107  ! ----- calculate search direction by limited BFGS method
108  do i=1,hecmesh%n_node*ndof
109  g_prev(i) = -hecmat%B(i)
110  enddo
111  call fstr_calc_direction_lbfgs(hecmesh, g_prev, s_k, y_k, rho_k, z_k, n_mem)
112 
113  ! ! ----- Set Boundary condition
114  call fstr_addbc_to_direction_vector(z_k, hecmesh,fstrsolid, cstep)
115 
116  !----- line search of step length
117  call fstr_line_search_along_direction(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k)
118 
119  ! ----- update the small displacement and the displacement for 1step
120  ! \delta u^k => solver's solution
121  ! \Delta u_{n+1}^{k} = \Delta u_{n+1}^{k-1} + \delta u^k
122  do i = 1, hecmesh%n_node*ndof
123  fstrsolid%dunode(i) = fstrsolid%dunode(i) + hecmat%X(i)
124  enddo
125 
126  !! set du for non-zero Dirichlet condition
127  ! call fstr_AddBC(cstep, hecMESH, hecMAT, fstrSOLID, fstrPARAM, hecLagMAT, 1, RHSvector=fstrSOLID%dunode)
128  call fstr_calc_residual_vector(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam)
129 
130  ! ----- check convergence
131  iterstatus = fstr_check_iteration_converged(hecmesh, hecmat, fstrsolid, ndof, iter, sub_step, cstep )
132  if (iterstatus == kitrconverged) exit
133  if (iterstatus == kitrdiverged .or. iterstatus==kitrfloatingerror) return
134  ! if (iterStatus == kitrDiverged) exit
135  ! if (iterStatus==kitrFloatingError) return
136 
137  n_mem = min(n_mem+1, n_mem_max)
138  do k=n_mem, 2, -1
139  s_k(:,k) = s_k(:,k-1)
140  y_k(:,k) = y_k(:,k-1)
141  rho_k(k) = rho_k(k-1)
142  enddo
143 
144  do i = 1, hecmesh%n_node*ndof
145  s_k(i,1) = hecmat%X(i)
146  y_k(i,1) = -hecmat%B(i) - g_prev(i)
147  enddo
148  call hecmw_innerproduct_r(hecmesh,ndof,s_k(:,1), y_k(:,1), sdoty)
149  if (abs(sdoty) < 1.0d-10) then
150  rho_k(1) = 0.0d0
151  else
152  rho_k(1) = 1.0d0/sdoty
153  endif
154  enddo
155  ! ----- end of inner loop
156 
157  fstrsolid%NRstat_i(knstmaxit) = max(fstrsolid%NRstat_i(knstmaxit),iter) ! logging newton iteration(maxtier)
158  fstrsolid%NRstat_i(knstsumit) = fstrsolid%NRstat_i(knstsumit) + iter ! logging newton iteration(sum of iter)
159 
160  ! ----- update the total displacement
161  ! u_{n+1} = u_{n} + \Delta u_{n+1}
162  do i=1,hecmesh%n_node*ndof
163  fstrsolid%unode(i) = fstrsolid%unode(i) + fstrsolid%dunode(i)
164  enddo
165 
166  call fstr_updatestate( hecmesh, fstrsolid, tincr )
167 
168  fstrsolid%CutBack_stat = 0
169  call hecmw_mpc_mat_finalize(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
170 
171  fstrsolid%step_ctrl(cstep)%max_iter = max_iter_bak
172  end subroutine fstr_quasi_newton
173 
174  subroutine fstr_calc_direction_lbfgs(hecMESH, g_prev, s_k, y_k, rho_k, z_k, n_mem)
175  implicit none
176 
177  type (hecmwST_local_mesh) :: hecMESH
178  real(kind=kreal) :: z_k(:), s_k(:,:), y_k(:,:), g_prev(:), rho_k(:)
179  integer :: n_mem
180 
181  real(kind=kreal), allocatable :: q(:)
182  real(kind=kreal) :: alpha(n_mem), beta
183  real(kind=kreal) :: sdotq, ysq, gamma, ydotz, g_max
184 
185  integer :: len_vector, ndof
186  integer(kind=kint) :: k,i
187 
188  ndof = hecmesh%n_dof
189  len_vector = hecmesh%n_node*hecmesh%n_dof
190  allocate(q(len_vector))
191 
192  q(1:len_vector) = g_prev(1:len_vector)
193 
194  do k=1, n_mem
195  call hecmw_innerproduct_r(hecmesh,ndof,s_k(:,k), q, sdotq)
196  alpha(k) = rho_k(k) * sdotq
197  do i=1, len_vector
198  q(i) = q(i) - alpha(k)*y_k(i,k)
199  enddo
200  enddo
201  call hecmw_innerproduct_r(hecmesh,ndof,y_k(:,1), y_k(:,1), ysq)
202  if (n_mem==1) then
203  call hecmw_absmax_r(hecmesh, ndof, g_prev, g_max)
204  ! if (g_max==0.0d0) then
205  ! write(6,*) 'gradient of potential is zero-vector'
206  ! stop
207  ! endif
208  ! gamma = 1.0d0/g_max
209  gamma = 1.0d0
210  else if (abs(rho_k(1)) < 1.0d-10) then
211  gamma = 1.0d0
212  else
213  gamma = 1.0d0/(rho_k(1)*ysq)
214  endif
215 
216  do i=1, len_vector
217  z_k(i) = gamma*q(i)
218  enddo
219 
220  do k=n_mem, 1, -1
221  call hecmw_innerproduct_r(hecmesh,ndof,y_k(:,k), z_k, ydotz)
222  beta = rho_k(k)*ydotz
223  do i=1, len_vector
224  z_k(i)=z_k(i) + s_k(i,k)*(alpha(k)-beta)
225  enddo
226  enddo
227  deallocate(q)
228  end subroutine fstr_calc_direction_lbfgs
229 
230  subroutine fstr_addbc_to_direction_vector(z_k, hecMESH,fstrSOLID, cstep)
231  implicit none
232  type (hecmwST_local_mesh) :: hecMESH
233  type (fstr_solid) :: fstrSOLID
234  integer(kind=kint) :: cstep
235  real(kind=kreal) :: z_k(:)
236 
237  integer(kind=kint) :: ig0, grpid, ig, ityp, idofS, idofE, iS0, iE0, ik, in, idof, ndof
238 
239  ndof = hecmesh%n_dof
240  ! ----- Prescibed displacement Boundary Conditions
241  do ig0 = 1, fstrsolid%BOUNDARY_ngrp_tot
242  grpid = fstrsolid%BOUNDARY_ngrp_GRPID(ig0)
243  if( .not. fstr_isboundaryactive( fstrsolid, grpid, cstep ) ) cycle
244  ig = fstrsolid%BOUNDARY_ngrp_ID(ig0)
245  ityp = fstrsolid%BOUNDARY_ngrp_type(ig0)
246  idofs = ityp/10
247  idofe = ityp - idofs*10
248  !
249  is0 = hecmesh%node_group%grp_index(ig-1) + 1
250  ie0 = hecmesh%node_group%grp_index(ig )
251  !
252  do ik = is0, ie0
253  in = hecmesh%node_group%grp_item(ik)
254  !
255  do idof = idofs, idofe
256  z_k(ndof*(in-1)+idof) = 0.0d0
257  enddo
258  enddo
259  enddo
260  end subroutine fstr_addbc_to_direction_vector
261 
262  subroutine fstr_apply_alpha0(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k, h_prime, pot)
263  implicit none
264  type (hecmwST_local_mesh) :: hecMESH
265  type (hecmwST_matrix) :: hecMAT
266  type (fstr_solid) :: fstrSOLID
267  real(kind=kreal), intent(in) :: ctime
268  real(kind=kreal), intent(in) :: tincr
269  integer(kind=kint) :: iter
270  integer, intent(in) :: cstep
271  real(kind=kreal), intent(in) :: dtime
272  type (fstr_param) :: fstrPARAM
273  real(kind=kreal), intent(in) :: z_k(:)
274  real(kind=kreal) :: h_prime, pot
275 
276  hecmat%X(:) = 0.0d0
277  call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, hecmat%B, z_k, h_prime)
278  pot = fstr_get_potential(cstep,hecmesh,hecmat,fstrsolid,pot_type)
279  end subroutine fstr_apply_alpha0
280 
281  subroutine fstr_apply_alpha(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k, alpha, h_prime, pot)
282  implicit none
283  type (hecmwST_local_mesh) :: hecMESH
284  type (hecmwST_matrix) :: hecMAT
285  type (fstr_solid) :: fstrSOLID
286  real(kind=kreal), intent(in) :: ctime
287  real(kind=kreal), intent(in) :: tincr
288  integer(kind=kint) :: iter
289  integer, intent(in) :: cstep
290  real(kind=kreal), intent(in) :: dtime
291  type (fstr_param) :: fstrPARAM
292  real(kind=kreal), intent(in) :: z_k(:)
293  real(kind=kreal), intent(in) :: alpha
294  real(kind=kreal) :: h_prime, pot
295 
296  hecmat%X(:) = -alpha*z_k(:)
297  call fstr_calc_residual_vector_with_x(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam)
298  call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, hecmat%B, z_k, h_prime)
299  pot = fstr_get_potential_with_x(cstep,hecmesh,hecmat,fstrsolid,pot_type)
300  end subroutine fstr_apply_alpha
301 
302  subroutine fstr_init_line_search_range(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k, &
303  h_prime_0, pot_0, alpha_S, h_prime_S, pot_S, alpha_E, h_prime_E, pot_E)
304  implicit none
305  type (hecmwST_local_mesh) :: hecMESH
306  type (hecmwST_matrix) :: hecMAT
307  type (fstr_solid) :: fstrSOLID
308  real(kind=kreal), intent(in) :: ctime
309  real(kind=kreal), intent(in) :: tincr
310  integer(kind=kint) :: iter
311  integer, intent(in) :: cstep
312  real(kind=kreal), intent(in) :: dtime
313  type (fstr_param) :: fstrPARAM
314  real(kind=kreal), intent(in) :: z_k(:)
315  real(kind=kreal), intent(in) :: h_prime_0, pot_0
316  real(kind=kreal) :: alpha_s, h_prime_s, pot_s
317  real(kind=kreal) :: alpha_e, h_prime_e, pot_e
318 
319  real(kind=kreal) :: alpha_s_new, h_prime_s_new, pot_s_new
320  real(kind=kreal) :: alpha_e_new, h_prime_e_new, pot_e_new
321  real(kind=kreal) :: alpha_tmp, h_prime_tmp, pot_tmp
322  real(kind=kreal) :: z_max
323  real(kind=kreal) :: pot_eps
324  pot_eps = eps_wolfe*c_wolfe
325 
326  alpha_s = 0.0d0
327  h_prime_s = h_prime_0
328  pot_s = pot_0
329 
330  if (iter==1) then
331  call hecmw_absmax_r(hecmesh, hecmat%NDOF, z_k, z_max)
332  alpha_e = psi0_line_search/z_max
333  else
334  alpha_e = 1.0d0
335  end if
336  call fstr_apply_alpha(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, alpha_e, h_prime_e, pot_e)
337 
338  do while (h_prime_e < 0.0d0)
339  ! if (pot_E <= pot_0 + pot_eps) then
340  ! alpha_S = alpha_E
341  ! h_prime_S = h_prime_E ! so h_prime_S < 0
342  ! pot_S = h_prime_S
343 
344  alpha_e = alpha_e * c_line_search
345  call fstr_apply_alpha(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, &
346  & fstrparam, z_k, alpha_e, h_prime_e, pot_e)
347  ! else
348  ! alpha_tmp = 2.0d0*alpha_E
349  ! h_prime_tmp = 0.0d0 ! dummy value
350  ! pot_tmp = 0.0d0 ! dummy value
351  ! call fstr_get_new_range_with_potential(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k, pot_0, &
352  ! alpha_S, h_prime_S, pot_S, alpha_tmp, h_prime_tmp, pot_tmp, alpha_E, h_prime_E, pot_E, &
353  ! alpha_S_new, h_prime_S_new, pot_S_new, alpha_E_new, h_prime_E_new, pot_E_new)
354 
355  ! alpha_S = alpha_S_new
356  ! h_prime_S = h_prime_S_new
357  ! pot_S = pot_S_new
358 
359  ! alpha_E = alpha_E_new
360  ! h_prime_E = h_prime_E_new
361  ! pot_E = pot_E_new
362  ! return
363  ! end if
364  enddo
365  end subroutine fstr_init_line_search_range
366 
367  subroutine fstr_update_range(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k, &
368  alpha_S, h_prime_S, pot_S, alpha_E, h_prime_E, pot_E, alpha_c, h_prime_c, pot_c, h_prime_0, pot_0)
369  implicit none
370  type (hecmwST_local_mesh) :: hecMESH
371  type (hecmwST_matrix) :: hecMAT
372  type (fstr_solid) :: fstrSOLID
373  real(kind=kreal), intent(in) :: ctime
374  real(kind=kreal), intent(in) :: tincr
375  integer(kind=kint) :: iter
376  integer, intent(in) :: cstep
377  real(kind=kreal), intent(in) :: dtime
378  type (fstr_param) :: fstrPARAM
379  real(kind=kreal) :: z_k(:)
380  real(kind=kreal) :: alpha_s, h_prime_s, pot_s
381  real(kind=kreal) :: alpha_e, h_prime_e, pot_e
382  real(kind=kreal) :: alpha_c, h_prime_c, pot_c
383  real(kind=kreal) :: h_prime_0, pot_0
384 
385  real(kind=kreal) :: alpha_a, h_prime_a, pot_a
386  real(kind=kreal) :: alpha_b, h_prime_b, pot_b
387  real(kind=kreal) :: alpha_a_bar, h_prime_a_bar, pot_a_bar
388  real(kind=kreal) :: alpha_b_bar, h_prime_b_bar, pot_b_bar
389  real(kind=kreal) :: alpha_c_bar, h_prime_c_bar, pot_c_bar
390 
391  call fstr_get_new_range_with_potential(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, pot_0, &
392  alpha_s, h_prime_s, pot_s, alpha_e, h_prime_e, pot_e, alpha_c, h_prime_c, pot_c, &
393  alpha_a, h_prime_a, pot_a, alpha_b, h_prime_b, pot_b)
394 
395  if (alpha_c == alpha_a) then
396  call fstr_get_secant(alpha_s, h_prime_s, alpha_a, h_prime_a, alpha_c_bar)
397  call fstr_apply_alpha(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, &
398  alpha_c_bar, h_prime_c_bar, pot_c_bar)
399  end if
400  if (alpha_c == alpha_b) then
401  call fstr_get_secant(alpha_b, h_prime_b, alpha_e, h_prime_e, alpha_c_bar)
402  call fstr_apply_alpha(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, &
403  alpha_c_bar, h_prime_c_bar, pot_c_bar)
404  end if
405 
406  if ((alpha_c == alpha_a) .or. (alpha_c == alpha_b)) then
407  call fstr_get_new_range_with_potential(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, pot_0, &
408  alpha_a, h_prime_a, pot_a, alpha_b, h_prime_b, pot_b, alpha_c_bar, h_prime_c_bar, pot_c_bar, &
409  alpha_a_bar, h_prime_a_bar, pot_a_bar, alpha_b_bar, h_prime_b_bar, pot_b_bar)
410  else
411  alpha_a_bar = alpha_a
412  h_prime_a_bar = h_prime_a
413  pot_a_bar = pot_a
414 
415  alpha_b_bar = alpha_b
416  h_prime_b_bar = h_prime_b
417  pot_b_bar = pot_b
418  end if
419 
420 
421  if ((alpha_b_bar - alpha_a_bar) > 0.66*(alpha_e - alpha_s)) then
422  alpha_c_bar = 0.5*(alpha_a_bar+alpha_b_bar)
423  call fstr_apply_alpha(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, &
424  alpha_c_bar, h_prime_c_bar, pot_c_bar)
425  call fstr_get_new_range_with_potential(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, pot_0, &
426  alpha_a, h_prime_a, pot_a, alpha_b, h_prime_b, pot_b, alpha_c_bar, h_prime_c_bar, pot_c_bar, &
427  alpha_s, h_prime_s, pot_s, alpha_e, h_prime_e, pot_e)
428  else
429  alpha_s = alpha_a_bar
430  h_prime_s = h_prime_a_bar
431  pot_s = pot_a_bar
432 
433  alpha_e = alpha_b_bar
434  h_prime_e = h_prime_b_bar
435  pot_e = pot_b_bar
436  end if
437  end subroutine fstr_update_range
438 
439  subroutine fstr_get_new_range_with_potential(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, &
440  fstrPARAM, z_k, pot_0, alpha_a, h_prime_a, pot_a, alpha_b, h_prime_b, pot_b, alpha_c, h_prime_c, pot_c, &
441  alpha_a_bar, h_prime_a_bar, pot_a_bar, alpha_b_bar, h_prime_b_bar, pot_b_bar)
442  implicit none
443  type (hecmwST_local_mesh) :: hecMESH
444  type (hecmwST_matrix) :: hecMAT
445  type (fstr_solid) :: fstrSOLID
446  real(kind=kreal), intent(in) :: ctime
447  real(kind=kreal), intent(in) :: tincr
448  integer(kind=kint) :: iter
449  integer, intent(in) :: cstep
450  real(kind=kreal), intent(in) :: dtime
451  type (fstr_param) :: fstrPARAM
452  real(kind=kreal), intent(in) :: z_k(:)
453  real(kind=kreal), intent(in) :: pot_0
454  real(kind=kreal), intent(in) :: alpha_a, h_prime_a, pot_a
455  real(kind=kreal), intent(in) :: alpha_b, h_prime_b, pot_b
456  real(kind=kreal), intent(in) :: alpha_c, h_prime_c, pot_c
457  real(kind=kreal) :: alpha_a_bar, h_prime_a_bar, pot_a_bar
458  real(kind=kreal) :: alpha_b_bar, h_prime_b_bar, pot_b_bar
459 
460  integer, parameter :: count_max=100
461  integer :: count_while
462  real(kind=kreal) :: alpha_d, h_prime_d, pot_d
463  real(kind=kreal) :: theta_ls = 0.5d0
464  real(kind=kreal) :: pot_eps
465  pot_eps = eps_wolfe*c_wolfe
466 
467  ! case of NOT (a < c < b)
468  if (alpha_c <= alpha_a .or. alpha_b <= alpha_c) then
469  alpha_a_bar = alpha_a
470  h_prime_a_bar = h_prime_a
471  pot_a_bar = pot_a
472 
473  alpha_b_bar = alpha_b
474  h_prime_b_bar = h_prime_b
475  pot_b_bar = pot_b
476  return
477  end if
478 
479  if (h_prime_c >= 0.0d0) then
480  alpha_a_bar = alpha_a
481  h_prime_a_bar = h_prime_a
482  pot_a_bar = pot_a
483 
484  alpha_b_bar = alpha_c
485  h_prime_b_bar = h_prime_c
486  pot_b_bar = pot_c
487  return
488  end if
489  ! below here, it can be assumed that h_prime_c<0
490 
491  if(pot_c <= pot_0 + pot_eps) then
492  alpha_a_bar = alpha_c
493  h_prime_a_bar = h_prime_c
494  pot_a_bar = pot_c
495 
496  alpha_b_bar = alpha_b
497  h_prime_b_bar = h_prime_b
498  pot_b_bar = pot_b
499  return
500  end if
501 
502  alpha_a_bar = alpha_a
503  h_prime_a_bar = h_prime_a
504  pot_a_bar = pot_a
505 
506  alpha_b_bar = alpha_b
507  h_prime_b_bar = h_prime_b
508  pot_b_bar = pot_b
509 
510  count_while = 0
511  do while(count_while < count_max)
512  count_while = count_while + 1
513 
514  alpha_d = (1.0d0-theta_ls)*alpha_a_bar + theta_ls*alpha_b_bar
515  call fstr_apply_alpha(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, alpha_d, h_prime_d, pot_d)
516 
517  if (h_prime_d >= 0.0d0) then
518  alpha_b_bar = alpha_d
519  h_prime_b_bar = h_prime_d
520  pot_b_bar = pot_d
521  return
522  end if
523 
524  if(pot_d <= pot_0 + pot_eps) then
525  alpha_a_bar = alpha_d
526  h_prime_a_bar = h_prime_d
527  pot_a_bar = pot_d
528  else
529  alpha_b_bar = alpha_d
530  h_prime_b_bar = h_prime_d
531  pot_b_bar = pot_d
532  end if
533  end do
534  write(6,*) 'fstr_get_new_range_with_potential reached loop count max.', hecmesh%my_rank, alpha_a_bar, alpha_b_bar
535  end subroutine fstr_get_new_range_with_potential
536 
537 
538  subroutine fstr_line_search_along_direction(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k)
539  implicit none
540  type (hecmwST_local_mesh) :: hecMESH
541  type (hecmwST_matrix) :: hecMAT
542  type (fstr_solid) :: fstrSOLID
543  real(kind=kreal), intent(in) :: ctime
544  real(kind=kreal), intent(in) :: tincr
545  integer(kind=kint) :: iter
546  integer, intent(in) :: cstep
547  real(kind=kreal), intent(in) :: dtime
548  type (fstr_param) :: fstrPARAM
549  real(kind=kreal) :: z_k(:)
550 
551  real(kind=kreal) :: alpha_s, h_prime_s, pot_s
552  real(kind=kreal) :: alpha_e, h_prime_e, pot_e
553  real(kind=kreal) :: alpha_c, h_prime_c, pot_c
554  real(kind=kreal) :: h_prime_0, pot_0
555  logical :: flag_converged
556  integer :: ndof, len_vector
557  real(kind=kreal) :: res
558 
559  integer(kind=kint) :: i, ierr, iter_ls
560  real(kind=kreal) :: z_max
561  integer :: dummy
562 
563  ndof = hecmat%NDOF
564  len_vector = hecmesh%n_node*hecmesh%n_dof
565 
566  call fstr_apply_alpha0(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, h_prime_0, pot_0)
567  if (h_prime_0 > 0.0d0) then
568  write(6,*) 'residual vector is not directed to potential decretion.', h_prime_0
569  stop
570  endif
571 
572  call fstr_init_line_search_range(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, &
573  fstrparam, z_k, h_prime_0, pot_0, alpha_s, h_prime_s, pot_s, alpha_e, h_prime_e, pot_e)
574  if( hecmesh%my_rank == 0 ) then
575  write(6,'(a, 6es27.16e3)') 'range initialized: alpha_S, alpha_E, h_prime_S, h_prime_E, pot_S, pot_E', &
576  & alpha_s, alpha_e, h_prime_s, h_prime_e, pot_s, pot_e
577  endif
578 
579  q_wolfe = 1 + delta_approx_wolfe * q_wolfe
580  c_wolfe = c_wolfe + (abs(pot_0)-c_wolfe) / q_wolfe
581 
582  flag_converged = .false.
583  iter_ls=0
584  do while (iter_ls<maxiter_ls)
585  call fstr_get_secant(alpha_s, h_prime_s, alpha_e, h_prime_e, alpha_c)
586  call fstr_apply_alpha(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, alpha_c, h_prime_c, pot_c)
587 
588  if (abs( pot_c - pot_0) <= (omega_wolfe * c_wolfe) ) then
589  flag_converged = fstr_approx_wolfe_condition(alpha_c, h_prime_c, pot_c, h_prime_0, pot_0)
590  else
591  flag_converged = fstr_wolfe_condition(alpha_c, h_prime_c, pot_c, h_prime_0, pot_0)
592  endif
593  if (flag_converged) exit
594 
595  call fstr_update_range(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam, z_k, &
596  alpha_s, h_prime_s, pot_s, alpha_e, h_prime_e, pot_e, alpha_c, h_prime_c, pot_c, h_prime_0, pot_0)
597  iter_ls = iter_ls +1
598  enddo
599  if( hecmesh%my_rank == 0 ) then
600  write(6,'(a, 1i8, 5es27.16e3)') 'converged: alpha_S, alpha_E, alpha_c, h_prime_c, pot_c', &
601  & iter_ls, &
602  & alpha_s, alpha_e, alpha_c, h_prime_c, pot_c
603  endif
604  end subroutine fstr_line_search_along_direction
605 
606  subroutine fstr_get_secant(a, Fa, b, Fb, c)
607  implicit none
608  real(kind=kreal), intent(in) :: a,b, fa,fb
609  real(kind=kreal), intent(out) :: c
610  if(fb /= fa) then
611  c = (a*fb - b*fa) / ( fb - fa )
612  else
613  c = 0.5*a + 0.5*b
614  endif
615  end subroutine fstr_get_secant
616 
617  function fstr_wolfe_condition(alpha_c, h_prime_c, pot_c, h_prime_0, pot_0) result(flag_converged)
618  implicit none
619  logical :: flag_converged
620 
621  real(kind=kreal), intent(in) :: alpha_c, h_prime_0, h_prime_c, pot_0, pot_c
622  real(kind=kreal) :: wolfe1_left, wolfe1_right, wolfe2_left, wolfe2_right
623 
624  wolfe1_left = pot_c - pot_0
625  wolfe1_right = delta_wolfe*h_prime_0*alpha_c
626 
627  wolfe2_left = h_prime_c
628  wolfe2_right = sigma_wolfe*h_prime_0
629 
630  flag_converged = (wolfe1_left<=wolfe1_right) .and. (wolfe2_left >= wolfe2_right)
631  if (flag_converged) write(6,'(a, 4es27.16e3)') 'oWolfe: ', wolfe1_left, wolfe1_right, wolfe2_left, wolfe2_right
632  end function fstr_wolfe_condition
633 
634  function fstr_approx_wolfe_condition(alpha_c, h_prime_c, pot_c, h_prime_0, pot_0) result(flag_converged)
635  implicit none
636  logical :: flag_converged
637 
638  real(kind=kreal), intent(in) :: alpha_c, h_prime_0, h_prime_c, pot_0, pot_c
639  real(kind=kreal) :: wolfe1_left, wolfe1_right, wolfe2_left, wolfe2_right
640 
641  real(kind=kreal) :: pot_eps
642  pot_eps = eps_wolfe*c_wolfe
643 
644  wolfe1_left = ( 2.0 * delta_wolfe - 1.0d0 ) * h_prime_0
645  wolfe1_right = h_prime_c
646 
647  wolfe2_left = h_prime_c
648  wolfe2_right = sigma_wolfe*h_prime_0
649 
650  flag_converged = &
651  (wolfe1_left>=wolfe1_right) &
652  .and. (wolfe2_left >= wolfe2_right) &
653  .and. (pot_c <= pot_0 + pot_eps)
654  if (flag_converged) write(6,'(a, 6es27.16e3)') 'aWolfe: ', wolfe1_left, wolfe1_right, &
655  & wolfe2_left, wolfe2_right, pot_c, pot_0 + (eps_wolfe*c_wolfe)
656  end function fstr_approx_wolfe_condition
657 end module m_fstr_quasinewton
m_fstr_quasinewton::fstr_apply_alpha
subroutine fstr_apply_alpha(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k, alpha, h_prime, pot)
Definition: fstr_solve_QuasiNewton.f90:282
m_fstr_quasinewton::fstr_quasi_newton
subroutine fstr_quasi_newton(cstep, hecMESH, hecMAT, fstrSOLID, fstrPARAM, restrt_step_num, sub_step, ctime, dtime)
This subroutine solve nonlinear solid mechanics problems by Newton-Raphson method.
Definition: fstr_solve_QuasiNewton.f90:27
m_fstr_quasinewton::fstr_calc_direction_lbfgs
subroutine fstr_calc_direction_lbfgs(hecMESH, g_prev, s_k, y_k, rho_k, z_k, n_mem)
Definition: fstr_solve_QuasiNewton.f90:175
m_fstr_quasinewton::fstr_line_search_along_direction
subroutine fstr_line_search_along_direction(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k)
Definition: fstr_solve_QuasiNewton.f90:539
m_fstr_quasinewton::psi0_line_search
real(kind=kreal), parameter psi0_line_search
Definition: fstr_solve_QuasiNewton.f90:13
m_fstr_quasinewton::eps_wolfe
real(kind=kreal), parameter eps_wolfe
Definition: fstr_solve_QuasiNewton.f90:14
m_fstr_quasinewton::omega_wolfe
real(kind=kreal), parameter omega_wolfe
Definition: fstr_solve_QuasiNewton.f90:15
m_fstr_quasinewton::delta_approx_wolfe
real(kind=kreal), parameter delta_approx_wolfe
Definition: fstr_solve_QuasiNewton.f90:15
m_fstr_quasinewton::fstr_init_line_search_range
subroutine fstr_init_line_search_range(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k, h_prime_0, pot_0, alpha_S, h_prime_S, pot_S, alpha_E, h_prime_E, pot_E)
Definition: fstr_solve_QuasiNewton.f90:304
m_fstr_quasinewton::fstr_wolfe_condition
logical function fstr_wolfe_condition(alpha_c, h_prime_c, pot_c, h_prime_0, pot_0)
Definition: fstr_solve_QuasiNewton.f90:618
m_fstr_quasinewton::fstr_addbc_to_direction_vector
subroutine fstr_addbc_to_direction_vector(z_k, hecMESH, fstrSOLID, cstep)
Definition: fstr_solve_QuasiNewton.f90:231
m_fstr_quasinewton::q_wolfe
real(kind=kreal) q_wolfe
Definition: fstr_solve_QuasiNewton.f90:16
m_fstr_quasinewton::n_mem_max
integer, parameter n_mem_max
Definition: fstr_solve_QuasiNewton.f90:17
m_fstr_quasinewton
This module provides functions on nonlinear analysis.
Definition: fstr_solve_QuasiNewton.f90:7
m_fstr_quasinewton::sigma_wolfe
real(kind=kreal), parameter sigma_wolfe
Definition: fstr_solve_QuasiNewton.f90:14
m_fstr_quasinewton::fstr_get_new_range_with_potential
subroutine fstr_get_new_range_with_potential(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k, pot_0, alpha_a, h_prime_a, pot_a, alpha_b, h_prime_b, pot_b, alpha_c, h_prime_c, pot_c, alpha_a_bar, h_prime_a_bar, pot_a_bar, alpha_b_bar, h_prime_b_bar, pot_b_bar)
Definition: fstr_solve_QuasiNewton.f90:442
m_fstr_quasinewton::pot_type
integer(kind=kint), parameter pot_type
Definition: fstr_solve_QuasiNewton.f90:19
m_fstr_quasinewton::maxiter_ls
integer(kind=kint), parameter maxiter_ls
Definition: fstr_solve_QuasiNewton.f90:18
m_fstr_quasinewton::fstr_apply_alpha0
subroutine fstr_apply_alpha0(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k, h_prime, pot)
Definition: fstr_solve_QuasiNewton.f90:263
m_fstr_quasinewton::fstr_get_secant
subroutine fstr_get_secant(a, Fa, b, Fb, c)
Definition: fstr_solve_QuasiNewton.f90:607
m_fstr_nonlinearmethod::fstr_init_newton
subroutine fstr_init_newton(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, hecLagMAT, ndof)
Definition: fstr_solve_NonLinear.f90:622
m_fstr_nonlinearmethod
This module provides functions on nonlinear analysis.
Definition: fstr_solve_NonLinear.f90:7
m_fstr_quasinewton::delta_wolfe
real(kind=kreal), parameter delta_wolfe
Definition: fstr_solve_QuasiNewton.f90:14
m_fstr_nonlinearmethod::fstr_check_iteration_converged
integer(kind=kint) function fstr_check_iteration_converged(hecMESH, hecMAT, fstrSOLID, ndof, iter, sub_step, cstep)
\breaf This function check iteration status
Definition: fstr_solve_NonLinear.f90:650
m_fstr_quasinewton::c_line_search
real(kind=kreal), parameter c_line_search
Definition: fstr_solve_QuasiNewton.f90:13
m_fstr_quasinewton::fstr_update_range
subroutine fstr_update_range(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM, z_k, alpha_S, h_prime_S, pot_S, alpha_E, h_prime_E, pot_E, alpha_c, h_prime_c, pot_c, h_prime_0, pot_0)
Definition: fstr_solve_QuasiNewton.f90:369
m_fstr_quasinewton::fstr_approx_wolfe_condition
logical function fstr_approx_wolfe_condition(alpha_c, h_prime_c, pot_c, h_prime_0, pot_0)
Definition: fstr_solve_QuasiNewton.f90:635
m_fstr_quasinewton::c_wolfe
real(kind=kreal) c_wolfe
Definition: fstr_solve_QuasiNewton.f90:16