[Commits] [svn:einsteintoolkit] incoming/WorldTube/ (Rev. 36)
reisswig at tapir.caltech.edu
reisswig at tapir.caltech.edu
Mon Jun 13 13:37:44 CDT 2011
User: reisswig
Date: 2011/06/13 01:37 PM
Added:
/WorldTube/
README, configuration.ccl, interface.ccl, param.ccl, schedule.ccl
/WorldTube/doc/
documentation.tex
/WorldTube/src/
areal_sphere.cc, decompose.cc, make.code.defn, paramcheck.cc, register.cc, register.h
Log:
A thorn capable of decomposing metric data on spheres.
File Changes:
Directory: /WorldTube/
======================
File [added]: README
Delta lines: +12 -0
===================================================================
--- WorldTube/README (rev 0)
+++ WorldTube/README 2011-06-13 18:37:44 UTC (rev 36)
@@ -0,0 +1,12 @@
+CVS info : $Header:$
+
+Cactus Code Thorn WorldTube
+Thorn Author(s) : Christian Reisswig <reisswig at aei.mpg.de>
+Thorn Maintainer(s) : Christian Reisswig <reisswig at aei.mpg.de>
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
+This thorn is ment as an interface between Cauchy and characteristic codes.
+
+It can also just output metric components and find spheres of constant areal radius.
\ No newline at end of file
File [added]: configuration.ccl
Delta lines: +1 -0
===================================================================
--- WorldTube/configuration.ccl (rev 0)
+++ WorldTube/configuration.ccl 2011-06-13 18:37:44 UTC (rev 36)
@@ -0,0 +1 @@
+REQUIRES THORNS: ADMBase ADMDerivatives
File [added]: interface.ccl
Delta lines: +65 -0
===================================================================
--- WorldTube/interface.ccl (rev 0)
+++ WorldTube/interface.ccl 2011-06-13 18:37:44 UTC (rev 36)
@@ -0,0 +1,65 @@
+# Interface definition for thorn WorldTube
+# $Header:$
+
+implements: WorldTube
+inherits: SphericalSlice ADMBase
+
+
+uses include: slices.hh
+
+
+
+# Register a Cactus variable for being sliced and return registration number
+CCTK_INT FUNCTION \
+ SphericalSlice_RegisterVariable \
+ (CCTK_POINTER_TO_CONST IN varname, \
+ CCTK_INT IN sn, \
+ CCTK_INT IN timelevels, \
+ CCTK_POINTER_TO_CONST IN distrib_method)
+REQUIRES FUNCTION SphericalSlice_RegisterVariable
+
+
+# Return an array of all values of scalar product with all spin-weighted spherical harmonics sYlm up to lmax
+CCTK_INT FUNCTION \
+ SphericalSlice_ContractVariableWithAllsYlm \
+ (CCTK_INT IN varno, \
+ CCTK_INT IN timelevel, \
+ CCTK_INT IN s, \
+ CCTK_INT IN lmin, \
+ CCTK_INT IN lmax, \
+ CCTK_COMPLEX OUT ARRAY coeffs)
+REQUIRES FUNCTION SphericalSlice_ContractVariableWithAllsYlm \
+
+
+# Synchronzie slice, i.e. interpolate a whole set of variables on the same slice
+CCTK_INT FUNCTION \
+ SphericalSlice_CollectiveSyncVariables \
+ (CCTK_POINTER_TO_CONST IN cctkGH_, \
+ CCTK_POINTER_TO_CONST IN varno, \
+ CCTK_INT IN number_of_vars)
+REQUIRES FUNCTION SphericalSlice_CollectiveSyncVariables \
+
+
+
+CCTK_INT FUNCTION Boundary_SelectVarForBC(CCTK_POINTER_TO_CONST IN GH, \
+ CCTK_INT IN faces, CCTK_INT IN boundary_width, CCTK_INT IN table_handle, \
+ CCTK_STRING IN var_name, CCTK_STRING IN bc_name)
+CCTK_INT FUNCTION Boundary_SelectGroupForBC(CCTK_POINTER_TO_CONST IN GH, \
+ CCTK_INT IN faces, CCTK_INT IN boundary_width, CCTK_INT IN table_handle, \
+ CCTK_STRING IN group_name, CCTK_STRING IN bc_name)
+
+REQUIRES FUNCTION Boundary_SelectVarForBC
+REQUIRES FUNCTION Boundary_SelectGroupForBC
+
+
+PUBLIC:
+
+COMPLEX extracted_vars[50] TYPE=ARRAY TIMELEVELS=1 DIM=2 SIZE=((lmax+1)*(lmax+1)),ntubes DISTRIB=CONSTANT
+{
+ extraction_vars
+} "Stores extracted lm-modes for each metric component"
+
+
+
+
+
File [added]: param.ccl
Delta lines: +67 -0
===================================================================
--- WorldTube/param.ccl (rev 0)
+++ WorldTube/param.ccl 2011-06-13 18:37:44 UTC (rev 36)
@@ -0,0 +1,67 @@
+# Parameter definitions for thorn WorldTube
+# $Header:$
+
+
+KEYWORD boundary_behavior "The type of boundary"
+{
+ "none" :: "don't extract/inject anything"
+ "extract metric" :: "extract metric components, lapse and shift"
+ "CCE" :: "extract metric components, lapse, shift and derivatives for CCE"
+ "CCE_cart" :: "extract metric components, lapse, shift and (cartesian) derivatives for CCE"
+ "CCM" :: "to be implemented"
+} "CCE"
+
+
+INT lmax "Maximum extracted l-mode"
+{
+ 0:* :: "any positive number"
+} 8
+
+
+CCTK_INT ntubes "The number of worldtubes to be placed"
+{
+ 0:42 :: "Any number between 0 and 42"
+} 1
+
+
+
+CCTK_INT which_slice_to_take[42] "Which SphericalSlice to use"
+{
+ 0:* :: "Any positive number"
+} 0
+
+
+BOOLEAN seek_const_areal_radius[42] "Set slice radius to a surface of constant areal radius before extracting/injecting"
+{
+} no
+
+
+REAL tolerance "Tolerance for finding areal radius"
+{
+ 0:* :: "Any positive number"
+} 1e-6
+
+
+INT max_it "Maximum number of iterations for finding areal radius"
+{
+ 1:* :: "any positive number greater than 1"
+} 10
+
+
+CCTK_INT radial_deriv_order "Order of radial derivatives"
+{
+ -1 :: "Use whatever is specified by SummationByParts"
+ 2:8 :: "Use this order"
+} -1
+
+
+
+SHARES: SphericalSlice
+USES KEYWORD type[42]
+USES REAL radius[42]
+USES INT nslices AS ss_nslices
+USES INT ntheta
+USES INT nphi
+USES INT nghostzones
+
+
File [added]: schedule.ccl
Delta lines: +45 -0
===================================================================
--- WorldTube/schedule.ccl (rev 0)
+++ WorldTube/schedule.ccl 2011-06-13 18:37:44 UTC (rev 36)
@@ -0,0 +1,45 @@
+# Schedule definitions for thorn WorldTube
+# $Header:$
+
+
+STORAGE: extracted_vars
+
+
+SCHEDULE WorldTube_ParamCheck AT PARAMCHECK
+{
+ LANG: C
+ OPTIONS: GLOBAL
+} "Check parameters"
+
+
+SCHEDULE WorldTube_RegisterSlices AT CCTK_INITIAL AFTER SphericalSlice_HasBeenSet
+{
+ LANG: C
+ OPTIONS: GLOBAL
+} "Register variables for extraction"
+
+
+SCHEDULE WorldTube_RegisterSlices AT CCTK_POST_RECOVER_VARIABLES AFTER SphericalSlice_HasBeenSet
+{
+ LANG: C
+ OPTIONS: GLOBAL
+} "Register variables for extraction"
+
+
+
+SCHEDULE WorldTube_ArealSphere AT CCTK_POSTSTEP AFTER SphericalSlice_HasBeenSet
+{
+ LANG: C
+ OPTIONS: GLOBAL
+} "Calculate surfaces of constant areal radius if selected"
+
+
+
+SCHEDULE WorldTube_Decompose AT CCTK_ANALYSIS AFTER (SphericalSlice_HasBeenSet, WorldTube_ArealSphere)
+{
+ LANG: C
+ OPTIONS: GLOBAL
+ TRIGGERS: extracted_vars
+} "Contract sliced variables with sYlm's"
+
+
Directory: /WorldTube/doc/
==========================
File [added]: documentation.tex
Delta lines: +144 -0
===================================================================
--- WorldTube/doc/documentation.tex (rev 0)
+++ WorldTube/doc/documentation.tex 2011-06-13 18:37:44 UTC (rev 36)
@@ -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{Christian Reisswig \textless reisswig at aei.mpg.de\textgreater}
+
+% The title of the document (not necessarily the name of the Thorn)
+\title{WorldTube}
+
+% the date your document was last changed, if your document is in CVS,
+% please use:
+% \date{$ $Date: 2004/01/07 20:12:39 $ $}
+\date{January 30 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}
Directory: /WorldTube/src/
==========================
File [added]: areal_sphere.cc
Delta lines: +226 -0
===================================================================
--- WorldTube/src/areal_sphere.cc (rev 0)
+++ WorldTube/src/areal_sphere.cc 2011-06-13 18:37:44 UTC (rev 36)
@@ -0,0 +1,226 @@
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "cctk_Functions.h"
+#include <assert.h>
+#include <iostream>
+
+#include "register.h"
+#include "slices.hh"
+
+
+
+using namespace WorldTube;
+using namespace SPS;
+
+double norm_inf(double const * const f, int const ni, int const nj);
+int get_1patch_metric_component(double g_sph[4][4], double const r, double const theta, double const phi,
+ const double lapse, const double shiftx, const double shifty, const double shiftz,
+ const double gxx, const double gxy, const double gxz,
+ const double gyy, const double gyz,
+ const double gzz);
+double calc_R(double const target_radius,
+ double const rp, double const tp, double const pp, const double g_sph[4][4],
+ double const diff_theta, double const diff_phi);
+
+extern "C" void WorldTube_ArealSphere(CCTK_ARGUMENTS)
+{
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+ int mode_dim = (lmax+1)*(lmax+1);
+
+
+ for (int ii=0; ii < ntubes; ++ii)
+ {
+ const int si = which_slice_to_take[ii];
+
+ if (seek_const_areal_radius[ii])
+ {
+ int it = 0;
+ double res_Linf = 2.0*tolerance;
+ while ((res_Linf > tolerance) && (it < max_it))
+ {
+ // first of all we do a collective interpolation for the current iterated radius onto the sphere
+ SphericalSlice_CollectiveSyncVariables(cctkGH, &sid[si].front(), 10);
+
+
+ if (CCTK_Equals(type[si], "1patch"))
+ {
+ vector<CCTK_REAL> R(slices_1patch(INDEX1P(ss_radius_id[si]), 0).lsh(0)[0] * slices_1patch(INDEX1P(ss_radius_id[si]), 0).lsh(0)[1], 0);
+ vector<CCTK_REAL> res(slices_1patch(INDEX1P(ss_radius_id[si]), 0).lsh(0)[0] * slices_1patch(INDEX1P(ss_radius_id[si]), 0).lsh(0)[1], 0);
+
+ for (spheredata_1patch<CCTK_REAL>::const_iter i=slices_1patch(INDEX1P(ss_radius_id[si]), 0).begin(); !i.done(); ++i)
+ {
+ // dont't integrate in ghostzones
+ if (i.ghostzone())
+ continue;
+
+ int const idx = i.idx().ij;
+
+ // get metric components at current point in polar-spherical basis
+ CCTK_REAL g_sph[4][4];
+ const CCTK_REAL lapse = slices_1patch(INDEX1P(sid[si][0]), 0)(i);
+ const CCTK_REAL shiftx = slices_1patch(INDEX1P(sid[si][1]), 0)(i);
+ const CCTK_REAL shifty = slices_1patch(INDEX1P(sid[si][2]), 0)(i);
+ const CCTK_REAL shiftz = slices_1patch(INDEX1P(sid[si][3]), 0)(i);
+ const CCTK_REAL gxx = slices_1patch(INDEX1P(sid[si][4]), 0)(i);
+ const CCTK_REAL gxy = slices_1patch(INDEX1P(sid[si][5]), 0)(i);
+ const CCTK_REAL gxz = slices_1patch(INDEX1P(sid[si][6]), 0)(i);
+ const CCTK_REAL gyy = slices_1patch(INDEX1P(sid[si][7]), 0)(i);
+ const CCTK_REAL gyz = slices_1patch(INDEX1P(sid[si][8]), 0)(i);
+ const CCTK_REAL gzz = slices_1patch(INDEX1P(sid[si][9]), 0)(i);
+ get_1patch_metric_component(g_sph, *i, i.idx().theta, i.idx().phi,
+ lapse,
+ shiftx, shifty, shiftz,
+ gxx, gxy, gxz,
+ gyy, gyz,
+ gzz);
+
+ // get angular derivatives of radius
+ CCTK_REAL diff_theta = slices_1patch(INDEX1P(ss_radius_id[si]), 0).dx(i);
+ CCTK_REAL diff_phi = slices_1patch(INDEX1P(ss_radius_id[si]), 0).dy(i);
+
+ R[idx] = calc_R(radius[si], *i, i.idx().theta, i.idx().phi, g_sph, diff_theta, diff_phi);
+ res[idx] = fabs(R[idx] - r[idx]);
+ }
+
+ res_Linf = norm_inf(&res.front(), slices_1patch(INDEX1P(ss_radius_id[si]), 0).lsh(0)[0],
+ slices_1patch(INDEX1P(ss_radius_id[si]), 0).lsh(0)[1]);
+
+ for (spheredata_1patch<CCTK_REAL>::iter i=slices_1patch(INDEX1P(ss_radius_id[si]), 0).begin(); !i.done(); ++i)
+ *i = R[i.idx().ij];
+
+ ++it;
+ }
+ if (CCTK_Equals(type[si], "6patch"))
+ CCTK_WARN(0, "Sorry: Finding surface of constant areal radius not implemented for 6-patches yet :(");
+
+
+ }
+ }
+
+ }
+}
+
+
+
+/*
+ * Infinity norm over the slice.
+ */
+double
+norm_inf(double const * const f, int const ni, int const nj)
+{
+ double norm = f[0];
+ for (int idx=0; idx<ni*nj; ++idx)
+ if (fabs(f[idx]) > norm)
+ norm = fabs(f[idx]);
+ double global_norm;
+ MPI_Allreduce (&norm, &global_norm, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
+ return global_norm;
+}
+
+
+/*
+ * Get the metric at a point on the sphere (1-patch coordinate basis).
+ *
+ * The metric components should be given in cartesian components to
+ * fill the g_cart[] array. This is transformed to spherical polar
+ * coordinates g_sph[] using the jacobian J[].
+ *
+ * In practice, the jacobian J[] will also need to recognise the possibility
+ * of alternate patch systems that are in use, eg. thornburg04 style 6-patch
+ * systems.
+ *
+ * The metric components should be obtained by interpolation from the
+ * numerically generated spacetime.
+ *
+ */
+int
+get_1patch_metric_component(double g_sph[4][4], double const r, double const theta, double const phi,
+ const double lapse, const double shiftx, const double shifty, const double shiftz,
+ const double gxx, const double gxy, const double gxz,
+ const double gyy, const double gyz,
+ const double gzz)
+{
+ double const J[4][4] = {
+ { 1, 0, 0, 0 },
+ { 0, sin(theta)*cos(phi), r*cos(theta)*cos(phi), -r*sin(theta)*sin(phi) },
+ { 0, sin(theta)*sin(phi), r*cos(theta)*sin(phi), r*sin(theta)*cos(phi) },
+ { 0, cos(theta), -r*sin(theta), 0 }
+ };
+
+ double const dt2 = -lapse*lapse
+ + gxx*shiftx*shiftx
+ + 2*gxy*shiftx*shifty
+ + 2*gxz*shiftx*shiftz
+ + gyy*shifty*shifty
+ + 2*gyz*shifty*shiftz
+ + gzz*shiftz*shiftz;
+
+ double const dtdx = gxx*shiftx+gxy*shifty+gxz*shiftz;
+ double const dtdy = gxy*shiftx+gyy*shifty+gyz*shiftz;
+ double const dtdz = gxz*shiftx+gyz*shifty+gzz*shiftz;
+
+ double const g_cart[4][4] = {
+ { dt2, dtdx, dtdy, dtdz },
+ { dtdx, gxx, gxy, gxz },
+ { dtdy, gxy, gyy, gyz },
+ { dtdz, gxz, gyz, gzz }
+ };
+
+ for (int a=0; a<4; ++a)
+ for (int b=0; b<4; ++b)
+ {
+ g_sph[a][b] = 0.0;
+ for (int ix=1; ix<4; ++ix)
+ for (int jx=1; jx<4; ++jx)
+ g_sph[a][b] += g_cart[ix][jx]*J[ix][a]*J[jx][b];
+ }
+
+
+ return 0;
+}
+
+
+
+/*
+ * Calculate the updated radius, R, at a point given an initial guess
+ * for the radius, r, and the metric on the sphere. See Nigel's notes
+ * for details on the computation.
+ */
+double
+calc_R(double const target_radius,
+ double const rp, double const tp, double const pp, const double g_sph[4][4],
+ double const diff_theta, double const diff_phi)
+{
+ double const grr = g_sph[1][1];
+ double const hr2 = g_sph[1][2] / rp;
+ double const hr3 = g_sph[1][3] / rp;
+ double const h22 = g_sph[2][2] / (rp*rp);
+ double const h23 = g_sph[2][3] / (rp*rp);
+ double const h33 = g_sph[3][3] / (rp*rp);
+
+ double const log_r_theta = diff_theta / rp;
+ double const log_r_phi = diff_phi / rp;
+
+ double const cg = sin(tp)*sin(tp) /
+ (h22*h33 - h23 * h23
+ + 2.0 * log_r_theta * hr2 * log_r_phi * hr3
+ + 2.0 * log_r_theta * hr2 * h33
+ + h22 * log_r_phi * log_r_phi * grr
+ + log_r_theta * log_r_theta * grr * h33
+ - 2.0 * log_r_theta * log_r_phi * grr * h23
+ - log_r_phi * log_r_phi * hr2 * hr2
+ + 2.0 * h22 * log_r_phi * hr3
+ - log_r_theta * log_r_theta * hr3 * hr3
+ - 2.0 * log_r_theta * hr3 * h23
+ - 2.0 * log_r_phi * hr2 * h23);
+
+ return target_radius * pow(cg, 0.25);
+}
+
+
+
+
+
File [added]: decompose.cc
Delta lines: +45 -0
===================================================================
--- WorldTube/src/decompose.cc (rev 0)
+++ WorldTube/src/decompose.cc 2011-06-13 18:37:44 UTC (rev 36)
@@ -0,0 +1,45 @@
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "cctk_Functions.h"
+#include <assert.h>
+#include <iostream>
+
+#include "register.h"
+
+
+
+
+
+namespace WorldTube {
+
+
+extern "C" void WorldTube_Decompose(CCTK_ARGUMENTS)
+{
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+ int mode_dim = (lmax+1)*(lmax+1);
+
+ // set all coeffs to zero initially
+ for (int i=0; i < 50*ntubes*(lmax+1)*(lmax+1); ++i)
+ extraction_vars[i] = CCTK_Cmplx(0,0);
+
+ for (int si=0; si < ntubes; ++si)
+ {
+
+ // first of all we do a collective interpolation onto the sphere
+ SphericalSlice_CollectiveSyncVariables(cctkGH, &sid[si].front(), extr_vars.size());
+
+ for (int j=0; j < extr_vars.size(); ++j)
+ {
+ SphericalSlice_ContractVariableWithAllsYlm(sid[si][j], 0, 0, 0, lmax, &extraction_vars[mode_dim*(si + ntubes*j)]);
+ }
+ }
+}
+
+
+}
+
+
+
File [added]: make.code.defn
Delta lines: +8 -0
===================================================================
--- WorldTube/src/make.code.defn (rev 0)
+++ WorldTube/src/make.code.defn 2011-06-13 18:37:44 UTC (rev 36)
@@ -0,0 +1,8 @@
+# Main make.code.defn file for thorn WorldTube
+# $Header:$
+
+# Source files in this directory
+SRCS = register.cc decompose.cc areal_sphere.cc paramcheck.cc
+
+# Subdirectories containing source files
+SUBDIRS =
File [added]: paramcheck.cc
Delta lines: +73 -0
===================================================================
--- WorldTube/src/paramcheck.cc (rev 0)
+++ WorldTube/src/paramcheck.cc 2011-06-13 18:37:44 UTC (rev 36)
@@ -0,0 +1,73 @@
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+
+
+
+extern "C" void WorldTube_ParamCheck(CCTK_ARGUMENTS)
+{
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+ // do this only on processor 0
+ if (CCTK_MyProc(cctkGH) != 0)
+ return;
+
+ for (int i=0; i < ntubes; ++i)
+ {
+ if (which_slice_to_take[i] >= ss_nslices)
+ CCTK_PARAMWARN("which_slice_to_take[i]: Slice number is greater than available slices!");
+
+ for (int j=0; j < ntubes; ++j)
+ if (i != j && which_slice_to_take[i] == which_slice_to_take[j])
+ CCTK_PARAMWARN("which_slice_to_take[i]: You have selected the same slice twice!");
+
+ // check, if the number of points required for the (lmax, mmax)-mode is met
+ // Eg: l=2, m=2 has 4 maxima -> we need at least twice as many points in order to resolve that
+ int required_points = lmax*4;
+ int recommended_points = lmax*8;
+
+ if (CCTK_Equals(type[which_slice_to_take[i]], "1patch"))
+ {
+ if (nphi[which_slice_to_take[i]] < required_points || 2*ntheta[which_slice_to_take[i]] < required_points)
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "For mode lmax=%d you need at least nphi=2*ntheta=%d to resolve it.", lmax, required_points);
+ if (nphi[which_slice_to_take[i]] < recommended_points || 2*ntheta[which_slice_to_take[i]] < recommended_points)
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "For mode lmax=%d I recommend nphi=2*ntheta=%d points.", lmax, recommended_points);
+ }
+ if (CCTK_Equals(type[which_slice_to_take[i]], "6patch"))
+ {
+ if (6*nphi[which_slice_to_take[i]] < required_points || 6*ntheta[which_slice_to_take[i]] < required_points)
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "For mode lmax=%d you need at least 6*nphi=6*ntheta=%d to resolve it.", lmax, required_points);
+ if (6*nphi[which_slice_to_take[i]] < recommended_points || 6*ntheta[which_slice_to_take[i]] < recommended_points)
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "For mode lmax=%d I recommend 6*nphi=6*ntheta=%d points.", lmax, recommended_points);
+ }
+
+ if (seek_const_areal_radius[i])
+ if (nghostzones[which_slice_to_take[i]] != 2)
+ CCTK_PARAMWARN("Ghostsize must be >= 2.");
+ }
+
+ if (CCTK_Equals(boundary_behavior, "CCE"))
+ {
+ int calc_dr = *((int*) CCTK_ParameterGet("store_radial_derivatives", "ADMDerivatives", NULL));
+ int calc_dt = *((int*) CCTK_ParameterGet("store_time_derivatives", "ADMDerivatives", NULL));
+
+ if (calc_dr == 0 || calc_dt == 0)
+ CCTK_WARN(0, "ADMDerivatives has to be set to calculate radial and time derivatives!");
+ }
+
+ if (CCTK_Equals(boundary_behavior, "CCE_cart"))
+ {
+ int calc_dx = *((int*) CCTK_ParameterGet("store_cartesian_derivatives", "ADMDerivatives", NULL));
+ int calc_dt = *((int*) CCTK_ParameterGet("store_time_derivatives", "ADMDerivatives", NULL));
+
+ if (calc_dx == 0 || calc_dt == 0)
+ CCTK_WARN(0, "ADMDerivatives has to be set to calculate cartesian and time derivatives!");
+ }
+
+}
File [added]: register.cc
Delta lines: +168 -0
===================================================================
--- WorldTube/src/register.cc (rev 0)
+++ WorldTube/src/register.cc 2011-06-13 18:37:44 UTC (rev 36)
@@ -0,0 +1,168 @@
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+
+#include "register.h"
+
+namespace WorldTube {
+
+
+using namespace std;
+
+
+// the slice ids for each slice and variable
+vector<vector<int> > sid;
+// all Cauchy extracted variables
+vector<string> extr_vars;
+
+
+
+extern "C" void WorldTube_RegisterSlices(CCTK_ARGUMENTS)
+{
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+
+ if (CCTK_Equals(boundary_behavior, "extract metric"))
+ {
+ extr_vars = vector<string>(10);
+ extr_vars[0] = "ADMBase::alp";
+ extr_vars[1] = "ADMBase::betax";
+ extr_vars[2] = "ADMBase::betay";
+ extr_vars[3] = "ADMBase::betaz";
+ extr_vars[4] = "ADMBase::gxx";
+ extr_vars[5] = "ADMBase::gxy";
+ extr_vars[6] = "ADMBase::gxz";
+ extr_vars[7] = "ADMBase::gyy";
+ extr_vars[8] = "ADMBase::gyz";
+ extr_vars[9] = "ADMBase::gzz";
+ }
+ if (CCTK_Equals(boundary_behavior, "CCE"))
+ {
+ extr_vars = vector<string>(30);
+ extr_vars[0] = "ADMBase::alp";
+ extr_vars[1] = "ADMBase::betax";
+ extr_vars[2] = "ADMBase::betay";
+ extr_vars[3] = "ADMBase::betaz";
+ extr_vars[4] = "ADMBase::gxx";
+ extr_vars[5] = "ADMBase::gxy";
+ extr_vars[6] = "ADMBase::gxz";
+ extr_vars[7] = "ADMBase::gyy";
+ extr_vars[8] = "ADMBase::gyz";
+ extr_vars[9] = "ADMBase::gzz";
+
+ // radial derivatives
+
+ extr_vars[10] = "ADMDerivatives::alp_dr";
+ extr_vars[11] = "ADMDerivatives::betax_dr";
+ extr_vars[12] = "ADMDerivatives::betay_dr";
+ extr_vars[13] = "ADMDerivatives::betaz_dr";
+ extr_vars[14] = "ADMDerivatives::gxx_dr";
+ extr_vars[15] = "ADMDerivatives::gxy_dr";
+ extr_vars[16] = "ADMDerivatives::gxz_dr";
+ extr_vars[17] = "ADMDerivatives::gyy_dr";
+ extr_vars[18] = "ADMDerivatives::gyz_dr";
+ extr_vars[19] = "ADMDerivatives::gzz_dr";
+
+ // time derivatives
+
+ // this can only be checked here (and not in parametercheck) because this is initialized in basegrid!
+ if (*dtlapse_state == 0 || *dtshift_state == 0)
+ CCTK_WARN(0, "You need to activate storage for ADMBase::dtalp and ADMBase::dtshift for CCE!");
+
+ extr_vars[20] = "ADMBase::dtalp";
+ extr_vars[21] = "ADMBase::dtbetax";
+ extr_vars[22] = "ADMBase::dtbetay";
+ extr_vars[23] = "ADMBase::dtbetaz";
+ extr_vars[24] = "ADMDerivatives::gxx_dt";
+ extr_vars[25] = "ADMDerivatives::gxy_dt";
+ extr_vars[26] = "ADMDerivatives::gxz_dt";
+ extr_vars[27] = "ADMDerivatives::gyy_dt";
+ extr_vars[28] = "ADMDerivatives::gyz_dt";
+ extr_vars[29] = "ADMDerivatives::gzz_dt";
+ }
+
+ if (CCTK_Equals(boundary_behavior, "CCE_cart"))
+ {
+ extr_vars = vector<string>(50);
+ extr_vars[0] = "ADMBase::alp";
+ extr_vars[1] = "ADMBase::betax";
+ extr_vars[2] = "ADMBase::betay";
+ extr_vars[3] = "ADMBase::betaz";
+ extr_vars[4] = "ADMBase::gxx";
+ extr_vars[5] = "ADMBase::gxy";
+ extr_vars[6] = "ADMBase::gxz";
+ extr_vars[7] = "ADMBase::gyy";
+ extr_vars[8] = "ADMBase::gyz";
+ extr_vars[9] = "ADMBase::gzz";
+
+ // cartesian derivatives
+
+ extr_vars[10] = "ADMDerivatives::alp_dx";
+ extr_vars[11] = "ADMDerivatives::betax_dx";
+ extr_vars[12] = "ADMDerivatives::betay_dx";
+ extr_vars[13] = "ADMDerivatives::betaz_dx";
+ extr_vars[14] = "ADMDerivatives::gxx_dx";
+ extr_vars[15] = "ADMDerivatives::gxy_dx";
+ extr_vars[16] = "ADMDerivatives::gxz_dx";
+ extr_vars[17] = "ADMDerivatives::gyy_dx";
+ extr_vars[18] = "ADMDerivatives::gyz_dx";
+ extr_vars[19] = "ADMDerivatives::gzz_dx";
+
+ extr_vars[20] = "ADMDerivatives::alp_dy";
+ extr_vars[21] = "ADMDerivatives::betax_dy";
+ extr_vars[22] = "ADMDerivatives::betay_dy";
+ extr_vars[23] = "ADMDerivatives::betaz_dy";
+ extr_vars[24] = "ADMDerivatives::gxx_dy";
+ extr_vars[25] = "ADMDerivatives::gxy_dy";
+ extr_vars[26] = "ADMDerivatives::gxz_dy";
+ extr_vars[27] = "ADMDerivatives::gyy_dy";
+ extr_vars[28] = "ADMDerivatives::gyz_dy";
+ extr_vars[29] = "ADMDerivatives::gzz_dy";
+
+ extr_vars[30] = "ADMDerivatives::alp_dz";
+ extr_vars[31] = "ADMDerivatives::betax_dz";
+ extr_vars[32] = "ADMDerivatives::betay_dz";
+ extr_vars[33] = "ADMDerivatives::betaz_dz";
+ extr_vars[34] = "ADMDerivatives::gxx_dz";
+ extr_vars[35] = "ADMDerivatives::gxy_dz";
+ extr_vars[36] = "ADMDerivatives::gxz_dz";
+ extr_vars[37] = "ADMDerivatives::gyy_dz";
+ extr_vars[38] = "ADMDerivatives::gyz_dz";
+ extr_vars[39] = "ADMDerivatives::gzz_dz";
+
+ // time derivatives
+
+ // this can only be checked here (and not in parametercheck) because this is initialized in basegrid!
+ if (*dtlapse_state == 0 || *dtshift_state == 0)
+ CCTK_WARN(0, "You need to activate storage for ADMBase::dtalp and ADMBase::dtshift for CCE!");
+
+ extr_vars[40] = "ADMBase::dtalp";
+ extr_vars[41] = "ADMBase::dtbetax";
+ extr_vars[42] = "ADMBase::dtbetay";
+ extr_vars[43] = "ADMBase::dtbetaz";
+ extr_vars[44] = "ADMDerivatives::gxx_dt";
+ extr_vars[45] = "ADMDerivatives::gxy_dt";
+ extr_vars[46] = "ADMDerivatives::gxz_dt";
+ extr_vars[47] = "ADMDerivatives::gyy_dt";
+ extr_vars[48] = "ADMDerivatives::gyz_dt";
+ extr_vars[49] = "ADMDerivatives::gzz_dt";
+ }
+
+
+ sid = vector<vector<int> >(ntubes);
+
+ for (int i=0; i < ntubes; ++i)
+ {
+ sid[i] = vector<int>(extr_vars.size(), 0);
+ for (int j=0; j < extr_vars.size(); ++j)
+ {
+ sid[i][j] = SphericalSlice_RegisterVariable(extr_vars[j].c_str(), which_slice_to_take[i], 1, "split");
+ }
+ }
+}
+
+
+
+}
+
File [added]: register.h
Delta lines: +22 -0
===================================================================
--- WorldTube/src/register.h (rev 0)
+++ WorldTube/src/register.h 2011-06-13 18:37:44 UTC (rev 36)
@@ -0,0 +1,22 @@
+#ifndef _WORLDTUBE_
+#define _WORLDTUBE_
+
+#include <vector>
+#include <iostream>
+
+namespace WorldTube {
+
+using namespace std;
+
+// the slice ids for each slice and variable
+extern vector<vector<int> > sid;
+
+// all Cauchy extracted variables
+extern vector<string> extr_vars;
+
+
+
+}
+
+
+#endif
More information about the Commits
mailing list