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