| 806 | } |
| 807 | |
| 808 | void RadSolve::levelDterm(int level, MultiFab& Dterm, MultiFab& Er, int igroup) |
| 809 | { |
| 810 | BL_PROFILE("RadSolve::levelDterm"); |
| 811 | const BoxArray& grids = parent->boxArray(level); |
| 812 | const DistributionMapping& dmap = parent->DistributionMap(level); |
| 813 | const Geometry& geom = parent->Geom(level); |
| 814 | const GeometryData& geomdata = geom.data(); |
| 815 | auto dx = parent->Geom(level).CellSizeArray(); |
| 816 | const Castro *castro = dynamic_cast<Castro*>(&parent->getLevel(level)); |
| 817 | |
| 818 | Array<MultiFab, AMREX_SPACEDIM> Dterm_face; |
| 819 | for (int idim=0; idim<AMREX_SPACEDIM; idim++) { |
| 820 | Dterm_face[idim].define(castro->getEdgeBoxArray(idim), dmap, 1, 0); |
| 821 | } |
| 822 | |
| 823 | // grow a larger MultiFab to hold Er so we can difference across faces |
| 824 | MultiFab Erborder(grids, dmap, 1, 1); |
| 825 | Erborder.setVal(0.0); |
| 826 | MultiFab::Copy(Erborder, Er, igroup, 0, 1, 0); |
| 827 | |
| 828 | Erborder.FillBoundary(parent->Geom(level).periodicity()); // zeroes left in off-level boundaries |
| 829 | |
| 830 | #ifdef _OPENMP |
| 831 | #pragma omp parallel |
| 832 | #endif |
| 833 | for (int n = 0; n < AMREX_SPACEDIM; n++) { |
| 834 | const MultiFab *dp; |
| 835 | |
| 836 | dp = &hem->d2Coefficients(level, n); |
| 837 | MultiFab &dcoef = *(MultiFab*)dp; |
| 838 | |
| 839 | for (MFIter fi(dcoef,true); fi.isValid(); ++fi) { |
| 840 | const Box& bx = fi.tilebox(); |
| 841 | |
| 842 | auto Er = Erborder[fi].array(); |
| 843 | auto dc = dcoef[fi].array(); |
| 844 | auto dtf = Dterm_face[n][fi].array(); |
| 845 | |
| 846 | amrex::ParallelFor(bx, |
| 847 | [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) |
| 848 | { |
| 849 | if (n == 0) { |
| 850 | dtf(i,j,k) = (Er(i,j,k) - Er(i-1,j,k)) / dx[0] * dc(i,j,k); |
| 851 | } |
| 852 | else if (n == 1) { |
| 853 | dtf(i,j,k) = (Er(i,j,k) - Er(i,j-1,k)) / dx[1] * dc(i,j,k); |
| 854 | } |
| 855 | else { |
| 856 | dtf(i,j,k) = (Er(i,j,k) - Er(i,j,k-1)) / dx[2] * dc(i,j,k); |
| 857 | } |
| 858 | }); |
| 859 | } |
| 860 | } |
| 861 | |
| 862 | // Correct D terms at physical and coarse-fine boundaries. |
| 863 | hem->boundaryDterm(level, &Dterm_face[0], Er, igroup); |
| 864 | |
| 865 | // Correct for metric terms (only has an effect in non-Cartesian geometries). |
no test coverage detected