tensorprod/m_tensorprod_linalg.F90
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | !> @brief The tensor product linear algebra module | ||
| 2 | !> | ||
| 3 | !> This module provides an inner product, L2 projection, integration, and L2 error computation for tensor product spaces | ||
| 4 | module m_tensorprod_linalg | ||
| 5 | use m_tensorprod_basis, only: TensorProdFun, TensorProdSpace | ||
| 6 | use m_common, only: wp, user_function_3d_interface | ||
| 7 | #include "petsc.fi" | ||
| 8 | |||
| 9 | private | ||
| 10 | public :: inner_product, l2_projection, integrate, l2_error | ||
| 11 | |||
| 12 | !> @brief Compute the (L2) inner product of all elements of a tensor product B-spline space with a user-defined function | ||
| 13 | interface inner_product | ||
| 14 | procedure inner_product_3d | ||
| 15 | end interface inner_product | ||
| 16 | |||
| 17 | !> @brief Compute the L2 projection of a user-defined function onto a tensor product B-spline space | ||
| 18 | interface l2_projection | ||
| 19 | procedure l2_projection_3d | ||
| 20 | end interface l2_projection | ||
| 21 | |||
| 22 | !> @brief Integrate a user-defined function | ||
| 23 | interface integrate | ||
| 24 | procedure integrate_3d | ||
| 25 | end interface integrate | ||
| 26 | |||
| 27 | !> @brief Compute the L2 norm of the difference between a tensor product B-spline function and a user-defined function | ||
| 28 | interface l2_error | ||
| 29 | procedure l2_error_3d | ||
| 30 | end interface l2_error | ||
| 31 | |||
| 32 | contains | ||
| 33 | |||
| 34 | !> @brief Compute the (L2) inner product of all elements of a tensor product B-spline space with a user-defined function | ||
| 35 | !> | ||
| 36 | !> @param[inout] v The resulting tensor product function containing the inner product values | ||
| 37 | !> @param[in] tp_space The tensor product space for which the inner product is computed | ||
| 38 | !> @param[in] userfun The user-defined function to compute the inner product with | ||
| 39 | !> @param[in] n_quad_extra (optional) The number of extra quadrature points on top of the number needed for exact integration of a | ||
| 40 | !> product of B-splines | ||
| 41 | !> | ||
| 42 | !> @note The resulting tensor product function `v` is initialized only for the response B-spline indices of the tensor product | ||
| 43 | !> space (if needed, call v%distribute() to distribute the data) | ||
| 44 | 2166 | subroutine inner_product_3d(v, tp_space, userfun, n_quad_extra) | |
| 45 | use m_bspline_basis, only: get_internal_and_actual_inds | ||
| 46 | use m_bspline_quadrature, only: BSplineQuadrature, n_quad_exact | ||
| 47 | use m_tensorprod_basis, only: TensorProdSpace, size, actv_interval_inds | ||
| 48 | use m_tensorprod_quadrature, only: TensorProdQuadrature, evaluate_at_nodes | ||
| 49 | implicit none | ||
| 50 | |||
| 51 | type(TensorProdFun), intent(inout) :: v | ||
| 52 | type(TensorProdSpace), intent(in) :: tp_space | ||
| 53 | procedure(user_function_3d_interface) :: userfun | ||
| 54 | integer, intent(in), optional :: n_quad_extra | ||
| 55 | |||
| 56 | integer :: i, j, k, i_internal, j_internal, k_internal | ||
| 57 | integer :: ii_index, jj_index, kk_index | ||
| 58 | integer :: ii, jj, kk, i_minus_ii, j_minus_jj, k_minus_kk | ||
| 59 | integer :: degree_x, degree_y, degree_z, n_quad1, n_quad2, n_quad3 | ||
| 60 | integer :: nr_intervals_x, nr_intervals_y, nr_intervals_z, qdx, qdy, qdz | ||
| 61 | integer :: n_quad_extra_ | ||
| 62 | logical :: bspline_mirrored_x, bspline_mirrored_y, bspline_mirrored_z | ||
| 63 | real(wp) :: val | ||
| 64 | |||
| 65 | 2166 | real(wp), allocatable :: userfun_vals(:, :, :) | |
| 66 | real(wp), allocatable :: bspline_vals_x(:), bspline_vals_y(:), bspline_vals_z(:) | ||
| 67 | real(wp), allocatable :: partial_sum_z(:, :), partial_sum_yz(:) | ||
| 68 | |||
| 69 | 2166 | integer, allocatable :: actv_interval_inds_x(:), actv_interval_inds_y(:), actv_interval_inds_z(:) | |
| 70 | |||
| 71 |
2/2✓ Branch 3 → 4 taken 6498 times.
✓ Branch 3 → 5 taken 2166 times.
|
8664 | type(TensorProdQuadrature) :: tp_space_quad |
| 72 | |||
| 73 |
2/2✓ Branch 5 → 6 taken 2121 times.
✓ Branch 5 → 7 taken 45 times.
|
2166 | if (present(n_quad_extra)) then |
| 74 | 2121 | n_quad_extra_ = n_quad_extra | |
| 75 | else | ||
| 76 | n_quad_extra_ = 0 | ||
| 77 | end if | ||
| 78 | |||
| 79 | 2166 | call actv_interval_inds(actv_interval_inds_x, actv_interval_inds_y, actv_interval_inds_z, tp_space) | |
| 80 | |||
| 81 | 2166 | nr_intervals_x = size(actv_interval_inds_x) | |
| 82 | 2166 | nr_intervals_y = size(actv_interval_inds_y) | |
| 83 | 2166 | nr_intervals_z = size(actv_interval_inds_z) | |
| 84 | |||
| 85 | 2166 | n_quad1 = n_quad_exact(tp_space%spaces(1), tp_space%spaces(1)) + n_quad_extra_ | |
| 86 | 2166 | n_quad2 = n_quad_exact(tp_space%spaces(2), tp_space%spaces(2)) + n_quad_extra_ | |
| 87 | 2166 | n_quad3 = n_quad_exact(tp_space%spaces(3), tp_space%spaces(3)) + n_quad_extra_ | |
| 88 | |||
| 89 | tp_space_quad = TensorProdQuadrature(BSplineQuadrature(tp_space%spaces(1), n_quad1), & | ||
| 90 | BSplineQuadrature(tp_space%spaces(2), n_quad2), & | ||
| 91 |
29/50✓ Branch 15 → 16 taken 6498 times.
✓ Branch 15 → 23 taken 2166 times.
✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 6498 times.
✗ Branch 18 → 19 not taken.
✓ Branch 18 → 20 taken 6498 times.
✗ Branch 20 → 21 not taken.
✓ Branch 20 → 22 taken 6498 times.
✓ Branch 23 → 24 taken 2166 times.
✗ Branch 23 → 25 not taken.
✓ Branch 25 → 26 taken 2166 times.
✗ Branch 25 → 27 not taken.
✓ Branch 27 → 28 taken 2166 times.
✗ Branch 27 → 29 not taken.
✓ Branch 29 → 30 taken 2166 times.
✗ Branch 29 → 31 not taken.
✓ Branch 31 → 32 taken 2166 times.
✗ Branch 31 → 33 not taken.
✓ Branch 33 → 34 taken 2166 times.
✗ Branch 33 → 35 not taken.
✓ Branch 35 → 36 taken 2166 times.
✗ Branch 35 → 37 not taken.
✓ Branch 37 → 38 taken 2166 times.
✗ Branch 37 → 39 not taken.
✓ Branch 39 → 40 taken 2166 times.
✗ Branch 39 → 41 not taken.
✓ Branch 42 → 43 taken 6498 times.
✓ Branch 42 → 53 taken 2166 times.
✓ Branch 43 → 44 taken 6498 times.
✗ Branch 43 → 45 not taken.
✓ Branch 46 → 47 taken 6498 times.
✗ Branch 46 → 48 not taken.
✓ Branch 49 → 50 taken 6498 times.
✗ Branch 49 → 51 not taken.
✓ Branch 53 → 54 taken 6498 times.
✓ Branch 53 → 61 taken 2166 times.
✗ Branch 54 → 55 not taken.
✓ Branch 54 → 56 taken 6498 times.
✗ Branch 56 → 57 not taken.
✓ Branch 56 → 58 taken 6498 times.
✗ Branch 58 → 59 not taken.
✓ Branch 58 → 60 taken 6498 times.
✓ Branch 62 → 63 taken 6498 times.
✓ Branch 62 → 70 taken 2166 times.
✓ Branch 63 → 64 taken 6498 times.
✗ Branch 63 → 65 not taken.
✓ Branch 65 → 66 taken 6498 times.
✗ Branch 65 → 67 not taken.
✓ Branch 67 → 68 taken 6498 times.
✗ Branch 67 → 69 not taken.
|
30324 | BSplineQuadrature(tp_space%spaces(3), n_quad3)) |
| 92 | |||
| 93 | 2166 | degree_x = tp_space%spaces(1)%degree | |
| 94 | 2166 | degree_y = tp_space%spaces(2)%degree | |
| 95 | 2166 | degree_z = tp_space%spaces(3)%degree | |
| 96 | |||
| 97 |
13/20✓ Branch 70 → 71 taken 80 times.
✓ Branch 70 → 72 taken 2086 times.
✓ Branch 72 → 71 taken 2086 times.
✗ Branch 72 → 73 not taken.
✓ Branch 73 → 74 taken 20 times.
✓ Branch 73 → 75 taken 2146 times.
✓ Branch 75 → 74 taken 2146 times.
✗ Branch 75 → 76 not taken.
✗ Branch 76 → 77 not taken.
✓ Branch 76 → 78 taken 2166 times.
✓ Branch 78 → 77 taken 2166 times.
✗ Branch 78 → 79 not taken.
✓ Branch 79 → 80 taken 2166 times.
✗ Branch 79 → 81 not taken.
✓ Branch 81 → 82 taken 2066 times.
✓ Branch 81 → 83 taken 100 times.
✗ Branch 83 → 84 not taken.
✓ Branch 83 → 85 taken 2166 times.
✗ Branch 85 → 86 not taken.
✓ Branch 85 → 87 taken 2166 times.
|
10830 | allocate (userfun_vals(1:n_quad3, 1:n_quad2, 1:n_quad1)) |
| 98 | |||
| 99 |
6/12✗ Branch 87 → 88 not taken.
✓ Branch 87 → 89 taken 2166 times.
✓ Branch 89 → 88 taken 2166 times.
✗ Branch 89 → 90 not taken.
✓ Branch 90 → 91 taken 2166 times.
✗ Branch 90 → 92 not taken.
✓ Branch 92 → 93 taken 2166 times.
✗ Branch 92 → 94 not taken.
✗ Branch 94 → 95 not taken.
✓ Branch 94 → 96 taken 2166 times.
✗ Branch 96 → 97 not taken.
✓ Branch 96 → 98 taken 2166 times.
|
6498 | allocate (bspline_vals_x(n_quad1)) |
| 100 |
8/12✓ Branch 98 → 99 taken 20 times.
✓ Branch 98 → 100 taken 2146 times.
✓ Branch 100 → 99 taken 2146 times.
✗ Branch 100 → 101 not taken.
✓ Branch 101 → 102 taken 2166 times.
✗ Branch 101 → 103 not taken.
✓ Branch 103 → 104 taken 2146 times.
✓ Branch 103 → 105 taken 20 times.
✗ Branch 105 → 106 not taken.
✓ Branch 105 → 107 taken 2166 times.
✗ Branch 107 → 108 not taken.
✓ Branch 107 → 109 taken 2166 times.
|
6498 | allocate (bspline_vals_y(n_quad2)) |
| 101 |
8/12✓ Branch 109 → 110 taken 80 times.
✓ Branch 109 → 111 taken 2086 times.
✓ Branch 111 → 110 taken 2086 times.
✗ Branch 111 → 112 not taken.
✓ Branch 112 → 113 taken 2166 times.
✗ Branch 112 → 114 not taken.
✓ Branch 114 → 115 taken 2086 times.
✓ Branch 114 → 116 taken 80 times.
✗ Branch 116 → 117 not taken.
✓ Branch 116 → 118 taken 2166 times.
✗ Branch 118 → 119 not taken.
✓ Branch 118 → 120 taken 2166 times.
|
6498 | allocate (bspline_vals_z(n_quad3)) |
| 102 | |||
| 103 |
10/16✓ Branch 120 → 121 taken 20 times.
✓ Branch 120 → 122 taken 2146 times.
✓ Branch 122 → 121 taken 2146 times.
✗ Branch 122 → 123 not taken.
✗ Branch 123 → 124 not taken.
✓ Branch 123 → 125 taken 2166 times.
✓ Branch 125 → 124 taken 2166 times.
✗ Branch 125 → 126 not taken.
✓ Branch 126 → 127 taken 2166 times.
✗ Branch 126 → 128 not taken.
✓ Branch 128 → 129 taken 2146 times.
✓ Branch 128 → 130 taken 20 times.
✗ Branch 130 → 131 not taken.
✓ Branch 130 → 132 taken 2166 times.
✗ Branch 132 → 133 not taken.
✓ Branch 132 → 134 taken 2166 times.
|
8664 | allocate (partial_sum_z(n_quad2, n_quad1)) |
| 104 |
4/8✗ Branch 134 → 135 not taken.
✓ Branch 134 → 136 taken 2166 times.
✓ Branch 136 → 135 taken 2166 times.
✗ Branch 136 → 137 not taken.
✗ Branch 137 → 138 not taken.
✓ Branch 137 → 139 taken 2166 times.
✗ Branch 139 → 140 not taken.
✓ Branch 139 → 141 taken 2166 times.
|
4332 | allocate (partial_sum_yz(n_quad1)) |
| 105 | |||
| 106 |
2/2✓ Branch 141 → 142 taken 2121 times.
✓ Branch 141 → 144 taken 45 times.
|
2166 | if (v%is_initialized()) then |
| 107 | 2121 | call v%reset() | |
| 108 | else | ||
| 109 | 45 | call v%init(tp_space) ! TODO if used inside DeRham, then it is alread initialized | |
| 110 | end if | ||
| 111 | |||
| 112 |
2/2✓ Branch 147 → 148 taken 8186 times.
✓ Branch 147 → 257 taken 2166 times.
|
10352 | do kk_index = 0, nr_intervals_z - 1 |
| 113 |
2/2✓ Branch 148 → 149 taken 95548 times.
✓ Branch 148 → 255 taken 8186 times.
|
105900 | do jj_index = 0, nr_intervals_y - 1 |
| 114 |
2/2✓ Branch 149 → 150 taken 1318894 times.
✓ Branch 149 → 253 taken 95548 times.
|
1422628 | do ii_index = 0, nr_intervals_x - 1 |
| 115 | 1318894 | ii = actv_interval_inds_x(ii_index) | |
| 116 | 1318894 | jj = actv_interval_inds_y(jj_index) | |
| 117 | 1318894 | kk = actv_interval_inds_z(kk_index) | |
| 118 | 1318894 | call evaluate_at_nodes(tp_space_quad, userfun, ii, jj, kk, userfun_vals) | |
| 119 | |||
| 120 |
2/2✓ Branch 152 → 153 taken 3593180 times.
✓ Branch 152 → 251 taken 1318894 times.
|
5007622 | do k_minus_kk = 0, degree_z |
| 121 | 3593180 | call get_internal_and_actual_inds(tp_space%spaces(3), kk, k_minus_kk, k, k_internal, bspline_mirrored_z) | |
| 122 |
2/4✓ Branch 154 → 155 taken 3593180 times.
✗ Branch 154 → 249 not taken.
✗ Branch 155 → 156 not taken.
✓ Branch 155 → 157 taken 3593180 times.
|
3593180 | if (k < tp_space%rank_resp_bspline%k0 .or. k > tp_space%rank_resp_bspline%k1) cycle |
| 123 | |||
| 124 |
2/2✓ Branch 157 → 158 taken 49260 times.
✓ Branch 157 → 166 taken 3543920 times.
|
3593180 | if (bspline_mirrored_z) then |
| 125 |
4/8✓ Branch 158 → 159 taken 49260 times.
✗ Branch 158 → 160 not taken.
✗ Branch 159 → 160 not taken.
✓ Branch 159 → 163 taken 49260 times.
✗ Branch 160 → 161 not taken.
✗ Branch 160 → 162 not taken.
✓ Branch 164 → 165 taken 269316 times.
✓ Branch 164 → 174 taken 49260 times.
|
367836 | bspline_vals_z = tp_space_quad%bspline_quad(3)%bspline_at_nodes(k_internal, degree_z - k_minus_kk)%data(n_quad3:1:-1) |
| 126 | else | ||
| 127 |
4/8✓ Branch 166 → 167 taken 3543920 times.
✗ Branch 166 → 168 not taken.
✗ Branch 167 → 168 not taken.
✓ Branch 167 → 171 taken 3543920 times.
✗ Branch 168 → 169 not taken.
✗ Branch 168 → 170 not taken.
✓ Branch 172 → 173 taken 17792728 times.
✓ Branch 172 → 174 taken 3543920 times.
|
24880568 | bspline_vals_z = tp_space_quad%bspline_quad(3)%bspline_at_nodes(k_internal, k_minus_kk)%data(1:n_quad3) |
| 128 | end if | ||
| 129 |
3/4✗ Branch 174 → 175 not taken.
✓ Branch 174 → 176 taken 3593180 times.
✓ Branch 177 → 178 taken 18062044 times.
✓ Branch 177 → 179 taken 3593180 times.
|
25248404 | bspline_vals_z = bspline_vals_z * tp_space_quad%bspline_quad(3)%weights |
| 130 | |||
| 131 |
2/2✓ Branch 180 → 181 taken 19022230 times.
✓ Branch 180 → 190 taken 3593180 times.
|
22615410 | do qdx = 1, n_quad1 |
| 132 |
2/2✓ Branch 182 → 183 taken 104835885 times.
✓ Branch 182 → 188 taken 19022230 times.
|
127451295 | do qdy = 1, n_quad2 |
| 133 | val = 0._wp | ||
| 134 |
2/2✓ Branch 184 → 185 taken 550425399 times.
✓ Branch 184 → 186 taken 104835885 times.
|
655261284 | do qdz = 1, n_quad3 |
| 135 | 655261284 | val = val + bspline_vals_z(qdz) * userfun_vals(qdz, qdy, qdx) | |
| 136 | end do | ||
| 137 | 123858115 | partial_sum_z(qdy, qdx) = val | |
| 138 | end do | ||
| 139 | end do | ||
| 140 | |||
| 141 |
2/2✓ Branch 192 → 193 taken 3593180 times.
✓ Branch 192 → 194 taken 13712508 times.
|
22217762 | do j_minus_jj = 0, degree_y |
| 142 | 13712508 | call get_internal_and_actual_inds(tp_space%spaces(2), jj, j_minus_jj, j, j_internal, bspline_mirrored_y) | |
| 143 |
2/4✓ Branch 195 → 196 taken 13712508 times.
✗ Branch 195 → 247 not taken.
✗ Branch 196 → 197 not taken.
✓ Branch 196 → 198 taken 13712508 times.
|
13712508 | if (j < tp_space%rank_resp_bspline%j0 .or. j > tp_space%rank_resp_bspline%j1) cycle |
| 144 | |||
| 145 |
2/2✓ Branch 198 → 199 taken 1377290 times.
✓ Branch 198 → 207 taken 12335218 times.
|
13712508 | if (bspline_mirrored_y) then |
| 146 |
4/8✓ Branch 199 → 200 taken 1377290 times.
✗ Branch 199 → 201 not taken.
✗ Branch 200 → 201 not taken.
✓ Branch 200 → 204 taken 1377290 times.
✗ Branch 201 → 202 not taken.
✗ Branch 201 → 203 not taken.
✓ Branch 205 → 206 taken 6108838 times.
✓ Branch 205 → 215 taken 1377290 times.
|
8863418 | bspline_vals_y = tp_space_quad%bspline_quad(2)%bspline_at_nodes(j_internal, degree_y - j_minus_jj)%data(n_quad2:1:-1) |
| 147 | else | ||
| 148 |
4/8✓ Branch 207 → 208 taken 12335218 times.
✗ Branch 207 → 209 not taken.
✗ Branch 208 → 209 not taken.
✓ Branch 208 → 212 taken 12335218 times.
✗ Branch 209 → 210 not taken.
✗ Branch 209 → 211 not taken.
✓ Branch 213 → 214 taken 69639602 times.
✓ Branch 213 → 215 taken 12335218 times.
|
94310038 | bspline_vals_y = tp_space_quad%bspline_quad(2)%bspline_at_nodes(j_internal, j_minus_jj)%data(1:n_quad2) |
| 149 | end if | ||
| 150 |
3/4✗ Branch 215 → 216 not taken.
✓ Branch 215 → 217 taken 13712508 times.
✓ Branch 218 → 219 taken 75748440 times.
✓ Branch 218 → 220 taken 13712508 times.
|
103173456 | bspline_vals_y = bspline_vals_y * tp_space_quad%bspline_quad(2)%weights |
| 151 | |||
| 152 |
2/2✓ Branch 221 → 222 taken 74792097 times.
✓ Branch 221 → 227 taken 13712508 times.
|
88504605 | do qdx = 1, n_quad1 |
| 153 | val = 0._wp | ||
| 154 |
2/2✓ Branch 223 → 224 taken 432878709 times.
✓ Branch 223 → 225 taken 74792097 times.
|
507670806 | do qdy = 1, n_quad2 |
| 155 | 507670806 | val = val + bspline_vals_y(qdy) * partial_sum_z(qdy, qdx) | |
| 156 | end do | ||
| 157 | 88504605 | partial_sum_yz(qdx) = val | |
| 158 | end do | ||
| 159 | |||
| 160 |
2/2✓ Branch 229 → 230 taken 13712508 times.
✓ Branch 229 → 231 taken 54944533 times.
|
85962729 | do i_minus_ii = 0, degree_x |
| 161 | 54944533 | call get_internal_and_actual_inds(tp_space%spaces(1), ii, i_minus_ii, i, i_internal, bspline_mirrored_x) | |
| 162 |
2/4✗ Branch 232 → 233 not taken.
✓ Branch 232 → 234 taken 54944533 times.
✗ Branch 234 → 235 not taken.
✓ Branch 234 → 236 taken 54944533 times.
|
54944533 | if (i < tp_space%rank_resp_bspline%i0 .or. i > tp_space%rank_resp_bspline%i1) cycle |
| 163 | |||
| 164 | val = 0._wp | ||
| 165 |
2/2✓ Branch 236 → 237 taken 5867222 times.
✓ Branch 236 → 240 taken 49077311 times.
|
54944533 | if (bspline_mirrored_x) then |
| 166 |
2/2✓ Branch 238 → 239 taken 36267239 times.
✓ Branch 238 → 244 taken 5867222 times.
|
42134461 | do qdx = 1, n_quad1 |
| 167 | val = val + partial_sum_yz(qdx) * & | ||
| 168 | & tp_space_quad%bspline_quad(1)%bspline_at_nodes(i_internal, degree_x - i_minus_ii)%data(n_quad1 + 1 - qdx) * & | ||
| 169 | 42134461 | & tp_space_quad%bspline_quad(1)%weights(qdx) | |
| 170 | end do | ||
| 171 | else | ||
| 172 |
2/2✓ Branch 241 → 242 taken 49077311 times.
✓ Branch 241 → 243 taken 275851050 times.
|
324928361 | do qdx = 1, n_quad1 |
| 173 | val = val + partial_sum_yz(qdx) * & | ||
| 174 | & tp_space_quad%bspline_quad(1)%bspline_at_nodes(i_internal, i_minus_ii)%data(qdx) * & | ||
| 175 | 324928361 | & tp_space_quad%bspline_quad(1)%weights(qdx) | |
| 176 | end do | ||
| 177 | end if | ||
| 178 | |||
| 179 | 123601574 | v%data(i, j, k) = v%data(i, j, k) + val | |
| 180 | |||
| 181 | end do | ||
| 182 | end do | ||
| 183 | end do | ||
| 184 | |||
| 185 | end do | ||
| 186 | end do | ||
| 187 | end do | ||
| 188 | |||
| 189 |
1/2✗ Branch 258 → 259 not taken.
✓ Branch 258 → 260 taken 2166 times.
|
2166 | deallocate (userfun_vals) |
| 190 | |||
| 191 |
1/2✗ Branch 260 → 261 not taken.
✓ Branch 260 → 262 taken 2166 times.
|
2166 | deallocate (partial_sum_z) |
| 192 |
1/2✗ Branch 262 → 263 not taken.
✓ Branch 262 → 264 taken 2166 times.
|
2166 | deallocate (partial_sum_yz) |
| 193 | |||
| 194 |
1/2✗ Branch 264 → 265 not taken.
✓ Branch 264 → 266 taken 2166 times.
|
2166 | deallocate (bspline_vals_x) |
| 195 |
1/2✗ Branch 266 → 267 not taken.
✓ Branch 266 → 268 taken 2166 times.
|
2166 | deallocate (bspline_vals_y) |
| 196 |
1/2✗ Branch 268 → 269 not taken.
✓ Branch 268 → 270 taken 2166 times.
|
2166 | deallocate (bspline_vals_z) |
| 197 | |||
| 198 |
1/2✗ Branch 270 → 271 not taken.
✓ Branch 270 → 272 taken 2166 times.
|
2166 | deallocate (actv_interval_inds_x) |
| 199 |
1/2✗ Branch 272 → 273 not taken.
✓ Branch 272 → 274 taken 2166 times.
|
2166 | deallocate (actv_interval_inds_y) |
| 200 |
1/2✗ Branch 274 → 275 not taken.
✓ Branch 274 → 276 taken 2166 times.
|
2166 | deallocate (actv_interval_inds_z) |
| 201 | |||
| 202 | 2166 | call tp_space_quad%destroy() | |
| 203 | |||
| 204 |
5/8✓ Branch 278 → 279 taken 6498 times.
✓ Branch 278 → 286 taken 2166 times.
✗ Branch 279 → 280 not taken.
✓ Branch 279 → 281 taken 6498 times.
✗ Branch 281 → 282 not taken.
✓ Branch 281 → 283 taken 6498 times.
✗ Branch 283 → 284 not taken.
✓ Branch 283 → 285 taken 6498 times.
|
8664 | end subroutine inner_product_3d |
| 205 | |||
| 206 | !> @brief Compute the L2 projection of a user-defined function onto a tensor product B-spline space | ||
| 207 | !> | ||
| 208 | !> @param[out] tp_fun The resulting tensor product function containing the L2 projection | ||
| 209 | !> @param[in] tp_space The tensor product space onto which the L2 projection is computed | ||
| 210 | !> @param[in] userfun The user-defined function to project onto the tensor product space | ||
| 211 | !> @param[in] coordfun (optional) A coordinate function used in the definition of the mass matrix | ||
| 212 | !> @param[in] rtol (optional) The relative tolerance for the linear solver | ||
| 213 | !> @param[out] l2_err (optional) An estimate of the L2 error of the projection | ||
| 214 | !> | ||
| 215 | !> @note It is not recommended to use this routine because each time it is called a new mass matrix and linear solver are created | ||
| 216 | 45 | subroutine l2_projection_3d(tp_fun, tp_space, userfun, coordfun, rtol, l2_err) | |
| 217 | use m_tensorprod_matrix, only: TensorProdMat, mass_matrix, get_petsc | ||
| 218 | use m_tensorprod_basis, only: get_petsc | ||
| 219 | use m_tensorprod_quadrature, only: TensorProdQuadrature, evaluate_at_nodes | ||
| 220 | implicit none | ||
| 221 | |||
| 222 | type(TensorProdFun), intent(out) :: tp_fun | ||
| 223 | type(TensorProdSpace), intent(in) :: tp_space | ||
| 224 | procedure(user_function_3d_interface) :: userfun | ||
| 225 | procedure(user_function_3d_interface), optional :: coordfun | ||
| 226 | real(wp), intent(in), optional :: rtol | ||
| 227 | real(wp), intent(out), optional :: l2_err | ||
| 228 | |||
| 229 | Vec :: c | ||
| 230 | Vec :: f | ||
| 231 | 45 | type(TensorProdMat) :: M_TP | |
| 232 | type(TensorProdFun) :: v_tp | ||
| 233 | Mat :: M | ||
| 234 | KSP :: ksp | ||
| 235 | integer :: ierr | ||
| 236 | real(wp) :: rtol_ | ||
| 237 | |||
| 238 | 45 | if (present(rtol)) then | |
| 239 | 30 | rtol_ = rtol | |
| 240 | else | ||
| 241 | 15 | rtol_ = 1.e-10_wp | |
| 242 | end if | ||
| 243 | |||
| 244 | 45 | call mass_matrix(M_TP, tp_space, tp_space, coordfun=coordfun) | |
| 245 | 45 | call inner_product(v_tp, tp_space, userfun) | |
| 246 | |||
| 247 | 45 | call get_petsc(M, M_TP) | |
| 248 | 45 | call get_petsc(f, v_tp) | |
| 249 | |||
| 250 | ! Create linear solver | ||
| 251 |
1/2✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 45 times.
|
45 | PetscCall(KSPCreate(tp_space%domain%comm, ksp, ierr)) |
| 252 |
1/2✗ Branch 13 → 14 not taken.
✓ Branch 13 → 15 taken 45 times.
|
45 | PetscCall(KSPSetOperators(ksp, M, M, ierr)) |
| 253 |
1/2✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 45 times.
|
45 | PetscCall(KSPSetFromOptions(ksp, ierr)) |
| 254 |
1/2✗ Branch 19 → 20 not taken.
✓ Branch 19 → 21 taken 45 times.
|
45 | PetscCall(KSPSetType(ksp, "cg", ierr)) |
| 255 | |||
| 256 |
1/2✗ Branch 22 → 23 not taken.
✓ Branch 22 → 24 taken 45 times.
|
45 | PetscCall(KSPSetTolerances(ksp, rtol_, 1.e-30_wp, 1.e10_wp, 1000, ierr)) |
| 257 | |||
| 258 |
1/2✗ Branch 25 → 26 not taken.
✓ Branch 25 → 27 taken 45 times.
|
45 | PetscCall(VecDuplicate(f, c, ierr)) |
| 259 |
1/2✗ Branch 28 → 29 not taken.
✓ Branch 28 → 30 taken 45 times.
|
45 | PetscCall(KSPSolve(ksp, f, c, ierr)) |
| 260 | |||
| 261 |
1/2✓ Branch 30 → 31 taken 45 times.
✗ Branch 30 → 36 not taken.
|
45 | if (present(l2_err)) then |
| 262 | ! L2-error = c' * M * c - 2 * c' * f + int F^2 dx = int f^2 dx - c' * F | ||
| 263 |
1/2✗ Branch 32 → 33 not taken.
✓ Branch 32 → 34 taken 45 times.
|
45 | PetscCall(VecDot(c, f, l2_err, ierr)) |
| 264 | |||
| 265 | 45 | l2_err = l2_err - integrate(tp_space, userfun_squared) | |
| 266 | |||
| 267 | 45 | l2_err = sqrt(abs(l2_err)) | |
| 268 | end if | ||
| 269 | |||
| 270 | 45 | call tp_fun%init(tp_space, vec=c) | |
| 271 | |||
| 272 |
1/2✗ Branch 38 → 39 not taken.
✓ Branch 38 → 40 taken 45 times.
|
45 | PetscCall(MatDestroy(M, ierr)) |
| 273 |
1/2✗ Branch 41 → 42 not taken.
✓ Branch 41 → 43 taken 45 times.
|
45 | PetscCall(VecDestroy(f, ierr)) |
| 274 |
1/2✗ Branch 44 → 45 not taken.
✓ Branch 44 → 46 taken 45 times.
|
45 | PetscCall(VecDestroy(c, ierr)) |
| 275 |
1/2✗ Branch 47 → 48 not taken.
✓ Branch 47 → 49 taken 45 times.
|
45 | PetscCall(KSPDestroy(ksp, ierr)) |
| 276 | |||
| 277 | 45 | call v_tp%destroy() | |
| 278 |
5/14✓ Branch 2 → 3 taken 30 times.
✓ Branch 2 → 4 taken 15 times.
✗ Branch 51 → 52 not taken.
✓ Branch 51 → 53 taken 45 times.
✗ Branch 53 → 54 not taken.
✓ Branch 53 → 55 taken 45 times.
✓ Branch 55 → 56 taken 45 times.
✗ Branch 55 → 68 not taken.
✗ Branch 57 → 58 not taken.
✗ Branch 57 → 59 not taken.
✗ Branch 61 → 62 not taken.
✗ Branch 61 → 63 not taken.
✗ Branch 65 → 66 not taken.
✗ Branch 65 → 69 not taken.
|
90 | call M_TP%destroy() |
| 279 | contains | ||
| 280 | 3367728 | pure real(wp) function userfun_squared(x, y, z) result(val) | |
| 281 | real(wp), intent(in) :: x, y, z | ||
| 282 | 3367728 | val = userfun(x, y, z)**2 | |
| 283 | 3367728 | end function userfun_squared | |
| 284 | end subroutine | ||
| 285 | |||
| 286 | !> @brief Integrate a user-defined function over the subintervals defined by a tensor product space using quadrature | ||
| 287 | !> | ||
| 288 | !> @param[in] tp_space The tensor product space over which to integrate | ||
| 289 | !> @param[in] userfun The user-defined function to integrate | ||
| 290 | !> @param[in] n_quad_extra (optional) The number of quadrature points to use in each direction, on top of the number needed for | ||
| 291 | !> exact integration of a product of B-splines | ||
| 292 | !> | ||
| 293 | !> @return The result of the integration | ||
| 294 | 770 | real(wp) function integrate_3d(tp_space, userfun, n_quad_extra) result(ans) | |
| 295 | use m_bspline_quadrature, only: BSplineQuadrature, n_quad_exact, DEFAULT_EXTRA_QUADRATURE_POINTS | ||
| 296 | use m_tensorprod_basis, only: TensorProdSpace | ||
| 297 | use m_tensorprod_quadrature, only: TensorProdQuadrature, evaluate_at_nodes | ||
| 298 | implicit none | ||
| 299 | |||
| 300 | type(TensorProdSpace), intent(in) :: tp_space | ||
| 301 | procedure(user_function_3d_interface) :: userfun | ||
| 302 | integer, intent(in), optional :: n_quad_extra | ||
| 303 | |||
| 304 | integer :: n_quad1, n_quad2, n_quad3, n_quad_extra_ | ||
| 305 | integer :: ii, jj, kk, qdx, qdy, qdz, ierr | ||
| 306 | real(wp) :: weight_x, weight_xy, local_ans | ||
| 307 | 770 | real(wp), allocatable :: userfun_vals(:, :, :) | |
| 308 |
2/2✓ Branch 3 → 4 taken 2310 times.
✓ Branch 3 → 5 taken 770 times.
|
3080 | type(TensorProdQuadrature) :: tp_space_quad |
| 309 | |||
| 310 | 770 | n_quad1 = n_quad_exact(tp_space%spaces(1), tp_space%spaces(1)) | |
| 311 | 770 | n_quad2 = n_quad_exact(tp_space%spaces(2), tp_space%spaces(2)) | |
| 312 | 770 | n_quad3 = n_quad_exact(tp_space%spaces(3), tp_space%spaces(3)) | |
| 313 |
2/2✓ Branch 5 → 6 taken 86 times.
✓ Branch 5 → 7 taken 684 times.
|
770 | if (present(n_quad_extra)) then |
| 314 | 86 | n_quad_extra_ = n_quad_extra | |
| 315 | else | ||
| 316 | n_quad_extra_ = DEFAULT_EXTRA_QUADRATURE_POINTS | ||
| 317 | end if | ||
| 318 | |||
| 319 | ! The number of quadrature points in each direction is increased | ||
| 320 | 770 | n_quad1 = n_quad1 + n_quad_extra_ | |
| 321 | 770 | n_quad2 = n_quad2 + n_quad_extra_ | |
| 322 | 770 | n_quad3 = n_quad3 + n_quad_extra_ | |
| 323 | |||
| 324 | tp_space_quad = TensorProdQuadrature( & | ||
| 325 | BSplineQuadrature(tp_space%spaces(1), n_quad1), & | ||
| 326 | BSplineQuadrature(tp_space%spaces(2), n_quad2), & | ||
| 327 |
29/50✓ Branch 14 → 15 taken 2310 times.
✓ Branch 14 → 22 taken 770 times.
✗ Branch 15 → 16 not taken.
✓ Branch 15 → 17 taken 2310 times.
✗ Branch 17 → 18 not taken.
✓ Branch 17 → 19 taken 2310 times.
✗ Branch 19 → 20 not taken.
✓ Branch 19 → 21 taken 2310 times.
✓ Branch 22 → 23 taken 770 times.
✗ Branch 22 → 24 not taken.
✓ Branch 24 → 25 taken 770 times.
✗ Branch 24 → 26 not taken.
✓ Branch 26 → 27 taken 770 times.
✗ Branch 26 → 28 not taken.
✓ Branch 28 → 29 taken 770 times.
✗ Branch 28 → 30 not taken.
✓ Branch 30 → 31 taken 770 times.
✗ Branch 30 → 32 not taken.
✓ Branch 32 → 33 taken 770 times.
✗ Branch 32 → 34 not taken.
✓ Branch 34 → 35 taken 770 times.
✗ Branch 34 → 36 not taken.
✓ Branch 36 → 37 taken 770 times.
✗ Branch 36 → 38 not taken.
✓ Branch 38 → 39 taken 770 times.
✗ Branch 38 → 40 not taken.
✓ Branch 41 → 42 taken 2310 times.
✓ Branch 41 → 52 taken 770 times.
✓ Branch 42 → 43 taken 2310 times.
✗ Branch 42 → 44 not taken.
✓ Branch 45 → 46 taken 2310 times.
✗ Branch 45 → 47 not taken.
✓ Branch 48 → 49 taken 2310 times.
✗ Branch 48 → 50 not taken.
✓ Branch 52 → 53 taken 2310 times.
✓ Branch 52 → 60 taken 770 times.
✗ Branch 53 → 54 not taken.
✓ Branch 53 → 55 taken 2310 times.
✗ Branch 55 → 56 not taken.
✓ Branch 55 → 57 taken 2310 times.
✗ Branch 57 → 58 not taken.
✓ Branch 57 → 59 taken 2310 times.
✓ Branch 61 → 62 taken 2310 times.
✓ Branch 61 → 69 taken 770 times.
✓ Branch 62 → 63 taken 2310 times.
✗ Branch 62 → 64 not taken.
✓ Branch 64 → 65 taken 2310 times.
✗ Branch 64 → 66 not taken.
✓ Branch 66 → 67 taken 2310 times.
✗ Branch 66 → 68 not taken.
|
10780 | BSplineQuadrature(tp_space%spaces(3), n_quad3)) |
| 328 | |||
| 329 | allocate (userfun_vals(1:tp_space_quad%bspline_quad(3)%n_quad, & | ||
| 330 | 1:tp_space_quad%bspline_quad(2)%n_quad, & | ||
| 331 |
10/20✗ Branch 69 → 70 not taken.
✓ Branch 69 → 71 taken 770 times.
✓ Branch 71 → 70 taken 770 times.
✗ Branch 71 → 72 not taken.
✗ Branch 72 → 73 not taken.
✓ Branch 72 → 74 taken 770 times.
✓ Branch 74 → 73 taken 770 times.
✗ Branch 74 → 75 not taken.
✗ Branch 75 → 76 not taken.
✓ Branch 75 → 77 taken 770 times.
✓ Branch 77 → 76 taken 770 times.
✗ Branch 77 → 78 not taken.
✓ Branch 78 → 79 taken 770 times.
✗ Branch 78 → 80 not taken.
✓ Branch 80 → 81 taken 770 times.
✗ Branch 80 → 82 not taken.
✗ Branch 82 → 83 not taken.
✓ Branch 82 → 84 taken 770 times.
✗ Branch 84 → 85 not taken.
✓ Branch 84 → 86 taken 770 times.
|
3850 | 1:tp_space_quad%bspline_quad(1)%n_quad)) |
| 332 | |||
| 333 | 770 | local_ans = 0._wp | |
| 334 |
2/2✓ Branch 87 → 88 taken 4240 times.
✓ Branch 87 → 109 taken 770 times.
|
5010 | do kk = tp_space%rank_resp_intervals%k0, tp_space%rank_resp_intervals%k1 |
| 335 |
2/2✓ Branch 89 → 90 taken 39017 times.
✓ Branch 89 → 107 taken 4240 times.
|
44027 | do jj = tp_space%rank_resp_intervals%j0, tp_space%rank_resp_intervals%j1 |
| 336 |
2/2✓ Branch 91 → 92 taken 515995 times.
✓ Branch 91 → 105 taken 39017 times.
|
559252 | do ii = tp_space%rank_resp_intervals%i0, tp_space%rank_resp_intervals%i1 |
| 337 | 515995 | call evaluate_at_nodes(tp_space_quad, userfun, ii, jj, kk, userfun_vals) | |
| 338 |
2/2✓ Branch 93 → 94 taken 3151862 times.
✓ Branch 93 → 103 taken 515995 times.
|
3706874 | do qdx = 1, tp_space_quad%bspline_quad(1)%n_quad |
| 339 | 3151862 | weight_x = tp_space_quad%bspline_quad(1)%weights(qdx) | |
| 340 |
2/2✓ Branch 95 → 96 taken 19677494 times.
✓ Branch 95 → 101 taken 3151862 times.
|
23345351 | do qdy = 1, tp_space_quad%bspline_quad(2)%n_quad |
| 341 | 19677494 | weight_xy = tp_space_quad%bspline_quad(2)%weights(qdy) * weight_x | |
| 342 |
2/2✓ Branch 97 → 98 taken 106997674 times.
✓ Branch 97 → 99 taken 19677494 times.
|
129827030 | do qdz = 1, tp_space_quad%bspline_quad(3)%n_quad |
| 343 | 126675168 | local_ans = local_ans + weight_xy * tp_space_quad%bspline_quad(3)%weights(qdz) * userfun_vals(qdz, qdy, qdx) | |
| 344 | end do | ||
| 345 | end do | ||
| 346 | end do | ||
| 347 | end do | ||
| 348 | end do | ||
| 349 | end do | ||
| 350 | |||
| 351 |
1/2✗ Branch 111 → 112 not taken.
✓ Branch 111 → 114 taken 770 times.
|
770 | PetscCallMPI(MPI_Allreduce(local_ans, ans, 1, MPI_WP, MPI_SUM, tp_space%domain%comm, ierr)) |
| 352 | |||
| 353 |
1/2✗ Branch 114 → 115 not taken.
✓ Branch 114 → 116 taken 770 times.
|
770 | deallocate (userfun_vals) |
| 354 | |||
| 355 |
6/10✗ Branch 117 → 118 not taken.
✓ Branch 117 → 119 taken 770 times.
✓ Branch 120 → 121 taken 2310 times.
✓ Branch 120 → 128 taken 770 times.
✓ Branch 121 → 122 taken 2310 times.
✗ Branch 121 → 123 not taken.
✓ Branch 123 → 124 taken 2310 times.
✗ Branch 123 → 125 not taken.
✓ Branch 125 → 126 taken 2310 times.
✗ Branch 125 → 127 not taken.
|
3850 | end function integrate_3d |
| 356 | |||
| 357 | !> @brief Compute the L2 norm of the difference between a tensor product B-spline function and a user-defined function | ||
| 358 | !> | ||
| 359 | !> @param[in] tp_fun The tensor product function representing the B-spline function | ||
| 360 | !> @param[in] userfun The user-defined function to compare against | ||
| 361 | !> @param[in] n_quad_extra (optional) The number of quadrature points to use in each direction, on top of the number needed for | ||
| 362 | !> exact integration of a product of B-splines | ||
| 363 | 45 | real(wp) function l2_error_3d(tp_fun, userfun, n_quad_extra) result(ans) | |
| 364 | use m_tensorprod_basis, only: evaluate | ||
| 365 | |||
| 366 | implicit none | ||
| 367 | |||
| 368 | type(TensorProdFun), intent(in) :: tp_fun | ||
| 369 | procedure(user_function_3d_interface) :: userfun | ||
| 370 | integer, intent(in), optional :: n_quad_extra | ||
| 371 | |||
| 372 | 45 | ans = sqrt(integrate(tp_fun%space, int_function, n_quad_extra=n_quad_extra)) | |
| 373 | contains | ||
| 374 | 3367728 | pure real(wp) function int_function(x, y, z) result(val) | |
| 375 | real(wp), intent(in) :: x, y, z | ||
| 376 | 3367728 | val = (userfun(x, y, z) - evaluate(tp_fun, x, y, z))**2 | |
| 377 | 3367728 | end function int_function | |
| 378 | end function l2_error_3d | ||
| 379 | |||
| 380 |
4/4✓ Branch 8 → 9 taken 2310 times.
✓ Branch 8 → 10 taken 770 times.
✓ Branch 9 → 10 taken 6498 times.
✓ Branch 9 → 11 taken 2166 times.
|
11744 | end module m_tensorprod_linalg |
| 381 |