LevelRK4.H

Go to the documentation of this file.

#ifdef CH_LANG_CC
/*
*      _______              __
*     / ___/ /  ___  __ _  / /  ___
*    / /__/ _ \/ _ \/  V \/ _ \/ _ \
*    \___/_//_/\___/_/_/_/_.__/\___/
*    Please refer to Copyright.txt, in Chombo's root directory.
*/
#endif

#ifndef _LEVELRK4_H__
#define _LEVELRK4_H__



template <class TSoln, class TFR, class TOp>
void RK4LevelAdvance(TSoln& a_newSoln,
                     TSoln& a_oldSoln,
                     const TSoln& a_oldCrseSoln,
                     Real a_oldCrseTime,
                     const TSoln& a_newCrseSoln,
                     Real a_newCrseTime,
                     TFR& a_crseFRPtr,
                     TFR& a_fineFRPtr,
                     Real a_time,
                     Real a_dt,
                     TOp& a_op)
{
  // Advance by one time step using 4th-order Runge-Kutta.

  // allocate and define temp storage
  TSoln RHSTmp;
  TSoln solnTmp;
  a_op.defineRHSData(RHSTmp,a_oldSoln);
  a_op.defineSolnData(solnTmp,a_oldSoln);
  a_op.copySolnData(solnTmp,a_oldSoln);

  // start by copying old solution into new solution
  a_op.copySolnData(a_newSoln, a_oldSoln);

  // RK Stage 1: compute RHS.at old time
  // RHSTmp gets RHS(tOld)
  // flux registers are also incremented if appropriate
  a_op.evalRHS(RHSTmp,a_oldSoln,
               a_fineFRPtr,a_crseFRPtr,
               a_oldCrseSoln, a_oldCrseTime,
               a_newCrseSoln, a_newCrseTime,
               a_time,a_dt/6);

  // RK Stage 1: compute update {phi}_1, partial update to new vals.

  // first part of update...
  a_op.updateODE(a_newSoln,RHSTmp,a_dt/6);


  // predict half-time value using LofPhiTmp
  a_op.updateODE(solnTmp, RHSTmp,a_dt/2);

  // RK Stage 2: compute RHS based on first predicted value
  // also increment flux registers
  a_op.evalRHS(RHSTmp,solnTmp,
               a_fineFRPtr,a_crseFRPtr,
               a_oldCrseSoln, a_oldCrseTime,
               a_newCrseSoln, a_newCrseTime,
               a_time+a_dt/2,a_dt/3);


  // RK Stage 2: compute update {phi}_2, partial update to new vals.
  a_op.updateODE(a_newSoln,RHSTmp,a_dt/3);

  // now predict half-time value using latest estimate of RHS
  a_op.copySolnData(solnTmp,a_oldSoln);
  a_op.updateODE(solnTmp,RHSTmp,a_dt/2);

  // RK Stage 3: compute new RHS.

  a_op.evalRHS(RHSTmp,solnTmp,
               a_fineFRPtr,a_crseFRPtr,
               a_oldCrseSoln, a_oldCrseTime,
               a_newCrseSoln, a_newCrseTime,
               a_time+a_dt/2,a_dt/3);


  // RK Stage 3: compute update {phi}_3, partial update to new vals.

  a_op.updateODE(a_newSoln,RHSTmp,a_dt/3);

  // predict value at new time
  a_op.copySolnData(solnTmp,a_oldSoln);
  a_op.updateODE(solnTmp,RHSTmp,a_dt);

  // RK Stage 4: compute RHS.

  a_op.evalRHS(RHSTmp,solnTmp,
               a_fineFRPtr,a_crseFRPtr,
               a_oldCrseSoln, a_oldCrseTime,
               a_newCrseSoln, a_newCrseTime,
               a_time+a_dt,a_dt/6);

  // RK Stage 4: compute final update of solution.

  a_op.updateODE(a_newSoln,RHSTmp,a_dt/6);
}




