heat_solve_TRAN.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_heat_solve_tran
contains
 
  subroutine heat_solve_tran ( hecMESH,hecMAT,fstrSOLID,fstrRESULT,fstrPARAM,fstrHEAT,ISTEP,total_step,current_time )
    use m_fstr
    use m_heat_mat_ass_conductivity
    use m_heat_mat_ass_capacity
    use m_heat_mat_ass_boundary
    use m_heat_init
    use m_heat_solve_main
    use m_solve_lineq
    use m_heat_io
    implicit none
    integer(kind=kint) :: ISTEP, iterALL, i, inod, mnod
    integer(kind=kint) :: total_step
    real(kind=kreal)   :: start_time, delta_time_base, delta_time, current_time, next_time, total_time, end_time, delmax, delmin
    real(kind=kreal)   :: tmpmax, dltmp, tmpmax_myrank, remain_time
    type(hecmwst_local_mesh)  :: hecmesh
    type(hecmwst_matrix)      :: hecMAT
    type(fstr_solid)          :: fstrSOLID
    type(hecmwst_result_data) :: fstrRESULT
    type(fstr_param)          :: fstrPARAM
    type(fstr_heat)           :: fstrHEAT
    type(hecmwst_local_mesh), pointer :: hecMESHmpc
    type(hecmwst_matrix), pointer :: hecMATmpc
    integer(kind=kint), parameter :: miniter = 4
    logical :: is_end, outflag
 
    call hecmw_mpc_mat_init(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
 
    start_time = 0.0d0
    do i = 1, istep - 1
      start_time = start_time + fstrheat%STEP_EETIME(i)
    enddo
 
    total_time = start_time + current_time
 
    delta_time_base = fstrheat%STEP_DLTIME(istep)
    end_time = fstrheat%STEP_EETIME(istep)
    delmin = fstrheat%STEP_DELMIN(istep)
    delmax = fstrheat%STEP_DELMAX(istep)
 
    is_end = .false.
    hecmat%NDOF = 1
    hecmat%Iarray(98) = 1 !Assembly complete
    hecmat%X = 0.0d0
 
    if(fstrheat%beta == -1.0d0)then
      if(fstrheat%is_steady == 1)then
        fstrheat%beta = 1.0d0
      else
        fstrheat%beta = 0.5d0
      endif
    endif
 
    if(fstrheat%is_steady /= 1 .and. total_step == 1) then
      call heat_output_result(hecmesh, fstrheat, fstrsolid, 0, total_time, .true.)
      call heat_output_visual(hecmesh, fstrresult, fstrheat, fstrsolid, 0, total_time, .true.)
    endif
 
    !C--------------------   START TRANSIENT LOOP   ------------------------
    tr_loop: do
 
      if(end_time <= current_time + delta_time_base + delta_time_base*1.0d-6) then
        delta_time_base = end_time - current_time
      endif
      if( 0.0d0 < delmin .and. fstrheat%timepoint_id > 0 ) then
        remain_time = get_remain_to_next_timepoints(total_time, 0.0d0, fstrparam%timepoints(fstrheat%timepoint_id))
        delta_time = dmin1(delta_time_base, remain_time)
      else
        delta_time = delta_time_base
      endif
      next_time = current_time + delta_time
      total_time = start_time + next_time
 
      if( fstrheat%is_steady == 1 ) then
        is_end = .true.
      else
        if( (end_time - next_time) / end_time < 1.d-12 ) is_end = .true.
      endif
 
      if( 0.0d0 < delmin .and. fstrheat%timepoint_id > 0 ) then
        outflag = is_end .or. is_at_timepoints(total_time, 0.0d0, fstrparam%timepoints(fstrheat%timepoint_id))
      else
        outflag = is_end
      endif
 
      if( hecmesh%my_rank.eq.0 ) then
        write(imsg,"(a,i8,a,1pe12.5,a,1pe12.5)") " ** Increment No. :", total_step, ", total time: ", &
        & total_time, ", delta t: ", delta_time
        write(*,   "(a,i8,a,1pe12.5,a,1pe12.5)") " ** Increment No. :", total_step, ", total time: ", &
        & total_time, ", delta t: ", delta_time
      endif
 
      if(delta_time_base < delmin .and. (.not. is_end))then
        if(hecmesh%my_rank == 0) write(imsg,*) ' !!! DELTA TIME EXCEEDED TOLERANCE OF TIME INCREMENT'
        call hecmw_abort(hecmw_comm_get_comm())
      endif
 
      call heat_solve_main(hecmesh, hecmat, hecmeshmpc, hecmatmpc, &
       & fstrsolid, fstrparam, fstrheat, istep, iterall, total_time, delta_time)
 
      if(0.0d0 < delmin)then
        tmpmax = 0.0d0
        do i = 1, hecmesh%nn_internal
          inod = fstrparam%global_local_id(1,i)
          dltmp = fstrheat%TEMP0(i) - fstrheat%TEMP(i)
          if(tmpmax < dabs(dltmp)) then
            mnod = inod
            tmpmax = dabs(dltmp)
          endif
        enddo
        tmpmax_myrank = tmpmax
        call hecmw_allreduce_r1(hecmesh, tmpmax, hecmw_max)
 
        if(tmpmax_myrank < tmpmax) mnod = -1
        call hecmw_allreduce_i1(hecmesh, mnod, hecmw_max)
 
        if(delmax < tmpmax .or. fstrheat%is_iter_max_limit)then
          if(hecmesh%my_rank == 0)then
            write(*,*) ' *** EXCEEDED TOLERANCE OF VARIATION IN TEMPERATUTE.'
            write(*,*) ' : NODE NUMBER  = ', mnod, ' : DELTA TEMP   = ', tmpmax
          endif
          delta_time_base = 0.5d0*delta_time_base
          cycle tr_loop
        endif
 
        if(iterall <= miniter) delta_time_base = delta_time_base*1.5d0
      else
        if(fstrheat%is_iter_max_limit) call hecmw_abort( hecmw_comm_get_comm() )
      endif
 
      do i = 1, hecmesh%n_node
        fstrheat%TEMP0(i) = fstrheat%TEMP(i)
      enddo
 
      call heat_output_log(hecmesh, fstrparam, fstrheat, total_step, total_time)
      call heat_output_result(hecmesh, fstrheat, fstrsolid, total_step, total_time, outflag)
      call heat_output_visual(hecmesh, fstrresult, fstrheat, fstrsolid, total_step, total_time, outflag)
      call heat_output_restart(hecmesh, fstrheat, istep, total_step, next_time, outflag)
 
      total_step = total_step + 1
      current_time = next_time
      if( is_end ) exit
    enddo tr_loop
    !C--------------------   END TRANSIENT LOOP   ------------------------
 
    call hecmw_mpc_mat_finalize(hecmesh, hecmat, hecmeshmpc, hecmatmpc)
 
  end subroutine heat_solve_tran
end module m_heat_solve_tran