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