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