GCC Code Coverage Report


Directory: src/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 59.5% 184 / 0 / 309
Functions: 56.2% 9 / 0 / 16
Branches: 33.1% 256 / 0 / 774

mform/m_mform_constraint_global.f90
Line Branch Exec Source
1 !> @brief Module containing an object representing a linear constraint that can be imposed on an m-form function (MFormFun)
2 !>
3 !> In essence, this object simply wraps multiple local constraints (MFormConstraintLocal) and provides a way to assemble them into a
4 !> global constraint operator (either via extraction or projection). It is assumed that the local constraints do not overlap, i.e.,
5 !> that they act on different parts of the m-form coefficients.
6 module m_mform_constraint_global
7 use m_mform_basis, only: MFormSpace
8 use m_mform_constraint_abstract, only: MFormConstraintLocal
9 use m_tensorprod_indices, only: TensorProdIndexList
10 #include "petsc.fi"
11 implicit none
12
13 private
14 public :: MFormConstraint, CONSTRAINT_PROJECTION, CONSTRAINT_EXTRACTION
15
16 !> Impose the constraint via a (square) projection operator
17 integer, parameter :: CONSTRAINT_PROJECTION = 1
18 !> Impose the constraint via an extraction operator
19 integer, parameter :: CONSTRAINT_EXTRACTION = 2
20 !> The maximum number of constraints that can be imposed on an m-form
21 integer, parameter :: MAX_NR_CONSTRAINTS = 5
22
23 !> An object that wraps a local constraint (MFormConstraintLocal). This is needed because we can't have an array of abstract
24 !> types, but we can have an array of pointers to abstract types.
25 type ConstraintLocalHolder
26 !> The local constraint
27 class(MFormConstraintLocal), allocatable :: elt
28 contains
29 procedure :: init => init_constraint_local_holder
30 procedure :: destroy => destroy_constraint_local_holder
31 end type
32
33 !> @brief The object that represents linear constraints imposed on the coefficients of an m-form
34 type MFormConstraint
35 !> The m-form space to which the constraint applies
36 type(MFormSpace) :: space
37
38 !> The mode of the constraint, either CONSTRAINT_PROJECTION or CONSTRAINT_EXTRACTION
39 integer :: mode
40 !> Whether the constraint is symmetric (i.e., the projection matrix is symmetric)
41 logical :: is_symmetric
42 !> Whether the constraint has been assembled
43 logical :: is_assembled
44
45 !> The local constraints that are imposed on the m-form space
46 type(ConstraintLocalHolder) :: constraints(MAX_NR_CONSTRAINTS)
47 !> The number of constraints that have been inserted (less than or equal to MAX_NR_CONSTRAINTS)
48 integer :: nr_constraints
49 !> The number of columns in the projection/extraction matrix
50 integer :: nr_cols
51 !> The number of columns in the projection/extraction matrix on this rank
52 integer :: my_nr_cols
53 !> Whether the unknowns corresponding to each constraint are present on this rank
54 logical :: unknowns_are_on_this_rank(0:MAX_NR_CONSTRAINTS)
55
56 !> The total number of nonzeros per row in the projection/extraction matrix P
57 integer, allocatable :: total_local_nz_per_row(:)
58 !> The total number of nonzeros per row in the projection/extraction matrix P (off rank part)
59 integer, allocatable :: total_nonlocal_nz_per_row(:)
60 !> The total number of nonzeros per row in the perpendicular projection matrix I - P
61 integer, allocatable :: total_local_nz_per_row_perp(:)
62 !> The type of the column, 0 means no constraint, > 0 means corresponding to constraints(col_type)
63 integer, allocatable :: col_type(:)
64 !> The mapping from the original matrix indices to the new matrix indices (this is non-trivial if mode==CONSTRAINT_EXTRACTION)
65 integer, allocatable :: col_mapping(:)
66
67 !> The responsible tensor product indices of the coefficients that are kept in the projection/extraction matrix (this is
68 !> non-trivial if mode==CONSTRAINT_EXTRACTION)
69 type(TensorProdIndexList), allocatable :: resp_inds(:)
70
71 !> The projection matrix (if mode == CONSTRAINT_PROJECTION) or the extraction matrix (if mode == CONSTRAINT_EXTRACTION)
72 Mat :: mat
73 !> The perpendicular projection matrix (I - P) (if mode == CONSTRAINT_PROJECTION)
74 Mat :: mat_perp
75
76 !> The mass matrix that is used to compute the constrained operator
77 ! Mat :: mass_mat
78
79 contains
80 procedure :: init => init_mform_constraint
81 procedure :: insert => insert_in_mform_constraint
82 procedure, private :: insert_constraint_in_projection_matrix
83 procedure, private :: insert_constraint_in_extraction_matrix
84 procedure :: assemble
85 ! procedure :: apply => apply_mform_constraint
86 procedure :: destroy => destroy_mform_constraint
87 end type
88
89 contains
90 !> @brief Initialize the MFormConstraint object
91 !>
92 !> @param[inout] this The MFormConstraint object to initialize
93 !> @param[in] space The MFormSpace to which the constraint applies
94 ! !> @param[in] mass _(optional)_ The mass matrix that is used to compute the constrained operator
95 !> @param[in] mode _(optional)_ The mode in which the constraint is applied, either CONSTRAINT_PROJECTION or CONSTRAINT_EXTRACTION
96 512 subroutine init_mform_constraint(this, space, mode)
97 use m_common, only: wp
98 use m_mform_basis, only: size
99 use m_mform_matrix, only: MassMatrix
100 implicit none
101
102 class(MFormConstraint), intent(inout) :: this
103 type(MFormSpace), intent(in) :: space
104 ! type(MassMatrix), optional, intent(in) :: mass
105 character(*), optional, intent(in) :: mode
106
107 ! type(MassMatrix) :: mass_
108 integer :: ierr
109
110 512 call this%destroy()
111
112
1/2
✗ Branch 3 → 4 not taken.
✓ Branch 3 → 7 taken 512 times.
512 if (present(mode)) then
113 select case (mode)
114 ! case ('projection')
115 ! this%mode = CONSTRAINT_PROJECTION
116 case ('extraction')
117 this%mode = CONSTRAINT_EXTRACTION
118 case default
119 error stop "AbstractMatrix::init_mform_constraint: invalid mode"
120 end select
121 else
122 512 this%mode = CONSTRAINT_EXTRACTION
123 end if
124
125
5/8
✓ Branch 8 → 9 taken 512 times.
✗ Branch 8 → 12 not taken.
✓ Branch 9 → 10 taken 512 times.
✗ Branch 9 → 11 not taken.
✓ Branch 12 → 13 taken 512 times.
✗ Branch 12 → 15 not taken.
✓ Branch 13 → 14 taken 130 times.
✓ Branch 13 → 15 taken 382 times.
512 this%space = space
126 512 this%is_assembled = .false.
127
128 ! if (this%mode == CONSTRAINT_PROJECTION) then
129 ! if (present(mass)) then
130 ! ! NOTE a copy is required because the destructor of mass will destroy the matrix
131 ! PetscCall(MatConvert(mass%mat, "same", MAT_INITIAL_MATRIX, this%mass_mat, ierr))
132 ! else
133 ! call mass_%init(space)
134 ! this%mass_mat = mass_%mat
135 ! ! NOTE: we don't destroy the mass_ object because we made a shallow copy of the matrix
136 ! end if
137 ! end if
138
139 512 this%is_symmetric = .true.
140 512 this%nr_constraints = 0
141
142
7/14
✗ Branch 15 → 16 not taken.
✓ Branch 15 → 17 taken 512 times.
✓ Branch 17 → 16 taken 512 times.
✗ Branch 17 → 18 not taken.
✓ Branch 18 → 19 taken 512 times.
✗ Branch 18 → 20 not taken.
✓ Branch 20 → 21 taken 512 times.
✗ Branch 20 → 22 not taken.
✗ Branch 22 → 23 not taken.
✓ Branch 22 → 24 taken 512 times.
✗ Branch 24 → 25 not taken.
✓ Branch 24 → 26 taken 512 times.
✗ Branch 26 → 27 not taken.
✓ Branch 26 → 28 taken 512 times.
1536 allocate (this%total_local_nz_per_row(space%rank_l0:space%rank_l1))
143
5/10
✗ Branch 28 → 29 not taken.
✓ Branch 28 → 30 taken 512 times.
✓ Branch 30 → 29 taken 512 times.
✗ Branch 30 → 31 not taken.
✗ Branch 31 → 32 not taken.
✓ Branch 31 → 33 taken 512 times.
✗ Branch 33 → 34 not taken.
✓ Branch 33 → 35 taken 512 times.
✗ Branch 35 → 36 not taken.
✓ Branch 35 → 37 taken 512 times.
1024 allocate (this%total_nonlocal_nz_per_row(space%rank_l0:space%rank_l1))
144
5/10
✗ Branch 37 → 38 not taken.
✓ Branch 37 → 39 taken 512 times.
✓ Branch 39 → 38 taken 512 times.
✗ Branch 39 → 40 not taken.
✗ Branch 40 → 41 not taken.
✓ Branch 40 → 42 taken 512 times.
✗ Branch 42 → 43 not taken.
✓ Branch 42 → 44 taken 512 times.
✗ Branch 44 → 45 not taken.
✓ Branch 44 → 46 taken 512 times.
1024 allocate (this%total_local_nz_per_row_perp(space%rank_l0:space%rank_l1))
145
5/10
✗ Branch 46 → 47 not taken.
✓ Branch 46 → 48 taken 512 times.
✓ Branch 48 → 47 taken 512 times.
✗ Branch 48 → 49 not taken.
✗ Branch 49 → 50 not taken.
✓ Branch 49 → 51 taken 512 times.
✗ Branch 51 → 52 not taken.
✓ Branch 51 → 53 taken 512 times.
✗ Branch 53 → 54 not taken.
✓ Branch 53 → 55 taken 512 times.
1024 allocate (this%col_type(space%rank_l0:space%rank_l1))
146
5/10
✗ Branch 55 → 56 not taken.
✓ Branch 55 → 57 taken 512 times.
✓ Branch 57 → 56 taken 512 times.
✗ Branch 57 → 58 not taken.
✗ Branch 58 → 59 not taken.
✓ Branch 58 → 60 taken 512 times.
✗ Branch 60 → 61 not taken.
✓ Branch 60 → 62 taken 512 times.
✗ Branch 62 → 63 not taken.
✓ Branch 62 → 64 taken 512 times.
1024 allocate (this%col_mapping(space%rank_l0:space%rank_l1))
147
9/16
✗ Branch 64 → 65 not taken.
✓ Branch 64 → 66 taken 512 times.
✓ Branch 66 → 65 taken 512 times.
✗ Branch 66 → 67 not taken.
✓ Branch 67 → 68 taken 512 times.
✗ Branch 67 → 69 not taken.
✓ Branch 69 → 70 taken 512 times.
✗ Branch 69 → 71 not taken.
✗ Branch 71 → 72 not taken.
✓ Branch 71 → 73 taken 512 times.
✗ Branch 73 → 74 not taken.
✓ Branch 73 → 75 taken 512 times.
✗ Branch 75 → 76 not taken.
✓ Branch 75 → 77 taken 512 times.
✓ Branch 78 → 79 taken 674 times.
✓ Branch 78 → 80 taken 512 times.
2210 allocate (this%resp_inds(size(space%tp_spaces)))
148
149 512 this%nr_cols = 0
150 512 this%my_nr_cols = size(space, my_rank=.true.)
151
2/2
✓ Branch 81 → 82 taken 1035274 times.
✓ Branch 81 → 83 taken 512 times.
1035786 this%total_local_nz_per_row = 1 ! The default extraction/projection operator is the identity matrix
152
2/2
✓ Branch 84 → 85 taken 1035274 times.
✓ Branch 84 → 86 taken 512 times.
1035786 this%total_nonlocal_nz_per_row = 0
153
2/2
✓ Branch 87 → 88 taken 1035274 times.
✓ Branch 87 → 89 taken 512 times.
1035786 this%total_local_nz_per_row_perp = 0
154
2/2
✓ Branch 90 → 91 taken 1035274 times.
✓ Branch 90 → 92 taken 512 times.
1035786 this%col_type = 0 ! 0 means no constraint, > 0 means corresponding to constraints(col_type)
155
2/2
✓ Branch 92 → 93 taken 3072 times.
✓ Branch 92 → 94 taken 512 times.
3584 this%unknowns_are_on_this_rank = .true.
156 512 end subroutine init_mform_constraint
157
158 !> @brief Initialize the ConstraintLocalHolder object
159 !>
160 !> @param[inout] this The ConstraintLocalHolder object to initialize
161 !> @param[in] local_constraint The local constraint to wrap
162 !> @param[in] space The MFormSpace to which the local constraint applies
163 !> @param[in] coord_transform _(optional)_ The coordinate transformation associated with the m-form space
164 1402 subroutine init_constraint_local_holder(this, local_constraint, space, coord_transform)
165 use m_coord_transform_abstract, only: CoordTransformAbstract
166 class(ConstraintLocalHolder), intent(inout) :: this
167 class(MFormConstraintLocal), intent(in) :: local_constraint
168 type(MFormSpace), intent(in) :: space
169 class(CoordTransformAbstract), optional, intent(in) :: coord_transform
170
171
1/2
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 701 times.
701 if (allocated(this%elt)) then
172 call this%destroy()
173 end if
174
175
2/4
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 701 times.
✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 701 times.
701 allocate (this%elt, source=local_constraint)
176 701 call this%elt%init(space, coord_transform=coord_transform)
177 701 end subroutine init_constraint_local_holder
178
179 !> @brief Destroy the ConstraintLocalHolder object
180 !>
181 !> @param[inout] this The ConstraintLocalHolder object to destroy
182 4465 subroutine destroy_constraint_local_holder(this)
183 class(ConstraintLocalHolder), intent(inout) :: this
184
185
2/2
✓ Branch 2 → 3 taken 699 times.
✓ Branch 2 → 10 taken 3766 times.
4465 if (allocated(this%elt)) then
186 699 call this%elt%destroy()
187
2/4
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 699 times.
✓ Branch 6 → 7 taken 699 times.
✗ Branch 6 → 9 not taken.
1398 deallocate (this%elt)
188 end if
189 4465 end subroutine destroy_constraint_local_holder
190
191 !> @brief Insert a local constraint into the MFormConstraint object
192 !>
193 !> This subroutine does not actually insert the constraint into the projection matrix yet, but keeps a pointer to the constraint.
194 !> This is because the PETSc matrix assembly is more efficient if all of the nonzeros per row are known before inserting the
195 !> values, thus requiring us to first compute the number of nonzeros per row for all constraints, and then actually insert the
196 !> values into the PETSc matrix.
197 !>
198 !> @param[inout] this The MFormConstraint object in which the local constraint is inserted
199 !> @param[in] local_constraint The local constraint to insert
200 !> @param[in] coord_transform _(optional)_ The coordinate transformation associated with the m-form space
201 701 subroutine insert_in_mform_constraint(this, local_constraint, coord_transform)
202 use m_mform_basis, only: size
203 use m_coord_transform_abstract, only: CoordTransformAbstract
204 implicit none
205
206 class(MFormConstraint), intent(inout) :: this
207 class(MFormConstraintLocal), intent(in) :: local_constraint
208 class(CoordTransformAbstract), optional, intent(in) :: coord_transform
209
210 integer :: nr_blocks, blk_row, blk_col, local_row, row, row_start, row_end, ierr
211
212
1/2
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 701 times.
701 if (this%is_assembled) then
213 error stop "AbstractMatrix::insert_in_mform_constraint: matrix is already assembled"
214 end if
215
216 ! We don't actually insert the constraint into the projection matrix just yet
217 ! Instead, we keep a pointer to the constraint and insert it in the matrix when
218 ! all of the nonzeros per row are known, such that the PETSc matrix assembly
219 ! is efficient.
220
221 701 this%nr_constraints = this%nr_constraints + 1
222 701 call this%constraints(this%nr_constraints)%init(local_constraint, this%space, coord_transform=coord_transform)
223
224 701 select case (this%mode)
225 case (CONSTRAINT_PROJECTION)
226 error stop "MFormConstraint::insert_in_mform_constraint: projection constraints are not supported"
227 ! if (.not. allocated(this%constraints(this%nr_constraints)%elt%projector)) then
228 ! error stop "MFormConstraint::insert_in_mform_constraint: local constraint does not have a projector defined"
229 ! end if
230 ! if (any(this%constraints(this%nr_constraints)%elt%is_mpi_shared .and. &
231 ! this%space%tp_spaces(1)%domain%nr_subintervals > 0)) then
232 ! error stop "MFormConstraint::insert_in_mform_constraint: projection constraints cannot be mpi shared."
233 ! end if
234 case (CONSTRAINT_EXTRACTION)
235
1/2
✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 701 times.
701 if (.not. allocated(this%constraints(this%nr_constraints)%elt%extractor)) then
236 error stop "MFormConstraint::insert_in_mform_constraint: local constraint does not have an extractor defined"
237 end if
238
239 this%unknowns_are_on_this_rank(this%nr_constraints) = all(.not. this%constraints(this%nr_constraints)%elt%is_mpi_shared .or. &
240
4/9
✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 701 times.
✗ Branch 5 → 15 not taken.
✓ Branch 10 → 11 taken 2103 times.
✓ Branch 10 → 14 taken 701 times.
✗ Branch 11 → 12 not taken.
✓ Branch 11 → 13 taken 2103 times.
✗ Branch 12 → 13 not taken.
✗ Branch 12 → 14 not taken.
3505 this%space%tp_spaces(1)%domain%my_subinterval_ijk == 0)
241 end select
242
243 701 nr_blocks = size(this%space%tp_spaces)
244
245 ! Compute the nr of nonzeros per row
246
2/2
✓ Branch 16 → 17 taken 887 times.
✓ Branch 16 → 69 taken 701 times.
1588 do blk_row = 1, nr_blocks
247
248
2/2
✓ Branch 18 → 19 taken 198306 times.
✓ Branch 18 → 48 taken 887 times.
199193 do local_row = 0, size(this%constraints(this%nr_constraints)%elt%index_list(blk_row)%i) - 1
249 198306 row = this%constraints(this%nr_constraints)%elt%index_list(blk_row)%lin(local_row)
250 198306 this%total_local_nz_per_row(row) = 0
251 198306 this%total_local_nz_per_row_perp(row) = 1
252
1/2
✗ Branch 19 → 20 not taken.
✓ Branch 19 → 47 taken 198306 times.
198306 if (this%col_type(row) > 0) then
253 print*,'blk, row, i, j, k', blk_row, row, this%constraints(this%nr_constraints)%elt%index_list(blk_row)%i(local_row), &
254 & this%constraints(this%nr_constraints)%elt%index_list(blk_row)%j(local_row), &
255 & this%constraints(this%nr_constraints)%elt%index_list(blk_row)%k(local_row)
256 print*,'lbound, ubound', lbound(this%col_type), ubound(this%col_type)
257 print*,'rank', this%space%tp_spaces(1)%domain%my_rank
258 error stop "MFormConstraint::insert_in_mform_constraint: when using multiple constraints it must be ensured"// &
259 & " that they do not overlap"
260 end if
261 199193 this%col_type(row) = this%nr_constraints
262 end do
263
264 701 select case (this%mode)
265 case (CONSTRAINT_PROJECTION)
266 do blk_col = 1, nr_blocks
267 if (this%constraints(this%nr_constraints)%elt%projector(blk_row, blk_col)%nr_nonzero == 0) cycle
268 do local_row = 0, this%constraints(this%nr_constraints)%elt%projector(blk_row, blk_col)%nr_rows - 1
269
270 row_start = this%constraints(this%nr_constraints)%elt%projector(blk_row, blk_col)%row_starts_at_nz(local_row)
271 row_end = this%constraints(this%nr_constraints)%elt%projector(blk_row, blk_col)%row_starts_at_nz(local_row + 1) - 1
272
273 row = this%constraints(this%nr_constraints)%elt%index_list(blk_row)%lin(local_row)
274 this%total_local_nz_per_row(row) = this%total_local_nz_per_row(row) + 1 + row_end - row_start
275 this%total_local_nz_per_row_perp(row) = this%total_local_nz_per_row_perp(row) + 1 + row_end - row_start
276 end do
277 end do
278 case (CONSTRAINT_EXTRACTION)
279
2/2
✓ Branch 61 → 62 taken 198306 times.
✓ Branch 61 → 63 taken 887 times.
199193 do local_row = 0, size(this%constraints(this%nr_constraints)%elt%index_list(blk_row)%i) - 1
280 198306 row_start = this%constraints(this%nr_constraints)%elt%extractor(blk_row)%row_starts_at_nz(local_row)
281 198306 row_end = this%constraints(this%nr_constraints)%elt%extractor(blk_row)%row_starts_at_nz(local_row + 1) - 1
282
283 198306 row = this%constraints(this%nr_constraints)%elt%index_list(blk_row)%lin(local_row)
284
285 199193 this%total_local_nz_per_row(row) = this%total_local_nz_per_row(row) + 1 + row_end - row_start
286 end do
287
288 ! NOTE: the number of columns must be equal per block!
289
3/4
✓ Branch 64 → 65 taken 701 times.
✓ Branch 64 → 67 taken 186 times.
✓ Branch 65 → 66 taken 701 times.
✗ Branch 65 → 67 not taken.
887 if (blk_row == 1 .and. this%unknowns_are_on_this_rank(this%nr_constraints)) then
290 701 this%my_nr_cols = this%my_nr_cols + this%constraints(this%nr_constraints)%elt%extractor(blk_row)%nr_cols
291 end if
292
1/3
✗ Branch 49 → 50 not taken.
✓ Branch 49 → 60 taken 887 times.
✗ Branch 49 → 68 not taken.
1774 this%my_nr_cols = this%my_nr_cols - this%constraints(this%nr_constraints)%elt%extractor(blk_row)%nr_rows
293 end select
294 end do
295
296
2/2
✓ Branch 70 → 71 taken 412 times.
✓ Branch 70 → 72 taken 289 times.
701 if (.not. this%constraints(this%nr_constraints)%elt%is_symmetric) this%is_symmetric = .false.
297 701 end subroutine insert_in_mform_constraint
298
299 !> @brief Assemble the MFormConstraint object (i.e., assemble the PETSc matrix)
300 !>
301 !> @param[inout] this The MFormConstraint object to assemble
302 !>
303 !> @note Calling this subroutine is usually not needed as it is called automatically by the GenericSolver
304 512 subroutine assemble(this)
305 use m_tensorprod, only: size, TensorProdIndices
306 use m_common, only: wp
307 use m_mform_basis, only: size
308 implicit none
309
310 class(MFormConstraint), intent(inout) :: this
311
312 integer :: ierr, is_ignored, nonlocal_nz_per_row, cdx, row, col, blk, blk_col, nr_blocks
313 integer :: nr_unknowns_visisted_per_constraint(MAX_NR_CONSTRAINTS), nr_unknowns_per_constraint(MAX_NR_CONSTRAINTS)
314 integer :: local_unknown_counter, global_unknown0, nr_active(3)
315 512 integer, allocatable :: their_unknowns(:)
316 integer :: i, j, k, sze_x, sze_y, sze_z
317 logical :: unknown_is_kept, active_on_block_per_constraint(3, MAX_NR_CONSTRAINTS)
318 integer, parameter :: COL_UNSET = -1, COL_OTHER_RANK = -2
319
320
1/2
✓ Branch 2 → 3 taken 512 times.
✗ Branch 2 → 48 not taken.
512 if (this%is_assembled) return
321
322 512 is_ignored = 0 ! Normally would use PETSC_NULL_INTEGER, but this has changed in different PETSc versions
323
324 512 nr_blocks = size(this%space%tp_spaces)
325 512 select case (this%mode)
326 case (CONSTRAINT_PROJECTION)
327 call init_resp_indexlist_and_colmapping_projector()
328 case (CONSTRAINT_EXTRACTION)
329 512 call init_resp_indexlist_and_colmapping_extractor()
330
331 ! At this point this%col_mapping is initialized, except for those columns for which another rank is responsible (by means of
332 ! a constraint which has is_mpi_shared, unknowns_are_on_this_rank(cdx) == .false.). Hence, we now obtain the mapping for those
333 ! columns by communicating with the other ranks.
334
2/2
✓ Branch 7 → 8 taken 701 times.
✓ Branch 7 → 16 taken 512 times.
1213 do cdx = 1, this%nr_constraints
335
4/8
✓ Branch 9 → 10 taken 2103 times.
✓ Branch 9 → 13 taken 701 times.
✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 2103 times.
✗ Branch 11 → 12 not taken.
✗ Branch 11 → 13 not taken.
✗ Branch 13 → 14 not taken.
✓ Branch 13 → 15 taken 701 times.
2804 if (any(this%constraints(cdx)%elt%is_mpi_shared .and. &
336 512 this%space%tp_spaces(1)%domain%nr_subintervals > 1)) then
337 call insert_offrank_colmapping_extractor(cdx)
338 end if
339 end do
340
341
5/9
✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 512 times.
✗ Branch 3 → 23 not taken.
✓ Branch 18 → 19 taken 1035274 times.
✓ Branch 18 → 21 taken 512 times.
✓ Branch 19 → 20 taken 1035274 times.
✗ Branch 19 → 21 not taken.
✗ Branch 21 → 22 not taken.
✓ Branch 21 → 23 taken 512 times.
1036298 if (any(this%col_mapping == COL_OTHER_RANK)) then
342 error stop "MFormConstraint::assemble: some columns still have COL_OTHER_RANK after communication"
343 end if
344 end select
345
346 512 call init_petsc_matrices()
347
348
2/2
✓ Branch 25 → 26 taken 701 times.
✓ Branch 25 → 32 taken 512 times.
1213 do cdx = 1, this%nr_constraints
349 512 select case (this%mode)
350 case (CONSTRAINT_PROJECTION)
351 call this%insert_constraint_in_projection_matrix(this%constraints(cdx)%elt)
352 case (CONSTRAINT_EXTRACTION)
353
1/3
✗ Branch 26 → 27 not taken.
✓ Branch 26 → 29 taken 701 times.
✗ Branch 26 → 31 not taken.
701 call this%insert_constraint_in_extraction_matrix(this%constraints(cdx)%elt, active_on_block_per_constraint(:, cdx))
354 end select
355 end do
356
357
1/2
✗ Branch 34 → 35 not taken.
✓ Branch 34 → 36 taken 512 times.
512 PetscCall(MatAssemblyBegin(this%mat, MAT_FINAL_ASSEMBLY, ierr))
358
1/2
✗ Branch 37 → 38 not taken.
✓ Branch 37 → 39 taken 512 times.
512 PetscCall(MatAssemblyEnd(this%mat, MAT_FINAL_ASSEMBLY, ierr))
359
360
1/2
✗ Branch 39 → 40 not taken.
✓ Branch 39 → 46 taken 512 times.
512 if (this%mode == CONSTRAINT_PROJECTION) then
361 PetscCall(MatAssemblyBegin(this%mat_perp, MAT_FINAL_ASSEMBLY, ierr))
362 PetscCall(MatAssemblyEnd(this%mat_perp, MAT_FINAL_ASSEMBLY, ierr))
363 end if
364
365
1/4
✗ Branch 46 → 47 not taken.
✓ Branch 46 → 52 taken 512 times.
✗ Branch 48 → 49 not taken.
✗ Branch 48 → 50 not taken.
512 this%is_assembled = .true.
366
367 contains
368 pure logical function col_was_set(col)
369
370 integer, intent(in) :: col
371
372 1035274 col_was_set = col > COL_UNSET
373 end function col_was_set
374
375 subroutine init_resp_indexlist_and_colmapping_projector()
376 implicit none
377
378 do row = this%space%rank_l0, this%space%rank_l1
379 this%col_mapping(row) = row
380 end do
381 do blk_col = 1, nr_blocks
382 call this%resp_inds(blk_col)%init(this%space%tp_spaces(blk_col)%rank_resp_bspline, &
383 this%space%tp_spaces(blk_col)%rank_resp_bspline, this%space%tp_spaces(blk_col)%rank_l0)
384 end do
385
386 this%nr_cols = size(this%space)
387
388 end subroutine init_resp_indexlist_and_colmapping_projector
389
390 512 subroutine init_resp_indexlist_and_colmapping_extractor()
391 implicit none
392
393 ! The vector this%col_mapping maps indices of the original matrix
394 ! to the indices of the new smaller matrix
395 ! TensorProdIndicesList are untouched by the extractor are kept (this%col_mapping > 0) and
396 ! for each constraint the first constraint%nr_cols indices are kept as well
397
398 ! Determine the number of unknowns corresponding to each constraint
399 512 nr_unknowns_per_constraint = 0
400
2/2
✓ Branch 3 → 4 taken 701 times.
✓ Branch 3 → 5 taken 512 times.
1213 do cdx = 1, this%nr_constraints
401 ! NOTE: the number of columns must be equal per block!
402 1213 nr_unknowns_per_constraint(cdx) = this%constraints(cdx)%elt%extractor(1)%nr_cols
403 end do
404
405 ! We gather the number of unknowns on all ranks such that we can define a global index
406
7/14
✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 512 times.
✓ Branch 8 → 7 taken 512 times.
✗ Branch 8 → 9 not taken.
✓ Branch 9 → 10 taken 512 times.
✗ Branch 9 → 11 not taken.
✓ Branch 11 → 12 taken 512 times.
✗ Branch 11 → 13 not taken.
✗ Branch 13 → 14 not taken.
✓ Branch 13 → 15 taken 512 times.
✗ Branch 15 → 16 not taken.
✓ Branch 15 → 17 taken 512 times.
✗ Branch 17 → 18 not taken.
✓ Branch 17 → 19 taken 512 times.
1536 allocate (their_unknowns(0:this%space%tp_spaces(1)%domain%nr_ranks - 1))
407 call MPI_Allgather(this%my_nr_cols, 1, MPI_INTEGER, &
408 & their_unknowns, 1, MPI_INTEGER, &
409
1/2
✗ Branch 20 → 21 not taken.
✓ Branch 20 → 22 taken 512 times.
512 & this%space%tp_spaces(1)%domain%comm, ierr); CHKERRMPI(ierr)
410
411
1/2
✗ Branch 23 → 24 not taken.
✓ Branch 23 → 25 taken 512 times.
512 global_unknown0 = sum(their_unknowns(0:this%space%tp_spaces(1)%domain%my_rank - 1))
412
2/2
✓ Branch 26 → 27 taken 512 times.
✓ Branch 26 → 28 taken 512 times.
1024 this%nr_cols = sum(their_unknowns)
413
1/2
✗ Branch 28 → 29 not taken.
✓ Branch 28 → 30 taken 512 times.
512 deallocate (their_unknowns)
414
415 ! We use the global index to define col_mapping
416 512 nr_active = 0
417
2/2
✓ Branch 31 → 32 taken 1035274 times.
✓ Branch 31 → 33 taken 512 times.
1035786 this%col_mapping = COL_UNSET
418 512 active_on_block_per_constraint = .false.
419 512 nr_unknowns_visisted_per_constraint = 0
420 512 local_unknown_counter = 0
421
2/2
✓ Branch 34 → 35 taken 674 times.
✓ Branch 34 → 48 taken 512 times.
1186 do blk_col = 1, nr_blocks
422
2/2
✓ Branch 36 → 37 taken 1035274 times.
✓ Branch 36 → 46 taken 674 times.
1036460 do col = this%space%tp_spaces(blk_col)%rank_l0, this%space%tp_spaces(blk_col)%rank_l1
423 1035274 cdx = this%col_type(col)
424
425 ! An unknowns is kept if it is not part of a constraint (cdx == 0) or if it is part of a constraint
426 ! and the number of unknowns visited for that constraint is no more than the number of unknowns per constraint
427 1035274 unknown_is_kept = cdx == 0
428
2/2
✓ Branch 37 → 38 taken 198306 times.
✓ Branch 37 → 39 taken 836968 times.
1035274 if (cdx > 0) then
429 198306 nr_unknowns_visisted_per_constraint(cdx) = nr_unknowns_visisted_per_constraint(cdx) + 1
430 198306 unknown_is_kept = nr_unknowns_visisted_per_constraint(cdx) <= nr_unknowns_per_constraint(cdx)
431 end if
432
433
2/2
✓ Branch 39 → 40 taken 854633 times.
✓ Branch 39 → 45 taken 180641 times.
1035948 if (unknown_is_kept) then
434
1/2
✓ Branch 40 → 41 taken 854633 times.
✗ Branch 40 → 42 not taken.
854633 if (this%unknowns_are_on_this_rank(cdx)) then
435 854633 this%col_mapping(col) = global_unknown0 + local_unknown_counter
436 854633 nr_active(blk_col) = nr_active(blk_col) + 1
437 854633 local_unknown_counter = local_unknown_counter + 1
438 else
439 this%col_mapping(col) = COL_OTHER_RANK ! Indicate that another rank is responsible for this unknown
440 end if
441
2/2
✓ Branch 43 → 44 taken 17665 times.
✓ Branch 43 → 45 taken 836968 times.
854633 if (cdx > 0) active_on_block_per_constraint(blk_col, cdx) = .true.
442 end if
443 end do
444 end do
445
446 ! Update the total_nonlocal_nz_per_row to account for the fact that some columns are not present on this rank
447
4/6
✓ Branch 50 → 51 taken 1035274 times.
✓ Branch 50 → 53 taken 512 times.
✓ Branch 51 → 52 taken 1035274 times.
✗ Branch 51 → 53 not taken.
✗ Branch 53 → 54 not taken.
✓ Branch 53 → 59 taken 512 times.
1035786 if (any(this%col_mapping == COL_OTHER_RANK)) then
448 ! TODO: optimize this
449 where (this%col_type > 0) this%total_nonlocal_nz_per_row = this%total_local_nz_per_row
450 end if
451
452 ! Make a list of responsible tensor product indices for the kept unknowns on this rank
453
2/2
✓ Branch 60 → 61 taken 674 times.
✓ Branch 60 → 145 taken 512 times.
1186 do blk = 1, size(this%resp_inds)
454
455
1/2
✗ Branch 61 → 62 not taken.
✓ Branch 61 → 63 taken 674 times.
674 if (allocated(this%resp_inds(blk)%lin)) deallocate (this%resp_inds(blk)%lin)
456
1/2
✗ Branch 63 → 64 not taken.
✓ Branch 63 → 65 taken 674 times.
674 if (allocated(this%resp_inds(blk)%i)) deallocate (this%resp_inds(blk)%i)
457
1/2
✗ Branch 65 → 66 not taken.
✓ Branch 65 → 67 taken 674 times.
674 if (allocated(this%resp_inds(blk)%j)) deallocate (this%resp_inds(blk)%j)
458
1/2
✗ Branch 67 → 68 not taken.
✓ Branch 67 → 69 taken 674 times.
674 if (allocated(this%resp_inds(blk)%k)) deallocate (this%resp_inds(blk)%k)
459
460
9/14
✓ Branch 69 → 70 taken 46 times.
✓ Branch 69 → 71 taken 628 times.
✓ Branch 71 → 70 taken 628 times.
✗ Branch 71 → 72 not taken.
✓ Branch 72 → 73 taken 674 times.
✗ Branch 72 → 74 not taken.
✓ Branch 74 → 75 taken 628 times.
✓ Branch 74 → 76 taken 46 times.
✗ Branch 76 → 77 not taken.
✓ Branch 76 → 78 taken 674 times.
✗ Branch 78 → 79 not taken.
✓ Branch 78 → 80 taken 674 times.
✗ Branch 80 → 81 not taken.
✓ Branch 80 → 82 taken 674 times.
2022 allocate (this%resp_inds(blk)%lin(0:nr_active(blk) - 1))
461
9/14
✓ Branch 82 → 83 taken 46 times.
✓ Branch 82 → 84 taken 628 times.
✓ Branch 84 → 83 taken 628 times.
✗ Branch 84 → 85 not taken.
✓ Branch 85 → 86 taken 674 times.
✗ Branch 85 → 87 not taken.
✓ Branch 87 → 88 taken 628 times.
✓ Branch 87 → 89 taken 46 times.
✗ Branch 89 → 90 not taken.
✓ Branch 89 → 91 taken 674 times.
✗ Branch 91 → 92 not taken.
✓ Branch 91 → 93 taken 674 times.
✗ Branch 93 → 94 not taken.
✓ Branch 93 → 95 taken 674 times.
2022 allocate (this%resp_inds(blk)%i(0:nr_active(blk) - 1))
462
9/14
✓ Branch 95 → 96 taken 46 times.
✓ Branch 95 → 97 taken 628 times.
✓ Branch 97 → 96 taken 628 times.
✗ Branch 97 → 98 not taken.
✓ Branch 98 → 99 taken 674 times.
✗ Branch 98 → 100 not taken.
✓ Branch 100 → 101 taken 628 times.
✓ Branch 100 → 102 taken 46 times.
✗ Branch 102 → 103 not taken.
✓ Branch 102 → 104 taken 674 times.
✗ Branch 104 → 105 not taken.
✓ Branch 104 → 106 taken 674 times.
✗ Branch 106 → 107 not taken.
✓ Branch 106 → 108 taken 674 times.
2022 allocate (this%resp_inds(blk)%j(0:nr_active(blk) - 1))
463
9/14
✓ Branch 108 → 109 taken 46 times.
✓ Branch 108 → 110 taken 628 times.
✓ Branch 110 → 109 taken 628 times.
✗ Branch 110 → 111 not taken.
✓ Branch 111 → 112 taken 674 times.
✗ Branch 111 → 113 not taken.
✓ Branch 113 → 114 taken 628 times.
✓ Branch 113 → 115 taken 46 times.
✗ Branch 115 → 116 not taken.
✓ Branch 115 → 117 taken 674 times.
✗ Branch 117 → 118 not taken.
✓ Branch 117 → 119 taken 674 times.
✗ Branch 119 → 120 not taken.
✓ Branch 119 → 121 taken 674 times.
2022 allocate (this%resp_inds(blk)%k(0:nr_active(blk) - 1))
464
465 674 sze_x = size(this%space%tp_spaces(blk), 1, my_rank=.true.)
466 674 sze_y = size(this%space%tp_spaces(blk), 2, my_rank=.true.)
467 674 sze_z = size(this%space%tp_spaces(blk), 3, my_rank=.true.)
468
469 674 local_unknown_counter = 0
470
2/2
✓ Branch 122 → 123 taken 3586 times.
✓ Branch 122 → 134 taken 674 times.
4260 do k = this%space%tp_spaces(blk)%rank_resp_bspline%k0, this%space%tp_spaces(blk)%rank_resp_bspline%k1
471
2/2
✓ Branch 124 → 125 taken 55574 times.
✓ Branch 124 → 132 taken 3586 times.
59834 do j = this%space%tp_spaces(blk)%rank_resp_bspline%j0, this%space%tp_spaces(blk)%rank_resp_bspline%j1
472
2/2
✓ Branch 126 → 127 taken 1035274 times.
✓ Branch 126 → 130 taken 55574 times.
1094434 do i = this%space%tp_spaces(blk)%rank_resp_bspline%i0, this%space%tp_spaces(blk)%rank_resp_bspline%i1
473 col = this%space%tp_spaces(blk)%rank_l0 + &
474 (i - this%space%tp_spaces(blk)%rank_resp_bspline%i0) + &
475 (j - this%space%tp_spaces(blk)%rank_resp_bspline%j0) * sze_x + &
476 1035274 (k - this%space%tp_spaces(blk)%rank_resp_bspline%k0) * sze_x * sze_y
477
2/2
✓ Branch 127 → 128 taken 854633 times.
✓ Branch 127 → 129 taken 180641 times.
1090848 if (col_was_set(this%col_mapping(col))) then
478 854633 this%resp_inds(blk)%lin(local_unknown_counter) = col
479 854633 this%resp_inds(blk)%i(local_unknown_counter) = i
480 854633 this%resp_inds(blk)%j(local_unknown_counter) = j
481 854633 this%resp_inds(blk)%k(local_unknown_counter) = k
482 854633 local_unknown_counter = local_unknown_counter + 1
483 end if
484 end do
485 end do
486 end do
487
1/2
✗ Branch 135 → 136 not taken.
✓ Branch 135 → 144 taken 674 times.
1186 if (local_unknown_counter /= nr_active(blk)) then
488 print*,'blk, sze, nr_active, local_unknown_counter', blk, sze_x * sze_y * sze_z, nr_active(blk), local_unknown_counter
489 error stop 'ERROR: MFormConstraint::assemble: number of active indices does not match the expected number'
490 end if
491 end do
492 end subroutine init_resp_indexlist_and_colmapping_extractor
493
494 subroutine insert_offrank_colmapping_extractor(cdx)
495 use m_common, only: cartesian_to_linear
496 implicit none
497
498 integer, intent(in) :: cdx
499
500 integer, allocatable :: subinterval_ijk(:, :)
501 integer :: leader_rank, leader_l0, leader_l1, leader_row, leader_ls(2)
502 integer, allocatable :: leader_col_mapping(:), leader_col_type(:)
503
504 integer, allocatable :: constraint_inds(:), constraint_size(:), non_constraint_inds(:), non_constraint_size(:)
505 integer :: constraint_comm, constraint_rank, constraint_nr_ranks, comm_color, comm_key
506
507 ! Create a sub communicator for the ranks involved in the constraint that share mpi unknowns
508 constraint_inds = pack(this%space%tp_spaces(1)%domain%my_subinterval_ijk, this%constraints(cdx)%elt%is_mpi_shared)
509 constraint_size = pack(this%space%tp_spaces(1)%domain%nr_subintervals, this%constraints(cdx)%elt%is_mpi_shared)
510 non_constraint_inds = pack(this%space%tp_spaces(1)%domain%my_subinterval_ijk,.not. this%constraints(cdx)%elt%is_mpi_shared)
511 non_constraint_size = pack(this%space%tp_spaces(1)%domain%nr_subintervals,.not. this%constraints(cdx)%elt%is_mpi_shared)
512
513 comm_color = cartesian_to_linear(non_constraint_inds, non_constraint_size) ! Same color ends up on same comm
514 comm_key = cartesian_to_linear(constraint_inds, constraint_size) ! Key to order ranks within comm
515
516 PetscCallMPI(MPI_Comm_split(this%space%tp_spaces(1)%domain%comm, comm_color, comm_key, constraint_comm, ierr))
517 PetscCallMPI(MPI_Comm_rank(constraint_comm, constraint_rank, ierr))
518 PetscCallMPI(MPI_Comm_size(constraint_comm, constraint_nr_ranks, ierr))
519
520 ! MPI communicate the subinterval indices as function of the rank
521 allocate (subinterval_ijk(1:size(constraint_inds), 0:constraint_nr_ranks - 1))
522
523 call MPI_Allgather(constraint_inds, size(constraint_inds), MPI_INTEGER, &
524 subinterval_ijk, size(constraint_inds), MPI_INTEGER, &
525 constraint_comm, ierr); CHKERRMPI(ierr)
526
527 ! Determine the leader rank for the constraint
528 leader_rank = -1
529 do i = 0, constraint_nr_ranks - 1
530 if (all(subinterval_ijk(:, i) == 0)) then
531 leader_rank = i
532 exit
533 end if
534 end do
535 if (leader_rank == -1) then
536 error stop "MFormConstraint::insert_offrank_colmapping_extractor: could not determine leader rank"
537 end if
538
539 if (constraint_rank == leader_rank) then
540 ! First, we send rank_l0 and rank_l1 to the other ranks
541 leader_l0 = this%space%rank_l0
542 leader_l1 = this%space%rank_l1
543
544 PetscCallMPI(MPI_Bcast([leader_l0, leader_l1], 2, MPI_INTEGER, leader_rank, constraint_comm, ierr))
545
546 ! Next, we send the col_mapping to the other ranks
547 PetscCallMPI(MPI_Bcast(this%col_mapping, leader_l1 - leader_l0 + 1, MPI_INTEGER, leader_rank, constraint_comm, ierr))
548
549 ! as well as the col_type
550 PetscCallMPI(MPI_Bcast(this%col_type, leader_l1 - leader_l0 + 1, MPI_INTEGER, leader_rank, constraint_comm, ierr))
551 else
552 ! First, we receive rank_l0 and rank_l1 from the leader rank
553 PetscCallMPI(MPI_Bcast(leader_ls, 2, MPI_INTEGER, leader_rank, constraint_comm, ierr))
554 leader_l0 = leader_ls(1); leader_l1 = leader_ls(2)
555
556 ! Next, we receive the col_mapping from the leader rank
557 allocate (leader_col_mapping(leader_l0:leader_l1))
558 PetscCallMPI(MPI_Bcast(leader_col_mapping, leader_l1 - leader_l0 + 1, MPI_INTEGER, leader_rank, constraint_comm, ierr))
559 ! as well as the col_type
560 allocate (leader_col_type(leader_l0:leader_l1))
561 PetscCallMPI(MPI_Bcast(leader_col_type, leader_l1 - leader_l0 + 1, MPI_INTEGER, leader_rank, constraint_comm, ierr))
562
563 ! Finally, we insert the received col_mapping into our own col_mapping: match the n-th constraint unknown missing on this
564 ! rank with the n-th constraint unknown on the leader rank
565 leader_row = leader_l0
566 do row = this%space%rank_l0, this%space%rank_l1
567 if (this%col_mapping(row) == COL_OTHER_RANK .and. this%col_type(row) == cdx) then
568 do while (leader_col_type(leader_row) /= cdx)
569 leader_row = leader_row + 1
570 if (leader_row > leader_l1) then
571 error stop "MFormConstraint::insert_offrank_colmapping_extractor: could not find matching constraint"// &
572 " unknown on leader rank"
573 end if
574 end do
575 this%col_mapping(row) = leader_col_mapping(leader_row)
576 leader_row = leader_row + 1
577 end if
578 end do
579
580 deallocate (leader_col_mapping)
581 deallocate (leader_col_type)
582 deallocate (subinterval_ijk)
583 end if
584
585 PetscCallMPI(MPI_Comm_free(constraint_comm, ierr))
586 end subroutine insert_offrank_colmapping_extractor
587
588 512 subroutine init_petsc_matrices()
589 implicit none
590
591 ! The projection matrix is initialized
592 call MatCreateAIJ(this%space%tp_spaces(1)%domain%comm, size(this%space, my_rank=.true.), this%my_nr_cols, size(this%space), &
593 512 & this%nr_cols, is_ignored, this%total_local_nz_per_row, is_ignored, this%total_nonlocal_nz_per_row, this%mat, ierr)
594
1/2
✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 512 times.
512 CHKERRQ(ierr)
595
596
1/2
✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 512 times.
512 PetscCall(MatSetFromOptions(this%mat, ierr))
597
1/2
✗ Branch 9 → 10 not taken.
✓ Branch 9 → 11 taken 512 times.
512 PetscCall(MatSetUp(this%mat, ierr))
598
1/2
✗ Branch 12 → 13 not taken.
✓ Branch 12 → 14 taken 512 times.
512 PetscCall(MatSetOption(this%mat, MAT_SYMMETRIC, merge(PETSC_TRUE, PETSC_FALSE, this%is_symmetric), ierr))
599
600
1/2
✗ Branch 14 → 15 not taken.
✓ Branch 14 → 27 taken 512 times.
512 if (this%mode == CONSTRAINT_PROJECTION) then
601 nonlocal_nz_per_row = 0 ! NOTE: the projector acts on local indices of this rank only
602 call MatCreateAIJ(this%space%tp_spaces(1)%domain%comm, size(this%space, my_rank=.true.), size(this%space, my_rank=.true.), &
603 & size(this%space), size(this%space), is_ignored, this%total_local_nz_per_row_perp, nonlocal_nz_per_row, &
604 & PETSC_NULL_INTEGER_ARRAY, this%mat_perp, ierr); CHKERRQ(ierr)
605
606 PetscCall(MatSetFromOptions(this%mat_perp, ierr))
607 PetscCall(MatSetUp(this%mat_perp, ierr))
608 PetscCall(MatSetOption(this%mat_perp, MAT_SYMMETRIC, merge(PETSC_TRUE, PETSC_FALSE, this%is_symmetric), ierr))
609 end if
610
611
2/2
✓ Branch 28 → 29 taken 1035274 times.
✓ Branch 28 → 51 taken 512 times.
1035786 do row = this%space%rank_l0, this%space%rank_l1
612 1035274 col = this%col_mapping(row)
613
2/2
✓ Branch 29 → 30 taken 854633 times.
✓ Branch 29 → 50 taken 180641 times.
1035786 if (col > -1) then
614
2/2
✓ Branch 30 → 31 taken 836968 times.
✓ Branch 30 → 40 taken 17665 times.
854633 if (this%col_type(row) == 0) then
615
5/6
✓ Branch 32 → 33 taken 836968 times.
✓ Branch 32 → 34 taken 836968 times.
✓ Branch 35 → 36 taken 836968 times.
✓ Branch 35 → 37 taken 836968 times.
✗ Branch 38 → 39 not taken.
✓ Branch 38 → 50 taken 836968 times.
2510904 PetscCall(MatSetValues(this%mat, 1, [row], 1, [col], [1._wp], INSERT_VALUES, ierr))
616
1/2
✗ Branch 40 → 41 not taken.
✓ Branch 40 → 50 taken 17665 times.
17665 else if (this%mode == CONSTRAINT_PROJECTION) then
617 PetscCall(MatSetValues(this%mat_perp, 1, [row], 1, [row], [1._wp], INSERT_VALUES, ierr))
618 end if
619 end if
620 end do
621 end subroutine init_petsc_matrices
622 end subroutine assemble
623
624 !> @brief Insert a local constraint into the projection matrix
625 !>
626 !> @param[inout] this The MFormConstraint object in which the local constraint is inserted
627 !> @param[in] local_constraint The local constraint to insert
628 subroutine insert_constraint_in_projection_matrix(this, local_constraint)
629 use m_common, only: wp
630 use m_mform_basis, only: size
631 implicit none
632
633 class(MFormConstraint), intent(inout) :: this
634 class(MFormConstraintLocal), intent(in) :: local_constraint
635
636 integer :: ierr, blk_row, blk_col, local_row, local_idx, local_col
637 integer :: row, nr_blocks, row_start, row_end, nr_values
638 integer, allocatable :: cols(:)
639 real(wp), allocatable :: local_vals(:)
640
641 nr_blocks = size(this%space%tp_spaces)
642
643 ! And we insert the local constraint projector
644 do blk_row = 1, nr_blocks
645 do blk_col = 1, nr_blocks
646
647 if (local_constraint%projector(blk_row, blk_col)%nr_nonzero == 0) cycle
648
649 allocate (cols(0:local_constraint%projector(blk_row, blk_col)%nr_rows - 1))
650 allocate (local_vals(0:local_constraint%projector(blk_row, blk_col)%nr_rows - 1))
651 do local_row = 0, local_constraint%projector(blk_row, blk_col)%nr_rows - 1
652 row = local_constraint%index_list(blk_row)%lin(local_row)
653
654 row_start = local_constraint%projector(blk_row, blk_col)%row_starts_at_nz(local_row)
655 row_end = local_constraint%projector(blk_row, blk_col)%row_starts_at_nz(local_row + 1) - 1
656
657 nr_values = 1 + row_end - row_start
658
659 if (nr_values == 0) cycle
660
661 do local_idx = 0, nr_values - 1
662 local_col = local_constraint%projector(blk_row, blk_col)%col_idx(row_start + local_idx)
663 cols(local_idx) = local_constraint%index_list(blk_col)%lin(local_col)
664 local_vals(local_idx) = local_constraint%projector(blk_row, blk_col)%values(row_start + local_idx)
665 end do
666
667 PetscCall(MatSetValues(this%mat, 1, [row], nr_values, cols(0:nr_values - 1), local_vals(0:nr_values - 1), INSERT_VALUES, \
668 ierr))
669
670 do local_idx = 0, nr_values - 1
671 if (cols(local_idx) == row) then
672 local_vals(local_idx) = 1._wp - local_vals(local_idx)
673 else
674 local_vals(local_idx) = -local_vals(local_idx)
675 end if
676 end do
677 PetscCall(MatSetValues(this%mat_perp, 1, [row], nr_values, cols(0:nr_values - 1), local_vals(0:nr_values - 1), \
678 INSERT_VALUES, ierr))
679
680 end do
681
682 deallocate (cols)
683 deallocate (local_vals)
684 end do
685 end do
686
687 end subroutine insert_constraint_in_projection_matrix
688
689 !> @brief Insert a local constraint into the extraction matrix
690 !>
691 !> @param[inout] this The MFormConstraint object in which the local constraint is inserted
692 !> @param[in] local_constraint The local constraint to insert
693 !> @param[in] active_on_block A logical array indicating whether the constraint is active on each block
694 701 subroutine insert_constraint_in_extraction_matrix(this, local_constraint, active_on_block)
695 use m_common, only: wp
696 use m_mform_basis, only: size
697 implicit none
698
699 class(MFormConstraint), intent(inout) :: this
700 class(MFormConstraintLocal), intent(in) :: local_constraint
701 logical, intent(in) :: active_on_block(3)
702
703 integer :: ierr, blk, local_row, local_idx, local_col, local_blk, local_blk0(3), count
704 integer :: row, nr_blocks, row_start, row_end, nr_values
705 integer, allocatable :: cols(:)
706 real(wp), allocatable :: local_vals(:)
707
708 701 nr_blocks = size(this%space%tp_spaces)
709
710 count = 0
711
2/2
✓ Branch 3 → 4 taken 887 times.
✓ Branch 3 → 7 taken 701 times.
1588 do blk = 1, nr_blocks
712
2/2
✓ Branch 4 → 5 taken 472 times.
✓ Branch 4 → 6 taken 415 times.
1588 if (active_on_block(blk)) then
713 472 count = count + 1
714 472 local_blk0(count) = blk
715 end if
716 end do
717
718 ! And we insert the local constraint projector
719
2/2
✓ Branch 8 → 9 taken 887 times.
✓ Branch 8 → 58 taken 701 times.
1588 do blk = 1, nr_blocks
720
721
2/2
✓ Branch 9 → 10 taken 335 times.
✓ Branch 9 → 11 taken 552 times.
887 if (local_constraint%extractor(blk)%nr_nonzero == 0) cycle
722
723
6/12
✗ Branch 11 → 12 not taken.
✓ Branch 11 → 13 taken 552 times.
✓ Branch 13 → 12 taken 552 times.
✗ Branch 13 → 14 not taken.
✓ Branch 14 → 15 taken 552 times.
✗ Branch 14 → 16 not taken.
✓ Branch 16 → 17 taken 552 times.
✗ Branch 16 → 18 not taken.
✗ Branch 18 → 19 not taken.
✓ Branch 18 → 20 taken 552 times.
✗ Branch 20 → 21 not taken.
✓ Branch 20 → 22 taken 552 times.
1656 allocate (cols(0:local_constraint%extractor(blk)%nr_cols - 1))
724
6/12
✗ Branch 22 → 23 not taken.
✓ Branch 22 → 24 taken 552 times.
✓ Branch 24 → 23 taken 552 times.
✗ Branch 24 → 25 not taken.
✓ Branch 25 → 26 taken 552 times.
✗ Branch 25 → 27 not taken.
✓ Branch 27 → 28 taken 552 times.
✗ Branch 27 → 29 not taken.
✗ Branch 29 → 30 not taken.
✓ Branch 29 → 31 taken 552 times.
✗ Branch 31 → 32 not taken.
✓ Branch 31 → 33 taken 552 times.
1656 allocate (local_vals(0:local_constraint%extractor(blk)%nr_cols - 1))
725
2/2
✓ Branch 35 → 36 taken 147848 times.
✓ Branch 35 → 51 taken 552 times.
148400 do local_row = 0, local_constraint%extractor(blk)%nr_rows - 1
726 147848 row = local_constraint%index_list(blk)%lin(local_row)
727
728 147848 row_start = local_constraint%extractor(blk)%row_starts_at_nz(local_row)
729 147848 row_end = local_constraint%extractor(blk)%row_starts_at_nz(local_row + 1) - 1
730
731 147848 nr_values = 1 + row_end - row_start
732
733
2/2
✓ Branch 36 → 34 taken 11658 times.
✓ Branch 36 → 37 taken 136190 times.
147848 if (nr_values == 0) cycle
734
735
2/2
✓ Branch 37 → 38 taken 994310 times.
✓ Branch 37 → 42 taken 136190 times.
1130500 do local_idx = 0, nr_values - 1
736 994310 local_col = local_constraint%extractor(blk)%col_idx(row_start + local_idx)
737
738 994310 local_blk = local_blk0(1)
739
1/2
✗ Branch 38 → 39 not taken.
✓ Branch 38 → 41 taken 994310 times.
994310 if (local_col > local_constraint%extractor(local_blk)%nr_rows - 1) then
740 local_blk = local_blk0(2)
741 local_col = local_col - local_constraint%extractor(local_blk0(1))%nr_rows
742 if (local_col > local_constraint%extractor(local_blk0(2))%nr_rows - 1) then
743 local_blk = local_blk0(3)
744 local_col = local_col - local_constraint%extractor(local_blk0(2))%nr_rows
745 end if
746 end if
747
748 994310 cols(local_idx) = this%col_mapping(local_constraint%index_list(local_blk)%lin(local_col))
749 1130500 local_vals(local_idx) = local_constraint%extractor(blk)%values(row_start + local_idx)
750 end do
751
752
3/4
✓ Branch 44 → 45 taken 136190 times.
✓ Branch 44 → 46 taken 136190 times.
✓ Branch 47 → 34 taken 136190 times.
✗ Branch 47 → 48 not taken.
272932 PetscCall(MatSetValues(this%mat, 1, [row], nr_values, cols(0:nr_values - 1), local_vals(0:nr_values - 1), INSERT_VALUES, \
753 ierr))
754 end do
755
756
1/2
✗ Branch 52 → 53 not taken.
✓ Branch 52 → 54 taken 552 times.
552 deallocate (cols)
757
1/2
✗ Branch 54 → 55 not taken.
✓ Branch 54 → 56 taken 552 times.
1588 deallocate (local_vals)
758 end do
759
760 end subroutine insert_constraint_in_extraction_matrix
761
762 !> @brief Destroy the MFormConstraint object
763 !>
764 !> @param[inout] this The MFormConstraint object to destroy
765 893 subroutine destroy_mform_constraint(this)
766 use m_common, only: wp
767 implicit none
768 class(MFormConstraint), intent(inout) :: this
769
770 integer :: ierr, s1
771
772
2/2
✓ Branch 2 → 3 taken 510 times.
✓ Branch 2 → 4 taken 383 times.
893 if (allocated(this%total_local_nz_per_row)) deallocate (this%total_local_nz_per_row)
773
2/2
✓ Branch 4 → 5 taken 510 times.
✓ Branch 4 → 6 taken 383 times.
893 if (allocated(this%total_nonlocal_nz_per_row)) deallocate (this%total_nonlocal_nz_per_row)
774
2/2
✓ Branch 6 → 7 taken 510 times.
✓ Branch 6 → 8 taken 383 times.
893 if (allocated(this%total_local_nz_per_row_perp)) deallocate (this%total_local_nz_per_row_perp)
775
2/2
✓ Branch 8 → 9 taken 510 times.
✓ Branch 8 → 10 taken 383 times.
893 if (allocated(this%col_type)) deallocate (this%col_type)
776
2/2
✓ Branch 10 → 11 taken 510 times.
✓ Branch 10 → 12 taken 383 times.
893 if (allocated(this%col_mapping)) deallocate (this%col_mapping)
777
2/2
✓ Branch 12 → 13 taken 510 times.
✓ Branch 12 → 33 taken 383 times.
893 if (allocated(this%resp_inds)) then
778
2/2
✓ Branch 14 → 15 taken 672 times.
✓ Branch 14 → 17 taken 510 times.
1182 do s1 = 1, size(this%resp_inds)
779 1182 call this%resp_inds(s1)%destroy()
780 end do
781
8/14
✓ Branch 18 → 19 taken 510 times.
✗ Branch 18 → 30 not taken.
✓ Branch 20 → 21 taken 672 times.
✓ Branch 20 → 30 taken 510 times.
✗ Branch 21 → 22 not taken.
✓ Branch 21 → 23 taken 672 times.
✗ Branch 23 → 24 not taken.
✓ Branch 23 → 25 taken 672 times.
✗ Branch 25 → 26 not taken.
✓ Branch 25 → 27 taken 672 times.
✗ Branch 27 → 28 not taken.
✓ Branch 27 → 29 taken 672 times.
✗ Branch 30 → 31 not taken.
✓ Branch 30 → 32 taken 510 times.
1182 deallocate (this%resp_inds)
782 end if
783
784
2/2
✓ Branch 33 → 34 taken 4 times.
✓ Branch 33 → 37 taken 889 times.
893 if (this%mode == CONSTRAINT_PROJECTION) then
785
1/2
✗ Branch 35 → 36 not taken.
✓ Branch 35 → 37 taken 4 times.
4 PetscCall(MatDestroy(this%mat, ierr))
786 end if
787
1/2
✗ Branch 38 → 39 not taken.
✓ Branch 38 → 40 taken 893 times.
893 PetscCall(MatDestroy(this%mat_perp, ierr))
788 ! PetscCall(MatDestroy(this%mass_mat, ierr))
789
790
2/2
✓ Branch 40 → 41 taken 4465 times.
✓ Branch 40 → 43 taken 893 times.
5358 do s1 = 1, MAX_NR_CONSTRAINTS
791 5358 call this%constraints(s1)%destroy()
792 end do
793 893 this%nr_constraints = 0
794 end subroutine destroy_mform_constraint
795 end module m_mform_constraint_global
796