User: reisswig Date: 2011/06/13 04:44 PM
Added: /PITTNullCode/WorldTube/ README, configuration.ccl, interface.ccl, param.ccl, schedule.ccl /PITTNullCode/WorldTube/doc/ documentation.tex /PITTNullCode/WorldTube/src/ areal_sphere.cc, decompose.cc, make.code.defn, paramcheck.cc, register.cc, register.h
Removed: /WorldTube/
Log: Place thorn for metric decomposition in PITTNullCode repository.
File Changes:
Directory: /PITTNullCode/WorldTube/ ===================================
File [added]: README Delta lines: +12 -0 =================================================================== --- PITTNullCode/WorldTube/README (rev 0) +++ PITTNullCode/WorldTube/README 2011-06-13 21:44:24 UTC (rev 37) @@ -0,0 +1,12 @@ +CVS info : $Header:$ + +Cactus Code Thorn WorldTube +Thorn Author(s) : Christian Reisswig reisswig@aei.mpg.de +Thorn Maintainer(s) : Christian Reisswig reisswig@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 =================================================================== --- PITTNullCode/WorldTube/configuration.ccl (rev 0) +++ PITTNullCode/WorldTube/configuration.ccl 2011-06-13 21:44:24 UTC (rev 37) @@ -0,0 +1 @@ +REQUIRES THORNS: ADMBase ADMDerivatives
File [added]: interface.ccl Delta lines: +65 -0 =================================================================== --- PITTNullCode/WorldTube/interface.ccl (rev 0) +++ PITTNullCode/WorldTube/interface.ccl 2011-06-13 21:44:24 UTC (rev 37) @@ -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 =================================================================== --- PITTNullCode/WorldTube/param.ccl (rev 0) +++ PITTNullCode/WorldTube/param.ccl 2011-06-13 21:44:24 UTC (rev 37) @@ -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 =================================================================== --- PITTNullCode/WorldTube/schedule.ccl (rev 0) +++ PITTNullCode/WorldTube/schedule.ccl 2011-06-13 21:44:24 UTC (rev 37) @@ -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: /PITTNullCode/WorldTube/doc/ =======================================
File [added]: documentation.tex Delta lines: +144 -0 =================================================================== --- PITTNullCode/WorldTube/doc/documentation.tex (rev 0) +++ PITTNullCode/WorldTube/doc/documentation.tex 2011-06-13 21:44:24 UTC (rev 37) @@ -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/%7D%7D +% +% *======================================================================* + +% 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@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: /PITTNullCode/WorldTube/src/ =======================================
File [added]: areal_sphere.cc Delta lines: +226 -0 =================================================================== --- PITTNullCode/WorldTube/src/areal_sphere.cc (rev 0) +++ PITTNullCode/WorldTube/src/areal_sphere.cc 2011-06-13 21:44:24 UTC (rev 37) @@ -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 =================================================================== --- PITTNullCode/WorldTube/src/decompose.cc (rev 0) +++ PITTNullCode/WorldTube/src/decompose.cc 2011-06-13 21:44:24 UTC (rev 37) @@ -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 =================================================================== --- PITTNullCode/WorldTube/src/make.code.defn (rev 0) +++ PITTNullCode/WorldTube/src/make.code.defn 2011-06-13 21:44:24 UTC (rev 37) @@ -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 =================================================================== --- PITTNullCode/WorldTube/src/paramcheck.cc (rev 0) +++ PITTNullCode/WorldTube/src/paramcheck.cc 2011-06-13 21:44:24 UTC (rev 37) @@ -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 =================================================================== --- PITTNullCode/WorldTube/src/register.cc (rev 0) +++ PITTNullCode/WorldTube/src/register.cc 2011-06-13 21:44:24 UTC (rev 37) @@ -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 =================================================================== --- PITTNullCode/WorldTube/src/register.h (rev 0) +++ PITTNullCode/WorldTube/src/register.h 2011-06-13 21:44:24 UTC (rev 37) @@ -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