722HZ7UFINNE3YKSYKP2NHZ5XEG5QQLQHSKC7PREJZR3EX6RDYUAC
#include <AMReX.hxx>
using namespace AMReX;
#include <cctk.h>
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>
#include <AMReX_PlotFileUtil.H>
using namespace amrex;
#include <atomic>
#include <cmath>
using namespace std;
namespace WaveToyAMReX {
// Linear interpolation between (i0, x0) and (i1, x1)
template <typename Y, typename X> Y linterp(Y y0, Y y1, X x0, X x1, X x) {
return Y(x - x0) / Y(x1 - x0) * y0 + Y(x - x1) / Y(x0 - x1) * y1;
}
// Square
template <typename T> T sqr(T x) { return x * x; }
// Standing wave
template <typename T> T standing(T t, T x, T y, T z) {
T Lx = 4, Ly = 4, Lz = 4;
T kx = 2 * M_PI / Lx, ky = 2 * M_PI / Ly, kz = 2 * M_PI / Lz;
T omega = sqrt(sqr(kx) + sqr(ky) + sqr(kz));
return cos(omega * t) * cos(kx * x) * cos(ky * y) * cos(kz * z);
}
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);
}
}
}
extern "C" void WaveToyAMReX_Initialize(CCTK_ARGUMENTS) {
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
CCTK_VINFO("Initialize iteration %d", cctk_iteration);
ghext->time = 0.0;
ghext->delta_time = 0.5 / ghext->ncells;
const CCTK_REAL t0 = ghext->time;
const CCTK_REAL dt = ghext->delta_time;
const CCTK_REAL *restrict const x0 = ghext->geom.ProbLo();
const CCTK_REAL *restrict const x1 = ghext->geom.ProbHi();
// Initialize phi
#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);
CCTK_REAL *restrict const phi = vars.ptr(0, 0, 0, 0);
CCTK_REAL *restrict const phi_p = vars.ptr(0, 0, 0, 1);
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;
CCTK_REAL x = linterp(x0[0], x1[0], -1, 2 * ghext->ncells - 1, 2 * i);
CCTK_REAL y = linterp(x0[1], x1[1], -1, 2 * ghext->ncells - 1, 2 * j);
CCTK_REAL z = linterp(x0[2], x1[2], -1, 2 * ghext->ncells - 1, 2 * k);
phi[idx] = standing(t0, x, y, z);
phi_p[idx] = standing(t0 - dt, x, y, z);
}
}
}
extern "C" void WaveToyAMReX_Evolve(CCTK_ARGUMENTS) {
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
CCTK_VINFO("Evolve iteration %d", cctk_iteration);
// Cycle time levels
MultiFab::Copy(ghext->mfab, ghext->mfab, 1, 2, 1, ghext->nghostzones);
MultiFab::Copy(ghext->mfab, ghext->mfab, 0, 1, 1, ghext->nghostzones);
const CCTK_REAL dt = ghext->delta_time;
const CCTK_REAL *restrict const dx = ghext->geom.CellSize();
// Evolve phi
#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.tilebox();
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);
const CCTK_REAL *restrict const phi_p = vars.ptr(0, 0, 0, 1);
const CCTK_REAL *restrict const phi_p_p = vars.ptr(0, 0, 0, 2);
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;
CCTK_REAL ddx_phi =
(phi_p[idx - di] - 2 * phi_p[idx] + phi_p[idx + di]) / sqr(dx[0]);
CCTK_REAL ddy_phi =
(phi_p[idx - dj] - 2 * phi_p[idx] + phi_p[idx + dj]) / sqr(dx[1]);
CCTK_REAL ddz_phi =
(phi_p[idx - dk] - 2 * phi_p[idx] + phi_p[idx + dk]) / sqr(dx[2]);
phi[idx] = -phi_p_p[idx] + 2 * phi_p[idx] +
sqr(dt) * (ddx_phi + ddy_phi + ddz_phi);
}
}
// Synchronize
ghext->mfab.FillBoundary(ghext->geom.periodicity());
}
extern "C" void WaveToyAMReX_Error(CCTK_ARGUMENTS) {
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
CCTK_VINFO("Error iteration %d", cctk_iteration);
const CCTK_REAL t0 = ghext->time;
const CCTK_REAL *restrict const x0 = ghext->geom.ProbLo();
const CCTK_REAL *restrict const x1 = ghext->geom.ProbHi();
// Initialize phi
#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);
CCTK_REAL *restrict const err = vars.ptr(0, 0, 0, 3);
const 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)
#pragma omp simd
for (int i = imin.x; i <= imax.x; ++i) {
const int idx = i * di + j * dj + k * dk;
CCTK_REAL x = linterp(x0[0], x1[0], -1, 2 * ghext->ncells - 1, 2 * i);
CCTK_REAL y = linterp(x0[1], x1[1], -1, 2 * ghext->ncells - 1, 2 * j);
CCTK_REAL z = linterp(x0[2], x1[2], -1, 2 * ghext->ncells - 1, 2 * k);
err[idx] = phi[idx] - standing(t0, x, y, z);
}
}
}
extern "C" void WaveToyAMReX_Output(CCTK_ARGUMENTS) {
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
CCTK_VINFO("Output iteration %d", cctk_iteration);
// Output phi
string filename = amrex::Concatenate("wavetoy/phi", cctk_iteration, 6);
WriteSingleLevelPlotfile(filename, ghext->mfab,
{"phi", "phi_p", "phi_p_p", "error"}, ghext->geom,
ghext->time, cctk_iteration);
}
} // namespace WaveToyAMReX
# Main make.code.defn file for thorn WaveToyAMReX
# Source files in this directory
SRCS = wavetoy.cxx
# Subdirectories containing source files
SUBDIRS =
# Schedule definitions for thorn WaveToyAMReX
SCHEDULE WaveToyAMReX_Check AT basegrid
{
LANG: C
} "Check looping constructs"
SCHEDULE WaveToyAMReX_Initialize AT initial
{
LANG: C
} "Set up initial conditions for the wave equation"
SCHEDULE WaveToyAMReX_Evolve AT evol
{
LANG: C
} "Evolve the wave equation"
SCHEDULE WaveToyAMReX_Error AT analysis BEFORE WaveToyAMReX_Output
{
LANG: C
} "Calculate error for the wave equation"
SCHEDULE WaveToyAMReX_Output AT analysis
{
LANG: C
} "Output wave equation state"
# Parameter definitions for thorn WaveToyAMReX
# Interface definition for thorn WaveToyAMReX
IMPLEMENTS: WaveToyAMReX
USES INCLUDE HEADER: AMReX.hxx
% *======================================================================*
% 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}
# Configuration definitions for thorn WaveToyAMReX
REQUIRES AMReX
Cactus Code Thorn WaveToyAMReX
Author(s) : Erik Schnetter <schnetter@gmail.com>
Maintainer(s): Erik Schnetter <schnetter@gmail.com>
Licence : LGPL
--------------------------------------------------------------------------
1. Purpose
A wave toy using the AMReX driver
# Main make.code.defn file for thorn AMReX
# Source files in this directory
SRCS = AMReX.cxx
# Subdirectories containing source files
SUBDIRS =
#ifndef AMREX_HXX
#define AMREX_HXX
#include <AMReX.H>
#include <AMReX_BoxArray.H>
#include <AMReX_DistributionMapping.H>
#include <AMReX_Geometry.H>
#include <AMReX_MultiFab.H>
using namespace amrex;
#include <cctk.h>
#include <memory>
namespace AMReX {
using namespace std;
static_assert(is_same<Real, CCTK_REAL>::value,
"AMReX's Real type must be the same as Cactus's CCTK_REAL");
struct GHExt {
const int ncells = 100; // grid size
const int nghostzones = 1; // number of ghost zones
CCTK_REAL time, delta_time;
BoxArray ba;
Geometry geom;
MultiFab mfab;
};
extern unique_ptr<GHExt> ghext;
} // namespace AMReX
#endif // #ifndef AMREX_HXX
#include <AMReX.hxx>
#include <AMReX.H>
#include <cctk.h>
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>
#include <mpi.h>
#include <string>
#include <type_traits>
#include <utility>
namespace AMReX {
using namespace std;
amrex::AMReX *pamrex = nullptr;
unique_ptr<GHExt> ghext;
extern "C" int AMReX_Startup() {
string banner = "AMR driver provided by AMReX " + amrex::Version();
CCTK_RegisterBanner(banner.c_str());
// Initialize AMReX
pamrex = amrex::Initialize(MPI_COMM_WORLD);
// Create grid structure
ghext = make_unique<GHExt>();
// Define box array
IntVect dom_lo(AMREX_D_DECL(0, 0, 0));
IntVect dom_hi(
AMREX_D_DECL(ghext->ncells - 1, ghext->ncells - 1, ghext->ncells - 1));
Box domain(dom_lo, dom_hi);
ghext->ba.define(domain);
// Break up box array into chunks no larger than max_grid_size along
// each direction
const int max_grid_size = 32;
ghext->ba.maxSize(max_grid_size);
// Define physical box
RealBox real_box({AMREX_D_DECL(-1.0, -1.0, -1.0)},
{AMREX_D_DECL(1.0, 1.0, 1.0)});
// Define geometry
Vector<int> is_periodic(AMREX_SPACEDIM, 1); // periodic in all directions
ghext->geom.define(domain, &real_box, CoordSys::cartesian,
is_periodic.data());
const int nvars = 4; // number of grid functions
// Distributed boxes
DistributionMapping dm(ghext->ba);
// Allocate grid hierarchy
ghext->mfab = MultiFab(ghext->ba, dm, nvars, ghext->nghostzones);
return 0;
}
extern "C" int AMReX_Shutdown() {
// Deallocate grid hierarchy
ghext = nullptr;
// Finalize AMReX
amrex::Finalize(pamrex);
pamrex = nullptr;
return 0;
}
} // namespace AMReX
# Schedule definitions for thorn AMReX
SCHEDULE AMReX_Startup AT startup AS Driver_Startup
{
LANG: C
} "Start up the driver"
SCHEDULE AMReX_Shutdown AT shutdown AS zzz_AMReX_Shutdown
{
LANG: C
} "Shut down the driver"
# Parameter definitions for thorn AMReX
# Interface definition for thorn AMReX
IMPLEMENTS: AMReX
INCLUDES HEADER: AMReX.hxx IN AMReX.hxx
#! /bin/bash
################################################################################
# Prepare
################################################################################
# Set up shell
if [ "$(echo ${VERBOSE} | tr '[:upper:]' '[:lower:]')" = 'yes' ]; then
set -x # Output commands
fi
set -e # Abort on errors
################################################################################
# Configure Cactus
################################################################################
if [ -z "${AMREX_DIR}" ]; then
echo "BEGIN ERROR"
echo "Configuration variable AMREX_DIR is not set"
echo "END ERROR"
exit 1
fi
# Set options
: ${AMREX_INC_DIRS="${AMREX_DIR}/include"}
: ${AMREX_LIB_DIRS="${AMREX_DIR}/lib"}
: ${AMREX_LIBS='amrex'}
AMREX_INC_DIRS="$(${CCTK_HOME}/lib/sbin/strip-incdirs.sh ${AMREX_INC_DIRS})"
AMREX_LIB_DIRS="$(${CCTK_HOME}/lib/sbin/strip-libdirs.sh ${AMREX_LIB_DIRS})"
# Pass options to Cactus
echo "BEGIN MAKE_DEFINITION"
echo "AMREX_DIR = ${AMREX_DIR}"
echo "AMREX_INC_DIRS = ${AMREX_INC_DIRS}"
echo "AMREX_LIB_DIRS = ${AMREX_LIB_DIRS}"
echo "AMREX_LIBS = ${AMREX_LIBS}"
echo "END MAKE_DEFINITION"
echo 'INCLUDE_DIRECTORY $(AMREX_INC_DIRS)'
echo 'LIBRARY_DIRECTORY $(AMREX_LIB_DIRS)'
echo 'LIBRARY $(AMREX_LIBS)'
# Configuration definitions for thorn AMReX
PROVIDES AMReX
{
SCRIPT configure.sh
LANG bash
OPTIONS AMREX_DIR
}
REQUIRES AMReX MPI
Cactus Code Thorn AMReX
Author(s) : Erik Schnetter <schnetter@gmail.com>
Maintainer(s): Erik Schnetter <schnetter@gmail.com>
Licence : LGPL
--------------------------------------------------------------------------
1. Purpose
Provide a Cactus driver based on AMReX; see
<https://amrex-codes.github.io>.
From the AMReX web site:
A Massively Block-Structured AMR Software Framework and Applications.
2. Building AMReX library
How to build the AMReX library so that this thorn can use it:
cd $HOME/src/amrex
mkdir -p build
cd build
rm -rf *
cmake -DCMAKE_C_COMPILER=/opt/local/bin/gcc -DCMAKE_CXX_COMPILER=/opt/local/bin/g++ -DCMAKE_Fortran_COMPILER=/opt/local/bin/gfortran -DENABLE_OMP=ON -DCMAKE_CXX_FLAGS='-march=native' -DCMAKE_Fortran_FLAGS='-march=native' -DCMAKE_INSTALL_PREFIX=$HOME/amrex ..
make -j12
make -j12 install
3. Visualizing output