Downward sweep: solves, layer by layer from the top of every column,
(D - M^{-1} e_N e_N^T) q^(k) = r^(k) - M^{-1} e_N q_in^(k) , q_in^(nz) = qTop , q_in^(k) = q^(k+1)_0 ,
on every vertical line (i,j) and for every variable, as q^(k) = q_in^(k) 1 + At^{-1} r^(k). This is the upwind DG discretization of dq/dzeta = r with the inflow at the top of each element, the adjoint of VerticalSweepUp. r is in reference scaling, as for VerticalSweepUp. Requires Gauss-Lobatto control nodes.
Input - r : right hand side, volume field - qTop : inflow at the top of every column, surface field Output - q : solution, volume field, overwritten. It must be a different object from r.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(ColumnOperators3D_t), | intent(in) | :: | this | |||
| class(Scalar3D), | intent(in) | :: | r | |||
| class(Scalar2D), | intent(in) | :: | qTop | |||
| class(Scalar3D), | intent(inout) | :: | q |
subroutine VerticalSweepDown_ColumnOperators3D_t(this,r,qTop,q)
!! Downward sweep: solves, layer by layer from the top of every column,
!!
!! (D - M^{-1} e_N e_N^T) q^(k) = r^(k) - M^{-1} e_N q_in^(k) ,
!! q_in^(nz) = qTop , q_in^(k) = q^(k+1)_0 ,
!!
!! on every vertical line (i,j) and for every variable, as q^(k) = q_in^(k) 1 + At^{-1} r^(k).
!! This is the upwind DG discretization of dq/dzeta = r with the inflow at the top of
!! each element, the adjoint of VerticalSweepUp. r is in reference scaling, as for
!! VerticalSweepUp. Requires Gauss-Lobatto control nodes.
!!
!! Input
!! - r : right hand side, volume field
!! - qTop : inflow at the top of every column, surface field
!! Output
!! - q : solution, volume field, overwritten. It must be a different object from r.
implicit none
class(ColumnOperators3D_t),intent(in) :: this
class(Scalar3D),intent(in) :: r
class(Scalar2D),intent(in) :: qTop
class(Scalar3D),intent(inout) :: q
! Local
integer :: i,j,l,m,k,e2,iel,ivar,nz,np
call CheckSweep_ColumnOperators3D_t(this,r,qTop,q,"VerticalSweepDown")
nz = this%nLayers
np = this%N+1
do ivar = 1,r%nVar
do e2 = 1,this%nColumns
do k = nz,1,-1
iel = (e2-1)*nz+k
do l = 1,np
do j = 1,np
do i = 1,np
q%interior(i,j,l,iel,ivar) = 0.0_prec
enddo
enddo
do m = 1,np
do j = 1,np
do i = 1,np
q%interior(i,j,l,iel,ivar) = q%interior(i,j,l,iel,ivar)+ &
this%downInverse(m,l)*r%interior(i,j,m,iel,ivar)
enddo
enddo
enddo
enddo
! Inflow: the top trace for the top layer, else the bottom node of the layer above
if(k == nz) then
do l = 1,np
do j = 1,np
do i = 1,np
q%interior(i,j,l,iel,ivar) = qTop%interior(i,j,e2,ivar)+q%interior(i,j,l,iel,ivar)
enddo
enddo
enddo
else
do l = 1,np
do j = 1,np
do i = 1,np
q%interior(i,j,l,iel,ivar) = q%interior(i,j,1,iel+1,ivar)+q%interior(i,j,l,iel,ivar)
enddo
enddo
enddo
endif
enddo
enddo
enddo
endsubroutine VerticalSweepDown_ColumnOperators3D_t