REMORA_TimeStep.cpp

Go to the documentation of this file.

#include <REMORA.H>
 
using namespace amrex;
 
/**
 * @param[in] lev            level of refinement
 * @param[in] time           simulation time at start of step
 * @param[in] iteration      iteration in subcycling, if using
 */
void
REMORA::timeStep (int lev, Real time, int iteration)
{
    if (regrid_int > 0)  // We may need to regrid
    {
        // help keep track of whether a level was already regridded
        // from a coarser level call to regrid
        static Vector<int> last_regrid_step(max_level+1, 0);
 
        // regrid changes level "lev+1" so we don't regrid on max_level
        // also make sure we don't regrid fine levels again if
        // it was taken care of during a coarser regrid
        if (lev < max_level && istep[lev] > last_regrid_step[lev])
        {
            if (istep[lev] % regrid_int == 0)
            {
                // regrid interpolates this level's data to build lev+1, so its ghost cells
                // must be current. The swap has not happened, so "new" is what regrid reads.
                FillPatchNoBC(lev, time, *cons_new[lev], cons_new, BdyVars::t,0,true,true);
                FillPatchNoBC(lev, time, *xvel_new[lev], xvel_new, BdyVars::u,0,true,true);
                FillPatchNoBC(lev, time, *yvel_new[lev], yvel_new, BdyVars::v,0,true,true);
                FillPatch(lev, time, *zvel_new[lev], zvel_new, zvel_bc(), BdyVars::null,0,true,true);
 
                // regrid could add newly refine levels (if finest_level < max_level)
                // so we save the previous finest level index
                int old_finest = finest_level;
                regrid(lev, time);
 
                // The refined footprint moved, so the level pairs are new
                check_mask_consistency();
 
#ifdef REMORA_USE_PARTICLES
                if (finest_level != old_finest) {
                    particleData.Redistribute();
                }
#endif
 
                // Mark that we have regridded this level already
                for (int k = lev; k <= finest_level; ++k) {
                    last_regrid_step[k] = istep[k];
                }
 
                // If there are newly created levels, set the time step. nsubsteps, not the
                // refinement ratio: without subcycling every level shares dt[0].
                for (int k = old_finest+1; k <= finest_level; ++k) {
                    dt[k] = dt[k-1] / nsubsteps[k];
                }
            }
        }
    }
 
    scale_rhs_vars(lev);
 
    // Update what we call "old" and "new" time
    t_old[lev] = t_new[lev];
    // A child's last substep lands exactly on its parent's new time. Summing dt[lev] would
    // drift when it is inexact (dt/3 in single precision), leaving the next step's t_old
    // just before the parent's t_old -- outside the interval the fill patchers bracket.
    if (lev > 0 && iteration == nsubsteps[lev]) {
        t_new[lev] = t_new[lev-1];
    } else {
        t_new[lev] += dt[lev];
    }
 
    // A child must stay inside its parent's step, or the fill patchers have nothing to
    // interpolate within. REMORAFillPatcher asserts this too, but cannot name the level.
    if (lev > 0) {
        const Real eps = Real(1.e-6) * dt[lev];
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
            t_old[lev] >= t_old[lev-1] - eps && t_new[lev] <= t_new[lev-1] + eps,
            "REMORA::timeStep: level " + std::to_string(lev) + " stepped outside the time "
            "interval of level " + std::to_string(lev-1));
    }
 
    if (Verbose()) {
        amrex::Print() << "[Level " << lev << " step " << istep[lev]+1 << "] ";
        amrex::Print() << "ADVANCE from time = " << t_old[lev] << " to " << t_new[lev]
                       << " with dt = " << dt[lev] << std::endl;
    }
 
    // We must swap the pointers so the previous step's "new" is now this step's "old"
    std::swap(cons_old[lev], cons_new[lev]);
    std::swap(xvel_old[lev], xvel_new[lev]);
    std::swap(yvel_old[lev], yvel_new[lev]);
    std::swap(zvel_old[lev], zvel_new[lev]);
 
    // Advance a single level for a single time step
    Advance(lev, time, dt[lev], iteration, nsubsteps[lev]);
 
    ++istep[lev];
 
    scale_rhs_vars_inv(lev);
 
    if (Verbose())
    {
        amrex::Print() << "[Level " << lev << " step " << istep[lev] << "] ";
        amrex::Print() << "Advanced " << CountCells(lev) << " cells" << std::endl;
    }
 
    if (lev < finest_level)
    {
        // recursive call for next-finer level
        for (int i = 1; i <= nsubsteps[lev+1]; ++i)
        {
            timeStep(lev+1, time+(i-1)*dt[lev+1], i);
        }
 
        // Before the average-down: refluxing writes coarse cells under the fine grid on the
        // assumption they are about to be overwritten from it.
        reflux_to(lev);
 
        if (solverChoice.coupling_type == CouplingType::two_way) {
            AverageDownTo(lev); // average lev+1 down to lev
        }
    }
}