template <class TSoln, class TInterp, class TFR, class TOp>
void RK4LevelAdvance(TSoln& a_newSoln,
                     TSoln& a_oldSoln,
                     TInterp& a_interpPtr,
                     const TSoln& a_oldCrseSoln,
                     Real a_oldCrseTime,
                     const TSoln& a_newCrseSoln,
                     Real a_newCrseTime,
                     TFR& a_crseFRPtr,
                     TFR& a_fineFRPtr,
                     Real a_time,
                     Real a_dt,
                     TOp& a_op)
{

  // Advance by one time step using 4th-order Runge-Kutta.

  // allocate and define temp storage
  TSoln RHSTmp;
  TSoln solnTmp;
  a_op.defineRHSData(RHSTmp,a_oldSoln);
  a_op.defineSolnData(solnTmp,a_oldSoln);
  a_op.copySolnData(solnTmp,a_oldSoln);

  // start by copying old solution into new solution
  a_op.copySolnData(a_newSoln, a_oldSoln);
  a_interpPtr.setDt(a_dt);
  a_interpPtr.saveInitialSoln(a_oldSoln);

  // RK Stage 1: compute RHS.at old time
  // RHSTmp gets RHS(tOld)
  // flux registers are also incremented if appropriate
  a_op.evalRHS(RHSTmp,a_oldSoln,
               a_fineFRPtr,a_crseFRPtr,
               a_oldCrseSoln, a_oldCrseTime,
               a_newCrseSoln, a_newCrseTime,
               a_time,a_dt/6);
  a_interpPtr.saveRHS(RHSTmp);

  // RK Stage 1: compute update {phi}_1, partial update to new vals.

  // first part of update...
  a_op.updateODE(a_newSoln,RHSTmp,a_dt/6);


  // predict half-time value using LofPhiTmp
  a_op.updateODE(solnTmp, RHSTmp,a_dt/2);

  // RK Stage 2: compute RHS based on first predicted value
  // also increment flux registers
  a_op.evalRHS(RHSTmp,solnTmp,
               a_fineFRPtr,a_crseFRPtr,
               a_oldCrseSoln, a_oldCrseTime,
               a_newCrseSoln, a_newCrseTime,
               a_time+a_dt/2,a_dt/3);
  a_interpPtr.saveRHS(RHSTmp);


  // RK Stage 2: compute update {phi}_2, partial update to new vals.
  a_op.updateODE(a_newSoln,RHSTmp,a_dt/3);

  // now predict half-time value using latest estimate of RHS
  a_op.copySolnData(solnTmp,a_oldSoln);
  a_op.updateODE(solnTmp,RHSTmp,a_dt/2);

  // RK Stage 3: compute new RHS.

  a_op.evalRHS(RHSTmp,solnTmp,
               a_fineFRPtr,a_crseFRPtr,
               a_oldCrseSoln, a_oldCrseTime,
               a_newCrseSoln, a_newCrseTime,
               a_time+a_dt/2,a_dt/3);
  a_interpPtr.saveRHS(RHSTmp);


  // RK Stage 3: compute update {phi}_3, partial update to new vals.

  a_op.updateODE(a_newSoln,RHSTmp,a_dt/3);

  // predict value at new time
  a_op.copySolnData(solnTmp,a_oldSoln);
  a_op.updateODE(solnTmp,RHSTmp,a_dt);

  // RK Stage 4: compute RHS.

  a_op.evalRHS(RHSTmp,solnTmp,
               a_fineFRPtr,a_crseFRPtr,
               a_oldCrseSoln, a_oldCrseTime,
               a_newCrseSoln, a_newCrseTime,
               a_time+a_dt,a_dt/6);
  a_interpPtr.saveRHS(RHSTmp);

  // RK Stage 4: compute final update of solution.

  a_op.updateODE(a_newSoln,RHSTmp,a_dt/6);
}

#endif // end multiple-include preventer

---