Cactus Code Thorn AMReXTest
Author(s)    : Erik Schnetter <schnetter@gmail.com>
Maintainer(s): Erik Schnetter <schnetter@gmail.com>
Licence      : LGPL
--------------------------------------------------------------------------
1. Purpose
Various self-tests of the AMReX driver

# Configuration definitions for thorn AMReXTest
REQUIRES AMReX

% *======================================================================*
%  Cactus Thorn template for ThornGuide documentation
%  Author: Ian Kelley
%  Date: Sun Jun 02, 2002
%  $Header$
%
%  Thorn documentation in the latex file doc/documentation.tex
%  will be included in ThornGuides built with the Cactus make system.
%  The scripts employed by the make system automatically include
%  pages about variables, parameters and scheduling parsed from the
%  relevant thorn CCL files.
%
%  This template contains guidelines which help to assure that your
%  documentation will be correctly added to ThornGuides. More
%  information is available in the Cactus UsersGuide.
%
%  Guidelines:
%   - Do not change anything before the line
%       % START CACTUS THORNGUIDE",
%     except for filling in the title, author, date, etc. fields.
%        - Each of these fields should only be on ONE line.
%        - Author names should be separated with a \\ or a comma.
%   - You can define your own macros, but they must appear after
%     the START CACTUS THORNGUIDE line, and must not redefine standard
%     latex commands.
%   - To avoid name clashes with other thorns, 'labels', 'citations',
%     'references', and 'image' names should conform to the following
%     convention:
%       ARRANGEMENT_THORN_LABEL
%     For example, an image wave.eps in the arrangement CactusWave and
%     thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps
%   - Graphics should only be included using the graphicx package.
%     More specifically, with the "\includegraphics" command.  Do
%     not specify any graphic file extensions in your .tex file. This
%     will allow us to create a PDF version of the ThornGuide
%     via pdflatex.
%   - References should be included with the latex "\bibitem" command.
%   - Use \begin{abstract}...\end{abstract} instead of \abstract{...}
%   - Do not use \appendix, instead include any appendices you need as
%     standard sections.
%   - For the benefit of our Perl scripts, and for future extensions,
%     please use simple latex.
%
% *======================================================================*
%
% Example of including a graphic image:
%    \begin{figure}[ht]
% 	\begin{center}
%    	   \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure}
% 	\end{center}
% 	\caption{Illustration of this and that}
% 	\label{MyArrangement_MyThorn_MyLabel}
%    \end{figure}
%
% Example of using a label:
%   \label{MyArrangement_MyThorn_MyLabel}
%
% Example of a citation:
%    \cite{MyArrangement_MyThorn_Author99}
%
% Example of including a reference
%   \bibitem{MyArrangement_MyThorn_Author99}
%   {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999),
%   1--16. {\tt http://www.nowhere.com/}}
%
% *======================================================================*
% If you are using CVS use this line to give version information
% $Header$
\documentclass{article}
% Use the Cactus ThornGuide style file
% (Automatically used from Cactus distribution, if you have a
%  thorn without the Cactus Flesh download this from the Cactus
%  homepage at www.cactuscode.org)
\usepackage{../../../../doc/latex/cactus}
\begin{document}
% The author of the documentation
\author{Erik Schnetter \textless schnetter@gmail.com\textgreater}
% The title of the document (not necessarily the name of the Thorn)
\title{WaveToyAMReX}
% the date your document was last changed, if your document is in CVS,
% please use:
%    \date{$ $Date$ $}
\date{July 10 2019}
\maketitle
% Do not delete next line
% START CACTUS THORNGUIDE
% Add all definitions used in this documentation here
%   \def\mydef etc
% Add an abstract for this thorn's documentation
\begin{abstract}
\end{abstract}
% The following sections are suggestive only.
% Remove them or add your own.
\section{Introduction}
\section{Physical System}
\section{Numerical Implementation}
\section{Using This Thorn}
\subsection{Obtaining This Thorn}
\subsection{Basic Usage}
\subsection{Special Behaviour}
\subsection{Interaction With Other Thorns}
\subsection{Examples}
\subsection{Support and Feedback}
\section{History}
\subsection{Thorn Source Code}
\subsection{Thorn Documentation}
\subsection{Acknowledgements}
\begin{thebibliography}{9}
\end{thebibliography}
% Do not delete next line
% END CACTUS THORNGUIDE
\end{document}

# Interface definition for thorn AMReXTest
IMPLEMENTS: AMReXTest
USES INCLUDE HEADER: AMReX.hxx

# Parameter definitions for thorn AMReXTest

# Schedule definitions for thorn AMReXTest
SCHEDULE AMReXTest_CheckLoops AT basegrid
{
  LANG: C
} "Check looping constructs"

#include <AMReX.hxx>
using namespace AMReX;
#include <cctk.h>
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>
#include <atomic>
using namespace std;
namespace AMReXTest {
extern "C" void AMReXTest_CheckLoops(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  CCTK_VINFO("Check loop iterators");
  // Set all grid points to zero: Loop over the interior, then
  // synchronize
  {
    int count = 0;
    for (MFIter mfi(ghext->mfab); mfi.isValid(); ++mfi) {
      const Box &fbx = mfi.fabbox();
      const Box &bx = mfi.validbox();
      const Dim3 amin = lbound(fbx);
      const Dim3 amax = ubound(fbx);
      const Dim3 imin = lbound(bx);
      const Dim3 imax = ubound(bx);
      constexpr int di = 1;
      const int dj = di * (amax.x - amin.x + 1);
      const int dk = dj * (amax.y - amin.y + 1);
      const Array4<CCTK_REAL> &vars = ghext->mfab.array(mfi);
      CCTK_REAL *restrict const phi = vars.ptr(0, 0, 0, 0);
      for (int k = imin.z; k <= imax.z; ++k)
        for (int j = imin.y; j <= imax.y; ++j)
          for (int i = imin.x; i <= imax.x; ++i) {
            const int idx = i * di + j * dj + k * dk;
            phi[idx] = 0;
            ++count;
          }
    }
    ParallelDescriptor::ReduceIntSum(count);
    assert(count == ghext->ncells * ghext->ncells * ghext->ncells);
  }
  ghext->mfab.FillBoundary(ghext->geom.periodicity());
  // Increase each grid point (both interior and ghost zones) by one
  // using atomic operations. This will catch cases where we either
  // omit points or traverse points twice.
#pragma omp parallel
  for (MFIter mfi(ghext->mfab,
                  MFItInfo().SetDynamic(true).EnableTiling({1024000, 16, 32}));
       mfi.isValid(); ++mfi) {
    const Box &fbx = mfi.fabbox();
    const Box &bx = mfi.growntilebox();
    const Dim3 amin = lbound(fbx);
    const Dim3 amax = ubound(fbx);
    const Dim3 imin = lbound(bx);
    const Dim3 imax = ubound(bx);
    constexpr int di = 1;
    const int dj = di * (amax.x - amin.x + 1);
    const int dk = dj * (amax.y - amin.y + 1);
    const Array4<CCTK_REAL> &vars = ghext->mfab.array(mfi);
    atomic<CCTK_REAL> *restrict const phi =
        (atomic<CCTK_REAL> *)vars.ptr(0, 0, 0, 0);
    for (int k = imin.z; k <= imax.z; ++k)
      for (int j = imin.y; j <= imax.y; ++j)
#pragma omp simd
        for (int i = imin.x; i <= imax.x; ++i) {
          const int idx = i * di + j * dj + k * dk;
          // phi[idx] += 1.0;
          CCTK_REAL expected = 0.0;
          bool success = phi[idx].compare_exchange_strong(expected, 1.0);
          // assert(success);
        }
  }
  // Check all grid points whether they are set to one
  for (MFIter mfi(ghext->mfab); mfi.isValid(); ++mfi) {
    const Box &fbx = mfi.fabbox();
    const Box &bx = mfi.fabbox();
    const Dim3 amin = lbound(fbx);
    const Dim3 amax = ubound(fbx);
    const Dim3 imin = lbound(bx);
    const Dim3 imax = ubound(bx);
    constexpr int di = 1;
    const int dj = di * (amax.x - amin.x + 1);
    const int dk = dj * (amax.y - amin.y + 1);
    const Array4<CCTK_REAL> &vars = ghext->mfab.array(mfi);
    CCTK_REAL *restrict const phi = vars.ptr(0, 0, 0, 0);
    for (int k = imin.z; k <= imax.z; ++k)
      for (int j = imin.y; j <= imax.y; ++j)
        for (int i = imin.x; i <= imax.x; ++i) {
          const int idx = i * di + j * dj + k * dk;
          assert(phi[idx] == 1);
        }
  }
}
} // namespace AMReXTest

# Main make.code.defn file for thorn AMReXTest
# Source files in this directory
SRCS = check_loops.cxx
# Subdirectories containing source files
SUBDIRS = 

SCHEDULE WaveToyAMReX_Check AT basegrid
{
  LANG: C
} "Check looping constructs"

#include <atomic>

}
extern "C" void WaveToyAMReX_Check(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  CCTK_VINFO("Check iterators");
  {
    int count = 0;
    for (MFIter mfi(ghext->mfab); mfi.isValid(); ++mfi) {
      const Box &fbx = mfi.fabbox();
      const Box &bx = mfi.validbox();
      const Dim3 amin = lbound(fbx);
      const Dim3 amax = ubound(fbx);
      const Dim3 imin = lbound(bx);
      const Dim3 imax = ubound(bx);
      constexpr int di = 1;
      const int dj = di * (amax.x - amin.x + 1);
      const int dk = dj * (amax.y - amin.y + 1);
      const Array4<CCTK_REAL> &vars = ghext->mfab.array(mfi);
      CCTK_REAL *restrict const phi = vars.ptr(0, 0, 0, 0);
      for (int k = imin.z; k <= imax.z; ++k)
        for (int j = imin.y; j <= imax.y; ++j)
          for (int i = imin.x; i <= imax.x; ++i) {
            const int idx = i * di + j * dj + k * dk;
            phi[idx] = 0;
            ++count;
          }
    }
    ParallelDescriptor::ReduceIntSum(count);
    assert(count == ghext->ncells * ghext->ncells * ghext->ncells);
  }
  ghext->mfab.FillBoundary(ghext->geom.periodicity());
#pragma omp parallel
  for (MFIter mfi(ghext->mfab,
                  MFItInfo().SetDynamic(true).EnableTiling({1024000, 16, 32}));
       mfi.isValid(); ++mfi) {
    const Box &fbx = mfi.fabbox();
    const Box &bx = mfi.growntilebox();
    const Dim3 amin = lbound(fbx);
    const Dim3 amax = ubound(fbx);
    const Dim3 imin = lbound(bx);
    const Dim3 imax = ubound(bx);
    constexpr int di = 1;
    const int dj = di * (amax.x - amin.x + 1);
    const int dk = dj * (amax.y - amin.y + 1);
    const Array4<CCTK_REAL> &vars = ghext->mfab.array(mfi);
    atomic<CCTK_REAL> *restrict const phi =
        (atomic<CCTK_REAL> *)vars.ptr(0, 0, 0, 0);
    for (int k = imin.z; k <= imax.z; ++k)
      for (int j = imin.y; j <= imax.y; ++j)
#pragma omp simd
        for (int i = imin.x; i <= imax.x; ++i) {
          const int idx = i * di + j * dj + k * dk;
          // phi[idx] += 1.0;
          CCTK_REAL expected = 0.0;
          bool success = phi[idx].compare_exchange_strong(expected, 1.0);
          // assert(success);
        }
  }
  for (MFIter mfi(ghext->mfab); mfi.isValid(); ++mfi) {
    const Box &fbx = mfi.fabbox();
    const Box &bx = mfi.fabbox();
    const Dim3 amin = lbound(fbx);
    const Dim3 amax = ubound(fbx);
    const Dim3 imin = lbound(bx);
    const Dim3 imax = ubound(bx);
    constexpr int di = 1;
    const int dj = di * (amax.x - amin.x + 1);
    const int dk = dj * (amax.y - amin.y + 1);
    const Array4<CCTK_REAL> &vars = ghext->mfab.array(mfi);
    CCTK_REAL *restrict const phi = vars.ptr(0, 0, 0, 0);
    for (int k = imin.z; k <= imax.z; ++k)
      for (int j = imin.y; j <= imax.y; ++j)
        for (int i = imin.x; i <= imax.x; ++i) {
          const int idx = i * di + j * dj + k * dk;
          assert(phi[idx] == 1);
        }
  }