wavetoy.cc

wavetoy.cc

#include "dgops.hh"

#include <cctk.h>
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>

#include <cassert>

namespace DGOps {

typedef some_dgop<4> dgop;

constexpr ptrdiff_t order = dgop::order;

constexpr ptrdiff_t npoints_i = order + 1;
constexpr ptrdiff_t npoints_j = order + 1;
constexpr ptrdiff_t npoints_k = order + 1;

extern "C" void DGOps_WaveToy(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

  const ptrdiff_t ni = cctk_lsh[0];
  const ptrdiff_t nj = cctk_lsh[1];
  const ptrdiff_t nk = cctk_lsh[2];

  assert((ni - 2) % npoints_i == 0);
  assert((nj - 2) % npoints_j == 0);
  assert((nk - 2) % npoints_k == 0);

  constexpr ptrdiff_t di = sizeof(CCTK_REAL);
  const ptrdiff_t dj = di * cctk_ash[0];
  const ptrdiff_t dk = dj * cctk_ash[1];

  constexpr ptrdiff_t vs = CCTK_REAL_VEC_SIZE;

  constexpr ptrdiff_t tni = alignup(npoints_i, vs);
  constexpr ptrdiff_t tnj = npoints_j;
  constexpr ptrdiff_t tnk = npoints_k;

  constexpr ptrdiff_t tdi = sizeof(CCTK_REAL);
  constexpr ptrdiff_t tdj = tdi * tni;
  constexpr ptrdiff_t tdk = tdj * tnj;
  constexpr ptrdiff_t tsz = tdk * tnk;

#pragma omp parallel
  {
    unsigned char utT[tsz] CCTK_ATTRIBUTE_ALIGNED(sizeof(CCTK_REAL_VEC));

    unsigned char ut_dxT[tsz] CCTK_ATTRIBUTE_ALIGNED(sizeof(CCTK_REAL_VEC));
    unsigned char ut_dyT[tsz] CCTK_ATTRIBUTE_ALIGNED(sizeof(CCTK_REAL_VEC));
    unsigned char ut_dzT[tsz] CCTK_ATTRIBUTE_ALIGNED(sizeof(CCTK_REAL_VEC));
    unsigned char ux_dxT[tsz] CCTK_ATTRIBUTE_ALIGNED(sizeof(CCTK_REAL_VEC));
    unsigned char uy_dyT[tsz] CCTK_ATTRIBUTE_ALIGNED(sizeof(CCTK_REAL_VEC));
    unsigned char uz_dzT[tsz] CCTK_ATTRIBUTE_ALIGNED(sizeof(CCTK_REAL_VEC));

    unsigned char u_rhsT[tsz] CCTK_ATTRIBUTE_ALIGNED(sizeof(CCTK_REAL_VEC));
    unsigned char ut_rhsT[tsz] CCTK_ATTRIBUTE_ALIGNED(sizeof(CCTK_REAL_VEC));
    unsigned char ux_rhsT[tsz] CCTK_ATTRIBUTE_ALIGNED(sizeof(CCTK_REAL_VEC));
    unsigned char uy_rhsT[tsz] CCTK_ATTRIBUTE_ALIGNED(sizeof(CCTK_REAL_VEC));
    unsigned char uz_rhsT[tsz] CCTK_ATTRIBUTE_ALIGNED(sizeof(CCTK_REAL_VEC));

#pragma omp for collapse(3)
    for (ptrdiff_t k1 = 1; k1 < nk - 1; k1 += npoints_k) {
      for (ptrdiff_t j1 = 1; j1 < nj - 1; j1 += npoints_j) {
        for (ptrdiff_t i1 = 1; i1 < ni - 1; i1 += npoints_i) {
          const ptrdiff_t off1 = i1 * di + j1 * dj + k1 * dk;

          // d/dt u = ut
          // d/dt ut = d/dx ux + d/dy uy + d/dz uz
          // d/dt ux = d/dx ut
          // d/dt uy = d/dy ut
          // d/dt uz = d/dz ut

          copy_dg_dim3<dgop>(&((const unsigned char *)ut)[off1], utT, dj, dk);
          stencil_dg_dim3_dir0<dgop>(&((const unsigned char *)ut)[off1], ut_dxT,
                                     dj, dk);
          stencil_dg_dim3_dir1<dgop>(&((const unsigned char *)ut)[off1], ut_dyT,
                                     dj, dk);
          stencil_dg_dim3_dir2<dgop>(&((const unsigned char *)ut)[off1], ut_dzT,
                                     dj, dk);
          stencil_dg_dim3_dir0<dgop>(&((const unsigned char *)ux)[off1], ux_dxT,
                                     dj, dk);
          stencil_dg_dim3_dir1<dgop>(&((const unsigned char *)uy)[off1], uy_dyT,
                                     dj, dk);
          stencil_dg_dim3_dir2<dgop>(&((const unsigned char *)uz)[off1], uz_dzT,
                                     dj, dk);

          for (ptrdiff_t k = 0; k < npoints_k; ++k) {
            for (ptrdiff_t j = 0; j < npoints_j; ++j) {
              for (ptrdiff_t i = 0; i < npoints_i; ++i) {
                const ptrdiff_t off = i * tdi + j * tdj + k * tdk;

                const CCTK_REAL_VEC utL = vec_load(getelt(utT, off));
                const CCTK_REAL_VEC ut_dxL = vec_load(getelt(ut_dxT, off));
                const CCTK_REAL_VEC ut_dyL = vec_load(getelt(ut_dyT, off));
                const CCTK_REAL_VEC ut_dzL = vec_load(getelt(ut_dzT, off));
                const CCTK_REAL_VEC ux_dxL = vec_load(getelt(ux_dxT, off));
                const CCTK_REAL_VEC uy_dyL = vec_load(getelt(uy_dyT, off));
                const CCTK_REAL_VEC uz_dzL = vec_load(getelt(uz_dzT, off));

                const CCTK_REAL_VEC u_rhsL = utL;
                const CCTK_REAL_VEC ut_rhsL =
                    kadd(kadd(ux_dxL, uy_dyL), uz_dzL);
                const CCTK_REAL_VEC ux_rhsL = ut_dxL;
                const CCTK_REAL_VEC uy_rhsL = ut_dyL;
                const CCTK_REAL_VEC uz_rhsL = ut_dzL;

                vec_store(getelt(u_rhsT, off), u_rhsL);
                vec_store(getelt(ut_rhsT, off), ut_rhsL);
                vec_store(getelt(ux_rhsT, off), ux_rhsL);
                vec_store(getelt(uy_rhsT, off), uy_rhsL);
                vec_store(getelt(uz_rhsT, off), uz_rhsL);
              }
            }
          }

          save_dg_dim3<dgop>(&((unsigned char *)ut_rhs)[off1], ut_rhsT, dj, dk);
          save_dg_dim3<dgop>(&((unsigned char *)ux_rhs)[off1], ux_rhsT, dj, dk);
          save_dg_dim3<dgop>(&((unsigned char *)uy_rhs)[off1], uy_rhsT, dj, dk);
          save_dg_dim3<dgop>(&((unsigned char *)uz_rhs)[off1], uz_rhsT, dj, dk);
        }
      }
    }
  }
}
}