GCC Code Coverage Report

Directory: src/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 94.8% 128 / 0 / 135
Functions: 100.0% 11 / 0 / 11
Branches: 60.2% 112 / 0 / 186
bspline/m_bspline_linalg.f90
Line Branch Exec Source
1 !> The B-spline linear algebra module
2 !>
3 !> This module provides an inner product, L2-error computation, as well as an L2-projection and a Poisson solver for the 1D
4 !> B-spline basis functions.
5 module m_bspline_linalg
6 #include "petsc.fi"
7
8 use m_bspline_basis, only: BSplineSpace, BSplineFun
9 use m_common, only: wp, user_function_1d_interface
10
11 implicit none
12
13 private
14
15 public :: solve_poisson_problem, inner_product, l2_projection, l2_error, orthogonalize, biorthogonalize, integral
16
17 !> @brief Solve the Poisson problem
18 interface solve_poisson_problem
19 module procedure solve_poisson_problem_1d
20 end interface
21
22 !> @brief Compute the inner product of a user function with the basis functions
23 interface inner_product
24 module procedure inner_product_1d, inner_product_bsplines_1d
25 end interface
26
27 !> @brief Compute the L2-projection of a user function onto the basis functions
28 interface l2_projection
29 module procedure l2_projection_1d
30 end interface
31
32 !> @brief Compute the L2-norm of a B-spline function minus a user function
33 interface l2_error
34 module procedure l2_error_1d
35 end interface
36
37 interface integral
38 module procedure integral_1d
39 end interface
40
41 interface orthogonalize
42 module procedure orthogonalize_single_1d
43 end interface
44
45 interface biorthogonalize
46 module procedure biorthogonalize_single_1d
47 end interface
48
49 contains
50
51 !> @brief Compute the inner product of a user function with each of the basis functions of a given B-spline space
52 !>
53 !> @param[out] bfun The resulting B-spline function
54 !> @param[in] bspline The B-spline basis
55 !> @param[in] userfun The user function
56 !> @param[in] n_quad_extra _(optional)_ The additional number of quadrature points on top of those needed for exact integration of
57 !> the B-spline space
58 131 subroutine inner_product_1d(bfun, bspline, userfun, n_quad_extra)
59 use m_bspline_basis, only: BSplineSpace, size, get_internal_and_actual_inds
60 use m_bspline_quadrature, only: BSplineQuadrature, quadrature, n_quad_exact
61 implicit none
62
63 type(BSplineFun), intent(out) :: bfun
64 type(BSplineSpace), intent(in) :: bspline
65
2/2
✓ Branch 2 → 3 taken 131 times.
✓ Branch 2 → 4 taken 7693 times.
 7824 type(BSplineQuadrature) :: bspline_quad
66 procedure(user_function_1d_interface) :: userfun
67 integer, intent(in), optional :: n_quad_extra
68
69 integer :: ierr, i, j, j_internal, j_minus_i, ndx, n_quad
70 real(wp) :: val, x
71 real(wp), allocatable :: userfun_vals(:)
72 logical :: bspline_mirrored
73
74 7824 n_quad = n_quad_exact(bspline, bspline)
75
1/2
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 7824 times.
 7824 if (present(n_quad_extra)) then
76 ✗ n_quad = n_quad + n_quad_extra
77 end if
78
79
3/6
✓ Branch 8 → 9 taken 7824 times.
✗ Branch 8 → 10 not taken.
✓ Branch 10 → 11 taken 7824 times.
✗ Branch 10 → 12 not taken.
✓ Branch 12 → 13 taken 7824 times.
✗ Branch 12 → 14 not taken.
 7824 bspline_quad = BSplineQuadrature(bspline, n_quad)
80
81
6/12
✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 7824 times.
✓ Branch 16 → 15 taken 7824 times.
✗ Branch 16 → 17 not taken.
✓ Branch 17 → 18 taken 7824 times.
✗ Branch 17 → 19 not taken.
✓ Branch 19 → 20 taken 7824 times.
✗ Branch 19 → 21 not taken.
✗ Branch 21 → 22 not taken.
✓ Branch 21 → 23 taken 7824 times.
✗ Branch 23 → 24 not taken.
✓ Branch 23 → 25 taken 7824 times.
 23472 allocate (userfun_vals(1:n_quad))
