Accessing flow quantities (viscosity, diffusion coefficients) in a derived quantity

[JayStandridge]

Hello,

I am trying to implement a local form of the entropy inequality into PeleC. Nick Wimer created a branch with a second order velocity gradient calculation which has been very useful. In Setup.cpp the routine (and mine) are added to the derive_lst with amrex::DeriveRec::GrowBoxByOne ( as opposed to amres:DeriveRec:TheSameBox for most other derive functions) in order to accomplish the central difference for the velocity gradient tensor. Now, I need the viscosity, diffusion coefficients, and thermal conductivity. There are derive functions already in place for them, but I am unable to get the code to compile when trying to access trans_parms.device_trans_parm() as they do because it wasn't declared in that scope. Is there a way for me to access these flow variables in my routine?

Thanks in advance.

[marchdf]

Hi, thanks for reaching out! The first thing I would recommend is for @nickwimer to make a PR so we can merge that code into mainline PeleC since it seems useful in general. The second thing is: can I see the code? It's hard to debug this without code. My first hunch is that maybe they aren't part of PeleC:: and so they can't access that member variable?

[JayStandridge]

Of course!
Here is my code from Setup.cpp:

//   // Entropy Inequality
   amrex::Vector<std::string> var_names_EI(9);
   var_names_EI[0] = "EITerm1";
   var_names_EI[1] = "EITerm2";
   var_names_EI[2] = "EITerm3";
   var_names_EI[3] = "EITerm4";
   var_names_EI[4] = "EI";
   var_names_EI[5] = "AUX1";
   var_names_EI[6] = "AUX2";
   var_names_EI[7] = "AUX3";
   var_names_EI[8] = "AUX4";

   derive_lst.add(
     "entropyInequality", amrex::IndexType::TheCellType(), 9, var_names_EI,
     pc_entropyInequality, amrex::DeriveRec::GrowBoxByOne);

     // derive_lst.add(
     // "entropyInequality", amrex::IndexType::TheCellType(), 1,
     // pc_entropyInequality, amrex::DeriveRec::TheSameBox);
   derive_lst.addComponent("entropyInequality", desc_lst, State_Type, Density, NVAR);

From IO.cpp:

bool plot_ei = true;
  pp.query("plot_ei", plot_ei);
  if (plot_ei) {
    amrex::Amr::addDerivePlotVar("entropyInequality");
  } else {
    amrex::Amr::deleteDerivePlotVar("entropyInequality");
  }

From Derive.H:

void pc_entropyInequality(
  const amrex::Box& bx,
  amrex::FArrayBox& derfab,
  int dcomp,
  int ncomp,
  const amrex::FArrayBox& datfab,
  const amrex::Geometry& geomdata,
  amrex::Real time,
  const int* bcrec,
  const int level);

And from Derive.cpp:

