1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5!> This module provides linear equation solver interface of MUMPS for
6!! contact problems using Lagrange multiplier.
7module m_solve_LINEQ_MUMPS_contact
8  use m_fstr
9  use m_sparse_matrix
10  use m_sparse_matrix_contact
11  use fstr_matrix_con_contact
12  use m_hecmw_MUMPS_wrapper
13
14  private
15  public :: solve_LINEQ_MUMPS_contact_init
16  public :: solve_LINEQ_MUMPS_contact
17
18  logical, save :: INITIALIZED = .false.
19  logical, save :: NEED_ANALYSIS = .true.
20  type (sparse_matrix), save :: spMAT
21
22contains
23
24  subroutine solve_LINEQ_MUMPS_contact_init(hecMESH,hecMAT,fstrMAT,is_sym)
25    implicit none
26    type (hecmwST_local_mesh), intent(in) :: hecMESH
27    type (hecmwST_matrix    ), intent(inout) :: hecMAT
28    type (fstrST_matrix_contact_lagrange), intent(in) :: fstrMAT !< type fstrST_matrix_contact_lagrange
29    logical, intent(in) :: is_sym
30    integer(kind=kint) :: spmat_type
31    integer(kind=kint) :: spmat_symtype
32    integer(kind=kint) :: mumps_job
33    integer(kind=kint) :: istat
34
35    if (INITIALIZED) then
36      mumps_job=-2
37      call hecmw_mumps_wrapper(spMAT, mumps_job, istat)
38      if (istat < 0) then
39        write(*,*) 'ERROR: MUMPS returned with error', istat
40        stop
41      endif
42      call sparse_matrix_finalize(spMAT)
43      INITIALIZED = .false.
44    endif
45
46    if (is_sym) then
47      !spmat_symtype = SPARSE_MATRIX_SYMTYPE_SPD
48      spmat_symtype = SPARSE_MATRIX_SYMTYPE_SYM
49    else
50      spmat_symtype = SPARSE_MATRIX_SYMTYPE_ASYM
51    endif
52    spmat_type = SPARSE_MATRIX_TYPE_COO
53    call sparse_matrix_set_type(spMAT, spmat_type, spmat_symtype)
54    mumps_job=-1
55    call hecmw_mumps_wrapper(spMAT, mumps_job, istat)
56    if (istat < 0) then
57      write(*,*) 'ERROR: MUMPS returned with error', istat
58      stop
59    endif
60
61    NEED_ANALYSIS = .true.
62    INITIALIZED = .true.
63  end subroutine solve_LINEQ_MUMPS_contact_init
64
65  subroutine solve_LINEQ_MUMPS_contact(hecMESH,hecMAT,fstrMAT,istat,conMAT)
66    implicit none
67    type (hecmwST_local_mesh), intent(in) :: hecMESH
68    type (hecmwST_matrix    ), intent(inout) :: hecMAT
69    type (fstrST_matrix_contact_lagrange), intent(inout) :: fstrMAT !< type fstrST_matrix_contact_lagrange
70    integer(kind=kint), intent(out) :: istat
71    type (hecmwST_matrix), intent(in),optional :: conMAT
72
73    integer(kind=kint) :: mumps_job
74
75    call hecmw_mat_dump(hecMAT, hecMESH)
76
77    ! ANALYSIS
78    if (NEED_ANALYSIS) then
79      call sparse_matrix_contact_init_prof(spMAT, hecMAT, fstrMAT, hecMESH)
80      mumps_job=1
81      call hecmw_mumps_wrapper(spMAT, mumps_job, istat)
82      if (istat < 0) then
83        write(*,*) 'ERROR: MUMPS returned with error', istat
84        stop
85      endif
86      if (myrank==0) write(*,*) ' [MUMPS]: Analysis completed.'
87      NEED_ANALYSIS = .false.
88    endif
89
90    ! FACTORIZATION and SOLUTION
91    !  ----  For Parallel Contact with Multi-Partition Domains
92    if(paraContactFlag.and.present(conMAT)) then
93      call sparse_matrix_para_contact_set_vals(spMAT, hecMAT, fstrMAT, conMAT)
94      call sparse_matrix_para_contact_set_rhs(spMAT, hecMAT, fstrMAT, conMAT)
95    else
96      call sparse_matrix_contact_set_vals(spMAT, hecMAT, fstrMAT)
97      !call sparse_matrix_dump(spMAT)
98      call sparse_matrix_contact_set_rhs(spMAT, hecMAT, fstrMAT)
99    endif
100    mumps_job=5
101    call hecmw_mumps_wrapper(spMAT, mumps_job, istat)
102    if (istat < 0) then
103      write(*,*) 'ERROR: MUMPS returned with error', istat
104      return
105    endif
106    call sparse_matrix_contact_get_rhs(spMAT, hecMAT, fstrMAT)
107    if (myrank==0) write(*,*) ' [MUMPS]: Factorization and Solution completed.'
108
109    call hecmw_mat_dump_solution(hecMAT)
110
111    !call sparse_matrix_finalize(spMAT)
112  end subroutine solve_LINEQ_MUMPS_contact
113
114end module m_solve_LINEQ_MUMPS_contact
115