82
83 7824 call bfun%init(bspline)
84
85
2/2
✓ Branch 27 → 28 taken 176302 times.
✓ Branch 27 → 45 taken 7824 times.
 184126 do i = 0, bspline%nr_intervals - 1
86
2/2
✓ Branch 28 → 29 taken 888422 times.
✓ Branch 28 → 31 taken 176302 times.
 1064724 do ndx = 1, n_quad
87 888422 x = (i + .5_wp + bspline_quad%nodes(ndx) / 2) * bspline_quad%bspline%h
88 1064724 userfun_vals(ndx) = userfun(x)
89 end do
90
2/2
✓ Branch 33 → 34 taken 888422 times.
✓ Branch 33 → 43 taken 176302 times.
 1072548 do j_minus_i = 0, bspline%degree
91 888422 call get_internal_and_actual_inds(bspline, i, j_minus_i, j, j_internal, bspline_mirrored)
92 val = 0._wp
93
2/2
✓ Branch 35 → 36 taken 55231 times.
✓ Branch 35 → 38 taken 833191 times.
 888422 if (bspline_mirrored) then
94
2/2
✓ Branch 36 → 37 taken 391511 times.
✓ Branch 36 → 41 taken 55231 times.
 446742 do ndx = 1, n_quad
95 val = val + bspline_quad%weights(ndx) * userfun_vals(ndx) * &
96 446742 & bspline_quad%bspline_at_nodes(j_internal, bspline%degree - j_minus_i)%data(n_quad + 1 - ndx)
97 end do
98 else
99
2/2
✓ Branch 38 → 39 taken 833191 times.
✓ Branch 38 → 40 taken 4398103 times.
 5231294 do ndx = 1, n_quad
100 val = val + bspline_quad%weights(ndx) * userfun_vals(ndx) * &
101 5231294 & bspline_quad%bspline_at_nodes(j_internal, j_minus_i)%data(ndx)
102 end do
103 end if
104
105 1953146 bfun%data(j) = bfun%data(j) + val
106 end do
107 end do
108
109
1/2
✗ Branch 46 → 47 not taken.
✓ Branch 46 → 48 taken 7824 times.
 7824 deallocate (userfun_vals)
110
3/6
✓ Branch 48 → 49 taken 7824 times.
✗ Branch 48 → 50 not taken.
✓ Branch 50 → 51 taken 7824 times.
✗ Branch 50 → 52 not taken.
✓ Branch 52 → 53 taken 7824 times.
✗ Branch 52 → 54 not taken.
 7824 end subroutine inner_product_1d
111
112 !> @brief Compute the L2-projection of a user function onto the basis functions of a given B-spline space
113 !>
114 !> @param[out] bfun The resulting B-spline function
115 !> @param[in] bspline The B-spline basis
116 !> @param[in] userfun The user function
117 !> @param[in] coordfun _(optional)_ The coordinate function
118 !> @param[in] rtol _(optional)_ The relative tolerance of the linear solver
119 !> @param[out] l2err _(optional)_ A cheap estimate of the L2-error
120 321 subroutine l2_projection_1d(bfun, bspline, userfun, coordfun, rtol, l2err)
121 use m_bspline_basis, only: BSplineSpace, get_petsc
122 use m_bspline_matrix, only: mass_matrix
123 use m_quadrature, only: quadrature
124
125 implicit none
126 type(BSplineFun), intent(out) :: bfun
127 type(BSplineSpace), intent(in) :: bspline
128 procedure(user_function_1d_interface) :: userfun
129 procedure(user_function_1d_interface), optional :: coordfun
130 real(wp), intent(in), optional :: rtol
131 real(wp), intent(out), optional :: l2err
132
133 Vec :: f, c
134 Mat :: M
135 KSP :: ksp
136 7668 type(BSplineFun) :: f_bfun
137 integer :: ierr
138 real(wp) :: rtol_
139
140
2/2
✓ Branch 4 → 5 taken 27 times.
✓ Branch 4 → 6 taken 7641 times.
 7668 if (present(rtol)) then
141 27 rtol_ = rtol
142 else
143 7641 rtol_ = 1.e-10_wp
144 end if
145
146
1/2
✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 7668 times.
 7668 if (present(coordfun)) then
