fstr_EIG_output.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_eig_output
 
contains
 
  subroutine fstr_eigen_residual(hecMESH, hecMAT, fstrEIG)
    use m_fstr
    use hecmw_solver_las
    implicit none
    type(hecmwst_local_mesh) :: hecMESH
    type(hecmwst_matrix)     :: hecMAT
    type(fstr_eigen)         :: fstrEIG
 
    integer(kind=kint) :: NNDOF, NPNDOF, i, j
    real(kind=kreal)   :: chk, shift
    real(kind=kreal), allocatable :: t(:), u(:)
 
    nndof  = hecmat%N *hecmat%NDOF
    npndof = hecmat%NP*hecmat%NDOF
 
    allocate(fstreig%resid(fstreig%nget))
    allocate(t(npndof))
    allocate(u(npndof))
 
    ! the free-free case works on [K] + sigma [M], see lanczos_set_initial_value
    shift = 0.0d0
    if(fstreig%is_free) shift = fstreig%sigma
 
    do j = 1, fstreig%nget
      u = 0.0d0
      do i = 1, nndof
        u(i) = fstreig%eigvec(i,j)
      enddo
      call hecmw_matvec(hecmesh, hecmat, u, t)
 
      chk = 0.0d0
      do i = 1, nndof
        chk = chk + (t(i) - (fstreig%eigval(j) + shift)*fstreig%mass(i)*u(i))**2
      enddo
      call hecmw_allreduce_r1(hecmesh, chk, hecmw_sum)
      fstreig%resid(j) = dsqrt(chk)
    enddo
 
    deallocate(t)
    deallocate(u)
  end subroutine fstr_eigen_residual
 
  subroutine fstr_eigen_output(hecMESH, hecMAT, fstrEIG)
    use m_fstr
    use m_fstr_eig_lanczos_util
    implicit none
    type(hecmwst_local_mesh) :: hecmesh
    type(hecmwst_matrix)     :: hecMAT
    type(fstr_eigen)         :: fstrEIG
 
    integer(kind=kint) :: N, NDOF
    integer(kind=kint) :: i, j, k, in, nget
    real(kind=kreal)   :: gm
    real(kind=kreal), allocatable :: r(:)
    real(kind=kreal), pointer     :: mass(:), eigval(:), eigvec(:,:)
 
    n      = hecmat%N
    ndof   = hecmesh%n_dof
 
    nget   =  fstreig%nget
    mass   => fstreig%mass
    eigval => fstreig%eigval
    eigvec => fstreig%eigvec
 
    allocate(fstreig%effmass(ndof*nget))
    allocate(fstreig%partfactor(ndof*nget))
    allocate(r(ndof))
    fstreig%effmass    = 0.0d0
    fstreig%partfactor = 0.0d0
 
    do i = 1, nget
      r  = 0.0d0
      gm = 0.0d0
      do j = 1, n
        do k = 1, ndof
          in = ndof*(j-1) + k
          r(k) = r(k) + mass(in)*eigvec(in,i)
          gm   = gm   + mass(in)*eigvec(in,i)*eigvec(in,i)
        enddo
      enddo
 
      call hecmw_allreduce_r(hecmesh, r, ndof, hecmw_sum)
      call hecmw_allreduce_r1(hecmesh, gm, hecmw_sum)
 
      gm = 1.0d0/gm
      do j = 1, ndof
        in = ndof*(i-1) + j
        fstreig%partfactor(in) = gm*r(j)
        fstreig%effmass(in)    = gm*r(j)*r(j)
      enddo
    enddo
 
    call eglist(hecmesh, hecmat, fstreig)
 
    if(myrank == 0)then
      write(imsg,*) ''
      write(imsg,*) '*----------------------------------------------*'
      write(imsg,*) '*--E I G E N V A L U E  C O N V E R G E N C E--*'
      write(imsg,*) 'Absolute residual =  |(||Kx - lambda*Mx||)|'
      write(imsg,*) '   Iter.#   Eigenvalue    Abs. Residual  '
      write(ilog,*) '   Iter.#   Eigenvalue    Abs. Residual  '
      write(imsg,*) '   *-----*  *---------*  *--------------*'
    endif
 
    do j = 1, fstreig%nget
      if(myrank == 0)then
        write(imsg,'(2x,i5,2x,1p5e15.6)') j, eigval(j), fstreig%resid(j)
        write(ilog,'(i5,1p5e12.4)') j, eigval(j), fstreig%resid(j)
      endif
    enddo
 
    if(myrank == 0)then
      write(imsg,*)'*    ---END Eigenvalue listing---     *'
    endif
 
    deallocate(r)
  end subroutine fstr_eigen_output
 
  subroutine fstr_eigen_make_result_step(hecMESH, fstrEIG, fstrRESULT, istep)
    use m_fstr
    use m_hecmw2fstr_mesh_conv
    use hecmw_util
    implicit none
    type(hecmwst_local_mesh)  :: hecmesh
    type(fstr_eigen)          :: fstrEIG
    type(hecmwst_result_data) :: fstrRESULT
    integer(kind=kint) :: istep, i
 
    integer(kind=kint) :: NP, NDOF, NPNDOF
    real(kind=kreal), pointer :: eigvec(:,:)
    real(kind=kreal), allocatable :: x(:)
 
    np     = hecmesh%n_node
    ndof   = hecmesh%n_dof
    npndof = np * ndof
 
    eigvec => fstreig%eigvec
    allocate(x(npndof))
    x = 0.0d0
 
    do i=1,npndof
      x(i) = eigvec(i,istep)
    enddo
    call hecmw_update_r(hecmesh, x, np, ndof)
 
    fstrresult%ng_component = 1
    fstrresult%ne_component = 0
    allocate(fstrresult%ng_dof(1))
    allocate(fstrresult%global_label(1))
    allocate(fstrresult%global_val_item(1))
    fstrresult%ng_dof(1) = 1
    fstrresult%global_label(1) = 'EIGENVALUE'
    fstrresult%global_val_item(1) = fstreig%eigval(istep)
 
    if( ndof == 6 ) then
      fstrresult%nn_component = 2
      allocate(fstrresult%nn_dof(2))
      allocate(fstrresult%node_label(2))
      allocate(fstrresult%node_val_item(6*np))
      fstrresult%nn_dof(1) = 3
      fstrresult%nn_dof(2) = 3
      fstrresult%node_label(1) = 'DISPLACEMENT'
      fstrresult%node_label(2) = 'ROTATION'
      do i = 1, np
        fstrresult%node_val_item((i-1)*6+1:(i-1)*6+3) = x((i-1)*ndof+1:(i-1)*ndof+3)
        fstrresult%node_val_item((i-1)*6+4:(i-1)*6+6) = x((i-1)*ndof+4:(i-1)*ndof+6)
      enddo
    else
      fstrresult%nn_component = 1
      allocate(fstrresult%nn_dof(1))
      allocate(fstrresult%node_label(1))
      allocate(fstrresult%node_val_item(ndof*np))
      fstrresult%nn_dof(1) = ndof
      fstrresult%node_label(1) = 'DISPLACEMENT'
      fstrresult%node_val_item = x
    endif
 
    deallocate(x)
  end subroutine
 
  subroutine fstr_eigen_output_visual(hecMESH, fstrEIG, fstrRESULT)
    use m_fstr
    use m_hecmw2fstr_mesh_conv
    use hecmw_util
    implicit none
    type(hecmwst_local_mesh)  :: hecmesh
    type(fstr_eigen)          :: fstrEIG
    type(hecmwst_result_data) :: fstrRESULT
    integer(kind=kint) :: istep, nget
    nget = fstreig%nget
    do istep = 1, nget
      call hecmw_nullify_result_data(fstrresult)
      call fstr_eigen_make_result_step(hecmesh, fstreig, fstrresult, istep)
      call fstr2hecmw_mesh_conv(hecmesh)
      call hecmw_visualize_init
      call hecmw_visualize( hecmesh, fstrresult, istep )
      call hecmw_visualize_finalize
      call hecmw2fstr_mesh_conv(hecmesh)
      call hecmw_result_free(fstrresult)
    end do
  end subroutine
 
  subroutine fstr_eigen_output_result(hecMESH, fstrEIG)
    use m_fstr
    use m_hecmw2fstr_mesh_conv
    use hecmw_util
    implicit none
    type(hecmwst_local_mesh)  :: hecMESH
    type(fstr_eigen)          :: fstrEIG
 
    integer(kind=kint) :: NP, NDOF, NPNDOF, istep, nget, i
    character(len=HECMW_HEADER_LEN) :: header
    character(len=HECMW_MSG_LEN)    :: comment
    character(len=HECMW_NAME_LEN)   :: label
    character(len=HECMW_NAME_LEN)   :: nameID
 
    real(kind=kreal), pointer :: eigvec(:,:)
    real(kind=kreal), allocatable :: x(:), egval(:)
    real(kind=kreal), allocatable :: disp3(:), rot3(:)
 
    nget   = fstreig%nget
    np     = hecmesh%n_node
    ndof   = hecmesh%n_dof
    npndof = np * ndof
 
    eigvec => fstreig%eigvec
    allocate(x(npndof))
    x = 0.0d0
    allocate(egval(1))
 
    do istep = 1, nget
      egval(1) = fstreig%eigval(istep)
      do i=1,npndof
        x(i) = eigvec(i,istep)
      enddo
      call hecmw_update_r(hecmesh, x, np, ndof)
 
      header = "*fstrresult"
      comment = "eigen_result"
      call hecmw_result_init(hecmesh,istep,header,comment)
      label = "EIGENVALUE"
      call hecmw_result_add(hecmw_result_dtype_global,1,label,egval)
      if( ndof == 6 ) then
        allocate(disp3(np*3))
        allocate(rot3(np*3))
        disp3 = 0.0d0
        rot3  = 0.0d0
        do i = 1, np
          disp3((i-1)*3+1:(i-1)*3+3) = x((i-1)*ndof+1:(i-1)*ndof+3)
          rot3((i-1)*3+1:(i-1)*3+3) = x((i-1)*ndof+4:(i-1)*ndof+6)
        enddo
        label = "DISPLACEMENT"
        call hecmw_result_add(hecmw_result_dtype_node,3,label,disp3)
        label = "ROTATION"
        call hecmw_result_add(hecmw_result_dtype_node,3,label,rot3)
        deallocate(disp3)
        deallocate(rot3)
      else
        label = "DISPLACEMENT"
        call hecmw_result_add(hecmw_result_dtype_node,ndof,label,x)
      endif
      nameid = "fstrRES"
      call hecmw_result_write_by_name(nameid)
      call hecmw_result_finalize
    end do
    deallocate(egval)
    deallocate(x)
  end subroutine
 
  subroutine eglist(hecMESH, hecMAT, fstrEIG)
    use m_fstr
    use hecmw_util
    use m_fstr_eig_lanczos_util
    implicit none
    type(hecmwst_local_mesh) :: hecmesh
    type(hecmwst_matrix)     :: hecMAT
    type(fstr_eigen)         :: fstrEIG
 
    integer(kind=kint) :: NDOF
    integer(kind=kint) :: i, j, in, iter ,nget
    real(kind=kreal)   :: pi, angle, freq, pf(3), em(3)
    real(kind=kreal), pointer :: eigval(:)
 
    ndof = hecmat%NDOF
    nget = fstreig%nget
    iter = fstreig%iter
    pi   = 4.0d0 * datan(1.0d0)
    eigval => fstreig%eigval
 
    if(myrank == 0)then
      write(ilog,*)""
      write(ilog,"(a)")"********************************"
      write(ilog,"(a)")"*RESULT OF EIGEN VALUE ANALYSIS*"
      write(ilog,"(a)")"********************************"
      write(ilog,"(a)")""
      write(ilog,"(a,i8)")"NUMBER OF ITERATIONS = ",iter
      write(ilog,"(a,1pe12.4)")"TOTAL MASS = ",fstreig%totalmass
      write(ilog,"(a)")""
      write(ilog,"(3a)")"                   ANGLE       FREQUENCY   ",&
        "PARTICIPATION FACTOR                EFFECTIVE MASS"
      write(ilog,"(3a)")"  NO.  EIGENVALUE  FREQUENCY   (HZ)        ",&
        "X           Y           Z           X           Y           Z"
      write(ilog,"(3a)")"  ---  ----------  ----------  ----------  ",&
        "----------  ----------  ----------  ----------  ----------  ----------"
      write(*,*)""
      write(*,"(a)")"#----------------------------------#"
      write(*,"(a)")"#  RESULT OF EIGEN VALUE ANALYSIS  #"
      write(*,"(a)")"#----------------------------------#"
      write(*,"(a)")""
      write(*,"(a,i8)")"### NUMBER OF ITERATIONS = ",iter
      write(*,"(a,1pe12.4)")"### TOTAL MASS = ",fstreig%totalmass
      write(*,"(a)")""
      write(*,"(3a)")"       PERIOD     FREQUENCY  ",&
        "PARTICIPATION FACTOR             EFFECTIVE MASS"
      write(*,"(3a)")"  NO.  [Sec]      [HZ]       ",&
        "X          Y          Z          X          Y          Z"
      write(*,"(3a)")"  ---  ---------  ---------  ",&
        "---------  ---------  ---------  ---------  ---------  ---------"
 
      in = 0
      do i = 1, nget
        in = in + 1
        if(eigval(i) < 0.0d0) eigval(i) = 0.0d0
        angle = dsqrt(eigval(i))
        freq = angle*0.5d0/pi
 
        pf = 0.0d0
        em = 0.0d0
        do j = 1, min(ndof, 3)
          pf(j) = fstreig%partfactor(ndof*(i-1) + j)
          em(j) = fstreig%effmass(ndof*(i-1) + j)
        enddo
 
        write(ilog,'(I5,1P9E12.4)') in, eigval(i), angle, freq, pf(1), pf(2), pf(3), em(1), em(2), em(3)
        write(*   ,'(I5,1P8E11.3)') in, 1.0d0/freq, freq, pf(1), pf(2), pf(3), em(1), em(2), em(3)
      enddo
      write(ilog,*)""
      write(*,*)""
    endif
 
  end subroutine eglist
end module m_fstr_eig_output