| 3 | using namespace amrex; |
| 4 | |
| 5 | void |
| 6 | Castro::pointmass_update(Real time, Real dt) |
| 7 | { |
| 8 | |
| 9 | amrex::ignore_unused(time); |
| 10 | amrex::ignore_unused(dt); |
| 11 | |
| 12 | int finest_level = parent->finestLevel(); |
| 13 | |
| 14 | if (level == finest_level && point_mass_fix_solution) |
| 15 | { |
| 16 | |
| 17 | MultiFab& S_old = get_old_data(State_Type); |
| 18 | MultiFab& S_new = get_new_data(State_Type); |
| 19 | |
| 20 | const auto dx = geom.CellSizeArray(); |
| 21 | const auto problo = geom.ProbLoArray(); |
| 22 | |
| 23 | ReduceOps<ReduceOpSum> reduce_op; |
| 24 | ReduceData<Real> reduce_data(reduce_op); |
| 25 | using ReduceTuple = typename decltype(reduce_data)::Type; |
| 26 | |
| 27 | #ifdef _OPENMP |
| 28 | #pragma omp parallel |
| 29 | #endif |
| 30 | for (MFIter mfi(S_new, TilingIfNotGPU()); mfi.isValid(); ++mfi) { |
| 31 | |
| 32 | const Box& bx = mfi.tilebox(); |
| 33 | |
| 34 | Array4<Real const> const uin = S_old.array(mfi); |
| 35 | Array4<Real const> const uout = S_new.array(mfi); |
| 36 | Array4<Real const> const vol = volume.array(mfi); |
| 37 | |
| 38 | reduce_op.eval(bx, reduce_data, |
| 39 | [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple |
| 40 | { |
| 41 | // This is just a small number to keep precision issues from making |
| 42 | // icen, jcen, kcen one cell too low. |
| 43 | const Real eps = 1.e-8_rt; |
| 44 | |
| 45 | // This should be the cell whose lower left corner is at center |
| 46 | int icen = static_cast<int>(std::floor((problem::center[0] - problo[0]) / dx[0] + eps)); |
| 47 | #if AMREX_SPACEDIM >= 2 |
| 48 | int jcen = static_cast<int>(std::floor((problem::center[1] - problo[1]) / dx[1] + eps)); |
| 49 | #endif |
| 50 | #if AMREX_SPACEDIM == 3 |
| 51 | int kcen = static_cast<int>(std::floor((problem::center[2] - problo[2]) / dx[2] + eps)); |
| 52 | #endif |
| 53 | |
| 54 | // Make sure we only count contributions from this grid |
| 55 | |
| 56 | const int box_size = 2; |
| 57 | |
| 58 | int istart = amrex::max(icen - box_size, bx.smallEnd(0)); |
| 59 | int iend = amrex::min(icen + box_size - 1, bx.bigEnd(0)); |
| 60 | #if AMREX_SPACEDIM >= 2 |
| 61 | int jstart = amrex::max(jcen - box_size, bx.smallEnd(1)); |
| 62 | int jend = amrex::min(jcen + box_size - 1, bx.bigEnd(1)); |
no test coverage detected