imexQ_oc_residual Subroutine

public subroutine imexQ_oc_residual(this, lev, dt, flags)

Arguments

Type IntentOptional AttributesName
class(pf_imexQ_oc_t), intent(inout) :: this
class(pf_level_t), intent(inout) :: lev
real(kind=pfdp), intent(in) :: dt
integer, intent(in), optional :: flags

Contents

Source Code


Source Code

  subroutine imexQ_oc_residual(this, lev, dt, flags)
    class(pf_imexQ_oc_t),  intent(inout) :: this
    class(pf_level_t), intent(inout)  :: lev
    real(pfdp),        intent(in)     :: dt
    integer, intent(in), optional     :: flags
    
    integer :: m, which

    which = 0
    if(present(flags)) then
      which = flags
    else
      print *, "flags not present in residual", which
      stop
    end if
!     print *, "IMEXQ_OC RESIDUAL ", which

    
    call this%integrate(lev, lev%Q, lev%F, dt, lev%I, which)

    ! add tau (which is 'node to node')
    if (allocated(lev%tauQ)) then
       do m = 1, lev%nnodes-1
          call lev%I(m)%axpy(1.0_pfdp, lev%tauQ(m), which)
       end do
    end if

    ! subtract out Q
    do m = 1, lev%nnodes-1
       if( (which .eq. 0) .or. (which .eq. 1) ) then
          call lev%R(m)%copy(lev%I(m), 1)
          call lev%R(m)%axpy(1.0_pfdp, lev%Q(1), 1)
          call lev%R(m)%axpy(-1.0_pfdp, lev%Q(m+1), 1)
       end if
       if( (which .eq. 0) .or. (which .eq. 2) ) then
          call lev%R(m)%copy(lev%I(m), 2)
          call lev%R(m)%axpy(1.0_pfdp,  lev%Q(lev%nnodes), 2)
          call lev%R(m)%axpy(-1.0_pfdp, lev%Q(m), 2)
       end if
    end do

  end subroutine imexQ_oc_residual