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

8 files changed, 493 insertions, 0 deletions

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 <schnetter@cct.lsu.edu>
+Maintainer(s): Erik Schnetter <schnetter@cct.lsu.edu>
+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 <math.h>
+
+#include <cctk.h>
+
+#ifdef CCTK_CXX_RESTRICT
+#  undef restrict
+#  define restrict CCTK_CXX_RESTRICT
+#endif
+
+#define KRANC_C
+extern "C" {
+#include <GenericFD.h>
+}
+
+
+
+// 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]; k<bmax[2]; ++k) {
+    for (int j=bmin[1]; j<bmax[1]; ++j) {
+      for (int i=bmin[0]; i<bmax[0]; ++i) {
+        int const ind = CCTK_GFINDEX3D(cctkGH, i,j,k);
+        
+        {
+          // The main part of the boundary condition assumes that we
+          // have an outgoing radial wave with some speed v0:
+          //
+          //    var  =  var0 + u(r-v0*t)/r
+          //
+          // This implies the following differential equation:
+          //
+          //    d_t var  =  - v^i d_i var  -  v0 (var - var0) / r
+          //
+          // where  vi = v0 xi/r
+          
+          // Find local wave speeds
+          CCTK_REAL const rp = r[ind];
+          CCTK_REAL const rpi = 1.0/rp;
+          
+          CCTK_REAL const vx = v0*x[ind]*rpi;
+          CCTK_REAL const vy = v0*y[ind]*rpi;
+          CCTK_REAL const vz = v0*z[ind]*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 (i>0 and i<ni-1) {
+            derivx = 0.5*(var[ind+di]-var[ind-di])*idx;
+          } else {
+            derivx = svi*0.5*(3*var[ind] - 4*var[ind-svi*di]
+                              + var[ind-2*svi*di])*idx;
+          }
+          
+          // Find y derivative
+          CCTK_REAL derivy;
+          if (j>0 and j<nj-1) {
+            derivy = 0.5*(var[ind+dj]-var[ind-dj])*idy;
+          } else {
+            derivy = svj*0.5*(3*var[ind] - 4*var[ind-svj*dj]
+                              + var[ind-2*svj*dj])*idy;
+          }
+          
+          // Find z derivative
+          CCTK_REAL derivz;
+          if (k>0 and k<nk-1) {
+            derivz = 0.5*(var[ind+dk]-var[ind-dk])*idz;
+          } else {
+            derivz = svk*0.5*(3*var[ind] - 4*var[ind-svk*dk]
+                              + var[ind-2*svk*dk])*idz;
+          }
+          
+          // Calculate source term
+          rhs[ind] =
+            - vx*derivx - vy*derivy - vz*derivz - v0*(var[ind] - var0)*rpi;
+          
+        }
+        
+        if (radpower >= 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 ip<ni-1) {
+            derivx = 0.5*(var[indp+di]-var[indp-di])*idx;
+          } else {
+            derivx = svi*0.5*(3*var[indp] - 4*var[indp-svi*di]
+                              + var[indp-2*svi*di])*idx;
+          }
+          
+          // Find y derivative
+          CCTK_REAL derivy;
+          if (jp>0 and jp<nj-1) {
+            derivy = 0.5*(var[indp+dj]-var[indp-dj])*idy;
+          } else {
+            derivy = svj*0.5*(3*var[indp] - 4*var[indp-svj*dj]
+                              + var[indp-2*svj*dj])*idy;
+          }
+          
+          // Find z derivative
+          CCTK_REAL derivz;
+          if (kp>0 and kp<nk-1) {
+            derivz = 0.5*(var[indp+dk]-var[indp-dk])*idz;
+          } else {
+            derivz = svk*0.5*(3*var[indp] - 4*var[indp-svk*dk]
+                              + var[indp-2*svk*dk])*idz;
+          }
+          
+          // Find difference in sources
+          CCTK_REAL const aux =
+            rhs[indp] +
+            vx*derivx + vy*derivy + vz*derivz + v0*(var[indp] - var0)*rpi;
+          
+          // Extrapolate difference and add it to source in boundary
+          rhs[ind] += aux*pow(rp/r[ind],radpower);
+          
+        } // if radpower>=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<cGH const*> (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_REAL const*> (CCTK_VarDataPtr (cctkGH, 0, "grid::x"));
+  CCTK_REAL const* restrict const y =
+    static_cast<CCTK_REAL const*> (CCTK_VarDataPtr (cctkGH, 0, "grid::y"));
+  CCTK_REAL const* restrict const z =
+    static_cast<CCTK_REAL const*> (CCTK_VarDataPtr (cctkGH, 0, "grid::z"));
+  CCTK_REAL const* restrict const r =
+    static_cast<CCTK_REAL const*> (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;
+}