quadrature/m_quadrature.f90
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | !> @brief Module for numerical quadrature. | ||
| 2 | !> | ||
| 3 | !> This module provides functions for numerical quadrature in 1D and 3D. | ||
| 4 | module m_quadrature | ||
| 5 | use m_common, only: wp, user_function_1d_interface, user_function_3d_interface | ||
| 6 | use m_quadrature_data | ||
| 7 | use m_quadrature_functions | ||
| 8 | |||
| 9 | implicit none | ||
| 10 | |||
| 11 | private | ||
| 12 | public :: quadrature, MAX_QUADRATURE_NODES | ||
| 13 | |||
| 14 | !> @brief Numerical quadrature | ||
| 15 | interface quadrature | ||
| 16 | module procedure quadrature_1d, quadrature_3d | ||
| 17 | module procedure quadrature_composite_1d, quadrature_composite_3d | ||
| 18 | end interface | ||
| 19 | |||
| 20 | contains | ||
| 21 | |||
| 22 | !> @brief Numerical quadrature in 3D | ||
| 23 | !> | ||
| 24 | !> @param[in] n_quad Number of quadrature nodes | ||
| 25 | !> @param[in] f Function to integrate | ||
| 26 | !> @param[in] x0,x1,y0,y1,z0,z1 Integration limits in the x, y and z directions | ||
| 27 | !> | ||
| 28 | !> @return Integral of f over the domain | ||
| 29 | 14584 | pure real(wp) function quadrature_3d(n_quad, f, x0, x1, y0, y1, z0, z1) result(ans) | |
| 30 | implicit none | ||
| 31 | |||
| 32 | integer, intent(in) :: n_quad | ||
| 33 | procedure(user_function_3d_interface) :: f | ||
| 34 | real(wp), intent(in) :: x0, x1, y0, y1, z0, z1 | ||
| 35 | |||
| 36 | 19268 | select case (n_quad) | |
| 37 | case (1) | ||
| 38 | 4684 | ans = quadrature_3d_arrayinput(gl_weights_1, gl_nodes_1, f, x0, x1, y0, y1, z0, z1) | |
| 39 | case (2) | ||
| 40 | 4688 | ans = quadrature_3d_arrayinput(gl_weights_2, gl_nodes_2, f, x0, x1, y0, y1, z0, z1) | |
| 41 | case (3) | ||
| 42 | 4692 | ans = quadrature_3d_arrayinput(gl_weights_3, gl_nodes_3, f, x0, x1, y0, y1, z0, z1) | |
| 43 | case (4) | ||
| 44 | 16 | ans = quadrature_3d_arrayinput(gl_weights_4, gl_nodes_4, f, x0, x1, y0, y1, z0, z1) | |
| 45 | case (5) | ||
| 46 | 20 | ans = quadrature_3d_arrayinput(gl_weights_5, gl_nodes_5, f, x0, x1, y0, y1, z0, z1) | |
| 47 | case (6) | ||
| 48 | 24 | ans = quadrature_3d_arrayinput(gl_weights_6, gl_nodes_6, f, x0, x1, y0, y1, z0, z1) | |
| 49 | case (7) | ||
| 50 | 28 | ans = quadrature_3d_arrayinput(gl_weights_7, gl_nodes_7, f, x0, x1, y0, y1, z0, z1) | |
| 51 | case (8) | ||
| 52 | 32 | ans = quadrature_3d_arrayinput(gl_weights_8, gl_nodes_8, f, x0, x1, y0, y1, z0, z1) | |
| 53 | case (9) | ||
| 54 | 36 | ans = quadrature_3d_arrayinput(gl_weights_9, gl_nodes_9, f, x0, x1, y0, y1, z0, z1) | |
| 55 | case (10) | ||
| 56 | 40 | ans = quadrature_3d_arrayinput(gl_weights_10, gl_nodes_10, f, x0, x1, y0, y1, z0, z1) | |
| 57 | case (11) | ||
| 58 | 44 | ans = quadrature_3d_arrayinput(gl_weights_11, gl_nodes_11, f, x0, x1, y0, y1, z0, z1) | |
| 59 | case (12) | ||
| 60 | 48 | ans = quadrature_3d_arrayinput(gl_weights_12, gl_nodes_12, f, x0, x1, y0, y1, z0, z1) | |
| 61 | case (13) | ||
| 62 | 52 | ans = quadrature_3d_arrayinput(gl_weights_13, gl_nodes_13, f, x0, x1, y0, y1, z0, z1) | |
| 63 | case (14) | ||
| 64 | 56 | ans = quadrature_3d_arrayinput(gl_weights_14, gl_nodes_14, f, x0, x1, y0, y1, z0, z1) | |
| 65 | case (15) | ||
| 66 | 60 | ans = quadrature_3d_arrayinput(gl_weights_15, gl_nodes_15, f, x0, x1, y0, y1, z0, z1) | |
| 67 | case (16) | ||
| 68 | 64 | ans = quadrature_3d_arrayinput(gl_weights_16, gl_nodes_16, f, x0, x1, y0, y1, z0, z1) | |
| 69 | case default | ||
| 70 |
16/17✓ Branch 2 → 3 taken 4684 times.
✓ Branch 2 → 4 taken 4688 times.
✓ Branch 2 → 5 taken 4692 times.
✓ Branch 2 → 6 taken 16 times.
✓ Branch 2 → 7 taken 20 times.
✓ Branch 2 → 8 taken 24 times.
✓ Branch 2 → 9 taken 28 times.
✓ Branch 2 → 10 taken 32 times.
✓ Branch 2 → 11 taken 36 times.
✓ Branch 2 → 12 taken 40 times.
✓ Branch 2 → 13 taken 44 times.
✓ Branch 2 → 14 taken 48 times.
✓ Branch 2 → 15 taken 52 times.
✓ Branch 2 → 16 taken 56 times.
✓ Branch 2 → 17 taken 60 times.
✓ Branch 2 → 18 taken 64 times.
✗ Branch 2 → 19 not taken.
|
14584 | error stop "quadrature parameter N incorrect" |
| 71 | end select | ||
| 72 | |||
| 73 | 14584 | end function | |
| 74 | |||
| 75 | !> @brief Numerical quadrature in 3D | ||
| 76 | !> | ||
| 77 | !> @param[in] weights Weights for the quadrature | ||
| 78 | !> @param[in] nodes Nodes for the quadrature | ||
| 79 | !> @param[in] f Function to integrate | ||
| 80 | !> @param[in] x0,x1,y0,y1,z0,z1 Integration limits in the x, y and z directions | ||
| 81 | !> | ||
| 82 | !> @return Integral of f over the domain | ||
| 83 |
2/4✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 14584 times.
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 14584 times.
|
14584 | pure real(wp) function quadrature_3d_arrayinput(weights, nodes, f, x0, x1, y0, y1, z0, z1) result(ans) |
| 84 | implicit none | ||
| 85 | |||
| 86 | real(wp), intent(in) :: weights(1:), nodes(1:) | ||
| 87 | procedure(user_function_3d_interface) :: f | ||
| 88 | real(wp), intent(in) :: x0, x1, y0, y1, z0, z1 | ||
| 89 | |||
| 90 | integer :: i, j, k | ||
| 91 | real(wp) :: hx, hy, hz, xmid, ymid, zmid | ||
| 92 | |||
| 93 | 14584 | hx = (x1 - x0) / 2 | |
| 94 | 14584 | xmid = (x0 + x1) / 2 | |
| 95 | 14584 | hy = (y1 - y0) / 2 | |
| 96 | 14584 | ymid = (y0 + y1) / 2 | |
| 97 | 14584 | hz = (z1 - z0) / 2 | |
| 98 | 14584 | zmid = (z0 + z1) / 2 | |
| 99 | |||
| 100 | ans = 0.0_wp | ||
| 101 |
2/2✓ Branch 7 → 8 taken 34064 times.
✓ Branch 7 → 18 taken 14584 times.
|
48648 | do i = 1, size(weights) |
| 102 |
2/2✓ Branch 9 → 10 taken 139504 times.
✓ Branch 9 → 16 taken 34064 times.
|
188152 | do j = 1, size(weights) |
| 103 |
2/2✓ Branch 11 → 12 taken 1143872 times.
✓ Branch 11 → 14 taken 139504 times.
|
1317440 | do k = 1, size(weights) |
| 104 | ans = ans + weights(i) * weights(j) * weights(k) * & | ||
| 105 | 1283376 | & f(xmid + hx * nodes(i), ymid + hy * nodes(j), zmid + hz * nodes(k)) | |
| 106 | end do | ||
| 107 | end do | ||
| 108 | end do | ||
| 109 | |||
| 110 | 14584 | ans = hx * hy * hz * ans | |
| 111 | |||
| 112 | 14584 | end function | |
| 113 | |||
| 114 | !> @brief Numerical quadrature in 1D | ||
| 115 | !> | ||
| 116 | !> @param[in] n_quad Number of quadrature nodes | ||
| 117 | !> @param[in] f Function to integrate | ||
| 118 | !> @param[in] xleft,xright Integration limits | ||
| 119 | !> | ||
| 120 | !> @return Integral of f over the domain | ||
| 121 | 3758770 | pure real(wp) function quadrature_1d(n_quad, f, xleft, xright) result(ans) | |
| 122 | implicit none | ||
| 123 | |||
| 124 | integer, intent(in) :: n_quad | ||
| 125 | procedure(user_function_1d_interface) :: f | ||
| 126 | real(wp), intent(in) :: xleft, xright | ||
| 127 | |||
| 128 | 3801334 | select case (n_quad) | |
| 129 | case (1) | ||
| 130 | 42564 | ans = quadrature_1(f, xleft, xright) | |
| 131 | case (2) | ||
| 132 | 70771 | ans = quadrature_2(f, xleft, xright) | |
| 133 | case (3) | ||
| 134 | 139646 | ans = quadrature_3(f, xleft, xright) | |
| 135 | case (4) | ||
| 136 | 638075 | ans = quadrature_4(f, xleft, xright) | |
| 137 | case (5) | ||
| 138 | 2505391 | ans = quadrature_5(f, xleft, xright) | |
| 139 | case (6) | ||
| 140 | 65161 | ans = quadrature_6(f, xleft, xright) | |
| 141 | case (7) | ||
| 142 | 76389 | ans = quadrature_7(f, xleft, xright) | |
| 143 | case (8) | ||
| 144 | 98079 | ans = quadrature_8(f, xleft, xright) | |
| 145 | case (9) | ||
| 146 | 121759 | ans = quadrature_9(f, xleft, xright) | |
| 147 | case (10) | ||
| 148 | 407 | ans = quadrature_10(f, xleft, xright) | |
| 149 | case (11) | ||
| 150 | 451 | ans = quadrature_11(f, xleft, xright) | |
| 151 | case (12) | ||
| 152 | 1 | ans = quadrature_12(f, xleft, xright) | |
| 153 | case (13) | ||
| 154 | 73 | ans = quadrature_13(f, xleft, xright) | |
| 155 | case (14) | ||
| 156 | 1 | ans = quadrature_14(f, xleft, xright) | |
| 157 | case (15) | ||
| 158 | 1 | ans = quadrature_15(f, xleft, xright) | |
| 159 | case (16) | ||
| 160 | 1 | ans = quadrature_16(f, xleft, xright) | |
| 161 | case default | ||
| 162 |
16/17✓ Branch 2 → 3 taken 42564 times.
✓ Branch 2 → 4 taken 70771 times.
✓ Branch 2 → 5 taken 139646 times.
✓ Branch 2 → 6 taken 638075 times.
✓ Branch 2 → 7 taken 2505391 times.
✓ Branch 2 → 8 taken 65161 times.
✓ Branch 2 → 9 taken 76389 times.
✓ Branch 2 → 10 taken 98079 times.
✓ Branch 2 → 11 taken 121759 times.
✓ Branch 2 → 12 taken 407 times.
✓ Branch 2 → 13 taken 451 times.
✓ Branch 2 → 14 taken 1 time.
✓ Branch 2 → 15 taken 73 times.
✓ Branch 2 → 16 taken 1 time.
✓ Branch 2 → 17 taken 1 time.
✓ Branch 2 → 18 taken 1 time.
✗ Branch 2 → 19 not taken.
|
3758770 | error stop "quadrature parameter N incorrect" |
| 163 | end select | ||
| 164 | 3758770 | end function | |
| 165 | |||
| 166 | !> @brief Composite numerical quadrature in 1D | ||
| 167 | !> | ||
| 168 | !> @param[in] n_quad Number of quadrature nodes per interval | ||
| 169 | !> @param[in] f Function to integrate | ||
| 170 | !> @param[in] xleft,xright Integration limits | ||
| 171 | !> @param[in] n_intervals Number of intervals for the composite quadrature | ||
| 172 | !> | ||
| 173 | !> @return Integral of f over the domain | ||
| 174 | 12 | pure real(wp) function quadrature_composite_1d(n_quad, f, xleft, xright, n_intervals) result(ans) | |
| 175 | implicit none | ||
| 176 | |||
| 177 | integer, intent(in) :: n_quad, n_intervals | ||
| 178 | procedure(user_function_1d_interface) :: f | ||
| 179 | real(wp), intent(in) :: xleft, xright | ||
| 180 | |||
| 181 | integer :: i | ||
| 182 | real(wp) :: h, x_i, x_ip1 | ||
| 183 | |||
| 184 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 12 times.
|
12 | if (n_intervals < 1) error stop "quadrature composite parameter n_intervals incorrect" |
| 185 | |||
| 186 | 12 | h = (xright - xleft) / n_intervals | |
| 187 | ans = 0.0_wp | ||
| 188 |
2/2✓ Branch 5 → 6 taken 90 times.
✓ Branch 5 → 7 taken 12 times.
|
102 | do i = 1, n_intervals |
| 189 | 90 | x_i = xleft + (i - 1) * h | |
| 190 | 90 | x_ip1 = x_i + h | |
| 191 | 102 | ans = ans + quadrature_1d(n_quad, f, x_i, x_ip1) | |
| 192 | end do | ||
| 193 | 12 | end function | |
| 194 | |||
| 195 | !> @brief Composite numerical quadrature in 3D | ||
| 196 | !> | ||
| 197 | !> @param[in] n_quad Number of quadrature nodes per interval in each direction | ||
| 198 | !> @param[in] f Function to integrate | ||
| 199 | !> @param[in] x0,x1,y0,y1,z0,z1 Integration limits in the x, y and z directions | ||
| 200 | !> @param[in] nx,ny,nz Number of intervals in x, y and z | ||
| 201 | !> | ||
| 202 | !> @return Integral of f over the domain | ||
| 203 | 12 | pure real(wp) function quadrature_composite_3d(n_quad, f, x0, x1, y0, y1, z0, z1, nx, ny, nz) result(ans) | |
| 204 | implicit none | ||
| 205 | |||
| 206 | integer, intent(in) :: n_quad, nx, ny, nz | ||
| 207 | procedure(user_function_3d_interface) :: f | ||
| 208 | real(wp), intent(in) :: x0, x1, y0, y1, z0, z1 | ||
| 209 | |||
| 210 | integer :: i, j, k | ||
| 211 | real(wp) :: hx, hy, hz | ||
| 212 | real(wp) :: x_i, x_ip1, y_j, y_jp1, z_k, z_kp1 | ||
| 213 | |||
| 214 |
3/6✓ Branch 2 → 3 taken 12 times.
✗ Branch 2 → 5 not taken.
✓ Branch 3 → 4 taken 12 times.
✗ Branch 3 → 5 not taken.
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 12 times.
|
12 | if (nx < 1 .or. ny < 1 .or. nz < 1) then |
| 215 | ✗ | error stop "quadrature composite parameters n{xyz} incorrect" | |
| 216 | end if | ||
| 217 | |||
| 218 | 12 | hx = (x1 - x0) / nx | |
| 219 | 12 | hy = (y1 - y0) / ny | |
| 220 | 12 | hz = (z1 - z0) / nz | |
| 221 | |||
| 222 | ans = 0.0_wp | ||
| 223 |
2/2✓ Branch 7 → 8 taken 90 times.
✓ Branch 7 → 17 taken 12 times.
|
102 | do i = 1, nx |
| 224 | 90 | x_i = x0 + (i - 1) * hx | |
| 225 | 90 | x_ip1 = x_i + hx | |
| 226 |
2/2✓ Branch 9 → 10 taken 1020 times.
✓ Branch 9 → 15 taken 90 times.
|
1122 | do j = 1, ny |
| 227 | 1020 | y_j = y0 + (j - 1) * hy | |
| 228 | 1020 | y_jp1 = y_j + hy | |
| 229 |
2/2✓ Branch 11 → 12 taken 14040 times.
✓ Branch 11 → 13 taken 1020 times.
|
15150 | do k = 1, nz |
| 230 | 14040 | z_k = z0 + (k - 1) * hz | |
| 231 | 14040 | z_kp1 = z_k + hz | |
| 232 | 15060 | ans = ans + quadrature_3d(n_quad, f, x_i, x_ip1, y_j, y_jp1, z_k, z_kp1) | |
| 233 | end do | ||
| 234 | end do | ||
| 235 | end do | ||
| 236 | 12 | end function | |
| 237 | |||
| 238 | end module m_quadrature | ||
| 239 |