mform/m_mform_constraint_polar.F90
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | !> @brief Module containing implementations of various local m-form constraints, such as: | ||
| 2 | !> - MFormPolarRegularity, which imposes a polar regularity constraint on an m-form | ||
| 3 | module m_mform_constraint_polar | ||
| 4 | use m_bspline, only: BSplineFun | ||
| 5 | use m_common, only: wp | ||
| 6 | use m_mform_basis, only: MFormSpace | ||
| 7 | use m_sparsemat, only: SparseMat | ||
| 8 | use m_tensorprod_indices, only: TensorProdIndexList, size, size_2d | ||
| 9 | use m_tensorprod_shared, only: SharedMemoryWindow | ||
| 10 | use m_mform_constraint_abstract, only: MFormConstraintLocal | ||
| 11 | #include "petsc.fi" | ||
| 12 | #include "timer.fi" | ||
| 13 | |||
| 14 | implicit none | ||
| 15 | |||
| 16 | private | ||
| 17 | public :: MFormPolarRegularity, MFormPolarBoundedness | ||
| 18 | |||
| 19 | !> The threshold under which extraction/projection coefficients are considered zero | ||
| 20 | real(wp), parameter :: COEFF_THR = 1e-14_wp | ||
| 21 | |||
| 22 | !> An object that represents a polar regularity constraint on an m-form space | ||
| 23 | type, extends(MFormConstraintLocal) :: MFormPolarRegularity | ||
| 24 | !> The regularity of the polar constraint | ||
| 25 | integer :: regularity | ||
| 26 | !> Enriched space | ||
| 27 | logical :: enriched | ||
| 28 | !> Orthogonal extraction (Euclidean inner product) | ||
| 29 | logical :: orthogonal | ||
| 30 | !> The radial B-splines for the polar regularity constraint | ||
| 31 | type(BSplineFun), allocatable :: radial_bsplines(:) | ||
| 32 | !> The angular B-splines for the polar regularity constraint | ||
| 33 | type(BSplineFun), allocatable :: angular_bsplines(:) | ||
| 34 | |||
| 35 | integer, private :: row0_2d(4) | ||
| 36 | |||
| 37 | !> The zero-form space | ||
| 38 | type(MFormSpace) :: space0 | ||
| 39 | |||
| 40 | !> Communicator for the shared memory window | ||
| 41 | integer :: comm_window = -1, my_rank_window = -1 | ||
| 42 | |||
| 43 | !> The dense 2D extractor operators (shared memory) | ||
| 44 | real(wp), pointer :: E1(:, :, :, :), E2(:, :, :, :) ! i, j, block, polar_index (0, 0, 1, 1)-based indexing | ||
| 45 | |||
| 46 | !> The shared memory windows for the extractor matrices | ||
| 47 | type(SharedMemoryWindow) :: shmem_E1, shmem_E2 | ||
| 48 | contains | ||
| 49 | procedure :: init => init_mform_polar | ||
| 50 | procedure :: apply => apply_mform_polar | ||
| 51 | procedure :: destroy_derived => destroy_mform_polarregularity | ||
| 52 | procedure, private :: init_shmem_window | ||
| 53 | procedure, private :: compute_1d_bsplines | ||
| 54 | procedure, private :: insert_polar_extractor_0forms | ||
| 55 | procedure, private :: insert_polar_extractor_1forms | ||
| 56 | procedure, private :: insert_polar_extractor_1forms_enrich | ||
| 57 | procedure, private :: insert_polar_extractor_2forms | ||
| 58 | procedure, private :: insert_polar_extractor_2forms_enrich | ||
| 59 | procedure, private :: insert_polar_extractor_3forms | ||
| 60 | procedure, private :: insert_polar_extractor_3forms_enrich | ||
| 61 | procedure, private :: finalize_polar_extractor_0forms | ||
| 62 | procedure, private :: finalize_polar_extractor_1forms | ||
| 63 | procedure, private :: finalize_polar_extractor_2forms | ||
| 64 | procedure, private :: finalize_polar_extractor_3forms | ||
| 65 | end type | ||
| 66 | |||
| 67 | interface MFormPolarRegularity | ||
| 68 | module procedure construct_mform_polar | ||
| 69 | end interface | ||
| 70 | |||
| 71 | !> An alias for MFormPolarRegularity(0), i.e., a boundedness constraint on an m-form resulting in conforming FEM spaces | ||
| 72 | interface MFormPolarBoundedness | ||
| 73 | module procedure construct_mform_bounded | ||
| 74 | end interface | ||
| 75 | contains | ||
| 76 | |||
| 77 | !> @brief Construct a boundedness constraint on an m-form | ||
| 78 | !> | ||
| 79 | !> @return The constructed MFormPolarRegularity object | ||
| 80 | 40 | pure type(MFormPolarRegularity) function construct_mform_bounded() result(this) | |
| 81 | implicit none | ||
| 82 | |||
| 83 | this = MFormPolarRegularity(0) | ||
| 84 | 40 | end function construct_mform_bounded | |
| 85 | |||
| 86 | !> @brief Construct a polar regularity constraint on an m-form | ||
| 87 | !> | ||
| 88 | !> @param[in] _(optional)_ regularity The regularity of the polar constraint (default: degree of the B-spline basis in the | ||
| 89 | !> x-direction) | ||
| 90 | !> @param[in] _(optional)_ enriched Whether to use an enriched polar regularity constraint (default: .true.) | ||
| 91 | !> @param[in] _(optional)_ orthogonal Whether to use orthogonal extraction operators (default: .true.) | ||
| 92 | !> | ||
| 93 | !> @return The constructed MFormPolarRegularity object | ||
| 94 | 804 | pure type(MFormPolarRegularity) function construct_mform_polar(regularity, enriched, orthogonal) result(this) | |
| 95 | implicit none | ||
| 96 | |||
| 97 | integer, intent(in), optional :: regularity | ||
| 98 | logical, intent(in), optional :: enriched | ||
| 99 | logical, intent(in), optional :: orthogonal | ||
| 100 | |||
| 101 |
2/2✓ Branch 2 → 3 taken 692 times.
✓ Branch 2 → 4 taken 152 times.
|
844 | if (present(regularity)) then |
| 102 | 652 | this%regularity = regularity | |
| 103 | else | ||
| 104 | this%regularity = -1 ! will be set properly in init method based on degree in x-direction | ||
| 105 | end if | ||
| 106 | |||
| 107 |
2/2✓ Branch 4 → 5 taken 700 times.
✓ Branch 4 → 6 taken 104 times.
|
804 | if (present(enriched)) then |
| 108 | 700 | this%enriched = enriched | |
| 109 | else | ||
| 110 | this%enriched = .true. | ||
| 111 | end if | ||
| 112 | |||
| 113 |
2/2✓ Branch 6 → 7 taken 280 times.
✓ Branch 6 → 8 taken 524 times.
|
804 | if (present(orthogonal)) then |
| 114 | 280 | this%orthogonal = orthogonal | |
| 115 | else | ||
| 116 | this%orthogonal = .true. | ||
| 117 | end if | ||
| 118 | 804 | end function construct_mform_polar | |
| 119 | |||
| 120 | !> @brief Initialize a polar regularity constraint on an m-form space | ||
| 121 | !> | ||
| 122 | !> @param[inout] this The MFormPolarRegularity object to initialize | ||
| 123 | !> @param[in] space The MFormSpace to which the polar regularity constraint applies | ||
| 124 | !> @param[in] _(optional)_ coord_transform The coordinate transformation of the m-form space | ||
| 125 | 1087 | subroutine init_mform_polar(this, space, coord_transform) | |
| 126 | use m_mform_basis, only: MFormSpace, get_zero_form_space | ||
| 127 | use m_qr_factorisation, only: orthogonalize | ||
| 128 | use m_tensorprod, only: size, TensorProdIndices | ||
| 129 | use m_coord_transform_abstract, only: CoordTransformAbstract | ||
| 130 | implicit none | ||
| 131 | |||
| 132 | class(MFormPolarRegularity), intent(inout) :: this | ||
| 133 | type(MFormSpace), intent(in) :: space | ||
| 134 | class(CoordTransformAbstract), optional, intent(in) :: coord_transform | ||
| 135 | |||
| 136 | integer :: s1, s2, nr_blocks | ||
| 137 | |||
| 138 | logical :: coord_has_polar_xy_singularity_ | ||
| 139 | |||
| 140 | type(TensorProdIndices) :: inds | ||
| 141 | |||
| 142 | TIMER_DECL(t_total) | ||
| 143 | TIMER_DECL(t_orthogonalize) | ||
| 144 | TIMER_INIT_START(t_total, "MFormPolarRegularity::init") | ||
| 145 | |||
| 146 | 1087 | call this%destroy() | |
| 147 | |||
| 148 |
5/8✓ Branch 3 → 4 taken 1087 times.
✗ Branch 3 → 7 not taken.
✓ Branch 4 → 5 taken 1087 times.
✗ Branch 4 → 6 not taken.
✓ Branch 7 → 8 taken 1087 times.
✗ Branch 7 → 10 not taken.
✓ Branch 8 → 9 taken 60 times.
✓ Branch 8 → 10 taken 1027 times.
|
1087 | this%space = space |
| 149 |
2/2✓ Branch 11 → 12 taken 60 times.
✓ Branch 11 → 13 taken 1027 times.
|
1087 | this%space0 = get_zero_form_space(space) |
| 150 | |||
| 151 | ! The constraint is shared amongst MPI ranks, however, most ranks contribute zeros to the extractor | ||
| 152 | 1087 | this%is_mpi_shared(1) = this%space%tp_spaces(1)%domain%nr_subintervals(1) > 1 | |
| 153 | 1087 | this%is_mpi_shared(2) = this%space%tp_spaces(1)%domain%nr_subintervals(2) > 1 | |
| 154 | 1087 | this%is_mpi_shared(3) = .false. | |
| 155 | |||
| 156 | coord_has_polar_xy_singularity_ = .false. | ||
| 157 |
3/4✓ Branch 13 → 14 taken 1067 times.
✓ Branch 13 → 16 taken 20 times.
✓ Branch 14 → 15 taken 1067 times.
✗ Branch 14 → 16 not taken.
|
1087 | if (present(coord_transform)) coord_has_polar_xy_singularity_ = coord_transform%has_polar_xy_singularity |
| 158 | |||
| 159 | 1087 | nr_blocks = size(space%tp_spaces) | |
| 160 |
9/16✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 1087 times.
✓ Branch 18 → 17 taken 1087 times.
✗ Branch 18 → 19 not taken.
✓ Branch 19 → 20 taken 1087 times.
✗ Branch 19 → 21 not taken.
✓ Branch 21 → 22 taken 1087 times.
✗ Branch 21 → 23 not taken.
✗ Branch 23 → 24 not taken.
✓ Branch 23 → 25 taken 1087 times.
✗ Branch 25 → 26 not taken.
✓ Branch 25 → 27 taken 1087 times.
✗ Branch 27 → 28 not taken.
✓ Branch 27 → 29 taken 1087 times.
✓ Branch 30 → 31 taken 2195 times.
✓ Branch 30 → 32 taken 1087 times.
|
5456 | allocate (this%index_list(nr_blocks)) |
| 161 |
9/16✗ Branch 32 → 33 not taken.
✓ Branch 32 → 34 taken 1087 times.
✓ Branch 34 → 33 taken 1087 times.
✗ Branch 34 → 35 not taken.
✓ Branch 35 → 36 taken 1087 times.
✗ Branch 35 → 37 not taken.
✓ Branch 37 → 38 taken 1087 times.
✗ Branch 37 → 39 not taken.
✗ Branch 39 → 40 not taken.
✓ Branch 39 → 41 taken 1087 times.
✗ Branch 41 → 42 not taken.
✓ Branch 41 → 43 taken 1087 times.
✗ Branch 43 → 44 not taken.
✓ Branch 43 → 45 taken 1087 times.
✓ Branch 46 → 47 taken 2195 times.
✓ Branch 46 → 48 taken 1087 times.
|
5456 | allocate (this%extractor(nr_blocks)) |
| 162 | |||
| 163 |
2/2✓ Branch 48 → 49 taken 20 times.
✓ Branch 48 → 58 taken 1067 times.
|
1087 | if (.not. coord_has_polar_xy_singularity_) then |
| 164 | 20 | call inds%empty() | |
| 165 | |||
| 166 |
2/2✓ Branch 51 → 52 taken 20 times.
✓ Branch 51 → 56 taken 20 times.
|
40 | do s1 = 1, nr_blocks |
| 167 | 20 | call this%index_list(s1)%init(inds, space%tp_spaces(s1)%rank_resp_bspline, space%tp_spaces(s1)%rank_l0) | |
| 168 | 20 | call this%extractor(s1)%init(nr_rows=0, nr_cols=0, nr_nonzero=0) | |
| 169 | 40 | call this%extractor(s1)%finalize() | |
| 170 | end do | ||
| 171 | |||
| 172 | TIMER_STOP(t_total) | ||
| 173 | 20 | return | |
| 174 | end if | ||
| 175 | |||
| 176 |
1/4✗ Branch 58 → 59 not taken.
✓ Branch 58 → 65 taken 1067 times.
✗ Branch 59 → 60 not taken.
✗ Branch 59 → 65 not taken.
|
1067 | if (space%tp_spaces(1)%domain%is_distributed_memory(1) .and. space%tp_spaces(1)%domain%my_subinterval_ijk(1) > 0) then |
| 177 | ✗ | do s1 = 1, size(space%tp_spaces) | |
| 178 | ✗ | if (space%tp_spaces(s1)%rank_resp_bspline%i1 < this%regularity) then | |
| 179 | error stop "MFormPolarRegularity::init_mform_polar: the first subinterval in the x-direction of the tensor product"// & | ||
| 180 | ✗ | & " space must have at least 'regularity + 1' B-splines whenever this direction is distributed in memory." | |
| 181 | end if | ||
| 182 | end do | ||
| 183 | end if | ||
| 184 |
1/2✗ Branch 65 → 66 not taken.
✓ Branch 65 → 67 taken 1067 times.
|
1067 | if (space%tp_spaces(1)%domain%is_distributed_memory(2)) then |
| 185 | error stop "MFormPolarRegularity::init_mform_polar: The polar regularity constraint is not compatible with a distributed "// & | ||
| 186 | ✗ | & "memory domain decomposition in the y-direction" | |
| 187 | end if | ||
| 188 | |||
| 189 |
2/4✓ Branch 67 → 68 taken 1067 times.
✗ Branch 67 → 69 not taken.
✗ Branch 68 → 69 not taken.
✓ Branch 68 → 70 taken 1067 times.
|
1067 | if (space%tp_spaces(1)%spaces(1)%is_periodic .or. .not. space%tp_spaces(1)%spaces(2)%is_periodic) then |
| 190 | error stop "MFormPolarRegularity::init_mform_polar: The polar regularity constraint requires the B-spline basis in the"// & | ||
| 191 | ✗ | & " x-direction to be non-periodic and the B-spline basis in the y-direction to be periodic" | |
| 192 | end if | ||
| 193 | |||
| 194 | 1067 | this%is_homogeneous = .true. | |
| 195 | 1067 | this%is_symmetric = .false. | |
| 196 | |||
| 197 |
2/2✓ Branch 70 → 71 taken 160 times.
✓ Branch 70 → 72 taken 907 times.
|
1067 | if (this%regularity <= -1) then |
| 198 | 160 | this%regularity = this%space0%tp_spaces(1)%spaces(1)%degree | |
| 199 | else | ||
| 200 |
1/2✗ Branch 72 → 73 not taken.
✓ Branch 72 → 74 taken 907 times.
|
907 | if (this%regularity < 0 .or. this%regularity > this%space0%tp_spaces(1)%spaces(1)%degree) then |
| 201 | error stop "MFormPolarRegularity::init_mform_polar: 'regularity' must be in [0, degree_x] where 'degree_x' is the "// & | ||
| 202 | ✗ | & "degree of the B-spline basis in the x-direction of the zero-form space" | |
| 203 | end if | ||
| 204 | end if | ||
| 205 | |||
| 206 | ! The forward part of commutativty (strong part of Hodge decomposition) | ||
| 207 | 1470 | select case (space%m) | |
| 208 | case (0) | ||
| 209 | 403 | call this%insert_polar_extractor_0forms() | |
| 210 | case (1) | ||
| 211 | 333 | call this%insert_polar_extractor_1forms() | |
| 212 | case (2) | ||
| 213 | 221 | call this%insert_polar_extractor_2forms() | |
| 214 | case (3) | ||
| 215 |
4/5✓ Branch 74 → 75 taken 403 times.
✓ Branch 74 → 76 taken 333 times.
✓ Branch 74 → 77 taken 221 times.
✓ Branch 74 → 78 taken 110 times.
✗ Branch 74 → 79 not taken.
|
1067 | call this%insert_polar_extractor_3forms() |
| 216 | end select | ||
| 217 | |||
| 218 |
4/4✓ Branch 79 → 80 taken 731 times.
✓ Branch 79 → 85 taken 336 times.
✓ Branch 80 → 81 taken 420 times.
✓ Branch 80 → 85 taken 311 times.
|
1067 | if (this%enriched .and. this%regularity > 0) then |
| 219 | ! The backward part of commutativity (weak part of Hodge decomposition) | ||
| 220 | 106 | select case (space%m) | |
| 221 | case (1) | ||
| 222 | 106 | call this%insert_polar_extractor_1forms_enrich() | |
| 223 | case (2) | ||
| 224 | 94 | call this%insert_polar_extractor_2forms_enrich() | |
| 225 | case (3) | ||
| 226 |
4/4✓ Branch 81 → 82 taken 106 times.
✓ Branch 81 → 83 taken 94 times.
✓ Branch 81 → 84 taken 49 times.
✓ Branch 81 → 85 taken 171 times.
|
420 | call this%insert_polar_extractor_3forms_enrich() |
| 227 | end select | ||
| 228 | end if | ||
| 229 | |||
| 230 |
1/2✓ Branch 85 → 86 taken 1067 times.
✗ Branch 85 → 106 not taken.
|
1067 | if (this%orthogonal) then |
| 231 | TIMER_INIT_START(t_orthogonalize, "MFormPolarRegularity::orthogonalize_extractors") | ||
| 232 | ! The 3D polar forms are Kronecker products of identity matrix (z) with 2D polar forms, so 2D orthogonalization is sufficient | ||
| 233 | ! NOTE: E1 and E2 are assumed to act on diffent components of the m-form, hence orthogonalization is independent | ||
| 234 |
1/2✓ Branch 86 → 87 taken 1067 times.
✗ Branch 86 → 96 not taken.
|
1067 | if (associated(this%E1)) then |
| 235 | 1067 | call this%shmem_E1%fence() | |
| 236 |
2/4✓ Branch 88 → 89 taken 1067 times.
✗ Branch 88 → 94 not taken.
✗ Branch 91 → 92 not taken.
✓ Branch 91 → 94 taken 1067 times.
|
1067 | if (this%shmem_E1%leader()) call orthogonalize(this%E1) |
| 237 | 1067 | call this%shmem_E1%fence() | |
| 238 | end if | ||
| 239 |
2/2✓ Branch 96 → 97 taken 337 times.
✓ Branch 96 → 106 taken 730 times.
|
1067 | if (associated(this%E2)) then |
| 240 | 337 | call this%shmem_E2%fence() | |
| 241 |
2/4✓ Branch 98 → 99 taken 337 times.
✗ Branch 98 → 104 not taken.
✗ Branch 101 → 102 not taken.
✓ Branch 101 → 104 taken 337 times.
|
337 | if (this%shmem_E2%leader()) call orthogonalize(this%E2) |
| 242 | 337 | call this%shmem_E2%fence() | |
| 243 | end if | ||
| 244 | TIMER_STOP(t_orthogonalize) | ||
| 245 | end if | ||
| 246 | |||
| 247 | 403 | select case (space%m) | |
| 248 | case (0) | ||
| 249 | 403 | call this%finalize_polar_extractor_0forms() | |
| 250 | case (1) | ||
| 251 | 333 | call this%finalize_polar_extractor_1forms() | |
| 252 | case (2) | ||
| 253 | 221 | call this%finalize_polar_extractor_2forms() | |
| 254 | case (3) | ||
| 255 |
4/5✓ Branch 106 → 107 taken 403 times.
✓ Branch 106 → 108 taken 333 times.
✓ Branch 106 → 109 taken 221 times.
✓ Branch 106 → 110 taken 110 times.
✗ Branch 106 → 111 not taken.
|
1067 | call this%finalize_polar_extractor_3forms() |
| 256 | end select | ||
| 257 | |||
| 258 |
1/2✓ Branch 111 → 112 taken 1067 times.
✗ Branch 111 → 114 not taken.
|
1067 | if (associated(this%E1)) then |
| 259 | 1067 | call this%shmem_E1%destroy() | |
| 260 | 1067 | nullify (this%E1) | |
| 261 | end if | ||
| 262 |
2/2✓ Branch 114 → 115 taken 337 times.
✓ Branch 114 → 117 taken 730 times.
|
1067 | if (associated(this%E2)) then |
| 263 | 337 | call this%shmem_E2%destroy() | |
| 264 | 337 | nullify (this%E2) | |
| 265 | end if | ||
| 266 | |||
| 267 | TIMER_STOP(t_total) | ||
| 268 | end subroutine init_mform_polar | ||
| 269 | |||
| 270 | !> @brief Initialize the shared memory windows for the extractor matrices | ||
| 271 | !> | ||
| 272 | !> @param[inout] this The MFormPolarRegularity object to initialize | ||
| 273 | !> @param[inout] E4d The 4D extractor matrix to be initialized as a shared memory window | ||
| 274 | !> @param[inout] E_window The SharedMemoryWindow object managing the shared memory for E4d | ||
| 275 | !> @param[in] nr_polar_splines The number of polar splines used in the extractor | ||
| 276 | 1568 | subroutine init_shmem_window(this, E4d, E_window, nr_polar_splines) | |
| 277 | use m_common, only: cartesian_to_linear | ||
| 278 | use m_tensorprod, only: size | ||
| 279 | |||
| 280 | implicit none | ||
| 281 | |||
| 282 | class(MFormPolarRegularity), intent(inout) :: this | ||
| 283 | real(wp), pointer :: E4d(:, :, :, :) | ||
| 284 | type(SharedMemoryWindow), intent(inout) :: E_window | ||
| 285 | integer, intent(in) :: nr_polar_splines | ||
| 286 | |||
| 287 | integer :: comm_color, comm_key | ||
| 288 | integer :: nr_x, nr_y, nr_blocks, nr_elements | ||
| 289 | 1568 | integer, allocatable :: cart_inds(:), cart_size(:), non_cart_inds(:), non_cart_size(:) | |
| 290 | integer :: ierr | ||
| 291 | logical :: mask(3) | ||
| 292 | |||
| 293 | TIMER_DECL(t_total) | ||
| 294 | TIMER_INIT_START(t_total, "MFormPolarRegularity::init_shmem_window") | ||
| 295 | |||
| 296 |
2/2✓ Branch 2 → 3 taken 1067 times.
✓ Branch 2 → 19 taken 501 times.
|
1568 | if (this%comm_window == -1) then |
| 297 | ! Create a shared memory sub communicator for the ranks involved in the constraint that share mpi unknowns | ||
| 298 | 1067 | mask(1) = this%space%tp_spaces(1)%domain%is_shared_memory(1) | |
| 299 | 1067 | mask(2) = this%space%tp_spaces(1)%domain%is_shared_memory(2) | |
| 300 | 1067 | mask(3) = .false. | |
| 301 | |||
| 302 | 1067 | cart_inds = pack(this%space%tp_spaces(1)%domain%my_subinterval_ijk, mask=mask) | |
| 303 | 1067 | cart_size = pack(this%space%tp_spaces(1)%domain%nr_subintervals, mask=mask) | |
| 304 |
2/2✓ Branch 6 → 7 taken 3201 times.
✓ Branch 6 → 8 taken 1067 times.
|
4268 | non_cart_inds = pack(this%space%tp_spaces(1)%domain%my_subinterval_ijk, mask=.not. mask) |
| 305 |
2/2✓ Branch 10 → 11 taken 3201 times.
✓ Branch 10 → 12 taken 1067 times.
|
4268 | non_cart_size = pack(this%space%tp_spaces(1)%domain%nr_subintervals, mask=.not. mask) |
| 306 | |||
| 307 | 1067 | comm_color = cartesian_to_linear(non_cart_inds, non_cart_size) ! Same color ends up on same comm | |
| 308 | 1067 | comm_key = cartesian_to_linear(cart_inds, cart_size) ! Key to order ranks within comm | |
| 309 | |||
| 310 |
1/2✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 1067 times.
|
1067 | PetscCallMPI(MPI_Comm_split(this%space%tp_spaces(1)%domain%comm, comm_color, comm_key, this%comm_window, ierr)) |
| 311 |
1/2✗ Branch 17 → 18 not taken.
✓ Branch 17 → 19 taken 1067 times.
|
1067 | PetscCallMPI(MPI_Comm_rank(this%comm_window, this%my_rank_window, ierr)) |
| 312 | end if | ||
| 313 | |||
| 314 | nr_x = merge(this%regularity + 1, 0, this%space%tp_spaces(1)%domain%is_shared_memory(1) .or. & | ||
| 315 |
2/4✗ Branch 19 → 20 not taken.
✓ Branch 19 → 21 taken 1568 times.
✓ Branch 21 → 20 taken 1568 times.
✗ Branch 21 → 22 not taken.
|
1568 | this%space%tp_spaces(1)%domain%my_subinterval_ijk(1) == 0) |
| 316 | 1568 | nr_y = size(this%space%tp_spaces(1), 2) ! All spaces have same size in y-direction | |
| 317 | 1568 | nr_blocks = size(this%space%tp_spaces) | |
| 318 | |||
| 319 | 1568 | nr_elements = nr_x * nr_y * nr_blocks * nr_polar_splines | |
| 320 | |||
| 321 | 1568 | call E_window%init(nr_elements, this%comm_window, this%my_rank_window) | |
| 322 | 1568 | E4d(0:nr_x - 1, 0:nr_y - 1, 1:nr_blocks, 1:nr_polar_splines) => E_window%window(:) | |
| 323 |
9/10✓ Branch 23 → 24 taken 1568 times.
✗ Branch 23 → 35 not taken.
✓ Branch 24 → 25 taken 4869 times.
✓ Branch 24 → 35 taken 1568 times.
✓ Branch 26 → 27 taken 10309 times.
✓ Branch 26 → 34 taken 4869 times.
✓ Branch 28 → 29 taken 291314 times.
✓ Branch 28 → 33 taken 10309 times.
✓ Branch 30 → 31 taken 1159546 times.
✓ Branch 30 → 32 taken 291314 times.
|
1467606 | if (E_window%leader()) E4d = 0._wp |
| 324 | 1568 | call E_window%fence() | |
| 325 | |||
| 326 | TIMER_STOP(t_total) | ||
| 327 |
8/16✓ Branch 36 → 37 taken 1067 times.
✓ Branch 36 → 38 taken 501 times.
✓ Branch 38 → 39 taken 1067 times.
✓ Branch 38 → 40 taken 501 times.
✓ Branch 40 → 41 taken 1067 times.
✓ Branch 40 → 42 taken 501 times.
✓ Branch 42 → 43 taken 1067 times.
✓ Branch 42 → 59 taken 501 times.
✗ Branch 44 → 45 not taken.
✗ Branch 44 → 46 not taken.
✗ Branch 48 → 49 not taken.
✗ Branch 48 → 50 not taken.
✗ Branch 52 → 53 not taken.
✗ Branch 52 → 54 not taken.
✗ Branch 56 → 57 not taken.
✗ Branch 56 → 60 not taken.
|
1568 | end subroutine init_shmem_window |
| 328 | |||
| 329 | !> @brief Apply the polar regularity constraint to an m-form (overload of MFormConstraintLocal::apply) | ||
| 330 | !> | ||
| 331 | !> @param[in] this The constraint to apply | ||
| 332 | !> @param[out] vout The MForm to which the constraint has been applied | ||
| 333 | !> @param[in] vin The MForm to which the constraint is applied | ||
| 334 |
4/4✓ Branch 2 → 3 taken 388 times.
✓ Branch 2 → 4 taken 253 times.
✓ Branch 4 → 5 taken 388 times.
✓ Branch 4 → 6 taken 253 times.
|
641 | subroutine apply_mform_polar(this, vout, vin) |
| 335 | use m_mform_basis, only: MFormFun | ||
| 336 | implicit none | ||
| 337 | |||
| 338 | class(MFormPolarRegularity), intent(in) :: this | ||
| 339 | type(MFormFun), intent(out) :: vout | ||
| 340 | type(MFormFun), intent(in) :: vin | ||
| 341 | |||
| 342 | integer :: blk, local_idx, local_row, local_col | ||
| 343 | integer :: i, j, k, row_start, row_end | ||
| 344 | real(wp) :: local_val | ||
| 345 | real(wp), allocatable :: reduced_coeffs(:), my_reduced_coeffs(:) | ||
| 346 | integer :: nr_cols, ierr | ||
| 347 | |||
| 348 | TIMER_DECL(t_total) | ||
| 349 | TIMER_DECL(t_mpi_reduce) | ||
| 350 | TIMER_INIT_START(t_total, "MFormPolarRegularity::apply") | ||
| 351 | |||
| 352 | ! The columns of the extraction operator are orthonormal, so we can use the following projection | ||
| 353 | !> vout = E * (E^T * vin) | ||
| 354 | !> | ||
| 355 | !> This is done in two steps: | ||
| 356 | !> 1. Compute reduced_coeffs = E^T * vin (project vin onto the constrained subspace) | ||
| 357 | !> 2. Compute vout = E * reduced_coeffs (reconstruct in the full space) | ||
| 358 | |||
| 359 | ! Make a deep copy | ||
| 360 | 641 | vout = vin | |
| 361 | |||
| 362 |
1/2✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 641 times.
|
641 | if (.not. allocated(this%extractor)) error stop "MFormPolarRegularity::apply_mform_polar: extractor not allocated" |
| 363 |
1/2✗ Branch 9 → 10 not taken.
✓ Branch 9 → 11 taken 641 times.
|
641 | if (.not. this%orthogonal) error stop "MFormPolarRegularity::apply_mform_polar: only orthogonal extraction can be used in apply" |
| 364 |
1/2✗ Branch 11 → 12 not taken.
✓ Branch 11 → 13 taken 641 times.
|
641 | if (this%comm_window == -1) error stop "MFormPolarRegularity::apply_mform_polar: shared memory window communicator not"// & |
| 365 | ✗ | " initialized" | |
| 366 | |||
| 367 | ! Apply projection using extraction operator: vout = E * (E^T * vin) | ||
| 368 | ! Only apply on the rank at the polar singularity (my_rank_i == 0) | ||
| 369 | ! TODO node-level-mpi: requires shared memory window with sparse extractor matrices | ||
| 370 | ! Set to zero where the constraint acts | ||
| 371 |
2/2✓ Branch 14 → 15 taken 1587 times.
✓ Branch 14 → 20 taken 641 times.
|
2228 | do blk = 1, size(vout%tp_funs) |
| 372 |
2/2✓ Branch 16 → 17 taken 142812 times.
✓ Branch 16 → 18 taken 1587 times.
|
145040 | do local_idx = 0, size(this%index_list(blk)) - 1 |
| 373 | 142812 | i = this%index_list(blk)%i(local_idx) | |
| 374 | 142812 | j = this%index_list(blk)%j(local_idx) | |
| 375 | 142812 | k = this%index_list(blk)%k(local_idx) | |
| 376 | |||
| 377 | 144399 | vout%tp_funs(blk)%data(i, j, k) = 0._wp | |
| 378 | end do | ||
| 379 | end do | ||
| 380 | |||
| 381 | 641 | nr_cols = this%extractor(1)%nr_cols | |
| 382 |
8/12✓ Branch 21 → 22 taken 369 times.
✓ Branch 21 → 23 taken 272 times.
✓ Branch 23 → 22 taken 272 times.
✗ Branch 23 → 24 not taken.
✓ Branch 24 → 25 taken 641 times.
✗ Branch 24 → 26 not taken.
✓ Branch 26 → 27 taken 272 times.
✓ Branch 26 → 28 taken 369 times.
✗ Branch 28 → 29 not taken.
✓ Branch 28 → 30 taken 641 times.
✗ Branch 30 → 31 not taken.
✓ Branch 30 → 32 taken 641 times.
|
1923 | allocate (my_reduced_coeffs(0:nr_cols - 1)) |
| 383 |
2/2✓ Branch 32 → 33 taken 4884 times.
✓ Branch 32 → 34 taken 641 times.
|
5525 | my_reduced_coeffs = 0._wp |
| 384 | |||
| 385 | ! Step 1: Compute E^T * vin (transpose multiply) | ||
| 386 | ! Loop over all input blocks that contribute | ||
| 387 |
2/2✓ Branch 35 → 36 taken 1587 times.
✓ Branch 35 → 44 taken 641 times.
|
2228 | do blk = 1, size(vin%tp_funs) |
| 388 |
2/2✓ Branch 37 → 38 taken 142812 times.
✓ Branch 37 → 42 taken 1587 times.
|
145040 | do local_row = 0, this%extractor(blk)%nr_rows - 1 |
| 389 | 142812 | row_start = this%extractor(blk)%row_starts_at_nz(local_row) | |
| 390 | 142812 | row_end = this%extractor(blk)%row_starts_at_nz(local_row + 1) - 1 | |
| 391 | |||
| 392 | 142812 | i = this%index_list(blk)%i(local_row) | |
| 393 | 142812 | j = this%index_list(blk)%j(local_row) | |
| 394 | 142812 | k = this%index_list(blk)%k(local_row) | |
| 395 | |||
| 396 |
2/2✓ Branch 39 → 40 taken 142812 times.
✓ Branch 39 → 41 taken 453180 times.
|
597579 | do local_idx = row_start, row_end |
| 397 | 453180 | local_col = this%extractor(blk)%col_idx(local_idx) | |
| 398 | 453180 | local_val = this%extractor(blk)%values(local_idx) | |
| 399 | |||
| 400 | ! E^T multiply: reduced_coeffs(col) += E(col, row) * vin(row) | ||
| 401 | 595992 | my_reduced_coeffs(local_col) = my_reduced_coeffs(local_col) + local_val * vin%tp_funs(blk)%data(i, j, k) | |
| 402 | end do | ||
| 403 | end do | ||
| 404 | end do | ||
| 405 | |||
| 406 | TIMER_INIT_START(t_mpi_reduce, "MPI_Allreduce") | ||
| 407 |
5/8✓ Branch 45 → 46 taken 369 times.
✓ Branch 45 → 47 taken 272 times.
✓ Branch 47 → 46 taken 272 times.
✗ Branch 47 → 48 not taken.
✗ Branch 48 → 49 not taken.
✓ Branch 48 → 50 taken 641 times.
✗ Branch 50 → 51 not taken.
✓ Branch 50 → 52 taken 641 times.
|
1282 | allocate (reduced_coeffs(0:nr_cols - 1)) |
| 408 | ! TODO shared memory equivalent? yes, use a fenced shmem window | ||
| 409 |
1/2✗ Branch 53 → 54 not taken.
✓ Branch 53 → 57 taken 641 times.
|
641 | PetscCallMPI(MPI_Allreduce(my_reduced_coeffs, reduced_coeffs, nr_cols, MPI_WP, MPI_SUM, this%comm_window, ierr)) |
| 410 | TIMER_STOP(t_mpi_reduce) | ||
| 411 | |||
| 412 | ! Step 2: Compute E * reduced_coeffs (forward multiply) | ||
| 413 |
2/2✓ Branch 58 → 59 taken 1587 times.
✓ Branch 58 → 67 taken 641 times.
|
2228 | do blk = 1, size(vout%tp_funs) |
| 414 |
2/2✓ Branch 60 → 61 taken 142812 times.
✓ Branch 60 → 65 taken 1587 times.
|
145040 | do local_row = 0, this%extractor(blk)%nr_rows - 1 |
| 415 | 142812 | row_start = this%extractor(blk)%row_starts_at_nz(local_row) | |
| 416 | 142812 | row_end = this%extractor(blk)%row_starts_at_nz(local_row + 1) - 1 | |
| 417 | |||
| 418 | 142812 | i = this%index_list(blk)%i(local_row) | |
| 419 | 142812 | j = this%index_list(blk)%j(local_row) | |
| 420 | 142812 | k = this%index_list(blk)%k(local_row) | |
| 421 | |||
| 422 |
2/2✓ Branch 62 → 63 taken 142812 times.
✓ Branch 62 → 64 taken 453180 times.
|
597579 | do local_idx = row_start, row_end |
| 423 | 453180 | local_col = this%extractor(blk)%col_idx(local_idx) | |
| 424 | 453180 | local_val = this%extractor(blk)%values(local_idx) | |
| 425 | |||
| 426 | ! E multiply: vout(row) += E(row, col) * reduced_coeffs(col) | ||
| 427 | vout%tp_funs(blk)%data(i, j, k) = & | ||
| 428 | 595992 | vout%tp_funs(blk)%data(i, j, k) + local_val * reduced_coeffs(local_col) | |
| 429 | end do | ||
| 430 | end do | ||
| 431 | end do | ||
| 432 | |||
| 433 |
1/2✗ Branch 68 → 69 not taken.
✓ Branch 68 → 70 taken 641 times.
|
641 | deallocate (reduced_coeffs) |
| 434 |
1/2✗ Branch 70 → 71 not taken.
✓ Branch 70 → 72 taken 641 times.
|
641 | deallocate (my_reduced_coeffs) |
| 435 | |||
| 436 | 641 | call vout%distribute() | |
| 437 | |||
| 438 | TIMER_STOP(t_total) | ||
| 439 | ✗ | end subroutine apply_mform_polar | |
| 440 | |||
| 441 | !> @brief Initialize the 0-form extractor for the polar regularity constraint | ||
| 442 | !> | ||
| 443 | !> @param[inout] this The MFormPolarRegularity object to initialize | ||
| 444 | 403 | subroutine insert_polar_extractor_0forms(this) | |
| 445 | use m_tensorprod, only: size | ||
| 446 | use m_tensorprod_indices, only: TensorProdIndices | ||
| 447 | implicit none | ||
| 448 | |||
| 449 | class(MFormPolarRegularity), intent(inout) :: this | ||
| 450 | |||
| 451 | type(TensorProdIndices) :: inds | ||
| 452 | |||
| 453 | integer :: i, j, i_min, i_max, j_min, j_max, k_min, k_max | ||
| 454 | integer :: gamma, alpha, polar_index, nr_polar_splines, gamma_min, gamma_max | ||
| 455 | integer :: s1 | ||
| 456 | real(wp) :: coeff | ||
| 457 | |||
| 458 | TIMER_DECL(t_total) | ||
| 459 | TIMER_INIT_START(t_total, "MFormPolarRegularity::insert0") | ||
| 460 | |||
| 461 | 403 | gamma_min = 0 | |
| 462 | 403 | gamma_max = this%regularity | |
| 463 | 403 | nr_polar_splines = get_nr_polar_zeroforms(gamma_max) | |
| 464 | |||
| 465 | !> Initialize the indices on which this extractor acts | ||
| 466 | 403 | i_min = max(0, this%space%tp_spaces(1)%rank_resp_bspline%i0) | |
| 467 | 403 | i_max = min(this%regularity, this%space%tp_spaces(1)%rank_resp_bspline%i1) | |
| 468 | 403 | j_min = this%space%tp_spaces(1)%rank_resp_bspline%j0; j_max = this%space%tp_spaces(1)%rank_resp_bspline%j1 | |
| 469 | 403 | k_min = this%space%tp_spaces(1)%rank_resp_bspline%k0; k_max = this%space%tp_spaces(1)%rank_resp_bspline%k1 | |
| 470 | 403 | call inds%init(i_min, i_max, j_min, j_max, k_min, k_max, this%space%tp_spaces(1)%spaces) | |
| 471 | 403 | call this%index_list(1)%init(inds, this%space%tp_spaces(1)%rank_resp_bspline, this%space%tp_spaces(1)%rank_l0) | |
| 472 |
2/2✓ Branch 5 → 6 taken 1612 times.
✓ Branch 5 → 7 taken 403 times.
|
2015 | this%row0_2d = 0 |
| 473 | |||
| 474 | !> Initialize the radial and angular B-splines | ||
| 475 | call this%compute_1d_bsplines(this%space%tp_spaces(1)%spaces(1), this%space%tp_spaces(1)%spaces(2), gamma_min, gamma_max, & | ||
| 476 | 403 | radial_truncate=this%regularity) | |
| 477 | |||
| 478 | 403 | call this%init_shmem_window(this%E1, this%shmem_E1, nr_polar_splines) | |
| 479 | s1 = 1 | ||
| 480 |
2/2✓ Branch 9 → 10 taken 11388 times.
✓ Branch 9 → 25 taken 403 times.
|
11791 | do j = j_min, j_max |
| 481 |
2/2✓ Branch 11 → 12 taken 26746 times.
✓ Branch 11 → 23 taken 11388 times.
|
38537 | do i = i_min, i_max |
| 482 | polar_index = 0 | ||
| 483 |
2/2✓ Branch 13 → 14 taken 86938 times.
✓ Branch 13 → 21 taken 26746 times.
|
125072 | do gamma = gamma_min, gamma_max |
| 484 |
1/2✓ Branch 14 → 15 taken 86938 times.
✗ Branch 14 → 19 not taken.
|
113684 | do alpha = -gamma, gamma, 2 |
| 485 |
2/2✓ Branch 16 → 17 taken 86938 times.
✓ Branch 16 → 18 taken 127824 times.
|
214762 | this%E1(i, j, s1, polar_index + 1) = this%radial_bsplines(gamma)%data(i) * this%angular_bsplines(alpha)%data(j) |
| 486 | 86938 | polar_index = polar_index + 1 | |
| 487 | end do | ||
| 488 | end do | ||
| 489 | end do | ||
| 490 | end do | ||
| 491 | |||
| 492 | TIMER_STOP(t_total) | ||
| 493 | 403 | end subroutine insert_polar_extractor_0forms | |
| 494 | |||
| 495 | !> @brief Finalize the 0-form extractor for the polar regularity constraint | ||
| 496 | !> | ||
| 497 | !> @param[inout] this The MFormPolarRegularity object to initialize | ||
| 498 | 403 | subroutine finalize_polar_extractor_0forms(this) | |
| 499 | implicit none | ||
| 500 | |||
| 501 | class(MFormPolarRegularity), intent(inout) :: this | ||
| 502 | |||
| 503 | integer :: nr_nonzeros | ||
| 504 | integer :: i, j, k, row_3d, col, polar_index, s1 | ||
| 505 | integer :: i_min, i_max, j_min, j_max, k_min, k_max, nr_polar_splines | ||
| 506 | real(wp) :: coeff | ||
| 507 | |||
| 508 | TIMER_DECL(t_total) | ||
| 509 | TIMER_INIT_START(t_total, "MFormPolarRegularity::finalize0") | ||
| 510 | |||
| 511 | 403 | call this%index_list(1)%get_bounds(i_min, i_max, j_min, j_max, k_min, k_max) | |
| 512 | |||
| 513 | 403 | nr_polar_splines = size(this%E1, 4) | |
| 514 | |||
| 515 | !> Insert the tensor product of radial and angular B-splines into the extractor | ||
| 516 | 403 | nr_nonzeros = size(this%index_list(1)) * nr_polar_splines | |
| 517 | call this%extractor(1)%init(nr_rows=size(this%index_list(1)), nr_cols=nr_polar_splines * (1 + k_max - k_min), & | ||
| 518 | 403 | nr_nonzero=nr_nonzeros) | |
| 519 | |||
| 520 | s1 = 1 | ||
| 521 |
2/2✓ Branch 5 → 6 taken 1183 times.
✓ Branch 5 → 22 taken 403 times.
|
1586 | do k = k_min, k_max |
| 522 |
2/2✓ Branch 7 → 8 taken 30416 times.
✓ Branch 7 → 20 taken 1183 times.
|
32002 | do j = j_min, j_max |
| 523 |
2/2✓ Branch 9 → 10 taken 94446 times.
✓ Branch 9 → 18 taken 30416 times.
|
126045 | do i = i_min, i_max |
| 524 | 94446 | row_3d = this%index_list(1)%get_local_index(i, j, k) | |
| 525 | |||
| 526 |
2/2✓ Branch 12 → 13 taken 942462 times.
✓ Branch 12 → 16 taken 94446 times.
|
1067324 | do polar_index = 0, nr_polar_splines - 1 |
| 527 | 942462 | col = polar_index + (k - k_min) * nr_polar_splines | |
| 528 | 942462 | coeff = this%E1(i, j, s1, polar_index + 1) | |
| 529 |
2/2✓ Branch 13 → 11 taken 247562 times.
✓ Branch 13 → 14 taken 694900 times.
|
1036908 | if (abs(coeff) > COEFF_THR) call this%extractor(1)%insert(row_3d, col, coeff) |
| 530 | end do | ||
| 531 | end do | ||
| 532 | end do | ||
| 533 | end do | ||
| 534 | |||
| 535 | 403 | call this%extractor(1)%finalize() | |
| 536 | |||
| 537 | TIMER_STOP(t_total) | ||
| 538 | 403 | end subroutine finalize_polar_extractor_0forms | |
| 539 | |||
| 540 | !> @brief Initialize the 1-form extractor for the polar regularity constraint using gradients of polar zero-forms | ||
| 541 | !> | ||
| 542 | !> @param[inout] this The MFormPolarRegularity object to initialize | ||
| 543 | 333 | subroutine insert_polar_extractor_1forms(this) | |
| 544 | use m_tensorprod, only: size | ||
| 545 | use m_qr_factorisation, only: orthogonalize | ||
| 546 | use m_tensorprod_indices, only: TensorProdIndices | ||
| 547 | implicit none | ||
| 548 | |||
| 549 | class(MFormPolarRegularity), intent(inout) :: this | ||
| 550 | |||
| 551 | type(TensorProdIndices) :: inds | ||
| 552 | |||
| 553 | integer :: nr_polar_oneforms_xy, nr_polar_splines_z | ||
| 554 | integer :: nr_bspline_y | ||
| 555 | integer :: i, j, k, row, s1, polar_index, col, row0_2d(4), j_plus_one | ||
| 556 | integer :: i_min(3), i_max(3), j_min(3), j_max(3), k_min(3), k_max(3) | ||
| 557 | integer :: gamma, alpha, gamma_min(3), gamma_max | ||
| 558 | real(wp) :: coeff | ||
| 559 | |||
| 560 | TIMER_DECL(t_total) | ||
| 561 | TIMER_INIT_START(t_total, "MFormPolarRegularity::insert1") | ||
| 562 | |||
| 563 | ! NOTE: the y,z spaces are periodic and therefore the number of basis functions is the same | ||
| 564 | 333 | nr_bspline_y = size(this%space%tp_spaces(1), 2) | |
| 565 | |||
| 566 | 333 | gamma_min(1) = 1 | |
| 567 | 333 | gamma_min(2) = 1 | |
| 568 | gamma_min(3) = 0 | ||
| 569 | 333 | gamma_max = this%regularity | |
| 570 | |||
| 571 | 333 | nr_polar_oneforms_xy = get_nr_polar_zeroforms(gamma_max) - get_nr_polar_zeroforms(0) | |
| 572 | 666 | nr_polar_splines_z = merge(get_nr_polar_zeroforms(gamma_max), 0, this%space%nr_components == 3) | |
| 573 | |||
| 574 | !> Initialize the indices on which this extractor acts | ||
| 575 | 333 | row0_2d(1) = 0 | |
| 576 |
2/2✓ Branch 5 → 6 taken 999 times.
✓ Branch 5 → 11 taken 333 times.
|
1332 | do s1 = 1, 3 |
| 577 | 999 | i_min(s1) = max(0, this%space%tp_spaces(s1)%rank_resp_bspline%i0) | |
| 578 |
2/2✓ Branch 6 → 7 taken 666 times.
✓ Branch 6 → 8 taken 333 times.
|
999 | i_max(s1) = this%regularity - merge(1, 0, s1 == 1) |
| 579 | 999 | i_max(s1) = min(i_max(s1), this%space%tp_spaces(s1)%rank_resp_bspline%i1) | |
| 580 | 999 | j_min(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%j0 | |
| 581 | 999 | j_max(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%j1 | |
| 582 | 999 | k_min(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%k0 | |
| 583 | 999 | k_max(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%k1 | |
| 584 | 999 | call inds%init(i_min(s1), i_max(s1), j_min(s1), j_max(s1), k_min(s1), k_max(s1), this%space%tp_spaces(s1)%spaces) | |
| 585 | 999 | call this%index_list(s1)%init(inds, this%space%tp_spaces(s1)%rank_resp_bspline, this%space%tp_spaces(s1)%rank_l0) | |
| 586 | 1332 | row0_2d(s1 + 1) = row0_2d(s1) + size_2d(this%index_list(s1)) | |
| 587 | end do | ||
| 588 |
2/2✓ Branch 13 → 14 taken 1332 times.
✓ Branch 13 → 15 taken 333 times.
|
1665 | this%row0_2d = row0_2d |
| 589 | |||
| 590 | !> Initialize the radial and angular B-splines (data members of this object, using space0!) | ||
| 591 | call this%compute_1d_bsplines(this%space0%tp_spaces(1)%spaces(1), this%space0%tp_spaces(1)%spaces(2), 0, this%regularity, & | ||
| 592 | 333 | radial_truncate=this%regularity) | |
| 593 | |||
| 594 | ! xy-components | ||
| 595 | 333 | call this%init_shmem_window(this%E1, this%shmem_E1, nr_polar_oneforms_xy) | |
| 596 |
2/2✓ Branch 17 → 18 taken 666 times.
✓ Branch 17 → 38 taken 333 times.
|
999 | do s1 = 1, 2 |
| 597 |
2/2✓ Branch 19 → 20 taken 19528 times.
✓ Branch 19 → 36 taken 666 times.
|
20527 | do j = j_min(s1), j_max(s1) |
| 598 |
2/2✓ Branch 21 → 22 taken 32072 times.
✓ Branch 21 → 34 taken 19528 times.
|
52266 | do i = i_min(s1), i_max(s1) |
| 599 | polar_index = 0 | ||
| 600 | ! NOTE: we skip gamma == 0 because it's gradient is zero | ||
| 601 |
2/2✓ Branch 23 → 24 taken 73126 times.
✓ Branch 23 → 32 taken 32072 times.
|
124726 | do gamma = gamma_min(s1), gamma_max |
| 602 |
1/2✗ Branch 24 → 25 not taken.
✓ Branch 24 → 26 taken 73126 times.
|
105198 | do alpha = -gamma, gamma, 2 |
| 603 | 97176 | select case (s1) | |
| 604 | case (1) | ||
| 605 | coeff = (this%radial_bsplines(gamma)%data(i + 1) - this%radial_bsplines(gamma)%data(i)) * & | ||
| 606 | 97176 | this%angular_bsplines(alpha)%data(j) | |
| 607 | case (2) | ||
| 608 | 129400 | j_plus_one = mod(j + 1, nr_bspline_y) | |
| 609 | coeff = this%radial_bsplines(gamma)%data(i) * & | ||
| 610 |
2/2✓ Branch 26 → 27 taken 97176 times.
✓ Branch 26 → 28 taken 129400 times.
|
226576 | (this%angular_bsplines(alpha)%data(j_plus_one) - this%angular_bsplines(alpha)%data(j)) |
| 611 | end select | ||
| 612 |
2/2✓ Branch 29 → 26 taken 153450 times.
✓ Branch 29 → 30 taken 73126 times.
|
226576 | this%E1(i, j, s1, polar_index + 1) = coeff |
| 613 | 73126 | polar_index = polar_index + 1 | |
| 614 | end do | ||
| 615 | end do | ||
| 616 | end do | ||
| 617 | end do | ||
| 618 | end do | ||
| 619 | |||
| 620 | ! z-component | ||
| 621 |
2/2✓ Branch 39 → 40 taken 116 times.
✓ Branch 39 → 56 taken 217 times.
|
333 | if (this%space%nr_components == 3) then |
| 622 | 116 | call this%init_shmem_window(this%E2, this%shmem_E2, nr_polar_splines_z) | |
| 623 | |||
| 624 | s1 = 3 | ||
| 625 |
2/2✓ Branch 42 → 43 taken 3324 times.
✓ Branch 42 → 55 taken 116 times.
|
3440 | do j = j_min(s1), j_max(s1) |
| 626 |
2/2✓ Branch 44 → 45 taken 8736 times.
✓ Branch 44 → 53 taken 3324 times.
|
12176 | do i = i_min(s1), i_max(s1) |
| 627 | polar_index = 0 | ||
| 628 |
2/2✓ Branch 45 → 46 taken 28996 times.
✓ Branch 45 → 51 taken 8736 times.
|
41056 | do gamma = gamma_min(s1), gamma_max |
| 629 |
1/2✗ Branch 46 → 47 not taken.
✓ Branch 46 → 48 taken 28996 times.
|
37732 | do alpha = -gamma, gamma, 2 |
| 630 | this%E2(i, j, s1, polar_index + 1) = & | ||
| 631 |
2/2✓ Branch 48 → 48 taken 41452 times.
✓ Branch 48 → 49 taken 28996 times.
|
70448 | this%radial_bsplines(gamma)%data(i) * this%angular_bsplines(alpha)%data(j) |
| 632 | 28996 | polar_index = polar_index + 1 | |
| 633 | end do | ||
| 634 | end do | ||
| 635 | end do | ||
| 636 | end do | ||
| 637 | end if | ||
| 638 | |||
| 639 | TIMER_STOP(t_total) | ||
| 640 | 333 | end subroutine insert_polar_extractor_1forms | |
| 641 | |||
| 642 | !> @brief Enrich 1-form basis for the polar regularity constraint using weak curls of polar two-forms | ||
| 643 | !> | ||
| 644 | !> @param[inout] this The MFormPolarEnriched object to initialize | ||
| 645 | 106 | subroutine insert_polar_extractor_1forms_enrich(this) | |
| 646 | use m_bspline_basis, only: BSplineFun, BSplineSpace, evaluate | ||
| 647 | use m_bspline_linalg, only: l2_projection, biorthogonalize | ||
| 648 | use m_common, only: zero_function, wp, pi, VERBOSITY_WARN_ALWAYS | ||
| 649 | use m_coord_transform, only: CylinderTransform | ||
| 650 | use m_mform_basis, only: MFormSpace, get_zero_form_space, evaluate, MFormFun | ||
| 651 | use m_mform_constraint_boundary, only: MFormDirichlet, MFormNeumann | ||
| 652 | use m_mform_derham, only: DeRhamSequence | ||
| 653 | use m_mform_linalg, only: inner_product, l2_norm, l2_error, solve_curlA_equals_B | ||
| 654 | use m_mform_solver, only: GenericSolver | ||
| 655 | use m_mform_matrix, only: DiffDiffMatrix, MassMatrix | ||
| 656 | use m_sparsemat, only: SparseMat | ||
| 657 | use m_tensorprod, only: size | ||
| 658 | use m_domain_decomp, only: DomainDecomp, memory_layout_convert | ||
| 659 | implicit none | ||
| 660 | |||
| 661 | class(MFormPolarRegularity), intent(inout) :: this | ||
| 662 | |||
| 663 | type(CylinderTransform) :: disk | ||
| 664 | 106 | type(GenericSolver) :: divgrad_solver_2d | |
| 665 | 106 | type(DiffDiffMatrix) :: divgrad_matrix_2d | |
| 666 | type(DomainDecomp) :: domain_2d | ||
| 667 |
6/6✓ Branch 3 → 4 taken 318 times.
✓ Branch 3 → 5 taken 106 times.
✓ Branch 280 → 281 taken 318 times.
✓ Branch 280 → 284 taken 106 times.
✓ Branch 281 → 282 taken 210 times.
✓ Branch 281 → 283 taken 108 times.
|
848 | type(MFormSpace) :: derham_2d(0:2) |
| 668 | 106 | type(MFormFun) :: f, tmp | |
| 669 | type(MFormFun), allocatable :: polar_oneforms(:) | ||
| 670 | 106 | type(MFormPolarRegularity) :: bounded_oneform_constraint | |
| 671 | |||
| 672 | type(BSplineSpace) :: bspline_x, bspline_y | ||
| 673 | |||
| 674 | 106 | real(wp), pointer :: E1_curl(:, :, :, :) | |
| 675 | type(SharedMemoryWindow) :: shmem_E1_curl | ||
| 676 | |||
| 677 | integer :: comm_2d, ierr, comm_color, comm_key | ||
| 678 | integer :: nr_enriched_basis, nr_gradient_basis, gamma, alpha, gamma_min, gamma_max | ||
| 679 | integer :: s1, nr_nonzeros(3), i, j, row, polar_index, nr_x, nr_bspline_y | ||
| 680 | integer :: i_min(3), i_max(3), j_min(3), j_max(3) | ||
| 681 | |||
| 682 | TIMER_DECL(t_total) | ||
| 683 | TIMER_INIT_START(t_total, "MFormPolarRegularity::enrich1") | ||
| 684 | |||
| 685 | gamma_min = 0 | ||
| 686 | |||
| 687 | ! The pushforward yields extra factor xp, and the 3rd component of the 2-form space has degree reduced by 1 | ||
| 688 | 106 | gamma_max = this%regularity - 2 | |
| 689 | 106 | if (gamma_max < gamma_min) then | |
| 690 | TIMER_STOP(t_total) | ||
| 691 | return | ||
| 692 | end if | ||
| 693 | |||
| 694 | 70 | disk = CylinderTransform(radius=1.0_wp) | |
| 695 | |||
| 696 | ! Create a 2D domain decomposition for the disk cross-section, respecting the original 3D decomposition | ||
| 697 | 70 | comm_color = this%space%tp_spaces(1)%domain%my_subinterval_ijk(3) | |
| 698 | comm_key = this%space%tp_spaces(1)%domain%my_subinterval_ijk(1) + this%space%tp_spaces(1)%domain%my_subinterval_ijk(2) * & | ||
| 699 | 70 | this%space%tp_spaces(1)%domain%nr_subintervals(1) | |
| 700 |
1/2✗ Branch 9 → 10 not taken.
✓ Branch 9 → 11 taken 70 times.
|
70 | PetscCallMPI(MPI_Comm_split(this%space%tp_spaces(1)%domain%comm, comm_color, comm_key, comm_2d, ierr)) |
| 701 | domain_2d = DomainDecomp(& | ||
| 702 | & Nx=this%space%tp_spaces(1)%domain%nr_subintervals(1), & | ||
| 703 | & Ny=this%space%tp_spaces(1)%domain%nr_subintervals(2), & | ||
| 704 | & memory_layout_x=memory_layout_convert(this%space%tp_spaces(1)%domain%memory_layout(1)), & | ||
| 705 | & memory_layout_y=memory_layout_convert(this%space%tp_spaces(1)%domain%memory_layout(2)), & | ||
| 706 | 70 | comm=comm_2d, flat_mpi=this%space%tp_spaces(1)%domain%flat_mpi, dimensionality=2) | |
| 707 | |||
| 708 | ! Get the 2D DeRham sequence corresponding to the disk cross-section | ||
| 709 |
3/4✓ Branch 15 → 16 taken 210 times.
✓ Branch 15 → 19 taken 70 times.
✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 210 times.
|
280 | derham_2d = DeRhamSequence(this%space0%tp_spaces(1)%spaces(1), this%space0%tp_spaces(1)%spaces(2), domain=domain_2d) |
| 710 | |||
| 711 | ! Set up the curl-curl and divgrad solver in 2D | ||
| 712 | 70 | call divgrad_matrix_2d%init(derham_2d(0), disk) | |
| 713 | call divgrad_solver_2d%init(divgrad_matrix_2d, constraint1=MFormNeumann(), constraint2=MFormPolarBoundedness(), & | ||
| 714 |
5/36✓ Branch 21 → 22 taken 70 times.
✗ Branch 21 → 23 not taken.
✗ Branch 24 → 25 not taken.
✓ Branch 24 → 26 taken 70 times.
✗ Branch 26 → 27 not taken.
✓ Branch 26 → 39 taken 70 times.
✗ Branch 28 → 29 not taken.
✗ Branch 28 → 38 not taken.
✗ Branch 29 → 30 not taken.
✗ Branch 29 → 31 not taken.
✗ Branch 31 → 32 not taken.
✗ Branch 31 → 33 not taken.
✗ Branch 33 → 34 not taken.
✗ Branch 33 → 35 not taken.
✗ Branch 35 → 36 not taken.
✗ Branch 35 → 37 not taken.
✗ Branch 39 → 40 not taken.
✓ Branch 39 → 50 taken 70 times.
✗ Branch 41 → 42 not taken.
✗ Branch 41 → 49 not taken.
✗ Branch 42 → 43 not taken.
✗ Branch 42 → 44 not taken.
✗ Branch 44 → 45 not taken.
✗ Branch 44 → 46 not taken.
✗ Branch 46 → 47 not taken.
✗ Branch 46 → 48 not taken.
✗ Branch 50 → 51 not taken.
✓ Branch 50 → 61 taken 70 times.
✗ Branch 52 → 53 not taken.
✗ Branch 52 → 60 not taken.
✗ Branch 53 → 54 not taken.
✗ Branch 53 → 55 not taken.
✗ Branch 55 → 56 not taken.
✗ Branch 55 → 57 not taken.
✗ Branch 57 → 58 not taken.
✗ Branch 57 → 59 not taken.
|
140 | coord_transform=disk) |
| 715 | |||
| 716 |
7/48✗ Branch 61 → 62 not taken.
✓ Branch 61 → 63 taken 70 times.
✗ Branch 63 → 64 not taken.
✓ Branch 63 → 76 taken 70 times.
✗ Branch 65 → 66 not taken.
✗ Branch 65 → 75 not taken.
✗ Branch 66 → 67 not taken.
✗ Branch 66 → 68 not taken.
✗ Branch 68 → 69 not taken.
✗ Branch 68 → 70 not taken.
✗ Branch 70 → 71 not taken.
✗ Branch 70 → 72 not taken.
✗ Branch 72 → 73 not taken.
✗ Branch 72 → 74 not taken.
✗ Branch 76 → 77 not taken.
✓ Branch 76 → 87 taken 70 times.
✗ Branch 78 → 79 not taken.
✗ Branch 78 → 86 not taken.
✗ Branch 79 → 80 not taken.
✗ Branch 79 → 81 not taken.
✗ Branch 81 → 82 not taken.
✗ Branch 81 → 83 not taken.
✗ Branch 83 → 84 not taken.
✗ Branch 83 → 85 not taken.
✗ Branch 87 → 88 not taken.
✓ Branch 87 → 98 taken 70 times.
✗ Branch 89 → 90 not taken.
✗ Branch 89 → 97 not taken.
✗ Branch 90 → 91 not taken.
✗ Branch 90 → 92 not taken.
✗ Branch 92 → 93 not taken.
✗ Branch 92 → 94 not taken.
✗ Branch 94 → 95 not taken.
✗ Branch 94 → 96 not taken.
✗ Branch 98 → 99 not taken.
✓ Branch 98 → 105 taken 70 times.
✗ Branch 100 → 101 not taken.
✗ Branch 100 → 104 not taken.
✗ Branch 101 → 102 not taken.
✗ Branch 101 → 103 not taken.
✗ Branch 105 → 106 not taken.
✓ Branch 105 → 112 taken 70 times.
✗ Branch 107 → 108 not taken.
✗ Branch 107 → 111 not taken.
✗ Branch 108 → 109 not taken.
✗ Branch 108 → 110 not taken.
✗ Branch 112 → 113 not taken.
✓ Branch 112 → 114 taken 70 times.
|
70 | bounded_oneform_constraint = MFormPolarBoundedness() |
| 717 | 70 | call bounded_oneform_constraint%init(derham_2d(1), coord_transform=disk) | |
| 718 | |||
| 719 | 70 | nr_gradient_basis = size(this%E1, 4) | |
| 720 | 70 | nr_enriched_basis = get_nr_polar_zeroforms(gamma_max) - get_nr_polar_zeroforms(gamma_min - 1) | |
| 721 |
6/10✓ Branch 117 → 116 taken 70 times.
✗ Branch 117 → 118 not taken.
✓ Branch 118 → 119 taken 70 times.
✗ Branch 118 → 120 not taken.
✗ Branch 120 → 121 not taken.
✓ Branch 120 → 122 taken 70 times.
✗ Branch 122 → 123 not taken.
✓ Branch 122 → 124 taken 70 times.
✓ Branch 125 → 126 taken 214 times.
✓ Branch 125 → 127 taken 70 times.
|
424 | allocate (polar_oneforms(0:nr_enriched_basis - 1)) |
| 722 | |||
| 723 |
2/2✓ Branch 128 → 129 taken 210 times.
✓ Branch 128 → 131 taken 70 times.
|
280 | do s1 = 1, 3 |
| 724 | 280 | call this%index_list(s1)%get_bounds(i_min(s1), i_max(s1), j_min(s1), j_max(s1)) | |
| 725 | end do | ||
| 726 | |||
| 727 | ! Compute L2 approximation of the radial and angular functions which make up the Zernike basis | ||
| 728 | ! (see also MFormPolarRegularity::compute_1d_bsplines) | ||
| 729 | 70 | bspline_x = derham_2d(2)%tp_spaces(3)%spaces(1) ! Radial B-spline space for 2-forms | |
| 730 | 70 | bspline_y = derham_2d(2)%tp_spaces(3)%spaces(2) ! Angular B-spline space for 2-forms | |
| 731 | |||
| 732 | ! NB we let gamma_max = gamma_max + 1 here because the radial part of the 2-form has an extra factor r | ||
| 733 | 70 | call this%compute_1d_bsplines(bspline_x, bspline_y, gamma_min=1, gamma_max=gamma_max + 1, radial_truncate=gamma_max + 1) | |
| 734 | |||
| 735 | polar_index = 0 | ||
| 736 |
2/2✓ Branch 134 → 135 taken 130 times.
✓ Branch 134 → 142 taken 70 times.
|
200 | do gamma = gamma_min, gamma_max |
| 737 |
1/2✗ Branch 135 → 136 not taken.
✓ Branch 135 → 137 taken 130 times.
|
200 | do alpha = -gamma, gamma, 2 |
| 738 | call solve_curlA_equals_B(tmp, this%radial_bsplines(gamma + 1), this%angular_bsplines(alpha), & | ||
| 739 | 214 | twoform_space=derham_2d(2), cylinder_transform=disk, zeroform_solver=divgrad_solver_2d) | |
| 740 | |||
| 741 | 214 | call bounded_oneform_constraint%apply(polar_oneforms(polar_index), tmp) | |
| 742 |
2/2✓ Branch 139 → 137 taken 84 times.
✓ Branch 139 → 140 taken 130 times.
|
214 | polar_index = polar_index + 1 |
| 743 | end do | ||
| 744 | end do | ||
| 745 | |||
| 746 | ! Set up the enriched extractor matrices; only for the ranks that are at the polar axis | ||
| 747 | 70 | call this%init_shmem_window(E1_curl, shmem_E1_curl, nr_enriched_basis) | |
| 748 |
2/2✓ Branch 144 → 145 taken 140 times.
✓ Branch 144 → 157 taken 70 times.
|
210 | do s1 = 1, 2 |
| 749 |
2/2✓ Branch 146 → 147 taken 3764 times.
✓ Branch 146 → 155 taken 140 times.
|
3974 | do j = j_min(s1), j_max(s1) |
| 750 |
2/2✓ Branch 148 → 149 taken 12638 times.
✓ Branch 148 → 153 taken 3764 times.
|
16542 | do i = i_min(s1), i_max(s1) |
| 751 |
2/2✓ Branch 149 → 150 taken 46852 times.
✓ Branch 149 → 151 taken 12638 times.
|
63254 | do polar_index = 0, nr_enriched_basis - 1 |
| 752 | 59490 | E1_curl(i, j, s1, polar_index + 1) = polar_oneforms(polar_index)%tp_funs(s1)%data(i, j, 0) | |
| 753 | end do | ||
| 754 | end do | ||
| 755 | end do | ||
| 756 | end do | ||
| 757 | |||
| 758 | ! Concatenate the shared memory windows | ||
| 759 | 70 | call this%shmem_E1%cat(shmem_E1_curl) | |
| 760 | 70 | nr_bspline_y = size(this%space%tp_spaces(1), 2) | |
| 761 | nr_x = merge(this%regularity + 1, 0, this%space%tp_spaces(1)%domain%is_shared_memory(1) .or. & | ||
| 762 |
2/4✗ Branch 159 → 160 not taken.
✓ Branch 159 → 161 taken 70 times.
✓ Branch 161 → 160 taken 70 times.
✗ Branch 161 → 162 not taken.
|
70 | this%space%tp_spaces(1)%domain%my_subinterval_ijk(1) == 0) |
| 763 | 70 | this%E1(0:nr_x - 1, 0:nr_bspline_y - 1, 1:3, 1:(nr_gradient_basis + nr_enriched_basis)) => this%shmem_E1%window(:) | |
| 764 | |||
| 765 | ! Clean up | ||
| 766 | 70 | call shmem_E1_curl%destroy() | |
| 767 | 70 | nullify (E1_curl) | |
| 768 | |||
| 769 |
1/2✗ Branch 164 → 165 not taken.
✓ Branch 164 → 166 taken 70 times.
|
70 | PetscCallMPI(MPI_Comm_free(comm_2d, ierr)) |
| 770 | 70 | call divgrad_matrix_2d%destroy() | |
| 771 | 70 | call divgrad_solver_2d%destroy() | |
| 772 | |||
| 773 | polar_index = 0 | ||
| 774 |
2/2✓ Branch 169 → 170 taken 130 times.
✓ Branch 169 → 176 taken 70 times.
|
200 | do gamma = gamma_min, gamma_max |
| 775 |
1/2✗ Branch 170 → 171 not taken.
✓ Branch 170 → 172 taken 130 times.
|
200 | do alpha = -gamma, gamma, 2 |
| 776 | 214 | call polar_oneforms(polar_index)%destroy() | |
| 777 |
2/2✓ Branch 173 → 172 taken 84 times.
✓ Branch 173 → 174 taken 130 times.
|
214 | polar_index = polar_index + 1 |
| 778 | end do | ||
| 779 | end do | ||
| 780 | |||
| 781 |
6/10✓ Branch 177 → 178 taken 70 times.
✗ Branch 177 → 184 not taken.
✓ Branch 178 → 179 taken 214 times.
✓ Branch 178 → 184 taken 70 times.
✓ Branch 179 → 180 taken 214 times.
✗ Branch 179 → 181 not taken.
✗ Branch 181 → 182 not taken.
✓ Branch 181 → 183 taken 214 times.
✗ Branch 184 → 185 not taken.
✓ Branch 184 → 186 taken 70 times.
|
284 | deallocate (polar_oneforms) |
| 782 | |||
| 783 | TIMER_STOP(t_total) | ||
| 784 |
46/138✓ Branch 5 → 6 taken 36 times.
✓ Branch 5 → 7 taken 70 times.
✓ Branch 187 → 188 taken 70 times.
✓ Branch 187 → 189 taken 36 times.
✓ Branch 189 → 190 taken 70 times.
✓ Branch 189 → 191 taken 36 times.
✗ Branch 191 → 192 not taken.
✓ Branch 191 → 199 taken 106 times.
✗ Branch 192 → 193 not taken.
✗ Branch 192 → 198 not taken.
✗ Branch 193 → 194 not taken.
✗ Branch 193 → 195 not taken.
✗ Branch 195 → 196 not taken.
✗ Branch 195 → 197 not taken.
✓ Branch 199 → 200 taken 70 times.
✓ Branch 199 → 201 taken 36 times.
✓ Branch 201 → 202 taken 70 times.
✓ Branch 201 → 203 taken 36 times.
✗ Branch 203 → 204 not taken.
✓ Branch 203 → 275 taken 106 times.
✗ Branch 204 → 205 not taken.
✗ Branch 204 → 206 not taken.
✗ Branch 207 → 208 not taken.
✗ Branch 207 → 251 not taken.
✗ Branch 208 → 209 not taken.
✗ Branch 208 → 211 not taken.
✗ Branch 209 → 210 not taken.
✗ Branch 209 → 211 not taken.
✗ Branch 211 → 212 not taken.
✗ Branch 211 → 250 not taken.
✗ Branch 212 → 213 not taken.
✗ Branch 212 → 214 not taken.
✗ Branch 214 → 215 not taken.
✗ Branch 214 → 227 not taken.
✗ Branch 216 → 217 not taken.
✗ Branch 216 → 226 not taken.
✗ Branch 217 → 218 not taken.
✗ Branch 217 → 219 not taken.
✗ Branch 219 → 220 not taken.
✗ Branch 219 → 221 not taken.
✗ Branch 221 → 222 not taken.
✗ Branch 221 → 223 not taken.
✗ Branch 223 → 224 not taken.
✗ Branch 223 → 225 not taken.
✗ Branch 227 → 228 not taken.
✗ Branch 227 → 238 not taken.
✗ Branch 229 → 230 not taken.
✗ Branch 229 → 237 not taken.
✗ Branch 230 → 231 not taken.
✗ Branch 230 → 232 not taken.
✗ Branch 232 → 233 not taken.
✗ Branch 232 → 234 not taken.
✗ Branch 234 → 235 not taken.
✗ Branch 234 → 236 not taken.
✗ Branch 238 → 239 not taken.
✗ Branch 238 → 249 not taken.
✗ Branch 240 → 241 not taken.
✗ Branch 240 → 248 not taken.
✗ Branch 241 → 242 not taken.
✗ Branch 241 → 243 not taken.
✗ Branch 243 → 244 not taken.
✗ Branch 243 → 245 not taken.
✗ Branch 245 → 246 not taken.
✗ Branch 245 → 247 not taken.
✗ Branch 251 → 252 not taken.
✗ Branch 251 → 253 not taken.
✗ Branch 253 → 254 not taken.
✗ Branch 253 → 255 not taken.
✗ Branch 255 → 256 not taken.
✗ Branch 255 → 257 not taken.
✗ Branch 257 → 258 not taken.
✗ Branch 257 → 259 not taken.
✗ Branch 259 → 260 not taken.
✗ Branch 259 → 261 not taken.
✗ Branch 261 → 262 not taken.
✗ Branch 261 → 274 not taken.
✗ Branch 263 → 264 not taken.
✗ Branch 263 → 273 not taken.
✗ Branch 264 → 265 not taken.
✗ Branch 264 → 266 not taken.
✗ Branch 266 → 267 not taken.
✗ Branch 266 → 268 not taken.
✗ Branch 268 → 269 not taken.
✗ Branch 268 → 270 not taken.
✗ Branch 270 → 271 not taken.
✗ Branch 270 → 272 not taken.
✓ Branch 275 → 276 taken 70 times.
✓ Branch 275 → 277 taken 36 times.
✓ Branch 277 → 278 taken 70 times.
✓ Branch 277 → 279 taken 36 times.
✓ Branch 284 → 285 taken 70 times.
✓ Branch 284 → 286 taken 36 times.
✓ Branch 286 → 287 taken 70 times.
✓ Branch 286 → 299 taken 36 times.
✓ Branch 288 → 289 taken 210 times.
✓ Branch 288 → 298 taken 70 times.
✓ Branch 289 → 290 taken 210 times.
✗ Branch 289 → 291 not taken.
✓ Branch 291 → 292 taken 210 times.
✗ Branch 291 → 293 not taken.
✓ Branch 293 → 294 taken 210 times.
✗ Branch 293 → 295 not taken.
✓ Branch 295 → 296 taken 210 times.
✗ Branch 295 → 297 not taken.
✗ Branch 299 → 300 not taken.
✓ Branch 299 → 310 taken 106 times.
✗ Branch 301 → 302 not taken.
✗ Branch 301 → 309 not taken.
✗ Branch 302 → 303 not taken.
✗ Branch 302 → 304 not taken.
✗ Branch 304 → 305 not taken.
✗ Branch 304 → 306 not taken.
✗ Branch 306 → 307 not taken.
✗ Branch 306 → 308 not taken.
✓ Branch 310 → 311 taken 70 times.
✓ Branch 310 → 321 taken 36 times.
✓ Branch 312 → 313 taken 210 times.
✓ Branch 312 → 320 taken 70 times.
✓ Branch 313 → 314 taken 210 times.
✗ Branch 313 → 315 not taken.
✓ Branch 315 → 316 taken 210 times.
✗ Branch 315 → 317 not taken.
✓ Branch 317 → 318 taken 210 times.
✗ Branch 317 → 319 not taken.
✓ Branch 321 → 322 taken 70 times.
✓ Branch 321 → 328 taken 36 times.
✓ Branch 323 → 324 taken 70 times.
✓ Branch 323 → 327 taken 70 times.
✓ Branch 324 → 325 taken 70 times.
✗ Branch 324 → 326 not taken.
✓ Branch 328 → 329 taken 70 times.
✓ Branch 328 → 335 taken 36 times.
✓ Branch 330 → 331 taken 70 times.
✓ Branch 330 → 334 taken 70 times.
✓ Branch 331 → 332 taken 70 times.
✗ Branch 331 → 333 not taken.
✓ Branch 335 → 336 taken 70 times.
✓ Branch 335 → 337 taken 36 times.
|
1018 | end subroutine insert_polar_extractor_1forms_enrich |
| 785 | |||
| 786 | 333 | subroutine finalize_polar_extractor_1forms(this) | |
| 787 | implicit none | ||
| 788 | |||
| 789 | class(MFormPolarRegularity), intent(inout) :: this | ||
| 790 | |||
| 791 | integer :: i, j, k, row_3d, s1, polar_index, col, nr_polar_splines | ||
| 792 | integer :: i_min(3), i_max(3), j_min(3), j_max(3), k_min(3), k_max(3), col0, nr_nonzeros(3) | ||
| 793 | real(wp) :: coeff | ||
| 794 | |||
| 795 | TIMER_DECL(t_total) | ||
| 796 | TIMER_INIT_START(t_total, "MFormPolarRegularity::finalize1") | ||
| 797 | |||
| 798 |
2/2✓ Branch 3 → 4 taken 999 times.
✓ Branch 3 → 6 taken 333 times.
|
1332 | do s1 = 1, 3 |
| 799 | 1332 | call this%index_list(s1)%get_bounds(i_min(s1), i_max(s1), j_min(s1), j_max(s1), k_min(s1), k_max(s1)) | |
| 800 | end do | ||
| 801 | |||
| 802 | 333 | nr_nonzeros(1) = size(this%index_list(1)) * size(this%E1, 4) | |
| 803 | 333 | nr_nonzeros(2) = size(this%index_list(2)) * size(this%E1, 4) | |
| 804 |
2/2✓ Branch 7 → 8 taken 217 times.
✓ Branch 7 → 9 taken 116 times.
|
333 | if (this%space%nr_components == 2) then |
| 805 | 217 | nr_nonzeros(3) = 0 | |
| 806 | nr_polar_splines = size(this%E1, 4) | ||
| 807 | else | ||
| 808 | 116 | nr_nonzeros(3) = size(this%index_list(3)) * size(this%E2, 4) | |
| 809 | 116 | nr_polar_splines = size(this%E1, 4) + size(this%E2, 4) | |
| 810 | end if | ||
| 811 |
2/2✓ Branch 11 → 12 taken 999 times.
✓ Branch 11 → 14 taken 333 times.
|
1332 | do s1 = 1, 3 |
| 812 | ! TODO is using k_max(1) OK with MPI? nr cols needs to be same for all blocks | ||
| 813 | call this%extractor(s1)%init(nr_rows=size(this%index_list(s1)), nr_cols=nr_polar_splines * & | ||
| 814 | 1332 | & size(this%index_list(1), dir=3), nr_nonzero=nr_nonzeros(s1)) | |
| 815 | end do | ||
| 816 | |||
| 817 |
2/2✓ Branch 15 → 16 taken 999 times.
✓ Branch 15 → 45 taken 333 times.
|
1332 | do s1 = 1, 3 |
| 818 |
4/4✓ Branch 16 → 17 taken 651 times.
✓ Branch 16 → 18 taken 348 times.
✓ Branch 17 → 18 taken 434 times.
✓ Branch 17 → 43 taken 217 times.
|
999 | if (this%space%nr_components == 2 .and. s1 == 3) cycle |
| 819 | |||
| 820 | 666 | select case (s1) | |
| 821 | case (1, 2) | ||
| 822 | 666 | nr_polar_splines = size(this%E1, 4) | |
| 823 | 666 | col0 = 0 | |
| 824 | case (3) | ||
| 825 | 116 | nr_polar_splines = size(this%E2, 4) | |
| 826 |
2/2✓ Branch 18 → 19 taken 666 times.
✓ Branch 18 → 20 taken 116 times.
|
782 | col0 = size(this%E1, 4) * size(this%index_list(1), dir=3) |
| 827 | end select | ||
| 828 | |||
| 829 |
2/2✓ Branch 22 → 23 taken 782 times.
✓ Branch 22 → 24 taken 2342 times.
|
3457 | do k = k_min(s1), k_max(s1) |
| 830 |
2/2✓ Branch 25 → 26 taken 60904 times.
✓ Branch 25 → 41 taken 2342 times.
|
64245 | do j = j_min(s1), j_max(s1) |
| 831 |
2/2✓ Branch 27 → 28 taken 137060 times.
✓ Branch 27 → 39 taken 60904 times.
|
200306 | do i = i_min(s1), i_max(s1) |
| 832 | 137060 | row_3d = this%index_list(s1)%get_local_index(i, j, k) | |
| 833 | polar_index = 0 | ||
| 834 | ! NOTE: we skip gamma == 0 because it's gradient is zero | ||
| 835 |
2/2✓ Branch 29 → 30 taken 1266040 times.
✓ Branch 29 → 37 taken 137060 times.
|
1464004 | do polar_index = 0, nr_polar_splines - 1 |
| 836 | 1266040 | col = polar_index + (k - k_min(s1)) * nr_polar_splines | |
| 837 | |||
| 838 | 881408 | select case (s1) | |
| 839 | case (1, 2) | ||
| 840 | 881408 | coeff = this%E1(i, j, s1, polar_index + 1) | |
| 841 | case (3) | ||
| 842 |
2/2✓ Branch 30 → 31 taken 881408 times.
✓ Branch 30 → 32 taken 384632 times.
|
1266040 | coeff = this%E2(i, j, s1, polar_index + 1) |
| 843 | end select | ||
| 844 | |||
| 845 |
2/2✓ Branch 33 → 34 taken 754205 times.
✓ Branch 33 → 36 taken 511835 times.
|
1403100 | if (abs(coeff) > COEFF_THR) call this%extractor(s1)%insert(row_3d, col0 + col, coeff) |
| 846 | end do | ||
| 847 | end do | ||
| 848 | end do | ||
| 849 | end do | ||
| 850 | end do | ||
| 851 | |||
| 852 |
2/2✓ Branch 46 → 47 taken 999 times.
✓ Branch 46 → 49 taken 333 times.
|
1332 | do s1 = 1, 3 |
| 853 | 1332 | call this%extractor(s1)%finalize() | |
| 854 | end do | ||
| 855 | TIMER_STOP(t_total) | ||
| 856 | 333 | end subroutine finalize_polar_extractor_1forms | |
| 857 | |||
| 858 | !> @brief Initialize the 2-form extractor for the polar regularity constraint with curls of polar one-forms | ||
| 859 | !> | ||
| 860 | !> @param[inout] this The MFormPolarRegularity object to initialize | ||
| 861 | 221 | subroutine insert_polar_extractor_2forms(this) | |
| 862 | use m_tensorprod, only: size | ||
| 863 | use m_tensorprod_indices, only: TensorProdIndices | ||
| 864 | use m_qr_factorisation, only: dense_column_space | ||
| 865 | implicit none | ||
| 866 | |||
| 867 | class(MFormPolarRegularity), intent(inout) :: this | ||
| 868 | |||
| 869 | type(TensorProdIndices) :: inds | ||
| 870 | |||
| 871 | integer :: nr_nonzeros(3), nr_polar_splines(3), nr_polar_splines_xy, nr_polar_splines_z, nr_bspline_y | ||
| 872 | integer :: i, j, j_plus_one, s1, polar_index, col0, col, row_xy, row_xy0, non_zeros_xy | ||
| 873 | integer :: i_min(3), i_max(3), j_min(3), j_max(3), k_min(3), k_max(3) | ||
| 874 | integer :: gamma, alpha, gamma_min, gamma_max | ||
| 875 | real(wp) :: coeff | ||
| 876 | |||
| 877 | TIMER_DECL(t_total) | ||
| 878 | TIMER_INIT_START(t_total, "MFormPolarRegularity::insert2") | ||
| 879 | |||
| 880 | 221 | nr_bspline_y = size(this%space%tp_spaces(1), 2) | |
| 881 | |||
| 882 | gamma_min = 1 | ||
| 883 | 221 | gamma_max = this%regularity | |
| 884 | |||
| 885 | !> Initialize the indices on which this extractor acts | ||
| 886 | 221 | i_max(1) = min(this%regularity, this%space%tp_spaces(1)%rank_resp_bspline%i1) | |
| 887 | 221 | i_max(2) = min(this%regularity - 1, this%space%tp_spaces(2)%rank_resp_bspline%i1) | |
| 888 | ! NOTE the 3rd component gets at least 0, because of the pushforward factor 1/xp | ||
| 889 | 221 | i_max(3) = min(this%regularity, this%space%tp_spaces(3)%spaces(1)%degree, this%space%tp_spaces(3)%rank_resp_bspline%i1) | |
| 890 | 221 | this%row0_2d(1) = 0 | |
| 891 |
2/2✓ Branch 3 → 4 taken 663 times.
✓ Branch 3 → 7 taken 221 times.
|
884 | do s1 = 1, 3 |
| 892 | 663 | i_min(s1) = max(0, this%space%tp_spaces(s1)%rank_resp_bspline%i0) | |
| 893 | 663 | j_min(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%j0; j_max(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%j1 | |
| 894 | 663 | k_min(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%k0; k_max(s1) = this%space%tp_spaces(s1)%rank_resp_bspline%k1 | |
| 895 | |||
| 896 | 663 | call inds%init(i_min(s1), i_max(s1), j_min(s1), j_max(s1), k_min(s1), k_max(s1), this%space%tp_spaces(s1)%spaces) | |
| 897 | 663 | call this%index_list(s1)%init(inds, this%space%tp_spaces(s1)%rank_resp_bspline, this%space%tp_spaces(s1)%rank_l0) | |
| 898 | 884 | this%row0_2d(s1 + 1) = this%row0_2d(s1) + size_2d(this%index_list(s1)) | |
| 899 | end do | ||
| 900 | |||
| 901 | ! Assume no enrichment: can be overridden later | ||
| 902 | 221 | call this%init_shmem_window(this%E2, this%shmem_E2, 0) | |
| 903 | |||
| 904 |
2/2✓ Branch 9 → 10 taken 111 times.
✓ Branch 9 → 12 taken 110 times.
|
221 | if (this%space%nr_components == 2) then |
| 905 | 111 | call this%init_shmem_window(this%E1, this%shmem_E1, 0) | |
| 906 | TIMER_STOP(t_total) | ||
| 907 | 111 | return | |
| 908 | end if | ||
| 909 | |||
| 910 | !> Initialize the radial and angular B-splines (data members of this object, using space0!) | ||
| 911 | call this%compute_1d_bsplines(this%space0%tp_spaces(1)%spaces(1), this%space0%tp_spaces(1)%spaces(2), 0, this%regularity, & | ||
| 912 | 110 | radial_truncate=this%regularity) | |
| 913 | |||
| 914 | ! Curls of the Z-component of polar one-forms; resulting in the xy-components of the 2-form space | ||
| 915 | 110 | nr_polar_splines_xy = get_nr_polar_zeroforms(gamma_max) - get_nr_polar_zeroforms(0) | |
| 916 | |||
| 917 | 110 | call this%init_shmem_window(this%E1, this%shmem_E1, nr_polar_splines_xy) | |
| 918 |
2/2✓ Branch 14 → 15 taken 220 times.
✓ Branch 14 → 36 taken 110 times.
|
330 | do s1 = 1, 2 |
| 919 |
2/2✓ Branch 16 → 17 taken 6504 times.
✓ Branch 16 → 34 taken 220 times.
|
6834 | do j = j_min(s1), j_max(s1) |
| 920 |
2/2✓ Branch 18 → 19 taken 13764 times.
✓ Branch 18 → 32 taken 6504 times.
|
20488 | do i = i_min(s1), i_max(s1) |
| 921 | polar_index = 0 | ||
| 922 |
2/2✓ Branch 19 → 20 taken 34184 times.
✓ Branch 19 → 30 taken 13764 times.
|
54452 | do gamma = gamma_min, gamma_max |
| 923 |
1/2✗ Branch 20 → 21 not taken.
✓ Branch 20 → 22 taken 34184 times.
|
47948 | do alpha = -gamma, gamma, 2 |
| 924 | |||
| 925 | 60200 | select case (s1) | |
| 926 | case (1) | ||
| 927 | 60200 | j_plus_one = mod(j + 1, nr_bspline_y) | |
| 928 | coeff = (this%radial_bsplines(gamma)%data(i) * & | ||
| 929 | 60200 | (this%angular_bsplines(alpha)%data(j_plus_one) - this%angular_bsplines(alpha)%data(j))) | |
| 930 | case (2) | ||
| 931 | coeff = -(this%radial_bsplines(gamma)%data(i + 1) - this%radial_bsplines(gamma)%data(i)) * & | ||
| 932 |
2/2✓ Branch 22 → 23 taken 60200 times.
✓ Branch 22 → 24 taken 45084 times.
|
105284 | this%angular_bsplines(alpha)%data(j) |
| 933 | end select | ||
| 934 |
2/2✓ Branch 25 → 26 taken 49811 times.
✓ Branch 25 → 27 taken 55473 times.
|
105284 | if (abs(coeff) > COEFF_THR) this%E1(i, j, s1, polar_index + 1) = coeff |
| 935 |
2/2✓ Branch 27 → 22 taken 71100 times.
✓ Branch 27 → 28 taken 34184 times.
|
105284 | polar_index = polar_index + 1 |
| 936 | end do | ||
| 937 | end do | ||
| 938 | end do | ||
| 939 | end do | ||
| 940 | end do | ||
| 941 | |||
| 942 | TIMER_STOP(t_total) | ||
| 943 | end subroutine insert_polar_extractor_2forms | ||
| 944 | |||
| 945 | !> @brief Enrich 2-form basis for the polar regularity constraint using weak gradients of polar 3-forms | ||
| 946 | !> | ||
| 947 | !> @param[inout] this The MFormPolarEnriched object to initialize | ||
| 948 | 94 | subroutine insert_polar_extractor_2forms_enrich(this) | |
| 949 | implicit none | ||
| 950 | |||
| 951 | class(MFormPolarRegularity), intent(inout) :: this | ||
| 952 | |||
| 953 | integer :: i, j, polar_index, alpha, gamma, gamma_min, gamma_max, nr_polar_splines_z | ||
| 954 | integer :: i_min(3), i_max(3), j_min(3), j_max(3), s1 | ||
| 955 | |||
| 956 | TIMER_DECL(t_total) | ||
| 957 | TIMER_INIT_START(t_total, "MFormPolarRegularity::enrich2") | ||
| 958 | |||
| 959 | ! TODO include weak gradients of 3-forms for enrichment | ||
| 960 | |||
| 961 | gamma_min = 0 | ||
| 962 | |||
| 963 | ! The pushforward yields extra factor xp, and the 3rd component of the 2-form space has degree reduced by 1 | ||
| 964 | 94 | gamma_max = this%regularity - 2 | |
| 965 |
2/2✓ Branch 2 → 3 taken 32 times.
✓ Branch 2 → 4 taken 62 times.
|
94 | if (gamma_max < gamma_min) then |
| 966 | TIMER_STOP(t_total) | ||
| 967 | 32 | return | |
| 968 | end if | ||
| 969 | |||
| 970 | 62 | nr_polar_splines_z = get_nr_polar_zeroforms(gamma_max) | |
| 971 | |||
| 972 |
2/2✓ Branch 5 → 6 taken 186 times.
✓ Branch 5 → 8 taken 62 times.
|
248 | do s1 = 1, 3 |
| 973 | 248 | call this%index_list(s1)%get_bounds(i_min(s1), i_max(s1), j_min(s1), j_max(s1)) | |
| 974 | end do | ||
| 975 | |||
| 976 | !> Initialize the radial and angular B-splines (data members of this object) | ||
| 977 | call this%compute_1d_bsplines(this%space%tp_spaces(3)%spaces(1), this%space%tp_spaces(3)%spaces(2), & | ||
| 978 | 62 | gamma_min + 1, gamma_max + 1, radial_truncate=gamma_max + 1) | |
| 979 | |||
| 980 | 62 | call this%init_shmem_window(this%E2, this%shmem_E2, nr_polar_splines_z) | |
| 981 | s1 = 3 | ||
| 982 |
2/2✓ Branch 12 → 13 taken 1684 times.
✓ Branch 12 → 25 taken 62 times.
|
1746 | do j = j_min(s1), j_max(s1) |
| 983 |
2/2✓ Branch 14 → 15 taken 5648 times.
✓ Branch 14 → 23 taken 1684 times.
|
7394 | do i = i_min(s1), i_max(s1) |
| 984 | polar_index = 0 | ||
| 985 |
2/2✓ Branch 15 → 16 taken 11860 times.
✓ Branch 15 → 21 taken 5648 times.
|
19192 | do gamma = gamma_min, gamma_max |
| 986 |
1/2✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 11860 times.
|
17508 | do alpha = -gamma, gamma, 2 |
| 987 | this%E2(i, j, s1, polar_index + 1) = & | ||
| 988 |
2/2✓ Branch 18 → 18 taken 9178 times.
✓ Branch 18 → 19 taken 11860 times.
|
21038 | this%radial_bsplines(gamma + 1)%data(i) * this%angular_bsplines(alpha)%data(j) |
| 989 | 11860 | polar_index = polar_index + 1 | |
| 990 | end do | ||
| 991 | end do | ||
| 992 | end do | ||
| 993 | end do | ||
| 994 | |||
| 995 | TIMER_STOP(t_total) | ||
| 996 | end subroutine insert_polar_extractor_2forms_enrich | ||
| 997 | |||
| 998 | 221 | subroutine finalize_polar_extractor_2forms(this) | |
| 999 | implicit none | ||
| 1000 | |||
| 1001 | class(MFormPolarRegularity), intent(inout) :: this | ||
| 1002 | |||
| 1003 | integer :: s1 | ||
| 1004 | integer :: i, j, k, row_3d, polar_index, col, nr_polar_splines | ||
| 1005 | integer :: i_min(3), i_max(3), j_min(3), j_max(3), k_min(3), k_max(3), col0, nr_nonzeros(3) | ||
| 1006 | real(wp) :: coeff | ||
| 1007 | |||
| 1008 | TIMER_DECL(t_total) | ||
| 1009 | TIMER_INIT_START(t_total, "MFormPolarRegularity::finalize2") | ||
| 1010 | |||
| 1011 |
2/2✓ Branch 3 → 4 taken 663 times.
✓ Branch 3 → 6 taken 221 times.
|
884 | do s1 = 1, 3 |
| 1012 | 884 | call this%index_list(s1)%get_bounds(i_min(s1), i_max(s1), j_min(s1), j_max(s1), k_min(s1), k_max(s1)) | |
| 1013 | end do | ||
| 1014 | |||
| 1015 | 221 | nr_nonzeros(1) = size(this%index_list(1)) * size(this%E1, 4) | |
| 1016 | 221 | nr_nonzeros(2) = size(this%index_list(2)) * size(this%E1, 4) | |
| 1017 | 221 | nr_nonzeros(3) = size(this%index_list(3)) * size(this%E2, 4) | |
| 1018 | 221 | nr_polar_splines = size(this%E1, 4) + size(this%E2, 4) | |
| 1019 | |||
| 1020 |
2/2✓ Branch 8 → 9 taken 663 times.
✓ Branch 8 → 11 taken 221 times.
|
884 | do s1 = 1, 3 |
| 1021 | call this%extractor(s1)%init(nr_rows=size(this%index_list(s1)), nr_cols=nr_polar_splines * & | ||
| 1022 | 884 | & size(this%index_list(3), dir=3), nr_nonzero=nr_nonzeros(s1)) | |
| 1023 | end do | ||
| 1024 | |||
| 1025 |
2/2✓ Branch 12 → 13 taken 663 times.
✓ Branch 12 → 41 taken 221 times.
|
884 | do s1 = 1, 3 |
| 1026 |
2/2✓ Branch 13 → 14 taken 442 times.
✓ Branch 13 → 16 taken 221 times.
|
663 | if (s1 < 3) then |
| 1027 |
2/2✓ Branch 14 → 15 taken 220 times.
✓ Branch 14 → 39 taken 222 times.
|
442 | if (this%space%nr_components == 2) cycle |
| 1028 | 220 | nr_polar_splines = size(this%E1, 4) | |
| 1029 | col0 = 0 | ||
| 1030 | else | ||
| 1031 | 221 | nr_polar_splines = size(this%E2, 4) | |
| 1032 | 221 | col0 = size(this%E1, 4) * size(this%index_list(1), dir=3) | |
| 1033 | end if | ||
| 1034 | |||
| 1035 |
2/2✓ Branch 18 → 19 taken 441 times.
✓ Branch 18 → 20 taken 1803 times.
|
2465 | do k = k_min(s1), k_max(s1) |
| 1036 |
2/2✓ Branch 21 → 22 taken 48406 times.
✓ Branch 21 → 37 taken 1803 times.
|
50872 | do j = j_min(s1), j_max(s1) |
| 1037 |
2/2✓ Branch 23 → 24 taken 111838 times.
✓ Branch 23 → 35 taken 48406 times.
|
162047 | do i = i_min(s1), i_max(s1) |
| 1038 | 111838 | row_3d = this%index_list(s1)%get_local_index(i, j, k) | |
| 1039 | polar_index = 0 | ||
| 1040 | |||
| 1041 |
2/2✓ Branch 25 → 26 taken 625570 times.
✓ Branch 25 → 33 taken 111838 times.
|
785814 | do polar_index = 0, nr_polar_splines - 1 |
| 1042 | 625570 | col = polar_index + (k - k_min(s1)) * nr_polar_splines | |
| 1043 |
2/2✓ Branch 26 → 27 taken 554088 times.
✓ Branch 26 → 28 taken 71482 times.
|
625570 | if (s1 < 3) then |
| 1044 | 554088 | coeff = this%E1(i, j, s1, polar_index + 1) | |
| 1045 | else | ||
| 1046 | 71482 | coeff = this%E2(i, j, s1, polar_index + 1) | |
| 1047 | end if | ||
| 1048 | |||
| 1049 |
2/2✓ Branch 29 → 30 taken 307984 times.
✓ Branch 29 → 32 taken 317586 times.
|
737408 | if (abs(coeff) > COEFF_THR) call this%extractor(s1)%insert(row_3d, col0 + col, coeff) |
| 1050 | end do | ||
| 1051 | end do | ||
| 1052 | end do | ||
| 1053 | end do | ||
| 1054 | end do | ||
| 1055 | |||
| 1056 |
2/2✓ Branch 42 → 43 taken 663 times.
✓ Branch 42 → 45 taken 221 times.
|
884 | do s1 = 1, 3 |
| 1057 | 884 | call this%extractor(s1)%finalize() | |
| 1058 | end do | ||
| 1059 | |||
| 1060 | TIMER_STOP(t_total) | ||
| 1061 | 221 | end subroutine finalize_polar_extractor_2forms | |
| 1062 | |||
| 1063 | !> @brief Initialize the 3-form extractor for the polar regularity constraint | ||
| 1064 | !> | ||
| 1065 | !> @param[inout] this The MFormPolarRegularity object to initialize | ||
| 1066 | 110 | subroutine insert_polar_extractor_3forms(this) | |
| 1067 | use m_tensorprod, only: size | ||
| 1068 | use m_tensorprod_indices, only: TensorProdIndices | ||
| 1069 | implicit none | ||
| 1070 | |||
| 1071 | class(MFormPolarRegularity), intent(inout) :: this | ||
| 1072 | |||
| 1073 | type(TensorProdIndices) :: inds | ||
| 1074 | |||
| 1075 | integer :: nr_bspline_y, i_min, i_max, j_min, j_max, k_min, k_max | ||
| 1076 | |||
| 1077 | TIMER_DECL(t_total) | ||
| 1078 | TIMER_INIT_START(t_total, "MFormPolarRegularity::insert3") | ||
| 1079 | |||
| 1080 | 110 | nr_bspline_y = size(this%space%tp_spaces(1), 2) | |
| 1081 | |||
| 1082 | !> Initialize the indices on which this extractor acts | ||
| 1083 | 110 | i_min = max(0, this%space%tp_spaces(1)%rank_resp_bspline%i0) | |
| 1084 | ! NOTE i_max is at least 0, because of the pushforward factor 1/xp | ||
| 1085 | 110 | i_max = min(this%regularity, this%space%tp_spaces(1)%spaces(1)%degree, this%space%tp_spaces(1)%rank_resp_bspline%i1) | |
| 1086 | 110 | j_min = this%space%tp_spaces(1)%rank_resp_bspline%j0; j_max = this%space%tp_spaces(1)%rank_resp_bspline%j1 | |
| 1087 | 110 | k_min = this%space%tp_spaces(1)%rank_resp_bspline%k0; k_max = this%space%tp_spaces(1)%rank_resp_bspline%k1 | |
| 1088 | |||
| 1089 | 110 | call inds%init(i_min, i_max, j_min, j_max, k_min, k_max, this%space%tp_spaces(1)%spaces) | |
| 1090 | 110 | call this%index_list(1)%init(inds, this%space%tp_spaces(1)%rank_resp_bspline, this%space%tp_spaces(1)%rank_l0) | |
| 1091 | |||
| 1092 | 110 | call this%init_shmem_window(this%E1, this%shmem_E1, 0) | |
| 1093 | |||
| 1094 | TIMER_STOP(t_total) | ||
| 1095 | 110 | end subroutine insert_polar_extractor_3forms | |
| 1096 | |||
| 1097 | !> @brief Enrich 3-form basis for the polar regularity constraint using regular 3-forms | ||
| 1098 | !> | ||
| 1099 | !> @param[inout] this The MFormPolarEnriched object to initialize | ||
| 1100 | 49 | subroutine insert_polar_extractor_3forms_enrich(this) | |
| 1101 | implicit none | ||
| 1102 | |||
| 1103 | class(MFormPolarRegularity), intent(inout) :: this | ||
| 1104 | |||
| 1105 | integer :: i, j, polar_index, alpha, gamma, gamma_min, gamma_max, nr_polar_splines_z | ||
| 1106 | integer :: i_min, i_max, j_min, j_max | ||
| 1107 | |||
| 1108 | TIMER_DECL(t_total) | ||
| 1109 | TIMER_INIT_START(t_total, "MFormPolarRegularity::enrich3") | ||
| 1110 | |||
| 1111 | gamma_min = 0 | ||
| 1112 | |||
| 1113 | ! The pushforward yields extra factor xp, and the 3rd component of the 2-form space has degree reduced by 1 | ||
| 1114 | 49 | gamma_max = this%regularity - 2 | |
| 1115 |
2/2✓ Branch 2 → 3 taken 17 times.
✓ Branch 2 → 4 taken 32 times.
|
49 | if (gamma_max < gamma_min) then |
| 1116 | TIMER_STOP(t_total) | ||
| 1117 | 17 | return | |
| 1118 | end if | ||
| 1119 | |||
| 1120 | 32 | nr_polar_splines_z = get_nr_polar_zeroforms(gamma_max) | |
| 1121 | |||
| 1122 | 32 | call this%index_list(1)%get_bounds(i_min, i_max, j_min, j_max) | |
| 1123 | |||
| 1124 | !> Initialize the radial and angular B-splines (data members of this object) | ||
| 1125 | call this%compute_1d_bsplines(this%space%tp_spaces(1)%spaces(1), this%space%tp_spaces(1)%spaces(2), & | ||
| 1126 | 32 | gamma_min + 1, gamma_max + 1, radial_truncate=gamma_max + 1) | |
| 1127 | |||
| 1128 | 32 | call this%init_shmem_window(this%E1, this%shmem_E1, nr_polar_splines_z) | |
| 1129 | |||
| 1130 |
2/2✓ Branch 8 → 9 taken 848 times.
✓ Branch 8 → 21 taken 32 times.
|
880 | do j = j_min, j_max |
| 1131 |
2/2✓ Branch 10 → 11 taken 2788 times.
✓ Branch 10 → 19 taken 848 times.
|
3668 | do i = i_min, i_max |
| 1132 | polar_index = 0 | ||
| 1133 |
2/2✓ Branch 11 → 12 taken 5624 times.
✓ Branch 11 → 17 taken 2788 times.
|
9260 | do gamma = gamma_min, gamma_max |
| 1134 |
1/2✗ Branch 12 → 13 not taken.
✓ Branch 12 → 14 taken 5624 times.
|
8412 | do alpha = -gamma, gamma, 2 |
| 1135 | this%E1(i, j, 1, polar_index + 1) = & | ||
| 1136 |
2/2✓ Branch 14 → 14 taken 4004 times.
✓ Branch 14 → 15 taken 5624 times.
|
9628 | this%radial_bsplines(gamma + 1)%data(i) * this%angular_bsplines(alpha)%data(j) |
| 1137 | 5624 | polar_index = polar_index + 1 | |
| 1138 | end do | ||
| 1139 | end do | ||
| 1140 | end do | ||
| 1141 | end do | ||
| 1142 | |||
| 1143 | TIMER_STOP(t_total) | ||
| 1144 | end subroutine insert_polar_extractor_3forms_enrich | ||
| 1145 | |||
| 1146 | 110 | subroutine finalize_polar_extractor_3forms(this) | |
| 1147 | implicit none | ||
| 1148 | |||
| 1149 | class(MFormPolarRegularity), intent(inout) :: this | ||
| 1150 | |||
| 1151 | integer :: i, j, k, row_2d, row_3d, polar_index, col, nr_polar_splines | ||
| 1152 | integer :: i_min, i_max, j_min, j_max, k_min, k_max, nr_nonzeros | ||
| 1153 | real(wp) :: coeff | ||
| 1154 | |||
| 1155 | TIMER_DECL(t_total) | ||
| 1156 | TIMER_INIT_START(t_total, "MFormPolarRegularity::finalize3") | ||
| 1157 | |||
| 1158 | 110 | call this%index_list(1)%get_bounds(i_min, i_max, j_min, j_max, k_min, k_max) | |
| 1159 | |||
| 1160 | 110 | nr_nonzeros = size(this%index_list(1)) * size(this%E1, 4) | |
| 1161 | nr_polar_splines = size(this%E1, 4) | ||
| 1162 | |||
| 1163 | call this%extractor(1)%init(nr_rows=size(this%index_list(1)), nr_cols=nr_polar_splines * & | ||
| 1164 | 110 | & size(this%index_list(1), dir=3), nr_nonzero=nr_nonzeros) | |
| 1165 | |||
| 1166 |
2/2✓ Branch 5 → 6 taken 564 times.
✓ Branch 5 → 22 taken 110 times.
|
674 | do k = k_min, k_max |
| 1167 |
2/2✓ Branch 7 → 8 taken 15008 times.
✓ Branch 7 → 20 taken 564 times.
|
15682 | do j = j_min, j_max |
| 1168 |
2/2✓ Branch 9 → 10 taken 35680 times.
✓ Branch 9 → 18 taken 15008 times.
|
51252 | do i = i_min, i_max |
| 1169 | 35680 | row_3d = this%index_list(1)%get_local_index(i, j, k) | |
| 1170 | row_2d = this%index_list(1)%get_local_index(i, j) | ||
| 1171 | polar_index = 0 | ||
| 1172 | |||
| 1173 |
2/2✓ Branch 12 → 13 taken 60072 times.
✓ Branch 12 → 16 taken 35680 times.
|
110760 | do polar_index = 0, nr_polar_splines - 1 |
| 1174 | 60072 | col = polar_index + (k - k_min) * nr_polar_splines | |
| 1175 | 60072 | coeff = this%E1(i, j, 1, polar_index + 1) | |
| 1176 | |||
| 1177 |
2/2✓ Branch 13 → 11 taken 33340 times.
✓ Branch 13 → 14 taken 26732 times.
|
95752 | if (abs(coeff) > COEFF_THR) call this%extractor(1)%insert(row_3d, col, coeff) |
| 1178 | end do | ||
| 1179 | end do | ||
| 1180 | end do | ||
| 1181 | end do | ||
| 1182 | |||
| 1183 | 110 | call this%extractor(1)%finalize() | |
| 1184 | |||
| 1185 | TIMER_STOP(t_total) | ||
| 1186 | 110 | end subroutine finalize_polar_extractor_3forms | |
| 1187 | |||
| 1188 | !> @brief Compute the 1D B-spline coefficients for the polar regularity constraint | ||
| 1189 | !> | ||
| 1190 | !> @param[inout] this The MFormPolarRegularity object for which to compute the B-splines | ||
| 1191 | !> @param[in] bspline_x The B-spline space in the x-direction | ||
| 1192 | !> @param[in] bspline_y The B-spline space in the y-direction | ||
| 1193 | !> @param[in] gamma_min The minimum radial degree | ||
| 1194 | !> @param[in] gamma_max The maximum radial degree | ||
| 1195 | !> @param[in] radial_truncate _(optional)_ The index at which to truncate the radial functions, 0:radial_truncate is kept | ||
| 1196 | 1010 | subroutine compute_1d_bsplines(this, bspline_x, bspline_y, gamma_min, gamma_max, radial_truncate) | |
| 1197 | use m_bspline, only: BSplineSpace, size | ||
| 1198 | use m_bspline_linalg, only: l2_projection, biorthogonalize | ||
| 1199 | use m_common, only: user_function_1d_interface | ||
| 1200 | |||
| 1201 | implicit none | ||
| 1202 | |||
| 1203 | class(MFormPolarRegularity), intent(inout) :: this | ||
| 1204 | type(BSplineSpace), intent(in) :: bspline_x, bspline_y | ||
| 1205 | integer, intent(in) :: gamma_min, gamma_max | ||
| 1206 | integer, optional, intent(in) :: radial_truncate | ||
| 1207 | |||
| 1208 | integer :: alpha, gamma, nr_bsplines_y | ||
| 1209 | |||
| 1210 | TIMER_DECL(t_total) | ||
| 1211 | TIMER_INIT_START(t_total, "MFormPolarRegularity::compute_1d_bsplines") | ||
| 1212 | |||
| 1213 |
2/4✓ Branch 2 → 3 taken 1010 times.
✗ Branch 2 → 4 not taken.
✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 1010 times.
|
1010 | if (gamma_min < 0 .or. gamma_max > bspline_x%degree) then |
| 1214 | error stop "MFormPolarRegularity::compute_1d_bsplines: regularity must be in [0, degree_x] where degree_x is "// & | ||
| 1215 | ✗ | & " the degree of the B-spline basis in the x-direction" | |
| 1216 | end if | ||
| 1217 | |||
| 1218 |
1/2✗ Branch 5 → 6 not taken.
✓ Branch 5 → 7 taken 1010 times.
|
1010 | if (2 * gamma_max + 1 >= size(bspline_y)) then |
| 1219 | error stop "MFormPolarRegularity::compute_1d_bsplines: 2*regularity + 1 must be less than the number of basis"// & | ||
| 1220 | ✗ | & " functions in the y-direction" | |
| 1221 | end if | ||
| 1222 | |||
| 1223 |
6/8✓ Branch 7 → 8 taken 102 times.
✓ Branch 7 → 16 taken 908 times.
✓ Branch 9 → 10 taken 392 times.
✓ Branch 9 → 13 taken 102 times.
✓ Branch 10 → 11 taken 392 times.
✗ Branch 10 → 12 not taken.
✗ Branch 13 → 14 not taken.
✓ Branch 13 → 15 taken 102 times.
|
1402 | if (allocated(this%radial_bsplines)) deallocate (this%radial_bsplines) |
| 1224 |
6/8✓ Branch 16 → 17 taken 102 times.
✓ Branch 16 → 25 taken 908 times.
✓ Branch 18 → 19 taken 682 times.
✓ Branch 18 → 22 taken 102 times.
✓ Branch 19 → 20 taken 682 times.
✗ Branch 19 → 21 not taken.
✗ Branch 22 → 23 not taken.
✓ Branch 22 → 24 taken 102 times.
|
1692 | if (allocated(this%angular_bsplines)) deallocate (this%angular_bsplines) |
| 1225 | |||
| 1226 |
9/16✗ Branch 25 → 26 not taken.
✓ Branch 25 → 27 taken 1010 times.
✓ Branch 27 → 26 taken 1010 times.
✗ Branch 27 → 28 not taken.
✓ Branch 28 → 29 taken 1010 times.
✗ Branch 28 → 30 not taken.
✓ Branch 30 → 31 taken 1010 times.
✗ Branch 30 → 32 not taken.
✗ Branch 32 → 33 not taken.
✓ Branch 32 → 34 taken 1010 times.
✗ Branch 34 → 35 not taken.
✓ Branch 34 → 36 taken 1010 times.
✗ Branch 36 → 37 not taken.
✓ Branch 36 → 38 taken 1010 times.
✓ Branch 39 → 40 taken 2281 times.
✓ Branch 39 → 41 taken 1010 times.
|
5311 | allocate (this%radial_bsplines(gamma_min:gamma_max)) |
| 1227 |
9/16✗ Branch 41 → 42 not taken.
✓ Branch 41 → 43 taken 1010 times.
✓ Branch 43 → 42 taken 1010 times.
✗ Branch 43 → 44 not taken.
✓ Branch 44 → 45 taken 1010 times.
✗ Branch 44 → 46 not taken.
✓ Branch 46 → 47 taken 1010 times.
✗ Branch 46 → 48 not taken.
✗ Branch 48 → 49 not taken.
✓ Branch 48 → 50 taken 1010 times.
✗ Branch 50 → 51 not taken.
✓ Branch 50 → 52 taken 1010 times.
✗ Branch 52 → 53 not taken.
✓ Branch 52 → 54 taken 1010 times.
✓ Branch 55 → 56 taken 3880 times.
✓ Branch 55 → 57 taken 1010 times.
|
6910 | allocate (this%angular_bsplines(-gamma_max:gamma_max)) |
| 1228 | |||
| 1229 | ! Compute L2 approximation of the radial and angular functions which make up the Zernike basis | ||
| 1230 |
2/2✓ Branch 58 → 59 taken 2281 times.
✓ Branch 58 → 84 taken 1010 times.
|
3291 | do gamma = gamma_min, gamma_max |
| 1231 |
2/2✓ Branch 59 → 60 taken 846 times.
✓ Branch 59 → 64 taken 1435 times.
|
2281 | if (gamma == 0) then |
| 1232 | 846 | call this%radial_bsplines(gamma)%init(bspline_x) | |
| 1233 |
2/2✓ Branch 62 → 63 taken 17286 times.
✓ Branch 62 → 68 taken 846 times.
|
18132 | this%radial_bsplines(gamma)%data(0:) = 1._wp |
| 1234 | else | ||
| 1235 | 1435 | call l2_projection(this%radial_bsplines(gamma), bspline_x, radial_fun) | |
| 1236 | ! Analytically the first gamma coefficients are zero | ||
| 1237 |
2/2✓ Branch 66 → 67 taken 2667 times.
✓ Branch 66 → 68 taken 1435 times.
|
4102 | this%radial_bsplines(gamma)%data(0:gamma - 1) = 0._wp |
| 1238 | end if | ||
| 1239 |
3/4✓ Branch 68 → 69 taken 2281 times.
✗ Branch 68 → 72 not taken.
✓ Branch 70 → 71 taken 37926 times.
✓ Branch 70 → 72 taken 2281 times.
|
40207 | if (present(radial_truncate)) this%radial_bsplines(gamma)%data(radial_truncate + 1:) = 0._wp |
| 1240 | |||
| 1241 | ! Optional: we normalize the coefficients for numerical stability (merely a change of basis) | ||
| 1242 | this%radial_bsplines(gamma)%data(0:) = this%radial_bsplines(gamma)%data(0:) / & | ||
| 1243 |
6/10✓ Branch 73 → 74 taken 2281 times.
✗ Branch 73 → 76 not taken.
✗ Branch 74 → 75 not taken.
✓ Branch 74 → 78 taken 2281 times.
✗ Branch 76 → 77 not taken.
✗ Branch 76 → 78 not taken.
✓ Branch 79 → 80 taken 45235 times.
✓ Branch 79 → 81 taken 2281 times.
✓ Branch 81 → 82 taken 45235 times.
✓ Branch 81 → 83 taken 2281 times.
|
96042 | maxval(abs(this%radial_bsplines(gamma)%data(0:))) |
| 1244 | end do | ||
| 1245 | |||
| 1246 | 1010 | nr_bsplines_y = size(bspline_y) | |
| 1247 |
2/2✓ Branch 86 → 87 taken 3880 times.
✓ Branch 86 → 100 taken 1010 times.
|
4890 | do alpha = -gamma_max, gamma_max |
| 1248 |
2/2✓ Branch 87 → 88 taken 1010 times.
✓ Branch 87 → 92 taken 2870 times.
|
4890 | if (alpha == 0) then |
| 1249 | 1010 | call this%angular_bsplines(alpha)%init(bspline_y) | |
| 1250 |
2/2✓ Branch 90 → 91 taken 32398 times.
✓ Branch 90 → 99 taken 1010 times.
|
33408 | this%angular_bsplines(alpha)%data(:) = 1._wp |
| 1251 | else | ||
| 1252 | 2870 | call l2_projection(this%angular_bsplines(alpha), bspline_y, angular_fun) | |
| 1253 | ! Force zero mean | ||
| 1254 | this%angular_bsplines(alpha)%data(:) = this%angular_bsplines(alpha)%data(:) - & | ||
| 1255 |
4/4✓ Branch 94 → 95 taken 80032 times.
✓ Branch 94 → 96 taken 2870 times.
✓ Branch 97 → 98 taken 90028 times.
✓ Branch 97 → 99 taken 2870 times.
|
172930 | sum(this%angular_bsplines(alpha)%data(0:nr_bsplines_y - 1)) / nr_bsplines_y |
| 1256 | end if | ||
| 1257 | end do | ||
| 1258 | |||
| 1259 | TIMER_STOP(t_total) | ||
| 1260 | contains | ||
| 1261 | 100600 | pure real(wp) function radial_fun(xp) result(ans) | |
| 1262 | implicit none | ||
| 1263 | real(wp), intent(in) :: xp | ||
| 1264 | |||
| 1265 | 100600 | ans = xp**gamma | |
| 1266 | 100600 | end function radial_fun | |
| 1267 | |||
| 1268 | 370272 | pure real(wp) function angular_fun(yp) result(ans) | |
| 1269 | use m_common, only: pi | ||
| 1270 | |||
| 1271 | implicit none | ||
| 1272 | real(wp), intent(in) :: yp | ||
| 1273 | |||
| 1274 |
1/2✓ Branch 2 → 3 taken 370272 times.
✗ Branch 2 → 6 not taken.
|
370272 | if (alpha == 0) then |
| 1275 | ans = 1._wp | ||
| 1276 |
2/2✓ Branch 3 → 4 taken 185136 times.
✓ Branch 3 → 5 taken 185136 times.
|
370272 | else if (alpha < 0) then |
| 1277 | 185136 | ans = sin(-2 * pi * alpha * yp) | |
| 1278 | else | ||
| 1279 | 185136 | ans = cos(2 * pi * alpha * yp) | |
| 1280 | end if | ||
| 1281 | |||
| 1282 | 370272 | end function angular_fun | |
| 1283 | end subroutine compute_1d_bsplines | ||
| 1284 | |||
| 1285 | !> @brief Destroy the MFormPolarRegularity object | ||
| 1286 | !> | ||
| 1287 | !> @param[inout] this The MFormPolarRegularity object to destroy | ||
| 1288 | 1450 | subroutine destroy_mform_polarregularity(this) | |
| 1289 | implicit none | ||
| 1290 | class(MFormPolarRegularity), intent(inout) :: this | ||
| 1291 | integer :: s1, regularity | ||
| 1292 | |||
| 1293 |
2/2✓ Branch 2 → 3 taken 393 times.
✓ Branch 2 → 21 taken 1057 times.
|
1450 | if (allocated(this%radial_bsplines)) then |
| 1294 |
4/6✓ Branch 3 → 4 taken 393 times.
✗ Branch 3 → 5 not taken.
✓ Branch 5 → 6 taken 393 times.
✗ Branch 5 → 7 not taken.
✓ Branch 8 → 9 taken 855 times.
✓ Branch 8 → 11 taken 393 times.
|
2034 | do s1 = lbound(this%radial_bsplines, 1), ubound(this%radial_bsplines, 1) |
| 1295 | 1248 | call this%radial_bsplines(s1)%destroy() | |
| 1296 | end do | ||
| 1297 |
5/8✓ Branch 12 → 13 taken 393 times.
✗ Branch 12 → 18 not taken.
✓ Branch 14 → 15 taken 855 times.
✓ Branch 14 → 18 taken 393 times.
✗ Branch 15 → 16 not taken.
✓ Branch 15 → 17 taken 855 times.
✗ Branch 18 → 19 not taken.
✓ Branch 18 → 20 taken 393 times.
|
1248 | deallocate (this%radial_bsplines) |
| 1298 | end if | ||
| 1299 | |||
| 1300 |
2/2✓ Branch 21 → 22 taken 393 times.
✓ Branch 21 → 40 taken 1057 times.
|
1450 | if (allocated(this%angular_bsplines)) then |
| 1301 |
4/6✓ Branch 22 → 23 taken 393 times.
✗ Branch 22 → 24 not taken.
✓ Branch 24 → 25 taken 393 times.
✗ Branch 24 → 26 not taken.
✓ Branch 27 → 28 taken 1425 times.
✓ Branch 27 → 30 taken 393 times.
|
2604 | do s1 = lbound(this%angular_bsplines, 1), ubound(this%angular_bsplines, 1) |
| 1302 | 1818 | call this%angular_bsplines(s1)%destroy() | |
| 1303 | end do | ||
| 1304 |
5/8✓ Branch 31 → 32 taken 393 times.
✗ Branch 31 → 37 not taken.
✓ Branch 33 → 34 taken 1425 times.
✓ Branch 33 → 37 taken 393 times.
✗ Branch 34 → 35 not taken.
✓ Branch 34 → 36 taken 1425 times.
✗ Branch 37 → 38 not taken.
✓ Branch 37 → 39 taken 393 times.
|
1818 | deallocate (this%angular_bsplines) |
| 1305 | end if | ||
| 1306 | |||
| 1307 | 1450 | call this%shmem_E1%destroy() | |
| 1308 | 1450 | call this%shmem_E2%destroy() | |
| 1309 |
2/2✓ Branch 42 → 43 taken 881 times.
✓ Branch 42 → 44 taken 569 times.
|
1450 | if (associated(this%E1)) nullify (this%E1) |
| 1310 |
2/2✓ Branch 44 → 45 taken 937 times.
✓ Branch 44 → 46 taken 513 times.
|
1450 | if (associated(this%E2)) nullify (this%E2) |
| 1311 | |||
| 1312 | 1450 | this%my_rank_window = -1 | |
| 1313 | 1450 | this%comm_window = -1 | |
| 1314 | 1450 | end subroutine destroy_mform_polarregularity | |
| 1315 | |||
| 1316 | pure integer function get_nr_polar_zeroforms(regularity) result(nr) | ||
| 1317 | implicit none | ||
| 1318 | integer, intent(in) :: regularity | ||
| 1319 | |||
| 1320 |
3/4✓ Branch 2 → 3 taken 217 times.
✓ Branch 2 → 4 taken 116 times.
✗ Branch 115 → 116 not taken.
✓ Branch 115 → 117 taken 70 times.
|
1010 | nr = (regularity + 1) * (regularity + 2) / 2 |
| 1321 | end function get_nr_polar_zeroforms | ||
| 1322 | |||
| 1323 |
74/122__m_mform_constraint_polar_MOD___copy_m_mform_constraint_polar_Mformpolarregularity:
✓ Branch 2 → 3 taken 365 times.
✗ Branch 2 → 70 not taken.
✓ Branch 3 → 4 taken 60 times.
✓ Branch 3 → 5 taken 305 times.
✓ Branch 6 → 7 taken 60 times.
✓ Branch 6 → 22 taken 305 times.
✓ Branch 8 → 9 taken 60 times.
✓ Branch 8 → 23 taken 60 times.
✓ Branch 9 → 10 taken 60 times.
✗ Branch 9 → 11 not taken.
✓ Branch 12 → 13 taken 60 times.
✗ Branch 12 → 14 not taken.
✓ Branch 15 → 16 taken 60 times.
✗ Branch 15 → 17 not taken.
✓ Branch 18 → 19 taken 60 times.
✗ Branch 18 → 20 not taken.
✗ Branch 23 → 24 not taken.
✓ Branch 23 → 36 taken 365 times.
✗ Branch 25 → 26 not taken.
✗ Branch 25 → 37 not taken.
✗ Branch 26 → 27 not taken.
✗ Branch 26 → 28 not taken.
✗ Branch 29 → 30 not taken.
✗ Branch 29 → 31 not taken.
✗ Branch 32 → 33 not taken.
✗ Branch 32 → 34 not taken.
✓ Branch 37 → 38 taken 60 times.
✓ Branch 37 → 50 taken 305 times.
✓ Branch 39 → 40 taken 60 times.
✓ Branch 39 → 51 taken 60 times.
✓ Branch 40 → 41 taken 60 times.
✗ Branch 40 → 42 not taken.
✓ Branch 43 → 44 taken 60 times.
✗ Branch 43 → 45 not taken.
✓ Branch 46 → 47 taken 60 times.
✗ Branch 46 → 48 not taken.
✓ Branch 51 → 52 taken 60 times.
✓ Branch 51 → 58 taken 305 times.
✓ Branch 53 → 54 taken 165 times.
✓ Branch 53 → 59 taken 60 times.
✓ Branch 54 → 55 taken 165 times.
✗ Branch 54 → 56 not taken.
✓ Branch 59 → 60 taken 60 times.
✓ Branch 59 → 66 taken 305 times.
✓ Branch 61 → 62 taken 270 times.
✓ Branch 61 → 67 taken 60 times.
✓ Branch 62 → 63 taken 270 times.
✗ Branch 62 → 64 not taken.
✓ Branch 67 → 68 taken 60 times.
✓ Branch 67 → 69 taken 305 times.
__m_mform_constraint_polar_MOD___final_m_mform_constraint_polar_Mformpolarregularity:
✗ Branch 3 → 4 not taken.
✓ Branch 3 → 7 taken 365 times.
✗ Branch 4 → 5 not taken.
✗ Branch 4 → 6 not taken.
✓ Branch 8 → 9 taken 365 times.
✓ Branch 8 → 76 taken 365 times.
✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 365 times.
✓ Branch 12 → 13 taken 365 times.
✗ Branch 12 → 15 not taken.
✓ Branch 13 → 14 taken 365 times.
✗ Branch 13 → 16 not taken.
✗ Branch 15 → 16 not taken.
✗ Branch 15 → 55 not taken.
✓ Branch 16 → 17 taken 2 times.
✓ Branch 16 → 28 taken 363 times.
✓ Branch 18 → 19 taken 2 times.
✓ Branch 18 → 28 taken 2 times.
✓ Branch 19 → 20 taken 2 times.
✗ Branch 19 → 21 not taken.
✓ Branch 21 → 22 taken 2 times.
✗ Branch 21 → 23 not taken.
✓ Branch 23 → 24 taken 2 times.
✗ Branch 23 → 25 not taken.
✓ Branch 25 → 26 taken 2 times.
✗ Branch 25 → 27 not taken.
✓ Branch 28 → 29 taken 2 times.
✓ Branch 28 → 30 taken 363 times.
✗ Branch 30 → 31 not taken.
✓ Branch 30 → 40 taken 365 times.
✗ Branch 32 → 33 not taken.
✗ Branch 32 → 40 not taken.
✗ Branch 33 → 34 not taken.
✗ Branch 33 → 35 not taken.
✗ Branch 35 → 36 not taken.
✗ Branch 35 → 37 not taken.
✗ Branch 37 → 38 not taken.
✗ Branch 37 → 39 not taken.
✗ Branch 40 → 41 not taken.
✓ Branch 40 → 42 taken 365 times.
✓ Branch 42 → 43 taken 2 times.
✓ Branch 42 → 52 taken 363 times.
✓ Branch 44 → 45 taken 2 times.
✓ Branch 44 → 52 taken 2 times.
✓ Branch 45 → 46 taken 2 times.
✗ Branch 45 → 47 not taken.
✓ Branch 47 → 48 taken 2 times.
✗ Branch 47 → 49 not taken.
✓ Branch 49 → 50 taken 2 times.
✗ Branch 49 → 51 not taken.
✓ Branch 52 → 53 taken 2 times.
✓ Branch 52 → 54 taken 363 times.
✓ Branch 55 → 56 taken 365 times.
✗ Branch 55 → 73 not taken.
✓ Branch 56 → 57 taken 1 time.
✓ Branch 56 → 62 taken 364 times.
✓ Branch 58 → 59 taken 6 times.
✓ Branch 58 → 62 taken 1 time.
✓ Branch 59 → 60 taken 6 times.
✗ Branch 59 → 61 not taken.
✓ Branch 62 → 63 taken 1 time.
✓ Branch 62 → 64 taken 364 times.
✓ Branch 64 → 65 taken 1 time.
✓ Branch 64 → 70 taken 364 times.
✓ Branch 66 → 67 taken 11 times.
✓ Branch 66 → 70 taken 1 time.
✓ Branch 67 → 68 taken 11 times.
✗ Branch 67 → 69 not taken.
✓ Branch 70 → 71 taken 1 time.
✓ Branch 70 → 72 taken 364 times.
✓ Branch 73 → 74 taken 365 times.
✗ Branch 73 → 75 not taken.
|
2036 | end module m_mform_constraint_polar |
| 1324 |