void
pc_entropyInequality(
  const amrex::Box& bx,
  amrex::FArrayBox& derfab,
  int /*dcomp*/,
  int /*ncomp*/,
  const amrex::FArrayBox& datfab,
  const amrex::Geometry& geomdata,
  amrex::Real /*time*/,
  const int* /*bcrec*/,
  const int /*level*/)
{

   
  auto entropyInequality = derfab.array();
  auto vel_ders = derfab.array();
  auto const dat = datfab.const_array();
  auto mu_arr = derfab.array();
  auto divu   = derfab.array();

  // This section stolen from pc_derviscosity, code does not compile with it in
  // auto const* ltransparm = trans_parms.device_trans_parm();
  // amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
  //   amrex::Real massfrac[NUM_SPECIES];
  //   const amrex::Real rho = dat(i, j, k, URHO);
  //   const amrex::Real rhoInv = 1.0 / rho;
  //   const amrex::Real T = dat(i, j, k, UTEMP);

  //   for (int n = 0; n < NUM_SPECIES; n++) {
  //     massfrac[n] = dat(i, j, k, UFS + n) * rhoInv;
  //   }
  //   auto trans = pele::physics::PhysicsType::transport();
  //   amrex::Real mu = 0.0, dum1 = 0.0, dum2 = 0.0;
  //   const bool get_xi = false, get_mu = true, get_lam = false,
  //              get_Ddiag = false;
  //   trans.transport(
  //     get_xi, get_mu, get_lam, get_Ddiag, T, rho, massfrac, nullptr, mu, dum1,
  //     dum2, ltransparm);
  //   mu_arr(i, j, k) = mu;
  // });
  
 // From @nickwimer, calculates velocity gradient tensor using central differences
  const amrex::Box& gbx = amrex::grow(bx, 1);

  amrex::FArrayBox local(gbx, 3, amrex::The_Async_Arena());
  auto larr = local.array();

  const auto& flag_fab = amrex::getEBCellFlagFab(datfab);
  const auto& typ = flag_fab.getType(bx);
  if (typ == amrex::FabType::covered) {
    derfab.setVal<amrex::RunOn::Device>(0.0, bx);
    return;
  }
  const auto& flags = flag_fab.const_array();
  const bool all_regular = typ == amrex::FabType::regular;
  // Convert momentum to velocity.
  amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
    const amrex::Real rhoInv = 1.0 / dat(i, j, k, URHO);
    larr(i, j, k, 0) = dat(i, j, k, UMX) * rhoInv;
    larr(i, j, k, 1) = dat(i, j, k, UMY) * rhoInv;
    larr(i, j, k, 2) = dat(i, j, k, UMZ) * rhoInv;
  });
  AMREX_D_TERM(const amrex::Real dx = geomdata.CellSize(0);
               , const amrex::Real dy = geomdata.CellSize(1);
               , const amrex::Real dz = geomdata.CellSize(2););

    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
  			   AMREX_D_TERM(int im; int ip;, int jm; int jp;, int km; int kp;);
  			   AMREX_D_TERM(get_idx(i, 0, all_regular, flags(i, j, k), im, ip);
  			   		, get_idx(j, 1, all_regular, flags(i, j, k), jm, jp);
  			   		, get_idx(k, 2, all_regular, flags(i, j, k), km, kp););
  			   AMREX_D_TERM(const amrex::Real wi = get_weight(im, ip);
  			   		, const amrex::Real wj = get_weight(jm, jp);
  			   		, const amrex::Real wk = get_weight(km, kp););

  			   vel_ders(i, j, k, 0) = wi * (larr(ip, j, k, 0) - larr(im, j, k, 0)) / dx; // dudx
  			   vel_ders(i, j, k, 1) = wj * (larr(i, jp, k, 0) - larr(i, jm, k, 0)) / dy; // dudy
  			   vel_ders(i, j, k, 2) = wk * (larr(i, j, kp, 0) - larr(i, j, km, 0)) / dz; // dudz

  			   vel_ders(i, j, k, 3) = wi * (larr(ip, j, k, 1) - larr(im, j, k, 1)) / dx; // dvdx
  			   vel_ders(i, j, k, 4) = wj * (larr(i, jp, k, 1) - larr(i, jm, k, 1)) / dy; // dvdy
  			   vel_ders(i, j, k, 5) = wk * (larr(i, j, kp, 1) - larr(i, j, km, 1)) / dz; // dvdz

  			   vel_ders(i, j, k, 6) = wi * (larr(ip, j, k, 2) - larr(im, j, k, 2)) / dx; // dwdx
  			   vel_ders(i, j, k, 7) = wj * (larr(i, jp, k, 2) - larr(i, jm, k, 2)) / dy; // dwdy
  			   vel_ders(i, j, k, 8) = wk * (larr(i, j, kp, 2) - larr(i, j, km, 2)) / dz; // dwdz
 
                           // Calculate one of the quantities of interest to me
			   divu(i,j,k,0) = vel_ders(i, j, k, 0) + vel_ders(i, j, k, 4) + vel_ders(i, j, k, 8);
			   entropyInequality(i,j,k,0) = (2/3*pow(divu(i,j,k,0),2)
                            - 2*( pow(vel_ders(i, j, k, 0),2) + pow(vel_ders(i, j, k, 4),2) +pow( vel_ders(i, j, k, 8),2))
                            - (pow(vel_ders(i, j, k, 2) + vel_ders(i, j, k, 3),2) + pow(vel_ders(i, j, k, 2) 
                            + vel_ders(i, j, k, 6),2) + pow(vel_ders(i, j, k, 7) + vel_ders(i, j, k, 5),2)));
    });
   }

Thank you for the help

[marchdf]

Ok yeah the problem is that they aren't member functions of PeleC class. You are going to have to follow the way (e.g.) diffusivity is done:https://github.com/AMReX-Combustion/PeleC/blob/development/Source/PeleC.H#L541 and https://github.com/AMReX-Combustion/PeleC/blob/development/Source/Derive.cpp#L847. Basically move the pc_entropyInequality declaration from Derive.H to PeleC.H and then add PeleC:: to the function definition.

[JayStandridge]

Okay, that makes sense. I hadn't noticed that those declarations weren't in Derive.H.

Thank you!