147 ✗ call mass_matrix(M, bspline, bspline, coordfun)
148 else
149 7668 call mass_matrix(M, bspline, bspline)
150 end if
151 7668 call inner_product_1d(f_bfun, bspline, userfun)
152 7668 call get_petsc(f, f_bfun)
153
154 ! Create linear solver
155
1/2
✗ Branch 13 → 14 not taken.
✓ Branch 13 → 15 taken 7668 times.
 7668 PetscCall(KSPCreate(PETSC_COMM_SELF, ksp, ierr))
156
1/2
✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 7668 times.
 7668 PetscCall(KSPSetOperators(ksp, M, M, ierr))
157
1/2
✗ Branch 19 → 20 not taken.
✓ Branch 19 → 21 taken 7668 times.
 7668 PetscCall(KSPSetFromOptions(ksp, ierr))
158
1/2
✗ Branch 22 → 23 not taken.
✓ Branch 22 → 24 taken 7668 times.
 7668 PetscCall(KSPSetType(ksp, "cg", ierr))
159
160
1/2
✗ Branch 25 → 26 not taken.
✓ Branch 25 → 27 taken 7668 times.
 7668 PetscCall(KSPSetTolerances(ksp, rtol_, 1.e-30_wp, 1.e10_wp, 1000, ierr))
161
162
1/2
✗ Branch 28 → 29 not taken.
✓ Branch 28 → 30 taken 7668 times.
 7668 PetscCall(VecDuplicate(f, c, ierr))
163
1/2
✗ Branch 31 → 32 not taken.
✓ Branch 31 → 33 taken 7668 times.
 7668 PetscCall(KSPSolve(ksp, f, c, ierr))
164
165
2/2
✓ Branch 33 → 34 taken 480 times.
✓ Branch 33 → 38 taken 7188 times.
 7668 if (present(l2err)) then
166 ! TODO more accurate (quadrature per interval)
167 ! L2-error = c' * M * c - 2 * c' * f + int F^2 dx = int f^2 dx - c' * F
168
1/2
✗ Branch 35 → 36 not taken.
✓ Branch 35 → 37 taken 480 times.
 480 PetscCall(VecDot(c, f, l2err, ierr))
169 480 l2err = quadrature(1 + int(bspline%degree * 1.5), userfun_squared, 0._wp, 1._wp) - l2err
170 480 l2err = sqrt(abs(l2err))
171 end if
172
173 7668 call bfun%init(bspline, vec=c)
174
175
1/2
✗ Branch 40 → 41 not taken.
✓ Branch 40 → 42 taken 7668 times.
 7668 PetscCall(MatDestroy(M, ierr))
176
1/2
✗ Branch 43 → 44 not taken.
✓ Branch 43 → 45 taken 7668 times.
 7668 PetscCall(VecDestroy(f, ierr))
177
1/2
✗ Branch 46 → 47 not taken.
✓ Branch 46 → 48 taken 7668 times.
 7668 PetscCall(VecDestroy(c, ierr))
178
4/8
✓ Branch 2 → 3 taken 321 times.
✓ Branch 2 → 4 taken 7347 times.
✗ Branch 49 → 50 not taken.
✓ Branch 49 → 51 taken 7668 times.
✓ Branch 51 → 52 taken 7668 times.
✗ Branch 51 → 56 not taken.
✗ Branch 53 → 54 not taken.
✗ Branch 53 → 57 not taken.
 15336 PetscCall(KSPDestroy(ksp, ierr))
179
180 contains
181 3880 pure real(wp) function userfun_squared(x) result(f2)
182 implicit none
183
184 real(wp), intent(in) :: x
185
186 3880 f2 = userfun(x)**2
187 3880 end function userfun_squared
188 end subroutine
189
190 !> @brief Compute the L2-norm of a B-spline function minus a user function
191 !>
192 !> @param[in] bfun The B-spline function
193 !> @param[in] userfun The user function
194 !> @param[in] n_quad _(optional)_ The number of quadrature points
195 !>
196 !> @return The L2-error
197 830 real(wp) function l2_error_1d(bfun, userfun, n_quad) result(l2)
198 use m_bspline_basis, only: BSplineSpace, size, evaluate
199 use m_quadrature, only: quadrature
200
201 implicit none
202
203 type(BSplineFun), intent(in) :: bfun
204 procedure(user_function_1d_interface) :: userfun
205 integer, intent(in), optional :: n_quad
206
207 integer :: n_quad_, i
208
209
2/2
✓ Branch 2 → 3 taken 78 times.
✓ Branch 2 → 4 taken 752 times.
 830 if (present(n_quad)) then
