LCOV - code coverage report
Current view: | top level - physics/cam - vdiff_lu_solver.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 15 | 39 | 38.5 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 1 | 2 | 50.0 %

              Line data    Source code
           1             : module vdiff_lu_solver
           2             :
           3             : ! This module provides a function returning the matrix decomposition for
           4             : ! an implicit finite volume solver for vertical diffusion. It accepts
           5             : ! diffusion coefficients, time/grid spacing, and boundary condition
           6             : ! objects, and returns a TriDiagDecomp object that can be used to diffuse
           7             : ! an array for one time step with the "left_div" method.
           8             :
           9             : use coords_1d, only: Coords1D
          10             : use linear_1d_operators, only: TriDiagOp, operator(+), TriDiagDecomp
          11             :
          12             : implicit none
          13             : private
          14             : save
          15             :
          16             : ! Public interfaces
          17             : public :: vd_lu_decomp
          18             : public :: fin_vol_lu_decomp
          19             :
          20             : ! 8-byte real.
          21             : integer, parameter :: r8 = selected_real_kind(12)
          22             :
          23             : contains
          24             :
          25             : ! ========================================================================!
          26             :
          27             : ! Designed to solve the equation:
          28             : ! dq/dt = c1 q'' + c2 q' + c q
          29             :
          30           0 : function vd_lu_decomp(dt, dp, coef_q,  coef_q_d, coef_q_d2, upper_bndry, &
          31             :      lower_bndry) result(decomp)
          32             :
          33             :   use linear_1d_operators, only: &
          34             :        identity_operator, &
          35             :        diagonal_operator, &
          36             :        first_derivative, &
          37             :        second_derivative, &
          38             :        BoundaryType
          39             :
          40             :   ! ---------------------- !
          41             :   ! Input-Output Arguments !
          42             :   ! ---------------------- !
          43             :
          44             :   ! Time step.
          45             :   real(r8), intent(in) :: dt
          46             :   ! Grid spacing (deltas).
          47             :   real(r8), USE_CONTIGUOUS intent(in) :: dp(:,:)
          48             :
          49             :   ! Coefficients for q, q', and q''.
          50             :   real(r8), USE_CONTIGUOUS intent(in), optional :: coef_q(:,:), &
          51             :        coef_q_d(:,:), coef_q_d2(:,:)
          52             :
          53             :   ! Boundary conditions (optional, default to 0 flux through boundary).
          54             :   class(BoundaryType), target, intent(in), optional :: &
          55             :        upper_bndry, lower_bndry
          56             :
          57             :   ! Output decomposition.
          58             :   type(TriDiagDecomp) :: decomp
          59             :
          60             :   ! --------------- !
          61             :   ! Local Variables !
          62             :   ! --------------- !
          63             :
          64             :   ! Operator objects.
          65           0 :   type(TriDiagOp) :: add_term
          66           0 :   type(TriDiagOp) :: net_operator
          67             :
          68             :   ! ----------------------- !
          69             :   ! Main Computation Begins !
          70             :   ! ----------------------- !
          71             :
          72           0 :   if (present(coef_q)) then
          73           0 :      net_operator = diagonal_operator(1._r8 - dt*coef_q)
          74             :   else
          75           0 :      net_operator = identity_operator(size(dp, 1), size(dp, 2) + 1)
          76             :   end if
          77             :
          78           0 :   if (present(coef_q_d)) then
          79           0 :      add_term = first_derivative(dp, upper_bndry, lower_bndry)
          80           0 :      call add_term%lmult_as_diag(-dt*coef_q_d)
          81           0 :      call net_operator%add(add_term)
          82             :   end if
          83             :
          84           0 :   if (present(coef_q_d2)) then
          85           0 :      add_term = second_derivative(dp, upper_bndry, lower_bndry)
          86           0 :      call add_term%lmult_as_diag(-dt*coef_q_d2)
          87           0 :      call net_operator%add(add_term)
          88             :   end if
          89             :
          90           0 :   decomp = TriDiagDecomp(net_operator)
          91             :
          92           0 :   call net_operator%finalize()
          93           0 :   call add_term%finalize()
          94             :
          95           0 : end function vd_lu_decomp
          96             :
          97             : ! ========================================================================!
          98             :
          99             : ! Designed to solve the equation:
         100             : !
         101             : ! w * dq/dt = d/dp (D q' - v q) + c q
         102             : !
         103             : ! where q is a grid-cell average, and p is the vertical coordinate
         104             : ! (presumably pressure).
         105             : !
         106             : ! In this function, coef_q_weight == w, coef_q_diff == D,
         107             : ! coef_q_adv == v, and coef_q == c. All these are optional; omitting a
         108             : ! coefficient is equivalent to setting the entire array to 0.
         109             : !
         110             : ! coef_q_diff and coef_q_adv are defined at the level interfaces, while
         111             : ! coef_q and coef_q_weight are grid-cell averages.
         112             :
         113     2978352 : function fin_vol_lu_decomp(dt, p, coef_q, coef_q_diff, coef_q_adv, &
         114     2978352 :      coef_q_weight, upper_bndry, lower_bndry, graft_decomp) result(decomp)
         115             :
         116             :   use linear_1d_operators, only: &
         117             :        zero_operator, &
         118             :        diagonal_operator, &
         119             :        diffusion_operator, &
         120             :        advection_operator, &
         121             :        BoundaryType
         122             :
         123             :   ! ---------------------- !
         124             :   ! Input-Output Arguments !
         125             :   ! ---------------------- !
         126             :
         127             :   ! Time step.
         128             :   real(r8), intent(in) :: dt
         129             :   ! Grid spacings.
         130             :   type(Coords1D), intent(in) :: p
         131             :
         132             :   ! Coefficients for diffusion and advection.
         133             :   !
         134             :   ! The sizes must be consistent among all the coefficients that are
         135             :   ! actually present, i.e. coef_q_diff and coef_q_adv should be one level
         136             :   ! bigger than coef_q and coef_q_weight, and have the same column number.
         137             :   real(r8), USE_CONTIGUOUS intent(in), optional :: coef_q(:,:), &
         138             :        coef_q_diff(:,:), coef_q_adv(:,:), coef_q_weight(:,:)
         139             :
         140             :   ! Boundary conditions (optional, default to 0 flux through boundary).
         141             :   class(BoundaryType), target, intent(in), optional :: &
         142             :        upper_bndry, lower_bndry
         143             :
         144             :   ! Decomposition to graft onto. If this is provided, you can pass in
         145             :   ! smaller coefficients.
         146             :   type(TriDiagDecomp), intent(in), optional :: graft_decomp
         147             :
         148             :   ! Output decomposition.
         149             :   type(TriDiagDecomp) :: decomp
         150             :
         151             :   ! --------------- !
         152             :   ! Local Variables !
         153             :   ! --------------- !
         154             :
         155             :   ! Operator objects.
         156     2978352 :   type(TriDiagOp) :: add_term
         157     2978352 :   type(TriDiagOp) :: net_operator
         158             :
         159             :   ! ----------------------- !
         160             :   ! Main Computation Begins !
         161             :   ! ----------------------- !
         162             :
         163             :   ! A diffusion term is probably present, so start with that. Otherwise
         164             :   ! start with an operator of all 0s.
         165             :
         166     2978352 :   if (present(coef_q_diff)) then
         167             :      net_operator = diffusion_operator(p, coef_q_diff, &
         168     2978352 :           upper_bndry, lower_bndry)
         169             :   else
         170           0 :      net_operator = zero_operator(p%n, p%d)
         171             :   end if
         172             :
         173             :   ! Constant term (damping).
         174     2978352 :   if (present(coef_q)) then
         175           0 :      add_term = diagonal_operator(coef_q)
         176           0 :      call net_operator%add(add_term)
         177             :   end if
         178             :
         179             :   ! Effective advection.
         180     2978352 :   if (present(coef_q_adv)) then
         181             :      add_term = advection_operator(p, coef_q_adv, &
         182           0 :           upper_bndry, lower_bndry)
         183           0 :      call net_operator%add(add_term)
         184             :   end if
         185             :
         186             :   ! We want I-dt*(w^-1)*A for a single time step, implicit method, where
         187             :   ! A is the right-hand-side operator (i.e. what net_operator is now).
         188     2978352 :   if (present(coef_q_weight)) then
         189           0 :      call net_operator%lmult_as_diag(-dt/coef_q_weight)
         190             :   else
         191     2978352 :      call net_operator%lmult_as_diag(-dt)
         192             :   end if
         193     2978352 :   call net_operator%add_to_diag(1._r8)
         194             :
         195             :   ! Decompose, grafting on an optional input decomp. The graft is a way to
         196             :   ! avoid re-calculating the ending (bottom) levels when the coefficients
         197             :   ! have only changed at the beginning (top), e.g. for different
         198             :   ! constituents in the molecular diffusion.
         199     2978352 :   decomp = TriDiagDecomp(net_operator, graft_decomp=graft_decomp)
         200             :
         201             :   ! Ensure local objects are deallocated.
         202     2978352 :   call net_operator%finalize()
         203     2978352 :   call add_term%finalize()
         204             :
         205     5956704 : end function fin_vol_lu_decomp
         206             :
         207             : end module vdiff_lu_solver

Generated by: LCOV version 1.14