Line data Source code
1 : module linear_1d_operators
2 :
3 : ! This module provides the type "TriDiagOp" to represent operators on a 1D
4 : ! grid as tridiagonal matrices, and related types to represent boundary
5 : ! conditions.
6 : !
7 : ! The focus is on solving diffusion equations with a finite volume method
8 : ! in one dimension, but other utility operators are provided, e.g. a second
9 : ! order approximation to the first derivative.
10 : !
11 : ! In order to allow vectorization to occur, as well as to avoid unnecessary
12 : ! copying/reshaping of data in CAM, TriDiagOp actually represents a
13 : ! collection of independent operators that can be applied to collections of
14 : ! independent data; the innermost index is over independent systems (e.g.
15 : ! CAM columns).
16 : !
17 : ! A simple example:
18 : ! ! First derivative operator
19 : ! op = first_derivative(coords)
20 : ! ! Convert data to its derivative (extrapolate at boundaries).
21 : ! call op%apply(data)
22 : !
23 : ! With explicit boundary conditions:
24 : ! op = first_derivative(coords, &
25 : ! l_bndry=BoundaryFixedFlux(), &
26 : ! r_bndry=BoundaryFixedLayer(layer_distance))
27 : ! call op%apply(data, &
28 : ! l_cond=BoundaryFlux(flux, dt, thickness), &
29 : ! r_cond=BoundaryData(boundary))
30 : !
31 : ! Implicit solution example:
32 : ! ! Construct diffusion matrix.
33 : ! op = diffusion_operator(coords, d)
34 : ! call op%lmult_as_diag(-dt)
35 : ! call op%add_to_diag(1._r8)
36 : ! ! Decompose in order to invert the operation.
37 : ! decomp = TriDiagDecomp(op)
38 : ! ! Diffuse data for one time step (fixed flux boundaries).
39 : ! call decomp%left_div(data)
40 :
41 : use shr_kind_mod, only: r8 => shr_kind_r8
42 : use shr_log_mod, only: errMsg => shr_log_errMsg
43 : use shr_sys_mod, only: shr_sys_abort
44 : use coords_1d, only: Coords1D
45 :
46 : implicit none
47 : private
48 : save
49 :
50 : ! Main type.
51 : public :: TriDiagOp
52 : public :: operator(+)
53 : public :: operator(-)
54 :
55 : ! Decomposition used for inversion (left division).
56 : public :: TriDiagDecomp
57 :
58 : ! Multiplies by 0.
59 : public :: zero_operator
60 :
61 : ! Construct identity.
62 : public :: identity_operator
63 :
64 : ! Produce a TriDiagOp that is simply a diagonal matrix.
65 : public :: diagonal_operator
66 :
67 : ! For solving the diffusion-advection equation with implicit Euler.
68 : public :: diffusion_operator
69 : public :: advection_operator
70 :
71 : ! Derivatives accurate to second order on a non-uniform grid.
72 : public :: first_derivative
73 : public :: second_derivative
74 :
75 : ! Boundary condition types.
76 : public :: BoundaryType
77 : public :: BoundaryZero
78 : public :: BoundaryFirstOrder
79 : public :: BoundaryExtrapolate
80 : public :: BoundaryFixedLayer
81 : public :: BoundaryFixedFlux
82 :
83 : ! Boundary data types.
84 : public :: BoundaryCond
85 : public :: BoundaryNoData
86 : public :: BoundaryData
87 : public :: BoundaryFlux
88 :
89 : ! TriDiagOp represents operators that can work between nearest neighbors,
90 : ! with some extra logic at the boundaries. The implementation is a
91 : ! tridiagonal matrix plus boundary info.
92 : type :: TriDiagOp
93 : private
94 : ! The number of independent systems.
95 : integer, public :: nsys
96 : ! The size of the matrix (number of grid cells).
97 : integer, public :: ncel
98 : ! Super-, sub-, and regular diagonals.
99 : real(r8), allocatable :: spr(:,:)
100 : real(r8), allocatable :: sub(:,:)
101 : real(r8), allocatable :: diag(:,:)
102 : ! Buffers to hold boundary data; Details depend on the type of boundary
103 : ! being used.
104 : real(r8), allocatable :: left_bound(:)
105 : real(r8), allocatable :: right_bound(:)
106 : contains
107 : ! Applies the operator to a set of data.
108 : procedure :: apply => apply_tridiag
109 : ! Given the off-diagonal elements, fills in the diagonal so that the
110 : ! operator will have the constant function as an eigenvector with
111 : ! eigenvalue 0. This is used internally as a utility for construction of
112 : ! derivative operators.
113 : procedure :: deriv_diag => make_tridiag_deriv_diag
114 : ! Add/substract another tridiagonal from this one in-place (without
115 : ! creating a temporary object).
116 : procedure :: add => add_in_place_tridiag_ops
117 : procedure :: subtract => subtract_in_place_tridiag_ops
118 : ! Add input vector or scalar to the diagonal.
119 : procedure :: scalar_add_tridiag
120 : procedure :: diagonal_add_tridiag
121 : generic :: add_to_diag => scalar_add_tridiag, diagonal_add_tridiag
122 : ! Treat input vector (or scalar) as if it was the diagonal of an
123 : ! operator, and multiply this operator on the left by that value.
124 : procedure :: scalar_lmult_tridiag
125 : procedure :: diagonal_lmult_tridiag
126 : generic :: lmult_as_diag => &
127 : scalar_lmult_tridiag, diagonal_lmult_tridiag
128 : ! Deallocate and reset.
129 : procedure :: finalize => tridiag_finalize
130 : end type TriDiagOp
131 :
132 : interface operator(+)
133 : module procedure add_tridiag_ops
134 : end interface operator(+)
135 :
136 : interface operator(-)
137 : module procedure subtract_tridiag_ops
138 : end interface operator(-)
139 :
140 : interface TriDiagOp
141 : module procedure new_TriDiagOp
142 : end interface TriDiagOp
143 :
144 : !
145 : ! Boundary condition types for the operators.
146 : !
147 : ! Note that BoundaryFixedLayer and BoundaryFixedFlux are the only options
148 : ! supported for backwards operation (i.e. decomp%left_div). The others are
149 : ! meant for direct application only (e.g. to find a derivative).
150 : !
151 : ! BoundaryZero means that the operator fixes boundaries to 0.
152 : ! BoundaryFirstOrder means a one-sided approximation for the first
153 : ! derivative.
154 : ! BoundaryExtrapolate means that a second order approximation will be used,
155 : ! even at the boundaries. Boundary points do this by using their next-
156 : ! nearest neighbor to extrapolate.
157 : ! BoundaryFixedLayer means that there's an extra layer outside of the given
158 : ! grid, which must be specified when applying/inverting the operator.
159 : ! BoundaryFixedFlux is intended to provide a fixed-flux condition for
160 : ! typical advection/diffusion operators. It tweaks the edge condition
161 : ! to work on an input current rather than a value.
162 : !
163 : ! The different types were originally implemented through polymorphism, but
164 : ! PGI required this to be done via enum instead.
165 : integer, parameter :: zero_bndry = 0
166 : integer, parameter :: first_order_bndry = 1
167 : integer, parameter :: extrapolate_bndry = 2
168 : integer, parameter :: fixed_layer_bndry = 3
169 : integer, parameter :: fixed_flux_bndry = 4
170 :
171 : type :: BoundaryType
172 : private
173 : integer :: bndry_type = fixed_flux_bndry
174 : real(r8), allocatable :: edge_width(:)
175 : contains
176 : procedure :: make_left
177 : procedure :: make_right
178 : procedure :: finalize => boundary_type_finalize
179 : end type BoundaryType
180 :
181 : abstract interface
182 : subroutine deriv_seed(del_minus, del_plus, sub, spr)
183 : import :: r8
184 : real(r8), USE_CONTIGUOUS intent(in) :: del_minus(:)
185 : real(r8), USE_CONTIGUOUS intent(in) :: del_plus(:)
186 : real(r8), USE_CONTIGUOUS intent(out) :: sub(:)
187 : real(r8), USE_CONTIGUOUS intent(out) :: spr(:)
188 : end subroutine deriv_seed
189 : end interface
190 :
191 : interface BoundaryZero
192 : module procedure new_BoundaryZero
193 : end interface BoundaryZero
194 :
195 : interface BoundaryFirstOrder
196 : module procedure new_BoundaryFirstOrder
197 : end interface BoundaryFirstOrder
198 :
199 : interface BoundaryExtrapolate
200 : module procedure new_BoundaryExtrapolate
201 : end interface BoundaryExtrapolate
202 :
203 : interface BoundaryFixedLayer
204 : module procedure new_BoundaryFixedLayer
205 : end interface BoundaryFixedLayer
206 :
207 : interface BoundaryFixedFlux
208 : module procedure new_BoundaryFixedFlux
209 : end interface BoundaryFixedFlux
210 :
211 : !
212 : ! Data for boundary conditions themselves.
213 : !
214 : ! "No data" conditions perform extrapolation, if BoundaryExtrapolate was
215 : ! the boundary type used to construct the operator.
216 : !
217 : ! "Data" conditions contain extra data, which effectively extends the
218 : ! system with an extra cell.
219 : !
220 : ! "Flux" conditions contain prescribed fluxes.
221 : !
222 : ! The condition you can use depends on the boundary type from above that
223 : ! was used in the operator's construction. For BoundaryFixedLayer use
224 : ! BoundaryData. For BoundaryFixedFlux use BoundaryFlux. For everything
225 : ! else, use BoundaryNoData.
226 :
227 : ! The switches using this enumeration used to be unnecessary due to use of
228 : ! polymorphism, but this had to be backed off due to insufficient PGI
229 : ! support for type extension.
230 : integer, parameter :: no_data_cond = 0
231 : integer, parameter :: data_cond = 1
232 : integer, parameter :: flux_cond = 2
233 :
234 : type :: BoundaryCond
235 : private
236 : integer :: cond_type = no_data_cond
237 : real(r8), allocatable :: edge_data(:)
238 : contains
239 : procedure :: apply_left
240 : procedure :: apply_right
241 : procedure :: finalize => boundary_cond_finalize
242 : end type BoundaryCond
243 :
244 : ! Constructors for different types of BoundaryCond.
245 : interface BoundaryNoData
246 : module procedure new_BoundaryNoData
247 : end interface BoundaryNoData
248 :
249 : interface BoundaryData
250 : module procedure new_BoundaryData
251 : end interface BoundaryData
252 :
253 : interface BoundaryFlux
254 : module procedure new_BoundaryFlux
255 : end interface BoundaryFlux
256 :
257 : ! Opaque type to hold a tridiagonal matrix decomposition.
258 : !
259 : ! Method used is similar to Richtmyer and Morton (1967,pp 198-201), but
260 : ! the order of iteration is reversed, leading to A and C being swapped, and
261 : ! some differences in the indexing.
262 : type :: TriDiagDecomp
263 : private
264 : integer :: nsys = 0
265 : integer :: ncel = 0
266 : ! These correspond to A_k, E_k, and 1 / (B_k - A_k * E_{k+1})
267 : real(r8), allocatable :: ca(:,:)
268 : real(r8), allocatable :: ze(:,:)
269 : real(r8), allocatable :: dnom(:,:)
270 : contains
271 : procedure :: left_div => decomp_left_div
272 : procedure :: finalize => decomp_finalize
273 : end type TriDiagDecomp
274 :
275 : interface TriDiagDecomp
276 : module procedure new_TriDiagDecomp
277 : end interface TriDiagDecomp
278 :
279 : contains
280 :
281 : ! Operator that sets to 0.
282 0 : function zero_operator(nsys, ncel) result(op)
283 : ! Sizes for operator.
284 : integer, intent(in) :: nsys, ncel
285 :
286 : type(TriDiagOp) :: op
287 :
288 0 : op = TriDiagOp(nsys, ncel)
289 :
290 0 : op%spr = 0._r8
291 0 : op%sub = 0._r8
292 0 : op%diag = 0._r8
293 0 : op%left_bound = 0._r8
294 0 : op%right_bound = 0._r8
295 :
296 0 : end function zero_operator
297 :
298 : ! Operator that does nothing.
299 0 : function identity_operator(nsys, ncel) result(op)
300 : ! Sizes for operator.
301 : integer, intent(in) :: nsys, ncel
302 :
303 : type(TriDiagOp) :: op
304 :
305 0 : op = TriDiagOp(nsys, ncel)
306 :
307 0 : op%spr = 0._r8
308 0 : op%sub = 0._r8
309 0 : op%diag = 1._r8
310 0 : op%left_bound = 0._r8
311 0 : op%right_bound = 0._r8
312 :
313 0 : end function identity_operator
314 :
315 : ! Create an operator that just does an element-wise product by some data.
316 5956704 : function diagonal_operator(diag) result(op)
317 : ! Data to multiply by.
318 : real(r8), USE_CONTIGUOUS intent(in) :: diag(:,:)
319 :
320 : type(TriDiagOp) :: op
321 :
322 2978352 : op = TriDiagOp(size(diag, 1), size(diag, 2))
323 :
324 4578281136 : op%spr = 0._r8
325 4578281136 : op%sub = 0._r8
326 4630991040 : op%diag = diag
327 49731552 : op%left_bound = 0._r8
328 49731552 : op%right_bound = 0._r8
329 :
330 2978352 : end function diagonal_operator
331 :
332 : ! Diffusion matrix operator constructor. Given grid coordinates, a set of
333 : ! diffusion coefficients, and boundaries, creates a matrix corresponding
334 : ! to a finite volume representation of the operation:
335 : !
336 : ! d/dx (d_coef * d/dx)
337 : !
338 : ! This differs from what you would get from combining the first and second
339 : ! derivative operations, which would be more appropriate for a finite
340 : ! difference scheme that does not use grid cell averages.
341 50631984 : function diffusion_operator(coords, d_coef, l_bndry, r_bndry) &
342 : result(op)
343 : ! Grid cell locations.
344 : type(Coords1D), intent(in) :: coords
345 : ! Diffusion coefficient defined on interfaces.
346 : real(r8), USE_CONTIGUOUS intent(in) :: d_coef(:,:)
347 : ! Objects representing the kind of boundary on each side.
348 : class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry
349 : ! Output operator.
350 : type(TriDiagOp) :: op
351 :
352 : ! Selectors to implement default boundary.
353 : class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc
354 : ! Fixed flux is default, no allocation/deallocation needed.
355 25315992 : type(BoundaryType), target :: bndry_default
356 :
357 : ! Level index.
358 : integer :: k
359 :
360 25315992 : if (present(l_bndry)) then
361 1489176 : l_bndry_loc => l_bndry
362 : else
363 23826816 : l_bndry_loc => bndry_default
364 : end if
365 :
366 25315992 : if (present(r_bndry)) then
367 0 : r_bndry_loc => r_bndry
368 : else
369 25315992 : r_bndry_loc => bndry_default
370 : end if
371 :
372 : ! Allocate the operator.
373 25315992 : op = TriDiagOp(coords%n, coords%d)
374 :
375 : ! d_coef over the distance to the next cell gives you the matrix term for
376 : ! flux of material between cells. Dividing by cell thickness translates
377 : ! this to a tendency on the concentration. Hence the basic pattern is
378 : ! d_coef*rdst*rdel.
379 : !
380 : ! Boundary conditions for a fixed layer simply extend this by calculating
381 : ! the distance to the midpoint of the extra edge layer.
382 :
383 26805168 : select case (l_bndry_loc%bndry_type)
384 : case (fixed_layer_bndry)
385 0 : op%left_bound = 2._r8*d_coef(:,1)*coords%rdel(:,1) / &
386 26354952 : (l_bndry_loc%edge_width+coords%del(:,1))
387 : case default
388 399341592 : op%left_bound = 0._r8
389 : end select
390 :
391 2240922219 : do k = 1, coords%d-1
392 36996435438 : op%spr(:,k) = d_coef(:,k+1)*coords%rdst(:,k)*coords%rdel(:,k)
393 37021751430 : op%sub(:,k) = d_coef(:,k+1)*coords%rdst(:,k)*coords%rdel(:,k+1)
394 : end do
395 :
396 25315992 : select case (r_bndry_loc%bndry_type)
397 : case (fixed_layer_bndry)
398 0 : op%right_bound = 2._r8*d_coef(:,coords%d+1)*coords%rdel(:,coords%d) / &
399 0 : (r_bndry_loc%edge_width+coords%del(:,coords%d))
400 : case default
401 422718192 : op%right_bound = 0._r8
402 : end select
403 :
404 : ! Above, we found all off-diagonals. Now get the diagonal.
405 25315992 : call op%deriv_diag()
406 :
407 50631984 : end function diffusion_operator
408 :
409 : ! Advection matrix operator constructor. Similar to diffusion_operator, it
410 : ! constructs an operator A corresponding to:
411 : !
412 : ! A y = d/dx (-v_coef * y)
413 : !
414 : ! Again, this is targeted at representing this operator acting on grid-cell
415 : ! averages in a finite volume scheme, rather than a literal representation.
416 0 : function advection_operator(coords, v_coef, l_bndry, r_bndry) &
417 : result(op)
418 : ! Grid cell locations.
419 : type(Coords1D), intent(in) :: coords
420 : ! Advection coefficient (effective velocity).
421 : real(r8), USE_CONTIGUOUS intent(in) :: v_coef(:,:)
422 : ! Objects representing the kind of boundary on each side.
423 : class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry
424 : ! Output operator.
425 : type(TriDiagOp) :: op
426 :
427 : ! Selectors to implement default boundary.
428 : class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc
429 : ! Fixed flux is default, no allocation/deallocation needed.
430 0 : type(BoundaryType), target :: bndry_default
431 :
432 : ! Negative derivative of v.
433 0 : real(r8) :: v_deriv(coords%n,coords%d)
434 :
435 0 : if (present(l_bndry)) then
436 0 : l_bndry_loc => l_bndry
437 : else
438 0 : l_bndry_loc => bndry_default
439 : end if
440 :
441 0 : if (present(r_bndry)) then
442 0 : r_bndry_loc => r_bndry
443 : else
444 0 : r_bndry_loc => bndry_default
445 : end if
446 :
447 : ! Allocate the operator.
448 0 : op = TriDiagOp(coords%n, coords%d)
449 :
450 : ! Construct the operator in two stages using the product rule. First
451 : ! create (-v * d/dx), then -dv/dx, and add the two.
452 : !
453 : ! For the first part, we want to interpolate to interfaces (weighted
454 : ! average involving del/2*dst), multiply by -v to get flux, then divide
455 : ! by cell thickness, which gives a concentration tendency:
456 : !
457 : ! (del/(2*dst))*(-v_coef)/del
458 : !
459 : ! Simplifying gives -v_coef*rdst*0.5, as seen below.
460 :
461 0 : select case (l_bndry_loc%bndry_type)
462 : case (fixed_layer_bndry)
463 0 : op%left_bound = v_coef(:,1) / &
464 0 : (l_bndry_loc%edge_width+coords%del(:,1))
465 : case default
466 0 : op%left_bound = 0._r8
467 : end select
468 :
469 0 : op%sub = v_coef(:,2:coords%d)*coords%rdst*0.5_r8
470 0 : op%spr = -op%sub
471 :
472 0 : select case (r_bndry_loc%bndry_type)
473 : case (fixed_layer_bndry)
474 0 : op%right_bound = v_coef(:,coords%d+1) / &
475 0 : (r_bndry_loc%edge_width+coords%del(:,coords%d))
476 : case default
477 0 : op%right_bound = 0._r8
478 : end select
479 :
480 : ! Above, we found all off-diagonals. Now get the diagonal. This must be
481 : ! done at this specific point, since the other half of the operator is
482 : ! not "derivative-like" in the sense of yielding 0 for a constant input.
483 0 : call op%deriv_diag()
484 :
485 : ! The second half of the operator simply involves taking a first-order
486 : ! derivative of v. Since v is on the interfaces, just use:
487 : ! (v(k+1) - v(k))*rdel(k)
488 0 : v_deriv(:,1) = v_coef(:,2)*coords%rdel(:,1)
489 :
490 0 : select case (l_bndry_loc%bndry_type)
491 : case (fixed_layer_bndry)
492 0 : v_deriv(:,1) = v_deriv(:,1) - v_coef(:,1)*coords%rdel(:,1)
493 : end select
494 :
495 0 : v_deriv(:,2:coords%d-1) = (v_coef(:,3:coords%d) - &
496 0 : v_coef(:,2:coords%d-1))*coords%rdel(:,2:coords%d-1)
497 :
498 0 : v_deriv(:,coords%d) = -v_coef(:,coords%d)*coords%rdel(:,coords%d)
499 :
500 0 : select case (r_bndry_loc%bndry_type)
501 : case (fixed_layer_bndry)
502 : v_deriv(:,coords%d) = v_deriv(:,coords%d) &
503 0 : + v_coef(:,coords%d+1)*coords%del(:,coords%d)
504 : end select
505 :
506 : ! Combine the two pieces.
507 0 : op%diag = op%diag - v_deriv
508 :
509 0 : end function advection_operator
510 :
511 : ! Second order approximation to the first and second derivatives on a non-
512 : ! uniform grid.
513 : !
514 : ! Both operators are constructed with the same method, except for a "seed"
515 : ! function that takes local distances between points to create the
516 : ! off-diagonal terms.
517 0 : function first_derivative(grid_spacing, l_bndry, r_bndry) result(op)
518 : ! Distances between points.
519 : real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:)
520 : ! Boundary conditions.
521 : class(BoundaryType), intent(in), optional :: l_bndry, r_bndry
522 : ! Output operator.
523 : type(TriDiagOp) :: op
524 :
525 : op = deriv_op_from_seed(grid_spacing, first_derivative_seed, &
526 0 : l_bndry, r_bndry)
527 :
528 0 : end function first_derivative
529 :
530 0 : subroutine first_derivative_seed(del_minus, del_plus, sub, spr)
531 : ! Distances to next and previous point.
532 : real(r8), USE_CONTIGUOUS intent(in) :: del_minus(:)
533 : real(r8), USE_CONTIGUOUS intent(in) :: del_plus(:)
534 : ! Off-diagonal matrix terms.
535 : real(r8), USE_CONTIGUOUS intent(out) :: sub(:)
536 : real(r8), USE_CONTIGUOUS intent(out) :: spr(:)
537 :
538 0 : real(r8) :: del_sum(size(del_plus))
539 :
540 0 : del_sum = del_plus + del_minus
541 :
542 0 : sub = - del_plus / (del_minus*del_sum)
543 0 : spr = del_minus / (del_plus*del_sum)
544 :
545 0 : end subroutine first_derivative_seed
546 :
547 0 : function second_derivative(grid_spacing, l_bndry, r_bndry) result(op)
548 : ! Distances between points.
549 : real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:)
550 : ! Boundary conditions.
551 : class(BoundaryType), intent(in), optional :: l_bndry, r_bndry
552 : ! Output operator.
553 : type(TriDiagOp) :: op
554 :
555 : op = deriv_op_from_seed(grid_spacing, second_derivative_seed, &
556 0 : l_bndry, r_bndry)
557 :
558 0 : end function second_derivative
559 :
560 0 : subroutine second_derivative_seed(del_minus, del_plus, sub, spr)
561 : ! Distances to next and previous point.
562 : real(r8), USE_CONTIGUOUS intent(in) :: del_minus(:)
563 : real(r8), USE_CONTIGUOUS intent(in) :: del_plus(:)
564 : ! Off-diagonal matrix terms.
565 : real(r8), USE_CONTIGUOUS intent(out) :: sub(:)
566 : real(r8), USE_CONTIGUOUS intent(out) :: spr(:)
567 :
568 0 : real(r8) :: del_sum(size(del_plus))
569 :
570 0 : del_sum = del_plus + del_minus
571 :
572 0 : sub = 2._r8 / (del_minus*del_sum)
573 0 : spr = 2._r8 / (del_plus*del_sum)
574 :
575 0 : end subroutine second_derivative_seed
576 :
577 : ! Brains behind the first/second derivative functions.
578 0 : function deriv_op_from_seed(grid_spacing, seed, l_bndry, r_bndry) result(op)
579 : ! Distances between points.
580 : real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:)
581 : ! Function to locally construct matrix elements.
582 : procedure(deriv_seed) :: seed
583 : ! Boundary conditions.
584 : class(BoundaryType), target, intent(in), optional :: l_bndry, r_bndry
585 : ! Output operator.
586 : type(TriDiagOp) :: op
587 :
588 : ! Selectors to implement default boundary.
589 : class(BoundaryType), pointer :: l_bndry_loc, r_bndry_loc
590 : ! Fixed flux is default, no allocation/deallocation needed.
591 0 : type(BoundaryType), target :: bndry_default
592 :
593 : integer :: k
594 :
595 0 : if (present(l_bndry)) then
596 0 : l_bndry_loc => l_bndry
597 : else
598 0 : l_bndry_loc => bndry_default
599 : end if
600 :
601 0 : if (present(r_bndry)) then
602 0 : r_bndry_loc => r_bndry
603 : else
604 0 : r_bndry_loc => bndry_default
605 : end if
606 :
607 : ! Number of grid points is one greater than the spacing.
608 0 : op = TriDiagOp(size(grid_spacing, 1), size(grid_spacing, 2) + 1)
609 :
610 : ! Left boundary condition.
611 : call l_bndry_loc%make_left(grid_spacing, seed, &
612 0 : op%left_bound, op%spr(:,1))
613 :
614 0 : do k = 2, op%ncel-1
615 : call seed(grid_spacing(:,k-1), grid_spacing(:,k), &
616 0 : op%sub(:,k-1), op%spr(:,k))
617 : end do
618 :
619 : ! Right boundary condition.
620 : call r_bndry_loc%make_right(grid_spacing, seed, &
621 0 : op%sub(:,op%ncel-1), op%right_bound)
622 :
623 : ! Above, we found all off-diagonals. Now get the diagonal.
624 0 : call op%deriv_diag()
625 :
626 0 : end function deriv_op_from_seed
627 :
628 : ! Boundary constructors. Most simply set an internal flag, but
629 : ! BoundaryFixedLayer accepts an argument representing the distance to the
630 : ! location where the extra layer is defined.
631 :
632 0 : function new_BoundaryZero() result(new_bndry)
633 : type(BoundaryType) :: new_bndry
634 :
635 0 : new_bndry%bndry_type = zero_bndry
636 :
637 0 : end function new_BoundaryZero
638 :
639 0 : function new_BoundaryFirstOrder() result(new_bndry)
640 : type(BoundaryType) :: new_bndry
641 :
642 0 : new_bndry%bndry_type = first_order_bndry
643 :
644 0 : end function new_BoundaryFirstOrder
645 :
646 0 : function new_BoundaryExtrapolate() result(new_bndry)
647 : type(BoundaryType) :: new_bndry
648 :
649 0 : new_bndry%bndry_type = extrapolate_bndry
650 :
651 0 : end function new_BoundaryExtrapolate
652 :
653 2978352 : function new_BoundaryFixedLayer(width) result(new_bndry)
654 : real(r8), USE_CONTIGUOUS intent(in) :: width(:)
655 : type(BoundaryType) :: new_bndry
656 :
657 1489176 : new_bndry%bndry_type = fixed_layer_bndry
658 24865776 : new_bndry%edge_width = width
659 :
660 1489176 : end function new_BoundaryFixedLayer
661 :
662 0 : function new_BoundaryFixedFlux() result(new_bndry)
663 : type(BoundaryType) :: new_bndry
664 :
665 : new_bndry%bndry_type = fixed_flux_bndry
666 :
667 0 : end function new_BoundaryFixedFlux
668 :
669 : ! The make_left and make_right methods implement the boundary conditions
670 : ! using an input seed.
671 :
672 0 : subroutine make_left(self, grid_spacing, seed, term1, term2)
673 : class(BoundaryType), intent(in) :: self
674 : real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:)
675 : procedure(deriv_seed) :: seed
676 : real(r8), USE_CONTIGUOUS intent(out) :: term1(:)
677 : real(r8), USE_CONTIGUOUS intent(out) :: term2(:)
678 :
679 0 : real(r8) :: del_plus(size(term1)), del_minus(size(term1))
680 :
681 0 : select case (self%bndry_type)
682 : case (zero_bndry)
683 0 : term1 = 0._r8
684 0 : term2 = 0._r8
685 : case (first_order_bndry)
686 : ! To calculate to first order, just use a really huge del_minus (i.e.
687 : ! pretend that there's a point so far away it doesn't matter).
688 0 : del_plus = grid_spacing(:,1)
689 0 : del_minus = del_plus * 4._r8 / epsilon(1._r8)
690 0 : call seed(del_minus, del_plus, term1, term2)
691 : case (extrapolate_bndry)
692 : ! To extrapolate from the boundary, use distance from the nearest
693 : ! neighbor (as usual) and the second nearest neighbor (with a negative
694 : ! sign, since we are using two points on the same side).
695 0 : del_plus = grid_spacing(:,1)
696 0 : del_minus = - (grid_spacing(:,1) + grid_spacing(:,2))
697 0 : call seed(del_minus, del_plus, term1, term2)
698 : case (fixed_layer_bndry)
699 : ! Use edge value to extend the grid.
700 0 : del_plus = grid_spacing(:,1)
701 0 : del_minus = self%edge_width
702 0 : call seed(del_minus, del_plus, term1, term2)
703 : case (fixed_flux_bndry)
704 : ! Treat grid as uniform, but then zero out the contribution from data
705 : ! on one side (since it will be prescribed).
706 0 : del_plus = grid_spacing(:,1)
707 0 : del_minus = del_plus
708 0 : call seed(del_minus, del_plus, term1, term2)
709 0 : term1 = 0._r8
710 : case default
711 : call shr_sys_abort("Invalid boundary type at "// &
712 0 : errMsg(__FILE__, __LINE__))
713 : end select
714 :
715 0 : end subroutine make_left
716 :
717 0 : subroutine make_right(self, grid_spacing, seed, term1, term2)
718 : class(BoundaryType), intent(in) :: self
719 : real(r8), USE_CONTIGUOUS intent(in) :: grid_spacing(:,:)
720 : procedure(deriv_seed) :: seed
721 : real(r8), USE_CONTIGUOUS intent(out) :: term1(:)
722 : real(r8), USE_CONTIGUOUS intent(out) :: term2(:)
723 :
724 0 : real(r8) :: del_plus(size(term1)), del_minus(size(term1))
725 :
726 0 : select case (self%bndry_type)
727 : case (zero_bndry)
728 0 : term1 = 0._r8
729 0 : term2 = 0._r8
730 : case (first_order_bndry)
731 : ! Use huge del_plus, analogous to how left boundary works.
732 0 : del_minus = grid_spacing(:,size(grid_spacing, 2))
733 0 : del_plus = del_minus * 4._r8 / epsilon(1._r8)
734 0 : call seed(del_minus, del_plus, term1, term2)
735 : case (extrapolate_bndry)
736 : ! Same strategy as left boundary, but reversed.
737 0 : del_plus = - (grid_spacing(:,size(grid_spacing, 2) - 1) + &
738 0 : grid_spacing(:,size(grid_spacing, 2)))
739 0 : del_minus = grid_spacing(:,size(grid_spacing, 2))
740 0 : call seed(del_minus, del_plus, term1, term2)
741 : case (fixed_layer_bndry)
742 : ! Use edge value to extend the grid.
743 0 : del_plus = self%edge_width
744 0 : del_minus = grid_spacing(:,size(grid_spacing, 2))
745 0 : call seed(del_minus, del_plus, term1, term2)
746 : case (fixed_flux_bndry)
747 : ! Uniform grid, but with edge zeroed.
748 0 : del_plus = grid_spacing(:,size(grid_spacing, 2))
749 0 : del_minus = del_plus
750 0 : call seed(del_minus, del_plus, term1, term2)
751 0 : term2 = 0._r8
752 : case default
753 : call shr_sys_abort("Invalid boundary type at "// &
754 0 : errMsg(__FILE__, __LINE__))
755 : end select
756 :
757 0 : end subroutine make_right
758 :
759 0 : subroutine boundary_type_finalize(self)
760 : class(BoundaryType), intent(inout) :: self
761 :
762 0 : self%bndry_type = fixed_flux_bndry
763 0 : if (allocated(self%edge_width)) deallocate(self%edge_width)
764 :
765 0 : end subroutine boundary_type_finalize
766 :
767 : ! Constructor for TriDiagOp; this just sets the size and allocates
768 : ! arrays.
769 28294344 : type(TriDiagOp) function new_TriDiagOp(nsys, ncel)
770 :
771 : integer, intent(in) :: nsys, ncel
772 :
773 28294344 : new_TriDiagOp%nsys = nsys
774 28294344 : new_TriDiagOp%ncel = ncel
775 :
776 : allocate(new_TriDiagOp%spr(nsys,ncel-1), &
777 : new_TriDiagOp%sub(nsys,ncel-1), &
778 : new_TriDiagOp%diag(nsys,ncel), &
779 : new_TriDiagOp%left_bound(nsys), &
780 339532128 : new_TriDiagOp%right_bound(nsys))
781 :
782 28294344 : end function new_TriDiagOp
783 :
784 : ! Deallocator for TriDiagOp.
785 50631984 : subroutine tridiag_finalize(self)
786 : class(TriDiagOp), intent(inout) :: self
787 :
788 50631984 : self%nsys = 0
789 50631984 : self%ncel = 0
790 :
791 50631984 : if (allocated(self%spr)) deallocate(self%spr)
792 50631984 : if (allocated(self%sub)) deallocate(self%sub)
793 50631984 : if (allocated(self%diag)) deallocate(self%diag)
794 50631984 : if (allocated(self%left_bound)) deallocate(self%left_bound)
795 50631984 : if (allocated(self%right_bound)) deallocate(self%right_bound)
796 :
797 50631984 : end subroutine tridiag_finalize
798 :
799 : ! Boundary condition constructors.
800 :
801 0 : function new_BoundaryNoData() result(new_cond)
802 : type(BoundaryCond) :: new_cond
803 :
804 : new_cond%cond_type = no_data_cond
805 : ! No edge data, so leave it unallocated.
806 :
807 0 : end function new_BoundaryNoData
808 :
809 2978352 : function new_BoundaryData(data) result(new_cond)
810 : real(r8), USE_CONTIGUOUS intent(in) :: data(:)
811 : type(BoundaryCond) :: new_cond
812 :
813 1489176 : new_cond%cond_type = data_cond
814 24865776 : new_cond%edge_data = data
815 :
816 1489176 : end function new_BoundaryData
817 :
818 0 : function new_BoundaryFlux(flux, dt, spacing) result(new_cond)
819 : real(r8), USE_CONTIGUOUS intent(in) :: flux(:)
820 : real(r8), intent(in) :: dt
821 : real(r8), USE_CONTIGUOUS intent(in) :: spacing(:)
822 : type(BoundaryCond) :: new_cond
823 :
824 0 : new_cond%cond_type = flux_cond
825 0 : new_cond%edge_data = flux*dt/spacing
826 :
827 0 : end function new_BoundaryFlux
828 :
829 : ! Application of input data.
830 : !
831 : ! When no data is input, assume that any bound term is applied to the
832 : ! third element in from the edge for extrapolation. Boundary conditions
833 : ! that don't need any edge data at all can then simply set the boundary
834 : ! terms to 0.
835 :
836 1489176 : function apply_left(self, bound_term, array) result(delta_edge)
837 : class(BoundaryCond), intent(in) :: self
838 : real(r8), USE_CONTIGUOUS intent(in) :: bound_term(:)
839 : real(r8), USE_CONTIGUOUS intent(in) :: array(:,:)
840 : real(r8) :: delta_edge(size(array, 1))
841 :
842 1489176 : select case (self%cond_type)
843 : case (no_data_cond)
844 0 : delta_edge = bound_term*array(:,3)
845 : case (data_cond)
846 24865776 : delta_edge = bound_term*self%edge_data
847 : case (flux_cond)
848 0 : delta_edge = self%edge_data
849 : case default
850 : call shr_sys_abort("Invalid boundary condition at "// &
851 1489176 : errMsg(__FILE__, __LINE__))
852 : end select
853 :
854 1489176 : end function apply_left
855 :
856 0 : function apply_right(self, bound_term, array) result(delta_edge)
857 : class(BoundaryCond), intent(in) :: self
858 : real(r8), USE_CONTIGUOUS intent(in) :: bound_term(:)
859 : real(r8), USE_CONTIGUOUS intent(in) :: array(:,:)
860 : real(r8) :: delta_edge(size(array, 1))
861 :
862 0 : select case (self%cond_type)
863 : case (no_data_cond)
864 0 : delta_edge = bound_term*array(:,size(array, 2)-2)
865 : case (data_cond)
866 0 : delta_edge = bound_term*self%edge_data
867 : case (flux_cond)
868 0 : delta_edge = self%edge_data
869 : case default
870 : call shr_sys_abort("Invalid boundary condition at "// &
871 0 : errMsg(__FILE__, __LINE__))
872 : end select
873 :
874 0 : end function apply_right
875 :
876 0 : subroutine boundary_cond_finalize(self)
877 : class(BoundaryCond), intent(inout) :: self
878 :
879 0 : self%cond_type = no_data_cond
880 0 : if (allocated(self%edge_data)) deallocate(self%edge_data)
881 :
882 0 : end subroutine boundary_cond_finalize
883 :
884 : ! Apply an operator and return the new data.
885 0 : function apply_tridiag(self, array, l_cond, r_cond) result(output)
886 : ! Operator to apply.
887 : class(TriDiagOp), intent(in) :: self
888 : ! Data to act on.
889 : real(r8), USE_CONTIGUOUS intent(in) :: array(:,:)
890 : ! Objects representing boundary conditions.
891 : class(BoundaryCond), target, intent(in), optional :: l_cond, r_cond
892 : ! Function result.
893 : real(r8) :: output(size(array, 1), size(array, 2))
894 :
895 : ! Local objects to implement default.
896 : class(BoundaryCond), pointer :: l_cond_loc, r_cond_loc
897 : ! Default state is no data, no allocation/deallocation needed.
898 0 : type(BoundaryCond), target :: cond_default
899 :
900 : ! Level index.
901 : integer :: k
902 :
903 0 : if (present(l_cond)) then
904 0 : l_cond_loc => l_cond
905 : else
906 0 : l_cond_loc => cond_default
907 : end if
908 :
909 0 : if (present(r_cond)) then
910 0 : r_cond_loc => r_cond
911 : else
912 0 : r_cond_loc => cond_default
913 : end if
914 :
915 : ! Left boundary.
916 0 : output(:,1) = self%diag(:,1)*array(:,1) + &
917 0 : self%spr(:,1)*array(:,2) + &
918 0 : l_cond_loc%apply_left(self%left_bound, array)
919 :
920 0 : do k = 2, self%ncel-1
921 0 : output(:,k) = &
922 0 : self%sub(:,k-1)*array(:,k-1) + &
923 0 : self%diag(:,k)*array(:,k ) + &
924 0 : self%spr(:,k)*array(:,k+1)
925 : end do
926 :
927 : ! Right boundary.
928 0 : output(:,self%ncel) = &
929 0 : self%sub(:,self%ncel-1)*array(:,self%ncel-1) + &
930 0 : self%diag(:,self%ncel)*array(:,self%ncel) + &
931 0 : r_cond_loc%apply_right(self%right_bound, array)
932 :
933 0 : end function apply_tridiag
934 :
935 : ! Fill in the diagonal for a TriDiagOp for a derivative operator, where
936 : ! the off diagonal elements are already filled in.
937 25315992 : subroutine make_tridiag_deriv_diag(self)
938 :
939 : class(TriDiagOp), intent(inout) :: self
940 :
941 : ! If a derivative operator operates on a constant function, it must
942 : ! return 0 everywhere. To force this, make sure that all rows add to
943 : ! zero in the matrix.
944 37021751430 : self%diag(:,:self%ncel-1) = - self%spr
945 422718192 : self%diag(:,self%ncel) = - self%right_bound
946 422718192 : self%diag(:,1) = self%diag(:,1) - self%left_bound
947 37021751430 : self%diag(:,2:) = self%diag(:,2:) - self%sub
948 :
949 25315992 : end subroutine make_tridiag_deriv_diag
950 :
951 : ! Sum two TriDiagOp objects into a new one; this is just the addition of
952 : ! all the entries.
953 0 : function add_tridiag_ops(op1, op2) result(new_op)
954 :
955 : type(TriDiagOp), intent(in) :: op1, op2
956 : type(TriDiagOp) :: new_op
957 :
958 0 : new_op = op1
959 :
960 0 : call new_op%add(op2)
961 :
962 0 : end function add_tridiag_ops
963 :
964 2978352 : subroutine add_in_place_tridiag_ops(self, other)
965 :
966 : class(TriDiagOp), intent(inout) :: self
967 : class(TriDiagOp), intent(in) :: other
968 :
969 4581259488 : self%spr = self%spr + other%spr
970 4581259488 : self%sub = self%sub + other%sub
971 4630991040 : self%diag = self%diag + other%diag
972 :
973 52709904 : self%left_bound = self%left_bound + other%left_bound
974 52709904 : self%right_bound = self%right_bound + other%right_bound
975 :
976 2978352 : end subroutine add_in_place_tridiag_ops
977 :
978 : ! Subtract two TriDiagOp objects.
979 0 : function subtract_tridiag_ops(op1, op2) result(new_op)
980 :
981 : type(TriDiagOp), intent(in) :: op1, op2
982 : type(TriDiagOp) :: new_op
983 :
984 0 : new_op = op1
985 :
986 0 : call new_op%subtract(op2)
987 :
988 0 : end function subtract_tridiag_ops
989 :
990 : ! Subtract two TriDiagOp objects.
991 0 : subroutine subtract_in_place_tridiag_ops(self, other)
992 :
993 : class(TriDiagOp), intent(inout) :: self
994 : class(TriDiagOp), intent(in) :: other
995 :
996 0 : self%spr = self%spr - other%spr
997 0 : self%sub = self%sub - other%sub
998 0 : self%diag = self%diag - other%diag
999 :
1000 0 : self%left_bound = self%left_bound - other%left_bound
1001 0 : self%right_bound = self%right_bound - other%right_bound
1002 :
1003 0 : end subroutine subtract_in_place_tridiag_ops
1004 :
1005 : ! Equivalent to adding a multiple of the identity.
1006 25315992 : subroutine scalar_add_tridiag(self, constant)
1007 :
1008 : class(TriDiagOp), intent(inout) :: self
1009 : real(r8), intent(in) :: constant
1010 :
1011 37444469622 : self%diag = self%diag + constant
1012 :
1013 25315992 : end subroutine scalar_add_tridiag
1014 :
1015 : ! Equivalent to adding the diagonal operator constructed from diag_array.
1016 0 : subroutine diagonal_add_tridiag(self, diag_array)
1017 :
1018 : class(TriDiagOp), intent(inout) :: self
1019 : real(r8), USE_CONTIGUOUS intent(in) :: diag_array(:,:)
1020 :
1021 0 : self%diag = self%diag + diag_array
1022 :
1023 0 : end subroutine diagonal_add_tridiag
1024 :
1025 : ! Multiply a scalar by an array.
1026 25315992 : subroutine scalar_lmult_tridiag(self, constant)
1027 :
1028 : class(TriDiagOp), intent(inout) :: self
1029 : real(r8), intent(in) :: constant
1030 :
1031 37021751430 : self%spr = self%spr * constant
1032 37021751430 : self%sub = self%sub * constant
1033 37444469622 : self%diag = self%diag * constant
1034 :
1035 422718192 : self%left_bound = self%left_bound * constant
1036 422718192 : self%right_bound = self%right_bound * constant
1037 :
1038 25315992 : end subroutine scalar_lmult_tridiag
1039 :
1040 : ! Multiply in an array as if it contained the entries of a diagonal matrix
1041 : ! being multiplied from the left.
1042 0 : subroutine diagonal_lmult_tridiag(self, diag_array)
1043 :
1044 : class(TriDiagOp), intent(inout) :: self
1045 : real(r8), USE_CONTIGUOUS intent(in) :: diag_array(:,:)
1046 :
1047 0 : self%spr = self%spr * diag_array(:,:self%ncel-1)
1048 0 : self%sub = self%sub * diag_array(:,2:)
1049 0 : self%diag = self%diag * diag_array(:,:)
1050 :
1051 0 : self%left_bound = self%left_bound * diag_array(:,1)
1052 0 : self%right_bound = self%right_bound * diag_array(:,self%ncel)
1053 :
1054 0 : end subroutine diagonal_lmult_tridiag
1055 :
1056 : ! Decomposition constructor
1057 : !
1058 : ! The equation to be solved later (with left_div) is:
1059 : ! - A(k)*q(k+1) + B(k)*q(k) - C(k)*q(k-1) = D(k)
1060 : !
1061 : ! The solution (effectively via LU decomposition) has the form:
1062 : ! E(k) = C(k) / (B(k) - A(k)*E(k+1))
1063 : ! F(k) = (D(k) + A(k)*F(k+1)) / (B(k) - A(k)*E(k+1))
1064 : ! q(k) = E(k) * q(k-1) + F(k)
1065 : !
1066 : ! Unlike Richtmyer and Morton, E and F are defined by iterating backward
1067 : ! down to level 1, and then q iterates forward.
1068 : !
1069 : ! E can be calculated and stored now, without knowing D.
1070 : ! To calculate F later, we store A and the denominator.
1071 25315992 : function new_TriDiagDecomp(op, graft_decomp) result(decomp)
1072 : type(TriDiagOp), intent(in) :: op
1073 : type(TriDiagDecomp), intent(in), optional :: graft_decomp
1074 :
1075 : type(TriDiagDecomp) :: decomp
1076 :
1077 : integer :: k
1078 :
1079 25315992 : if (present(graft_decomp)) then
1080 0 : decomp%nsys = graft_decomp%nsys
1081 0 : decomp%ncel = graft_decomp%ncel
1082 : else
1083 25315992 : decomp%nsys = op%nsys
1084 25315992 : decomp%ncel = op%ncel
1085 : end if
1086 :
1087 : ! Simple allocation with no error checking.
1088 101263968 : allocate(decomp%ca(decomp%nsys,decomp%ncel))
1089 75947976 : allocate(decomp%dnom(decomp%nsys,decomp%ncel))
1090 75947976 : allocate(decomp%ze(decomp%nsys,decomp%ncel))
1091 :
1092 : ! decomp%ca is simply the negative of the tridiagonal's superdiagonal.
1093 37021751430 : decomp%ca(:,:op%ncel-1) = -op%spr
1094 422718192 : decomp%ca(:,op%ncel) = -op%right_bound
1095 :
1096 25315992 : if (present(graft_decomp)) then
1097 : ! Copy in graft_decomp beyond op%ncel.
1098 0 : decomp%ca(:,op%ncel+1:) = graft_decomp%ca(:,op%ncel+1:)
1099 0 : decomp%dnom(:,op%ncel+1:) = graft_decomp%dnom(:,op%ncel+1:)
1100 0 : decomp%ze(:,op%ncel+1:) = graft_decomp%ze(:,op%ncel+1:)
1101 : ! Fill in dnom edge value.
1102 0 : decomp%dnom(:,op%ncel) = 1._r8 / (op%diag(:,op%ncel) - &
1103 0 : decomp%ca(:,op%ncel)*decomp%ze(:,op%ncel+1))
1104 : else
1105 : ! If no grafting, the edge value of dnom comes from the diagonal.
1106 422718192 : decomp%dnom(:,op%ncel) = 1._r8 / op%diag(:,op%ncel)
1107 : end if
1108 :
1109 2240922219 : do k = op%ncel - 1, 1, -1
1110 36996435438 : decomp%ze(:,k+1) = - op%sub(:,k) * decomp%dnom(:,k+1)
1111 : decomp%dnom(:,k) = 1._r8 / &
1112 37021751430 : (op%diag(:,k) - decomp%ca(:,k)*decomp%ze(:,k+1))
1113 : end do
1114 :
1115 : ! Don't multiply edge level by denom, because we want to leave it up to
1116 : ! the BoundaryCond object to decide what this means in left_div.
1117 448034184 : decomp%ze(:,1) = -op%left_bound
1118 :
1119 25315992 : end function new_TriDiagDecomp
1120 :
1121 : ! Left-division (multiplication by inverse) using a decomposed operator.
1122 : !
1123 : ! See the comment above for the constructor for a quick explanation of the
1124 : ! intermediate variables. The "q" argument is "D(k)" on input and "q(k)" on
1125 : ! output.
1126 253159920 : subroutine decomp_left_div(decomp, q, l_cond, r_cond)
1127 :
1128 : ! Decomposed matrix.
1129 : class(TriDiagDecomp), intent(in) :: decomp
1130 : ! Data to left-divide by the matrix.
1131 : real(r8), USE_CONTIGUOUS intent(inout) :: q(:,:)
1132 : ! Objects representing boundary conditions.
1133 : class(BoundaryCond), intent(in), optional :: l_cond, r_cond
1134 :
1135 : ! "F" from the equation above.
1136 506319840 : real(r8) :: zf(decomp%nsys,decomp%ncel)
1137 :
1138 : ! Level index.
1139 : integer :: k
1140 :
1141 : ! Include boundary conditions.
1142 253159920 : if (present(l_cond)) then
1143 26354952 : q(:,1) = q(:,1) + l_cond%apply_left(decomp%ze(:,1), q)
1144 : end if
1145 :
1146 253159920 : if (present(r_cond)) then
1147 0 : q(:,decomp%ncel) = q(:,decomp%ncel) + &
1148 0 : r_cond%apply_right(decomp%ca(:,decomp%ncel), q)
1149 : end if
1150 :
1151 4227181920 : zf(:,decomp%ncel) = q(:,decomp%ncel) * decomp%dnom(:,decomp%ncel)
1152 :
1153 22182292116 : do k = decomp%ncel - 1, 1, -1
1154 >36643*10^7 : zf(:,k) = (q(:,k) + decomp%ca(:,k)*zf(:,k+1)) * decomp%dnom(:,k)
1155 : end do
1156 :
1157 : ! Perform back substitution
1158 :
1159 4480341840 : q(:,1) = zf(:,1)
1160 :
1161 22182292116 : do k = 2, decomp%ncel
1162 >36643*10^7 : q(:,k) = zf(:,k) + decomp%ze(:,k)*q(:,k-1)
1163 : end do
1164 :
1165 253159920 : end subroutine decomp_left_div
1166 :
1167 : ! Decomposition deallocation.
1168 25315992 : subroutine decomp_finalize(decomp)
1169 : class(TriDiagDecomp), intent(inout) :: decomp
1170 :
1171 25315992 : decomp%nsys = 0
1172 25315992 : decomp%ncel = 0
1173 :
1174 25315992 : if (allocated(decomp%ca)) deallocate(decomp%ca)
1175 25315992 : if (allocated(decomp%dnom)) deallocate(decomp%dnom)
1176 25315992 : if (allocated(decomp%ze)) deallocate(decomp%ze)
1177 :
1178 25315992 : end subroutine decomp_finalize
1179 :
1180 0 : end module linear_1d_operators
1181 :
|