210 78 n_quad_ = n_quad
211 else
212 752 n_quad_ = bfun%bspline%degree + 1
213 end if
214
215 830 l2 = sqrt(integral(error_fun, bfun%bspline, n_quad=n_quad_))
216
217 contains
218 85908 pure real(wp) function error_fun(x) result(err)
219 implicit none
220 real(wp), intent(in) :: x
221
222 85908 err = (userfun(x) - evaluate(bfun, x))**2
223 85908 end function error_fun
224
225 end function l2_error_1d
226
227 !> @brief Compute the integral of a user function over the domain [0, 1] using interval-wise quadrature
228 !>
229 !> @param[in] userfun The user function
230 !> @param[in] space The B-spline space defining the intervals
231 !> @param[in] n_quad _(optional)_ The number of quadrature points per interval
232 !>
233 2236 real(wp) function integral_1d(userfun, space, n_quad) result(int_val)
234 use m_bspline_basis, only: BSplineSpace
235 use m_quadrature, only: quadrature
236
237 implicit none
238
239 procedure(user_function_1d_interface) :: userfun
240 type(BSplineSpace), intent(in) :: space
241 integer, intent(in), optional :: n_quad
242
243 integer :: n_quad_, i
244
245
1/2
✓ Branch 2 → 3 taken 2236 times.
✗ Branch 2 → 4 not taken.
 2236 if (present(n_quad)) then
246 2236 n_quad_ = n_quad
247 else
248 ✗ n_quad_ = space%degree + 1
249 end if
250
251 int_val = 0._wp
252
2/2
✓ Branch 6 → 7 taken 26008 times.
✓ Branch 6 → 8 taken 2236 times.
 28244 do i = 0, space%nr_intervals - 1
253 28244 int_val = int_val + quadrature(n_quad_, userfun, i * space%h, (i + 1) * space%h)
254 end do
255 2236 end function integral_1d
256
257 !> @brief Compute the inner product of two B-splines functions
258 !>
259 !> @param[out] ip The resulting inner product
260 !> @param[in] bfun1 The first B-spline function
261 !> @param[in] bfun2 The second B-spline function
262 !> @param[in] userfun _(optional)_ The weight function defining the inner product
263 !> @param[in] n_quad_extra _(optional)_ The additional number of quadrature points on top of those needed for exact integration of
264 !> the B-spline spaces
265 1406 subroutine inner_product_bsplines_1d(ip, bfun1, bfun2, userfun, n_quad_extra)
266 use m_bspline_quadrature, only: n_quad_exact
267 implicit none
268
269 real(wp), intent(out) :: ip
270 type(BSplineFun), intent(in) :: bfun1, bfun2
271 procedure(user_function_1d_interface), optional :: userfun
272 integer, intent(in), optional :: n_quad_extra
273
274 integer :: n_quad
275
276 1406 n_quad = n_quad_exact(bfun1%bspline, bfun2%bspline)
277
1/2
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 1406 times.
 1406 if (present(n_quad_extra) .and. present(userfun)) then
278 ✗ n_quad = n_quad + n_quad_extra
279 end if
280
281 1406 ip = integral_1d(integrand, bfun1%bspline, n_quad=n_quad)
282 contains
283 69300 pure real(wp) function integrand(x) result(val)
284 use m_bspline_basis, only: evaluate
285 implicit none
286
287 real(wp), intent(in) :: x
288
289 69300 val = evaluate(bfun1, x) * evaluate(bfun2, x)
290
2/2
✓ Branch 2 → 3 taken 896 times.
✓ Branch 2 → 5 taken 68404 times.
 69300 if (present(userfun)) then
