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