Modifying checkpoint data for restart

[dashp7]

Hi,
Is there a way of modifying checkpoint data in PeleC? The need arises out of the fact that the function pc_initdata() is not called when a simulation is resumed from a checkpoint using amr.restart=chk#####. Basically, any approach towards reading in checkpoint data as a state Array4 like it appears in pc_initdata() can work. Any suggestion on a starting point for this functionality using AMReX features like ReadCheckpointFile() and WriteCheckpointFile() will be useful.
Thanks.

[marchdf]

Hi, the user can access the PeleC member function problem_post_restart(). It is defined empty in almost all problems right now. Presumably, one could grab the state multifab in there, and then write an MFIter loop and modify the state. That might work.

[dashp7]

@marchdf thanks for the reply.
Can the function problem_post_restart() be defined along the lines of pc_initdata() where the state variable(Array4) can be accessed using integral indices?

[RSuryaNarayan]

@marchdf would this method of grabbing the state work

PeleC/Source/PeleC.cpp

Line 627 in 0fbdc11

amrex::MultiFab& S_new = get_new_data(State_Type);

[marchdf]

I don't think it's necessary to change the API. I would define that function that does something like

PeleC::problem_post_restart(){ 
   amrex::MultiFab& S_new = get_new_data(State_Type);
   
#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
  for (amrex::MFIter mfi(S_new, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
    const amrex::Box& bx = mfi.tilebox();
auto const& Sarr = S_new.array(mfi);
    amrex::ParallelFor(
      bx, NVAR, [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
        Sarr(i, j, k, n) = 0.0; // set to whatever
      });
  }
}

You can grab things like geomdata in there if you want too.

That code probably doesn't compile... but something like that should do it

[RSuryaNarayan]

@marchdf thanks a ton! was very helpful

[RSuryaNarayan]

Hi @marchdf,

Based on your inputs, I worked out a more digestible and general version of this function. Just posting it for the clarity of other users.

void
PeleC::problem_post_restart()
{
  
  amrex::MultiFab& S_new = get_new_data(State_Type);
  for (amrex::MFIter mfi(S_new, amrex::TilingIfNotGPU()); mfi.isValid();
       ++mfi) {
    const amrex::Box& box = mfi.tilebox();
    auto state = S_new.array(mfi);
    const auto geomdata = geom.data();
    
    const amrex::Real* prob_lo = geomdata.ProbLo();
    const amrex::Real* dx = geomdata.CellSize();
    
    const ProbParmDevice* lprobparm = d_prob_parm_device;
#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
    amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
	// compute the local coordinates	
	amrex::Real x[3];
	x[0] = prob_lo[0] + (i + 0.5) * dx[0];
    	x[1] = prob_lo[1] + (j + 0.5) * dx[1];
	x[2] = prob_lo[2] + (k + 0.5) * dx[2];
       pc_post_restart(i, j, k, x, state, lprobparm);
    });
  }
}

Where we define the function pc_post_restart in prob.H like pc_initdata as:

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
void
pc_post_restart(int i, int j, int k, 
  const amrex::Real x[AMREX_SPACEDIM],
  const amrex::Array4<amrex::Real>& state,
  const ProbParmDevice* const prob_parm)
{
// custom function on the lines of pc_initdata
}

Example: let's say in EB-Bluffbody I want to define a rectangular region of high-temperature, I do:

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
void
pc_post_restart(int i, int j, int k, 
  const amrex::Real x[AMREX_SPACEDIM],
  const amrex::Array4<amrex::Real>& state,
  const ProbParmDevice* const prob_parm)
{
    if (x[0]>7.0 && x[0]<10.0 && x[1]>2 && x[1]<4)
      {
        std::cout<<"\n I am here!";
        amrex::Real p = 1013250.0;
        amrex::Real T = 1000; //set a large temperature value
        amrex::Real rho, e;
        auto eos = pele::physics::PhysicsType::eos();
  	eos.PYT2RE(
    	p, PeleC::h_prob_parm_device->massfrac.begin(),
    	T, rho, e);
    	state(i, j, k, UTEMP) = T;
    	state(i, j, k, URHO) = rho;
    	state(i, j, k, UEINT) = rho*e;
    	state(i, j, k, UEDEN) = rho*
  			  	(e + 0.5 * (prob_parm->vx_in * prob_parm->vx_in +
                             	prob_parm->vy_in * prob_parm->vy_in));
	}
}

Thanks once again!

[marchdf]

Super cool! I am glad it is working for you. Thanks for posting a better/working example.

[dashp7]

Once again, thanks for the valuable insights!

[Nadexterbrown]

Would this be able to be used for domain modification from the checkpoint file, such as increasing the length of rectangular domain? I understand that there is further nuances that would required to ensure the domain modification is physically reasonable and doesn't place the original solution in undesired locations.

[baperry2]

Domain modification when restarting from a checkpoint file is not supported. The physical size and base resolution must match for the data structures to be created and data loaded from the file correctly. The init from plotfile capability (pelec.init_plotfile = plt00100) offers a bit more flexibility and may be a better starting point. But based on some assumptions buried within the PltFileManager utility from PelePhysics, I'm not sure if that would let you go to bigger domains (especially if you have EB), and there's still the issue of needing to populate physically correct data into the "new" parts of the domain if it is being enlarged.

This would definitely be very useful capability, but I'm not sure if anyone has figured out a solution for it yet.

[Nadexterbrown]

thank you so much for your quick response. I was hoping maybe there was something using the checkpoint file, since it has the complete solution before attempting with the plot file approach. I will look into the plot file approach and see if there is something I can figure out from that approach. Thank you again and have a wonderful day.