Method levelDCoeffs

Source/radiation/RadSolve.cpp:413–502

}

void RadSolve::levelDCoeffs(int level, Array<MultiFab, AMREX_SPACEDIM>& lambda,
                            MultiFab& vel, MultiFab& dcf)
{
    BL_PROFILE("RadSolve::levelDCoeffs");
    const Castro *castro = dynamic_cast<Castro*>(&parent->getLevel(level));
    const DistributionMapping& dm = castro->DistributionMap();
    const Geometry& geom = parent->Geom(level);
    const auto dx = geom.CellSizeArray();
    const auto geomdata = geom.data();

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

        MultiFab dcoefs(castro->getEdgeBoxArray(idim), dm, 1, 0);

#ifdef _OPENMP
#pragma omp parallel
#endif
        for (MFIter mfi(dcoefs, TilingIfNotGPU()); mfi.isValid(); ++mfi) {

            const Box& bx = mfi.tilebox();

            auto dcoefs_arr = dcoefs[mfi].array();
            auto lambda_arr = lambda[idim][mfi].array();
            auto vel_arr = vel[mfi].array();
            auto dcf_arr = dcf[mfi].array();

            amrex::ParallelFor(bx,
            [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
            {
                Real r, s;

                if (idim == 0) {

                    edge_center_metric(i, j, k, idim, geomdata, r, s);

                    if (vel_arr(i-1,j,k,0) + vel_arr(i,j,k,0) > 0.e0_rt) {
                        dcoefs_arr(i,j,k) = dcf_arr(i-1,j,k) * vel_arr(i-1,j,k,0) * lambda_arr(i,j,k);
                    }
                    else if (vel_arr(i-1,j,k,0) + vel_arr(i,j,k,0) < 0.e0_rt) {
                        dcoefs_arr(i,j,k) = dcf_arr(i,j,k) * vel_arr(i,j,k,0) * lambda_arr(i,j,k);
                    }
                    else {
                        dcoefs_arr(i,j,k) = 0.e0_rt;
                    }

                    dcoefs_arr(i,j,k) = dcoefs_arr(i,j,k) * r;

                }
                else if (idim == 1) {

                    edge_center_metric(i, j, k, idim, geomdata, r, s);

                    if (vel_arr(i,j-1,k,1) + vel_arr(i,j,k,1) > 0.e0_rt) {
                        dcoefs_arr(i,j,k) = dcf_arr(i,j-1,k) * vel_arr(i,j-1,k,1) * lambda_arr(i,j,k);
                    }
                    else if (vel_arr(i,j-1,k,1) + vel_arr(i,j,k,1) < 0.e0_rt) {
                        dcoefs_arr(i,j,k) = dcf_arr(i,j,k) * vel_arr(i,j,k,1) * lambda_arr(i,j,k);
                    }