Method gray_accel

Source/radiation/MGFLD.cpp:440–701

void Radiation::gray_accel(MultiFab& Er_new, MultiFab& Er_pi,
                           MultiFab& kappa_p, MultiFab& kappa_r,
                           MultiFab& etaT, MultiFab& eta1,
                           MultiFab& mugT,
                           Array<MultiFab, AMREX_SPACEDIM>& lambda,
                           RadSolve* solver, MGRadBndry& mgbd,
                           const BoxArray& grids, int level, Real time,
                           Real delta_t, Real ptc_tau)
{
  const Geometry& geom = parent->Geom(level);
  auto dx = parent->Geom(level).CellSizeArray();
  const Castro *castro = dynamic_cast<Castro*>(&parent->getLevel(level));
  const DistributionMapping& dmap = castro->DistributionMap();

  if (nGroups > 1) {
    solver->setHypreMulti(1.0);
  }
  else {
    solver->setHypreMulti(0.0);
  }
  mgbd.setCorrection();

  MultiFab Er_zero(grids, dmap, 1, 0);
  Er_zero.setVal(0.0);
  getBndryDataMG_ga(mgbd, Er_zero, level);

  MultiFab spec(grids, dmap, nGroups, 1);

#ifdef _OPENMP
#pragma omp parallel
#endif
  for (MFIter mfi(spec, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
    const Box& bx = mfi.tilebox();

    auto kappa_p_arr = kappa_p[mfi].array();
    auto mugT_arr = mugT[mfi].array();
    auto spec_arr = spec[mfi].array();

    Real cdt1 = 1.e0_rt / (C::c_light * delta_t);

    amrex::ParallelFor(bx,
    [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
    {
        Real epsilon[NGROUPS];
        Real sumeps = 0.0;
        for (int g = 0; g < NGROUPS; ++g) {
            Real kapt = kappa_p_arr(i,j,k,g) + (1.e0_rt + ptc_tau) * cdt1;
            epsilon[g] = mugT_arr(i,j,k,g) / kapt;
            sumeps += epsilon[g];
        }

        if (sumeps == 0.e0_rt) {
            for (int g = 0; g < NGROUPS; ++g) {
                spec_arr(i,j,k,g) = 0.e0_rt;
            }
        } else {
            for (int g = 0; g < NGROUPS; ++g) {
                spec_arr(i,j,k,g) = epsilon[g] / sumeps;