box_gs_clpl Subroutine

private subroutine box_gs_clpl(mg, id, nc, redblack_cntr)

Perform Gauss-Seidel relaxation on box for a Laplacian operator in cylindrical geometry. TODO: in 3D this does not converge well, maybe it will for a stretched grid.

Arguments

Type IntentOptional Attributes Name
type(mg_t), intent(inout) :: mg
integer, intent(in) :: id
integer, intent(in) :: nc
integer, intent(in) :: redblack_cntr

Iteration counter


Source Code

  subroutine box_gs_clpl(mg, id, nc, redblack_cntr)
    type(mg_t), intent(inout) :: mg
    integer, intent(in)       :: id
    integer, intent(in)       :: nc
    integer, intent(in)       :: redblack_cntr
    !< Iteration counter
    integer                   :: i, j, i0, di
    logical                   :: redblack
    real(dp)                  :: idr2(2), dr(2), dr2(2), fac
    real(dp)                  :: r_face(nc+1), r_inv(nc)

    dr     = mg%dr(:, mg%boxes(id)%lvl)
    dr2    = dr**2
    idr2   = 1/dr**2
    fac    = 0.5_dp / sum(idr2)
    r_face = mg%boxes(id)%r_min(1) + dr(1) * [(i, i=0,nc)]
    r_inv  = 1/(mg%boxes(id)%r_min(1) + dr(1) * [(i-0.5_dp, i=1,nc)])

    i0  = 1
    redblack = (mg%smoother_type == mg_smoother_gsrb)
    if (redblack) then
       di = 2
    else
       di = 1
    end if

    ! The parity of redblack_cntr determines which cells we use. If
    ! redblack_cntr is even, we use the even cells and vice versa.
    associate (cc => mg%boxes(id)%cc, n => mg_iphi)
      do j = 1, nc
         if (redblack) i0 = 2 - iand(ieor(redblack_cntr, j), 1)

         do i = i0, nc, di
            cc(i, j, n) = fac * (idr2(1) * (r_face(i+1) * r_inv(i) * cc(i+1, j,&
                n) + r_face(i) * r_inv(i) * cc(i-1, j, n)) + idr2(2) * (cc(i,&
                j+1, n) + cc(i, j-1, n)) - cc(i, j, mg_irhs))
         end do
      end do
    end associate
  end subroutine box_gs_clpl