From 836af57b85364cac331a6879f81ca855065c15ba Mon Sep 17 00:00:00 2001 From: schnetter Date: Fri, 4 Sep 2009 18:19:51 +0000 Subject: Add thorn implementing the NewRad boundary condition (known from BSSN_MoL) git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/NewRad/trunk@2 16f5cafb-89ad-4b9f-9d0b-9d7a0a37d03b --- README | 10 ++ configuration.ccl | 3 + doc/documentation.tex | 144 +++++++++++++++++++++++ interface.ccl | 16 +++ param.ccl | 1 + schedule.ccl | 1 + src/make.code.defn | 7 ++ src/newrad.cc | 311 ++++++++++++++++++++++++++++++++++++++++++++++++++ 8 files changed, 493 insertions(+) create mode 100644 README create mode 100644 configuration.ccl create mode 100644 doc/documentation.tex create mode 100644 interface.ccl create mode 100644 param.ccl create mode 100644 schedule.ccl create mode 100644 src/make.code.defn create mode 100644 src/newrad.cc diff --git a/README b/README new file mode 100644 index 0000000..b6aebd8 --- /dev/null +++ b/README @@ -0,0 +1,10 @@ +Cactus Code Thorn NewRad +Author(s) : Erik Schnetter +Maintainer(s): Erik Schnetter +Licence : GPL +-------------------------------------------------------------------------- + +1. Purpose + +Implement the "NewRad" radiative boundary conditions known from +BSSN_MoL. diff --git a/configuration.ccl b/configuration.ccl new file mode 100644 index 0000000..b8c5c16 --- /dev/null +++ b/configuration.ccl @@ -0,0 +1,3 @@ +# Configuration definition for thorn NewRad + +REQUIRES THORNS: GenericFD diff --git a/doc/documentation.tex b/doc/documentation.tex new file mode 100644 index 0000000..089de90 --- /dev/null +++ b/doc/documentation.tex @@ -0,0 +1,144 @@ +% *======================================================================* +% 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@cct.lsu.edu{\textgreater} \textless n\textgreater} + +% The title of the document (not necessarily the name of the Thorn) +\title{NewRad} + +% the date your document was last changed, if your document is in CVS, +% please use: +% \date{$ $Date$ $} +\date{September 03 2009} + +\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} diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..57f6416 --- /dev/null +++ b/interface.ccl @@ -0,0 +1,16 @@ +# Interface definition for thorn NewRad + +IMPLEMENTS: NewRad + +USES INCLUDE HEADER: GenericFD.h + +CCTK_INT FUNCTION \ + NewRad_Apply \ + (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_REAL ARRAY IN var, \ + CCTK_REAL ARRAY INOUT rhs, \ + CCTK_REAL IN var0, \ + CCTK_REAL IN v0, \ + CCTK_REAL IN radpower, \ + CCTK_INT IN width) +PROVIDES FUNCTION NewRad_Apply WITH NewRad_Apply1 LANGUAGE C diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..46bf4ed --- /dev/null +++ b/param.ccl @@ -0,0 +1 @@ +# Parameter definitions for thorn NewRad diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..0814a06 --- /dev/null +++ b/schedule.ccl @@ -0,0 +1 @@ +# Schedule definitions for thorn NewRad diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..244b721 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,7 @@ +# Main make.code.defn file for thorn NewRad + +# Source files in this directory +SRCS = newrad.cc + +# Subdirectories containing source files +SUBDIRS = diff --git a/src/newrad.cc b/src/newrad.cc new file mode 100644 index 0000000..a1edcf0 --- /dev/null +++ b/src/newrad.cc @@ -0,0 +1,311 @@ +#include + +#include + +#ifdef CCTK_CXX_RESTRICT +# undef restrict +# define restrict CCTK_CXX_RESTRICT +#endif + +#define KRANC_C +extern "C" { +#include +} + + + +// Adapted from BSSN_MoL's files NewRad.F and newrad.h +static +void newrad_kernel (cGH const* restrict const cctkGH, + int const* restrict const bmin, + int const* restrict const bmax, + int const* restrict const dir, + CCTK_REAL const* restrict const var, + CCTK_REAL * restrict const rhs, + CCTK_REAL const* restrict const x, + CCTK_REAL const* restrict const y, + CCTK_REAL const* restrict const z, + CCTK_REAL const* restrict const r, + CCTK_REAL const& var0, + CCTK_REAL const& v0, + CCTK_REAL const& radpower) +{ + int const ni = cctkGH->cctk_lsh[0]; + int const nj = cctkGH->cctk_lsh[1]; + int const nk = cctkGH->cctk_lsh[2]; + + int const si = -dir[0]; + int const sj = -dir[1]; + int const sk = -dir[2]; + + int const di = 1; + int const dj = ni; + int const dk = ni*nj; + + CCTK_REAL const dx = cctkGH->cctk_delta_space[0] / cctkGH->cctk_levfac[0]; + CCTK_REAL const dy = cctkGH->cctk_delta_space[1] / cctkGH->cctk_levfac[1]; + CCTK_REAL const dz = cctkGH->cctk_delta_space[2] / cctkGH->cctk_levfac[2]; + CCTK_REAL const idx = 1.0/dx; + CCTK_REAL const idy = 1.0/dy; + CCTK_REAL const idz = 1.0/dz; + + for (int k=bmin[2]; k0 ? +1 : -1; + int const svj = j==0 ? +1 : j==nj-1 ? -1 : vy>0 ? +1 : -1; + int const svk = k==0 ? +1 : k==nk-1 ? -1 : vz>0 ? +1 : -1; + + // Find x derivative + CCTK_REAL derivx; + if (i>0 and i0 and j0 and k= 0) { + // ***************************************** + // *** EXTRAPOLATION OF MISSING PART *** + // ***************************************** + // + // Here we try to extrapolate for the part of the boundary + // that does not behave as a pure wave (i.e. Coulomb type + // terms caused by infall of the coordinate lines). + // + // This we do by comparing the source term one grid point + // away from the boundary (which we already have), to what + // we would have obtained if we had used the boundary + // condition there. The difference gives us an idea of the + // missing part and we extrapolate that to the boundary + // assuming a power-law decay. + + int const ip = i-si; + int const jp = j-sj; + int const kp = k-sk; + + // Find local wave speeds + int const indp = CCTK_GFINDEX3D(cctkGH, ip,jp,kp); + + CCTK_REAL const rp = r[indp]; + CCTK_REAL const rpi = 1.0/rp; + + CCTK_REAL const vx = v0*x[indp]*rpi; + CCTK_REAL const vy = v0*y[indp]*rpi; + CCTK_REAL const vz = v0*z[indp]*rpi; + + int const svi = i==0 ? +1 : i==ni-1 ? -1 : vx>0 ? +1 : -1; + int const svj = j==0 ? +1 : j==nj-1 ? -1 : vy>0 ? +1 : -1; + int const svk = k==0 ? +1 : k==nk-1 ? -1 : vz>0 ? +1 : -1; + + // Find x derivative + CCTK_REAL derivx; + if (ip>0 and ip0 and jp0 and kp=0 + + } // for i j k + } + } +} + + + +// Adapted from Kranc's KrancNumericalTools/GenericFD's file +// GenericFD.c +static +void newrad_loop (cGH const* restrict const cctkGH, + CCTK_REAL const* restrict const var, + CCTK_REAL * restrict const rhs, + CCTK_REAL const* restrict const x, + CCTK_REAL const* restrict const y, + CCTK_REAL const* restrict const z, + CCTK_REAL const* restrict const r, + CCTK_REAL const& var0, + CCTK_REAL const& v0, + CCTK_REAL const& radpower, + int const width) +{ + int imin[3], imax[3], is_symbnd[6], is_physbnd[6], is_ipbnd[6]; + GenericFD_GetBoundaryInfo + (cctkGH, cctkGH->cctk_lsh, cctkGH->cctk_lssh, cctkGH->cctk_bbox, + cctkGH->cctk_nghostzones, + imin, imax, is_symbnd, is_physbnd, is_ipbnd); + + // Loop over all faces + for (int dir2=-1; dir2<=+1; ++dir2) { + for (int dir1=-1; dir1<=+1; ++dir1) { + for (int dir0=-1; dir0<=+1; ++dir0) { + int const dir[3] = { dir0, dir1, dir2 }; + + // one of tahe faces is a boundary + bool have_bnd = false; + // all boundary faces are physical boundaries + bool all_physbnd = true; + + int bmin[3], bmax[3]; + for (int d=0; d<3; ++d) { + switch (dir[d]) { + case -1: + bmin[d] = 0; + bmax[d] = imin[d]; + have_bnd = true; + all_physbnd = all_physbnd and is_physbnd[2*d+0]; + break; + case 0: + bmin[d] = imin[d]; + bmax[d] = imax[d]; + break; + case +1: + bmin[d] = imax[d]; + bmax[d] = cctkGH->cctk_lssh[CCTK_LSSH_IDX(0,d)]; + have_bnd = true; + all_physbnd = all_physbnd and is_physbnd[2*d+1]; + break; + } + } + + if (have_bnd and all_physbnd) { + newrad_kernel (cctkGH, bmin, bmax, dir, + var, rhs, x,y,z,r, var0, v0, radpower); + } + + } // for dir0 dir1 dir2 + } + } +} + + + +extern "C" +CCTK_INT NewRad_Apply1 (CCTK_POINTER_TO_CONST const cctkGH_, + CCTK_REAL const* restrict const var, + CCTK_REAL * restrict const rhs, + CCTK_REAL const var0, + CCTK_REAL const v0, + CCTK_REAL const radpower, + CCTK_INT const width) +{ + cGH const* restrict const cctkGH = static_cast (cctkGH_); + if (not cctkGH) { + CCTK_WARN (CCTK_WARN_ABORT, + "cctkGH is NULL"); + } + +#if 0 + CCTK_REAL const* restrict const var = CCTK_VarDataPtr (cctkGH, 0, varname); + if (not var) { + CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot access variable \"%s\"", varname); + } + CCTK_REAL * restrict const rhs = CCTK_VarDataPtr (cctkGH, 0, rhsname); + if (not rhs) { + CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot access RHS variable \"%s\"", rhsname); + } +#endif + + if (not var) { + CCTK_WARN (CCTK_WARN_ABORT, + "Pointer to variable is NULL"); + } + if (not rhs) { + CCTK_WARN (CCTK_WARN_ABORT, + "Pointer to RHS is NULL"); + } + + CCTK_REAL const* restrict const x = + static_cast (CCTK_VarDataPtr (cctkGH, 0, "grid::x")); + CCTK_REAL const* restrict const y = + static_cast (CCTK_VarDataPtr (cctkGH, 0, "grid::y")); + CCTK_REAL const* restrict const z = + static_cast (CCTK_VarDataPtr (cctkGH, 0, "grid::z")); + CCTK_REAL const* restrict const r = + static_cast (CCTK_VarDataPtr (cctkGH, 0, "grid::r")); + if (not x or not y or not z or not z) { + CCTK_WARN (CCTK_WARN_ABORT, + "Cannot access coordinate variables x, y, z, and r"); + } + + newrad_loop (cctkGH, var, rhs, x,y,z,r, var0, v0, radpower, width); + + return 0; +} -- cgit v1.2.3