Method levelFlux

Source/radiation/RadSolve.cpp:688–777

}

void RadSolve::levelFlux(int level,
                         Array<MultiFab, AMREX_SPACEDIM>& Flux,
                         MultiFab& Er, int igroup)
{
  BL_PROFILE("RadSolve::levelFlux");
  const BoxArray& grids = parent->boxArray(level);
  const DistributionMapping& dmap = parent->DistributionMap(level);

  // grow a larger MultiFab to hold Er so we can difference across faces
  MultiFab Erborder(grids, dmap, 1, 1);
  Erborder.setVal(0.0);
  MultiFab::Copy(Erborder, Er, igroup, 0, 1, 0);

  Erborder.FillBoundary(parent->Geom(level).periodicity()); // zeroes left in off-level boundaries

  auto dx = parent->Geom(level).CellSizeArray();

  for (int n = 0; n < AMREX_SPACEDIM; n++) {

      const MultiFab *bp;

      if (hd) {
          bp = &hd->bCoefficients(n);
      }
      else if (hm) {
          bp = &hm->bCoefficients(level, n);
      }
      else if (hem) {
          bp = &hem->bCoefficients(level, n);
      }

      MultiFab &bcoef = *(MultiFab*)bp;

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

          auto Erborder_arr = Erborder[mfi].array();
          auto bcoef_arr = bcoef[mfi].array();
          auto Flux_arr = Flux[n][mfi].array();

          Real beta = radsolve::beta;

          amrex::ParallelFor(bx,
          [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
          {
              if (n == 0) {

                  const Real fac = -beta / dx[0];

                  Flux_arr(i,j,k) = bcoef_arr(i,j,k) * (Erborder_arr(i,j,k) - Erborder_arr(i-1,j,k)) * fac;

              }
              else if (n == 1) {

                  const Real fac = -beta / dx[1];