291 896 val = val * userfun(x)
292 end if
293
294 69300 end function integrand
295
296 end subroutine inner_product_bsplines_1d
297
298 !> @brief Solve the Poisson problem
299 !> \f$\nabla^2 u = f\f$
300 !> with Dirichlet boundary conditions
301 !>
302 !> @param[out] bfun The resulting B-spline function
303 !> @param[in] bspline The B-spline basis
304 !> @param[in] L The Laplace matrix
305 !> @param[in] f The right-hand side
306 !> @param[in] bdy_fun _(optional)_ The boundary function (homogeneous if not present)
307 !> @param[in] rtol _(optional)_ The relative tolerance of the linear solver
308 ✗ subroutine solve_poisson_problem_1d(bfun, bspline, L, f, bdy_fun, rtol)
309 use m_bspline_basis, only: BSplineSpace, size, get_petsc
310 use m_common, only: user_function_1d_interface
311 implicit none
312
313 type(BSplineFun), intent(out) :: bfun
314 type(BSplineSpace), intent(in) :: bspline
315 Mat, intent(in) :: L
316 type(BSplineFun), intent(in) :: f
317 procedure(user_function_1d_interface), optional :: bdy_fun
318 real(wp), intent(in), optional :: rtol
319
320 KSP :: ksp
321 Mat :: L_sub
322 Vec :: f_sub, c_sub, c, f_vec
323 IS :: is_sub
324 real(wp) :: val, val1(1)
325 integer :: ierr
326 integer, allocatable :: inds_sub(:)
327 integer :: i
328 real(wp) :: rtol_
329
330
1/2
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 156 times.
 156 if (present(rtol)) then
331 ✗ rtol_ = rtol
332 else
333 156 rtol_ = 1.e-10_wp
334 end if
335
336 156 call get_petsc(f_vec, f)
337
338 ! Impose Dirichlet BC
339
6/12
✗ Branch 8 → 9 not taken.
✓ Branch 8 → 10 taken 156 times.
✓ Branch 10 → 9 taken 156 times.
✗ Branch 10 → 11 not taken.
✓ Branch 11 → 12 taken 156 times.
✗ Branch 11 → 13 not taken.
✓ Branch 13 → 14 taken 156 times.
✗ Branch 13 → 15 not taken.
✗ Branch 15 → 16 not taken.
✓ Branch 15 → 17 taken 156 times.
✗ Branch 17 → 18 not taken.
✓ Branch 17 → 19 taken 156 times.
 468 allocate (inds_sub(1:size(bspline) - 2))
340
2/2
✓ Branch 20 → 21 taken 6060 times.
✓ Branch 20 → 22 taken 156 times.
 6216 do i = 1, size(bspline) - 2
341 6216 inds_sub(i) = i
342 end do
343
1/2
✗ Branch 24 → 25 not taken.
✓ Branch 24 → 26 taken 156 times.
 156 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size(bspline) - 2, inds_sub, PETSC_COPY_VALUES, is_sub, ierr))
344
1/2
✗ Branch 26 → 27 not taken.
✓ Branch 26 → 28 taken 156 times.
 156 deallocate (inds_sub)
345
1/2
✗ Branch 29 → 30 not taken.
✓ Branch 29 → 32 taken 156 times.
 156 PetscCall(MatCreateSubMatrix(L, is_sub, is_sub, MAT_INITIAL_MATRIX, L_sub, ierr))
346
347 ! Set boundary values
348
1/2
✗ Branch 33 → 34 not taken.
✓ Branch 33 → 35 taken 156 times.
 156 PetscCall(VecCreate(PETSC_COMM_SELF, c, ierr))
349
1/2
✗ Branch 36 → 37 not taken.
✓ Branch 36 → 38 taken 156 times.
 156 PetscCall(VecSetSizes(c, size(bspline), size(bspline), ierr))
350
1/2
✗ Branch 39 → 40 not taken.
✓ Branch 39 → 41 taken 156 times.
 156 PetscCall(VecSetFromOptions(c, ierr))
351
2/2
✓ Branch 42 → 43 taken 6372 times.
✓ Branch 42 → 55 taken 156 times.
 6528 do i = 0, size(bspline) - 1
352
2/2
✓ Branch 43 → 44 taken 156 times.
✓ Branch 43 → 46 taken 6216 times.
 6372 if (i == 0 .and. present(bdy_fun)) then
