fstr_solve_NLGEOM.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
 
module m_fstr_solve_nlgeom
  use m_fstr
  use m_static_lib
  use m_static_output
  use m_fstr_nonlinearmethod
  use m_fstr_quasinewton
  use m_fstr_nodalkinematics, only: fstr_begin_nodal_kinematics_step
  use m_fstr_restart
  use fstr_matrix_con_contact
  use m_fstr_timeinc
  use m_fstr_cutback
  use mcontact
  use m_solve_lineq_contact
 
  implicit none
 
contains
 
  !======================================================================!
  subroutine fstr_solve_nlgeom(hecMESH,hecMAT,fstrSOLID,hecLagMAT,fstrPARAM,conMAT)
    type (hecmwST_local_mesh)              :: hecMESH
    type (hecmwST_matrix    )              :: hecMAT
    type (fstr_param       )               :: fstrPARAM
    type (fstr_solid       )               :: fstrSOLID
    type (hecmwST_matrix_lagrange)         :: hecLagMAT
    type (fstr_info_contactChange)         :: infoCTChange, infoCTChange_bak
    type (hecmwST_matrix    )              :: conMAT
 
    integer(kind=kint) :: ndof, nn
    integer(kind=kint) :: j, i, tot_step, step_count, tot_step_print, CBbound
    integer(kind=kint) :: sub_step
    integer(kind=kint) :: restart_step_num, restart_substep_num
    real(kind=kreal)   :: ctime, dtime, endtime, factor
    real(kind=kreal)   :: time_1, time_2
    logical            :: ctchanged, is_OutPoint, is_interaction_active
 
    if(hecmesh%my_rank==0) call fstr_timeinc_printstatus_init
 
    hecmat%NDOF = hecmesh%n_dof
 
    ndof = hecmat%NDOF
    nn = ndof*ndof
 
    is_interaction_active = ( associated( fstrsolid%contacts ) .or. associated( fstrsolid%embeds ) )
 
    if( fstrsolid%TEMP_ngrp_tot>0 .and. hecmesh%hecmw_flag_initcon==1 ) then
      fstrsolid%last_temp = 0.0d0
      fstrsolid%temperature = 0.0d0
      do j=1, size(hecmesh%node_init_val_item)
        i = hecmesh%node_init_val_index(j)
        fstrsolid%last_temp(j) = hecmesh%node_init_val_item(i)
        fstrsolid%temperature(j) = hecmesh%node_init_val_item(i)
      end do
    endif
    if( fstrsolid%TEMP_ngrp_tot>0 .and. associated(g_initialcnd) ) then
      fstrsolid%last_temp = 0.0d0
      fstrsolid%temperature = 0.0d0
      do j=1,size(g_initialcnd)
          if( g_initialcnd(j)%cond_name=="temperature" ) then
            if( .not. associated(fstrsolid%temperature) ) then
                allocate( fstrsolid%temperature( hecmesh%n_node ) )
                allocate( fstrsolid%temp_bak( hecmesh%n_node ) )
                allocate( fstrsolid%last_temp( hecmesh%n_node ) )
            endif
            do i= 1, hecmesh%n_node
              fstrsolid%last_temp(i) = g_initialcnd(j)%realval(i)
              fstrsolid%temperature(i) = fstrsolid%last_temp(i)
            enddo
          endif
      end do
    endif
 
    if( associated( fstrsolid%contacts ) ) then
      call initialize_contact_output_vectors(fstrsolid,hecmat)
      call setup_contact_elesurf_for_area( 1, hecmesh, fstrsolid )
    endif
    if( fstrsolid%n_embeds > 0 ) call initialize_embed_vectors(fstrsolid,hecmat)
 
    restart_step_num    = 1
    restart_substep_num = 1
    fstrsolid%unode = 0.0d0
    step_count = 0 !**
    infoctchange%contactNode_previous = 0
    infoctchange%contactNode_current  = 0
    if( fstrsolid%restart_nout < 0 ) then
      call fstr_read_restart(restart_step_num,restart_substep_num,step_count,ctime,dtime,hecmesh,fstrsolid, &
        fstrparam,infoctchange%contactNode_previous)
      hecmat%Iarray(98) = 1
      call fstr_set_time( ctime )
      call fstr_set_timeinc_base( dtime )
      fstrsolid%restart_nout = - fstrsolid%restart_nout
    else
      call fstr_static_output( 1, 0, 0.d0, hecmesh, fstrsolid, fstrparam, fstrpr%solution_type, .true., 0.d0 )
    endif
 
    fstrsolid%FACTOR = 0.0d0
    call fstr_begin_nodal_kinematics_step( hecmesh, fstrsolid, hecmat%NDOF )
    call fstr_cutback_init( hecmesh, fstrsolid, fstrparam )
    call fstr_cutback_save( fstrsolid, infoctchange, infoctchange_bak )
 
    do tot_step=1, fstrsolid%nstep_tot
      tot_step_print = tot_step+restart_step_num-1
      if(hecmesh%my_rank==0) write(*,*) ''
      if(hecmesh%my_rank==0) write(*,'(a,i5)') ' loading step=',tot_step_print
 
      if( fstrsolid%TEMP_ngrp_tot>0 ) then
        do j=1, hecmesh%n_node
          fstrsolid%temp_bak(j) = fstrsolid%temperature(j)
        end do
      endif
      call fstr_updatestate( hecmesh, fstrsolid, 0.0d0 )
 
      call fstr_begin_nodal_kinematics_step( hecmesh, fstrsolid, hecmat%NDOF )
      fstrsolid%unode_bak(:) = fstrsolid%unode(:)
 
      ! -------------------------------------------------------------------------
      !      STEP LOOP
      ! -------------------------------------------------------------------------
      sub_step = restart_substep_num
      do while(.true.)
 
        ! ----- time history of factor
        call fstr_timeinc_settimeincrement( fstrsolid%step_ctrl(tot_step), fstrparam, sub_step, &
          &  fstrsolid%NRstat_i, fstrsolid%NRstat_r, fstrsolid%AutoINC_stat, fstrsolid%CutBack_stat )
        if( fstrsolid%TEMP_irres > 0 ) then
          fstrsolid%FACTOR(1) = 0.d0
          fstrsolid%FACTOR(2) = 1.d0
          call table_nlsta(hecmesh,fstrsolid,tot_step,fstr_get_time()+fstr_get_timeinc(), factor)
          fstrsolid%TEMP_FACTOR = factor
        else
          call table_nlsta(hecmesh,fstrsolid,tot_step,fstr_get_time(),factor)
          fstrsolid%FACTOR(1) = factor
          call table_nlsta(hecmesh,fstrsolid,tot_step,fstr_get_time()+fstr_get_timeinc(), factor)
          fstrsolid%FACTOR(2) = factor
        endif
 
        if(hecmesh%my_rank==0) then
          write(*,'(A,I0,2(A,E12.4))') ' sub_step= ',sub_step,', &
            &  current_time=',fstr_get_time(), ', time_inc=',fstr_get_timeinc()
          write(*,'(A,2f12.7)') ' loading_factor= ', fstrsolid%FACTOR
          if( fstrsolid%TEMP_irres > 0 ) write(*,'(A,2f12.7)') ' readtemp_factor= ', fstrsolid%TEMP_FACTOR
        endif
 
        time_1 = hecmw_wtime()
 
        ! analysis algorithm ( Newton-Rapshon Method )
        if( .not. is_interaction_active ) then
          if( fstrparam%nlsolver_method == knsmnewton ) then
            call fstr_newton( tot_step, hecmesh, hecmat, fstrsolid, fstrparam,   &
            restart_step_num, sub_step, fstr_get_time(), fstr_get_timeinc() )
          else if( fstrparam%nlsolver_method == knsmquasinewton ) then
            call fstr_quasi_newton( tot_step, hecmesh, hecmat, fstrsolid, fstrparam,   &
            restart_step_num, sub_step, fstr_get_time(), fstr_get_timeinc() )
          endif
        else
          if( fstrparam%contact_algo == kcaslagrange ) then
            call fstr_newton_contactslag( tot_step, hecmesh, hecmat, fstrsolid, fstrparam, heclagmat,  &
              restart_step_num, restart_substep_num, sub_step, fstr_get_time(), fstr_get_timeinc(), infoctchange, conmat )
          else if( fstrparam%contact_algo == kcaalagrange ) then
            call fstr_newton_contactalag( tot_step, hecmesh, hecmat, fstrsolid, fstrparam,            &
              restart_step_num, restart_substep_num, sub_step, fstr_get_time(), fstr_get_timeinc(), infoctchange, conmat )
          endif
        endif
 
        ! Time Increment
        if( hecmesh%my_rank == 0 ) call fstr_timeinc_printstatus( fstrsolid%step_ctrl(tot_step), fstrparam, &
          &  tot_step_print, sub_step, fstrsolid%NRstat_i, fstrsolid%NRstat_r,   &
          &  fstrsolid%AutoINC_stat, fstrsolid%CutBack_stat )
        if( fstr_cutback_active() ) then
 
          if( fstrsolid%CutBack_stat == 0 ) then ! converged
            call fstr_cutback_save( fstrsolid, infoctchange, infoctchange_bak )  ! save analysis state
            call fstr_proceed_time()             ! current time += time increment
 
          else                                   ! not converged
            cbbound = fstrparam%ainc(fstrsolid%step_ctrl(tot_step)%AincParam_id)%CBbound
            if( fstrsolid%CutBack_stat == cbbound ) then
              if( hecmesh%my_rank == 0 ) then
                write(*,*) 'Number of successive cutback reached max number: ',cbbound
                call fstr_timeinc_printstatus_final(.false.)
              endif
              call hecmw_abort( hecmw_comm_get_comm() )
            endif
            call fstr_cutback_load( fstrsolid, infoctchange, infoctchange_bak )  ! load analysis state
            call fstr_set_contact_active( infoctchange%contactNode_current > 0 )
 
            ! restore matrix structure for slagrange contact analysis
            if( is_interaction_active .and. fstrparam%contact_algo == kcaslagrange ) then
              call fstr_mat_con_contact( tot_step, fstrparam%contact_algo, hecmat, fstrsolid, heclagmat, &
                &  infoctchange, conmat, fstr_is_contact_active())
              conmat%B(:) = 0.0d0
              call solve_lineq_contact_init(hecmesh, hecmat, heclagmat, fstr_is_matrixstruct_symmetric(fstrsolid, hecmesh))
            endif
            if( hecmesh%my_rank == 0 ) write(*,*) '### State has been restored at time =',fstr_get_time()
 
            !stop if # of substeps reached upper bound.
            if( sub_step == fstrsolid%step_ctrl(tot_step)%num_substep ) then
              if( hecmesh%my_rank == 0 ) then
                write(*,'(a,i5,a,f6.3)') '### Number of substeps reached max number: at total_step=', &
                  & tot_step_print, '  time=', fstr_get_time()
              endif
              call hecmw_abort( hecmw_comm_get_comm())
            endif
 
            ! output time
            time_2 = hecmw_wtime()
            if( hecmesh%my_rank==0) write(imsg,'(a,",",2(I8,","),f10.2)') &
              &  'step, substep, solve (sec) :', tot_step_print, sub_step, time_2 - time_1
            cycle
          endif
        else
          if( fstrsolid%CutBack_stat > 0 ) stop
          call fstr_proceed_time() ! current time += time increment
        endif
 
        step_count = step_count + 1
 
        ! ----- Restart
        if( fstrsolid%restart_nout > 0) then
          if( mod(step_count,fstrsolid%restart_nout) == 0 ) then
            call fstr_write_restart(tot_step,tot_step_print,sub_step,step_count,fstr_get_time(),  &
              & fstr_get_timeinc_base(), hecmesh,fstrsolid,fstrparam,.false.,infoctchange%contactNode_current)
          endif
        endif
 
        ! ----- Result output (include visualize output)
        is_outpoint = fstr_timeinc_istimepoint( fstrsolid%step_ctrl(tot_step), fstrparam ) &
          & .or. fstr_timeinc_isstepfinished( fstrsolid%step_ctrl(tot_step) )
        call fstr_static_output( tot_step, step_count, fstr_get_time(), hecmesh, fstrsolid, fstrparam, &
          &                      fstrpr%solution_type, is_outpoint, fstr_get_timeinc() )
 
        time_2 = hecmw_wtime()
        if( hecmesh%my_rank==0 ) then
          write(imsg,'(A,",",2(I8,","),f10.2)') 'step, substep, solve (sec) :', tot_step_print, sub_step, time_2 - time_1
          write(imsg,'(A,I0,",",1pE15.8)') '### stepcount (for output), time :', step_count, fstr_get_time()
        endif
 
        !if time reached the end time of step, exit loop.
        if( fstr_timeinc_isstepfinished( fstrsolid%step_ctrl(tot_step) ) ) exit
 
        if( sub_step == fstrsolid%step_ctrl(tot_step)%num_substep ) then
          if( hecmesh%my_rank == 0 ) then
            write(*,'(a,i5,a,f6.3)') '### Number of substeps reached max number: at total_step=', &
              & tot_step_print, '  time=', fstr_get_time()
          endif
          if( hecmesh%my_rank == 0 ) call fstr_timeinc_printstatus_final(.false.)
          stop !stop if # of substeps reached upper bound.
        endif
 
        sub_step = sub_step + 1
      enddo    !--- end of substep  loop
 
      ! ----- Restart at the end of step
      if( fstrsolid%restart_nout > 0 ) then
        call fstr_write_restart(tot_step,tot_step_print,sub_step,step_count,fstr_get_time(),fstr_get_timeinc_base(), &
          &  hecmesh,fstrsolid,fstrparam,.true.,infoctchange%contactNode_current)
      endif
      restart_substep_num = 1
      if( fstrsolid%TEMP_irres > 0 ) exit
    enddo      !--- end of tot_step loop
 
    call fstr_cutback_finalize( fstrsolid )
 
    !  message
    if(myrank == 0)then
      call fstr_timeinc_printstatus_final(.true.)
      write(imsg,'("### FSTR_SOLVE_NLGEOM FINISHED!")')
      write(*,'("### FSTR_SOLVE_NLGEOM FINISHED!")')
    endif
 
  end subroutine fstr_solve_nlgeom
 
  !C================================================================C
  !C================================================================C
  subroutine table_nlsta(hecMESH, fstrSOLID, cstep, time, f_t)
    type ( hecmwST_local_mesh ), intent(in) :: hecMESH
    type ( fstr_solid         ), intent(in) :: fstrSOLID
    integer(kind=kint), intent(in)          :: cstep
    real(kind=kreal), intent(in)            :: time       
    real(kind=kreal), intent(out)           :: f_t        
 
    integer(kind=kint) :: jj_n_amp
 
    jj_n_amp = fstrsolid%step_ctrl( cstep )%amp_id
 
    if( jj_n_amp <= 0 ) then  ! Amplitude not defined
      f_t = (time-fstrsolid%step_ctrl(cstep)%starttime)/fstrsolid%step_ctrl(cstep)%elapsetime
      if( f_t>1.d0 ) f_t=1.d0
    else
      call table_amp(hecmesh, fstrsolid, cstep, jj_n_amp, time, f_t)
    endif
 
  end subroutine table_nlsta
 
end module m_fstr_solve_nlgeom