CoarsenIO.H

Go to the documentation of this file.

#ifndef WARPX_COARSEN_IO_H_
#define WARPX_COARSEN_IO_H_

#include "WarpX.H"
#include <AMReX_Math.H>

namespace CoarsenIO{

    using namespace amrex;

    AMREX_GPU_DEVICE
    AMREX_FORCE_INLINE
    Real Interp ( Array4<Real const> const& arr_src,
                  GpuArray<int,3> const& sf,
                  GpuArray<int,3> const& sc,
                  GpuArray<int,3> const& cr,
                  const int i,
                  const int j,
                  const int k,
                  const int comp )
    {
        // Indices of destination array (coarse)
        const int ic[3] = { i, j, k };

        // Number of points and starting indices of source array (fine)
        int np[3], idx_min[3];

        // Compute number of points
        for ( int l = 0; l < 3; ++l ) {
            if ( cr[l] == 1 ) np[l] = 1+amrex::Math::abs(sf[l]-sc[l]); // no coarsening
            else              np[l] = 2-sf[l];
        }

        // Compute starting indices of source array (fine)
        for ( int l = 0; l < 3; ++l ) {
            if ( cr[l] == 1 ) idx_min[l] = ic[l]-sc[l]*(1-sf[l]); // no coarsening
            else              idx_min[l] = ic[l]*cr[l]+static_cast<int>(cr[l]/2)*(1-sc[l])-(1-sf[l]);
        }

        // Auxiliary integer variables
        const int numx = np[0];
        const int numy = np[1];
        const int numz = np[2];
        const int imin = idx_min[0];
        const int jmin = idx_min[1];
        const int kmin = idx_min[2];
        int  ii, jj, kk;
        Real wx, wy, wz;

        // Interpolate over points computed above
        Real c = 0.0_rt;
        for         (int kref = 0; kref < numz; ++kref) {
            for     (int jref = 0; jref < numy; ++jref) {
                for (int iref = 0; iref < numx; ++iref) {
                    ii = imin+iref;
                    jj = jmin+jref;
                    kk = kmin+kref;
                    wx = 1.0_rt/static_cast<Real>(numx);
                    wy = 1.0_rt/static_cast<Real>(numy);
                    wz = 1.0_rt/static_cast<Real>(numz);
                    c += wx*wy*wz*arr_src(ii,jj,kk,comp);
                }
            }
        }
        return c;
    }

    void Loop ( MultiFab& mf_dst,
                const MultiFab& mf_src,
                const int dcomp,
                const int scomp,
                const int ncomp,
                const int ngrow,
                const IntVect crse_ratio=IntVect(1) );

    void Coarsen ( MultiFab& mf_dst,
                   const MultiFab& mf_src,
                   const int dcomp,
                   const int scomp,
                   const int ncomp,
                   const int ngrow,
                   const IntVect crse_ratio=IntVect(1) );
}

#endif // WARPX_COARSEN_IO_H_