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