mform/m_mform_solver.f90
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | !> @brief Module containing the linear solver for linear systems of equations where the unknown is an m-form | ||
| 2 | module m_mform_solver | ||
| 3 | use m_common, only: wp, to_lowercase | ||
| 4 | use m_mform_basis, only: MFormSpace | ||
| 5 | use m_mform_constraint_global, only: MFormConstraint | ||
| 6 | #include "petsc.fi" | ||
| 7 | |||
| 8 | implicit none | ||
| 9 | private | ||
| 10 | |||
| 11 | public :: AbstractSolver, GenericSolver, get_constrained_petsc_matrix, create_ksp | ||
| 12 | public :: SolverInitOptions, SolverSolveOptions, SolverInitTraits | ||
| 13 | |||
| 14 | type SolverSolveOptions | ||
| 15 | !> Relative tolerance for the solver | ||
| 16 | real(wp) :: rtol = -1._wp | ||
| 17 | !> Absolute tolerance for the solver | ||
| 18 | real(wp) :: atol = -1._wp | ||
| 19 | !> Maximum number of iterations for the solver | ||
| 20 | integer :: max_iter = -1 | ||
| 21 | !> Verbosity level for the solver (default: VERBOSITY_WARN_ON_FAILURE) | ||
| 22 | integer :: verbosity = -1 | ||
| 23 | contains | ||
| 24 | procedure :: apply_default_solve_options | ||
| 25 | end type SolverSolveOptions | ||
| 26 | |||
| 27 | type, extends(SolverSolveOptions) :: SolverInitOptions | ||
| 28 | !> KSP type for the PETSc solver (e.g., 'cg', 'gmres', 'minres', etc.) | ||
| 29 | character(len=80) :: ksp_type = '' | ||
| 30 | !> PC type for the PETSc preconditioner (e.g., 'jacobi', 'ilu', 'hypre', etc.) | ||
| 31 | character(len=80) :: pc_type = '' | ||
| 32 | !> HYPRE preconditioner type (e.g., 'boomeramg', 'ams', 'ilu', etc.) | ||
| 33 | character(len=80) :: hypre_type = '' | ||
| 34 | !> ILU level for the HYPRE ILU preconditioner (e.g., '0', '1', '2', '3', or 'ilut') | ||
| 35 | character(len=80) :: ilu_level = '' | ||
| 36 | contains | ||
| 37 | procedure :: init_from_traits => init_solver_options_from_traits | ||
| 38 | end type SolverInitOptions | ||
| 39 | |||
| 40 | type SolverInitTraits | ||
| 41 | logical :: is_symmetric = .false. | ||
| 42 | logical :: is_spd = .false. | ||
| 43 | logical :: is_parallel = .false. | ||
| 44 | logical :: is_diffdiff = .false. | ||
| 45 | logical :: allow_ams = .true. | ||
| 46 | integer :: m_form = -1 | ||
| 47 | end type SolverInitTraits | ||
| 48 | |||
| 49 | interface SolverSolveOptions | ||
| 50 | module procedure make_solver_solve_options | ||
| 51 | end interface | ||
| 52 | |||
| 53 | interface SolverInitOptions | ||
| 54 | module procedure make_solver_options | ||
| 55 | end interface | ||
| 56 | |||
| 57 | !> @brief Abstract type for a generic linear solver for m-form systems of equations | ||
| 58 | type, abstract :: AbstractSolver | ||
| 59 | !> The matrix to be inverted | ||
| 60 | Mat :: mat | ||
| 61 | |||
| 62 | !> KSP object for the PETSc matrix | ||
| 63 | KSP :: mat_ksp | ||
| 64 | |||
| 65 | !> The number of iterations in the last solve | ||
| 66 | integer :: last_nr_iterations | ||
| 67 | !> The residual norm in the last solve | ||
| 68 | real(wp) :: last_residual_norm | ||
| 69 | !> The convergence reason in the last solve | ||
| 70 | KSPConvergedReason :: last_converged_reason | ||
| 71 | |||
| 72 | !> Resolved init options used to configure the solver backend | ||
| 73 | type(SolverInitOptions) :: options | ||
| 74 | |||
| 75 | contains | ||
| 76 | procedure :: get_warning_message | ||
| 77 | procedure :: destroy => destroy_abstract_solver | ||
| 78 | end type | ||
| 79 | |||
| 80 | !> @brief Generic solver for m-form systems of equations | ||
| 81 | type, extends(AbstractSolver) :: GenericSolver | ||
| 82 | !> The m-form space of the solution vector | ||
| 83 | type(MFormSpace) :: space_col | ||
| 84 | |||
| 85 | !> The m-form space of the residual | ||
| 86 | type(MFormSpace) :: space_row | ||
| 87 | |||
| 88 | !> The optional constraint | ||
| 89 | type(MFormConstraint), allocatable :: constraint | ||
| 90 | |||
| 91 | contains | ||
| 92 | procedure :: init => init_generic_solver | ||
| 93 | procedure :: solve => solve_generic_solver | ||
| 94 | procedure :: destroy => destroy_generic_solver | ||
| 95 | end type | ||
| 96 | |||
| 97 | contains | ||
| 98 | |||
| 99 | !> @brief Create a SolverSolveOptions object from the provided arguments | ||
| 100 | !> | ||
| 101 | !> @param[in] rtol _(optional)_ The relative tolerance for the solver (default is 1.e-10, or the value set in the init() method) | ||
| 102 | !> @param[in] atol _(optional)_ The absolute tolerance for the solver (default is 1.e-20, or the value set in the init() method) | ||
| 103 | !> @param[in] max_iter _(optional)_ The maximum number of iterations for the solver (default is the value set in the init() method) | ||
| 104 | !> @param[in] verbosity _(optional)_ The verbosity level for the solver (default is the value set in the init() method) | ||
| 105 | !> @param[in] verbosity_flag _(optional)_ The verbosity flag for the solver (default is the value set in the init() method) | ||
| 106 | 134 | function make_solver_solve_options(rtol, atol, max_iter, verbosity, verbosity_flag) result(options) | |
| 107 | type(SolverSolveOptions) :: options | ||
| 108 | real(wp), optional, intent(in) :: rtol, atol | ||
| 109 | integer, optional, intent(in) :: max_iter, verbosity_flag | ||
| 110 | character(*), optional, intent(in) :: verbosity | ||
| 111 | |||
| 112 | 134 | if (present(rtol)) options%rtol = rtol | |
| 113 |
1/2✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 134 times.
|
134 | if (present(atol)) options%atol = atol |
| 114 |
1/2✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 134 times.
|
134 | if (present(max_iter)) options%max_iter = max_iter |
| 115 |
1/2✗ Branch 8 → 9 not taken.
✓ Branch 8 → 10 taken 134 times.
|
134 | if (present(verbosity)) options%verbosity = solver_option_verbosity_to_flag(verbosity, options%verbosity) |
| 116 |
1/2✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 134 times.
|
134 | if (present(verbosity_flag)) options%verbosity = verbosity_flag |
| 117 |
1/2✓ Branch 2 → 3 taken 134 times.
✗ Branch 2 → 4 not taken.
|
134 | end function make_solver_solve_options |
| 118 | |||
| 119 | !> @brief Create a SolverInitOptions object from the provided arguments | ||
| 120 | !> | ||
| 121 | !> @param[in] ksp_type _(optional)_ The type of the KSP solver to be used (default is 'gmres' for non-symmetric matrices, 'cg' for | ||
| 122 | !> symmetric positive definite matrices, and 'minres' for symmetric indefinite matrices) | ||
| 123 | !> @param[in] pc_type _(optional)_ The type of the preconditioner to be used (default is 'hypre' for parallel problems, 'ilu' for | ||
| 124 | !> serial problems, and 'hypre' for any curl-curl problems) | ||
| 125 | !> @param[in] hypre_type _(optional)_ The type of the Hypre preconditioner to be used (default is 'boomeramg' for parallel | ||
| 126 | !> problems, 'ams' for curl-curl problems; note that 'pc_type' is set to 'hypre' if 'hypre_type' is present | ||
| 127 | !> and 'pc_type' is not). For ILU the following options are supported: 'ilu' (for ILUT), 'ilu(0)', 'ilu(1)', | ||
| 128 | !> 'ilu(2)', 'ilu(3)' (for ILUK). | ||
| 129 | !> @param[in] rtol _(optional)_ The relative tolerance for the solver (default is 1.e-10, or the value set in the init() method) | ||
| 130 | !> @param[in] atol _(optional)_ The absolute tolerance for the solver (default is 1.e-20, or the value set in the init() method) | ||
| 131 | !> @param[in] max_iter _(optional)_ The maximum number of iterations for the solver (default is the value set in the init() method) | ||
| 132 | !> @param[in] verbosity _(optional)_ The verbosity level for the solver (default is the value set in the init() method) | ||
| 133 | !> @param[in] verbosity_flag _(optional)_ The verbosity flag for the solver (default is the value set in the init() method) | ||
| 134 | 879 | function make_solver_options(ksp_type, pc_type, hypre_type, rtol, atol, max_iter, verbosity, & | |
| 135 | verbosity_flag) result(options) | ||
| 136 | type(SolverInitOptions) :: options | ||
| 137 | character(*), optional, intent(in) :: ksp_type, pc_type, hypre_type | ||
| 138 | real(wp), optional, intent(in) :: rtol, atol | ||
| 139 | integer, optional, intent(in) :: max_iter, verbosity_flag | ||
| 140 | character(*), optional, intent(in) :: verbosity | ||
| 141 | |||
| 142 | character(len=80) :: hypre_type_ | ||
| 143 | integer :: ierr | ||
| 144 | |||
| 145 |
1/2✓ Branch 3 → 4 taken 6 times.
✗ Branch 3 → 5 not taken.
|
6 | if (present(ksp_type)) options%ksp_type = ksp_type |
| 146 |
1/4✗ Branch 6 → 7 not taken.
✓ Branch 6 → 10 taken 879 times.
✗ Branch 7 → 8 not taken.
✗ Branch 7 → 9 not taken.
|
879 | if (present(pc_type)) options%pc_type = pc_type |
| 147 |
2/2✓ Branch 10 → 11 taken 26 times.
✓ Branch 10 → 55 taken 853 times.
|
879 | if (present(hypre_type)) then |
| 148 |
1/2✓ Branch 11 → 12 taken 26 times.
✗ Branch 11 → 13 not taken.
|
26 | if (.not. present(pc_type)) options%pc_type = 'hypre' |
| 149 | |||
| 150 |
2/4✓ Branch 16 → 17 taken 26 times.
✗ Branch 16 → 18 not taken.
✓ Branch 18 → 19 taken 26 times.
✗ Branch 18 → 20 not taken.
|
26 | hypre_type_ = to_lowercase(trim(adjustl(hypre_type))) |
| 151 | 26 | select case (trim(hypre_type_)) | |
| 152 | case ('ilu', 'ilu(0)') | ||
| 153 | ✗ | options%ilu_level = '0' | |
| 154 | case ('ilu(1)') | ||
| 155 | ✗ | options%ilu_level = '1' | |
| 156 | case ('ilu(2)') | ||
| 157 | ✗ | options%ilu_level = '2' | |
| 158 | case ('ilu(3)') | ||
| 159 | ✗ | options%ilu_level = '3' | |
| 160 | case default | ||
| 161 |
1/5✗ Branch 21 → 22 not taken.
✗ Branch 21 → 23 not taken.
✗ Branch 21 → 24 not taken.
✗ Branch 21 → 25 not taken.
✓ Branch 21 → 26 taken 26 times.
|
26 | options%ilu_level = 'ilut' |
| 162 | end select | ||
| 163 |
1/2✗ Branch 27 → 28 not taken.
✓ Branch 27 → 29 taken 26 times.
|
26 | if (hypre_type_(1:3) == 'ilu') hypre_type_ = 'ilu' |
| 164 |
1/2✓ Branch 29 → 30 taken 26 times.
✗ Branch 29 → 31 not taken.
|
26 | options%hypre_type = hypre_type |
| 165 | |||
| 166 |
2/4✓ Branch 32 → 33 taken 26 times.
✗ Branch 32 → 55 not taken.
✗ Branch 33 → 34 not taken.
✓ Branch 33 → 55 taken 26 times.
|
26 | if (options%pc_type == 'hypre' .and. options%hypre_type == 'ilu') then |
| 167 | ✗ | PetscCall(PetscOptionsSetValue(PETSC_NULL_OPTIONS, '-pc_type', options%pc_type, ierr)) | |
| 168 | ✗ | PetscCall(PetscOptionsSetValue(PETSC_NULL_OPTIONS, '-pc_hypre_type', options%hypre_type, ierr)) | |
| 169 | |||
| 170 | ✗ | if (options%ilu_level == 'ilut') then | |
| 171 | ✗ | PetscCall(PetscOptionsSetValue(PETSC_NULL_OPTIONS, '-pc_hypre_ilu_type', 'Block-Jacobi-ILUT', ierr)) | |
| 172 | ! PetscCall(PetscOptionsSetValue(PETSC_NULL_OPTIONS, '-pc_hypre_ilu_drop_threshold', '1E-3', ierr)) | ||
| 173 | else | ||
| 174 | ✗ | PetscCall(PetscOptionsSetValue(PETSC_NULL_OPTIONS, '-pc_hypre_ilu_type', 'Block-Jacobi-ILUK', ierr)) | |
| 175 | ✗ | PetscCall(PetscOptionsSetValue(PETSC_NULL_OPTIONS, '-pc_hypre_ilu_level', options%ilu_level, ierr)) | |
| 176 | end if | ||
| 177 | end if | ||
| 178 | end if | ||
| 179 | |||
| 180 |
2/2✓ Branch 55 → 56 taken 205 times.
✓ Branch 55 → 57 taken 674 times.
|
879 | if (present(rtol)) options%rtol = rtol |
| 181 |
1/2✗ Branch 57 → 58 not taken.
✓ Branch 57 → 59 taken 879 times.
|
879 | if (present(atol)) options%atol = atol |
| 182 |
1/2✗ Branch 59 → 60 not taken.
✓ Branch 59 → 61 taken 879 times.
|
879 | if (present(max_iter)) options%max_iter = max_iter |
| 183 |
1/2✗ Branch 61 → 62 not taken.
✓ Branch 61 → 63 taken 879 times.
|
879 | if (present(verbosity)) options%verbosity = solver_option_verbosity_to_flag(verbosity, options%verbosity) |
| 184 |
1/2✗ Branch 63 → 64 not taken.
✓ Branch 63 → 65 taken 879 times.
|
879 | if (present(verbosity_flag)) options%verbosity = verbosity_flag |
| 185 |
2/2✓ Branch 2 → 3 taken 6 times.
✓ Branch 2 → 6 taken 873 times.
|
1758 | end function make_solver_options |
| 186 | |||
| 187 | !> @brief Initialize the generic solver for a linear system of equations defined on m-forms | ||
| 188 | !> | ||
| 189 | !> @param[inout] this The generic solver object to be initialized | ||
| 190 | !> @param[in] amat The matrix to be inverted (any type derived from AbstractMatrix) | ||
| 191 | !> @param[in] constraint1 _(optional)_ An optional constraint to be applied to the solution | ||
| 192 | !> @param[in] constraint2 _(optional)_ An optional constraint to be applied to the solution | ||
| 193 | !> @param[in] constraint3 _(optional)_ An optional constraint to be applied to the solution | ||
| 194 | !> @param[in] coord_transform _(optional)_ The coordinate transformation associated with the m-form spaces | ||
| 195 | !> @param[in] options _(optional)_ Options for initializing the solver | ||
| 196 | 618 | subroutine init_generic_solver(this, amat, constraint1, constraint2, constraint3, coord_transform, options) | |
| 197 | use m_common, only: VERBOSITY_WARN_ON_FAILURE | ||
| 198 | use m_tensorprod_matrix, only: get_petsc | ||
| 199 | use m_mform_basis, only: size | ||
| 200 | use m_mform_constraint_abstract, only: MFormConstraintLocal | ||
| 201 | use m_mform_matrix, only: AbstractMatrix, DiffDiffMatrix | ||
| 202 | use m_coord_transform_abstract, only: CoordTransformAbstract | ||
| 203 | |||
| 204 | implicit none | ||
| 205 | |||
| 206 | class(GenericSolver), intent(inout) :: this | ||
| 207 | class(AbstractMatrix), intent(in) :: amat | ||
| 208 | class(MFormConstraintLocal), optional, intent(in) :: constraint1, constraint2, constraint3 | ||
| 209 | class(CoordTransformAbstract), optional, intent(in) :: coord_transform | ||
| 210 | type(SolverInitOptions), optional, intent(in) :: options | ||
| 211 | |||
| 212 | Mat :: mat_tmp | ||
| 213 | type(SolverInitTraits) :: traits | ||
| 214 | integer :: ierr | ||
| 215 | logical :: is_diffdiff | ||
| 216 | |||
| 217 | 618 | call this%destroy() | |
| 218 | ! PETSCViewer :: viewer | ||
| 219 | |||
| 220 |
1/2✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 618 times.
|
618 | if (amat%space_col%m /= amat%space_row%m) then |
| 221 | ✗ | error stop "GenericSolver::init_generic_solver: Non-square matrix provided which is not supported by this solver." | |
| 222 | end if | ||
| 223 | |||
| 224 |
5/12✓ Branch 5 → 6 taken 350 times.
✓ Branch 5 → 7 taken 268 times.
✗ Branch 6 → 7 not taken.
✓ Branch 6 → 11 taken 350 times.
✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 268 times.
✗ Branch 8 → 9 not taken.
✗ Branch 8 → 11 not taken.
✗ Branch 9 → 10 not taken.
✓ Branch 9 → 33 taken 268 times.
✗ Branch 10 → 11 not taken.
✗ Branch 10 → 33 not taken.
|
618 | if (present(constraint1) .or. present(constraint2) .or. present(constraint3)) then |
| 225 |
4/6✗ Branch 11 → 12 not taken.
✓ Branch 11 → 13 taken 350 times.
✗ Branch 13 → 14 not taken.
✓ Branch 13 → 15 taken 350 times.
✓ Branch 16 → 17 taken 1750 times.
✓ Branch 16 → 18 taken 350 times.
|
2100 | allocate (this%constraint) |
| 226 | 350 | call this%constraint%init(amat%space_col) | |
| 227 | |||
| 228 |
2/4✓ Branch 19 → 20 taken 350 times.
✗ Branch 19 → 23 not taken.
✓ Branch 20 → 21 taken 350 times.
✗ Branch 20 → 23 not taken.
|
350 | if (present(constraint1)) call this%constraint%insert(constraint1, coord_transform=coord_transform) |
| 229 |
3/4✓ Branch 23 → 24 taken 159 times.
✓ Branch 23 → 27 taken 191 times.
✓ Branch 24 → 25 taken 159 times.
✗ Branch 24 → 27 not taken.
|
350 | if (present(constraint2)) call this%constraint%insert(constraint2, coord_transform=coord_transform) |
| 230 |
1/4✗ Branch 27 → 28 not taken.
✓ Branch 27 → 31 taken 350 times.
✗ Branch 28 → 29 not taken.
✗ Branch 28 → 31 not taken.
|
350 | if (present(constraint3)) call this%constraint%insert(constraint3, coord_transform=coord_transform) |
| 231 | |||
| 232 | 350 | call this%constraint%assemble() | |
| 233 | end if | ||
| 234 | |||
| 235 |
5/8✓ Branch 33 → 34 taken 618 times.
✗ Branch 33 → 37 not taken.
✓ Branch 34 → 35 taken 618 times.
✗ Branch 34 → 36 not taken.
✓ Branch 37 → 38 taken 618 times.
✗ Branch 37 → 40 not taken.
✓ Branch 38 → 39 taken 403 times.
✓ Branch 38 → 40 taken 215 times.
|
618 | this%space_col = amat%space_col |
| 236 |
5/8✓ Branch 40 → 41 taken 618 times.
✗ Branch 40 → 44 not taken.
✓ Branch 41 → 42 taken 618 times.
✗ Branch 41 → 43 not taken.
✓ Branch 44 → 45 taken 618 times.
✗ Branch 44 → 47 not taken.
✓ Branch 45 → 46 taken 403 times.
✓ Branch 45 → 47 taken 215 times.
|
618 | this%space_row = amat%space_row |
| 237 | |||
| 238 | is_diffdiff = .false. | ||
| 239 | select type (amat) | ||
| 240 | type is (DiffDiffMatrix) | ||
| 241 | is_diffdiff = .true. | ||
| 242 | end select | ||
| 243 | |||
| 244 | 618 | this%options = SolverInitOptions() | |
| 245 |
2/2✓ Branch 50 → 51 taken 223 times.
✓ Branch 50 → 52 taken 395 times.
|
618 | if (present(options)) this%options = options |
| 246 | |||
| 247 | 618 | traits%is_symmetric = amat%is_symmetric | |
| 248 | 618 | traits%is_spd = amat%is_spd | |
| 249 | 618 | traits%is_parallel = amat%space_row%tp_spaces(1)%domain%nr_ranks > 1 | |
| 250 | 618 | traits%is_diffdiff = is_diffdiff | |
| 251 | 618 | traits%m_form = amat%space_col%m | |
| 252 |
4/4✓ Branch 52 → 53 taken 259 times.
✓ Branch 52 → 54 taken 359 times.
✓ Branch 53 → 54 taken 237 times.
✓ Branch 53 → 55 taken 22 times.
|
618 | traits%allow_ams = is_diffdiff .and. amat%space_col%m == 1 |
| 253 | |||
| 254 | 618 | call this%options%init_from_traits(traits) | |
| 255 | 618 | call this%options%apply_default_solve_options(size(amat%space_col), VERBOSITY_WARN_ON_FAILURE) | |
| 256 | |||
| 257 |
2/2✓ Branch 57 → 58 taken 268 times.
✓ Branch 57 → 62 taken 350 times.
|
618 | if (.not. allocated(this%constraint)) then |
| 258 | ! NOTE a copy is required because the destructor of amat will destroy the matrix | ||
| 259 |
1/2✗ Branch 59 → 60 not taken.
✓ Branch 59 → 65 taken 268 times.
|
268 | PetscCall(MatConvert(amat%mat, "same", MAT_INITIAL_MATRIX, this%mat, ierr)) |
| 260 | else | ||
| 261 | 350 | call get_constrained_petsc_matrix(this%mat, amat%mat, this%constraint, coeff=amat%constraint_coeff()) | |
| 262 | end if | ||
| 263 | |||
| 264 |
2/2✓ Branch 65 → 66 taken 350 times.
✓ Branch 65 → 67 taken 268 times.
|
618 | if (allocated(this%constraint)) then |
| 265 | 350 | call create_ksp(this%mat_ksp, this%mat, this%space_col, this%options, constraint=this%constraint) | |
| 266 | else | ||
| 267 | 268 | call create_ksp(this%mat_ksp, this%mat, this%space_col, this%options) | |
| 268 | end if | ||
| 269 | |||
| 270 | 618 | end subroutine init_generic_solver | |
| 271 | |||
| 272 | 626 | subroutine create_ksp(ksp, mat, space, options, constraint) | |
| 273 | use m_mform_constraint_global, only: MFormConstraint | ||
| 274 | |||
| 275 | KSP, intent(out) :: ksp | ||
| 276 | Mat, intent(in) :: mat | ||
| 277 | type(MFormSpace), intent(in) :: space | ||
| 278 | type(SolverInitOptions), intent(in) :: options | ||
| 279 | type(MFormConstraint), intent(in), optional :: constraint | ||
| 280 | |||
| 281 | PC :: pc | ||
| 282 | integer :: ierr | ||
| 283 | |||
| 284 |
1/2✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 626 times.
|
626 | PetscCall(KSPCreate(space%tp_spaces(1)%domain%comm, ksp, ierr)) |
| 285 |
1/2✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 626 times.
|
626 | PetscCall(KSPSetOperators(ksp, mat, mat, ierr)) |
| 286 |
1/2✗ Branch 9 → 10 not taken.
✓ Branch 9 → 11 taken 626 times.
|
626 | PetscCall(KSPSetFromOptions(ksp, ierr)) |
| 287 |
1/2✗ Branch 12 → 13 not taken.
✓ Branch 12 → 14 taken 626 times.
|
626 | PetscCall(KSPSetType(ksp, options%ksp_type, ierr)) |
| 288 | |||
| 289 |
1/2✗ Branch 15 → 16 not taken.
✓ Branch 15 → 17 taken 626 times.
|
626 | PetscCall(KSPGetPC(ksp, pc, ierr)) |
| 290 | |||
| 291 |
3/4✓ Branch 17 → 18 taken 30 times.
✓ Branch 17 → 19 taken 596 times.
✓ Branch 18 → 19 taken 30 times.
✗ Branch 18 → 32 not taken.
|
626 | if (.not. (options%pc_type == 'hypre' .and. options%hypre_type == 'ilu')) then |
| 292 | ! TODO: two ways of settings options seem incompatible? | ||
| 293 |
1/2✗ Branch 20 → 21 not taken.
✓ Branch 20 → 22 taken 626 times.
|
626 | PetscCall(PCSetType(pc, options%pc_type, ierr)) |
| 294 |
2/2✓ Branch 22 → 23 taken 30 times.
✓ Branch 22 → 26 taken 596 times.
|
626 | if (options%pc_type == 'hypre') then |
| 295 |
1/2✗ Branch 24 → 25 not taken.
✓ Branch 24 → 26 taken 30 times.
|
30 | PetscCall(PCHYPRESetType(pc, options%hypre_type, ierr)) |
| 296 | end if | ||
| 297 | if ((options%pc_type == 'cholesky' .or. options%pc_type == 'lu') & | ||
| 298 |
2/6✓ Branch 26 → 27 taken 626 times.
✗ Branch 26 → 28 not taken.
✗ Branch 27 → 28 not taken.
✓ Branch 27 → 32 taken 626 times.
✗ Branch 28 → 29 not taken.
✗ Branch 28 → 32 not taken.
|
626 | & .and. space%tp_spaces(1)%domain%nr_ranks > 1) then |
| 299 | ✗ | PetscCall(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS, ierr)) | |
| 300 | end if | ||
| 301 | end if | ||
| 302 | |||
| 303 |
4/4✓ Branch 32 → 33 taken 30 times.
✓ Branch 32 → 35 taken 596 times.
✓ Branch 33 → 34 taken 22 times.
✓ Branch 33 → 35 taken 8 times.
|
626 | if (options%pc_type == 'hypre' .and. options%hypre_type == 'ams') call init_ams_pc(pc, space, constraint) |
| 304 | |||
| 305 |
1/2✗ Branch 36 → 37 not taken.
✓ Branch 36 → 38 taken 626 times.
|
626 | PetscCall(KSPSetPC(ksp, pc, ierr)) |
| 306 | 626 | end subroutine create_ksp | |
| 307 | |||
| 308 | 626 | subroutine init_solver_options_from_traits(this, traits) | |
| 309 | class(SolverInitOptions), intent(inout) :: this | ||
| 310 | type(SolverInitTraits), intent(in) :: traits | ||
| 311 | |||
| 312 |
2/2✓ Branch 2 → 3 taken 620 times.
✓ Branch 2 → 12 taken 6 times.
|
626 | if (len_trim(this%ksp_type) == 0) then |
| 313 |
5/6✓ Branch 3 → 4 taken 253 times.
✓ Branch 3 → 7 taken 367 times.
✓ Branch 4 → 5 taken 16 times.
✓ Branch 4 → 7 taken 237 times.
✓ Branch 5 → 6 taken 16 times.
✗ Branch 5 → 7 not taken.
|
620 | if (traits%is_diffdiff .and. traits%m_form == 1 .and. traits%allow_ams) then |
| 314 | 16 | this%ksp_type = 'minres' | |
| 315 |
2/2✓ Branch 7 → 8 taken 596 times.
✓ Branch 7 → 11 taken 8 times.
|
604 | else if (traits%is_symmetric) then |
| 316 |
1/2✓ Branch 8 → 9 taken 596 times.
✗ Branch 8 → 10 not taken.
|
596 | if (traits%is_spd) then |
| 317 | 596 | this%ksp_type = 'cg' | |
| 318 | else | ||
| 319 | ✗ | this%ksp_type = 'minres' | |
| 320 | end if | ||
| 321 | else | ||
| 322 | 8 | this%ksp_type = 'gmres' | |
| 323 | end if | ||
| 324 | end if | ||
| 325 | |||
| 326 |
2/2✓ Branch 12 → 13 taken 600 times.
✓ Branch 12 → 27 taken 26 times.
|
626 | if (len_trim(this%pc_type) == 0) then |
| 327 |
1/2✗ Branch 13 → 14 not taken.
✓ Branch 13 → 15 taken 600 times.
|
600 | if (len_trim(this%hypre_type) > 0) then |
| 328 | ✗ | this%pc_type = 'hypre' | |
| 329 |
1/2✓ Branch 15 → 16 taken 600 times.
✗ Branch 15 → 24 not taken.
|
600 | else if (traits%is_symmetric) then |
| 330 |
6/8✓ Branch 16 → 17 taken 600 times.
✗ Branch 16 → 20 not taken.
✓ Branch 17 → 18 taken 241 times.
✓ Branch 17 → 21 taken 359 times.
✓ Branch 18 → 19 taken 4 times.
✓ Branch 18 → 21 taken 237 times.
✓ Branch 19 → 20 taken 4 times.
✗ Branch 19 → 21 not taken.
|
600 | if (traits%is_parallel .or. (traits%is_diffdiff .and. traits%m_form == 1 .and. traits%allow_ams)) then |
| 331 | 4 | this%pc_type = 'hypre' | |
| 332 |
1/2✓ Branch 21 → 22 taken 596 times.
✗ Branch 21 → 23 not taken.
|
596 | else if (traits%is_spd) then |
| 333 | 596 | this%pc_type = 'icc' | |
| 334 | else | ||
| 335 | ✗ | this%pc_type = 'ilu' | |
| 336 | end if | ||
| 337 | else | ||
| 338 | ✗ | if (traits%is_parallel) then | |
| 339 | ✗ | this%pc_type = 'hypre' | |
| 340 | else | ||
| 341 | ✗ | this%pc_type = 'ilu' | |
| 342 | end if | ||
| 343 | end if | ||
| 344 | end if | ||
| 345 | |||
| 346 |
2/2✓ Branch 27 → 28 taken 600 times.
✓ Branch 27 → 33 taken 26 times.
|
626 | if (len_trim(this%hypre_type) == 0) then |
| 347 |
5/6✓ Branch 28 → 29 taken 241 times.
✓ Branch 28 → 32 taken 359 times.
✓ Branch 29 → 30 taken 4 times.
✓ Branch 29 → 32 taken 237 times.
✓ Branch 30 → 31 taken 4 times.
✗ Branch 30 → 32 not taken.
|
600 | if (traits%is_diffdiff .and. traits%m_form == 1 .and. traits%allow_ams) then |
| 348 | 4 | this%hypre_type = 'ams' | |
| 349 | else | ||
| 350 | 596 | this%hypre_type = 'boomeramg' | |
| 351 | end if | ||
| 352 | end if | ||
| 353 | 626 | end subroutine init_solver_options_from_traits | |
| 354 | |||
| 355 | 626 | subroutine apply_default_solve_options(options, default_max_iter, default_verbosity) | |
| 356 | class(SolverSolveOptions), intent(inout) :: options | ||
| 357 | integer, intent(in) :: default_max_iter | ||
| 358 | integer, intent(in) :: default_verbosity | ||
| 359 | |||
| 360 |
2/2✓ Branch 2 → 3 taken 421 times.
✓ Branch 2 → 4 taken 205 times.
|
626 | if (options%rtol < 0._wp) options%rtol = 1.e-10_wp |
| 361 |
1/2✓ Branch 4 → 5 taken 626 times.
✗ Branch 4 → 6 not taken.
|
626 | if (options%atol < 0._wp) options%atol = 1.e-20_wp |
| 362 |
1/2✓ Branch 6 → 7 taken 626 times.
✗ Branch 6 → 8 not taken.
|
626 | if (options%max_iter < 0) options%max_iter = default_max_iter |
| 363 |
1/2✓ Branch 8 → 9 taken 626 times.
✗ Branch 8 → 10 not taken.
|
626 | if (options%verbosity < 0) options%verbosity = default_verbosity |
| 364 | 626 | end subroutine apply_default_solve_options | |
| 365 | |||
| 366 | ✗ | integer function solver_option_verbosity_to_flag(verbosity, default_flag) result(flag) | |
| 367 | use m_common, only: VERBOSITY_WARN_NEVER, VERBOSITY_WARN_ON_FAILURE, VERBOSITY_WARN_ALWAYS, VERBOSITY_ERROR_ON_FAILURE | ||
| 368 | |||
| 369 | character(*), intent(in) :: verbosity | ||
| 370 | integer, intent(in) :: default_flag | ||
| 371 | |||
| 372 | character(len=64) :: v | ||
| 373 | |||
| 374 | ✗ | v = to_lowercase(trim(adjustl(verbosity))) | |
| 375 | ✗ | select case (trim(v)) | |
| 376 | case ('warn never', 'never', 'none') | ||
| 377 | flag = VERBOSITY_WARN_NEVER | ||
| 378 | case ('warn on failure', 'on failure', 'warn') | ||
| 379 | flag = VERBOSITY_WARN_ON_FAILURE | ||
| 380 | case ('warn always', 'always') | ||
| 381 | flag = VERBOSITY_WARN_ALWAYS | ||
| 382 | case ('error on failure', 'error') | ||
| 383 | ✗ | flag = VERBOSITY_ERROR_ON_FAILURE | |
| 384 | case default | ||
| 385 | ✗ | flag = default_flag | |
| 386 | end select | ||
| 387 | ✗ | end function solver_option_verbosity_to_flag | |
| 388 | |||
| 389 | 984 | subroutine destroy_generic_solver(this) | |
| 390 | class(GenericSolver), intent(inout) :: this | ||
| 391 | |||
| 392 | 984 | call destroy_abstract_solver(this) | |
| 393 | |||
| 394 |
2/2✓ Branch 3 → 4 taken 350 times.
✓ Branch 3 → 78 taken 634 times.
|
984 | if (allocated(this%constraint)) then |
| 395 | 350 | call this%constraint%destroy() | |
| 396 |
12/68✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 350 times.
✓ Branch 7 → 8 taken 350 times.
✗ Branch 7 → 9 not taken.
✓ Branch 10 → 11 taken 1750 times.
✓ Branch 10 → 54 taken 350 times.
✗ Branch 11 → 12 not taken.
✓ Branch 11 → 14 taken 1750 times.
✗ Branch 12 → 13 not taken.
✗ Branch 12 → 14 not taken.
✗ Branch 14 → 15 not taken.
✓ Branch 14 → 53 taken 1750 times.
✗ Branch 15 → 16 not taken.
✗ Branch 15 → 17 not taken.
✗ Branch 17 → 18 not taken.
✗ Branch 17 → 30 not taken.
✗ Branch 19 → 20 not taken.
✗ Branch 19 → 29 not taken.
✗ Branch 20 → 21 not taken.
✗ Branch 20 → 22 not taken.
✗ Branch 22 → 23 not taken.
✗ Branch 22 → 24 not taken.
✗ Branch 24 → 25 not taken.
✗ Branch 24 → 26 not taken.
✗ Branch 26 → 27 not taken.
✗ Branch 26 → 28 not taken.
✗ Branch 30 → 31 not taken.
✗ Branch 30 → 41 not taken.
✗ Branch 32 → 33 not taken.
✗ Branch 32 → 40 not taken.
✗ Branch 33 → 34 not taken.
✗ Branch 33 → 35 not taken.
✗ Branch 35 → 36 not taken.
✗ Branch 35 → 37 not taken.
✗ Branch 37 → 38 not taken.
✗ Branch 37 → 39 not taken.
✗ Branch 41 → 42 not taken.
✗ Branch 41 → 52 not taken.
✗ Branch 43 → 44 not taken.
✗ Branch 43 → 51 not taken.
✗ Branch 44 → 45 not taken.
✗ Branch 44 → 46 not taken.
✗ Branch 46 → 47 not taken.
✗ Branch 46 → 48 not taken.
✗ Branch 48 → 49 not taken.
✗ Branch 48 → 50 not taken.
✗ Branch 54 → 55 not taken.
✓ Branch 54 → 56 taken 350 times.
✗ Branch 56 → 57 not taken.
✓ Branch 56 → 58 taken 350 times.
✗ Branch 58 → 59 not taken.
✓ Branch 58 → 60 taken 350 times.
✗ Branch 60 → 61 not taken.
✓ Branch 60 → 62 taken 350 times.
✗ Branch 62 → 63 not taken.
✓ Branch 62 → 64 taken 350 times.
✗ Branch 64 → 65 not taken.
✓ Branch 64 → 77 taken 350 times.
✗ Branch 66 → 67 not taken.
✗ Branch 66 → 76 not taken.
✗ Branch 67 → 68 not taken.
✗ Branch 67 → 69 not taken.
✗ Branch 69 → 70 not taken.
✗ Branch 69 → 71 not taken.
✗ Branch 71 → 72 not taken.
✗ Branch 71 → 73 not taken.
✗ Branch 73 → 74 not taken.
✗ Branch 73 → 75 not taken.
|
2450 | deallocate (this%constraint) |
| 397 | end if | ||
| 398 | 984 | end subroutine destroy_generic_solver | |
| 399 | |||
| 400 | !> @brief Solve the linear system of equations defined by the initialized generic solver and the right-hand side vector | ||
| 401 | !> | ||
| 402 | !> @param[in] this The generic solver object | ||
| 403 | !> @param[inout] x The solution vector (m-form function, x%distribute() is called) | ||
| 404 | !> @param[in] b The right-hand side vector (m-form function, b%distribute() is not needed) | ||
| 405 | !> @param[in] options _(optional)_ Options for solving the linear system | ||
| 406 | 1209 | subroutine solve_generic_solver(this, x, b, options) | |
| 407 | use m_common, only: wp, VERBOSITY_WARN_NEVER, VERBOSITY_WARN_ON_FAILURE, VERBOSITY_WARN_ALWAYS, VERBOSITY_ERROR_ON_FAILURE | ||
| 408 | use m_mform_basis, only: size, MFormFun, get_petsc | ||
| 409 | use m_mform_constraint_global, only: CONSTRAINT_EXTRACTION, CONSTRAINT_PROJECTION | ||
| 410 | |||
| 411 | implicit none | ||
| 412 | |||
| 413 | class(GenericSolver), intent(in) :: this | ||
| 414 | type(MFormFun), intent(inout) :: x | ||
| 415 | type(MFormFun), intent(in) :: b | ||
| 416 | type(SolverSolveOptions), optional, intent(in) :: options | ||
| 417 | |||
| 418 | integer :: ierr, verbosity_ | ||
| 419 | real(wp) :: rtol_, atol_ | ||
| 420 | integer :: max_iter_ | ||
| 421 | character(len=512) :: warn_msg | ||
| 422 | logical :: failure | ||
| 423 | |||
| 424 | Vec :: x_vec, b_vec, b_proj_vec, x_proj_vec | ||
| 425 | |||
| 426 | 1209 | rtol_ = this%options%rtol | |
| 427 | 1209 | atol_ = this%options%atol | |
| 428 | 1209 | max_iter_ = this%options%max_iter | |
| 429 | 1209 | if (present(options)) then | |
| 430 |
1/2✓ Branch 3 → 4 taken 134 times.
✗ Branch 3 → 5 not taken.
|
134 | if (options%rtol >= 0._wp) rtol_ = options%rtol |
| 431 |
1/2✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 134 times.
|
134 | if (options%atol >= 0._wp) atol_ = options%atol |
| 432 |
1/2✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 134 times.
|
134 | if (options%max_iter >= 0) max_iter_ = options%max_iter |
| 433 | end if | ||
| 434 | |||
| 435 |
1/2✓ Branch 9 → 10 taken 1209 times.
✗ Branch 9 → 13 not taken.
|
1209 | if (this%space_row%tp_spaces(1)%domain%my_rank > 0) then |
| 436 | verbosity_ = VERBOSITY_WARN_NEVER | ||
| 437 | else | ||
| 438 | 1209 | verbosity_ = this%options%verbosity | |
| 439 |
2/2✓ Branch 10 → 11 taken 1075 times.
✓ Branch 10 → 12 taken 134 times.
|
1209 | if (present(options)) then |
| 440 |
1/2✓ Branch 12 → 11 taken 134 times.
✗ Branch 12 → 13 not taken.
|
134 | if (options%verbosity >= 0) verbosity_ = options%verbosity |
| 441 | end if | ||
| 442 | end if | ||
| 443 | |||
| 444 |
1/2✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 1209 times.
|
1209 | PetscCall(KSPSetTolerances(this%mat_ksp, rtol_, atol_, 1.e10_wp, max_iter_, ierr)) |
| 445 | |||
| 446 |
1/2✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 1209 times.
|
1209 | if (this%space_row%m /= b%space%m) then |
| 447 | ✗ | error stop "GenericSolver::solve_generic_solver: m-form spaces do not match" | |
| 448 | end if | ||
| 449 | |||
| 450 | 1209 | call get_petsc(b_vec, b) | |
| 451 |
2/2✓ Branch 19 → 20 taken 535 times.
✓ Branch 19 → 21 taken 674 times.
|
1209 | if (.not. allocated(this%constraint)) then |
| 452 | 535 | b_proj_vec = b_vec | |
| 453 | else | ||
| 454 | ! TODO if inhomoheneous then update b_vec with this%rhs_constraint_vec | ||
| 455 | 1348 | select case (this%constraint%mode) | |
| 456 | case (CONSTRAINT_EXTRACTION) | ||
| 457 |
1/2✗ Branch 23 → 24 not taken.
✓ Branch 23 → 25 taken 674 times.
|
674 | PetscCall(VecCreate(this%space_row%tp_spaces(1)%domain%comm, b_proj_vec, ierr)) |
| 458 |
1/2✗ Branch 26 → 27 not taken.
✓ Branch 26 → 28 taken 674 times.
|
674 | PetscCall(VecSetSizes(b_proj_vec, this%constraint%my_nr_cols, this%constraint%nr_cols, ierr)) |
| 459 |
1/2✗ Branch 29 → 30 not taken.
✓ Branch 29 → 31 taken 674 times.
|
674 | PetscCall(VecSetFromOptions(b_proj_vec, ierr)) |
| 460 | |||
| 461 |
1/2✗ Branch 32 → 33 not taken.
✓ Branch 32 → 34 taken 674 times.
|
674 | PetscCall(VecAssemblyBegin(b_proj_vec, ierr)) |
| 462 |
1/2✗ Branch 35 → 36 not taken.
✓ Branch 35 → 40 taken 674 times.
|
674 | PetscCall(VecAssemblyEnd(b_proj_vec, ierr)) |
| 463 | case (CONSTRAINT_PROJECTION) | ||
| 464 |
1/5✓ Branch 21 → 22 taken 674 times.
✗ Branch 21 → 37 not taken.
✗ Branch 21 → 40 not taken.
✗ Branch 38 → 39 not taken.
✗ Branch 38 → 40 not taken.
|
674 | PetscCall(VecDuplicate(b_vec, b_proj_vec, ierr)) |
| 465 | end select | ||
| 466 |
1/2✗ Branch 41 → 42 not taken.
✓ Branch 41 → 43 taken 674 times.
|
674 | PetscCall(MatMultTranspose(this%constraint%mat, b_vec, b_proj_vec, ierr)) |
| 467 | end if | ||
| 468 | |||
| 469 |
1/2✗ Branch 44 → 45 not taken.
✓ Branch 44 → 46 taken 1209 times.
|
1209 | PetscCall(VecDuplicate(b_proj_vec, x_proj_vec, ierr)) |
| 470 |
1/2✗ Branch 47 → 48 not taken.
✓ Branch 47 → 49 taken 1209 times.
|
1209 | PetscCall(KSPSolve(this%mat_ksp, b_proj_vec, x_proj_vec, ierr)) |
| 471 | |||
| 472 |
1/2✗ Branch 50 → 51 not taken.
✓ Branch 50 → 52 taken 1209 times.
|
1209 | PetscCall(KSPGetIterationNumber(this%mat_ksp, this%last_nr_iterations, ierr)) |
| 473 |
1/2✗ Branch 53 → 54 not taken.
✓ Branch 53 → 55 taken 1209 times.
|
1209 | PetscCall(KSPGetResidualNorm(this%mat_ksp, this%last_residual_norm, ierr)) |
| 474 |
1/2✗ Branch 56 → 57 not taken.
✓ Branch 56 → 58 taken 1209 times.
|
1209 | PetscCall(KSPGetConvergedReason(this%mat_ksp, this%last_converged_reason, ierr)) |
| 475 | |||
| 476 |
1/2✓ Branch 58 → 59 taken 1209 times.
✗ Branch 58 → 72 not taken.
|
1209 | if (verbosity_ > VERBOSITY_WARN_NEVER) then |
| 477 | 1209 | warn_msg = this%get_warning_message(failure) | |
| 478 |
1/2✗ Branch 60 → 61 not taken.
✓ Branch 60 → 67 taken 1209 times.
|
1209 | if (failure) then |
| 479 | ✗ | if (verbosity_ == VERBOSITY_ERROR_ON_FAILURE) then | |
| 480 | ✗ | error stop trim(warn_msg) | |
| 481 | else | ||
| 482 | ✗ | write (*, '(A)') trim(warn_msg) | |
| 483 | end if | ||
| 484 |
1/2✗ Branch 67 → 68 not taken.
✓ Branch 67 → 72 taken 1209 times.
|
1209 | else if (verbosity_ == VERBOSITY_WARN_ALWAYS) then |
| 485 | ✗ | write (*, '(A)') trim(warn_msg) | |
| 486 | end if | ||
| 487 | end if | ||
| 488 | |||
| 489 |
2/2✓ Branch 72 → 73 taken 674 times.
✓ Branch 72 → 83 taken 535 times.
|
1209 | if (allocated(this%constraint)) then |
| 490 |
1/2✓ Branch 73 → 74 taken 674 times.
✗ Branch 73 → 80 not taken.
|
674 | if (this%constraint%mode == CONSTRAINT_EXTRACTION) then |
| 491 |
1/2✗ Branch 75 → 76 not taken.
✓ Branch 75 → 77 taken 674 times.
|
674 | PetscCall(VecDuplicate(b_vec, x_vec, ierr)) |
| 492 |
1/2✗ Branch 78 → 79 not taken.
✓ Branch 78 → 84 taken 674 times.
|
674 | PetscCall(MatMult(this%constraint%mat, x_proj_vec, x_vec, ierr)) |
| 493 | else | ||
| 494 | ✗ | PetscCall(VecCopy(x_proj_vec, x_vec, ierr)) | |
| 495 | end if | ||
| 496 | else | ||
| 497 | 535 | x_vec = x_proj_vec | |
| 498 | end if | ||
| 499 | |||
| 500 |
2/2✓ Branch 84 → 85 taken 125 times.
✓ Branch 84 → 87 taken 1084 times.
|
1209 | if (x%is_initialized()) then |
| 501 | 125 | call x%reset(x_vec) | |
| 502 | else | ||
| 503 | 1084 | call x%init(this%space_col, x_vec) | |
| 504 | end if | ||
| 505 | |||
| 506 |
2/2✓ Branch 89 → 90 taken 674 times.
✓ Branch 89 → 96 taken 535 times.
|
1209 | if (allocated(this%constraint)) then |
| 507 |
1/2✗ Branch 91 → 92 not taken.
✓ Branch 91 → 93 taken 674 times.
|
674 | PetscCall(VecDestroy(b_proj_vec, ierr)) |
| 508 |
1/2✗ Branch 94 → 95 not taken.
✓ Branch 94 → 96 taken 674 times.
|
674 | PetscCall(VecDestroy(x_proj_vec, ierr)) |
| 509 | end if | ||
| 510 |
1/2✗ Branch 97 → 98 not taken.
✓ Branch 97 → 99 taken 1209 times.
|
1209 | PetscCall(VecDestroy(b_vec, ierr)) |
| 511 |
1/2✗ Branch 100 → 101 not taken.
✓ Branch 100 → 102 taken 1209 times.
|
1209 | PetscCall(VecDestroy(x_vec, ierr)) |
| 512 | |||
| 513 | ! PetscCall(KSPView(this%mat_ksp, PETSC_VIEWER_STDOUT_WORLD, ierr)) | ||
| 514 | |||
| 515 |
2/2✓ Branch 2 → 3 taken 134 times.
✓ Branch 2 → 9 taken 1075 times.
|
1209 | end subroutine solve_generic_solver |
| 516 | |||
| 517 | 1217 | character(512) function get_warning_message(this, failure) result(warn_msg) | |
| 518 | implicit none | ||
| 519 | class(AbstractSolver), intent(in) :: this | ||
| 520 | |||
| 521 | logical, intent(out) :: failure | ||
| 522 | character(len=30) :: conv_msg, solv_msg, reason_msg | ||
| 523 | character(len=100) :: reason_detail_msg | ||
| 524 | |||
| 525 | 1217 | failure = PETSC_CONVERGED_REASON_INT(this%last_converged_reason) <= 0 | |
| 526 | |||
| 527 |
2/2✓ Branch 2 → 3 taken 34 times.
✓ Branch 2 → 13 taken 1183 times.
|
1217 | if (this%options%pc_type == 'hypre') then |
| 528 | 34 | write (solv_msg, '(7A)') 'KSP(', trim(this%options%ksp_type), ', ', trim(this%options%pc_type), ', ', trim(this%options%hypre_type), ')' | |
| 529 | else | ||
| 530 | 1183 | write (solv_msg, '(5A)') 'KSP(', trim(this%options%ksp_type), ', ', trim(this%options%pc_type), ')' | |
| 531 | end if | ||
| 532 | |||
| 533 |
1/2✗ Branch 21 → 22 not taken.
✓ Branch 21 → 26 taken 1217 times.
|
1217 | if (failure) then |
| 534 | ✗ | write (conv_msg, '(A)') ' did not converge' | |
| 535 | else | ||
| 536 | 1217 | write (conv_msg, '(A)') ' converged' | |
| 537 | end if | ||
| 538 | |||
| 539 | 1217 | select case (PETSC_CONVERGED_REASON_INT(this%last_converged_reason)) | |
| 540 | case (1) | ||
| 541 | ✗ | reason_msg = 'KSP_CONVERGED_RTOL_NORMAL' | |
| 542 | ✗ | reason_detail_msg = 'requested decrease in the residual for the normal equations' | |
| 543 | case (9) | ||
| 544 | ✗ | reason_msg = 'KSP_CONVERGED_ATOL_NORMAL' | |
| 545 | ✗ | reason_detail_msg = 'requested absolute value in the residual for the normal equations' | |
| 546 | case (2) | ||
| 547 | 983 | reason_msg = 'KSP_CONVERGED_RTOL' | |
| 548 | 983 | reason_detail_msg = 'requested decrease in the residual' | |
| 549 | case (3) | ||
| 550 | 234 | reason_msg = 'KSP_CONVERGED_ATOL' | |
| 551 | 234 | reason_detail_msg = 'requested absolute value in the residual' | |
| 552 | case (4) | ||
| 553 | ✗ | reason_msg = 'KSP_CONVERGED_ITS' | |
| 554 | ✗ | reason_detail_msg = 'requested number of iterations' | |
| 555 | case (5) | ||
| 556 | ✗ | reason_msg = 'KSP_CONVERGED_NEG_CURVE' | |
| 557 | ✗ | reason_detail_msg = '' | |
| 558 | case (6) | ||
| 559 | ✗ | reason_msg = 'KSP_CONVERGED_STEP_LENGTH' | |
| 560 | ✗ | reason_detail_msg = '' | |
| 561 | case (7) | ||
| 562 | ✗ | reason_msg = 'KSP_CONVERGED_HAPPY_BREAKDOWN' | |
| 563 | ✗ | reason_detail_msg = 'happy breakdown (meaning early convergence of the `KSPType` occurred).' | |
| 564 | case (-2) | ||
| 565 | ✗ | reason_msg = 'KSP_DIVERGED_NULL' | |
| 566 | ✗ | reason_detail_msg = 'breakdown when solving the Hessenberg system within GMRES' | |
| 567 | case (-3) | ||
| 568 | ✗ | reason_msg = 'KSP_DIVERGED_ITS' | |
| 569 | ✗ | reason_detail_msg = 'requested number of iterations' | |
| 570 | case (-4) | ||
| 571 | ✗ | reason_msg = 'KSP_DIVERGED_DTOL' | |
| 572 | ✗ | reason_detail_msg = 'large increase in the residual norm' | |
| 573 | case (-5) | ||
| 574 | ✗ | reason_msg = 'KSP_DIVERGED_BREAKDOWN' | |
| 575 | ✗ | reason_detail_msg = 'breakdown in the Krylov method' | |
| 576 | case (-6) | ||
| 577 | ✗ | reason_msg = 'KSP_DIVERGED_BREAKDOWN_BICG' | |
| 578 | ✗ | reason_detail_msg = 'breakdown in the `KSPBGCS` Krylov method' | |
| 579 | case (-7) | ||
| 580 | ✗ | reason_msg = 'KSP_DIVERGED_NONSYMMETRIC' | |
| 581 | ✗ | reason_detail_msg = 'the operator or preonditioner was not symmetric for a `KSPType` that requires symmetry' | |
| 582 | case (-8) | ||
| 583 | ✗ | reason_msg = 'KSP_DIVERGED_INDEFINITE_PC' | |
| 584 | ✗ | reason_detail_msg = 'the preconditioner was indefinite for a `KSPType` that requires it be definite' | |
| 585 | case (-9) | ||
| 586 | ✗ | reason_msg = 'KSP_DIVERGED_NANORINF' | |
| 587 | ✗ | reason_detail_msg = 'a not a number of infinity was detected in a vector during the computation' | |
| 588 | case (-10) | ||
| 589 | ✗ | reason_msg = 'KSP_DIVERGED_INDEFINITE_MAT' | |
| 590 | ✗ | reason_detail_msg = 'the operator was indefinite for a `KSPType` that requires it be definite' | |
| 591 | case (-11) | ||
| 592 | ✗ | reason_msg = 'KSP_DIVERGED_PC_FAILED' | |
| 593 | ✗ | reason_detail_msg = 'the action of the preconditioner failed for some reason' | |
| 594 | case (0) | ||
| 595 | ✗ | reason_msg = 'KSP_CONVERGED_ITERATING' | |
| 596 |
2/20✗ Branch 30 → 31 not taken.
✗ Branch 30 → 32 not taken.
✓ Branch 30 → 33 taken 983 times.
✓ Branch 30 → 34 taken 234 times.
✗ Branch 30 → 35 not taken.
✗ Branch 30 → 36 not taken.
✗ Branch 30 → 37 not taken.
✗ Branch 30 → 38 not taken.
✗ Branch 30 → 39 not taken.
✗ Branch 30 → 40 not taken.
✗ Branch 30 → 41 not taken.
✗ Branch 30 → 42 not taken.
✗ Branch 30 → 43 not taken.
✗ Branch 30 → 44 not taken.
✗ Branch 30 → 45 not taken.
✗ Branch 30 → 46 not taken.
✗ Branch 30 → 47 not taken.
✗ Branch 30 → 48 not taken.
✗ Branch 30 → 49 not taken.
✗ Branch 30 → 50 not taken.
|
1217 | reason_detail_msg = 'the solver is still iterating' |
| 597 | end select | ||
| 598 | |||
| 599 | 1217 | write (warn_msg, '(4A,I0,A,E10.3,5A)') 'AbstractSolver::solve: ', trim(solv_msg), trim(conv_msg), ' after ', & | |
| 600 | 1217 | & this%last_nr_iterations, ' iterations with residual norm ', this%last_residual_norm, ' and convergence reason ', & | |
| 601 | 2434 | & trim(reason_msg), ' (', trim(reason_detail_msg), ')' | |
| 602 | 1217 | end function get_warning_message | |
| 603 | |||
| 604 | ! !> @brief Subtract the gradient-divergence matrix from the curl-curl matrix | ||
| 605 | ! !> | ||
| 606 | ! !> @param[out] mat_out The output matrix (curl-curl matrix minus the gradient-divergence matrix) | ||
| 607 | ! !> @param[in] mat_curlcurl The input curl-curl matrix | ||
| 608 | ! !> @param[in] space The 1-form space | ||
| 609 | ! subroutine subtract_graddiv_from_curlcul_matrix(mat_out, mat_curlcurl, space) | ||
| 610 | ! use m_mform_basis, only: MFormSpace, previous | ||
| 611 | ! use m_mform_matrix, only: MassMatrix, gradient_matrix | ||
| 612 | ! implicit none | ||
| 613 | |||
| 614 | ! Mat, intent(out) :: mat_out | ||
| 615 | ! Mat, intent(in) :: mat_curlcurl | ||
| 616 | ! type(MFormSpace), intent(in) :: space | ||
| 617 | |||
| 618 | ! Mat :: grad | ||
| 619 | ! type(MassMatrix) :: mass0 | ||
| 620 | ! type(MFormSpace) :: space0 | ||
| 621 | ! real(wp), parameter :: EXPECTED_FILL_RATIO = 1.0_wp | ||
| 622 | ! integer :: ierr | ||
| 623 | |||
| 624 | ! space0 = previous(space) | ||
| 625 | ! call gradient_matrix(grad, space0, space) | ||
| 626 | |||
| 627 | ! call mass0%init(space0) | ||
| 628 | ! PetscCall(MatRARt(mass0%mat, grad, MAT_INITIAL_MATRIX, EXPECTED_FILL_RATIO, mat_out, ierr)) | ||
| 629 | ! call mass0%destroy() | ||
| 630 | |||
| 631 | ! ! PetscCall(MatMatTransposeMult(grad, grad, MAT_INITIAL_MATRIX, EXPECTED_FILL_RATIO, mat_out, ierr)) | ||
| 632 | |||
| 633 | ! PetscCall(MatAXPY(mat_out, +1._wp, mat_curlcurl, DIFFERENT_NONZERO_PATTERN, ierr)) | ||
| 634 | |||
| 635 | ! PetscCall(MatDestroy(grad, ierr)) | ||
| 636 | |||
| 637 | ! end subroutine subtract_graddiv_from_curlcul_matrix | ||
| 638 | |||
| 639 | !> @brief Initialize the AMS (Auxilliary-space Maxwell Solver) preconditioner from HYPRE for the generic solver | ||
| 640 | !> | ||
| 641 | !> @param[inout] pc The preconditioner object to be initialized | ||
| 642 | !> @param[in] space The 1-form space associated with the solver | ||
| 643 | !> @param[in] constraint _(optional)_ An optional constraint to be applied to the solution | ||
| 644 | 22 | subroutine init_ams_pc(pc, space, constraint) | |
| 645 | use m_mform_matrix, only: DiffDiffMatrix, DiffMatrix, gradient_matrix | ||
| 646 | use m_mform_basis, only: previous, size | ||
| 647 | use m_tensorprod, only: size, TensorProdIndices, TensorProdIndexList | ||
| 648 | use m_domain_decomp, only: get_petsc | ||
| 649 | implicit none | ||
| 650 | |||
| 651 | PC, intent(inout) :: pc | ||
| 652 | type(MFormSpace), intent(in) :: space | ||
| 653 | type(MFormConstraint), intent(in), optional :: constraint | ||
| 654 | |||
| 655 | 22 | type(MFormSpace) :: space0 | |
| 656 | Vec :: constant_x, constant_y, constant_z | ||
| 657 | Vec :: constant_x_super, constant_y_super, constant_z_super | ||
| 658 | Vec :: zero_conductivity | ||
| 659 | Mat :: mat_tmp, grad | ||
| 660 |
17/26✓ Branch 3 → 4 taken 22 times.
✓ Branch 3 → 5 taken 22 times.
✓ Branch 6 → 7 taken 66 times.
✓ Branch 6 → 8 taken 22 times.
✓ Branch 126 → 127 taken 66 times.
✓ Branch 126 → 136 taken 22 times.
✓ Branch 127 → 128 taken 66 times.
✗ Branch 127 → 129 not taken.
✓ Branch 129 → 130 taken 66 times.
✗ Branch 129 → 131 not taken.
✓ Branch 131 → 132 taken 66 times.
✗ Branch 131 → 133 not taken.
✓ Branch 133 → 134 taken 66 times.
✗ Branch 133 → 135 not taken.
✓ Branch 137 → 138 taken 22 times.
✓ Branch 137 → 146 taken 22 times.
✓ Branch 138 → 139 taken 22 times.
✗ Branch 138 → 140 not taken.
✓ Branch 140 → 141 taken 22 times.
✗ Branch 140 → 142 not taken.
✓ Branch 142 → 143 taken 22 times.
✗ Branch 142 → 144 not taken.
✗ Branch 144 → 137 not taken.
✓ Branch 144 → 145 taken 22 times.
✓ Branch 146 → 147 taken 22 times.
✗ Branch 146 → 148 not taken.
|
264 | type(TensorProdIndexList) :: resp_inds(3), col_inds(1) |
| 661 | IS :: is_resp_inds_rows, is_resp_inds_cols | ||
| 662 | integer :: ierr, blk, nr_unknowns, my_nr_unknowns | ||
| 663 | real(wp), parameter :: EXPECTED_FILL_RATIO = 1.0_wp | ||
| 664 | |||
| 665 | 22 | if (.not. present(constraint)) then | |
| 666 | ✗ | do blk = 1, 3 | |
| 667 | call resp_inds(blk)%init(space%tp_spaces(blk)%rank_resp_bspline, space%tp_spaces(blk)%rank_resp_bspline, & | ||
| 668 | ✗ | & space%tp_spaces(blk)%rank_l0) | |
| 669 | end do | ||
| 670 | ✗ | nr_unknowns = size(space) | |
| 671 | ✗ | my_nr_unknowns = size(space, my_rank=.true.) | |
| 672 | else | ||
| 673 |
12/22✓ Branch 15 → 16 taken 66 times.
✓ Branch 15 → 39 taken 22 times.
✓ Branch 16 → 17 taken 66 times.
✗ Branch 16 → 29 not taken.
✓ Branch 17 → 18 taken 66 times.
✗ Branch 17 → 19 not taken.
✓ Branch 20 → 21 taken 66 times.
✗ Branch 20 → 22 not taken.
✓ Branch 23 → 24 taken 66 times.
✗ Branch 23 → 25 not taken.
✓ Branch 26 → 27 taken 66 times.
✗ Branch 26 → 28 not taken.
✓ Branch 29 → 30 taken 66 times.
✗ Branch 29 → 38 not taken.
✗ Branch 30 → 31 not taken.
✓ Branch 30 → 32 taken 66 times.
✗ Branch 32 → 33 not taken.
✓ Branch 32 → 34 taken 66 times.
✗ Branch 34 → 35 not taken.
✓ Branch 34 → 36 taken 66 times.
✗ Branch 36 → 37 not taken.
✓ Branch 36 → 38 taken 66 times.
|
88 | resp_inds = constraint%resp_inds |
| 674 | 22 | nr_unknowns = constraint%nr_cols | |
| 675 | 22 | my_nr_unknowns = constraint%my_nr_cols | |
| 676 | end if | ||
| 677 | |||
| 678 | ! It aims to solve the following curl-curl problem: | ||
| 679 | ! (\alpha curl v, curl u) + (\beta v, u) = (v, f) | ||
| 680 | ! where \alpha > 0, \beta >= 0 | ||
| 681 | ! in our case we have \beta = 0 | ||
| 682 | ! TODO more general solver with \beta > 0 | ||
| 683 | |||
| 684 | ! This solver requires | ||
| 685 | ! - the gradient matrix (whose columns span part of the null-space of the curl operator) | ||
| 686 | 22 | space0 = previous(space) | |
| 687 | 22 | call gradient_matrix(grad, space0, space) | |
| 688 |
1/2✓ Branch 42 → 43 taken 22 times.
✗ Branch 42 → 54 not taken.
|
22 | if (present(constraint)) then |
| 689 | ! PetscCall(MatTransposeMatMult(projector%mat, grad, MAT_INITIAL_MATRIX, EXPECTED_FILL_RATIO, mat_tmp, )ierr) | ||
| 690 | 22 | call get_petsc(is_resp_inds_rows, resp_inds, space%tp_spaces(1)%domain) | |
| 691 | |||
| 692 | call col_inds(1)%init(space0%tp_spaces(1)%rank_resp_bspline, space0%tp_spaces(1)%rank_resp_bspline, & | ||
| 693 | 22 | space0%tp_spaces(1)%rank_l0) | |
| 694 | 22 | call get_petsc(is_resp_inds_cols, col_inds, space0%tp_spaces(1)%domain) | |
| 695 | |||
| 696 |
1/2✗ Branch 47 → 48 not taken.
✓ Branch 47 → 50 taken 22 times.
|
22 | PetscCall(MatCreateSubMatrix(grad, is_resp_inds_rows, is_resp_inds_cols, MAT_INITIAL_MATRIX, mat_tmp, ierr)) |
| 697 |
1/2✗ Branch 51 → 52 not taken.
✓ Branch 51 → 53 taken 22 times.
|
22 | PetscCall(MatDestroy(grad, ierr)) |
| 698 | 22 | grad = mat_tmp | |
| 699 | |||
| 700 | end if | ||
| 701 | |||
| 702 |
1/2✗ Branch 55 → 56 not taken.
✓ Branch 55 → 57 taken 22 times.
|
22 | PetscCall(PCHYPRESetDiscreteGradient(pc, grad, ierr)) |
| 703 | |||
| 704 | 22 | call interior_vector(zero_conductivity, space, my_nr_unknowns, nr_unknowns) | |
| 705 |
1/2✗ Branch 59 → 60 not taken.
✓ Branch 59 → 61 taken 22 times.
|
22 | PetscCall(PCHYPREAMSSetInteriorNodes(pc, zero_conductivity, ierr)) |
| 706 | |||
| 707 | ! - the constant vectors in each direction (if periodic in the corresponding direction then such vectors are in the | ||
| 708 | ! null-space of the curl operator but not in the span of the columns of the gradient matrix) | ||
| 709 |
1/2✓ Branch 61 → 62 taken 22 times.
✗ Branch 61 → 74 not taken.
|
22 | if (present(constraint)) then |
| 710 | 22 | call constant_vector(constant_x_super, space, 1) | |
| 711 |
1/2✗ Branch 64 → 65 not taken.
✓ Branch 64 → 66 taken 22 times.
|
22 | PetscCall(VecGetSubVector(constant_x_super, is_resp_inds_rows, constant_x, ierr)) |
| 712 | |||
| 713 | 22 | call constant_vector(constant_y_super, space, 2) | |
| 714 |
1/2✗ Branch 68 → 69 not taken.
✓ Branch 68 → 70 taken 22 times.
|
22 | PetscCall(VecGetSubVector(constant_y_super, is_resp_inds_rows, constant_y, ierr)) |
| 715 | |||
| 716 | 22 | call constant_vector(constant_z_super, space, 3) | |
| 717 |
1/2✗ Branch 72 → 73 not taken.
✓ Branch 72 → 77 taken 22 times.
|
22 | PetscCall(VecGetSubVector(constant_z_super, is_resp_inds_rows, constant_z, ierr)) |
| 718 | else | ||
| 719 | ✗ | call constant_vector(constant_x, space, 1) | |
| 720 | ✗ | call constant_vector(constant_y, space, 2) | |
| 721 | ✗ | call constant_vector(constant_z, space, 3) | |
| 722 | end if | ||
| 723 |
1/2✗ Branch 78 → 79 not taken.
✓ Branch 78 → 80 taken 22 times.
|
22 | PetscCall(PCHYPRESetEdgeConstantVectors(pc, constant_x, constant_y, constant_z, ierr)) |
| 724 | |||
| 725 | ! Since \beta = 0 we provide a null matrix | ||
| 726 |
1/2✗ Branch 81 → 82 not taken.
✓ Branch 81 → 83 taken 22 times.
|
22 | PetscCall(PCHYPRESetBetaPoissonMatrix(pc, PETSC_NULL_MAT, ierr)) |
| 727 | |||
| 728 | ! TODO provide the \alpha matrix? | ||
| 729 | |||
| 730 | ! TODO consider the following from the Hypre documentation of the AMS solver: | ||
| 731 | ! In the high-order AMS interface, the user does not need to provide the coordinates of the vertices (or the representations of | ||
| 732 | ! the constant vector fields in the edge basis), but instead should construct and pass the Nedelec interpolation matrix which | ||
| 733 | ! maps (high-order) vector nodal finite elements into the (high-order) Nedelec space. | ||
| 734 | |||
| 735 |
1/2✓ Branch 83 → 84 taken 22 times.
✗ Branch 83 → 108 not taken.
|
22 | if (present(constraint)) then |
| 736 |
1/2✗ Branch 85 → 86 not taken.
✓ Branch 85 → 87 taken 22 times.
|
22 | PetscCall(VecRestoreSubVector(constant_x_super, is_resp_inds_rows, constant_x, ierr)) |
| 737 |
1/2✗ Branch 88 → 89 not taken.
✓ Branch 88 → 90 taken 22 times.
|
22 | PetscCall(VecRestoreSubVector(constant_y_super, is_resp_inds_rows, constant_y, ierr)) |
| 738 |
1/2✗ Branch 91 → 92 not taken.
✓ Branch 91 → 93 taken 22 times.
|
22 | PetscCall(VecRestoreSubVector(constant_z_super, is_resp_inds_rows, constant_z, ierr)) |
| 739 | |||
| 740 |
1/2✗ Branch 94 → 95 not taken.
✓ Branch 94 → 96 taken 22 times.
|
22 | PetscCall(VecDestroy(constant_x_super, ierr)) |
| 741 |
1/2✗ Branch 97 → 98 not taken.
✓ Branch 97 → 99 taken 22 times.
|
22 | PetscCall(VecDestroy(constant_y_super, ierr)) |
| 742 |
1/2✗ Branch 100 → 101 not taken.
✓ Branch 100 → 102 taken 22 times.
|
22 | PetscCall(VecDestroy(constant_z_super, ierr)) |
| 743 | |||
| 744 |
1/2✗ Branch 103 → 104 not taken.
✓ Branch 103 → 105 taken 22 times.
|
22 | PetscCall(ISDestroy(is_resp_inds_rows, ierr)) |
| 745 |
1/2✗ Branch 106 → 107 not taken.
✓ Branch 106 → 117 taken 22 times.
|
22 | PetscCall(ISDestroy(is_resp_inds_cols, ierr)) |
| 746 | else | ||
| 747 | ✗ | PetscCall(VecDestroy(constant_x, ierr)) | |
| 748 | ✗ | PetscCall(VecDestroy(constant_y, ierr)) | |
| 749 | ✗ | PetscCall(VecDestroy(constant_z, ierr)) | |
| 750 | end if | ||
| 751 |
1/2✗ Branch 118 → 119 not taken.
✓ Branch 118 → 120 taken 22 times.
|
22 | PetscCall(VecDestroy(zero_conductivity, ierr)) |
| 752 | |||
| 753 |
1/2✗ Branch 121 → 122 not taken.
✓ Branch 121 → 123 taken 22 times.
|
22 | PetscCall(MatDestroy(grad, ierr)) |
| 754 |
2/4✗ Branch 8 → 9 not taken.
✓ Branch 8 → 14 taken 22 times.
✓ Branch 123 → 124 taken 22 times.
✗ Branch 123 → 125 not taken.
|
66 | end subroutine init_ams_pc |
| 755 | |||
| 756 | !> @brief Get the constrained PETSc matrix for the given input matrix and constraint | ||
| 757 | !> | ||
| 758 | !> If projection is used, this yields the matrix | ||
| 759 | !> L <- P^T * L * P + coeff * (I - P)^T * M * (I - P) | ||
| 760 | !> and if extraction is used, this yields the matrix | ||
| 761 | !> L <- E^T * L * E | ||
| 762 | !> | ||
| 763 | !> @param[out] mat_out The output matrix (the constrained PETSc matrix) | ||
| 764 | !> @param[in] mat_in The input matrix (PETSc matrix) | ||
| 765 | !> @param[in] constraint_sol The constraint object acting on the solution space | ||
| 766 | !> @param[in] coeff _(optional)_ The coefficient for the (I - P)^T * M * (I - P) part (default is 0.0) | ||
| 767 | !> @param[in] constraint_res _(optional)_ The constraint object acting on the residual space (default is the same as constraint_sol) | ||
| 768 | 410 | subroutine get_constrained_petsc_matrix(mat_out, mat_in, constraint_sol, coeff, constraint_res) | |
| 769 | use m_common, only: wp | ||
| 770 | use m_mform_matrix, only: AbstractMatrix | ||
| 771 | use m_mform_constraint_global, only: MFormConstraint, CONSTRAINT_EXTRACTION, CONSTRAINT_PROJECTION | ||
| 772 | |||
| 773 | implicit none | ||
| 774 | |||
| 775 | Mat, intent(out) :: mat_out | ||
| 776 | Mat, intent(in) :: mat_in | ||
| 777 | type(MFormConstraint), intent(in) :: constraint_sol | ||
| 778 | real(wp), optional, intent(in) :: coeff | ||
| 779 | type(MFormConstraint), optional, intent(in) :: constraint_res | ||
| 780 | |||
| 781 | real(wp), parameter :: EXPECTED_FILL_RATIO = 1.0_wp | ||
| 782 | real(wp) :: coeff_ | ||
| 783 | integer :: ierr | ||
| 784 | Mat :: mat_tmp | ||
| 785 | |||
| 786 | if (present(coeff)) then | ||
| 787 | coeff_ = coeff | ||
| 788 | else | ||
| 789 | coeff_ = 0.0_wp | ||
| 790 | end if | ||
| 791 | |||
| 792 | 410 | if (.not. present(constraint_res)) then | |
| 793 | ! We assume constraint_res = constraint_sol | ||
| 794 | |||
| 795 | ! Compute | ||
| 796 | ! L <- P^T * L * P, | ||
| 797 | ! for L = mat_in, P = constraint_sol%mat | ||
| 798 |
1/2✗ Branch 4 → 5 not taken.
✓ Branch 4 → 15 taken 378 times.
|
378 | PetscCall(MatPtAP(mat_in, constraint_sol%mat, MAT_INITIAL_MATRIX, EXPECTED_FILL_RATIO, mat_out, ierr)) |
| 799 | |||
| 800 | ! if (constraint_sol%mode == CONSTRAINT_PROJECTION .and. coeff_ /= 0.0_wp) then | ||
| 801 | ! ! We add the (I - P)^T * M * (I - P) part, where M = constraint_sol%mass_mat | ||
| 802 | ! PetscCall(MatPtAP(constraint_sol%mass_mat, constraint_sol%mat_perp, MAT_INITIAL_MATRIX, EXPECTED_FILL_RATIO, mat_tmp, ierr)) | ||
| 803 | ! PetscCall(MatAXPY(mat_out, coeff_, mat_tmp, DIFFERENT_NONZERO_PATTERN, ierr)) | ||
| 804 | ! ! We now have | ||
| 805 | ! ! L <- P^T * L * P + \gamma * (I - P)^T * M * (I - P) | ||
| 806 | ! ! where \gamma = coeff_ is implemented by the derived AbstractMatrix type | ||
| 807 | ! PetscCall(MatDestroy(mat_tmp, ierr)) | ||
| 808 | ! end if | ||
| 809 | else | ||
| 810 | ! Compute | ||
| 811 | ! L <- P1^T * L * P0, | ||
| 812 | ! for L = mat_in, P0 = constraint_sol%mat, P1 = constraint_res%mat | ||
| 813 |
1/2✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 32 times.
|
32 | PetscCall(MatMatMult(mat_in, constraint_sol%mat, MAT_INITIAL_MATRIX, EXPECTED_FILL_RATIO, mat_tmp, ierr)) |
| 814 |
1/2✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 32 times.
|
32 | PetscCall(MatTransposeMatMult(constraint_res%mat, mat_tmp, MAT_INITIAL_MATRIX, EXPECTED_FILL_RATIO, mat_out, ierr)) |
| 815 |
1/2✗ Branch 13 → 14 not taken.
✓ Branch 13 → 15 taken 32 times.
|
32 | PetscCall(MatDestroy(mat_tmp, ierr)) |
| 816 | |||
| 817 | if (constraint_sol%mode == CONSTRAINT_PROJECTION .and. coeff_ /= 0.0_wp) then | ||
| 818 | ! TODO | ||
| 819 | end if | ||
| 820 | end if | ||
| 821 |
2/2✓ Branch 2 → 3 taken 378 times.
✓ Branch 2 → 6 taken 32 times.
|
410 | end subroutine get_constrained_petsc_matrix |
| 822 | |||
| 823 | !> @brief Create a vector with all entries set to 1.0 in the interior of the m-form space | ||
| 824 | !> | ||
| 825 | !> @param[out] vec The output PETSc vector | ||
| 826 | !> @param[in] space The m-form space for which the vector is created | ||
| 827 | !> @param[in] my_nr_unknowns The number of unknowns owned by the current rank in the m-form space | ||
| 828 | !> @param[in] nr_unknowns The number of unknowns in the m-form space | ||
| 829 | 22 | subroutine interior_vector(vec, space, my_nr_unknowns, nr_unknowns) | |
| 830 | use m_tensorprod, only: TensorProdIndexList | ||
| 831 | use m_common, only: wp | ||
| 832 | use m_mform_basis, only: MFormSpace, size | ||
| 833 | |||
| 834 | implicit none | ||
| 835 | |||
| 836 | Vec, intent(out) :: vec | ||
| 837 | type(MFormSpace), intent(in) :: space | ||
| 838 | integer, intent(in) :: my_nr_unknowns, nr_unknowns | ||
| 839 | |||
| 840 | integer :: ierr | ||
| 841 | |||
| 842 |
1/2✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 22 times.
|
22 | PetscCall(VecCreate(space%tp_spaces(1)%domain%comm, vec, ierr)) |
| 843 | |||
| 844 |
1/2✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 22 times.
|
22 | PetscCall(VecSetSizes(vec, my_nr_unknowns, nr_unknowns, ierr)) |
| 845 |
1/2✗ Branch 9 → 10 not taken.
✓ Branch 9 → 11 taken 22 times.
|
22 | PetscCall(VecSetFromOptions(vec, ierr)) |
| 846 | |||
| 847 |
1/2✗ Branch 12 → 13 not taken.
✓ Branch 12 → 14 taken 22 times.
|
22 | PetscCall(VecSet(vec, 1._wp, ierr)) |
| 848 | |||
| 849 |
1/2✗ Branch 15 → 16 not taken.
✓ Branch 15 → 17 taken 22 times.
|
22 | PetscCall(VecAssemblyBegin(vec, ierr)) |
| 850 |
1/2✗ Branch 18 → 19 not taken.
✓ Branch 18 → 20 taken 22 times.
|
22 | PetscCall(VecAssemblyEnd(vec, ierr)) |
| 851 | |||
| 852 | end subroutine interior_vector | ||
| 853 | |||
| 854 | !> @brief Create a constant vector with all entries in the given direction set to 1.0 in the interior of the m-form space | ||
| 855 | !> | ||
| 856 | !> @param[out] vec The output PETSc vector | ||
| 857 | !> @param[in] space The m-form space for which the vector is created | ||
| 858 | !> @param[in] dir The direction in which the constant vector is created (1 for 'x', 2 for 'y', or 3 for 'z') | ||
| 859 | 66 | subroutine constant_vector(vec, space, dir) | |
| 860 | use m_tensorprod, only: size, TensorProdIndexList | ||
| 861 | use m_common, only: wp | ||
| 862 | use m_mform_basis, only: MFormSpace, size, MFormFun, get_petsc | ||
| 863 | |||
| 864 | implicit none | ||
| 865 | |||
| 866 | Vec, intent(out) :: vec | ||
| 867 | type(MFormSpace), intent(in) :: space | ||
| 868 | integer, intent(in) :: dir | ||
| 869 | |||
| 870 | 66 | type(MFormFun) :: const_fun | |
| 871 | |||
| 872 | 66 | call const_fun%init(space) | |
| 873 | ! Only the leader writes to the shared memory window to avoid race conditions | ||
| 874 |
1/2✓ Branch 3 → 4 taken 66 times.
✗ Branch 3 → 13 not taken.
|
66 | if (const_fun%tp_funs(dir)%shmem_window%leader()) then |
| 875 |
6/6✓ Branch 5 → 6 taken 1094 times.
✓ Branch 5 → 13 taken 66 times.
✓ Branch 7 → 8 taken 10852 times.
✓ Branch 7 → 12 taken 1094 times.
✓ Branch 9 → 10 taken 218036 times.
✓ Branch 9 → 11 taken 10852 times.
|
230048 | const_fun%tp_funs(dir)%data(:, :, :) = 1.0_wp |
| 876 | end if | ||
| 877 | ! Synchronize so all ranks see the updated data before get_petsc reads it | ||
| 878 | 66 | call const_fun%tp_funs(dir)%shmem_window%fence() | |
| 879 | 66 | call get_petsc(vec, const_fun) | |
| 880 | 66 | call const_fun%destroy() | |
| 881 |
2/4✓ Branch 16 → 17 taken 66 times.
✗ Branch 16 → 18 not taken.
✗ Branch 18 → 19 not taken.
✓ Branch 18 → 20 taken 66 times.
|
66 | end subroutine constant_vector |
| 882 | |||
| 883 | !> @brief Destroy the generic solver object | ||
| 884 | !> | ||
| 885 | !> @param[inout] this The generic solver object to be destroyed | ||
| 886 | 984 | subroutine destroy_abstract_solver(this) | |
| 887 | use m_common, only: wp | ||
| 888 | implicit none | ||
| 889 | |||
| 890 | class(AbstractSolver), intent(inout) :: this | ||
| 891 | |||
| 892 | integer :: ierr | ||
| 893 | |||
| 894 |
1/2✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 984 times.
|
984 | PetscCall(MatDestroy(this%mat, ierr)) |
| 895 |
1/2✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 984 times.
|
984 | PetscCall(KSPDestroy(this%mat_ksp, ierr)) |
| 896 | |||
| 897 | end subroutine destroy_abstract_solver | ||
| 898 |
2/120None:
✓ Branch 47 → 48 taken 259 times.
✓ Branch 47 → 49 taken 359 times.
__m_mform_solver_MOD___copy_m_mform_solver_Genericsolver:
✗ Branch 2 → 3 not taken.
✗ Branch 2 → 52 not taken.
✗ Branch 3 → 4 not taken.
✗ Branch 3 → 5 not taken.
✗ Branch 6 → 7 not taken.
✗ Branch 6 → 8 not taken.
✗ Branch 9 → 10 not taken.
✗ Branch 9 → 51 not taken.
✗ Branch 10 → 11 not taken.
✗ Branch 10 → 12 not taken.
✗ Branch 14 → 15 not taken.
✗ Branch 14 → 19 not taken.
✗ Branch 15 → 16 not taken.
✗ Branch 15 → 17 not taken.
✗ Branch 19 → 20 not taken.
✗ Branch 19 → 21 not taken.
✗ Branch 22 → 23 not taken.
✗ Branch 22 → 24 not taken.
✗ Branch 25 → 26 not taken.
✗ Branch 25 → 27 not taken.
✗ Branch 28 → 29 not taken.
✗ Branch 28 → 30 not taken.
✗ Branch 31 → 32 not taken.
✗ Branch 31 → 33 not taken.
✗ Branch 34 → 35 not taken.
✗ Branch 34 → 50 not taken.
✗ Branch 36 → 37 not taken.
✗ Branch 36 → 52 not taken.
✗ Branch 37 → 38 not taken.
✗ Branch 37 → 39 not taken.
✗ Branch 40 → 41 not taken.
✗ Branch 40 → 42 not taken.
✗ Branch 43 → 44 not taken.
✗ Branch 43 → 45 not taken.
✗ Branch 46 → 47 not taken.
✗ Branch 46 → 48 not taken.
__m_mform_solver_MOD___final_m_mform_solver_Genericsolver:
✗ Branch 3 → 4 not taken.
✗ Branch 3 → 7 not taken.
✗ Branch 4 → 5 not taken.
✗ Branch 4 → 6 not taken.
✗ Branch 8 → 9 not taken.
✗ Branch 8 → 91 not taken.
✗ Branch 10 → 11 not taken.
✗ Branch 10 → 12 not taken.
✗ Branch 12 → 13 not taken.
✗ Branch 12 → 14 not taken.
✗ Branch 14 → 15 not taken.
✗ Branch 14 → 16 not taken.
✗ Branch 16 → 17 not taken.
✗ Branch 16 → 90 not taken.
✗ Branch 17 → 18 not taken.
✗ Branch 17 → 89 not taken.
✗ Branch 18 → 19 not taken.
✗ Branch 18 → 20 not taken.
✗ Branch 21 → 22 not taken.
✗ Branch 21 → 65 not taken.
✗ Branch 22 → 23 not taken.
✗ Branch 22 → 25 not taken.
✗ Branch 23 → 24 not taken.
✗ Branch 23 → 25 not taken.
✗ Branch 25 → 26 not taken.
✗ Branch 25 → 64 not taken.
✗ Branch 26 → 27 not taken.
✗ Branch 26 → 28 not taken.
✗ Branch 28 → 29 not taken.
✗ Branch 28 → 41 not taken.
✗ Branch 30 → 31 not taken.
✗ Branch 30 → 40 not taken.
✗ Branch 31 → 32 not taken.
✗ Branch 31 → 33 not taken.
✗ Branch 33 → 34 not taken.
✗ Branch 33 → 35 not taken.
✗ Branch 35 → 36 not taken.
✗ Branch 35 → 37 not taken.
✗ Branch 37 → 38 not taken.
✗ Branch 37 → 39 not taken.
✗ Branch 41 → 42 not taken.
✗ Branch 41 → 52 not taken.
✗ Branch 43 → 44 not taken.
✗ Branch 43 → 51 not taken.
✗ Branch 44 → 45 not taken.
✗ Branch 44 → 46 not taken.
✗ Branch 46 → 47 not taken.
✗ Branch 46 → 48 not taken.
✗ Branch 48 → 49 not taken.
✗ Branch 48 → 50 not taken.
✗ Branch 52 → 53 not taken.
✗ Branch 52 → 63 not taken.
✗ Branch 54 → 55 not taken.
✗ Branch 54 → 62 not taken.
✗ Branch 55 → 56 not taken.
✗ Branch 55 → 57 not taken.
✗ Branch 57 → 58 not taken.
✗ Branch 57 → 59 not taken.
✗ Branch 59 → 60 not taken.
✗ Branch 59 → 61 not taken.
✗ Branch 65 → 66 not taken.
✗ Branch 65 → 67 not taken.
✗ Branch 67 → 68 not taken.
✗ Branch 67 → 69 not taken.
✗ Branch 69 → 70 not taken.
✗ Branch 69 → 71 not taken.
✗ Branch 71 → 72 not taken.
✗ Branch 71 → 73 not taken.
✗ Branch 73 → 74 not taken.
✗ Branch 73 → 75 not taken.
✗ Branch 75 → 76 not taken.
✗ Branch 75 → 88 not taken.
✗ Branch 77 → 78 not taken.
✗ Branch 77 → 87 not taken.
✗ Branch 78 → 79 not taken.
✗ Branch 78 → 80 not taken.
✗ Branch 80 → 81 not taken.
✗ Branch 80 → 82 not taken.
✗ Branch 82 → 83 not taken.
✗ Branch 82 → 84 not taken.
✗ Branch 84 → 85 not taken.
✗ Branch 84 → 86 not taken.
|
618 | end module m_mform_solver |
| 899 |