Setting up a new PeleLM Case

In order to set up and run a new case in PeleLM, the user must provide problem-specific code for two main tasks

Initial conditions
Boundary conditions

These functions are typically collected into a single subfolder in ${PELELM_HOME}/Exec, such as FlameSheet. The user can organize these tasks in any way that is convenient - the examples distributed with PeleLM represent a certain style for managing this with some level of flexibility, but the basic requirement is simply that source be linked into the build for the functions pelelm_initdata for initial conditions and bcnormal for boundary conditions.

Initial Conditions

At the beginning of a PeleLM run, for each level, after grids are generated, the cell-centered values of the state must be initialized. In the code, this is done in an MFIter loop over grids, and a call to the user’s initialization function, pelelm_initdata, that must be provided:

for (MFIter mfi(S_new,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
    const Box& box = mfi.validbox();
    auto sfab = S_new.array(mfi);

    amrex::ParallelFor(box,
    [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
    {
      pelelm_initdata(i, j, k, sfab, geomdata, *lprobparm, lpmfdata);
    });
}

where (i,j,k) are the cell indices, sfab is a light data pointer to the initial state MultiFab and geomdata, lprobparm and lpmfdata are container for the geometry, user-define input and PMF data. The associated user function (in pelelm_prob.H) will provide a value for each entry of the state (velocity, density, mass fraction, …) :

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
void
pelelm_initdata (int i, int j, int k,
                 amrex::Array4<amrex::Real> const& state,
                 amrex::GeometryData const& geomdata,
                 ProbParm const& prob_parm,
                 PmfData const *pmf_data)
{

   const amrex::Real* prob_lo = geomdata.ProbLo();
   const amrex::Real* prob_hi = geomdata.ProbHi();
   const amrex::Real* dx      = geomdata.CellSize();

   const amrex::Real z = prob_lo[2] + (k+0.5)*dx[2];
   const amrex::Real y = prob_lo[1] + (j+0.5)*dx[1];
   const amrex::Real x = prob_lo[0] + (i+0.5)*dx[0];

   constexpr amrex::Real Pi = 3.14159265358979323846264338327950288;
   const amrex::Real L_x = prob_hi[0] - prob_lo[0];
   const amrex::Real L_y = prob_hi[1] - prob_lo[1];


   ...


   state(i,j,k,DEF_Temp) = prob_parm.Temp;

   for (int n = 0; n < NUM_SPECIES; n++){
     massfrac[n] = 1.0/NUM_SPECIES;
   }


   state(i,j,k,Xvel) = 0.0;
   state(i,j,k,Yvel) =  prob_parm.Vel;
#elif ( AMREX_SPACEDIM == 3 )
   state(i,j,k,Zvel) = 0.0;
#endif

   amrex::Real rho_cgs, P_cgs;
   P_cgs = prob_parm.P_mean * 10.0;

   auto eos = pele::physics::PhysicsType::eos();
   eos.PYT2R(P_cgs, massfrac, state(i,j,k,DEF_Temp), rho_cgs);
   state(i,j,k,Density) = rho_cgs * 1.0e3;            // CGS -> MKS conversion

   eos.TY2H(state(i,j,k,DEF_Temp), massfrac, state(i,j,k,DEF_RhoH));
   state(i,j,k,DEF_RhoH) = state(i,j,k,DEF_RhoH) * 1.0e-4 * state(i,j,k,Density);   // CGS -> MKS conversion

   for (int n = 0; n < NUM_SPECIES; n++) {
     state(i,j,k,DEF_first_spec+n) = massfrac[n] * state(i,j,k,Density);
   }
}

Note home the geometry data are retrived from the geomdata object to obtain the coordinate of each cell center. The state data indices (Xvel, Yvel, Density, … ) are prescribed in the $(PELELM_HOME)/Source/IndexDefines.H file. Note that the conserved states are stored for species and enthalpy (i.e., $\rho Y_i$ and $\rho h$); these are the variables that the user must fill in the initial and boundary condition routines. Typically, however, the primitive state (i.e., $Y_i$ and $T$) is known directly. If that is the case, the user can make use of the compiled-in model-specific equation-of-state routines (eos.) to translate primitive to conserved state values. Consult the example setups provided to see how to call these routines, and how to load the final values required for initial data.

The runtime option (such as initial temperature, inlet velocity, …) are gathered in the ProbParm C++ structure defined in pelelm_prob_parm.H and filled from the input file in pelelm_prob.cpp using AMReX parser. This structure can be modified by the user to hold any data necessary for initial or boundary conditions.

Boundary Conditions

In PeleLM, a single function is used to fill all the state component at physical boundaries. The function bcnormal is in the pelelm_prob.H file. The main objective of this function is to fill the s_ext array fill boundary state data. The function bcnormal is called on each side (lo or hi) for each spatial dimension and will be used to fill the ghost cells of the state variables for which the PeleLM internal boundary condition is EXT_DIR (external Dirichlet) on that face. For example, specifying a PeleLM Inflow boundary condition on the lower face in the y direction in the input file leads to an EXT_DIR for species mass fraction, which then need to be provided in bcnormal. An example of the bcnormal of the FlameSheet is presented here:

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
void
bcnormal(
  const amrex::Real x[AMREX_SPACEDIM],
  amrex::Real s_ext[DEF_NUM_STATE],
  const int idir,
  const int sgn,
  const amrex::Real time,
  amrex::GeometryData const& geomdata,
  ProbParm const& prob_parm,
  ACParm const& ac_parm,
  PmfData const *pmf_data)
{
  const amrex::Real* prob_lo = geomdata.ProbLo();
  const amrex::Real* prob_hi = geomdata.ProbHi();
  amrex::GpuArray<amrex::Real, NUM_SPECIES + 4> pmf_vals = {0.0};
  amrex::Real molefrac[NUM_SPECIES] = {0.0};
  amrex::Real massfrac[NUM_SPECIES] = {0.0};

  if (sgn == 1) {
    PMF::pmf(pmf_data,prob_lo[idir], prob_lo[idir], pmf_vals);

    s_ext[Xvel] = 0.0;
#if ( AMREX_SPACEDIM == 2 )
    s_ext[Yvel] = pmf_vals[1]*1e-2;
#elif (AMREX_SPACEDIM == 3)
    s_ext[Yvel] = 0.0;
    s_ext[Zvel] = pmf_vals[1]*1e-2;
#endif

    s_ext[DEF_Temp] = pmf_vals[0];

    for (int n = 0; n < NUM_SPECIES; n++){
      molefrac[n] = pmf_vals[3 + n];
    }
    auto eos = pele::physics::PhysicsType::eos();
    eos.X2Y(molefrac, massfrac);

    amrex::Real rho_cgs, P_cgs, RhoH_temp;
    P_cgs = prob_parm.P_mean * 10.0;

    eos.PYT2R(P_cgs, massfrac, s_ext[DEF_Temp], rho_cgs);
    s_ext[Density] = rho_cgs * 1.0e3;

    eos.TY2H(s_ext[DEF_Temp], massfrac, RhoH_temp);
    s_ext[DEF_RhoH] = RhoH_temp * 1.0e-4 * s_ext[Density];   // CGS -> MKS conversion

    for (int n = 0; n < NUM_SPECIES; n++) {
      s_ext[DEF_first_spec+n] = massfrac[n] * s_ext[Density];
    }
  }
}

The sgn input takes a value of 1 on the low face and -1 on the high face, while ìdir provide the spatial direction (0, 1 or 2 corresponding to X, Y or Z, respectively). This allow to differentiate between the various boundary conditions when more than 1 ÈXT_DIR is needed. In this example, the boundary conditions are extracted from a pre-computed premixed flame which data are stored in the pmf_data structure.

Here, we’ve made use of a local convenience function, bcnormal endowed with the knowledge of all boundary values, and extract the appropriate quantity from the results of that call. This was done to localize all boundary condition calculations to a single routine in the code, and helps to preserve consistency. This is only one style though, and as long as appropriate Dirichlet values are set for this state, it makes no difference how the work is organized. For example, data may be provided by interpolating “live data” being actively generated by a co-running separate code, by interpolating data files, evaluating functional forms, etc.

Note that although the array structure to be filled contains valid cell-centered state data where it overlaps the valid domain, the values set in the grow cells of the container will be applied on the boundary face of the corresponding cells. Internally, all PeleLM code understands to apply Dirichlet conditions on the boundary faces.

Setting up a new PeleLM Case

Initial Conditions

Boundary Conditions

Refinement Criteria