Method gwstrain

Source/driver/sum_utils.cpp:317–605

#ifdef GRAVITY
void
Castro::gwstrain (Real time,
                  Real& h_plus_1, Real& h_cross_1,
                  Real& h_plus_2, Real& h_cross_2,
                  Real& h_plus_3, Real& h_cross_3,
                  bool local) {

    amrex::ignore_unused(time);

    BL_PROFILE("Castro::gwstrain()");

    // We have nothing to do if the user did not request the gravitational wave
    // strain (inferred from whether the observation distance is positive).
    if (castro::gw_dist <= 0.0_rt) {
        return;
    }

    constexpr amrex::Real c_light_4 = amrex::Math::powi<4>(C::c_light);

    GeometryData geomdata = geom.data();

    MultiFab& S_new = get_new_data(State_Type);
    MultiFab& grav_new = get_new_data(Gravity_Type);

    bool mask_available = level < parent->finestLevel();

    MultiFab tmp_mf;
    const MultiFab& mask_mf = mask_available ? getLevel(level+1).build_fine_mask() : tmp_mf;

    // Qtt stores the second time derivative of the quadrupole moment.
    // We calculate it directly rather than computing the quadrupole moment
    // and differentiating it in time, because the latter method is less accurate
    // and requires the state at other timesteps. See, e.g., Equation 5 of
    // Loren-Aguilar et al. 2005.

    ReduceOps<ReduceOpSum, ReduceOpSum, ReduceOpSum,
              ReduceOpSum, ReduceOpSum, ReduceOpSum,
              ReduceOpSum, ReduceOpSum, ReduceOpSum> reduce_op;
    ReduceData<Real, Real, Real,
               Real, Real, Real,
               Real, Real, Real> reduce_data(reduce_op);
    using ReduceTuple = typename decltype(reduce_data)::Type;

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

        const Box& box = mfi.tilebox();

        auto rho = S_new[mfi].array(URHO);
        auto vol = volume.array(mfi);
        auto xmom = S_new[mfi].array(UMX);
        auto ymom = S_new[mfi].array(UMY);
        auto zmom = S_new[mfi].array(UMZ);
        auto gravx = grav_new[mfi].array(0);
        auto gravy = grav_new[mfi].array(1);
        auto gravz = grav_new[mfi].array(2);