353 156 val = bdy_fun(0._wp)
354
3/4
✓ Branch 46 → 47 taken 156 times.
✓ Branch 46 → 50 taken 6060 times.
✓ Branch 47 → 48 taken 156 times.
✗ Branch 47 → 50 not taken.
 6216 else if (i == size(bspline) - 1 .and. present(bdy_fun)) then
355 156 val = bdy_fun(1._wp)
356 else
357 6060 val = 0._wp
358 end if
359
1/2
✗ Branch 52 → 53 not taken.
✓ Branch 52 → 54 taken 6372 times.
 6528 PetscCall(VecSetValue(c, i, val, INSERT_VALUES, ierr))
360 end do
361
362
1/2
✗ Branch 57 → 58 not taken.
✓ Branch 57 → 59 taken 156 times.
 156 PetscCall(VecAssemblyBegin(c, ierr))
363
1/2
✗ Branch 60 → 61 not taken.
✓ Branch 60 → 62 taken 156 times.
 156 PetscCall(VecAssemblyEnd(c, ierr))
364
365 ! Apply Dirichlet BC
366
1/2
✗ Branch 63 → 64 not taken.
✓ Branch 63 → 65 taken 156 times.
 156 PetscCall(VecScale(c, -1._wp, ierr))
367
1/2
✗ Branch 66 → 67 not taken.
✓ Branch 66 → 68 taken 156 times.
 156 PetscCall(MatMultAdd(L, c, f_vec, f_vec, ierr))
368
1/2
✗ Branch 69 → 70 not taken.
✓ Branch 69 → 71 taken 156 times.
 156 PetscCall(VecScale(c, -1._wp, ierr))
369
370
1/2
✗ Branch 72 → 73 not taken.
✓ Branch 72 → 74 taken 156 times.
 156 PetscCall(VecGetSubVector(f_vec, is_sub, f_sub, ierr))
371
372 ! TODO KSP should also be re-used; consider an object based on the matrix
373 ! TODO which has: the full matrix, the submatrix, the IS, the KSP (based on the submatrix)
374
375 ! Create linear solver
376
1/2
✗ Branch 75 → 76 not taken.
✓ Branch 75 → 77 taken 156 times.
 156 PetscCall(KSPCreate(PETSC_COMM_SELF, ksp, ierr))
377
1/2
✗ Branch 78 → 79 not taken.
✓ Branch 78 → 80 taken 156 times.
 156 PetscCall(KSPSetOperators(ksp, L_sub, L_sub, ierr))
378
1/2
✗ Branch 81 → 82 not taken.
✓ Branch 81 → 83 taken 156 times.
 156 PetscCall(KSPSetFromOptions(ksp, ierr))
379
1/2
✗ Branch 84 → 85 not taken.
✓ Branch 84 → 86 taken 156 times.
 156 PetscCall(KSPSetType(ksp, "cg", ierr))
380
381
1/2
✗ Branch 87 → 88 not taken.
✓ Branch 87 → 89 taken 156 times.
 156 PetscCall(KSPSetTolerances(ksp, rtol_, 1.e-30_wp, 1.e10_wp, 1000, ierr))
382
383
1/2
✗ Branch 90 → 91 not taken.
✓ Branch 90 → 92 taken 156 times.
 156 PetscCall(VecDuplicate(f_sub, c_sub, ierr))
384
1/2
✗ Branch 93 → 94 not taken.
✓ Branch 93 → 95 taken 156 times.
 156 PetscCall(KSPSolve(ksp, f_sub, c_sub, ierr))
385
386 ! Go back to full basis (with BC in place)
387
2/2
✓ Branch 96 → 97 taken 6060 times.
✓ Branch 96 → 107 taken 156 times.
 6216 do i = 1, size(bspline) - 2
388
3/4
✓ Branch 98 → 99 taken 6060 times.
✓ Branch 98 → 100 taken 6060 times.
✗ Branch 101 → 102 not taken.
✓ Branch 101 → 103 taken 6060 times.
 12120 PetscCall(VecGetValues(c_sub, 1, [i - 1], val1, ierr))
389
1/2
✗ Branch 104 → 105 not taken.
✓ Branch 104 → 106 taken 6060 times.
 6216 PetscCall(VecSetValue(c, i, val1(1), INSERT_VALUES, ierr))
390 end do
391
392
1/2
✗ Branch 109 → 110 not taken.
✓ Branch 109 → 111 taken 156 times.
 156 PetscCall(VecDestroy(f_sub, ierr))
