releasedExamples/AMRINS/NavierStokes/AMRNavierStokesAdvance.cpp
// ---------------------------------------------------------------
// advance solution by one timestep
// return max safe timestep
advance solution by one timestep
return max safe timestep
advance solution by one timestep return max safe timestep
Real
AMRNavierStokes::advance()
{
  if (s_verbosity >= 2)
  {
    pout () << "AMRNavierStokes::advance " << m_level
            << ", starting time = "
            << setiosflags(ios::fixed) << setprecision(6)
            << setw(12) << m_time
            << ", dt = "
            << setiosflags(ios::fixed) << setprecision(6) << setw(12) << dt()
            << endl;
  }
#ifdef DEBUG
  if (m_level == 0)
    {
      AMRNavierStokes* levelPtr = this;
      while (!levelPtr->finestLevel())
        {
          levelPtr = levelPtr->finerNSPtr();
        }
    int finest_level = levelPtr->m_level;
    levelPtr = this;
    for (int lev=0; lev<= finest_level; lev++)
      {
        levelPtr->m_saved_time = m_time;
        levelPtr->newVel().copyTo(levelPtr->newVel().interval(),
                                  *levelPtr->m_vel_save_ptr,
                                  levelPtr->m_vel_save_ptr->interval());
        levelPtr->newLambda().copyTo(levelPtr->newLambda().interval(),
                                     *levelPtr->m_lambda_save_ptr,
                                     levelPtr->m_lambda_save_ptr->interval());
        for (int comp=0; comp<s_num_scal_comps; comp++)
          {
            levelPtr->newScal(comp).copyTo(levelPtr->newScal(comp).interval(),
                                           *levelPtr->m_scal_save[comp],
                                           levelPtr->m_scal_save[comp]->interval());
          }
        levelPtr = levelPtr->finerNSPtr();
      }
    }
#endif
  Real old_time = m_time;
  Real new_time = m_time + m_dt;
  m_dt_save = m_dt;
  swapOldAndNewStates();
  const DisjointBoxLayout levelGrids = newVel().getBoxes();
  if (!finestLevel())
    {
      m_flux_register.setToZero();
      for (int comp=0; comp<s_num_scal_comps; comp++)
        {
          m_scal_fluxreg_ptrs[comp]->setToZero();
        }
      m_lambda_flux_reg.setToZero();
    }
  IntVect advVelGhost(IntVect::Unit);
  LevelData<FluxBox> advectionVel(levelGrids, 1, advVelGhost);
  if (s_set_bogus_values)
    {
      setValLevel(advectionVel, s_bogus_value);
    }
  computeAdvectionVelocities(advectionVel);
  LevelData<FArrayBox> lambdaTemp;
  fillLambda(lambdaTemp, old_time);
  LevelFluxRegister* crseLambdaFluxRegPtr = NULL;
  if (m_level > 0)
    {
      crseLambdaFluxRegPtr = &(crseNSPtr()->m_lambda_flux_reg);
    }
  advectScalar(newLambda(), lambdaTemp, advectionVel,
               crseLambdaFluxRegPtr, m_lambda_flux_reg,
               m_patchGodLambda,
               m_dt);
  if (s_num_scal_comps > 0)
    {
      for (int comp=0; comp<s_num_scal_comps; comp++)
        {
          BCHolder scalPhysBCs = m_physBCPtr->scalarTraceFuncBC(comp);
          LevelData<FArrayBox> scalarTemp;
          fillScalars(scalarTemp, old_time,comp);
          LevelData<FArrayBox>* newCrseScalPtr = NULL;
          LevelData<FArrayBox>* oldCrseScalPtr = NULL;
          Real oldCrseTime = -s_bogus_value;
          Real newCrseTime = s_bogus_value;
          LevelFluxRegister* crseScalFluxRegPtr = NULL;
          if (m_level > 0)
            {
              newCrseScalPtr = &(crseNSPtr()->newScal(comp));
              oldCrseScalPtr = &(crseNSPtr()->oldScal(comp));
              newCrseTime = crseNSPtr()->time();
              oldCrseTime = newCrseTime - crseNSPtr()->dt();
              crseScalFluxRegPtr = (crseNSPtr()->m_scal_fluxreg_ptrs[comp]);
            }
          advectDiffuseScalar(newScal(comp), scalarTemp,
                              advectionVel,
                              s_scal_coeffs[comp],
                              oldCrseScalPtr, newCrseScalPtr,
                              oldCrseTime, newCrseTime,
                              crseScalFluxRegPtr,
                              *(m_scal_fluxreg_ptrs[comp]),
                              *(m_patchGodScalars[comp]),
                              scalPhysBCs,
                              m_dt, comp);
        } 
    } 
  LevelData<FArrayBox>& uStar = newVel();
  predictVelocities(uStar, advectionVel);
  computeUStar(uStar);
  LevelData<FArrayBox>* crseVelPtr = NULL;
  if (m_level > 0)
    {
      const DisjointBoxLayout& crseGrids = crseNSPtr()->newVel().getBoxes();
      crseVelPtr = new LevelData<FArrayBox>(crseGrids, SpaceDim);
      crseNSPtr()->fillVelocity(*crseVelPtr, new_time);
    }
  bool isViscous = (s_nu > 0);
  VelBCHolder velBC(m_physBCPtr->uStarFuncBC(isViscous));
  velBC.applyBCs(uStar, levelGrids,
                 m_problem_domain, m_dx,
                 false); 
  m_projection.LevelProject(uStar, crseVelPtr, new_time, m_dt);
  Real newDt = computeDt();
  if (crseVelPtr != NULL)
    {
      delete crseVelPtr;
      crseVelPtr = NULL;
    }
  ++m_level_steps;
  if (s_verbosity >= 2)
    {
      pout () << "AMRNavierStokes::advance " << m_level
              << ",      end time = "
              << setiosflags(ios::fixed) << setprecision(6) << setw(12) << m_time
              << ", dt = "
              << setiosflags(ios::fixed) << setprecision(6)
              << setw(12) << dt()
              << endl;
    }
  return newDt;
}
Real AMRNavierStokes :: advance ( ) { if ( s_verbosity >= 2 ) { pout ( ) << " " << m_level << " " << setiosflags ( ios :: fixed ) << setprecision ( 6 ) << setw ( 12 ) << m_time << " " << setiosflags ( ios :: fixed ) << setprecision ( 6 ) << setw ( 12 ) << dt ( ) << endl ; } #ifdef DEBUG if ( m_level == 0 ) { AMRNavierStokes * levelPtr = this ; while ( ! levelPtr -> finestLevel ( ) ) { levelPtr = levelPtr -> finerNSPtr ( ) ; } int finest_level = levelPtr -> m_level ; levelPtr = this ; for ( int lev = 0 ; lev <= finest_level ; lev ++ ) { levelPtr -> m_saved_time = m_time ; levelPtr -> newVel ( ) . copyTo ( levelPtr -> newVel ( ) . interval ( ) , * levelPtr -> m_vel_save_ptr , levelPtr -> m_vel_save_ptr -> interval ( ) ) ; levelPtr -> newLambda ( ) . copyTo ( levelPtr -> newLambda ( ) . interval ( ) , * levelPtr -> m_lambda_save_ptr , levelPtr -> m_lambda_save_ptr -> interval ( ) ) ; for ( int comp = 0 ; comp < s_num_scal_comps ; comp ++ ) { levelPtr -> newScal ( comp ) . copyTo ( levelPtr -> newScal ( comp ) . interval ( ) , * levelPtr -> m_scal_save [ comp ] , levelPtr -> m_scal_save [ comp ] -> interval ( ) ) ; } levelPtr = levelPtr -> finerNSPtr ( ) ; } } #endif Real old_time = m_time ; Real new_time = m_time + m_dt ; m_dt_save = m_dt ; swapOldAndNewStates ( ) ; const DisjointBoxLayout levelGrids = newVel ( ) . getBoxes ( ) ; if ( ! finestLevel ( ) ) { m_flux_register . setToZero ( ) ; for ( int comp = 0 ; comp < s_num_scal_comps ; comp ++ ) { m_scal_fluxreg_ptrs [ comp ] -> setToZero ( ) ; } m_lambda_flux_reg . setToZero ( ) ; } IntVect advVelGhost ( IntVect :: Unit ) ; LevelData < FluxBox > advectionVel ( levelGrids , 1 , advVelGhost ) ; if ( s_set_bogus_values ) { setValLevel ( advectionVel , s_bogus_value ) ; } computeAdvectionVelocities ( advectionVel ) ; LevelData < FArrayBox > lambdaTemp ; fillLambda ( lambdaTemp , old_time ) ; LevelFluxRegister * crseLambdaFluxRegPtr = NULL ; if ( m_level > 0 ) { crseLambdaFluxRegPtr = & ( crseNSPtr ( ) -> m_lambda_flux_reg ) ; } advectScalar ( newLambda ( ) , lambdaTemp , advectionVel , crseLambdaFluxRegPtr , m_lambda_flux_reg , m_patchGodLambda , m_dt ) ; if ( s_num_scal_comps > 0 ) { for ( int comp = 0 ; comp < s_num_scal_comps ; comp ++ ) { BCHolder scalPhysBCs = m_physBCPtr -> scalarTraceFuncBC ( comp ) ; LevelData < FArrayBox > scalarTemp ; fillScalars ( scalarTemp , old_time , comp ) ; LevelData < FArrayBox > * newCrseScalPtr = NULL ; LevelData < FArrayBox > * oldCrseScalPtr = NULL ; Real oldCrseTime = - s_bogus_value ; Real newCrseTime = s_bogus_value ; LevelFluxRegister * crseScalFluxRegPtr = NULL ; if ( m_level > 0 ) { newCrseScalPtr = & ( crseNSPtr ( ) -> newScal ( comp ) ) ; oldCrseScalPtr = & ( crseNSPtr ( ) -> oldScal ( comp ) ) ; newCrseTime = crseNSPtr ( ) -> time ( ) ; oldCrseTime = newCrseTime - crseNSPtr ( ) -> dt ( ) ; crseScalFluxRegPtr = ( crseNSPtr ( ) -> m_scal_fluxreg_ptrs [ comp ] ) ; } advectDiffuseScalar ( newScal ( comp ) , scalarTemp , advectionVel , s_scal_coeffs [ comp ] , oldCrseScalPtr , newCrseScalPtr , oldCrseTime , newCrseTime , crseScalFluxRegPtr , * ( m_scal_fluxreg_ptrs [ comp ] ) , * ( m_patchGodScalars [ comp ] ) , scalPhysBCs , m_dt , comp ) ; } } LevelData < FArrayBox > & uStar = newVel ( ) ; predictVelocities ( uStar , advectionVel ) ; computeUStar ( uStar ) ; LevelData < FArrayBox > * crseVelPtr = NULL ; if ( m_level > 0 ) { const DisjointBoxLayout & crseGrids = crseNSPtr ( ) -> newVel ( ) . getBoxes ( ) ; crseVelPtr = new LevelData < FArrayBox > ( crseGrids , SpaceDim ) ; crseNSPtr ( ) -> fillVelocity ( * crseVelPtr , new_time ) ; } bool isViscous = ( s_nu > 0 ) ; VelBCHolder velBC ( m_physBCPtr -> uStarFuncBC ( isViscous ) ) ; velBC . applyBCs ( uStar , levelGrids , m_problem_domain , m_dx , false ) ; m_projection . LevelProject ( uStar , crseVelPtr , new_time , m_dt ) ; Real newDt = computeDt ( ) ; if ( crseVelPtr != NULL ) { delete crseVelPtr ; crseVelPtr = NULL ; } ++ m_level_steps ; if ( s_verbosity >= 2 ) { pout ( ) << " " << m_level << " " << setiosflags ( ios :: fixed ) << setprecision ( 6 ) << setw ( 12 ) << m_time << " " << setiosflags ( ios :: fixed ) << setprecision ( 6 ) << setw ( 12 ) << dt ( ) << endl ; } return newDt ; }
advance solution by one timestep
return max safe timestep
advance solution by one timestep return max safe timestep
// if we're at the beginning of a composite timestep, // save current multilevel state // figure out finest level // initialize flux registers // compute advection velocities -- if we're using the // patchGodunov approach, these will have ghost cells. // now that advection velocities have been computed, // update advected scalars // do scalar boundary conditions on lambda // (also grow to appropriate size) // now do advected-diffused scalars // loop over scalar components // now get coarse-level scalars for BC's if necessary // end loop over components // end advected-diffused scalars // now predict velocities -- this returns the advection term // put this in "new_vel" // if a coarser level exists, will need coarse-level data for proj // coarse velocity BC data is interpolated in time // need to do physical boundary conditions and exchanges // inhomogeneous // noel -- all time in level project // as things stand now, physical BC's are re-set in LevelProjection // compute maximum safe timestep for next iteration // clean up temp storage