Assign level pointer
Loop over sweeps
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(pf_imexQ_oc_t), | intent(inout) | :: | this | |||
| type(pf_pfasst_t), | intent(inout), | target | :: | pf | PFASST structure |
|
| integer, | intent(in) | :: | level_index | which level this is |
||
| real(kind=pfdp), | intent(in) | :: | t0 | Time at beginning of time step |
||
| real(kind=pfdp), | intent(in) | :: | dt | time step size |
||
| integer, | intent(in) | :: | nsweeps | number of sweeps to do |
||
| integer, | intent(in), | optional | :: | flags |
subroutine imexQ_oc_sweep(this, pf, level_index, t0, dt, nsweeps, flags)
use pf_mod_timer
use pf_mod_hooks
class(pf_imexQ_oc_t), intent(inout) :: this
type(pf_pfasst_t), intent(inout),target :: pf !! PFASST structure
real(pfdp), intent(in ) :: t0 !! Time at beginning of time step
real(pfdp), intent(in ) :: dt !! time step size
integer, intent(in) :: level_index !! which level this is
integer, intent(in) :: nsweeps !! number of sweeps to do
integer, optional, intent(in ) :: flags
class(pf_level_t), pointer :: lev !! points to current level
! indicate if sweep on both (0, default; might skip y or p if tolerance satisfied), just y (1), just p (2)
integer :: k, m, n,Nnodes, which
real(pfdp) :: t,tend
logical :: sweep_y, sweep_p
real(pfdp), allocatable :: norms_y(:) !, norms_p(Lev%nnodes-1)
integer ::step
lev => pf%levels(level_index) !! Assign level pointer
step = pf%state%step+1
! print *, 'sweep on step', step
which = 0
if (present(flags)) which = flags
if (.not.present(flags)) stop "IMEXQ_OC SWEEPER WITHOUT FLAGS"
! print *, "IMEXQ_OC SWEEP", which
Nnodes = lev%nnodes
tend = t0+dt
call start_timer(pf, TLEVEL+lev%index-1)
if (which .eq. 1) then
sweep_y = .true.
sweep_p = .false.
else if (which .eq. 2) then
sweep_y = .false.
sweep_p = .true.
else
sweep_y = .true.
sweep_p = .true.
allocate(norms_y(lev%nnodes-1))
do m = 1, Nnodes-1
norms_y(m) = lev%R(m)%norm(1)
end do
if ( maxval(abs(norms_y)) < pf%abs_res_tol ) then
sweep_y = .false.
if (level_index == pf%nlevels) pf%state%skippedy = pf%state%skippedy + 1
end if
deallocate(norms_y)
!if ( maxval(abs(norms_p)) < pf%abs_res_tol ) sweep_p = .false.
end if
! if( sweep_p .and. pf%rank == 0) print *, "sweep on p with which = ", which
do k = 1,nsweeps !! Loop over sweeps
call call_hooks(pf, level_index, PF_PRE_SWEEP)
! compute integrals from previous iteration and add fas correction
! do m = 1, Nnodes-1
! call Lev%encap%setval(Lev%S(m), 0.0_pfdp, 1)
! call Lev%encap%setval(Lev%S(m), 0.0_pfdp, 2)
if( sweep_y ) then
do m = 1, Nnodes-1
call lev%I(m)%setval(0.0_pfdp, 1)
! Forward in y
if (this%explicit) then
do n = 1, Nnodes
call lev%I(m)%axpy(dt*this%QdiffE(m,n), lev%F(n,1),1)
! call Lev%encap%axpy(Lev%S(m), dt*imexQ_oc%QdiffE(m,n), Lev%F(n,1),1)
end do
end if
if (this%implicit) then
do n = 1, lev%nnodes
call lev%I(m)%axpy(dt*this%QdiffI(m,n), lev%F(n,2), 1)
! call Lev%encap%axpy(Lev%S(m), dt*imexQ_oc%QdiffI(m,n), Lev%F(n,2),1)
end do
end if
if (allocated(lev%tauQ)) then
call lev%I(m)%axpy(1.0_pfdp, lev%tauQ(m), 1)
! call Lev%encap%axpy(Lev%S(m), 1.0_pfdp, Lev%tauQ(m),1)
end if
end do
end if
if( sweep_p ) then
do m = Nnodes-1,1,-1
call lev%I(m)%setval(0.0_pfdp, 2)
! Backward in p, note S(m) goes backward now
!2 do n = 1,Nnodes
!2 call Lev%encap%axpy(Lev%S(Nnodes-m), dt*imexQ_oc%QdiffE(m,n), Lev%F(Nnodes+1-n,1),2)
!2 call Lev%encap%axpy(Lev%S(Nnodes-m), dt*imexQ_oc%QdiffI(m,n), Lev%F(Nnodes+1-n,2),2)
!2 end do
if (this%explicit) then
do n = Nnodes,1,-1
call lev%I(m)%axpy(dt*this%QdiffE(Nnodes-m,Nnodes+1-n), lev%F(n,1), 2)
!call Lev%encap%axpy(Lev%S(m), dt*imexQ_oc%QdiffE(Nnodes-m,Nnodes+1-n), Lev%F(n,1),2)
end do
end if
if (this%implicit) then
do n = Nnodes,1,-1
call lev%I(m)%axpy(dt*this%QdiffI(Nnodes-m,Nnodes+1-n), lev%F(n,2), 2)
! call Lev%encap%axpy(Lev%S(m), dt*imexQ_oc%QdiffI(Nnodes-m,Nnodes+1-n), Lev%F(n,2),2)
end do
end if
if (allocated(lev%tauQ)) then
call lev%I(m)%axpy(1.0_pfdp, lev%tauQ(m), 2)
end if
end do
end if
! Reload the newest initial values
! Recompute first function values
if (sweep_y) then
if (k .eq. 1) then
call lev%Q(1)%copy(lev%q0, 1)
if (this%explicit) &
call this%f_eval(lev%Q(1), t0, lev%index, lev%F(1,1), 1, 1, 1, step)
if (this%implicit) &
call this%f_eval(lev%Q(1), t0, lev%index, lev%F(1,2), 2, 1, 1, step)
end if
! call Lev%encap%copy(Lev%Q(1), Lev%q0, 1)
! call imexQ_oc%f1eval(Lev%Q(1), t0, Lev%level, Lev%ctx, Lev%F(1,1), 1, 1, step)
! call imexQ_oc%f2eval(Lev%Q(1), t0, Lev%level, Lev%ctx, Lev%F(1,2), 1)
end if
!else
! if( sweep_p ) then
! if (k .eq. 1) then
! call lev%Q(Nnodes)%copy(lev%qend, 2)
! if (this%explicit) &
! call this%f_eval(lev%Q(Nnodes), tend, lev%index, lev%F(Nnodes,1), 1, 2, Nnodes, step)
! if (this%implicit) &
! call this%f_eval(lev%Q(Nnodes), tend, lev%index, lev%F(Nnodes,2), 2, 2, Nnodes, step)
! end if
! ! call Lev%encap%copy(Lev%Q(Nnodes), Lev%qend, 2)
! ! call imexQ_oc%f1eval(Lev%Q(Nnodes), tend, Lev%level, Lev%ctx, Lev%F(Nnodes,1), 2, Nnodes, step)
! ! call imexQ_oc%f2eval(Lev%Q(Nnodes), tend, Lev%level, Lev%ctx, Lev%F(Nnodes,2), 2)
! end if
! if (sweep_p) then
! ! Backward sweep on p
! t = tend
! do m = Nnodes-1,1,-1
! t = t - dt*this%dtsdc(m)
!
! ! Do the dirk parts
! call this%rhs%setval(0.0_pfdp, 2)
! do n = Nnodes, m+1,-1
! if (this%explicit) &
! call this%rhs%axpy(dt*this%QtilE(Nnodes-m,Nnodes-n+1), lev%F(n,1), 2)
! if (this%implicit) &
! call this%rhs%axpy(dt*this%QtilI(Nnodes-m,Nnodes-n+1), lev%F(n,2), 2)
! end do
!
! call this%rhs%axpy(1.0_pfdp, lev%I(m), 2)
! call this%rhs%axpy(1.0_pfdp, lev%Q(Nnodes), 2)
!
! ! Do implicit solve
! if (this%implicit) then
! call this%f_comp(lev%Q(m), t, dt*this%QtilI(Nnodes-m,Nnodes-m+1), this%rhs, lev%index, lev%F(m,2), 2, 2)
! else
! call lev%Q(m)%copy(this%rhs,2)
! end if
! if (this%explicit) &
! call this%f_eval(lev%Q(m), t, lev%index, lev%F(m,1), 1, 2, m, step)
! end do
! ! reset first value
! call lev%q0%copy(lev%Q(1), 2)
! call pf_residual(pf, lev, dt, 2)
! ! call pf_residual(pf, lev, dt, which)
! end if
! Make some space
! call Lev%encap%create(rhs, Lev%level, SDC_KIND_SOL_FEVAL, Lev%nvars, Lev%shape, Lev%ctx)
if (sweep_y) then
! Forward sweep on y
t = t0
do m = 1, Nnodes-1
t = t + dt*this%dtsdc(m) ! forward running time
! Form rhs with all explicit terms
call this%rhs%setval(0.0_pfdp,1)
do n = 1, m
if (this%explicit) &
call this%rhs%axpy(dt*this%QtilE(m,n), lev%F(n,1), 1)
if (this%implicit) &
call this%rhs%axpy(dt*this%QtilI(m,n), lev%F(n,2), 1)
end do
call this%rhs%axpy(1.0_pfdp, lev%I(m), 1)
call this%rhs%axpy(1.0_pfdp, lev%Q(1), 1)
! Do implicit solve
if (this%implicit) then
call this%f_comp(lev%Q(m+1), t, dt*this%QtilI(m,m+1), this%rhs, lev%index, lev%F(m+1,2), 2, 1)
else
call lev%Q(m+1)%copy(this%rhs,1)
end if
! Compute explicit piece on new value
if (this%explicit) &
call this%f_eval(lev%Q(m+1), t, lev%index, lev%F(m+1,1), 1, 1, m+1, step)
end do
! Reset last values
call lev%qend%copy(lev%Q(Nnodes), 1)
! call pf_residual(pf, lev, dt, 1)
! call pf_residual(pf, lev, dt, which)
end if
if( sweep_p ) then
! do m=1, Nnodes-1
! call lev%I(m)%setval(0.0_pfdp, 2)
!
! ! Backward in p, note S(m) goes backward now
! !2 do n = 1,Nnodes
! !2 call Lev%encap%axpy(Lev%S(Nnodes-m), dt*imexQ_oc%QdiffE(m,n), Lev%F(Nnodes+1-n,1),2)
! !2 call Lev%encap%axpy(Lev%S(Nnodes-m), dt*imexQ_oc%QdiffI(m,n), Lev%F(Nnodes+1-n,2),2)
! !2 end do
! if (this%explicit) then
! do n = Nnodes,1,-1
! call lev%I(m)%axpy(dt*this%QdiffE(Nnodes-m,Nnodes+1-n), lev%F(n,1), 2)
! !call Lev%encap%axpy(Lev%S(m), dt*imexQ_oc%QdiffE(Nnodes-m,Nnodes+1-n), Lev%F(n,1),2)
! end do
! end if
! if (this%implicit) then
! do n = Nnodes,1,-1
! call lev%I(m)%axpy(dt*this%QdiffI(Nnodes-m,Nnodes+1-n), lev%F(n,2), 2)
! ! call Lev%encap%axpy(Lev%S(m), dt*imexQ_oc%QdiffI(Nnodes-m,Nnodes+1-n), Lev%F(n,2),2)
! end do
! end if
! if (allocated(lev%tauQ)) then
! call lev%I(m)%axpy(1.0_pfdp, lev%tauQ(m), 2)
! end if
! end do
if (k .eq. 1) then
call lev%Q(Nnodes)%copy(lev%qend, 2)
if (this%explicit) &
call this%f_eval(lev%Q(Nnodes), tend, lev%index, lev%F(Nnodes,1), 1, 2, Nnodes, step)
if (this%implicit) &
call this%f_eval(lev%Q(Nnodes), tend, lev%index, lev%F(Nnodes,2), 2, 2, Nnodes, step)
end if
! call Lev%encap%copy(Lev%Q(Nnodes), Lev%qend, 2)
! call imexQ_oc%f1eval(Lev%Q(Nnodes), tend, Lev%level, Lev%ctx, Lev%F(Nnodes,1), 2, Nnodes, step)
! call imexQ_oc%f2eval(Lev%Q(Nnodes), tend, Lev%level, Lev%ctx, Lev%F(Nnodes,2), 2)
end if
if (sweep_p) then
! Backward sweep on p
t = tend
do m = Nnodes-1,1,-1
t = t - dt*this%dtsdc(m)
! Do the dirk parts
call this%rhs%setval(0.0_pfdp, 2)
do n = Nnodes, m+1,-1
if (this%explicit) &
call this%rhs%axpy(dt*this%QtilE(Nnodes-m,Nnodes-n+1), lev%F(n,1), 2)
if (this%implicit) &
call this%rhs%axpy(dt*this%QtilI(Nnodes-m,Nnodes-n+1), lev%F(n,2), 2)
end do
call this%rhs%axpy(1.0_pfdp, lev%I(m), 2)
call this%rhs%axpy(1.0_pfdp, lev%Q(Nnodes), 2)
! Do implicit solve
if (this%implicit) then
call this%f_comp(lev%Q(m), t, dt*this%QtilI(Nnodes-m,Nnodes-m+1), this%rhs, lev%index, lev%F(m,2), 2, 2)
else
call lev%Q(m)%copy(this%rhs,2)
end if
if (this%explicit) &
call this%f_eval(lev%Q(m), t, lev%index, lev%F(m,1), 1, 2, m, step)
end do
! reset first value
call lev%q0%copy(lev%Q(1), 2)
! call pf_residual(pf, lev, dt, 2)
end if
! if(sweep_p) &
! call this%evaluate_all(lev, t0 + dt*lev%nodes, 2, step)
if( sweep_p .and. sweep_y ) then
call pf_residual(pf, lev, dt, 0)
else if( sweep_y ) then
call pf_residual(pf, lev, dt, 1)
else if (sweep_p ) then
call pf_residual(pf, lev, dt, 2)
else
stop "neither sweep on p nor on y : that should not happen"
end if
! done
call call_hooks(pf, level_index, PF_POST_SWEEP)
end do !nsweeps
call end_timer(pf, TLEVEL+lev%index-1)
end subroutine imexQ_oc_sweep