393
1/2
✗ Branch 112 → 113 not taken.
✓ Branch 112 → 114 taken 156 times.
 156 PetscCall(VecDestroy(f_vec, ierr))
394
1/2
✗ Branch 115 → 116 not taken.
✓ Branch 115 → 117 taken 156 times.
 156 PetscCall(VecDestroy(c_sub, ierr))
395
396
1/2
✗ Branch 118 → 119 not taken.
✓ Branch 118 → 120 taken 156 times.
 156 PetscCall(MatDestroy(L_sub, ierr))
397
398 156 call bfun%init(bspline, vec=c)
399
1/2
✗ Branch 122 → 123 not taken.
✓ Branch 122 → 124 taken 156 times.
 156 PetscCall(VecDestroy(c, ierr))
400
1/2
✗ Branch 125 → 126 not taken.
✓ Branch 125 → 130 taken 156 times.
 156 PetscCall(KSPDestroy(ksp, ierr))
401
1/4
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 156 times.
✗ Branch 127 → 128 not taken.
✗ Branch 127 → 131 not taken.
 156 end subroutine solve_poisson_problem_1d
402
403 !> @brief Orthogonalize a B-spline function with respect to another B-spline function
404 !>
405 !> @param[inout] vout The orthogonalized B-spline function
406 !> @param[in] vin The B-spline function to orthogonalize against
407 !> @param[in] weight_fun _(optional)_ The weight function defining the inner product
408 560 subroutine orthogonalize_single_1d(vout, vin, weight_fun)
409 use m_bspline_basis, only: BSplineSpace, BSplineFun
410 implicit none
411
412 type(BSplineFun), intent(inout) :: vout
413 type(BSplineFun), intent(in) :: vin
414 procedure(user_function_1d_interface), optional :: weight_fun
415
416 real(wp) :: ip_vin_vout, ip_vin_vin, scale_factor
417
418 280 call inner_product_bsplines_1d(ip_vin_vout, vout, vin, weight_fun)
419 280 call inner_product_bsplines_1d(ip_vin_vin, vin, vin, weight_fun)
420 280 scale_factor = ip_vin_vout / ip_vin_vin
421
3/4
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 280 times.
✓ Branch 7 → 8 taken 3752 times.
✓ Branch 7 → 9 taken 280 times.
 4312 vout%data = vout%data - scale_factor * vin%data
422
423 280 end subroutine orthogonalize_single_1d
424
425 !> @brief Biorthogonalize a B-spline function with respect to another B-spline function
426 !>
427 !> vout <- vout - ( <vout, win> / <vin, win> ) * vin
428 !>
429 !> such that <vout, win> = 0
430 !>
431 !> @param[inout] vout The biorthogonalized B-spline function
432 !> @param[in] vin The B-spline function to biorthogonalize with
433 !> @param[in] win The B-spline function to biorthogonalize against
434 !> @param[in] weight_fun _(optional)_ The weight function defining the inner product
435 6 subroutine biorthogonalize_single_1d(vout, vin, win, weight_fun)
436 use m_bspline_basis, only: BSplineSpace, BSplineFun
437 implicit none
438
439 type(BSplineFun), intent(inout) :: vout
440 type(BSplineFun), intent(in) :: vin
441 type(BSplineFun), intent(in) :: win
442 procedure(user_function_1d_interface), optional :: weight_fun
443
444 real(wp) :: ip_win_vout, ip_vin_win, scale_factor
445
446 3 call inner_product_bsplines_1d(ip_win_vout, vout, win, weight_fun)
447 3 call inner_product_bsplines_1d(ip_vin_win, vin, win, weight_fun)
448
1/2
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 3 times.
 3 if (ip_vin_win == 0._wp) then
449 ✗ error stop 'biorthogonalize_single_1d: division by zero in inner product <vin, win>'
450 end if
451 3 scale_factor = ip_win_vout / ip_vin_win
452
3/4
✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 3 times.
✓ Branch 9 → 10 taken 107 times.
✓ Branch 9 → 11 taken 3 times.
 113 vout%data = vout%data - scale_factor * vin%data
453
454 3 end subroutine biorthogonalize_single_1d
455
456 end module
457