Step 1. Getting the initial condition on the finest level at each processor If we are doing multiple levels, then we need to coarsen to fine level Step 2: Proceed fine to coarse levels coarsening the fine solution and computing tau correction If RK_pred is true, just do some RK_steps
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(pf_pfasst_t), | intent(inout), | target | :: | pf | PFASST main data structure |
|
| real(kind=pfdp), | intent(in) | :: | t0 | Initial time of this processor |
||
| real(kind=pfdp), | intent(in) | :: | dt | time step |
||
| integer, | intent(in), | optional | :: | flags(:) | User defined flags |
subroutine pf_predictor_oc(pf, t0, dt, flags)
type(pf_pfasst_t), intent(inout), target :: pf !! PFASST main data structure
real(pfdp), intent(in ) :: t0 !! Initial time of this processor
real(pfdp), intent(in ) :: dt !! time step
integer, intent(in ), optional :: flags(:) !! User defined flags
class(pf_level_t), pointer :: c_lev_p
class(pf_level_t), pointer :: f_lev_p
integer :: k !! Loop indices
integer :: level_index !! Local variable for looping over levels
real(pfdp) :: t0k !! Initial time at time step k
integer :: which, dir, send_tag
which = 1 ! standard: predict and sweep forward-in-time
dir = 1 ! for MPI communication, standard is forward-in-time
if(present(flags)) then
if(flags(1)==2) then
which = 2 ! if we are computing an adjoint, predict and sweep backward-in-time
dir = 2 ! communication has to be backwards as well
end if
if(flags(1)==0) which = 0 ! sweep forward and backward simultaneously on two components, communication only forwards
end if
call call_hooks(pf, 1, PF_PRE_PREDICTOR)
call start_timer(pf, TPREDICTOR)
if (pf%debug) print*, 'DEBUG --', pf%rank, 'beginning predictor'
!! Step 1. Getting the initial condition on the finest level at each processor
!! If we are doing multiple levels, then we need to coarsen to fine level
f_lev_p => pf%levels(pf%nlevels)
if (pf%q0_style < 2) then ! Spread q0 to all the nodes
if( (which == 0) .or. (which == 1)) call f_lev_p%ulevel%sweeper%spreadq0(f_lev_p, t0, 1, pf%state%step+1)
if( (which == 0) .or. (which == 2)) call f_lev_p%ulevel%sweeper%spreadq0(f_lev_p, t0+dt, 2, pf%state%step+1)
endif
!! Step 2: Proceed fine to coarse levels coarsening the fine solution and computing tau correction
if (pf%debug) print*, 'DEBUG --', pf%rank, 'do coarsen in predictor'
if (pf%nlevels > 1) then
do level_index = pf%nlevels, 2, -1
f_lev_p => pf%levels(level_index);
c_lev_p => pf%levels(level_index-1)
call pf_residual(pf, f_lev_p, dt, which)
if( (which == 0) .or. (which == 1)) call f_lev_p%ulevel%restrict(f_lev_p, c_lev_p, f_lev_p%q0, c_lev_p%q0, t0, 1)
if( (which == 0) .or. (which == 2)) call f_lev_p%ulevel%restrict(f_lev_p, c_lev_p, f_lev_p%qend, c_lev_p%qend, t0+dt, 2)
call restrict_time_space_fas(pf, t0, dt, level_index, which) ! Restrict
call save(c_lev_p, which)
end do ! level_index = pf%nlevels, 2, -1
end if
level_index = 1
c_lev_p => pf%levels(1)
! Step 3. Do the "Burn in" step on the coarse level to make the coarse values consistent
! (this is skipped if the fine initial conditions are already consistent)
! The first processor does nothing, the second does one set of sweeps, the 2nd two, etc
! Hence, this is skipped completely if nprocs=1
if (pf%q0_style .eq. 0) then ! The coarse level needs burn in
if (pf%debug) print*, 'DEBUG --', pf%rank, 'do burnin in pred', ' RK_pred', pf%RK_pred, ' PFASST_pred', pf%PFASST_pred
!! If RK_pred is true, just do some RK_steps
if (pf%RK_pred) then ! Use Runge-Kutta to get the coarse initial data
! Get new initial conditions
call pf_recv(pf, c_lev_p, 100000+pf%rank, .true., dir)
! Do a RK_step
call c_lev_p%ulevel%sweeper%sweep(pf, level_index, t0, dt, 1, which )
! Send forward/backward
if (dir == 1) send_tag = 100000+pf%rank+1
if (dir == 2) send_tag = 100000+pf%rank-1
call pf_send(pf, c_lev_p, send_tag, .false., dir)
else ! Normal PFASST burn in
do k = 1, pf%rank + 1
pf%state%iter = -k
t0k = t0-(pf%rank)*dt + (k-1)*dt ! Remember t0=pf%rank*dt is the beginning of this time slice so t0-(pf%rank)*dt is 0
! and we iterate up to the correct time step.
! for optimal control problem t, t0k has no influence on f_eval, so there this does something else
! Get new initial value (skip on first iteration)
if (k > 1) then
if ((which == 0) .or. (which == 1)) call c_lev_p%q0%copy(c_lev_p%qend, 1)
! if ((which == 0) .or. (which == 2)) call c_lev_p%qend%copy(c_lev_p%q0, 2) ! for which==0, we solve with zero terminal conditions,
! but q0,2 is not zero (source term due to state sweeps)
if (which == 2) call c_lev_p%qend%copy(c_lev_p%q0, 2)
! If we are doing PFASST_pred, we use the old values at nodes, otherwise spread q0
if (.not. pf%PFASST_pred) then
if( (which == 0) .or. (which == 1)) call c_lev_p%ulevel%sweeper%spreadq0(c_lev_p, t0k, 1, pf%state%step+1)
! if( (which == 0) .or. (which == 2)) call c_lev_p%ulevel%sweeper%spreadq0(c_lev_p, t0k+dt, 2, pf%state%step+1)
if( which == 2) call c_lev_p%ulevel%sweeper%spreadq0(c_lev_p, t0k+dt, 2, pf%state%step+1)
end if
end if
! Do some sweeps
if( which == 0 .or. which == 1 ) call c_lev_p%ulevel%sweeper%sweep(pf, level_index, t0k, dt, pf%nsweeps_burn, 1)
if( which == 2 ) call c_lev_p%ulevel%sweeper%sweep(pf, level_index, t0k, dt, pf%nsweeps_burn, 2)
end do
endif ! RK_pred
end if ! (q0_style .eq. 0)
! Step 4: Now we have everyone burned in, so do some coarse sweeps
if (pf%nlevels > 1) then
pf%state%pstatus = PF_STATUS_ITERATING
pf%state%status = PF_STATUS_ITERATING
if (pf%debug) print*, 'DEBUG --', pf%rank, 'do sweeps in predictor', ' Pipeline_pred', pf%Pipeline_pred
level_index=1
c_lev_p => pf%levels(level_index)
if (pf%Pipeline_pred) then
do k = 1, c_lev_p%nsweeps_pred
pf%state%iter =-(pf%rank + 1) -k
! Get new initial conditions
call pf_recv(pf, c_lev_p, c_lev_p%index*110000+pf%rank+k, .true., dir)
! Do a sweep
call c_lev_p%ulevel%sweeper%sweep(pf, c_lev_p%index, t0, dt, 1, which)
! Send forward/backward
if (dir == 1) send_tag = c_lev_p%index*1110000+pf%rank+1+k
if (dir == 2) send_tag = c_lev_p%index*1110000+pf%rank-1+k
call pf_send(pf, c_lev_p, send_tag, .false., dir)
end do ! k = 1, c_lev_p%nsweeps_pred-1
else ! Don't pipeline
! Get new initial conditions
call pf_recv(pf, c_lev_p, c_lev_p%index*100000+pf%rank, .true., dir)
! Do sweeps
if(which == 0 .or. which == 1) call c_lev_p%ulevel%sweeper%sweep(pf, level_index, t0, dt, c_lev_p%nsweeps_pred, 1) !which ! why only state?
if(which == 2) call c_lev_p%ulevel%sweeper%sweep(pf, level_index, t0, dt, c_lev_p%nsweeps_pred, 2) !which
! Send forward/backward
if (dir == 1) send_tag = c_lev_p%index*100000+pf%rank+1
if (dir == 2) send_tag = c_lev_p%index*100000+pf%rank-1
call pf_send(pf, c_lev_p, send_tag, .false., dir)
endif ! (Pipeline_pred .eq. .true) then
end if ! pf%nlevels > 1
! Step 5: Return to fine level sweeping on any level in between coarsest and finest
if (pf%debug) print*, 'DEBUG --', pf%rank, 'returning to fine level in predictor'
do level_index = 2, pf%nlevels ! Will do nothing with one level
f_lev_p => pf%levels(level_index);
c_lev_p => pf%levels(level_index-1)
call interpolate_time_space(pf, t0, dt, level_index, c_lev_p%Finterp, which)
if ((which == 0) .or. (which == 1)) call interpolate_q0(pf, f_lev_p, c_lev_p, 1)
if (which == 2) call interpolate_qend(pf, f_lev_p, c_lev_p) ! for which==0, qend never changes, so don't need to interpolate
! Do sweeps on level unless we are at the finest level
if (level_index < pf%nlevels) then
call f_lev_p%ulevel%sweeper%sweep(pf, level_index, t0, dt, f_lev_p%nsweeps_pred, which)
end if
end do
call end_timer(pf, TPREDICTOR)
call call_hooks(pf, -1, PF_POST_PREDICTOR)
pf%state%iter = 0
pf%state%status = PF_STATUS_ITERATING
pf%state%pstatus = PF_STATUS_ITERATING
if (pf%debug) print*, 'DEBUG --', pf%rank, 'ending predictor'
end subroutine pf_predictor_oc