box_gs_ahelmh Subroutine

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

Perform Gauss-Seidel relaxation on box for a Helmholtz operator

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_ahelmh(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, i0, di
    logical                   :: redblack
    real(dp)                  :: idr2(2*1), u(2*1)
    real(dp)                  :: a0(2*1), a(2*1), c(2*1)

    ! Duplicate 1/dr^2 array to multiply neighbor values
    idr2(1:2*1:2) = 1/mg%dr(:, mg%boxes(id)%lvl)**2
    idr2(2:2*1:2) = idr2(1:2*1:2)
    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, i_eps1 => mg_iveps1)
      if (redblack) i0 = 2 - iand(redblack_cntr, 1)

      do i = i0, nc, di
         a0(1:2) = cc(i, i_eps1)
         u(1:2) = cc(i-1:i+1:2, n)
         a(1:2) = cc(i-1:i+1:2, i_eps1)
         c(:)   = 2 * a0(:) * a(:) / (a0(:) + a(:)) * idr2

         cc(i, n) = (sum(c(:) * u(:)) - cc(i,&
             mg_irhs)) / (sum(c(:)) + ahelmholtz_lambda)
      end do
    end associate
  end subroutine box_gs_ahelmh