mform/m_mform_coupled_solver.f90
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | !> @brief Module containing the linear solver for coupled block systems of m-form equations | ||
| 2 | module m_mform_coupled_solver | ||
| 3 | use m_common, only: wp, VERBOSITY_WARN_ON_FAILURE | ||
| 4 | use petscis | ||
| 5 | use m_mform_basis, only: MFormSpace | ||
| 6 | use m_mform_constraint_global, only: MFormConstraint | ||
| 7 | use m_mform_matrix, only: AbstractMatrix | ||
| 8 | use m_mform_solver, only: AbstractSolver, get_constrained_petsc_matrix | ||
| 9 | use m_mform_solver, only: SolverInitOptions, SolverSolveOptions, SolverInitTraits | ||
| 10 | #include "petsc.fi" | ||
| 11 | |||
| 12 | implicit none | ||
| 13 | private | ||
| 14 | |||
| 15 | public :: MatrixHolder, CoupledSolver | ||
| 16 | |||
| 17 | !> @brief Wrapper type that allows passing heterogeneous polymorphic matrix blocks. | ||
| 18 | !> | ||
| 19 | !> The wrapped matrix is a non-owning copy of the object state and should not call | ||
| 20 | !> matrix-specific destroy; the original matrix object owns the PETSc resources. | ||
| 21 | type MatrixHolder | ||
| 22 | class(AbstractMatrix), allocatable :: elt | ||
| 23 | contains | ||
| 24 | procedure :: init => init_matrix_holder | ||
| 25 | procedure :: destroy => destroy_matrix_holder | ||
| 26 | end type | ||
| 27 | |||
| 28 | !> @brief Solver for coupled block systems assembled into a PETSc nested matrix | ||
| 29 | type, extends(AbstractSolver) :: CoupledSolver | ||
| 30 | integer :: nr_blocks = 0 | ||
| 31 | type(MFormSpace), allocatable :: spaces(:) | ||
| 32 | type(MFormConstraint), allocatable :: constraints(:) | ||
| 33 | ! Stored row-major as [(1,1), (1,2), ..., (n,n)] for MatCreateNest. | ||
| 34 | Mat, allocatable :: block_mats(:) | ||
| 35 | contains | ||
| 36 | procedure :: init => init_coupled_solver | ||
| 37 | procedure :: solve => solve_coupled_solver | ||
| 38 | procedure :: destroy => destroy_coupled_solver | ||
| 39 | end type | ||
| 40 | |||
| 41 | contains | ||
| 42 | |||
| 43 | !> @brief Initialize a matrix holder from any AbstractMatrix-derived object | ||
| 44 | 64 | subroutine init_matrix_holder(this, mat) | |
| 45 | class(MatrixHolder), intent(inout) :: this | ||
| 46 | class(AbstractMatrix), intent(in) :: mat | ||
| 47 | |||
| 48 | 32 | call this%destroy() | |
| 49 |
2/4✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 32 times.
✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 32 times.
|
32 | allocate (this%elt, source=mat) |
| 50 | 32 | end subroutine init_matrix_holder | |
| 51 | |||
| 52 | !> @brief Destroy a matrix holder | ||
| 53 | 64 | subroutine destroy_matrix_holder(this) | |
| 54 | class(MatrixHolder), intent(inout) :: this | ||
| 55 | |||
| 56 |
3/4✓ Branch 2 → 3 taken 32 times.
✓ Branch 2 → 7 taken 32 times.
✓ Branch 3 → 4 taken 32 times.
✗ Branch 3 → 6 not taken.
|
96 | if (allocated(this%elt)) deallocate (this%elt) |
| 57 | 64 | end subroutine destroy_matrix_holder | |
| 58 | |||
| 59 | !> @brief Initialize the coupled solver for a square block system | ||
| 60 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 8 times.
|
8 | subroutine init_coupled_solver(this, amat_blocks, constraint1, constraint2, constraint3, coord_transform, options) |
| 61 | use m_mform_basis, only: size | ||
| 62 | use m_mform_constraint_abstract, only: MFormConstraintLocal | ||
| 63 | use m_coord_transform_abstract, only: CoordTransformAbstract | ||
| 64 | use m_mform_solver, only: create_ksp | ||
| 65 | |||
| 66 | class(CoupledSolver), intent(inout) :: this | ||
| 67 | type(MatrixHolder), intent(in) :: amat_blocks(:, :) | ||
| 68 | class(MFormConstraintLocal), optional, intent(in) :: constraint1, constraint2, constraint3 | ||
| 69 | class(CoordTransformAbstract), optional, intent(in) :: coord_transform | ||
| 70 | type(SolverInitOptions), optional, intent(in) :: options | ||
| 71 | |||
| 72 | PC :: pc | ||
| 73 | integer :: i, j, n_blocks, ierr, idx, default_max_iter | ||
| 74 | logical :: has_constraints, coupled_is_symmetric, coupled_is_spd | ||
| 75 | type(SolverInitTraits) :: traits | ||
| 76 | |||
| 77 | 8 | call this%destroy() | |
| 78 | |||
| 79 |
1/2✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 8 times.
|
8 | if (size(amat_blocks, 1) /= size(amat_blocks, 2)) then |
| 80 | ✗ | error stop "CoupledSolver::init_coupled_solver: coupled block matrix must be square." | |
| 81 | end if | ||
| 82 | |||
| 83 | 8 | n_blocks = size(amat_blocks, 1) | |
| 84 |
1/2✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 8 times.
|
8 | if (n_blocks <= 0) then |
| 85 | ✗ | error stop "CoupledSolver::init_coupled_solver: coupled block matrix must be non-empty." | |
| 86 | end if | ||
| 87 | |||
| 88 |
2/2✓ Branch 9 → 10 taken 16 times.
✓ Branch 9 → 16 taken 8 times.
|
24 | do i = 1, n_blocks |
| 89 |
2/2✓ Branch 10 → 11 taken 32 times.
✓ Branch 10 → 14 taken 16 times.
|
56 | do j = 1, n_blocks |
| 90 |
1/2✗ Branch 11 → 12 not taken.
✓ Branch 11 → 13 taken 32 times.
|
48 | if (.not. allocated(amat_blocks(i, j)%elt)) then |
| 91 | ✗ | error stop "CoupledSolver::init_coupled_solver: all matrix holders must be initialized." | |
| 92 | end if | ||
| 93 | end do | ||
| 94 | end do | ||
| 95 | |||
| 96 | 8 | this%nr_blocks = n_blocks | |
| 97 |
8/14✗ Branch 17 → 18 not taken.
✓ Branch 17 → 19 taken 8 times.
✓ Branch 19 → 18 taken 8 times.
✗ Branch 19 → 20 not taken.
✓ Branch 20 → 21 taken 8 times.
✗ Branch 20 → 22 not taken.
✗ Branch 22 → 23 not taken.
✓ Branch 22 → 24 taken 8 times.
✗ Branch 24 → 25 not taken.
✓ Branch 24 → 26 taken 8 times.
✗ Branch 26 → 27 not taken.
✓ Branch 26 → 28 taken 8 times.
✓ Branch 29 → 30 taken 16 times.
✓ Branch 29 → 31 taken 8 times.
|
40 | allocate (this%spaces(n_blocks)) |
| 98 | |||
| 99 |
2/2✓ Branch 32 → 33 taken 16 times.
✓ Branch 32 → 41 taken 8 times.
|
24 | do i = 1, n_blocks |
| 100 |
4/8✓ Branch 33 → 34 taken 16 times.
✗ Branch 33 → 37 not taken.
✓ Branch 34 → 35 taken 16 times.
✗ Branch 34 → 36 not taken.
✓ Branch 37 → 38 taken 16 times.
✗ Branch 37 → 40 not taken.
✗ Branch 38 → 39 not taken.
✓ Branch 38 → 40 taken 16 times.
|
24 | this%spaces(i) = amat_blocks(i, i)%elt%space_col |
| 101 | end do | ||
| 102 | |||
| 103 |
2/2✓ Branch 42 → 43 taken 16 times.
✓ Branch 42 → 53 taken 8 times.
|
24 | do i = 1, n_blocks |
| 104 |
2/2✓ Branch 43 → 44 taken 32 times.
✓ Branch 43 → 51 taken 16 times.
|
56 | do j = 1, n_blocks |
| 105 |
1/2✗ Branch 44 → 45 not taken.
✓ Branch 44 → 46 taken 32 times.
|
32 | if (amat_blocks(i, j)%elt%space_row%m /= this%spaces(i)%m) then |
| 106 | ✗ | error stop "CoupledSolver::init_coupled_solver: inconsistent row m-form space across a block row." | |
| 107 | end if | ||
| 108 |
1/2✗ Branch 46 → 47 not taken.
✓ Branch 46 → 48 taken 32 times.
|
32 | if (amat_blocks(i, j)%elt%space_col%m /= this%spaces(j)%m) then |
| 109 | ✗ | error stop "CoupledSolver::init_coupled_solver: inconsistent column m-form space across a block column." | |
| 110 | end if | ||
| 111 | |||
| 112 |
1/2✗ Branch 48 → 49 not taken.
✓ Branch 48 → 50 taken 32 times.
|
48 | if (amat_blocks(i, j)%elt%space_col%m /= amat_blocks(j, i)%elt%space_row%m) then |
| 113 | ✗ | error stop "CoupledSolver::init_coupled_solver: coupled block dimensions are not mutually compatible." | |
| 114 | end if | ||
| 115 | end do | ||
| 116 | end do | ||
| 117 | |||
| 118 |
2/12✓ Branch 54 → 55 taken 8 times.
✗ Branch 54 → 56 not taken.
✗ Branch 55 → 56 not taken.
✓ Branch 55 → 60 taken 8 times.
✗ Branch 56 → 57 not taken.
✗ Branch 56 → 58 not taken.
✗ Branch 57 → 58 not taken.
✗ Branch 57 → 60 not taken.
✗ Branch 58 → 59 not taken.
✗ Branch 58 → 104 not taken.
✗ Branch 59 → 60 not taken.
✗ Branch 59 → 104 not taken.
|
8 | has_constraints = present(constraint1) .or. present(constraint2) .or. present(constraint3) |
| 119 | if (has_constraints) then | ||
| 120 |
14/22✗ Branch 60 → 61 not taken.
✓ Branch 60 → 62 taken 8 times.
✓ Branch 62 → 61 taken 8 times.
✗ Branch 62 → 63 not taken.
✓ Branch 63 → 64 taken 8 times.
✗ Branch 63 → 65 not taken.
✗ Branch 65 → 66 not taken.
✓ Branch 65 → 67 taken 8 times.
✗ Branch 67 → 68 not taken.
✓ Branch 67 → 69 taken 8 times.
✗ Branch 69 → 70 not taken.
✓ Branch 69 → 71 taken 8 times.
✓ Branch 72 → 73 taken 40 times.
✓ Branch 72 → 74 taken 8 times.
✓ Branch 75 → 76 taken 16 times.
✓ Branch 75 → 86 taken 8 times.
✗ Branch 76 → 77 not taken.
✓ Branch 76 → 78 taken 16 times.
✓ Branch 80 → 81 taken 80 times.
✓ Branch 80 → 85 taken 16 times.
✗ Branch 81 → 82 not taken.
✓ Branch 81 → 83 taken 80 times.
|
160 | allocate (this%constraints(n_blocks)) |
| 121 |
2/2✓ Branch 87 → 88 taken 16 times.
✓ Branch 87 → 103 taken 8 times.
|
24 | do j = 1, n_blocks |
| 122 | 16 | call this%constraints(j)%init(this%spaces(j)) | |
| 123 |
2/4✓ Branch 89 → 90 taken 16 times.
✗ Branch 89 → 93 not taken.
✓ Branch 90 → 91 taken 16 times.
✗ Branch 90 → 93 not taken.
|
16 | if (present(constraint1)) call this%constraints(j)%insert(constraint1, coord_transform=coord_transform) |
| 124 |
2/4✓ Branch 93 → 94 taken 16 times.
✗ Branch 93 → 97 not taken.
✓ Branch 94 → 95 taken 16 times.
✗ Branch 94 → 97 not taken.
|
16 | if (present(constraint2)) call this%constraints(j)%insert(constraint2, coord_transform=coord_transform) |
| 125 |
1/4✗ Branch 97 → 98 not taken.
✓ Branch 97 → 101 taken 16 times.
✗ Branch 98 → 99 not taken.
✗ Branch 98 → 101 not taken.
|
16 | if (present(constraint3)) call this%constraints(j)%insert(constraint3, coord_transform=coord_transform) |
| 126 | 24 | call this%constraints(j)%assemble() | |
| 127 | end do | ||
| 128 | end if | ||
| 129 | |||
| 130 | 8 | this%options = SolverInitOptions() | |
| 131 |
1/2✓ Branch 105 → 106 taken 8 times.
✗ Branch 105 → 107 not taken.
|
8 | if (present(options)) this%options = options |
| 132 | coupled_is_symmetric = .true. | ||
| 133 | coupled_is_spd = .false. | ||
| 134 |
2/2✓ Branch 108 → 109 taken 16 times.
✓ Branch 108 → 116 taken 8 times.
|
24 | do i = 1, n_blocks |
| 135 |
2/2✓ Branch 109 → 110 taken 32 times.
✓ Branch 109 → 114 taken 16 times.
|
56 | do j = 1, n_blocks |
| 136 |
4/4✓ Branch 110 → 111 taken 16 times.
✓ Branch 110 → 112 taken 16 times.
✓ Branch 111 → 112 taken 8 times.
✓ Branch 111 → 113 taken 8 times.
|
48 | coupled_is_symmetric = coupled_is_symmetric .and. amat_blocks(i, j)%elt%is_symmetric |
| 137 | end do | ||
| 138 | end do | ||
| 139 | |||
| 140 | 8 | traits%is_symmetric = coupled_is_symmetric | |
| 141 | traits%is_spd = coupled_is_spd | ||
| 142 | 8 | traits%is_parallel = this%spaces(1)%tp_spaces(1)%domain%nr_ranks > 1 | |
| 143 | traits%is_diffdiff = .false. | ||
| 144 | traits%m_form = -1 | ||
| 145 | 8 | traits%allow_ams = .false. | |
| 146 | |||
| 147 | 8 | call this%options%init_from_traits(traits) | |
| 148 | |||
| 149 | 8 | default_max_iter = this%options%max_iter | |
| 150 |
1/2✓ Branch 118 → 119 taken 8 times.
✗ Branch 118 → 120 not taken.
|
8 | if (default_max_iter < 0) default_max_iter = 0 |
| 151 |
1/2✓ Branch 120 → 121 taken 8 times.
✗ Branch 120 → 128 not taken.
|
8 | if (default_max_iter == 0) then |
| 152 | 8 | default_max_iter = 0 | |
| 153 |
2/2✓ Branch 122 → 123 taken 16 times.
✓ Branch 122 → 127 taken 8 times.
|
24 | do i = 1, n_blocks |
| 154 |
1/2✓ Branch 123 → 124 taken 16 times.
✗ Branch 123 → 125 not taken.
|
24 | if (allocated(this%constraints)) then |
| 155 | 16 | default_max_iter = default_max_iter + this%constraints(i)%nr_cols | |
| 156 | else | ||
| 157 | ✗ | default_max_iter = default_max_iter + size(this%spaces(i)) | |
| 158 | end if | ||
| 159 | end do | ||
| 160 | end if | ||
| 161 | |||
| 162 | 8 | call this%options%apply_default_solve_options(default_max_iter, VERBOSITY_WARN_ON_FAILURE) | |
| 163 | |||
| 164 |
8/14✗ Branch 129 → 130 not taken.
✓ Branch 129 → 131 taken 8 times.
✓ Branch 131 → 130 taken 8 times.
✗ Branch 131 → 132 not taken.
✓ Branch 132 → 133 taken 8 times.
✗ Branch 132 → 134 not taken.
✗ Branch 134 → 135 not taken.
✓ Branch 134 → 136 taken 8 times.
✗ Branch 136 → 137 not taken.
✓ Branch 136 → 138 taken 8 times.
✗ Branch 138 → 139 not taken.
✓ Branch 138 → 140 taken 8 times.
✓ Branch 141 → 142 taken 32 times.
✓ Branch 141 → 143 taken 8 times.
|
56 | allocate (this%block_mats(n_blocks*n_blocks)) |
| 165 | idx = 0 | ||
| 166 |
2/2✓ Branch 143 → 144 taken 16 times.
✓ Branch 143 → 155 taken 8 times.
|
24 | do i = 1, n_blocks |
| 167 |
2/2✓ Branch 144 → 145 taken 32 times.
✓ Branch 144 → 153 taken 16 times.
|
56 | do j = 1, n_blocks |
| 168 | 32 | idx = idx + 1 | |
| 169 |
1/2✓ Branch 145 → 146 taken 32 times.
✗ Branch 145 → 149 not taken.
|
48 | if (allocated(this%constraints)) then |
| 170 | call get_constrained_petsc_matrix(this%block_mats(idx), amat_blocks(i, j)%elt%mat, this%constraints(j), & | ||
| 171 | 32 | coeff=amat_blocks(i, j)%elt%constraint_coeff(), constraint_res=this%constraints(i)) | |
| 172 | else | ||
| 173 | ✗ | PetscCall(MatConvert(amat_blocks(i, j)%elt%mat, "same", MAT_INITIAL_MATRIX, this%block_mats(idx), ierr)) | |
| 174 | end if | ||
| 175 | end do | ||
| 176 | end do | ||
| 177 | |||
| 178 |
1/2✗ Branch 157 → 158 not taken.
✓ Branch 157 → 159 taken 8 times.
|
8 | PetscCall(MatCreateNest(this%spaces(1)%tp_spaces(1)%domain%comm, n_blocks, PETSC_NULL_IS_ARRAY, n_blocks, \ |
| 179 | PETSC_NULL_IS_ARRAY, this%block_mats, this%mat, ierr)) | ||
| 180 |
1/2✗ Branch 160 → 161 not taken.
✓ Branch 160 → 162 taken 8 times.
|
8 | PetscCall(MatAssemblyBegin(this%mat, MAT_FINAL_ASSEMBLY, ierr)) |
| 181 |
1/2✗ Branch 163 → 164 not taken.
✓ Branch 163 → 165 taken 8 times.
|
8 | PetscCall(MatAssemblyEnd(this%mat, MAT_FINAL_ASSEMBLY, ierr)) |
| 182 | |||
| 183 | 8 | call create_ksp(this%mat_ksp, this%mat, this%spaces(1), this%options) | |
| 184 | |||
| 185 | 8 | end subroutine init_coupled_solver | |
| 186 | |||
| 187 | !> @brief Solve the coupled system using nested PETSc vectors | ||
| 188 |
2/4✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 8 times.
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 8 times.
|
8 | subroutine solve_coupled_solver(this, x, b, options) |
| 189 | use m_common, only: VERBOSITY_WARN_NEVER, VERBOSITY_WARN_ON_FAILURE, VERBOSITY_WARN_ALWAYS, VERBOSITY_ERROR_ON_FAILURE | ||
| 190 | use m_mform_basis, only: MFormFun, get_petsc | ||
| 191 | use m_mform_constraint_global, only: CONSTRAINT_EXTRACTION, CONSTRAINT_PROJECTION | ||
| 192 | |||
| 193 | class(CoupledSolver), intent(inout) :: this | ||
| 194 | type(MFormFun), intent(inout) :: x(:) | ||
| 195 | type(MFormFun), intent(in) :: b(:) | ||
| 196 | type(SolverSolveOptions), optional, intent(in) :: options | ||
| 197 | |||
| 198 | integer :: i, n_blocks, ierr, verbosity_ | ||
| 199 | real(wp) :: rtol_, atol_ | ||
| 200 | integer :: max_iter_ | ||
| 201 | character(len=512) :: warn_msg | ||
| 202 | logical :: failure | ||
| 203 | |||
| 204 | Vec, allocatable :: b_vec(:), b_proj_vec(:), x_proj_vec(:), x_vec(:) | ||
| 205 | Vec :: b_nest, x_nest | ||
| 206 | |||
| 207 | 8 | n_blocks = this%nr_blocks | |
| 208 | |||
| 209 |
1/2✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 8 times.
|
8 | if (size(x) /= n_blocks .or. size(b) /= n_blocks) then |
| 210 | ✗ | error stop "CoupledSolver::solve_coupled_solver: x and b must match number of block unknowns." | |
| 211 | end if | ||
| 212 | |||
| 213 | 8 | rtol_ = this%options%rtol | |
| 214 | 8 | atol_ = this%options%atol | |
| 215 | 8 | max_iter_ = this%options%max_iter | |
| 216 |
1/2✗ Branch 9 → 10 not taken.
✓ Branch 9 → 16 taken 8 times.
|
8 | if (present(options)) then |
| 217 | ✗ | if (options%rtol >= 0._wp) rtol_ = options%rtol | |
| 218 | ✗ | if (options%atol >= 0._wp) atol_ = options%atol | |
| 219 | ✗ | if (options%max_iter >= 0) max_iter_ = options%max_iter | |
| 220 | end if | ||
| 221 | |||
| 222 |
1/2✓ Branch 16 → 17 taken 8 times.
✗ Branch 16 → 20 not taken.
|
8 | if (this%spaces(1)%tp_spaces(1)%domain%my_rank > 0) then |
| 223 | verbosity_ = VERBOSITY_WARN_NEVER | ||
| 224 | else | ||
| 225 | 8 | verbosity_ = this%options%verbosity | |
| 226 |
1/2✓ Branch 17 → 18 taken 8 times.
✗ Branch 17 → 19 not taken.
|
8 | if (present(options)) then |
| 227 | ✗ | if (options%verbosity >= 0) verbosity_ = options%verbosity | |
| 228 | end if | ||
| 229 | end if | ||
| 230 | |||
| 231 |
1/2✗ Branch 21 → 22 not taken.
✓ Branch 21 → 23 taken 8 times.
|
8 | PetscCall(KSPSetTolerances(this%mat_ksp, rtol_, atol_, 1.e10_wp, max_iter_, ierr)) |
| 232 | |||
| 233 |
26/44✗ Branch 23 → 24 not taken.
✓ Branch 23 → 25 taken 8 times.
✓ Branch 25 → 24 taken 8 times.
✗ Branch 25 → 26 not taken.
✓ Branch 26 → 27 taken 8 times.
✗ Branch 26 → 28 not taken.
✓ Branch 28 → 29 taken 8 times.
✗ Branch 28 → 30 not taken.
✗ Branch 30 → 31 not taken.
✓ Branch 30 → 32 taken 8 times.
✗ Branch 32 → 33 not taken.
✓ Branch 32 → 34 taken 8 times.
✓ Branch 35 → 36 taken 16 times.
✓ Branch 35 → 37 taken 8 times.
✗ Branch 37 → 38 not taken.
✓ Branch 37 → 39 taken 8 times.
✓ Branch 39 → 38 taken 8 times.
✗ Branch 39 → 40 not taken.
✗ Branch 40 → 41 not taken.
✓ Branch 40 → 42 taken 8 times.
✗ Branch 42 → 43 not taken.
✓ Branch 42 → 44 taken 8 times.
✓ Branch 44 → 45 taken 16 times.
✓ Branch 44 → 46 taken 8 times.
✗ Branch 46 → 47 not taken.
✓ Branch 46 → 48 taken 8 times.
✓ Branch 48 → 47 taken 8 times.
✗ Branch 48 → 49 not taken.
✗ Branch 49 → 50 not taken.
✓ Branch 49 → 51 taken 8 times.
✗ Branch 51 → 52 not taken.
✓ Branch 51 → 53 taken 8 times.
✓ Branch 53 → 54 taken 16 times.
✓ Branch 53 → 55 taken 8 times.
✗ Branch 55 → 56 not taken.
✓ Branch 55 → 57 taken 8 times.
✓ Branch 57 → 56 taken 8 times.
✗ Branch 57 → 58 not taken.
✗ Branch 58 → 59 not taken.
✓ Branch 58 → 60 taken 8 times.
✗ Branch 60 → 61 not taken.
✓ Branch 60 → 62 taken 8 times.
✓ Branch 62 → 63 taken 16 times.
✓ Branch 62 → 64 taken 8 times.
|
112 | allocate (b_vec(n_blocks), b_proj_vec(n_blocks), x_proj_vec(n_blocks), x_vec(n_blocks)) |
| 234 | |||
| 235 |
2/2✓ Branch 64 → 65 taken 16 times.
✓ Branch 64 → 101 taken 8 times.
|
24 | do i = 1, n_blocks |
| 236 |
1/2✗ Branch 65 → 66 not taken.
✓ Branch 65 → 67 taken 16 times.
|
16 | if (this%spaces(i)%m /= b(i)%space%m) then |
| 237 | ✗ | error stop "CoupledSolver::solve_coupled_solver: m-form spaces of rhs blocks do not match row spaces." | |
| 238 | end if | ||
| 239 | |||
| 240 | 16 | call get_petsc(b_vec(i), b(i)) | |
| 241 | |||
| 242 |
1/2✗ Branch 68 → 69 not taken.
✓ Branch 68 → 75 taken 16 times.
|
16 | if (.not. allocated(this%constraints)) then |
| 243 | ✗ | PetscCall(VecDuplicate(b_vec(i), b_proj_vec(i), ierr)) | |
| 244 | ✗ | PetscCall(VecCopy(b_vec(i), b_proj_vec(i), ierr)) | |
| 245 | else | ||
| 246 | 32 | select case (this%constraints(i)%mode) | |
| 247 | case (CONSTRAINT_EXTRACTION) | ||
| 248 |
1/2✗ Branch 77 → 78 not taken.
✓ Branch 77 → 79 taken 16 times.
|
16 | PetscCall(VecCreate(this%spaces(1)%tp_spaces(1)%domain%comm, b_proj_vec(i), ierr)) |
| 249 |
1/2✗ Branch 80 → 81 not taken.
✓ Branch 80 → 82 taken 16 times.
|
16 | PetscCall(VecSetSizes(b_proj_vec(i), this%constraints(i)%my_nr_cols, this%constraints(i)%nr_cols, ierr)) |
| 250 |
1/2✗ Branch 83 → 84 not taken.
✓ Branch 83 → 85 taken 16 times.
|
16 | PetscCall(VecSetFromOptions(b_proj_vec(i), ierr)) |
| 251 |
1/2✗ Branch 86 → 87 not taken.
✓ Branch 86 → 88 taken 16 times.
|
16 | PetscCall(VecAssemblyBegin(b_proj_vec(i), ierr)) |
| 252 |
1/2✗ Branch 89 → 90 not taken.
✓ Branch 89 → 94 taken 16 times.
|
16 | PetscCall(VecAssemblyEnd(b_proj_vec(i), ierr)) |
| 253 | case (CONSTRAINT_PROJECTION) | ||
| 254 |
1/5✓ Branch 75 → 76 taken 16 times.
✗ Branch 75 → 91 not taken.
✗ Branch 75 → 94 not taken.
✗ Branch 92 → 93 not taken.
✗ Branch 92 → 94 not taken.
|
16 | PetscCall(VecDuplicate(b_vec(i), b_proj_vec(i), ierr)) |
| 255 | end select | ||
| 256 |
1/2✗ Branch 95 → 96 not taken.
✓ Branch 95 → 97 taken 16 times.
|
16 | PetscCall(MatMultTranspose(this%constraints(i)%mat, b_vec(i), b_proj_vec(i), ierr)) |
| 257 | end if | ||
| 258 | |||
| 259 |
1/2✗ Branch 98 → 99 not taken.
✓ Branch 98 → 100 taken 16 times.
|
24 | PetscCall(VecDuplicate(b_proj_vec(i), x_proj_vec(i), ierr)) |
| 260 | end do | ||
| 261 | |||
| 262 |
1/2✗ Branch 103 → 104 not taken.
✓ Branch 103 → 105 taken 8 times.
|
8 | PetscCall(VecCreateNest(this%spaces(1)%tp_spaces(1)%domain%comm, n_blocks, PETSC_NULL_IS_ARRAY, b_proj_vec, b_nest, ierr)) |
| 263 |
1/2✗ Branch 106 → 107 not taken.
✓ Branch 106 → 108 taken 8 times.
|
8 | PetscCall(VecCreateNest(this%spaces(1)%tp_spaces(1)%domain%comm, n_blocks, PETSC_NULL_IS_ARRAY, x_proj_vec, x_nest, ierr)) |
| 264 | |||
| 265 |
1/2✗ Branch 109 → 110 not taken.
✓ Branch 109 → 111 taken 8 times.
|
8 | PetscCall(KSPSolve(this%mat_ksp, b_nest, x_nest, ierr)) |
| 266 | |||
| 267 |
1/2✗ Branch 112 → 113 not taken.
✓ Branch 112 → 114 taken 8 times.
|
8 | PetscCall(KSPGetIterationNumber(this%mat_ksp, this%last_nr_iterations, ierr)) |
| 268 |
1/2✗ Branch 115 → 116 not taken.
✓ Branch 115 → 117 taken 8 times.
|
8 | PetscCall(KSPGetResidualNorm(this%mat_ksp, this%last_residual_norm, ierr)) |
| 269 |
1/2✗ Branch 118 → 119 not taken.
✓ Branch 118 → 120 taken 8 times.
|
8 | PetscCall(KSPGetConvergedReason(this%mat_ksp, this%last_converged_reason, ierr)) |
| 270 | |||
| 271 |
1/2✓ Branch 120 → 121 taken 8 times.
✗ Branch 120 → 134 not taken.
|
8 | if (verbosity_ > VERBOSITY_WARN_NEVER) then |
| 272 | 8 | warn_msg = this%get_warning_message(failure) | |
| 273 |
1/2✗ Branch 122 → 123 not taken.
✓ Branch 122 → 129 taken 8 times.
|
8 | if (failure) then |
| 274 | ✗ | if (verbosity_ == VERBOSITY_ERROR_ON_FAILURE) then | |
| 275 | ✗ | error stop trim(warn_msg) | |
| 276 | else | ||
| 277 | ✗ | write (*, '(A)') trim(warn_msg) | |
| 278 | end if | ||
| 279 |
1/2✗ Branch 129 → 130 not taken.
✓ Branch 129 → 134 taken 8 times.
|
8 | else if (verbosity_ == VERBOSITY_WARN_ALWAYS) then |
| 280 | ✗ | write (*, '(A)') trim(warn_msg) | |
| 281 | end if | ||
| 282 | end if | ||
| 283 | |||
| 284 |
2/2✓ Branch 135 → 136 taken 16 times.
✓ Branch 135 → 156 taken 8 times.
|
24 | do i = 1, n_blocks |
| 285 |
1/2✗ Branch 137 → 138 not taken.
✓ Branch 137 → 139 taken 16 times.
|
16 | PetscCall(VecDuplicate(b_vec(i), x_vec(i), ierr)) |
| 286 | |||
| 287 |
1/2✓ Branch 139 → 140 taken 16 times.
✗ Branch 139 → 147 not taken.
|
16 | if (allocated(this%constraints)) then |
| 288 |
1/2✓ Branch 140 → 141 taken 16 times.
✗ Branch 140 → 144 not taken.
|
16 | if (this%constraints(i)%mode == CONSTRAINT_EXTRACTION) then |
| 289 |
1/2✗ Branch 142 → 143 not taken.
✓ Branch 142 → 150 taken 16 times.
|
16 | PetscCall(MatMult(this%constraints(i)%mat, x_proj_vec(i), x_vec(i), ierr)) |
| 290 | else | ||
| 291 | ✗ | PetscCall(VecCopy(x_proj_vec(i), x_vec(i), ierr)) | |
| 292 | end if | ||
| 293 | else | ||
| 294 | ✗ | PetscCall(VecCopy(x_proj_vec(i), x_vec(i), ierr)) | |
| 295 | end if | ||
| 296 | |||
| 297 |
1/2✗ Branch 150 → 151 not taken.
✓ Branch 150 → 153 taken 16 times.
|
24 | if (x(i)%is_initialized()) then |
| 298 | ✗ | call x(i)%reset(x_vec(i)) | |
| 299 | else | ||
| 300 | 16 | call x(i)%init(this%spaces(i), x_vec(i)) | |
| 301 | end if | ||
| 302 | end do | ||
| 303 | |||
| 304 |
1/2✗ Branch 158 → 159 not taken.
✓ Branch 158 → 160 taken 8 times.
|
8 | PetscCall(VecDestroy(b_nest, ierr)) |
| 305 |
1/2✗ Branch 161 → 162 not taken.
✓ Branch 161 → 163 taken 8 times.
|
8 | PetscCall(VecDestroy(x_nest, ierr)) |
| 306 | |||
| 307 |
2/2✓ Branch 163 → 164 taken 16 times.
✓ Branch 163 → 177 taken 8 times.
|
24 | do i = 1, n_blocks |
| 308 |
1/2✗ Branch 165 → 166 not taken.
✓ Branch 165 → 167 taken 16 times.
|
16 | PetscCall(VecDestroy(b_vec(i), ierr)) |
| 309 |
1/2✗ Branch 168 → 169 not taken.
✓ Branch 168 → 170 taken 16 times.
|
16 | PetscCall(VecDestroy(b_proj_vec(i), ierr)) |
| 310 |
1/2✗ Branch 171 → 172 not taken.
✓ Branch 171 → 173 taken 16 times.
|
16 | PetscCall(VecDestroy(x_proj_vec(i), ierr)) |
| 311 |
1/2✗ Branch 174 → 175 not taken.
✓ Branch 174 → 176 taken 16 times.
|
24 | PetscCall(VecDestroy(x_vec(i), ierr)) |
| 312 | end do | ||
| 313 |
5/18✓ Branch 6 → 7 taken 8 times.
✗ Branch 6 → 8 not taken.
✓ Branch 178 → 179 taken 8 times.
✗ Branch 178 → 180 not taken.
✓ Branch 180 → 181 taken 8 times.
✗ Branch 180 → 182 not taken.
✓ Branch 182 → 183 taken 8 times.
✗ Branch 182 → 184 not taken.
✓ Branch 184 → 185 taken 8 times.
✗ Branch 184 → 202 not taken.
✗ Branch 187 → 188 not taken.
✗ Branch 187 → 189 not taken.
✗ Branch 191 → 192 not taken.
✗ Branch 191 → 193 not taken.
✗ Branch 195 → 196 not taken.
✗ Branch 195 → 197 not taken.
✗ Branch 199 → 200 not taken.
✗ Branch 199 → 203 not taken.
|
16 | end subroutine solve_coupled_solver |
| 314 | |||
| 315 | !> @brief Destroy the coupled solver | ||
| 316 | 16 | subroutine destroy_coupled_solver(this) | |
| 317 | class(CoupledSolver), intent(inout) :: this | ||
| 318 | |||
| 319 | integer :: i, ierr | ||
| 320 | |||
| 321 |
1/2✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 16 times.
|
16 | PetscCall(MatDestroy(this%mat, ierr)) |
| 322 |
1/2✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 16 times.
|
16 | PetscCall(KSPDestroy(this%mat_ksp, ierr)) |
| 323 | |||
| 324 |
2/2✓ Branch 8 → 9 taken 8 times.
✓ Branch 8 → 19 taken 8 times.
|
16 | if (allocated(this%block_mats)) then |
| 325 |
2/2✓ Branch 10 → 11 taken 32 times.
✓ Branch 10 → 15 taken 8 times.
|
40 | do i = 1, size(this%block_mats) |
| 326 |
1/2✗ Branch 12 → 13 not taken.
✓ Branch 12 → 14 taken 32 times.
|
40 | PetscCall(MatDestroy(this%block_mats(i), ierr)) |
| 327 | end do | ||
| 328 |
1/2✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 8 times.
|
8 | deallocate (this%block_mats) |
| 329 | end if | ||
| 330 | |||
| 331 |
2/2✓ Branch 19 → 20 taken 8 times.
✓ Branch 19 → 102 taken 8 times.
|
16 | if (allocated(this%constraints)) then |
| 332 |
2/2✓ Branch 21 → 22 taken 16 times.
✓ Branch 21 → 24 taken 8 times.
|
24 | do i = 1, size(this%constraints) |
| 333 | 24 | call this%constraints(i)%destroy() | |
| 334 | end do | ||
| 335 |
15/72✓ Branch 25 → 26 taken 8 times.
✗ Branch 25 → 99 not taken.
✓ Branch 27 → 28 taken 16 times.
✓ Branch 27 → 99 taken 8 times.
✓ Branch 28 → 29 taken 16 times.
✗ Branch 28 → 30 not taken.
✓ Branch 31 → 32 taken 80 times.
✓ Branch 31 → 75 taken 16 times.
✗ Branch 32 → 33 not taken.
✓ Branch 32 → 35 taken 80 times.
✗ Branch 33 → 34 not taken.
✗ Branch 33 → 35 not taken.
✗ Branch 35 → 36 not taken.
✓ Branch 35 → 74 taken 80 times.
✗ Branch 36 → 37 not taken.
✗ Branch 36 → 38 not taken.
✗ Branch 38 → 39 not taken.
✗ Branch 38 → 51 not taken.
✗ Branch 40 → 41 not taken.
✗ Branch 40 → 50 not taken.
✗ Branch 41 → 42 not taken.
✗ Branch 41 → 43 not taken.
✗ Branch 43 → 44 not taken.
✗ Branch 43 → 45 not taken.
✗ Branch 45 → 46 not taken.
✗ Branch 45 → 47 not taken.
✗ Branch 47 → 48 not taken.
✗ Branch 47 → 49 not taken.
✗ Branch 51 → 52 not taken.
✗ Branch 51 → 62 not taken.
✗ Branch 53 → 54 not taken.
✗ Branch 53 → 61 not taken.
✗ Branch 54 → 55 not taken.
✗ Branch 54 → 56 not taken.
✗ Branch 56 → 57 not taken.
✗ Branch 56 → 58 not taken.
✗ Branch 58 → 59 not taken.
✗ Branch 58 → 60 not taken.
✗ Branch 62 → 63 not taken.
✗ Branch 62 → 73 not taken.
✗ Branch 64 → 65 not taken.
✗ Branch 64 → 72 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 75 → 76 not taken.
✓ Branch 75 → 77 taken 16 times.
✗ Branch 77 → 78 not taken.
✓ Branch 77 → 79 taken 16 times.
✗ Branch 79 → 80 not taken.
✓ Branch 79 → 81 taken 16 times.
✗ Branch 81 → 82 not taken.
✓ Branch 81 → 83 taken 16 times.
✗ Branch 83 → 84 not taken.
✓ Branch 83 → 85 taken 16 times.
✗ Branch 85 → 86 not taken.
✓ Branch 85 → 98 taken 16 times.
✗ Branch 87 → 88 not taken.
✗ Branch 87 → 97 not taken.
✗ Branch 88 → 89 not taken.
✗ Branch 88 → 90 not taken.
✗ Branch 90 → 91 not taken.
✗ Branch 90 → 92 not taken.
✗ Branch 92 → 93 not taken.
✗ Branch 92 → 94 not taken.
✗ Branch 94 → 95 not taken.
✗ Branch 94 → 96 not taken.
✗ Branch 99 → 100 not taken.
✓ Branch 99 → 101 taken 8 times.
|
104 | deallocate (this%constraints) |
| 336 | end if | ||
| 337 | |||
| 338 |
6/8✓ Branch 102 → 103 taken 8 times.
✓ Branch 102 → 111 taken 8 times.
✓ Branch 104 → 105 taken 16 times.
✓ Branch 104 → 108 taken 8 times.
✓ Branch 105 → 106 taken 16 times.
✗ Branch 105 → 107 not taken.
✗ Branch 108 → 109 not taken.
✓ Branch 108 → 110 taken 8 times.
|
32 | if (allocated(this%spaces)) deallocate (this%spaces) |
| 339 | 16 | this%nr_blocks = 0 | |
| 340 | end subroutine destroy_coupled_solver | ||
| 341 | |||
| 342 | ✗ | end module m_mform_coupled_solver | |
| 343 |