CarpetMask: New thorn

CarpetMask allows removing certain regions from Carpet's reduction
operators.  This can be used to excise punctures or other unwanted
features.

darcs-hash:20070421045211-dae7b-90791ae8ab35f639b8873144fd66f2587cc5e068.gz

10 files changed, 416 insertions, 0 deletions

diff --git a/Carpet/CarpetMask/README b/Carpet/CarpetMask/README
new file mode 100644
index 000000000..9c2fde769
--- /dev/null
+++ b/Carpet/CarpetMask/README
@@ -0,0 +1,9 @@
+Cactus Code Thorn CarpetMask
+Thorn Author(s)     : Erik Schnetter <schnetter@cct.lsu.edu>
+Thorn Maintainer(s) : Erik Schnetter <schnetter@cct.lsu.edu>
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
+Remove unwanted regions from Carpet's reduction operations.  This can
+be used e.g. to excise the interior of horizons.
diff --git a/Carpet/CarpetMask/doc/documentation.tex b/Carpet/CarpetMask/doc/documentation.tex
new file mode 100644
index 000000000..41c9885e9
--- /dev/null
+++ b/Carpet/CarpetMask/doc/documentation.tex
@@ -0,0 +1,144 @@
+% *======================================================================*
+%  Cactus Thorn template for ThornGuide documentation
+%  Author: Ian Kelley
+%  Date: Sun Jun 02, 2002
+%  $Header: /cactusdevcvs/Cactus/doc/ThornGuide/template.tex,v 1.12 2004/01/07 20:12:39 rideout Exp $
+%
+%  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: /cactusdevcvs/Cactus/doc/ThornGuide/template.tex,v 1.12 2004/01/07 20:12:39 rideout Exp $
+
+\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}
+
+% The title of the document (not necessarily the name of the Thorn)
+\title{CarpetMask}
+
+% the date your document was last changed, if your document is in CVS,
+% please use:
+%    \date{$ $Date: 2004/01/07 20:12:39 $ $}
+\date{April 20 2007}
+
+\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/Carpet/CarpetMask/interface.ccl b/Carpet/CarpetMask/interface.ccl
new file mode 100644
index 000000000..7297a1c46
--- /dev/null
+++ b/Carpet/CarpetMask/interface.ccl
@@ -0,0 +1,5 @@
+# Interface definition for thorn CarpetMask
+
+IMPLEMENTS: CarpetMask
+
+INHERITS: grid SphericalSurface
diff --git a/Carpet/CarpetMask/param.ccl b/Carpet/CarpetMask/param.ccl
new file mode 100644
index 000000000..bc8fe1979
--- /dev/null
+++ b/Carpet/CarpetMask/param.ccl
@@ -0,0 +1,53 @@
+# Parameter definitions for thorn CarpetMask
+
+
+
+# The exclusion parameters are re-interpreted after each regridding,
+# and are therefore marked as "always steerable".
+
+
+
+# Exclude coordinate spheres
+
+CCTK_REAL excluded_centre_x[10] "x coordinate of excluded region" STEERABLE=always
+{
+  *:* :: ""
+} 0.0
+
+CCTK_REAL excluded_centre_y[10] "y coordinate of excluded region" STEERABLE=always
+{
+  *:* :: ""
+} 0.0
+
+CCTK_REAL excluded_centre_z[10] "z coordinate of excluded region" STEERABLE=always
+{
+  *:* :: ""
+} 0.0
+
+CCTK_REAL excluded_radius[10] "radius of excluded region" STEERABLE=always
+{
+  -1.0  :: "do not exlude this region"
+  0.0:* :: ""
+} -1.0
+
+
+
+# Exclude spherical surfaces
+
+CCTK_INT excluded_surface[10] "index of excluded surface" STEERABLE=always
+{
+  -1  :: "do not exclude"
+  0:* :: ""
+} -1
+
+CCTK_REAL excluded_surface_factor[10] "shrink factor of excluded surface" STEERABLE=always
+{
+  0:* :: ""
+} 1.0
+
+SHARES: SphericalSurface
+
+USES CCTK_INT nsurfaces
+
+USES CCTK_INT maxntheta
+USES CCTK_INT maxnphi
diff --git a/Carpet/CarpetMask/schedule.ccl b/Carpet/CarpetMask/schedule.ccl
new file mode 100644
index 000000000..5e468b098
--- /dev/null
+++ b/Carpet/CarpetMask/schedule.ccl
@@ -0,0 +1,7 @@
+# Schedule definitions for thorn CarpetMask
+
+SCHEDULE CarpetExcludedSetup IN SetupMask
+{
+  LANG: C
+  OPTIONS: global loop-local
+} "Set up the weight function for the restriction regions"
diff --git a/Carpet/CarpetMask/src/make.code.defn b/Carpet/CarpetMask/src/make.code.defn
new file mode 100644
index 000000000..92d82ef80
--- /dev/null
+++ b/Carpet/CarpetMask/src/make.code.defn
@@ -0,0 +1,7 @@
+# Main make.code.defn file for thorn CarpetMask
+
+# Source files in this directory
+SRCS = mask_excluded.cc mask_surface.cc
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/Carpet/CarpetMask/src/mask_excluded.cc b/Carpet/CarpetMask/src/mask_excluded.cc
new file mode 100644
index 000000000..97268ef42
--- /dev/null
+++ b/Carpet/CarpetMask/src/mask_excluded.cc
@@ -0,0 +1,69 @@
+#include <cmath>
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include "mask_excluded.hh"
+
+
+
+namespace CarpetMask {
+  
+  using namespace std;
+  
+  
+  
+  /**
+   * Set the weight in the excluded regions to zero.
+   */
+  
+  void
+  CarpetExcludedSetup (CCTK_ARGUMENTS)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+    
+    CCTK_REAL * const weight =
+      static_cast <CCTK_REAL *>
+      (CCTK_VarDataPtr (cctkGH, 0, "CarpetReduce::weight"));
+    
+    if (not weight) {
+      CCTK_WARN (CCTK_WARN_ABORT,
+                 "CarpetReduce is not active, or CarpetReduce::mask does not have storage");
+    }
+    
+    for (int n = 0; n < 10; ++ n) {
+      
+      CCTK_REAL const r0 = excluded_radius[n];
+      if (r0 >= 0.0) {
+        
+        CCTK_REAL const x0 = excluded_centre_x[n];
+        CCTK_REAL const y0 = excluded_centre_y[n];
+        CCTK_REAL const z0 = excluded_centre_z[n];
+        
+        CCTK_REAL const r2 = pow (r0, 2);
+        
+        for (int k = 0; k < cctk_lsh[2]; ++ k) {
+          for (int j = 0; j < cctk_lsh[1]; ++ j) {
+            for (int i = 0; i < cctk_lsh[0]; ++ i) {
+              int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+              
+              CCTK_REAL const dx2 = pow (x[ind] - x0, 2);
+              CCTK_REAL const dy2 = pow (y[ind] - y0, 2);
+              CCTK_REAL const dz2 = pow (z[ind] - z0, 2);
+              
+              if (dx2 + dy2 + dz2 <= r2) {
+                weight[ind] = 0.0;
+              }
+              
+            }
+          }
+        }
+        
+      } // if r>=0
+    }   // for n
+    
+  }
+  
+} // namespace CarpetMask
diff --git a/Carpet/CarpetMask/src/mask_excluded.hh b/Carpet/CarpetMask/src/mask_excluded.hh
new file mode 100644
index 000000000..70edd8e8f
--- /dev/null
+++ b/Carpet/CarpetMask/src/mask_excluded.hh
@@ -0,0 +1,11 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+
+namespace CarpetMask {
+
+  extern "C" {
+    void
+    CarpetExcludedSetup (CCTK_ARGUMENTS);
+  }
+  
+} // namespace CarpetMask
diff --git a/Carpet/CarpetMask/src/mask_surface.cc b/Carpet/CarpetMask/src/mask_surface.cc
new file mode 100644
index 000000000..667a38643
--- /dev/null
+++ b/Carpet/CarpetMask/src/mask_surface.cc
@@ -0,0 +1,100 @@
+#include <cassert>
+#include <cmath>
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include "mask_surface.hh"
+
+
+
+namespace CarpetMask {
+  
+  using namespace std;
+  
+  
+  
+  /**
+   * Set the weight in the excluded regions to zero.
+   */
+  
+  void
+  CarpetSurfaceSetup (CCTK_ARGUMENTS)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+    
+    CCTK_REAL * const weight =
+      static_cast <CCTK_REAL *>
+      (CCTK_VarDataPtr (cctkGH, 0, "CarpetReduce::weight"));
+    
+    if (not weight) {
+      CCTK_WARN (CCTK_WARN_ABORT,
+                 "CarpetReduce is not active, or CarpetReduce::mask does not have storage");
+    }
+    
+    for (int n = 0; n < 10; ++ n) {
+      
+      int const sn = excluded_surface[n];
+      if (sn >= 0) {
+        
+        assert (sn < nsurfaces);
+        
+        CCTK_REAL const shrink_factor = excluded_surface_factor[n];
+        
+        if (sf_valid[sn] > 0) {
+          
+          CCTK_REAL const x0 = sf_origin_x[sn];
+          CCTK_REAL const y0 = sf_origin_y[sn];
+          CCTK_REAL const z0 = sf_origin_z[sn];
+          
+          CCTK_REAL const theta0 = sf_origin_theta[sn];
+          CCTK_REAL const phi0   = sf_origin_phi  [sn];
+          CCTK_REAL const dtheta = sf_delta_theta[sn];
+          CCTK_REAL const dphi   = sf_delta_phi  [sn];
+          
+          int const ntheta = sf_ntheta[sn];
+          int const nphi   = sf_nphi  [sn];
+          
+          for (int k = 0; k < cctk_lsh[2]; ++ k) {
+            for (int j = 0; j < cctk_lsh[1]; ++ j) {
+              for (int i = 0; i < cctk_lsh[0]; ++ i) {
+                int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+                
+                CCTK_REAL const dx = x[ind] - x0;
+                CCTK_REAL const dy = y[ind] - y0;
+                CCTK_REAL const dz = z[ind] - z0;
+                
+                CCTK_REAL const rho =
+                  sqrt (pow (dx, 2) + pow (dy, 2) + pow (dz, 2));
+                if (rho < 1.0e-12) {
+                  // Always excise the surface origin
+                  weight[ind] = 0.0;
+                } else {
+                  CCTK_REAL const theta = acos (dz / rho);
+                  CCTK_REAL const phi   = atan2 (dy, dx);
+                  int const a = floor ((theta - theta0) / dtheta + 0.5);
+                  int const b = floor ((phi   - phi0  ) / dphi   + 0.5);
+                  
+                  assert (a >= 0 and a < ntheta);
+                  assert (b >= 0 and b < nphi  );
+                  
+                  CCTK_REAL const dr = sf_radius[a + maxntheta * b];
+                  if (rho <= dr * shrink_factor) {
+                    weight[ind] = 0.0;
+                  }
+                }
+                
+              }
+            }
+          }
+          
+        } // if surface is valid
+        
+      } // if r>=0
+    }   // for n
+    
+  }
+  
+} // namespace CarpetMask
diff --git a/Carpet/CarpetMask/src/mask_surface.hh b/Carpet/CarpetMask/src/mask_surface.hh
new file mode 100644
index 000000000..a3e1f19d0
--- /dev/null
+++ b/Carpet/CarpetMask/src/mask_surface.hh
@@ -0,0 +1,11 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+
+namespace CarpetMask {
+
+  extern "C" {
+    void
+    CarpetSurfaceSetup (CCTK_ARGUMENTS);
+  }
+  
+} // namespace CarpetMask