diff options
49 files changed, 2190 insertions, 3062 deletions
diff --git a/Auxiliary/Cactus/IOASCII.m b/Auxiliary/Cactus/IOASCII.m new file mode 100644 index 0000000..118db8d --- /dev/null +++ b/Auxiliary/Cactus/IOASCII.m @@ -0,0 +1,26 @@ + +BeginPackage["IOASCII`"]; + +ReadIOASCII::usage = "ReadIOASCII[file] reads the IOASCII file and parses it into list format"; +MGraph::usage = "MGraph[file] plots the IOASCII file"; + +Begin["`Private`"]; + +ReadIOASCII[file_] := + Module[{data1, data2, data3}, + data1 = Import[file, "Table"]; + data2 = Select[SplitBy[data1, Length[#] == 0 &], #[[1]] != {} &]; + data3 = Map[{First[#][[3]], Drop[#, 1]} &, data2]]; + +MGraph[file_String] := + Module[{data = ReadIOASCII[file]}, + MGraph[data]]; + +MGraph[data_List] := + Module[{}, + Manipulate[ + ListLinePlot[data[[it, 2]], PlotLabel -> data[[it, 1]]], {{it,1,"it"}, 1, Length[data], 1}]]; + +End[]; + +EndPackage[]; diff --git a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/interface.ccl b/Auxiliary/Cactus/KrancNumericalTools/Ceiling/interface.ccl deleted file mode 100644 index 8a9294b..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/interface.ccl +++ /dev/null @@ -1,8 +0,0 @@ -# file produced by user shusa, 31/3/2004 - -# $Id$ - -implements: Ceiling - -inherits: Grid - diff --git a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/param.ccl b/Auxiliary/Cactus/KrancNumericalTools/Ceiling/param.ccl deleted file mode 100644 index a6dc073..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/param.ccl +++ /dev/null @@ -1,32 +0,0 @@ -# file produced by user shusa, 31/3/2004 -# Produced with Mathematica Version 5.0 for Linux (June 9, 2003) - -# Mathematica script written by Ian Hinder and Sascha Husa - -# $Id$ - -private: -BOOLEAN check_active "whether to check ceiling values" -{ -} "false" - -REAL ceiling_value "what value we use for the ceiling" -{ -*:* :: "with negative values no cutoff is set" -} -1 - -KEYWORD type "what type of checking to apply" -{ - "absolute" :: "check for absolute value of GF" - "differential" :: "check for difference between min & max" -} "absolute" - -STRING vars "List of evolved grid functions that should have dissipation added" STEERABLE=always -{ - .* :: "must be a valid list of grid functions" -} "" - -BOOLEAN verbose "produce log output" STEERABLE=always -{ -} "no" - diff --git a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/schedule.ccl b/Auxiliary/Cactus/KrancNumericalTools/Ceiling/schedule.ccl deleted file mode 100644 index b9d167c..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/schedule.ccl +++ /dev/null @@ -1,21 +0,0 @@ -# file produced by user shusa, 31/3/2004 -# Produced with Mathematica Version 5.0 for Linux (June 9, 2003) - -# Mathematica script written by Ian Hinder and Sascha Husa - -# $Id$ - - -schedule Ceiling_Startup at startup -{ -LANG: C -} "ceiling startup message" - -if (check_active) -{ - schedule check_ceiling at PostStep - { - LANG: C - - } "check ceiling" -} diff --git a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/Ceiling.F90 b/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/Ceiling.F90 deleted file mode 100644 index 9e1d2c1..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/Ceiling.F90 +++ /dev/null @@ -1,68 +0,0 @@ -! file written by s. husa, 5/6/2004 - -! $Id$ - -#include "cctk.h" - -subroutine apply_check_abs(var, ni, nj, nk, ceiling_value) - -implicit none - -CCTK_INT, intent(in) :: ni, nj, nk -CCTK_REAL, dimension (ni, nj, nk), intent(in) :: var(ni, nj, nk) -CCTK_REAL, intent(in) :: ceiling_value - -CCTK_REAL :: criterion -CCTK_REAL, save :: initial_value - -CCTK_INT, save :: counter - -counter = counter + 1 -criterion = maxval(abs(var) + epsilon(1.0d0)) - -if (counter == 1) then - initial_value = criterion - write (*,*) "<<<<<< using ceiling initial value", initial_value -else - criterion = criterion / initial_value - if ((ceiling_value > 0).AND.(criterion > ceiling_value)) then - - call CCTK_INFO("Ceiling thorn terminates evolution") - call CCTK_TerminateNext(var) - endif -endif - -end subroutine apply_check_abs - - -subroutine apply_check_diff(var, ni, nj, nk, ceiling_value) - -implicit none - -CCTK_INT, intent(in) :: ni, nj, nk -CCTK_REAL, dimension (ni, nj, nk), intent(in) :: var(ni, nj, nk) -CCTK_REAL, intent(in) :: ceiling_value - -CCTK_REAL :: criterion -CCTK_REAL, save :: initial_value - -CCTK_INT, save :: counter - -counter = counter + 1 - -criterion = maxval(var) - minval(var) - -if (counter == 1) then - initial_value = criterion - write (*,*) "<<<<<< using ceiling initial value", initial_value -else - criterion = criterion / initial_value - if ((ceiling_value > 0).AND.(criterion > ceiling_value)) then - - call CCTK_INFO("Ceiling thorn terminates evolution.") - call CCTK_TerminateNext (var) - endif -endif - -end subroutine apply_check_diff - diff --git a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/Startup.c b/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/Startup.c deleted file mode 100644 index b42b2a1..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/Startup.c +++ /dev/null @@ -1,15 +0,0 @@ -/* file produced by user shusa, 31/3/2004 */ -/* Produced with Mathematica Version 5.0 for Linux (June 9, 2003) */ - -/* Mathematica script written by Ian Hinder and Sascha Husa */ - -/* $Id$ */ - -#include "cctk.h" - -int Ceiling_Startup(void) -{ - const char * banner = "Ceiling: abort when solution grows through the roof"; - CCTK_RegisterBanner(banner); - return 0; -} diff --git a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/make.code.defn b/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/make.code.defn deleted file mode 100644 index dd1377a..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/make.code.defn +++ /dev/null @@ -1,8 +0,0 @@ -# file produced by user shusa, 31/3/2004 -# Produced with Mathematica Version 5.0 for Linux (June 9, 2003) - -# Mathematica script written by Ian Hinder and Sascha Husa - -# $Id$ - -SRCS = Startup.c Ceiling.F90 selectGFs.c diff --git a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/selectGFs.c b/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/selectGFs.c deleted file mode 100644 index 86e1196..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/Ceiling/src/selectGFs.c +++ /dev/null @@ -1,94 +0,0 @@ -/* $Header$ */ - -/* this code is based on Erik Schnetter's dissipation thorn */ - -#include <assert.h> -#include <stdlib.h> - -#include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Parameters.h" - -void CCTK_FCALL -CCTK_FNAME(apply_check_abs) (CCTK_REAL const * const var, - int const * const ni, - int const * const nj, - int const * const nk, - CCTK_REAL const * const ceiling_value); - - -void CCTK_FCALL -CCTK_FNAME(apply_check_diff) (CCTK_REAL const * const var, - int const * const ni, - int const * const nj, - int const * const nk, - CCTK_REAL const * const ceiling_value); - -static void -call_apply_check (int const varindex, char const * const optstring, void * const arg); - -void -check_ceiling (CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - CCTK_TraverseString (vars, call_apply_check, cctkGH, CCTK_GROUP_OR_VAR); -} - - -void -call_apply_check (int const varindex, char const * const optstring, void * const arg) -{ - cGH const * const cctkGH = (cGH const *) arg; - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - int vargroup; - cGroup vardata; - - CCTK_REAL const * varptr; - int ierr /* , terminate */ ; - - assert (varindex >= 0); - - if (verbose) { - char * const fullvarname = CCTK_FullName (varindex); - assert (fullvarname); - CCTK_VInfo (CCTK_THORNSTRING, - "Applying ceiling check to \"%s\" ", - fullvarname); - free (fullvarname); - } - - vargroup = CCTK_GroupIndexFromVarI (varindex); - assert (vargroup >= 0); - - ierr = CCTK_GroupData (vargroup, &vardata); - assert (!ierr); - - assert (vardata.grouptype == CCTK_GF); - assert (vardata.vartype == CCTK_VARIABLE_REAL); - assert (vardata.dim == cctk_dim); - - varptr = CCTK_VarDataPtrI (cctkGH, 0, varindex); - assert (varptr); - - if (CCTK_Equals (type, "absolute")) - { - CCTK_FNAME(apply_check_abs) - (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2], &ceiling_value); - } - else if (CCTK_Equals (type, "differential")) - { - CCTK_FNAME(apply_check_diff) - (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2], &ceiling_value); - } - else - { - CCTK_INFO("keyword ceiling::type only allows values 'absolute' and 'differential'"); - } - - /* if (terminate > 0) {CCTK_TerminateNext (cctkGH);} */ -} - diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/param.ccl b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/param.ccl index 430c2f2..8e410f7 100644 --- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/param.ccl +++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/param.ccl @@ -4,38 +4,46 @@ # $Header$ -private: - -KEYWORD FDscheme "Type of finite differencing to use" -{ - "2nd order centered macro" :: "centered 2nd order implemented with macros" - "4th order centered macro" :: "centered 4th order implemented with macros" -} "2nd order centered macro" - restricted: -CCTK_INT stencil_width "stencil width used near boundary" +CCTK_INT stencil_width "stencil width used near boundary DEPRECATED" { -1:* :: "outgoing characteristic speed > 0, default of -1 is intentionally invalid" } -1 -CCTK_INT stencil_width_x "stencil width used near boundary" +CCTK_INT stencil_width_x "stencil width used near boundary DEPRECATED" { -1:* :: "outgoing characteristic speed > 0, default of -1 is intentionally invalid" } -1 -CCTK_INT stencil_width_y "stencil width used near boundary" +CCTK_INT stencil_width_y "stencil width used near boundary DEPRECATED" { -1:* :: "outgoing characteristic speed > 0, default of -1 is intentionally invalid" } -1 -CCTK_INT stencil_width_z "stencil width used near boundary" +CCTK_INT stencil_width_z "stencil width used near boundary DEPRECATED" { -1:* :: "outgoing characteristic speed > 0, default of -1 is intentionally invalid" } -1 -CCTK_INT boundary_width "width of boundary (fix later to use Cactus boundary calls)" +CCTK_INT boundary_width "width of boundary (fix later to use Cactus boundary calls) DEPRECATED" { -1:* :: "Any integer" } 1 +restricted: +CCTK_STRING jacobian_group "Name of group containing Jacobian" STEERABLE=RECOVER +{ + "" :: "String of the form <implementation>::<groupname>" +} "" + +restricted: +CCTK_STRING jacobian_derivative_group "Name of group containing Jacobian derivative" STEERABLE=RECOVER +{ + "" :: "String of the form <implementation>::<groupname>" +} "" + +CCTK_INT jacobian_identity_map "Map number on which the Jacobian should not be applied" +{ + -1:* :: "Any integer" +} -1 diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/schedule.ccl b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/schedule.ccl index a013b00..05ca0e1 100644 --- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/schedule.ccl +++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/schedule.ccl @@ -3,16 +3,3 @@ # author: S. Husa # $Header$ - - - -schedule GenericFD_Startup at STARTUP -{ - LANG: C -} "Register Banner" - -# schedule GenericFD_ParamCheck at ParamCheck -# { -# LANG: Fortran -# } "check stencil width parameters for consistency" - diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c index 4e1e774..8c06940 100644 --- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c +++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.c @@ -59,18 +59,28 @@ int sgn(CCTK_REAL x) return 0; } -int GenericFD_GetBoundaryWidth(cGH const * restrict const cctkGH) +void GenericFD_GetBoundaryWidths(cGH const * restrict const cctkGH, int nboundaryzones[6]) { int is_internal[6]; int is_staggered[6]; - int nboundaryzones[6]; int shiftout[6]; int ierr = -1; if (CCTK_IsFunctionAliased ("MultiPatch_GetBoundarySpecification")) { int const map = MultiPatch_GetMap (cctkGH); + /* This doesn't make sense in level mode */ if (map < 0) - CCTK_WARN(0, "Could not determine current map"); + { + static int have_warned = 0; + if (!have_warned) + { + CCTK_WARN(1, "GenericFD_GetBoundaryWidths: Could not determine current map (can be caused by calling in LEVEL mode)"); + have_warned = 1; + } + for (int i = 0; i < 6; i++) + nboundaryzones[i] = 0; + return; + } ierr = MultiPatch_GetBoundarySpecification (map, 6, nboundaryzones, is_internal, is_staggered, shiftout); if (ierr != 0) @@ -83,6 +93,12 @@ int GenericFD_GetBoundaryWidth(cGH const * restrict const cctkGH) } else { CCTK_WARN(0, "Could not obtain boundary specification"); } +} + +int GenericFD_GetBoundaryWidth(cGH const * restrict const cctkGH) +{ + int nboundaryzones[6]; + GenericFD_GetBoundaryWidths(cctkGH, nboundaryzones); int bw = nboundaryzones[0]; @@ -93,20 +109,6 @@ int GenericFD_GetBoundaryWidth(cGH const * restrict const cctkGH) return bw; } -/* int GenericFD_BoundaryWidthTable(cGH const * restrict const cctkGH) */ -/* { */ -/* int nboundaryzones[6]; */ -/* GenericFD_GetBoundaryWidth(cctkGH, nboundaryzones); */ - -/* int table = Util_TableCreate(0); */ -/* if (table < 0) CCTK_WARN(0, "Could not create table"); */ - -/* if (Util_TableSetIntArray(table, 6, nboundaryzones, "BOUNDARY_WIDTH") < 0) */ -/* CCTK_WARN(0, "Could not set table"); */ -/* return table; */ -/* } */ - - /* Return the array indices in imin and imax for looping over the interior of the grid. imin is the index of the first grid point. imax is the index of the last grid point plus 1. So a loop over @@ -219,7 +221,7 @@ void GenericFD_GetBoundaryInfo(cGH const * restrict const cctkGH, imin[dir] = npoints; break; case 1: /* Upper boundary */ - imax[dir] = cctk_lssh[CCTK_LSSH_IDX(0,dir)] - npoints; + imax[dir] = CCTK_LSSH(0,dir) - npoints; break; default: CCTK_WARN(0, "internal error"); @@ -239,7 +241,7 @@ void GenericFD_LoopOverEverything(cGH const * restrict const cctkGH, Kranc_Calcu CCTK_REAL tangentA[] = {0,0,0}; CCTK_REAL tangentB[] = {0,0,0}; int bmin[] = {0,0,0}; - int bmax[] = {cctk_lssh[CCTK_LSSH_IDX(0,0)],cctk_lssh[CCTK_LSSH_IDX(0,1)],cctk_lssh[CCTK_LSSH_IDX(0,2)]}; + int bmax[] = {CCTK_LSSH(0,0),CCTK_LSSH(0,1),CCTK_LSSH(0,2)}; calc(cctkGH, dir, face, normal, tangentA, tangentB, bmin, bmax, 0, NULL); return; @@ -301,7 +303,7 @@ void GenericFD_LoopOverBoundary(cGH const * restrict const cctkGH, Kranc_Calcula break; case +1: bmin[d] = imax[d]; - bmax[d] = cctk_lssh[CCTK_LSSH_IDX(0,d)]; + bmax[d] = CCTK_LSSH(0,d); have_bnd = 1; all_physbnd = all_physbnd && is_physbnd[2*d+1]; break; @@ -394,7 +396,7 @@ void GenericFD_LoopOverBoundaryWithGhosts(cGH const * restrict const cctkGH, Kra break; case +1: bmin[d] = imax[d]; - bmax[d] = cctk_lssh[CCTK_LSSH_IDX(0,d)]; + bmax[d] = CCTK_LSSH(0,d); have_bnd = 1; have_physbnd = have_physbnd || is_physbnd[2*d+1]; break; @@ -474,7 +476,7 @@ void GenericFD_PenaltyPrim2Char(cGH const * restrict const cctkGH, int const dir CCTK_REAL tangentA[] = {0,0,0}; CCTK_REAL tangentB[] = {0,0,0}; int bmin[] = {0,0,0}; - int bmax[] = {cctk_lssh[CCTK_LSSH_IDX(0,0)],cctk_lssh[CCTK_LSSH_IDX(0,1)],cctk_lssh[CCTK_LSSH_IDX(0,2)]}; + int bmax[] = {cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]}; CCTK_REAL **all_vars; int i = 0; @@ -499,3 +501,100 @@ void GenericFD_PenaltyPrim2Char(cGH const * restrict const cctkGH, int const dir return; } + +void GenericFD_AssertGroupStorage(cGH const * restrict const cctkGH, const char *calc, + int ngroups, const char *group_names[ngroups]) +{ + for (int i = 0; i < ngroups; i++) + { + int result = CCTK_QueryGroupStorage(cctkGH, group_names[i]); + if (result == 0) + { + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Error in %s: Group \"%s\" does not have storage", calc, group_names[i]); + } + else if (result < 0) + { + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Error in %s: Invalid group name \"%s\"", calc, group_names[i]); + } + } +} + +/* Return a list of pointers to the members of a named group */ +void GenericFD_GroupDataPointers(cGH const * restrict const cctkGH, const char *group_name, + int nvars, CCTK_REAL const *restrict *ptrs) +{ + int group_index, status; + cGroup group_info; + + group_index = CCTK_GroupIndex(group_name); + if (group_index < 0) + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Error return %d trying to get group index for group \'%s\'", + group_index, + group_name); + + status = CCTK_GroupData(group_index, &group_info); + if (status < 0) + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Error return %d trying to get info for group \'%s\'", + status, + group_name); + + if (group_info.numvars != nvars) + { + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Group \'%s\' has %d variables but %d were expected", + group_name, group_info.numvars, nvars); + } + + int v1 = CCTK_FirstVarIndex(group_name); + + for (int v = 0; v < nvars; v++) + { + ptrs[v] = (CCTK_REAL const *) CCTK_VarDataPtrI(cctkGH, 0 /* timelevel */, v1+v); + } +} + +void GenericFD_EnsureStencilFits(cGH const * restrict const cctkGH, const char *calc, int ni, int nj, int nk) +{ + DECLARE_CCTK_ARGUMENTS + + int bws[6]; + GenericFD_GetBoundaryWidths(cctkGH, bws); + + int ns[] = {ni, nj, nk}; + const char *dirs[] = {"x", "y", "z"}; + const char *faces[] = {"lower", "upper"}; + int abort = 0; + + for (int dir = 0; dir < 3; dir++) + { + for (int face = 0; face < 2; face++) + { + int bw = bws[2*dir+face]; + if (bw < ns[dir]) + { + CCTK_VInfo(CCTK_THORNSTRING, + "The stencil for %s requires %d points, but the %s %s boundary has only %d points.", + calc, ns[dir], faces[face], dirs[dir], bw); + abort = 1; + } + } + int gz = cctk_nghostzones[dir]; + if (gz < ns[dir]) + { + CCTK_VInfo(CCTK_THORNSTRING, + "The stencil for %s requires %d points, but there are only %d ghost zones in the %s direction.", + calc, ns[dir], gz, dirs[dir]); + abort = 1; + } + } + + if (abort) + { + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Insufficient ghost or boundary points for %s", calc); + } +} diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h index 041347d..401b4e0 100644 --- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h +++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h @@ -27,655 +27,33 @@ along with Kranc; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ - -#ifndef NOPRECOMPUTE -#define PRECOMPUTE -#endif - -#if defined(FD_C0) || defined(FD_C2) || defined(FD_C4) || defined(FD_C2C4) -#define FD_SET_BY_USER -#endif -#ifndef FD_SET_BY_USER -#define FD_C2 -#endif +#include "cctk.h" -#if defined(FD_C0) -#define FD_METHOD_DESC "FD method: replace derivatives by zero" +#ifdef __cplusplus +extern "C" { #endif -#if defined(FD_C2) -#define FD_METHOD_DESC "FD method: second order centered finite differences" -#endif - -#if defined(FD_C4) -#define FD_METHOD_DESC "FD method: fourth order centered finite differences" -#endif - -#if defined(FD_C2C4) -#define FD_METHOD_DESC "FD method: weighted 2nd/4th order centered finite differences" -#endif - - - -/* utility functions */ -#if defined(KRANC_C) -#define string(d,f) d ## f -#else -#define string(d,f) d/**/f +#ifdef __cplusplus +# ifdef CCTK_CXX_RESTRICT +# define restrict CCTK_CXX_RESTRICT +# endif #endif -#if defined(KRANC_C) -#define IMAX(int1, int2) ((int1) > (int2) ? (int1) : (int2)) -#endif - - #include "MathematicaCompat.h" -/* finite differencing macros */ - -/* */ -/* add method argument to shorthands */ -/* */ - -/* third derivatives */ -#define D111x(gf) string(D111,gf) -#define D211x(gf) string(D211,gf) -#define D311x(gf) string(D311,gf) -#define D221x(gf) string(D221,gf) -#define D321x(gf) string(D321,gf) -#define D331x(gf) string(D331,gf) -#define D222x(gf) string(D222,gf) -#define D322x(gf) string(D322,gf) -#define D332x(gf) string(D332,gf) -#define D333x(gf) string(D333,gf) - -/* second derivatives */ -#define D11x(gf) string(D11,gf) -#define D22x(gf) string(D22,gf) -#define D33x(gf) string(D33,gf) -#define D21x(gf) string(D21,gf) -#define D32x(gf) string(D32,gf) -#define D31x(gf) string(D31,gf) - -/* first derivatives */ -#define D1x(gf) string(D1,gf) -#define D2x(gf) string(D2,gf) -#define D3x(gf) string(D3,gf) - -#ifdef PRECOMPUTE - -/* third derivatives */ -#define D111(gf) string(D111,gf) -#define D211(gf) string(D211,gf) -#define D311(gf) string(D311,gf) -#define D221(gf) string(D221,gf) -#define D321(gf) string(D321,gf) -#define D331(gf) string(D331,gf) -#define D222(gf) string(D222,gf) -#define D322(gf) string(D322,gf) -#define D332(gf) string(D332,gf) -#define D333(gf) string(D333,gf) - -/* second derivatives */ -#define D11(gf,i,j,k) string(D11,gf) -#define D22(gf,i,j,k) string(D22,gf) -#define D33(gf,i,j,k) string(D33,gf) -#define D21(gf,i,j,k) string(D21,gf) -#define D32(gf,i,j,k) string(D32,gf) -#define D31(gf,i,j,k) string(D31,gf) - -/* first derivatives */ -#define D1(gf,i,j,k) string(D1,gf) -#define D2(gf,i,j,k) string(D2,gf) -#define D3(gf,i,j,k) string(D3,gf) - -#else - -/* third derivatives */ -#define D111(gf) D111gf(gf,i,j,k) -#define D211(gf) D211gf(gf,i,j,k) -#define D311(gf) D311gf(gf,i,j,k) -#define D221(gf) D221gf(gf,i,j,k) -#define D321(gf) D321gf(gf,i,j,k) -#define D331(gf) D331gf(gf,i,j,k) -#define D222(gf) D222gf(gf,i,j,k) -#define D322(gf) D322gf(gf,i,j,k) -#define D332(gf) D332gf(gf,i,j,k) -#define D333(gf) D333gf(gf,i,j,k) - -/* second derivatives */ -#define D11(gf,i,j,k) D11gf(gf,i,j,k) -#define D22(gf,i,j,k) D22gf(gf,i,j,k) -#define D33(gf,i,j,k) D33gf(gf,i,j,k) -#define D21(gf,i,j,k) D21gf(gf,i,j,k) -#define D32(gf,i,j,k) D32gf(gf,i,j,k) -#define D31(gf,i,j,k) D31gf(gf,i,j,k) - -/* first derivatives */ -#define D1(gf,i,j,k) D1gf(gf, i,j,k) -#define D2(gf,i,j,k) D2gf(gf, i,j,k) -#define D3(gf,i,j,k) D3gf(gf, i,j,k) - -#endif - - -#ifdef FD_C0 -/* third derivatives */ -#define D111gf(gf,i,j,k) D111_c0(gf,i,j,k) -#define D211gf(gf,i,j,k) D211_c0(gf,i,j,k) -#define D311gf(gf,i,j,k) D311_c0(gf,i,j,k) -#define D221gf(gf,i,j,k) D221_c0(gf,i,j,k) -#define D321gf(gf,i,j,k) D321_c0(gf,i,j,k) -#define D331gf(gf,i,j,k) D331_c0(gf,i,j,k) -#define D222gf(gf,i,j,k) D222_c0(gf,i,j,k) -#define D322gf(gf,i,j,k) D322_c0(gf,i,j,k) -#define D332gf(gf,i,j,k) D332_c0(gf,i,j,k) -#define D333gf(gf,i,j,k) D333_c0(gf,i,j,k) - -/* second derivatives */ -#define D11gf(gf,i,j,k) D11_c0(gf,i,j,k) -#define D22gf(gf,i,j,k) D22_c0(gf,i,j,k) -#define D33gf(gf,i,j,k) D33_c0(gf,i,j,k) -#define D21gf(gf,i,j,k) D21_c0(gf,i,j,k) -#define D32gf(gf,i,j,k) D32_c0(gf,i,j,k) -#define D31gf(gf,i,j,k) D31_c0(gf,i,j,k) - -/* first derivatives */ -#define D1gf(gf,i,j,k) D1_c0(gf, i,j,k) -#define D2gf(gf,i,j,k) D2_c0(gf, i,j,k) -#define D3gf(gf,i,j,k) D3_c0(gf, i,j,k) -#endif - - - -#ifdef FD_C2 -/* third derivatives */ -#define D111gf(gf,i,j,k) D111_c2(gf,i,j,k) -#define D211gf(gf,i,j,k) D211_c2(gf,i,j,k) -#define D311gf(gf,i,j,k) D311_c2(gf,i,j,k) -#define D221gf(gf,i,j,k) D221_c2(gf,i,j,k) -#define D321gf(gf,i,j,k) D321_c2(gf,i,j,k) -#define D331gf(gf,i,j,k) D331_c2(gf,i,j,k) -#define D222gf(gf,i,j,k) D222_c2(gf,i,j,k) -#define D322gf(gf,i,j,k) D322_c2(gf,i,j,k) -#define D332gf(gf,i,j,k) D332_c2(gf,i,j,k) -#define D333gf(gf,i,j,k) D333_c2(gf,i,j,k) - -/* second derivatives */ -#define D11gf(gf,i,j,k) D11_c2(gf,i,j,k) -#define D22gf(gf,i,j,k) D22_c2(gf,i,j,k) -#define D33gf(gf,i,j,k) D33_c2(gf,i,j,k) -#define D21gf(gf,i,j,k) D21_c2(gf,i,j,k) -#define D32gf(gf,i,j,k) D32_c2(gf,i,j,k) -#define D31gf(gf,i,j,k) D31_c2(gf,i,j,k) - -/* first derivatives */ -#define D1gf(gf,i,j,k) D1_c2(gf, i,j,k) -#define D2gf(gf,i,j,k) D2_c2(gf, i,j,k) -#define D3gf(gf,i,j,k) D3_c2(gf, i,j,k) -#endif - - - -#ifdef FD_C4 -/* third derivatives */ -#define D111gf(gf,i,j,k) D111_c4(gf,i,j,k) -#define D211gf(gf,i,j,k) D211_c4(gf,i,j,k) -#define D311gf(gf,i,j,k) D311_c4(gf,i,j,k) -#define D221gf(gf,i,j,k) D221_c4(gf,i,j,k) -#define D321gf(gf,i,j,k) D321_c4(gf,i,j,k) -#define D331gf(gf,i,j,k) D331_c4(gf,i,j,k) -#define D222gf(gf,i,j,k) D222_c4(gf,i,j,k) -#define D322gf(gf,i,j,k) D322_c4(gf,i,j,k) -#define D332gf(gf,i,j,k) D332_c4(gf,i,j,k) -#define D333gf(gf,i,j,k) D333_c4(gf,i,j,k) - -/* second derivatives */ -#define D11gf(gf,i,j,k) D11_c4(gf,i,j,k) -#define D22gf(gf,i,j,k) D22_c4(gf,i,j,k) -#define D33gf(gf,i,j,k) D33_c4(gf,i,j,k) -#define D21gf(gf,i,j,k) D21_c4(gf,i,j,k) -#define D32gf(gf,i,j,k) D32_c4(gf,i,j,k) -#define D31gf(gf,i,j,k) D31_c4(gf,i,j,k) - -/* first derivatives */ -#define D1gf(gf,i,j,k) D1_c4(gf, i,j,k) -#define D2gf(gf,i,j,k) D2_c4(gf, i,j,k) -#define D3gf(gf,i,j,k) D3_c4(gf, i,j,k) -#endif - - -#ifdef FD_C2C4 -/* third derivatives */ -#define D111gf(gf,i,j,k) D111_c2c4(gf,i,j,k) -#define D211gf(gf,i,j,k) D211_c2c4(gf,i,j,k) -#define D311gf(gf,i,j,k) D311_c2c4(gf,i,j,k) -#define D221gf(gf,i,j,k) D221_c2c4(gf,i,j,k) -#define D321gf(gf,i,j,k) D321_c2c4(gf,i,j,k) -#define D331gf(gf,i,j,k) D331_c2c4(gf,i,j,k) -#define D222gf(gf,i,j,k) D222_c2c4(gf,i,j,k) -#define D322gf(gf,i,j,k) D322_c2c4(gf,i,j,k) -#define D332gf(gf,i,j,k) D332_c2c4(gf,i,j,k) -#define D333gf(gf,i,j,k) D333_c2c4(gf,i,j,k) - -/* second derivatives */ -#define D11gf(gf,i,j,k) D11_c2c4(gf,i,j,k) -#define D22gf(gf,i,j,k) D22_c2c4(gf,i,j,k) -#define D33gf(gf,i,j,k) D33_c2c4(gf,i,j,k) -#define D21gf(gf,i,j,k) D21_c2c4(gf,i,j,k) -#define D32gf(gf,i,j,k) D32_c2c4(gf,i,j,k) -#define D31gf(gf,i,j,k) D31_c2c4(gf,i,j,k) - -/* first derivatives */ -#define D1gf(gf,i,j,k) D1_c2c4(gf, i,j,k) -#define D2gf(gf,i,j,k) D2_c2c4(gf, i,j,k) -#define D3gf(gf,i,j,k) D3_c2c4(gf, i,j,k) -#endif - - - - -/*****************************************************/ -/* */ -/* METHODS */ -/* */ -/*****************************************************/ - -/* c0 */ - -/* set all derivatives = 0 */ -/* for debugging and benchmarking */ - -/* third derivatives */ - -#define D111_c0(gf,i,j,k) 0. -#define D211_c0(gf,i,j,k) 0. -#define D311_c0(gf,i,j,k) 0. -#define D221_c0(gf,i,j,k) 0. -#define D321_c0(gf,i,j,k) 0. -#define D331_c0(gf,i,j,k) 0. -#define D222_c0(gf,i,j,k) 0. -#define D322_c0(gf,i,j,k) 0. -#define D332_c0(gf,i,j,k) 0. -#define D333_c0(gf,i,j,k) 0. - -/* second derivatives */ - -#define D11_c0(gf,i,j,k) 0. -#define D22_c0(gf,i,j,k) 0. -#define D33_c0(gf,i,j,k) 0. -#define D21_c0(gf,i,j,k) 0. -#define D32_c0(gf,i,j,k) 0. -#define D31_c0(gf,i,j,k) 0. - -/* first derivatives */ - -#define D1_c0(gf,i,j,k) 0. -#define D2_c0(gf,i,j,k) 0. -#define D3_c0(gf,i,j,k) 0. - - - -#ifndef KRANC_C - -/* c2 */ -/* */ -/* 2nd order centered */ -/* */ - -/* third derivatives, centered, 2nd order */ - -#define D111_c2(gf,i,j,k) ((- gf(i+2,j,k) + 2*gf(i+1,j,k) - 2*gf(i-1,j,k) + gf(i-2,j,k)) * dxi*dxi*dxi * (1.0/2.0)) -#define D211_c2(gf,i,j,k) ((gf(i+1,j+1,k) - 2*gf(i,j+1,k) + gf(i-1,j+1,k) - gf(i+1,j-1,k) + 2*gf(i,j-1,k) - gf(i-1,j-1,k)) * dxi*dxi*dyi * (1.0/2.0)) -#define D311_c2(gf,i,j,k) ((gf(i+1,j,k+1) - 2*gf(i,j,k+1) + gf(i-1,j,k+1) - gf(i+1,j,k-1) + 2*gf(i,j,k-1) - gf(i-1,j,k-1)) * dxi*dxi*dzi * (1.0/2.0)) -#define D221_c2(gf,i,j,k) ((gf(i+1,j+1,k) - 2*gf(i+1,j,k) + gf(i+1,j-1,k) - gf(i-1,j+1,k) + 2*gf(i-1,j,k) - gf(i-1,j-1,k)) * dxi*dyi*dyi * (1.0/2.0)) -#define D321_c2(gf,i,j,k) ((gf(i+1,j+1,k+1) - gf(i-1,j+1,k+1) - gf(i+1,j-1,k+1) + gf(i-1,j-1,k+1) - gf(i+1,j+1,k-1) + gf(i-1,j+1,k-1) + gf(i+1,j-1,k-1) - gf(i-1,j-1,k-1)) * dxi*dyi*dzi * (1.0/8.0)) -#define D331_c2(gf,i,j,k) ((gf(i+1,j,k+1) - 2*gf(i+1,j,k) + gf(i+1,j,k-1) - gf(i-1,j,k+1) + 2*gf(i-1,j,k) - gf(i-1,j,k-1)) * dxi*dzi*dzi * (1.0/2.0)) -#define D222_c2(gf,i,j,k) ((- gf(i,j+2,k) + 2*gf(i,j+1,k) - 2*gf(i,j-1,k) + gf(i,j-2,k)) * dyi*dyi*dyi * (1.0/2.0)) -#define D322_c2(gf,i,j,k) ((gf(i,j+1,k+1) - 2*gf(i,j,k+1) + gf(i,j-1,k+1) - gf(i,j+1,k-1) + 2*gf(i,j,k-1) - gf(i,j-1,k-1)) * dyi*dyi*dzi * (1.0/2.0)) -#define D332_c2(gf,i,j,k) ((gf(i,j+1,k+1) - 2*gf(i,j+1,k) + gf(i,j+1,k-1) - gf(i,j-1,k+1) + 2*gf(i,j-1,k) - gf(i,j-1,k-1)) * dyi*dzi*dzi * (1.0/2.0)) -#define D333_c2(gf,i,j,k) ((- gf(i,j,k+2) + 2*gf(i,j,k+1) - 2*gf(i,j,k-1) + gf(i,j,k-2)) * dzi*dzi*dzi * (1.0/2.0)) - -/* second derivatives, centered, 2nd order */ - -#define D11_c2(gf,i,j,k) \ - (( gf(i+1,j,k) \ - - 2.*gf(i, j,k) \ - + gf(i-1,j,k)) * dxi * dxi ) - -#define D22_c2(gf,i,j,k) \ - (( gf(i,j+1,k) \ - - 2.*gf(i,j, k) \ - + gf(i,j-1,k)) * dyi * dyi ) - -#define D33_c2(gf,i,j,k) \ - (( gf(i,j,k+1) \ - - 2.*gf(i,j,k ) \ - + gf(i,j,k-1)) * dzi * dzi ) - -#define D21_c2(gf,i,j,k) \ - ((gf(i+1,j+1,k) \ - + gf(i-1,j-1,k) \ - - gf(i+1,j-1,k) \ - - gf(i-1,j+1,k)) * hdxi * hdyi ) - -#define D32_c2(gf,i,j,k) \ - ((gf(i,j+1,k+1) \ - + gf(i,j-1,k-1) \ - - gf(i,j+1,k-1) \ - - gf(i,j-1,k+1)) * hdyi * hdzi ) - -#define D31_c2(gf,i,j,k) \ - ((gf(i+1,j,k+1) \ - + gf(i-1,j,k-1) \ - - gf(i+1,j,k-1) \ - - gf(i-1,j,k+1)) * hdxi * hdzi ) - -/* first derivatives, centered, 2nd order */ - -#define D1_c2(gf,i,j,k) \ - ((gf(i+1,j,k) \ - - gf(i-1,j,k)) * hdxi) - -#define D2_c2(gf,i,j,k) \ - ((gf(i,j+1,k) \ - - gf(i,j-1,k)) * hdyi) - -#define D3_c2(gf,i,j,k) \ - ((gf(i,j,k+1) \ - - gf(i,j,k-1)) * hdzi) - -#else - - -#define D11_c2(gf,i,j,k) \ - (( gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k)] \ - - 2.*gf[CCTK_GFINDEX3D(cctkGH,i, j,k)] \ - + gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k)]) * dxi * dxi ) - -#define D22_c2(gf,i,j,k) \ - (( gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k)] \ - - 2.*gf[CCTK_GFINDEX3D(cctkGH,i,j, k)] \ - + gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k)]) * dyi * dyi ) - -#define D33_c2(gf,i,j,k) \ - (( gf[CCTK_GFINDEX3D(cctkGH,i,j,k+1)] \ - - 2.*gf[CCTK_GFINDEX3D(cctkGH,i,j,k )] \ - + gf[CCTK_GFINDEX3D(cctkGH,i,j,k-1)]) * dzi * dzi ) - -#define D21_c2(gf,i,j,k) \ - ((gf[CCTK_GFINDEX3D(cctkGH,i+1,j+1,k)] \ - + gf[CCTK_GFINDEX3D(cctkGH,i-1,j-1,k)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i+1,j-1,k)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i-1,j+1,k)]) * hdxi * hdyi ) - -#define D32_c2(gf,i,j,k) \ - ((gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k+1)] \ - + gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k-1)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k-1)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k+1)]) * hdyi * hdzi ) - -#define D31_c2(gf,i,j,k) \ - ((gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k+1)] \ - + gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k-1)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k-1)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k+1)]) * hdxi * hdzi ) - -/* first derivatives, centered, 2nd order */ - -#define D1_c2(gf,i,j,k) \ - ((gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k)]) * hdxi) - -#define D2_c2(gf,i,j,k) \ - ((gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k)]) * hdyi) - -#define D3_c2(gf,i,j,k) \ - ((gf[CCTK_GFINDEX3D(cctkGH,i,j,k+1)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i,j,k-1)]) * hdzi) - -#endif - - -/* c4 */ -/* */ -/* 4th order centered */ -/* */ - -/* second derivatives, centered, 4th order */ - -#ifndef KRANC_C - -#define D11_c4(gf,i,j,k) \ - ((- gf(i+2,j,k) \ - + 16. * gf(i+1,j,k) \ - - 30. * gf(i, j,k) \ - + 16. * gf(i-1,j,k) \ - - gf(i-2,j,k)) * dxi * dxi / 12.) - - -#define D22_c4(gf,i,j,k) \ - ((- gf(i,j+2,k) \ - + 16. * gf(i,j+1,k) \ - - 30. * gf(i,j, k) \ - + 16. * gf(i,j-1,k) \ - - gf(i,j-2,k)) * dyi * dyi / 12.) - - -#define D33_c4(gf,i,j,k) \ - ((- gf(i,j,k+2) \ - + 16. * gf(i,j,k+1) \ - - 30. * gf(i,j,k ) \ - + 16. * gf(i,j,k-1) \ - - gf(i,j,k-2)) * dzi * dzi / 12.) - -#define D21_c4(gf,i,j,k) \ - ((- gf(i+2,j+2,k) \ - + gf(i+2,j-2,k) \ - + gf(i-2,j+2,k) \ - - gf(i-2,j-2,k) \ - + 16. * gf(i+1,j+1,k) \ - - 16. * gf(i+1,j-1,k) \ - - 16. * gf(i-1,j+1,k) \ - + 16. * gf(i-1,j-1,k)) * dxi * dyi / 48.) - -#define D31_c4(gf,i,j,k) \ - ((- gf(i+2,j,k+2) \ - + gf(i+2,j,k-2) \ - + gf(i-2,j,k+2) \ - - gf(i-2,j,k-2) \ - + 16. * gf(i+1,j,k+1) \ - - 16. * gf(i+1,j,k-1) \ - - 16. * gf(i-1,j,k+1) \ - + 16. * gf(i-1,j,k-1)) * dxi * dzi / 48.) - - -#define D32_c4(gf,i,j,k) \ - ((- gf(i,j+2,k+2) \ - + gf(i,j+2,k-2) \ - + gf(i,j-2,k+2) \ - - gf(i,j-2,k-2) \ - + 16. * gf(i,j+1,k+1) \ - - 16. * gf(i,j+1,k-1) \ - - 16. * gf(i,j-1,k+1) \ - + 16. * gf(i,j-1,k-1)) * dzi * dyi / 48.) - - -/* first derivatives, centered, 4th order */ - -#define D1_c4(gf,i,j,k) \ - ((- gf(i+2,j,k) \ - + 8. * gf(i+1,j,k) \ - - 8. * gf(i-1,j,k) \ - + gf(i-2,j,k)) * (dxi / 12.)) - -#define D2_c4(gf,i,j,k) \ - ((- gf(i,j+2,k) \ - + 8. * gf(i,j+1,k) \ - - 8. * gf(i,j-1,k) \ - + gf(i,j-2,k)) * (dyi / 12.)) - -#define D3_c4(gf,i,j,k) \ - ((- gf(i,j,k+2) \ - + 8. * gf(i,j,k+1) \ - - 8. * gf(i,j,k-1) \ - + gf(i,j,k-2)) * (dzi / 12.)) - - -#else - -#define D11_c4(gf,i,j,k) \ - ((- gf[CCTK_GFINDEX3D(cctkGH,i+2,j,k)] \ - + 16. * gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k)] \ - - 30. * gf[CCTK_GFINDEX3D(cctkGH,i, j,k)] \ - + 16. * gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i-2,j,k)]) * dxi * dxi / 12.) - - -#define D22_c4(gf,i,j,k) \ - ((- gf[CCTK_GFINDEX3D(cctkGH,i,j+2,k)] \ - + 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k)] \ - - 30. * gf[CCTK_GFINDEX3D(cctkGH,i,j, k)] \ - + 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i,j-2,k)]) * dyi * dyi / 12.) - - -#define D33_c4(gf,i,j,k) \ - ((- gf[CCTK_GFINDEX3D(cctkGH,i,j,k+2)] \ - + 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j,k+1)] \ - - 30. * gf[CCTK_GFINDEX3D(cctkGH,i,j,k )] \ - + 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j,k-1)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i,j,k-2)]) * dzi * dzi / 12.) - -#define D21_c4(gf,i,j,k) \ - ((- gf[CCTK_GFINDEX3D(cctkGH,i+2,j+2,k)] \ - + gf[CCTK_GFINDEX3D(cctkGH,i+2,j-2,k)] \ - + gf[CCTK_GFINDEX3D(cctkGH,i-2,j+2,k)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i-2,j-2,k)] \ - + 16. * gf[CCTK_GFINDEX3D(cctkGH,i+1,j+1,k)] \ - - 16. * gf[CCTK_GFINDEX3D(cctkGH,i+1,j-1,k)] \ - - 16. * gf[CCTK_GFINDEX3D(cctkGH,i-1,j+1,k)] \ - + 16. * gf[CCTK_GFINDEX3D(cctkGH,i-1,j-1,k)]) * dxi * dyi / 48.) - -#define D31_c4(gf,i,j,k) \ - ((- gf[CCTK_GFINDEX3D(cctkGH,i+2,j,k+2)] \ - + gf[CCTK_GFINDEX3D(cctkGH,i+2,j,k-2)] \ - + gf[CCTK_GFINDEX3D(cctkGH,i-2,j,k+2)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i-2,j,k-2)] \ - + 16. * gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k+1)] \ - - 16. * gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k-1)] \ - - 16. * gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k+1)] \ - + 16. * gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k-1)]) * dxi * dzi / 48.) - - -#define D32_c4(gf,i,j,k) \ - ((- gf[CCTK_GFINDEX3D(cctkGH,i,j+2,k+2)] \ - + gf[CCTK_GFINDEX3D(cctkGH,i,j+2,k-2)] \ - + gf[CCTK_GFINDEX3D(cctkGH,i,j-2,k+2)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i,j-2,k-2)] \ - + 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k+1)] \ - - 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k-1)] \ - - 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k+1)] \ - + 16. * gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k-1)]) * dzi * dyi / 48.) - - -/* first derivatives, centered, 4th order */ - -#define D1_c4(gf,i,j,k) \ - ((- gf[CCTK_GFINDEX3D(cctkGH,i+2,j,k)] \ - + 8. * gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k)] \ - - 8. * gf[CCTK_GFINDEX3D(cctkGH,i-1,j,k)] \ - + gf[CCTK_GFINDEX3D(cctkGH,i-2,j,k)]) * (dxi / 12.)) - -#define D2_c4(gf,i,j,k) \ - ((- gf[CCTK_GFINDEX3D(cctkGH,i,j+2,k)] \ - + 8. * gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k)] \ - - 8. * gf[CCTK_GFINDEX3D(cctkGH,i,j-1,k)] \ - + gf[CCTK_GFINDEX3D(cctkGH,i,j-2,k)]) * (dyi / 12.)) - -#define D3_c4(gf,i,j,k) \ - ((- gf[CCTK_GFINDEX3D(cctkGH,i,j,k+2)] \ - + 8. * gf[CCTK_GFINDEX3D(cctkGH,i,j,k+1)] \ - - 8. * gf[CCTK_GFINDEX3D(cctkGH,i,j,k-1)] \ - + gf[CCTK_GFINDEX3D(cctkGH,i,j,k-2)]) * (dzi / 12.)) - -#endif - -/*****************************************************/ -/* */ -/* DERIVED METHODS */ -/* */ -/******************************************************/ - - -/* blend c2 and c4 */ -/* second derivatives */ -#define D11_c2c4(gf,i,j,k) (fdweight_c2*D11_c2(gf,i,j,k) + fdweight_c4*D11_c4(gf,i,j,k)) -#define D22_c2c4(gf,i,j,k) (fdweight_c2*D22_c2(gf,i,j,k) + fdweight_c4*D22_c4(gf,i,j,k)) -#define D33_c2c4(gf,i,j,k) (fdweight_c2*D33_c2(gf,i,j,k) + fdweight_c4*D33_c4(gf,i,j,k)) -#define D21_c2c4(gf,i,j,k) (fdweight_c2*D21_c2(gf,i,j,k) + fdweight_c4*D21_c4(gf,i,j,k)) -#define D32_c2c4(gf,i,j,k) (fdweight_c2*D32_c2(gf,i,j,k) + fdweight_c4*D32_c4(gf,i,j,k)) -#define D31_c2c4(gf,i,j,k) (fdweight_c2*D31_c2(gf,i,j,k) + fdweight_c4*D31_c4(gf,i,j,k)) - -/* first derivatives */ -#define D1_c2c4(gf,i,j,k) (fdweight_c2*D1_c2(gf, i,j,k) + fdweight_c4*D1_c4(gf,i,j,k)) -#define D2_c2c4(gf,i,j,k) (fdweight_c2*D2_c2(gf, i,j,k) + fdweight_c4*D2_c4(gf,i,j,k)) -#define D3_c2c4(gf,i,j,k) (fdweight_c2*D3_c2(gf, i,j,k) + fdweight_c4*D3_c4(gf,i,j,k)) - - -/*****************************************************/ -/* */ -/* Poor man's one-sided derivatives */ -/* */ -/******************************************************/ - - -#define Dplus1(gf,i,j,k) string(Dplus1,gf) -#define Dplus2(gf,i,j,k) string(Dplus2,gf) -#define Dplus3(gf,i,j,k) string(Dplus3,gf) - - -#define Dplus1gf(gf,i,j,k) Dplus1x(gf, i,j,k) -#define Dplus2gf(gf,i,j,k) Dplus2x(gf, i,j,k) -#define Dplus3gf(gf,i,j,k) Dplus3x(gf, i,j,k) - -#ifdef KRANC_C - -#define Dplus1x(gf,i,j,k) \ - ((gf[CCTK_GFINDEX3D(cctkGH,i+1,j,k)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i,j,k)]) * dxi) - -#define Dplus2x(gf,i,j,k) \ - ((gf[CCTK_GFINDEX3D(cctkGH,i,j+1,k)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i,j,k)]) * dyi) - -#define Dplus3x(gf,i,j,k) \ - ((gf[CCTK_GFINDEX3D(cctkGH,i,j,k+1)] \ - - gf[CCTK_GFINDEX3D(cctkGH,i,j,k)]) * dzi) - -#else - -#define Dplus1x(gf,i,j,k) ( ( gf(i+1, j, k) - gf(i,j,k) ) * dxi ) -#define Dplus2x(gf,i,j,k) ( ( gf(i, j + 1, k) - gf(i,j,k) ) * dxi ) -#define Dplus3x(gf,i,j,k) ( ( gf(i, j, k + 1) - gf(i,j,k) ) * dxi ) - -#endif - #ifdef KRANC_C int sgn(CCTK_REAL x); -#define Dupwind1(gf,dir,i,j,k) ((dir * gf[CCTK_GFINDEX3D(cctkGH,i+dir,j,k)] \ - - dir * gf[CCTK_GFINDEX3D(cctkGH,i,j,k)]) * dxi) -#define Dupwind2(gf,dir,i,j,k) ((dir * gf[CCTK_GFINDEX3D(cctkGH,i,j+dir,k)] \ - - dir * gf[CCTK_GFINDEX3D(cctkGH,i,j,k)]) * dxi) -#define Dupwind3(gf,dir,i,j,k) ((dir * gf[CCTK_GFINDEX3D(cctkGH,i,j,k+dir)] \ - - dir * gf[CCTK_GFINDEX3D(cctkGH,i,j,k)]) * dxi) - int GenericFD_GetBoundaryWidth(cGH const * restrict const cctkGH); -/* int GenericFD_BoundaryWidthTable(cGH const * restrict const cctkGH); */ +#ifdef __cplusplus +// Define the restrict qualifier +# ifdef CCTK_CXX_RESTRICT +# undef restrict +# define restrict CCTK_CXX_RESTRICT +# endif +#endif void GenericFD_GetBoundaryInfo(cGH const * restrict cctkGH, int const * restrict cctk_lsh, @@ -688,24 +66,8 @@ void GenericFD_GetBoundaryInfo(cGH const * restrict cctkGH, int * restrict is_physbnd, int * restrict is_ipbnd); -#if 0 -/* Finite differencing near boundaries */ - -/* The array var is to be accessed at the location - [i+ioff,j+joff,k+koff]. idir,jdir,kdir specify whether there is a - lower (dir<0), upper (dir>0), or no boundary nearby. If a boundary - is in the way, the value 0 is returned instead of the array - content. */ -#define CCTK_GFACCESS3D(cctkGH, var, i,j,k, ioff,joff,koff, idir,jdir,kdir) \ - (((idir)<0 && (ioff)<0) || \ - ((jdir)<0 && (joff)<0) || \ - ((kdir)<0 && (koff)<0) || \ - ((idir)>0 && (ioff)>0) || \ - ((jdir)>0 && (joff)>0) || \ - ((kdir)>0 && (koff)>0) || \ - ? 0 \ - : (var)[CCTK_GFINDEX3D((cctkGH), (i)+(ioff),(j)+(joff),(k)+(koff))]) -#endif +void GenericFD_AssertGroupStorage(cGH const * restrict const cctkGH, const char *calc, + int ngroups, const char *group_names[]); /* Summation by parts */ @@ -790,40 +152,13 @@ void GenericFD_LoopOverBoundary(cGH const * restrict cctkGH, Kranc_Calculation c void GenericFD_LoopOverBoundaryWithGhosts(cGH const * restrict cctkGH, Kranc_Calculation calc); void GenericFD_LoopOverInterior(cGH const * restrict cctkGH, Kranc_Calculation calc); +void GenericFD_GroupDataPointers(cGH const * restrict const cctkGH, const char *group_name, + int nvars, CCTK_REAL const *restrict *ptrs); +void GenericFD_EnsureStencilFits(cGH const * restrict const cctkGH, const char *calc, int ni, int nj, int nk); -/* Vectorisation of memory accesses */ - -#include <stdlib.h> -#include <cctk.h> - -#if defined(__SSE2__) && defined(CCTK_REAL_PRECISION_8) - -#include <emmintrin.h> - -/* A vector type corresponding to CCTK_REAL */ -typedef __m128d CCTK_REAL_VEC; - -/* Really only SSE is required, but there doesn't seem to be a - preprocessing flag to check for this */ -#elif defined(__SSE2__) && defined(CCTK_REAL_PRECISION_4) - -#include <emmintrin.h> - -/* A vector type corresponding to CCTK_REAL */ -typedef __m128 CCTK_REAL_VEC; - -#else - -/* There is no vector type corresponding to CCTK_REAL */ -typedef CCTK_REAL CCTK_REAL_VEC; - +#ifdef __cplusplus +} /* extern "C" */ #endif -/* The number of vector elements in a CCTK_REAL_VEC */ -static -size_t const CCTK_REAL_VEC_SIZE = sizeof(CCTK_REAL_VEC) / sizeof(CCTK_REAL); - - - #endif diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/MathematicaCompat.h b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/MathematicaCompat.h index 86b70eb..e48daa2 100644 --- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/MathematicaCompat.h +++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/MathematicaCompat.h @@ -1,19 +1,19 @@ -#define Power(x, y) (pow((x), (y))) +#define Power(x, y) (pow(x,y)) #define Sqrt(x) (sqrt(x)) #ifdef KRANC_C -#define Abs(x) (fabs(x)) -#define Min(x, y) (fmin((x), (y))) -#define Min3(x, y, z) (fmin(fmin((x), (y)), (z))) -#define Max(x, y) (fmax((x), (y))) -#define IfThen(x,y,z) ((x) ? (y) : (z)) +# define Abs(x) (fabs(x)) +# define Min(x, y) (fmin(x,y)) +# define Min3(x, y, z) (fmin(fmin((x), (y)), (z))) +# define Max(x, y) (fmax(x,y)) +# define IfThen(x,y,z) ((x) ? (y) : (z)) #else -#define Abs(x) (abs(x)) -#define Min(x, y) (min((x), (y))) -#define Max(x, y) (max((x), (y))) -/* IfThen cannot be expressed in Fortran */ +# define Abs(x) (abs(x)) +# define Min(x, y) (min(x,y)) +# define Max(x, y) (max(x,y)) +# define IfThen(x,y,z) ((x)*(y) + (1-(x))*(z)) #endif #ifdef KRANC_C @@ -42,11 +42,52 @@ #define Tanh(x) (tanh(x)) #ifdef KRANC_C -#define E M_E -#define Pi M_PI +# define Sign(x) (copysign(1.0,(x))) +# define ToReal(x) ((CCTK_REAL)(x)) #else -#define E 2.71828182845904523536029d0 -#define Pi 3.14159265358979323846264d0 +# define Sign(x) (sgn(x)) +# define ToReal(x) (real((x),kind(khalf))) #endif -#define StepFunction(x) ( (x) > 0 ? 1 : 0 ) +#if 0 + +/* TODO: use fma(x,y,z) to implement fmadd and friends? Note that fma + may be unsupported, or may be slow. */ + +/* #define fmadd(x,y,z) ((x)*(y)+(z)) */ +/* #define fmsub(x,y,z) ((x)*(y)-(z)) */ +/* #define fnmadd(x,y,z) (-(z)-(x)*(y)) */ +/* #define fnmsub(x,y,z) (+(z)-(x)*(y)) */ + +#define fpos(x) (+(x)) +#define fneg(x) (-(x)) +#define fmul(x,y) ((x)*(y)) +#define fdiv(x,y) ((x)/(y)) +#define fadd(x,y) ((x)+(y)) +#define fsub(x,y) ((x)-(y)) + +#define fmadd(x,y,z) (fadd(fmul(x,y),z)) +#define fmsub(x,y,z) (fsub(fmul(x,y),z)) +#define fnmadd(x,y,z) (fsub(fneg(z),fmul(x,y))) +#define fnmsub(x,y,z) (fsub(z,fmul(x,y))) + +#define kexp(x) (exp(x)) +#define kfabs(x) (fabs(x)) +#define kfmax(x,y) (fmax(x,y)) +#define kfmin(x,y) (fmin(x,y)) +#define klog(x) (log(x)) +#define kpow(x,y) (pow(x,y)) +#define ksqrt(x) (sqrt(x)) + +#endif + +#ifdef KRANC_C +# define E M_E +# define Pi M_PI + +#else +# define E 2.71828182845904523536029d0 +# define Pi 3.14159265358979323846264d0 +#endif + +#define StepFunction(x) ((x)>0) diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/ParamCheck.F90 b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/ParamCheck.F90 deleted file mode 100644 index c2baf82..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/ParamCheck.F90 +++ /dev/null @@ -1,96 +0,0 @@ -/*@@ - @file GenericFD/src/ParamCheck.F90 - @date October 20 2004 - @author S. Husa - @desc - Check consistency of parameters associated with stencil widths - - $Header$ - - @enddesc - @@*/ - -/* Copyright 2004 Sascha Husa, Ian Hinder, Christiane Lechner - - This file is part of Kranc. - - Kranc is free software; you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation; either version 2 of the License, or - (at your option) any later version. - - Kranc is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with Kranc; if not, write to the Free Software - Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -*/ - -#define KRANC_FORTRAN - -#include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Parameters.h" -#include "GenericFD.h" - -/*@@ - @routine GenericFD_ParamCheck - @date October 20 2004 - @author S. Husa - @desc - Check consistency of parameters associated with stencil widths - - @enddesc - @calls - @calledby - @history - @endhistory - @@*/ - - - -SUBROUTINE GenericFD_ParamCheck(CCTK_ARGUMENTS) - -implicit none - -DECLARE_CCTK_ARGUMENTS -DECLARE_CCTK_PARAMETERS - -if (stencil_width < 0) then - call CCTK_WARN(0, "stencil_width < 0 - set GenericFD::stencil_width > 0 in par file!") -endif - -if ((stencil_width_x < 0).AND.(stencil_width < 0)) then - call CCTK_WARN(0, "stencil_width_x < 0, set GenericFD::stencil_width_x (or stencil_width) > 0!") -endif - -if ((stencil_width_y < 0).AND.(stencil_width < 0)) then - call CCTK_WARN(0, "stencil_width_y < 0, set GenericFD::stencil_width_x (or stencil_width) > 0!") -endif - -if ((stencil_width_z < 0).AND.(stencil_width < 0)) then - call CCTK_WARN(0, "stencil_width_z < 0, set GenericFD::stencil_width_x (or stencil_width) > 0!") -endif - - -if (stencil_width > maxval(cctk_nghostzones)) then - call CCTK_WARN(0, "stencil_width is larger than max(cctk_nghostzones)!") -endif - -if (stencil_width_x > cctk_nghostzones(1)) then - call CCTK_WARN(0, "stencil_width is smaller than cctk_nghostzones(1)!") -endif - -if (stencil_width_y > cctk_nghostzones(2)) then - call CCTK_WARN(0, "stencil_width is smaller than cctk_nghostzones(2)!") -endif - -if (stencil_width_z > cctk_nghostzones(3)) then - call CCTK_WARN(0, "stencil_width is smaller than cctk_nghostzones(3)!") -endif - - -END SUBROUTINE GenericFD_ParamCheck diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/Startup.c b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/Startup.c deleted file mode 100644 index 78e9915..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/Startup.c +++ /dev/null @@ -1,60 +0,0 @@ -/*@@ - @file GenericFD/src/Startup.c - @date June 16 2002 - @author S. Husa - @desc - Register Banner - straight copy of WaveToy - - $Header$ - - @enddesc - @@*/ - -/* Copyright 2004 Sascha Husa, Ian Hinder, Christiane Lechner - - This file is part of Kranc. - - Kranc is free software; you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation; either version 2 of the License, or - (at your option) any later version. - - Kranc is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with Kranc; if not, write to the Free Software - Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -*/ - -#include "cctk.h" -#include "GenericFD.h" - -int GenericFD_Startup(void); - -/*@@ - @routine GenericFD_Startup - @date June 16 2002 - @author S. Husa - @desc - - @enddesc - @calls - @calledby - @history - @endhistory - @@*/ - - - -int GenericFD_Startup(void) -{ - const char *banner = "GenericFD: generic finite differencing"; - CCTK_RegisterBanner(banner); - - - CCTK_INFO(FD_METHOD_DESC); - return 0; - } diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/make.code.defn b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/make.code.defn index fde3275..c9ef2f1 100644 --- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/make.code.defn +++ b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/make.code.defn @@ -2,8 +2,7 @@ # $Header$ # Source files in this directory -SRCS = Startup.c GenericFD.c -#ParamCheck.F90 +SRCS = GenericFD.c # Subdirectories containing source files SUBDIRS = diff --git a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/testmacros.c b/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/testmacros.c deleted file mode 100644 index 3208ea0..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/testmacros.c +++ /dev/null @@ -1,23 +0,0 @@ -/* process me with something like - * cpp -DFD_C2 testmacros.c - */ - - -#include "GenericFD.h" - - -dummy = D11loc(testfunc,i,j,k); - -dummy = D11gf(testfunc,i,j,k); - -dummy = D11_c2c4(testfunc,i,j,k); - -dummy = D11(testfunc,i,j,k); - -dummy = D31(testfunc,i,j,k); - -dummy = D1(testfunc,i,j,k); - - -DECLAREVARS - diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/README b/Auxiliary/Cactus/KrancNumericalTools/Perturb/README deleted file mode 100644 index 899015f..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/README +++ /dev/null @@ -1,5 +0,0 @@ - -Purpose of the thorn: - -This thorn adds a perturbation to grid variables at initial data and/or -boundaries, in the style of Denis' noise thorn. diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/interface.ccl b/Auxiliary/Cactus/KrancNumericalTools/Perturb/interface.ccl deleted file mode 100644 index a745e5c..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/interface.ccl +++ /dev/null @@ -1,7 +0,0 @@ -# Interface definition for thorn Perturb -# $Header$ - -IMPLEMENTS: perturb - -INHERITS: grid - diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/param.ccl b/Auxiliary/Cactus/KrancNumericalTools/Perturb/param.ccl deleted file mode 100644 index 8e87dc6..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/param.ccl +++ /dev/null @@ -1,45 +0,0 @@ -# Parameter definitions for thorn Perturb -# $Header$ - -#------------------------------------------------------------------------------ -# Private: -#------------------------------------------------------------------------------ -private: - -BOOLEAN apply_id_perturb "Add perturbation to initial data" -{ -} "no" - -BOOLEAN apply_bc_perturb "Add perturbation to boundary data" -{ -} "no" - -STRING id_vars "Initial data variables to modify with perturbation" -{ - .* :: "A regex which matches everything" -} "" - -STRING bc_vars "Variables to modify with perturbation at boundary" -{ - .* :: "A regex which matches everything" -} "" - -BOOLEAN perturb_boundaries[6] "At which boundaries to apply perturbation" -{ -} "yes" - -INT perturb_stencil[3] "Number of boundary points" -{ - 0:* :: "0:*" -} 1 - -REAL amplitude "Amplitude of perturbation data" -{ - 0: :: "Positive number" -} 0.000001 - -REAL period "period of perturbation data" -{ - 0: :: "Positive number" -} 1.0 - diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/schedule.ccl b/Auxiliary/Cactus/KrancNumericalTools/Perturb/schedule.ccl deleted file mode 100644 index 268f81d..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/schedule.ccl +++ /dev/null @@ -1,21 +0,0 @@ -# Schedule definitions for thorn Perturb -# $Header$ - -if (apply_id_perturb) { - SCHEDULE id_perturb AT CCTK_POSTINITIAL - { - LANG: C - } "Add perturb to initial data" -} - -if (apply_bc_perturb) { - SCHEDULE bc_perturb AT CCTK_POSTSTEP - { - LANG: C - } "Add perturb to boundary condition" - - SCHEDULE bc_perturb AT CCTK_POSTRESTRICT - { - LANG: C - } "Add perturb to boundary condition" -} diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/bc_perturb.c b/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/bc_perturb.c deleted file mode 100644 index bf56039..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/bc_perturb.c +++ /dev/null @@ -1,322 +0,0 @@ -/* $Header$ */ - -#include <assert.h> -#include <stdio.h> -#include <stdlib.h> -#include <string.h> - - - -#include "perturb.h" - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "cctk_Arguments.h" -#include "cctk_FortranString.h" - -#include "Symmetry.h" - - - -/* #define DEBUG_BOUNDARY 1 */ - -static int ApplyBndperturb (const cGH *GH, - int stencil_dir, - const int *stencil, - int dir, - int first_var, - int num_vars); - - - -int -BndperturbVI (const cGH *GH, const int *stencil, int vi) -{ - int retval; - retval = ApplyBndperturb (GH, -1, stencil, 0, vi, 1); - return retval; -} - -void -CCTK_FCALL CCTK_FNAME (BndperturbVI) (int *ierr, const cGH **GH, - const int *stencil, const int *vi) -{ - *ierr = BndperturbVI (*GH, stencil, *vi); -} - -int -BndperturbVN (const cGH *GH, const int *stencil, const char *vn) -{ - int vi, retval; - vi = CCTK_VarIndex(vn); - retval = BndperturbVI (GH, stencil, vi); - return retval; -} - -void -CCTK_FCALL CCTK_FNAME (BndperturbVN) (int *ierr, const cGH **GH, - const int *stencil, ONE_FORTSTRING_ARG) -{ - ONE_FORTSTRING_CREATE (vn); - *ierr = BndperturbVN (*GH, stencil, vn); - free (vn); -} - - - -int -BndperturbGI (const cGH *GH, const int *stencil, int gi) -{ - int first_vi, retval; - - first_vi = CCTK_FirstVarIndexI (gi); - if (first_vi >= 0) - { - retval = ApplyBndperturb (GH, -1, stencil, 0, first_vi, - CCTK_NumVarsInGroupI (gi)); - } - else - { - CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, - "Invalid group index %d in BndperturbGI", gi); - retval = -1; - } - - return (retval); -} - - - -void -CCTK_FCALL CCTK_FNAME (BndperturbGI) (int *ierr, const cGH **GH, - const int *stencil, const int *gi) -{ - *ierr = BndperturbGI (*GH, stencil, *gi); -} - - -int -BndperturbGN (const cGH *GH, const int *stencil, const char *gn) -{ - int gi, retval; - - gi = CCTK_GroupIndex (gn); - if (gi >= 0) - { - retval = BndperturbGI (GH, stencil, gi); - } - else - { - CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, - "Invalid group name '%s' in BndperturbGN", gn); - retval = -1; - } - - return (retval); -} - -void CCTK_FCALL CCTK_FNAME (BndperturbGN) - (int *ierr, - const cGH **GH, - const int *stencil, - ONE_FORTSTRING_ARG) -{ - ONE_FORTSTRING_CREATE (gn) - *ierr = BndperturbGN (*GH, stencil, gn); - free (gn); -} - - - - - -static int ApplyBndperturb (const cGH *GH, - int stencil_dir, - const int *stencil, - int dir, - int first_var, - int num_vars) -{ - DECLARE_CCTK_PARAMETERS; - int i, j, k; - int var, vtypesize, gindex, gdim, timelvl; - int doBC[2*MAXDIM], dstag[MAXDIM], lsh[MAXDIM], lssh[MAXDIM]; - SymmetryGHex *sGHex; - int type; - - /* This argument is unused an undocumented; better make sure people - don't try to use it for something. */ - assert (stencil_dir == -1); - - /* get the group index of the variables */ - gindex = CCTK_GroupIndexFromVarI (first_var); - - /* get the number of dimensions and the size of the variables' type */ - gdim = CCTK_GroupDimI (gindex); - vtypesize = CCTK_VarTypeSize (CCTK_VarTypeI (first_var)); - - /* make sure we can deal with this number of dimensions */ - if (gdim > MAXDIM) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "ApplyBndperturb: Variable dimension of %d not supported", gdim); - return (-1); - } - - /* check the direction parameter */ - if (abs (dir) > gdim) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "ApplyBndperturb: direction %d greater than dimension %d", - dir, gdim); - return (-2); - } - - /* initialize arrays for variables with less dimensions than MAXDIM - so that we can use the INDEX_3D macro later on */ - for (i = gdim; i < MAXDIM; i++) - { - lsh[i] = 1; - lssh[i] = 1; - } - - /* get the directional staggering of the group */ - CCTK_GroupStaggerDirArrayGI (dstag, gdim, gindex); - - /* get the current timelevel */ - timelvl = 0; - - /* see if we have a symmetry array */ - sGHex = (SymmetryGHex *) CCTK_GHExtension (GH, "Symmetry"); - - - /* now loop over all variables */ - for (var = first_var; var < first_var + num_vars; var++) - { - /* Apply condition if: - + boundary is not a symmetry boundary (no symmetry or unset(=unsed)) - + boundary is a physical boundary - + have enough grid points - */ - memset (doBC, 1, sizeof (doBC)); - if (sGHex) - { - for (i = 0; i < 2 * gdim; i++) - { - doBC[i] = sGHex->GFSym[var][i] == GFSYM_NOSYM || - sGHex->GFSym[var][i] == GFSYM_UNSET; - } - } - for (i = 0; i < gdim; i++) - { - lsh[i] = GH->cctk_lsh[i]; - lssh[i] = GH->cctk_lssh[CCTK_LSSH_IDX (dstag[i], i)]; - doBC[i*2] &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2]; - doBC[i*2+1] &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2+1]; - if (dir != 0) - { - doBC[i*2] &= i+1 == -dir; - doBC[i*2+1] &= i+1 == dir; - } - } - - /* now apply the boundaries face by face */ - if (gdim > 0) - { -#ifdef DEBUG_BOUNDARY - if (doBC[0]) - { - printf("Boundary: Applying lower x perturb to boundary\n"); - } - if (doBC[1]) - { - printf("Boundary: Applying upper x perturb to boundary\n"); - } -#endif /* DEBUG_BOUNDARY */ - /* lower x */ - BOUNDARY_PERTURB (doBC[0], stencil[0], lssh[1], lssh[2], - i, j, k); - /* upper x */ - BOUNDARY_PERTURB (doBC[1], stencil[0], lssh[1], lssh[2], - lssh[0]-i-1, j, k); - - } - if (gdim > 1) - - { -#ifdef DEBUG_BOUNDARY - if (doBC[2]) - { - printf("Boundary: Applying lower y perturb to boundary\n"); - } - if (doBC[3]) - { - printf("Boundary: Applying upper y perturb to boundary\n"); - } -#endif /* DEBUG_BOUNDARY */ - /* lower y */ - BOUNDARY_PERTURB (doBC[2], lssh[0], stencil[1], lssh[2], - i, j, k); - /* upper y */ - BOUNDARY_PERTURB (doBC[3], lssh[0], stencil[1], lssh[2], - i, lssh[1]-j-1, k); - } - if (gdim > 2) - { -#ifdef DEBUG_BOUNDARY - if (doBC[4]) - { - printf("Boundary: Applying lower z perturb to boundary\n"); - } - if (doBC[5]) - { - printf("Boundary: Applying upper z perturb to boundary\n"); - } -#endif /* DEBUG_BOUNDARY */ - /* lower z */ - BOUNDARY_PERTURB (doBC[4], lssh[0], lssh[1], stencil[2], - i, j, k); - /* upper z */ - BOUNDARY_PERTURB (doBC[5], lssh[0], lssh[1], stencil[2], - i, j, lssh[2]-k-1); - } - } - - return(0); -} - - -static void -add_bc_perturb_to_var (int idx, const char* optstring, void* cctkGH) -{ - DECLARE_CCTK_PARAMETERS; - cGH* GH = cctkGH; - int sw[3]; - - /* Change type from CCTK_INT to int */ - sw[0] = perturb_stencil[0]; - sw[1] = perturb_stencil[1]; - sw[2] = perturb_stencil[2]; - - if (perturb_boundaries[0]) ApplyBndperturb (GH, -1, sw, -1, idx, 1); - if (perturb_boundaries[1]) ApplyBndperturb (GH, -1, sw, +1, idx, 1); - if (perturb_boundaries[2]) ApplyBndperturb (GH, -1, sw, -2, idx, 1); - if (perturb_boundaries[3]) ApplyBndperturb (GH, -1, sw, +2, idx, 1); - if (perturb_boundaries[4]) ApplyBndperturb (GH, -1, sw, -3, idx, 1); - if (perturb_boundaries[5]) ApplyBndperturb (GH, -1, sw, +3, idx, 1); -} - - -void -bc_perturb(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - -/* Boundary_MakeSureThatTheSelectionIsEmpty(); */ - if (CCTK_TraverseString(bc_vars, add_bc_perturb_to_var, cctkGH, - CCTK_GROUP_OR_VAR) < 0) - { - CCTK_WARN (1, "Failed to parse 'perturb::bc_vars' parameter"); - } -} diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/id_perturb.c b/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/id_perturb.c deleted file mode 100644 index 07bd106..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/id_perturb.c +++ /dev/null @@ -1,49 +0,0 @@ -/* - $Header$ -*/ - -#include <stdlib.h> - -#include "perturb.h" - -#include "cctk_Arguments.h" -#include "cctk_Parameters.h" - -void -add_perturb_to_var (int idx, const char* optstring, void* cctkGH) -{ - DECLARE_CCTK_PARAMETERS; - int i, j, k, ijk; - CCTK_REAL* data; - cGH* GH = cctkGH; - int type; - - data = (CCTK_REAL*) CCTK_VarDataPtrI(GH, 0, idx); - - for (k=1; k< GH->cctk_lsh[2]-1; ++k) - { - for (j=1; j< GH->cctk_lsh[1]-1; ++j) - { - for (i=1; i< GH->cctk_lsh[0]-1; ++i) - { - ijk = CCTK_GFINDEX3D(GH, i, j, k); - - data[ijk] += amplitude * pow(-1, i/period) * pow(-1, j/period) * pow(-1, k/period); - } - } - } -} - - -void -id_perturb(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - if (CCTK_TraverseString(id_vars, add_perturb_to_var, cctkGH, - CCTK_GROUP_OR_VAR) < 0) - { - CCTK_WARN (1, "Failed to parse 'IDRandom::noisy_id_vars' parameter"); - } -} diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/make.code.defn b/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/make.code.defn deleted file mode 100644 index 795dd51..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/make.code.defn +++ /dev/null @@ -1,8 +0,0 @@ -# Main make.code.defn file for thorn perturb -# $Header$ - -# Source files in this directory -SRCS = id_perturb.c bc_perturb.c - -# Subdirectories containing source files -SUBDIRS = diff --git a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/perturb.h b/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/perturb.h deleted file mode 100644 index e7da060..0000000 --- a/Auxiliary/Cactus/KrancNumericalTools/Perturb/src/perturb.h +++ /dev/null @@ -1,45 +0,0 @@ -/* $Header$ */ - -#ifndef PERTURB_H -#define PERTURB_H - -#include <assert.h> -#include <stdlib.h> -#include <string.h> -#include <time.h> -#include "cctk.h" - - -/* constants */ - -#define MAXDIM 3 - - -/* macros */ - -#define RAND_VAL ((random()*(1.0/RAND_MAX)-0.5)*amplitude) - -#define BOUNDARY_PERTURB(doBC, \ - iend, jend, kend, \ - ii, jj, kk) \ -{ \ - if (doBC) \ - { \ - CCTK_REAL* v= CCTK_VarDataPtrI(GH, timelvl, var); \ - for (k = 0; k < kend; k++) \ - { \ - for (j = 0; j < jend; j++) \ - { \ - for (i = 0; i < iend; i++) \ - { \ - const int _index = CCTK_GFINDEX3D(GH, (ii), (jj), (kk)); \ - v[_index] += RAND_VAL; \ - } \ - } \ - } \ - } \ -} - - - -#endif /* !define(PERTURB_H) */ diff --git a/Bin/krancvars.sh b/Bin/krancvars.sh new file mode 100755 index 0000000..0c91aae --- /dev/null +++ b/Bin/krancvars.sh @@ -0,0 +1,8 @@ + +export KRANCDIR=$(ls -d $(dirname ${BASH_ARGV[0]})/..) +if ! expr "$KRANCDIR">/dev/null : '/'; then + # Relative path + KRANCDIR="$PWD/$KRANCDIR" +fi + +export PATH="$KRANCDIR/Bin:$PATH" diff --git a/Bin/thorn-branch b/Bin/thorn-branch new file mode 100755 index 0000000..e0f9520 --- /dev/null +++ b/Bin/thorn-branch @@ -0,0 +1,82 @@ +#!/bin/bash + +function usage() +{ + cat <<EOF +Usage: thorn-branch <command> + +Available commands are: + + init Set up a new <branch>_thorns branch in the current repository. + Currently only Git repositories are supported, and there + must be at least one commit in the repository. + + push Clone a temporary copy of the current branch, run "make" in + it, and commit the contents of the resulting thorns + directory to the <branch>_thorns branch. Push this branch + and the current branch to origin. + +EOF + +} + +set -e + +if [ $# = 0 ]; then + usage + exit 1 +fi + +command="$1" + +case "$command" in + "help") + usage + ;; + + "push") + BRANCH=$(git branch|grep '^\*'|awk '{print $2}') + if [ -r ${BRANCH}_thorns ]; then + echo "Existing directory ${BRANCH}_thorns: please remove before running thorn-branch init" + exit 1 + fi + if [ -r ${BRANCH} ]; then + echo "Existing directory ${BRANCH}: please remove before running thorn-branch init" + exit 1 + fi + git fetch origin + git branch -f ${BRANCH}_thorns remotes/origin/${BRANCH}_thorns + git clone -b ${BRANCH} . ${BRANCH} + make -C ${BRANCH} + git clone -b ${BRANCH}_thorns . ${BRANCH}_thorns + rsync -av --exclude=".git/*" --delete ${BRANCH}/thorns/ ${BRANCH}_thorns/ + git --git-dir ${PWD}/${BRANCH}_thorns/.git --work-tree ${PWD}/${BRANCH}_thorns add --all + if git --git-dir ${PWD}/${BRANCH}_thorns/.git --work-tree ${PWD}/${BRANCH}_thorns commit -a -m "Updated from source"; then + true + fi + git fetch ${BRANCH}_thorns "${BRANCH}_thorns:${BRANCH}_thorns" + git push origin ${BRANCH} && git push origin ${BRANCH}_thorns + rm -rf ${BRANCH} ${BRANCH}_thorns + ;; + + "init") + BRANCH=$(git branch|grep '^\*'|awk '{print $2}') + if [ -r ${BRANCH}_thorns ]; then + echo "Existing directory ${BRANCH}_thorns: please remove before running thorn-branch init" + exit 1 + fi + git clone -b ${BRANCH} . ${BRANCH}_thorns + git --git-dir ${PWD}/${BRANCH}_thorns/.git --work-tree ${PWD}/${BRANCH}_thorns symbolic-ref HEAD refs/heads/master_thorns + rm ${BRANCH}_thorns/.git/index + git --git-dir ${PWD}/${BRANCH}_thorns/.git --work-tree ${PWD}/${BRANCH}_thorns clean -fdxq + echo "This is the initial commit" > ${PWD}/${BRANCH}_thorns/README + git --git-dir ${PWD}/${BRANCH}_thorns/.git --work-tree ${PWD}/${BRANCH}_thorns add ${PWD}/${BRANCH}_thorns/README + git --git-dir ${PWD}/${BRANCH}_thorns/.git --work-tree ${PWD}/${BRANCH}_thorns commit -a -m "Initial thorns commit" + git fetch ${BRANCH}_thorns "${BRANCH}_thorns:${BRANCH}_thorns" + rm -rf ${BRANCH}_thorns + ;; + + *) + echo "thorn-branch: Command \'$command\' not recognised" + exit 1 +esac diff --git a/Doc/KrancDoc.tex b/Doc/KrancDoc.tex index d6c9262..c6e1821 100644 --- a/Doc/KrancDoc.tex +++ b/Doc/KrancDoc.tex @@ -3,6 +3,7 @@ \usepackage{tabularx} \usepackage{graphicx} \usepackage{alltt} +\usepackage{hyperref} \addtolength{\oddsidemargin}{-0.25in} \addtolength{\textwidth}{1in} @@ -138,111 +139,112 @@ by the Kranc system, such as a name for the calculation. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\chapter{Getting Started} +%% \chapter{Getting Started} -\section{Tutorial for new Cactus users} +%% \section{Tutorial for new Cactus users} -This section provides a step-by-step tutorial for Kranc for people who -have not used Cactus before. It will lead you through downloading -Cactus and Kranc, building one of the example Kranc thorns, compiling -the thorn into a Cactus executable, running the executable and -analysing the output. People who have used Cactus before and already -have a working setup should skip this section and use -Sec.~\ref{sec:tutcur} instead. +%% This section provides a step-by-step tutorial for Kranc for people who +%% have not used Cactus before. It will lead you through downloading +%% Cactus and Kranc, building one of the example Kranc thorns, compiling +%% the thorn into a Cactus executable, running the executable and +%% analysing the output. People who have used Cactus before and already +%% have a working setup should skip this section and use +%% Sec.~\ref{sec:tutcur} instead. -The first step is to create a Cactus directory. Everything will then -be downloaded into this directory and compiled there. +%% The first step is to create a Cactus directory. Everything will then +%% be downloaded into this directory and compiled there. -\begin{verbatim} - mkdir Cactus - cd Cactus -\end{verbatim} +%% \begin{verbatim} +%% mkdir Cactus +%% cd Cactus +%% \end{verbatim} -Download Kranc into the newly created Cactus directory: -\begin{verbatim} - git clone git://github.com/ianhinder/kranc -\end{verbatim} +%% Download Kranc into the newly created Cactus directory: +%% \begin{verbatim} +%% git clone git://github.com/ianhinder/kranc +%% \end{verbatim} -Download the Cactus {\em GetComponents} script and make it executable. -This script allows you to download Cactus thorns from many different -repositories listed in a ``ThornList''. -\begin{verbatim} - wget http://www.cactuscode.org/download/GetComponents - chmod u+x GetComponents -\end{verbatim} +%% Download the Cactus {\em GetComponents} script and make it executable. +%% This script allows you to download Cactus thorns from many different +%% repositories listed in a ``ThornList''. +%% \begin{verbatim} +%% wget http://www.cactuscode.org/download/GetComponents +%% chmod u+x GetComponents +%% \end{verbatim} -Download Cactus and several thorns specified in the Kranc example thornlist: -\begin{verbatim} - ./GetComponents --root=. kranc/Auxiliary/Cactus/thornlist.th -\end{verbatim} +%% Download Cactus and several thorns specified in the Kranc example thornlist: +%% \begin{verbatim} +%% ./GetComponents --root=. kranc/Auxiliary/Cactus/thornlist.th +%% \end{verbatim} -Configure Cactus to run on your local machine. You will need a C -compiler. See the Cactus documentation for more information on this -step if Cactus does not configure automatically on your system. -\begin{verbatim} - make wave-config ThornList=kranc/Auxiliary/Cactus/thornlist.th -\end{verbatim} +%% Configure Cactus to run on your local machine. You will need a C +%% compiler. See the Cactus documentation for more information on this +%% step if Cactus does not configure automatically on your system. +%% \begin{verbatim} +%% make wave-config ThornList=kranc/Auxiliary/Cactus/thornlist.th +%% \end{verbatim} -Build a Cactus executable: -\begin{verbatim} - make wave -\end{verbatim} +%% Build a Cactus executable: +%% \begin{verbatim} +%% make wave +%% \end{verbatim} -Run the Cactus executable with the example wave equation parameter file: -\begin{verbatim} - cd Kranc/Examples - mpirun ../../exe/cactus_wave wave.par -\end{verbatim} +%% Run the Cactus executable with the example wave equation parameter file: +%% \begin{verbatim} +%% cd Kranc/Examples +%% mpirun ../../exe/cactus_wave wave.par +%% \end{verbatim} -Analyse the output of the simulation: -\begin{verbatim} - ygraph wave/phi.x.asc -\end{verbatim} +%% Analyse the output of the simulation: +%% \begin{verbatim} +%% ygraph wave/phi.x.asc +%% \end{verbatim} -\section{Tutorial for Current Cactus Users} -\label{sec:tutcur} +%% \section{Tutorial for Current Cactus Users} +%% \label{sec:tutcur} -Change into your Cactus directory and download the current version of Kranc: -\begin{verbatim} - cd Cactus - git clone git://github.com/ianhinder/kranc -\end{verbatim} +%% Change into your Cactus directory and download the current version of Kranc: +%% \begin{verbatim} +%% cd Cactus +%% git clone git://github.com/ianhinder/kranc +%% \end{verbatim} -Make a symbolic link from the Kranc examples to the arrangements directory: -\begin{verbatim} - ln -s ../kranc/Examples arrangements/KrancExamples -\end{verbatim} +%% Make a symbolic link from the Kranc examples to the arrangements directory: +%% \begin{verbatim} +%% ln -s ../kranc/Examples arrangements/KrancExamples +%% \end{verbatim} -Make a symbolic link from the Kranc support arrangement into your Cactus tree: -\begin{verbatim} - ln -s ../kranc/Auxiliary/Cactus/KrancNumericalTools arrangements/KrancNumericalTools -\end{verbatim} +%% Make a symbolic link from the Kranc support arrangement into your Cactus tree: +%% \begin{verbatim} +%% ln -s ../kranc/Auxiliary/Cactus/KrancNumericalTools arrangements/KrancNumericalTools +%% \end{verbatim} -Go into the KrancExamples arrangement and build the Wave thorn: -\begin{verbatim} - cd arrangements/KrancExamples - ../../kranc/Bin/kranc Wave.m -\end{verbatim} +%% Go into the KrancExamples arrangement and build the Wave thorn: +%% \begin{verbatim} +%% cd arrangements/KrancExamples +%% ../../kranc/Bin/kranc Wave.m +%% \end{verbatim} -It is assumed that Mathematica is available on your path with the name - "math", or at the location - \verb|/Applications/Mathematica.app/Contents/MacOS/MathKernel| for Mac OS. +%% It is assumed that Mathematica is available on your path with the name +%% "math", or at the location +%% \verb|/Applications/Mathematica.app/Contents/MacOS/MathKernel| for Mac OS. -You should get a thorn generated in the current directory called Wave. +%% You should get a thorn generated in the current directory called Wave. -Add \verb|KrancExamples/Wave| and \verb|KrancNumericalTools/GenericFD| to the - thornlist of one of your configurations, or make one from scratch. +%% Add \verb|KrancExamples/Wave| and \verb|KrancNumericalTools/GenericFD| to the +%% thornlist of one of your configurations, or make one from scratch. -Recompile Cactus. +%% Recompile Cactus. -Run Cactus with the parameter file wave.par found in - \verb|arrangements/KrancExamples|. +%% Run Cactus with the parameter file wave.par found in +%% \verb|arrangements/KrancExamples|. -Examine the 1D output phi.x.asc and pi.x.asc. +%% Examine the 1D output phi.x.asc and pi.x.asc. -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Using Kranc} +\label{chp:usingkranc} %% \section{Types of arguments} @@ -925,6 +927,111 @@ same arguments, but they can be tensorial in nature. \end{tabularx} \end{center} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\chapter{Additional Features} + +In addition to the basic functionality described in Chapter +\ref{chp:usingkranc}, Kranc provides a number of additional features. +Each of these features can be used independently, and are typically +enabled by an option of the form \verb|Use* -> True| in the call to +\verb|CreateThorn|. + +\section{Jacobians} + +Kranc allows the user to apply an arbitrary user-defined Jacobian +transformation (provided in a grid function) to all derivatives. This +feature is enabled by setting \verb|UseJacobian -> True| in the +options to \verb|CreateThorn| or \verb|CreateThornTT|. This can be +used to generate codes which work with multi-block schemes. The use +of the Jacobian is determined at run time, so a single code can be +generated which will work both with and without a Jacobian being +present. + +Kranc does not provide the Jacobian grid function; it might be +provided by an external infrastructure (for example for multi-block +schemes), or could be provided easily in the user's thorn, or another +Kranc-generated thorn. Wherever the Jacobian is provided, it must +adhere to the following conventions. There should be one Cactus group +for the components of the Jacobian matrix and another for the +components of its first spatial derivative (this is necessary for +systems containing second spatial derivatives). The components should +be real-valued grid functions declared in a similar manner to the +following: + +\begin{verbatim} +CCTK_REAL jac type=GF timelevels=1 +{ + J11, J12, J13, J21, J22, J23, J31, J32, J33 +} "Jacobian of the coordinate transformation" + +CCTK_REALd djac type=GF timelevels=1 +{ + dJ111, dJ112, dJ113, dJ122, dJ123, dJ133, + dJ211, dJ212, dJ213, dJ222, dJ223, dJ233, + dJ311, dJ312, dJ313, dJ322, dJ323, dJ333, +} "Derivative of the Jacobian" +\end{verbatim} + +The names of the groups and variables are not important, but the order +of the variables within the groups is critical. + +The GenericFD thorn provides two parameters, \verb|jacobian_group| and +\verb|jacobian_derivative_group| which should be set by the user in their +parameter file to the names of the Jacobian and Jacobian derivative +groups. With the above Jacobian definition, provided by a thorn with +implementation \verb|MyCoordTransform| the user would set + +\begin{verbatim} +GenericFD::jacobian = "MyCoordTransform::jac" +GenericFD::jacobian_derivative = "MyCoordTransform::djac" +\end{verbatim} + +The partial derivatives, associated with certain finite difference +operators, specified in the user's calculation, will then be +multiplied by the Jacobian. If the user specifies the following, + +\begin{center} +\begin{minipage}{0.8 \textwidth} +\begin{verbatim} +derivs = { + PDstandard2nd[i_] -> DZero[i]} + +... + +dot[v] -> PDstandard2nd[v,1] +\end{verbatim} +\end{minipage} +\end{center} + +the code that will actually be generated will be + +\begin{center} +\begin{minipage}{0.8 \textwidth} +\begin{verbatim} +dot[v] -> J11 PDstandard2nd[v,1] + J21 PDstandard2nd[v,2] + + J31 PDstandard2nd[v,3] +\end{verbatim} +\end{minipage} +\end{center} + +Note: + +\begin{itemize} +\item The Jacobian multiplication introduces an additional performance + cost to the simulation, so it should not be enabled unless + necessary. +\item If the Jacobian parameters are left unset, the partial + derivatives will not be multiplied by any Jacobian, and the + performance will be only slightly worse than if UseJacobian was not + set and the Jacobian code was not generated (the extra cost comes + from the use of run-time conditionals). +\item The parameter \verb|GenericFD::jacobian_identity_map| can be + used with multi-block systems to prevent the use of the Jacobian on + one of the maps. Set this parameter to the map number (e.g.~0) + where the Jacobian should not be applied. This can lead to a + significant performance improvement. +\end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/Examples/Analysis.nb b/Examples/Analysis.nb new file mode 100644 index 0000000..2aa39a2 --- /dev/null +++ b/Examples/Analysis.nb @@ -0,0 +1,61 @@ +(* Content-type: application/mathematica *) + +(*** Wolfram Notebook File ***) +(* http://www.wolfram.com/nb *) + +(* CreatedBy='Mathematica 7.0' *) + +(*CacheID: 234*) +(* Internal cache information: +NotebookFileLineBreakTest +NotebookFileLineBreakTest +NotebookDataPosition[ 145, 7] +NotebookDataLength[ 1395, 52] +NotebookOptionsPosition[ 1075, 37] +NotebookOutlinePosition[ 1415, 52] +CellTagsIndexPosition[ 1372, 49] +WindowFrame->Normal*) + +(* Beginning of Notebook Content *) +Notebook[{ +Cell[BoxData[ + RowBox[{ + RowBox[{"SetDirectory", "[", + RowBox[{"NotebookDirectory", "[", "]"}], "]"}], ";"}]], "Input", + CellChangeTimes->{{3.4951847946945972`*^9, 3.4951848185680017`*^9}}], + +Cell[BoxData[ + RowBox[{"Get", "[", "\"\<../Auxiliary/Cactus/IOASCII.m\>\"", "]"}]], "Input",\ + + CellChangeTimes->{{3.49518480029657*^9, 3.4951848005403*^9}}], + +Cell[BoxData[ + RowBox[{"MGraph", "[", + RowBox[{"ReadIOASCII", "[", "\"\<simplewave_sine/phi.xl\>\"", "]"}], + "]"}]], "Input", + CellChangeTimes->{3.49518480818021*^9}] +}, +WindowSize->{640, 750}, +WindowMargins->{{36, Automatic}, {Automatic, 84}}, +FrontEndVersion->"7.0 for Mac OS X x86 (32-bit) (February 18, 2009)", +StyleDefinitions->"Default.nb" +] +(* End of Notebook Content *) + +(* Internal cache information *) +(*CellTagsOutline +CellTagsIndex->{} +*) +(*CellTagsIndex +CellTagsIndex->{} +*) +(*NotebookFileOutline +Notebook[{ +Cell[545, 20, 193, 4, 27, "Input"], +Cell[741, 26, 157, 3, 27, "Input"], +Cell[901, 31, 170, 4, 27, "Input"] +} +] +*) + +(* End of internal cache information *) diff --git a/Examples/EM-xTensor.m b/Examples/EM-xTensor.m new file mode 100644 index 0000000..19837a5 --- /dev/null +++ b/Examples/EM-xTensor.m @@ -0,0 +1,156 @@ +(* ::Package:: *) + +(* Copyright 2004 Sascha Husa, Ian Hinder, Christiane Lechner, Barry Wardell + + This file is part of Kranc. + + Kranc is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + Kranc is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Foobar; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*) + +$KrancTensorPackage="xTensorKranc`"; +Get["KrancTensor`"]; +Get["Metrics`LoadMetric`"]; +$DefInfoQ=False; +$CVVerbose=False; +$PrePrint=ScreenDollarIndices; + +(**************************************************************************************) +(* Derivatives *) +(**************************************************************************************) + +derivatives = +{ + PDstandard2nd[i_] -> StandardCenteredDifferenceOperator[1,1,i], + PDstandard2nd[i_, i_] -> StandardCenteredDifferenceOperator[2,1,i], + PDstandard2nd[i_, j_] -> StandardCenteredDifferenceOperator[1,1,i] * + StandardCenteredDifferenceOperator[1,1,j], + + PDstandard4th[i_] -> StandardCenteredDifferenceOperator[1,2,i], + PDstandard4th[i_, i_] -> StandardCenteredDifferenceOperator[2,2,i], + PDstandard4th[i_, j_] -> StandardCenteredDifferenceOperator[1,2,i] * + StandardCenteredDifferenceOperator[1,2,j] +}; + +(**************************************************************************************) +(* Tensors *) +(**************************************************************************************) +(* Load Euclidean metric and set components of g, epsilon and gamma *) +LoadMetric["Euclidean"]; +MetricCompute[g, Euclidean, "Christoffel"[1,-1,-1], CVSimplify -> Simplify]; +AllComponentValues[ToBasis[Euclidean][epsilong[-a,e,f]], + Table[Signature[{i,j,k}],{i,3},{j,3},{k,3}]]; +RuleToSet[ChristoffelCDPDEuclidean]; +RuleToSet[g]; +RuleToSet[epsilong]; + +(* Register the tensor quantities with the TensorTools package *) +DefTensor[El[-a],EuclideanM]; +DefTensor[B[-a],EuclideanM]; + +(**************************************************************************************) +(* Groups *) +(**************************************************************************************) + +(* NB This will give B the symmetries of a VECTOR, which is not correct *) +evolvedGroups = Map[CreateGroupFromTensor, ToBasis[Euclidean][{El[-a], B[-a]}]]; +evaluatedGroups = + {{"constraints", {CEl, CB}}, + {"endens",{rho}}}; + +declaredGroups = Join[evolvedGroups, evaluatedGroups]; +declaredGroupNames = Map[First, declaredGroups]; + +groups = Join[declaredGroups]; + +(**************************************************************************************) +(* Initial data *) +(**************************************************************************************) + +initialCalc = +{ + Name -> "EM_initial", + Schedule -> {"at CCTK_INITIAL"}, + Equations -> + { + El1 -> sigma*Cos[2 Pi (x + y)] , + El2 -> - (1 - sigma)*Cos[2 Pi x] - sigma*Cos[2 Pi (x + y)], + El3 -> 0, + B1 -> 0, + B2 -> 0, + B3 -> (1 - sigma)*Cos[2 Pi x] + sigma*Cos[2 Pi (x + y)] + } +}; + +(**************************************************************************************) +(* Evolution equations *) +(**************************************************************************************) + +evolCalc = +{ + Name -> "EM_evol", + Schedule -> {"in MoL_CalcRHS"}, + Equations -> + { + dot[ToBasis[Euclidean][El[-a]]] -> ToBasis[Euclidean][(epsilong[-a,e,f] CD[-e][B[-f]])], + dot[ToBasis[Euclidean][ B[-a]]] -> ToBasis[Euclidean][-(epsilong[-a,e,f] CD[-e][El[-f]])] + } +}; + +(**************************************************************************************) +(* Constraint equations *) +(**************************************************************************************) + +constraintsCalc = +{ + Name -> "EM_constraints", + Equations -> + { + CEl -> ToBasis[Euclidean][g[a,b]CD[-b][El[-a]]], + CB -> ToBasis[Euclidean][ g[a,b]CD[-b][B[-a]]] + } +}; + +(**************************************************************************************) +(* Energy equation *) +(**************************************************************************************) + +energyCalc = +{ + Name -> "EM_energy", + Equations -> + { + rho -> ToBasis[Euclidean][El[-a] El[a]/2 + B[-a] B[a]/2] + } +}; + +realParameters = {sigma}; + +(**************************************************************************************) +(* Construct the thorn *) +(**************************************************************************************) + +calculations = +{ + initialCalc, + evolCalc, + constraintsCalc, + energyCalc +}; + +CreateKrancThornTT[groups, ".", "EM-xTensor", + Calculations -> calculations, + DeclaredGroups -> declaredGroupNames, + PartialDerivatives -> derivatives, + RealParameters -> realParameters]; diff --git a/Examples/Wave.m b/Examples/Wave.m index 0c19475..94dbbda 100644 --- a/Examples/Wave.m +++ b/Examples/Wave.m @@ -36,10 +36,10 @@ zerodims = {}; PD = PDstandard2nd; PDnorm = PDplus; -evolvedGroup = {"evolved", {phi, pi}}; -normGroup = {"norms", {VL2, VDP, EL2}}; -exactGroup = {"exact", {phiExact, piExact}}; -errorGroup = {"errors", {phiError, piError}}; +evolvedGroup = {"evolved", {phi, pi}, Timelevels -> 3}; +normGroup = {"norms", {VL2, VDP, EL2}, Timelevels -> 3}; +exactGroup = {"exact", {phiExact, piExact}, Timelevels -> 3}; +errorGroup = {"errors", {phiError, piError}, Timelevels -> 3}; groups = {evolvedGroup, normGroup, exactGroup, errorGroup}; @@ -105,7 +105,7 @@ importerCalc = evolveEquations = { dot[phi] -> pi, - dot[pi] -> PD[phi,li,ui] + dot[pi] -> PD[phi,li,lj] Euc[ui,uj] } evolveCalc = diff --git a/Examples/simplewave_sine.par b/Examples/simplewave_sine.par index f552f0f..362a8e5 100644 --- a/Examples/simplewave_sine.par +++ b/Examples/simplewave_sine.par @@ -1,22 +1,18 @@ -#!/usr/bin/perl -W ActiveThorns = " Boundary -Carpet -CarpetIOASCII -CarpetIOBasic -CarpetIOScalar -CarpetLib -CarpetReduce -CarpetSlab CartGrid3d CoordBase GenericFD +IOASCII IOUtil -LoopControl +IOBasic MoL NanChecker Periodic +PUGH +PUGHReduce +PUGHSlab SimpleWave Slab SymBase @@ -27,49 +23,32 @@ Time # Grid ############################################################# -CoordBase::domainsize = minmax - -CoordBase::xmin = 0 -CoordBase::ymin = 0 -CoordBase::zmin = 0 - -CoordBase::xmax = 1 -CoordBase::ymax = 0.1 -CoordBase::zmax = 0.1 - -CoordBase::dx = 0.05 -CoordBase::dy = 0.1 -CoordBase::dz = 0.1 - CoordBase::boundary_size_x_lower = 1 CoordBase::boundary_size_y_lower = 1 CoordBase::boundary_size_z_lower = 1 -CoordBase::boundary_shiftout_x_lower = 1 -CoordBase::boundary_shiftout_y_lower = 1 -CoordBase::boundary_shiftout_z_lower = 1 CoordBase::boundary_size_x_upper = 1 CoordBase::boundary_size_y_upper = 1 CoordBase::boundary_size_z_upper = 1 -CoordBase::boundary_shiftout_x_upper = 0 -CoordBase::boundary_shiftout_y_upper = 0 -CoordBase::boundary_shiftout_z_upper = 0 -CartGrid3D::type = "coordbase" -CartGrid3D::domain = "full" -CartGrid3D::avoid_origin = "no" +# Size of the grid (including boundary points) +PUGH::global_nx = 12 +PUGH::global_ny = 3 +PUGH::global_nz = 3 -Periodic::periodic = "yes" +PUGH::ghost_size = 1 -############################################################# -# Carpet -############################################################# +CartGrid3D::type = "byrange" +CartGrid3D::avoid_origin = "no" -Carpet::ghost_size = 1 -Carpet::domain_from_coordbase = "yes" -Carpet::max_refinement_levels = 1 -#Carpet::init_each_timelevel = "yes" -Carpet::num_integrator_substeps = 4 +CartGrid3D::xmin = -0.1 +CartGrid3D::ymin = -0.1 +CartGrid3D::zmin = -0.1 +CartGrid3D::xmax = 1 +CartGrid3D::ymax = 0.1 +CartGrid3D::zmax = 0.1 + +Periodic::periodic = "yes" ############################################################# # Time integration @@ -95,16 +74,15 @@ SimpleWave::evolved_group_bound = "none" ############################################################# IO::out_dir = $parfile -IO::out_fileinfo = "all" +IO::out_fileinfo = "none" +IO::new_filename_scheme = "no" -CarpetIOBasic::outInfo_every = 1 -CarpetIOBasic::outInfo_vars = "SimpleWave::phi" +IOBasic::outInfo_every = 1 +IOBasic::outInfo_vars = "SimpleWave::phi" IOASCII::out1D_every = 1 IOASCII::out1D_x = "yes" -IOASCII::out1D_y = "yes" -IOASCII::out1D_z = "yes" +IOASCII::out1D_y = "no" +IOASCII::out1D_z = "no" +IOASCII::out1D_d = "no" IOASCII::out1D_vars = "SimpleWave::phi SimpleWave::pi" - -CarpetIOASCII::out_precision = 19 -CarpetIOASCII::out3D_ghosts = "yes" diff --git a/Tools/CodeGen/CactusBoundary.m b/Tools/CodeGen/CactusBoundary.m index 090151c..610987b 100644 --- a/Tools/CodeGen/CactusBoundary.m +++ b/Tools/CodeGen/CactusBoundary.m @@ -69,7 +69,7 @@ GetScheduledGroups[thornName_] := {Name -> "ApplyBCs", Language -> "None", (* groups do not have a language *) SchedulePoint -> "as " <> thornName <> "_ApplyBCs in MoL_PostStep " - <> " after " <> boundariesName[thornName], + <> "after " <> boundariesName[thornName], Comment -> "Apply boundary conditions " <> "controlled by thorn Boundary" } @@ -95,10 +95,10 @@ GetScheduledFunctions[thornName_, evolvedGroups_] := } }; -createBoundTypeParam[groupOrGF_] := { +createBoundTypeParam[groupOrGF_, def_] := { Name -> ToString@groupOrGF <> "_bound", Type -> "KEYWORD", - Default -> "skip", + Default -> def, Description -> "Boundary condition to implement", Visibility -> "private", AllowedValues -> { @@ -108,8 +108,9 @@ createBoundTypeParam[groupOrGF_] := { {Value -> "radiative", Description -> "Radiation boundary condition"}, {Value -> "scalar", Description -> "Dirichlet boundary condition"}, {Value -> "newrad", Description -> "Improved radiative boundary condition"}, - {Value -> "skip", Description -> "skip boundary condition code"} -}}; + {Value -> "skip", Description -> "skip boundary condition code"}}, + Steerable -> Always +}; createBoundSpeedParam[groupOrGF_] := { @@ -119,7 +120,8 @@ createBoundSpeedParam[groupOrGF_] := { Description -> "characteristic speed at boundary", Visibility -> "private", AllowedValues -> {{Value -> "0:*" , - Description -> "outgoing characteristic speed > 0"}} + Description -> "outgoing characteristic speed > 0"}}, + Steerable -> Always }; createBoundLimitParam[groupOrGF_] := { @@ -129,7 +131,8 @@ createBoundLimitParam[groupOrGF_] := { Description -> "limit value for r -> infinity", Visibility -> "private", AllowedValues -> {{Value -> "*:*" , - Description -> "value of limit value is unrestricted"}} + Description -> "value of limit value is unrestricted"}}, + Steerable -> Always }; createBoundScalarParam[groupOrGF_] := { @@ -139,12 +142,13 @@ createBoundScalarParam[groupOrGF_] := { Description -> "Dirichlet boundary value", Visibility -> "private", AllowedValues -> {{Value -> "*:*" , - Description -> "unrestricted"}} + Description -> "unrestricted"}}, + Steerable -> Always }; GetParameters[evolvedGFs_, evolvedGroups_] := - Join[Map[createBoundTypeParam, evolvedGFs], - Map[createBoundTypeParam, Map[unqualifiedGroupName,evolvedGroups]], + Join[Map[createBoundTypeParam[#,"skip"] &, evolvedGFs], + Map[createBoundTypeParam[#,"none"] &, Map[unqualifiedGroupName,evolvedGroups]], Map[createBoundSpeedParam, evolvedGFs], Map[createBoundSpeedParam, Map[unqualifiedGroupName,evolvedGroups]], @@ -169,7 +173,7 @@ GetSources[evolvedGroups_, groups_, implementation_, thornName_] := ExcisionGFs -> evolvedGFs }; - Return[{{Filename -> "Boundaries.c", + Return[{{Filename -> "Boundaries.cc", Contents -> CreateMoLBoundariesSource[boundarySpec]}}]]; End[]; diff --git a/Tools/CodeGen/CalculationFunction.m b/Tools/CodeGen/CalculationFunction.m index 58e7a5f..4f236ef 100644 --- a/Tools/CodeGen/CalculationFunction.m +++ b/Tools/CodeGen/CalculationFunction.m @@ -20,7 +20,7 @@ BeginPackage["CalculationFunction`", {"CodeGen`", "MapLookup`", "KrancGroups`", "Differencing`", "Errors`", - "Helpers`", "Kranc`"}]; + "Helpers`", "Kranc`", "Optimize`", "Jacobian`"}]; CreateCalculationFunction::usage = ""; VerifyCalculation::usage = ""; @@ -159,6 +159,26 @@ removeUnusedShorthands[calc_] := removeUnusedShorthands[newCalc], newCalc]]; +(* Return all the groups that are used in a given calculation *) +groupsInCalculation[calc_, imp_] := + Module[{groups,gfs,eqs,gfsUsed, groupNames}, + groups = lookup[calc, Groups]; + gfs = allGroupVariables[groups]; + eqs = lookup[calc, Equations]; + gfsUsed = Union[Cases[eqs, _ ? (MemberQ[gfs,#] &), Infinity]]; + groupNames = containingGroups[gfsUsed, groups]; + Map[qualifyGroupName[#, imp] &, groupNames]]; + +CheckGroupStorage[groupNames_, calcName_] := + Module[{ignoreGroups, groupsNames2}, + ignoreGroups = {"TmunuBase::stress_energy_scalar", "TmunuBase::stress_energy_vector", + "TmunuBase::stress_energy_tensor"}; + groupNames2 = Select[groupNames, !MemberQ[ignoreGroups, #] &]; + {"\nconst char *groups[] = {", + Riffle[Map[Quote,groupNames2], ","], + "};\n", + "GenericFD_AssertGroupStorage(cctkGH, ", Quote[calcName],", ", Length[groupNames2], ", groups);\n"}]; + (* -------------------------------------------------------------------------- Variables -------------------------------------------------------------------------- *) @@ -180,7 +200,7 @@ localName[x_] := definePreDefinitions[pDefs_] := CommentedBlock["Initialize predefined quantities", - Map[DeclareAssignVariable["CCTK_REAL", #[[1]], #[[2]]] &, pDefs]]; + Map[DeclareAssignVariable[DataType[], #[[1]], #[[2]]] &, pDefs]]; (* -------------------------------------------------------------------------- Equations @@ -203,7 +223,7 @@ printEq[eq_] := split[ x_ + y___] := { x, " + ..."}; split[-x_ + y___] := {-x, " + ..."}; split[ x_ ] := { x, ""}; - rhsSplit = split[Expand@ReplacePowers@rhs]; + rhsSplit = split[Expand@ReplacePowers[rhs,False]]; rhsString = ToString@CForm[rhsSplit[[1]]] <> rhsSplit[[2]]; InfoMessage[InfoFull, " " <> ToString@lhs <> " -> " <> rhsString]]; @@ -228,11 +248,11 @@ simpCollect[collectList_, eqrhs_, localvar_, debug_] := all]; (* Return a CodeGen block which assigns dest by evaluating expr *) -assignVariableFromExpression[dest_, expr_, declare_] := +assignVariableFromExpression[dest_, expr_, declare_, vectorise_] := Module[{tSym, type, cleanExpr, code}, tSym = Unique[]; - type = If[StringMatchQ[ToString[dest], "dir*"], "int", "CCTK_REAL"]; - cleanExpr = ReplacePowers[expr] /. Kranc`t -> tSym; + type = If[StringMatchQ[ToString[dest], "dir*"], "ptrdiff_t", DataType[]]; + cleanExpr = ReplacePowers[expr, vectorise] /. Kranc`t -> tSym; If[SOURCELANGUAGE == "C", code = If[declare, type <> " ", ""] <> ToString[dest] <> " = " <> @@ -258,6 +278,34 @@ assignVariableFromExpression[dest_, expr_, declare_] := {code}]; +(* This and assignVariableFromExpression should be combined *) +generateCodeFromExpression[expr_, vectorise_] := + Module[{tSym, type, cleanExpr, code}, + tSym = Unique[]; + cleanExpr = ReplacePowers[expr, vectorise] /. Kranc`t -> tSym; + + If[SOURCELANGUAGE == "C", + code = + ToString[cleanExpr, CForm, PageWidth -> Infinity], + code = ToString[cleanExpr, FortranForm, PageWidth -> 120]]; + + If[SOURCELANGUAGE != "C", + code = StringReplace[code, "\n " -> " &\n"]; + code = StringReplace[code, " - " -> " & "]; + code = StringReplace[code, ".eq." -> " = "]; + code = StringReplace[code, "= " -> "="]; + code = StringReplace[code, "\\" -> ""]; + code = StringReplace[code, "(index)" -> "(i,j,k)"]]; + + code = StringReplace[code, "normal1" -> "normal[0]"]; + code = StringReplace[code, "normal2" -> "normal[1]"]; + code = StringReplace[code, "normal3" -> "normal[2]"]; + code = StringReplace[code, "BesselJ"-> "gsl_sf_bessel_Jn"]; + code = StringReplace[code, ToString@tSym -> "cctk_time"]; + code = StringReplace[code, "\"" -> ""]; + + {code}]; + (* -------------------------------------------------------------------------- Shorthands -------------------------------------------------------------------------- *) @@ -326,24 +374,32 @@ pdCanonicalOrdering[name_[inds___] -> x_] := Options[CreateCalculationFunction] = ThornOptions; -CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] := +CreateCalculationFunction[calcp_, debug_, imp_, opts:OptionsPattern[]] := Module[{gfs, allSymbols, knownSymbols, - shorts, eqs, parameters, + shorts, eqs, parameters, parameterRules, functionName, dsUsed, groups, pddefs, cleancalc, eqLoop, where, - addToStencilWidth, pDefs, haveCondTextuals, condTextuals}, + addToStencilWidth, pDefs, haveCondTextuals, condTextuals, calc}, + + calc = If[OptionValue[UseJacobian], InsertJacobian[calcp, opts], calcp]; cleancalc = removeUnusedShorthands[calc]; + If[OptionValue[CSE], + cleancalc = EliminateCommonSubexpressions[cleancalc, opts]]; + shorts = lookupDefault[cleancalc, Shorthands, {}]; eqs = lookup[cleancalc, Equations]; parameters = lookupDefault[cleancalc, Parameters, {}]; groups = lookup[cleancalc, Groups]; + If[OptionValue[UseJacobian], groups = Join[groups, JacobianGroups[]]]; pddefs = lookupDefault[cleancalc, PartialDerivatives, {}]; where = lookupDefault[cleancalc, Where, Everywhere]; addToStencilWidth = lookupDefault[cleancalc, AddToStencilWidth, 0]; pDefs = lookup[cleancalc, PreDefinitions]; haveCondTextuals = mapContains[cleancalc, ConditionalOnTextuals]; + SetDataType[If[OptionValue[UseVectors],"CCTK_REAL_VEC", "CCTK_REAL"]]; + VerifyCalculation[cleancalc]; gfs = allGroupVariables[groups]; @@ -368,10 +424,14 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] := If[!lookupDefault[cleancalc, NoSimplify, False], InfoMessage[InfoFull, "Simplifying equations", eqs]; - eqs = Simplify[eqs, {r>0}]]]; + eqs = Simplify[eqs, {r>=0}]]]; InfoMessage[InfoFull, "Equations:"]; + (* Wrap parameters with ToReal *) + parameterRules = Map[(#->ToReal[#])&, parameters]; + eqs = eqs /. parameterRules; + Map[printEq, eqs]; (* Check all the function names *) @@ -393,7 +453,8 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] := allSymbols = calculationSymbols[cleancalc]; knownSymbols = Join[lookupDefault[cleancalc, AllowedSymbols, {}], gfs, shorts, parameters, {dx,dy,dz,idx,idy,idz,t, Pi, E, Symbol["i"], Symbol["j"], Symbol["k"], normal1, normal2, - normal3, tangentA1, tangentA2, tangentA3, tangentB1, tangentB2, tangentB3}]; + normal3, tangentA1, tangentA2, tangentA3, tangentB1, tangentB2, tangentB3}, + If[OptionValue[UseJacobian], JacobianSymbols[], {}]]; unknownSymbols = Complement[allSymbols, knownSymbols]; @@ -402,7 +463,7 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] := "Calculation is:", cleancalc]]; { - DefineFunction[bodyFunctionName, "void", "cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]", + DefineFunction[bodyFunctionName, "static void", "cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]", { "DECLARE_CCTK_ARGUMENTS;\n", "DECLARE_CCTK_PARAMETERS;\n\n", @@ -415,14 +476,21 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] := ConditionalOnParameterTextual["cctk_iteration % " <> functionName <> "_calc_every != " <> functionName <> "_calc_offset", "return;\n"], + CheckGroupStorage[groupsInCalculation[cleancalc, imp], functionName], + + "\n", + CheckStencil[pddefs, eqs, functionName], + If[haveCondTextuals, Map[ConditionalOnParameterTextual["!(" <> # <> ")", "return;\n"] &,condTextuals], {}], CommentedBlock["Include user-supplied include files", Map[IncludeFile, lookupDefault[cleancalc, DeclarationIncludes, {}]]], - InitialiseFDVariables[], + InitialiseFDVariables[OptionValue[UseVectors]], definePreDefinitions[pDefs], + If[OptionValue[UseJacobian], CreateJacobianVariables[], {}], + If[Cases[{pddefs}, SBPDerivative[_], Infinity] != {}, CommentedBlock["Compute Summation By Parts derivatives", IncludeFile["sbp_calc_coeffs.h"]], @@ -431,7 +499,7 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] := If[gfs != {}, { eqLoop = equationLoop[eqs, cleancalc, gfs, shorts, {}, groups, - pddefs, where, addToStencilWidth, useCSE, opts]}, + pddefs, where, addToStencilWidth, opts]}, ConditionalOnParameterTextual["verbose > 1", "CCTK_VInfo(CCTK_THORNSTRING,\"Leaving " <> bodyFunctionName <> "\");\n"], @@ -464,12 +532,13 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] := Options[equationLoop] = ThornOptions; equationLoop[eqs_, cleancalc_, gfs_, shorts_, incs_, groups_, pddefs_, - where_, addToStencilWidth_, useCSE_, + where_, addToStencilWidth_, opts:OptionsPattern[]] := Module[{rhss, lhss, gfsInRHS, gfsInLHS, gfsOnlyInRHS, localGFs, localMap, eqs2, derivSwitch, code, functionName, calcCode, loopFunction, gfsInBoth, gfsDifferentiated, - gfsDifferentiatedAndOnLHS, declare, eqsReplaced}, + gfsDifferentiatedAndOnLHS, declare, eqsReplaced, + generateEquationsCode}, rhss = Map[#[[2]] &, eqs]; lhss = Map[#[[1]] &, eqs]; @@ -520,9 +589,6 @@ equationLoop[eqs_, cleancalc_, gfs_, shorts_, incs_, groups_, pddefs_, (* Replace grid functions with their local forms *) eqsReplaced = eqs2 /. localMap; - If[useCSE, - eqsReplaced = CSE[eqsReplaced]]; - (* Construct a list, corresponding to the list of equations, marking those which need their LHS variables declared. We declare variables at the same time as assigning to them as it @@ -531,22 +597,80 @@ equationLoop[eqs_, cleancalc_, gfs_, shorts_, incs_, groups_, pddefs_, functions which appear in the RHSs have been declared and set already (DeclareMaybeAssignVariableInLoop below), so assignments to these do not generate declarations here. *) - declare = markFirst[First /@ eqsReplaced, Map[localName, gfsInRHS]]; + declare = Block[{$RecursionLimit=Infinity},markFirst[First /@ eqsReplaced, Map[localName, gfsInRHS]]]; + + (* Generate the code for the equations, factoring out any + sequential IfThen statements with the same condition. Try to + declare variables as they are assigned, but it is only possible + to do this outside all if(){} statements. "declare2" is a list + of booleans indicating if the LHS variables should be declared + as they are assigned. *) + generateEquationsCode[eqs2_, declare2_] := + Module[{ifs, ind, cond, num, preDeclare}, + ifs = Position[eqs2, _ -> IfThen[__], {1}]; + If[Length[ifs] <= 1, + Riffle[MapThread[assignVariableFromExpression[#1[[1]], #1[[2]], #2, OptionValue[UseVectors]] &, + {eqs2, declare2}],"\n"], + (* else *) + ind = ifs[[1,1]]; + If[ind > 1, + {generateEquationsCode[Take[eqs2, ind-1], Take[declare2, ind-1]], "\n", + generateEquationsCode[Drop[eqs2, ind-1], Drop[declare2, ind-1]]}, + (* else *) + cond = eqs2[[1,2,1]]; + num = LengthWhile[eqs2, MatchQ[#, _ -> IfThen[cond, _, _]] &]; + If[num == 1, + {generateEquationsCode[Take[eqs2,1], Take[declare2,1]], "\n", + generateEquationsCode[Drop[eqs2,1], Drop[declare2,1]]}, + (* else *) + e1 = Take[eqs2, num]; + e2 = Drop[eqs2, num]; + preDeclare = #[[2,1]] & /@ Select[MapThread[List, {Take[declare2,num], e1}], #[[1]] &]; + {Map[DeclareVariableNoInit[#, DataType[]] &, Complement[Union[preDeclare], localName/@gfsInRHS]], {"\n"}, + {Conditional[generateCodeFromExpression[cond, OptionValue[UseVectors]], + generateEquationsCode[Map[#[[1]] -> #[[2,2]]&, e1], ConstantArray[False, Length[e1]]], + generateEquationsCode[Map[#[[1]] -> #[[2,3]]&, e1], ConstantArray[False, Length[e1]]]], + "\n", + generateEquationsCode[e2,Drop[declare2,Length[e1]]]}}]]]]; + + calcCode = generateEquationsCode[eqsReplaced, declare]; + + + assignLocalGridFunctions[gs_, useVectors_, useJacobian_] := + Module[{conds, varPatterns, varsInConds, simpleVars, code}, + conds = + {{"eT" ~~ _ ~~ _, "*stress_energy_state", "ToReal(0.0)"}}; (* This should be passed as an option *) + If[useJacobian, + conds = Append[conds, JacobianConditionalGridFunctions[]]]; + + varPatterns = Map[First, conds]; + varsInConds = Map[Function[pattern, Select[gs,StringMatchQ[ToString[#], pattern] &]], varPatterns]; + simpleVars = Complement[gs, Flatten[varsInConds]]; + code = {"\n", + Map[DeclareMaybeAssignVariableInLoop[ + DataType[], localName[#], GridName[#], + False,"", useVectors] &, simpleVars], + {"\n", + Riffle[ + MapThread[ + If[Length[#2] > 0, + {DeclareVariables[localName/@#2, DataType[]],"\n", + Conditional[#1, + Table[AssignVariableInLoop[localName[var], GridName[var], useVectors], {var, #2}], + Sequence@@If[#3 =!= None, {Table[AssignVariableInLoop[localName[var], #3, False (*useVectors*)], {var, #2}]}, {}]]}, + (* else *) + {}] &, + {Map[#[[2]]&, conds], varsInConds, Map[#[[3]]&, conds]}], "\n"]}}; + code + ]; - calcCode = - MapThread[{assignVariableFromExpression[#1[[1]], #1[[2]], #2], "\n"} &, - {eqsReplaced, declare}]; GenericGridLoop[functionName, { - DeclareDerivatives[defsWithoutShorts, eqsOrdered], + (* DeclareDerivatives[defsWithoutShorts, eqsOrdered], *) CommentedBlock["Assign local copies of grid functions", - Map[DeclareMaybeAssignVariableInLoop[ - "CCTK_REAL", localName[#], GridName[#], - StringMatchQ[ToString[GridName[#]], "eT" ~~ _ ~~ _ ~~ "[" ~~ __ ~~ "]"], - "*stress_energy_state"] &, - gfsInRHS]], + assignLocalGridFunctions[gfsInRHS, OptionValue[UseVectors], OptionValue[UseJacobian]]], CommentedBlock["Include user supplied include files", Map[IncludeFile, incs]], @@ -560,8 +684,34 @@ equationLoop[eqs_, cleancalc_, gfs_, shorts_, incs_, groups_, pddefs_, Map[InfoVariable[#[[1]]] &, (eqs2 /. localMap)], ""], + If[OptionValue[UseVectors], { + CommentedBlock["If necessary, store only partial vectors after the first iteration", + ConditionalOnParameterTextual["CCTK_REAL_VEC_SIZE > 2 && CCTK_BUILTIN_EXPECT(i < lc_imin && i+CCTK_REAL_VEC_SIZE > lc_imax, 0)", + { + DeclareAssignVariable["ptrdiff_t", "elt_count_lo", "lc_imin-i"], + DeclareAssignVariable["ptrdiff_t", "elt_count_hi", "lc_imax-i"], + Map[StoreMiddlePartialVariableInLoop[GridName[#], localName[#], "elt_count_lo", "elt_count_hi"] &, + gfsInLHS], + "break;\n" + }]], + CommentedBlock["If necessary, store only partial vectors after the first iteration", + ConditionalOnParameterTextual["CCTK_REAL_VEC_SIZE > 1 && CCTK_BUILTIN_EXPECT(i < lc_imin, 0)", + { + DeclareAssignVariable["ptrdiff_t", "elt_count", "lc_imin-i"], + Map[StoreHighPartialVariableInLoop[GridName[#], localName[#], "elt_count"] &, + gfsInLHS], + "continue;\n" + }]], + CommentedBlock["If necessary, store only partial vectors after the last iteration", + ConditionalOnParameterTextual["CCTK_REAL_VEC_SIZE > 1 && CCTK_BUILTIN_EXPECT(i+CCTK_REAL_VEC_SIZE > lc_imax, 0)", + { + DeclareAssignVariable["ptrdiff_t", "elt_count", "lc_imax-i"], + Map[StoreLowPartialVariableInLoop[GridName[#], localName[#], "elt_count"] &, + gfsInLHS], + "break;\n" + }]]}, {}], CommentedBlock["Copy local copies back to grid functions", - Map[AssignVariableInLoop[GridName[#], localName[#]] &, + Map[(If[OptionValue[UseVectors], StoreVariableInLoop, AssignVariableInLoop][GridName[#], localName[#]]) &, gfsInLHS]], If[debugInLoop, Map[InfoVariable[GridName[#]] &, gfsInLHS], ""]}, opts]]; diff --git a/Tools/CodeGen/CodeGen.m b/Tools/CodeGen/CodeGen.m index 7cae5fb..5e8cd6a 100644 --- a/Tools/CodeGen/CodeGen.m +++ b/Tools/CodeGen/CodeGen.m @@ -43,6 +43,9 @@ IncludeSystemFile::usage = "IncludeFile[name] returns a block of code" <> DeclareVariable::usage = "DeclareVariable[name, type] returns a block of code " <> "that declares a variable of given name and type. 'name' and 'type' should be " <> "strings."; +DeclareVariableNoInit::usage = "DeclareVariableNoInit[name, type] returns a block of code " <> + "that declares a variable of given name and type without initialising it. 'name' and 'type' should be " <> + "strings."; DeclareVariables::usage = "DeclareVariables[names, type] returns a block of code " <> "that declares a list of variables of given name and type. 'names' should be a list" <> " of strings and 'type' should be a string string."; @@ -60,14 +63,20 @@ DeclareAssignVariable::usage = "DeclareAssignVariable[type_, dest_, src_] return "that declares and sets a constant variable of given name and type."; AssignVariableInLoop::usage = "AssignVariableInLoop[dest_, src_] returns a block of code " <> "that assigns 'src' to 'dest'."; +StoreVariableInLoop::usage = "StoreVariableInLoop[dest_, src_] returns a block of code " <> + "that assigns 'src' to 'dest'."; +StoreLowPartialVariableInLoop::usage = "StoreLowPartialVariableInLoop[dest_, src_, count_] returns a block of code " <> + "that assigns 'src' to 'dest'."; +StoreHighPartialVariableInLoop::usage = "StoreHighPartialVariableInLoop[dest_, src_, count_] returns a block of code " <> + "that assigns 'src' to 'dest'."; +StoreMiddlePartialVariableInLoop::usage = "StoreMiddlePartialVariableInLoop[dest_, src_, countLow_, countHigh_] returns a block of code " <> + "that assigns 'src' to 'dest'."; DeclareAssignVariableInLoop::usage = "DeclareAssignVariableInLoop[type_, dest_, src_] returns a block of code " <> "that assigns 'src' to 'dest'."; MaybeAssignVariableInLoop::usage = "MaybeAssignVariableInLoop[dest_, src_, cond_] returns a block of code " <> "that assigns 'src' to 'dest'."; DeclareMaybeAssignVariableInLoop::usage = "DeclareMaybeAssignVariableInLoop[type_, dest_, src_, cond_] returns a block of code " <> "that assigns 'src' to 'dest'."; -DeclareVariablesInLoopVectorised::usage = ""; -AssignVariablesInLoopVectorised::usage = ""; TestForNaN::usage = "TestForNaN[expr_] returns a block of code " <> "that tests 'expr' for nan."; CommentedBlock::usage = "CommentedBlock[comment, block] returns a block consisting " <> @@ -117,7 +126,6 @@ CommaNewlineSeparated::usage = ""; (* This should not really be in CodeGen *) CommaSeparated::usage = ""; ReplacePowers::usage = ""; CFormHideStrings::usage = ""; -CSE::usage = ""; BoundaryLoop::usage = ""; BoundaryWithGhostsLoop::usage = ""; GenericGridLoop::usage = ""; @@ -125,15 +133,18 @@ GenericGridLoop::usage = ""; NameRoot::usage = ""; PartitionVarList::usage = ""; Quote::usage = "Quote[x] returns x surrounded by quotes"; +DataType::usage = "DataType[] returns a string for the grid function data type (e.g. CCTK_REAL)"; +SetDataType::usage = "SetDataType[type] sets a string for the grid function data type (e.g. CCTK_REAL)"; +Conditional; Begin["`Private`"]; SOURCELANGUAGE = "C"; -SOURCESUFFIX = ".c"; +SOURCESUFFIX = ".cc"; setSourceSuffix[lang_] := If[ (lang == "C"), - SOURCESUFFIX = ".c"; + SOURCESUFFIX = ".cc"; , SOURCESUFFIX = ".F90"; ]; @@ -146,10 +157,18 @@ If[ (lang == "C" || lang == "Fortran"), InfoMessage[Terse, "User set source language to " <> lang], SOURCELANGUAGE = "C"; - setSourceSuffix[".c"]; + setSourceSuffix[".cc"]; InfoMessage[Terse, "Setting Source Language to C"]; ]; +SetDataType[type_String] := + dataType = type; + +DataType[] := + If[dataType === Symbol["datatype"], + Throw["DataType: Have not set a data type"], + dataType]; + (* Code generation utilities; not specific to any language *) FlattenBlock[b_] := Apply[StringJoin,Map[ToString,If[! AtomQ[b], Flatten[b, Infinity], b]]]; @@ -234,6 +253,12 @@ If[SOURCELANGUAGE == "C", {type, " :: ", name, EOL[]} (* no value init here to avoid implicit SAVE attribute *) ]; +DeclareVariableNoInit[name_, type_] := +If[SOURCELANGUAGE == "C", + {type, " ", name, EOL[]}, + {type, " :: ", name, EOL[]} (* no value init here to avoid implicit SAVE attribute *) + ]; + DeclareVariables[names_?ListQ, type_] := If[SOURCELANGUAGE == "C", @@ -273,69 +298,41 @@ AssignVariable[dest_, src_] := DeclareAssignVariable[type_, dest_, src_] := {type, " const ", dest, " = ", src, EOL[]}; -AssignVariableInLoop[dest_, src_] := - {dest, " = ", src, EOL[]}; +AssignVariableInLoop[dest_, src_, vectorise_:False] := + Module[{loader}, + loader[x_] := If[vectorise, {"vec_load(", x, ")"}, x]; + {dest, " = ", loader[src], EOL[]}]; (* {dest, " = ", src, EOL[], TestForNaN[dest]}; *) -DeclareAssignVariableInLoop[type_, dest_, src_] := - {type, " const ", dest, " = ", src, EOL[]}; +StoreVariableInLoop[dest_, src_] := + {"vec_store_nta(", dest, ",", src, ")", EOL[]}; -MaybeAssignVariableInLoop[dest_, src_, cond_] := - If [cond, - {dest, " = useMatter ? ", src, " : 0.0", EOL[]}, - {dest, " = ", src, EOL[]}]; +StoreLowPartialVariableInLoop[dest_, src_, count_] := + {"vec_store_nta_partial_lo(", dest, ",", src, ",", count, ")", EOL[]}; -DeclareMaybeAssignVariableInLoop[type_, dest_, src_, mmaCond_, codeCond_] := - If [mmaCond, - {type, " ", dest, " = (", codeCond, ") ? (", src, ") : 0.0", EOL[]}, - {type, " ", dest, " = ", src, EOL[]}]; +StoreHighPartialVariableInLoop[dest_, src_, count_] := + {"vec_store_nta_partial_hi(", dest, ",", src, ",", count, ")", EOL[]}; -(* TODO: move these into OpenMP loop *) -DeclareVariablesInLoopVectorised[dests_, temps_, srcs_] := - { - {"#undef LC_PRELOOP_STATEMENTS", "\n"}, - {"#define LC_PRELOOP_STATEMENTS", " \\\n"}, - {"int const GFD_imin = lc_imin + ((lc_imin + cctk_lsh[0] * (j + cctk_lsh[1] * k)) & (CCTK_REAL_VEC_SIZE-1))", "; \\\n"}, - {"int const GFD_imax = lc_imax + ((lc_imax + cctk_lsh[0] * (j + cctk_lsh[1] * k)) & (CCTK_REAL_VEC_SIZE-1)) - CCTK_REAL_VEC_SIZE", "; \\\n"}, - Map[Function[x, Module[{dest, temp, src}, - {dest, temp, src} = x; - {"CCTK_REAL_VEC ", temp, "; \\\n"}]], - Transpose[{dests, temps, srcs}]], - {"\n"} - }; +StoreMiddlePartialVariableInLoop[dest_, src_, countLow_, countHigh_] := + {"vec_store_nta_partial_mid(", dest, ",", src, ",", countLow, ",", countHigh, ")", EOL[]}; -AssignVariablesInLoopVectorised[dests_, temps_, srcs_] := - { - {"{\n"}, - {" if (i < GFD_imin || i >= GFD_imax) {\n"}, - Map[Function[x, Module[{dest, temp, src}, - {dest, temp, src} = x; - {" ", dest, "[index] = ", src, EOL[]}]], - Transpose[{dests, temps, srcs}]], - {" } else {\n"}, - {" size_t const index0 = index & (CCTK_REAL_VEC_SIZE-1)", EOL[]}, - Map[Function[x, Module[{dest, temp, src}, - {dest, temp, src} = x; - {" ((CCTK_REAL*)&", temp, ")[index0] = ", - src, EOL[]}]], - Transpose[{dests, temps, srcs}]], - {" if (index0 == CCTK_REAL_VEC_SIZE-1) {\n"}, - {" size_t const index1 = index - (CCTK_REAL_VEC_SIZE-1)", EOL[]}, - Map[Function[x, Module[{dest, temp, src}, - {dest, temp, src} = x; - {" _mm_stream_pd (&", dest, "[index1], ", - temp, ")", EOL[]}]], - Transpose[{dests, temps, srcs}]], - {" }\n"}, - {" }\n"}, - {"}\n"} - }; +DeclareAssignVariableInLoop[type_, dest_, src_] := + {type, " const ", dest, " = vec_load(", src, ")", EOL[]}; + +MaybeAssignVariableInLoop[dest_, src_, cond_] := + If [cond, + {dest, " = useMatter ? vec_load(", src, ") : ToReal(0.0)", EOL[]}, + {dest, " = vec_load(", src, ")", EOL[]}]; -AssignVariableInLoopsVectorised[dest_, temp_, src_] := - {"GFD_save_and_store(", dest, ",", "index", ",", "&", temp, ",", src, ")", EOL[]}; +DeclareMaybeAssignVariableInLoop[type_, dest_, src_, mmaCond_, codeCond_, vectorise_:False] := + Module[{loader}, + loader[x_] := If[vectorise, {"vec_load(", x, ")"}, x]; + If [mmaCond, + {type, " ", dest, " = (", codeCond, ") ? ", loader[src], " : ToReal(0.0)", EOL[]}, + {type, " ", dest, " = ", loader[src], EOL[]}]]; TestForNaN[expr_] := {"if (isnan(", expr, ")) {\n", @@ -343,7 +340,7 @@ TestForNaN[expr_] := " CCTK_VInfo(CCTK_THORNSTRING, \"ipos: %d %d %d\", i, j, k);\n", " CCTK_VInfo(CCTK_THORNSTRING, \"lbnd: %d %d %d\", cctk_lbnd[0], cctk_lbnd[1], cctk_lbnd[2]);\n", " CCTK_VInfo(CCTK_THORNSTRING, \"lsh: %d %d %d\", cctk_lsh[0], cctk_lsh[1], cctk_lsh[2]);\n", - " CCTK_VInfo(CCTK_THORNSTRING, \"lssh: %d %d %d\", cctk_lssh[CCTK_LSSH_IDX(0,0)], cctk_lssh[CCTK_LSSH_IDX(0,1)], cctk_lssh[CCTK_LSSH_IDX(0,2)]);\n", + " CCTK_VInfo(CCTK_THORNSTRING, \"LSSH: %d %d %d\", CCTK_LSSH(0,0), CCTK_LSSH(0,1), CCTK_LSSH(0,2));\n", " CCTK_VInfo(CCTK_THORNSTRING, \"", expr, ": %.17g\", (double)", expr, ");\n", "}\n"}; @@ -406,7 +403,7 @@ defineSubroutine[name_, args_, contents_] := defineSubroutineC[name_, args_, contents_] := SeparatedBlock[ - {"void ", name, "(", args, ")", "\n", + {"extern \"C\" void ", name, "(", args, ")", "\n", CBlock[contents]}]; defineSubroutineF[name_, args_, contents_] := @@ -426,7 +423,7 @@ defineSubroutineF[name_, args_, contents_] := (* This is a Cactus-callable function *) DefineCCTKFunction[name_, type_, contents_] := - DefineFunction[name, type, "CCTK_ARGUMENTS", + DefineFunction[name, "extern \"C\" " <> type, "CCTK_ARGUMENTS", { "DECLARE_CCTK_ARGUMENTS;\n", "DECLARE_CCTK_PARAMETERS;\n\n", @@ -451,47 +448,65 @@ DeclareFDVariables[] := {khalf,kthird,ktwothird,kfourthird,keightthird}, {"hdxi", "hdyi", "hdzi"}}], "\n"}, - {Map[DeclareVariables[#, "int"] &, {{"di", "dj", "dk"}}], + {Map[DeclareVariables[#, "ptrdiff_t"] &, {{"di", "dj", "dk"}}], "\n"}]; *) CommentedBlock["Declare finite differencing variables", {}]; -InitialiseFDSpacingVariablesC[] := +InitialiseFDSpacingVariablesC[vectorise_:False] := { - DeclareAssignVariable["CCTK_REAL", "dx", "CCTK_DELTA_SPACE(0)"], - DeclareAssignVariable["CCTK_REAL", "dy", "CCTK_DELTA_SPACE(1)"], - DeclareAssignVariable["CCTK_REAL", "dz", "CCTK_DELTA_SPACE(2)"], - (* DeclareAssignVariable["int", "di", "CCTK_GFINDEX3D(cctkGH,1,0,0) - CCTK_GFINDEX3D(cctkGH,0,0,0)"], *) - DeclareAssignVariable["int", "di", "1"], - DeclareAssignVariable["int", "dj", "CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0)"], - DeclareAssignVariable["int", "dk", "CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0)"] + (* DeclareAssignVariable["ptrdiff_t", "di", "CCTK_GFINDEX3D(cctkGH,1,0,0) - CCTK_GFINDEX3D(cctkGH,0,0,0)"], *) + DeclareAssignVariable["ptrdiff_t", "di", "1"], + DeclareAssignVariable["ptrdiff_t", "dj", "CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0)"], + DeclareAssignVariable["ptrdiff_t", "dk", "CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0)"], + If[vectorise, + { + DeclareAssignVariable["ptrdiff_t", "cdi", "sizeof(CCTK_REAL) * di"], + DeclareAssignVariable["ptrdiff_t", "cdj", "sizeof(CCTK_REAL) * dj"], + DeclareAssignVariable["ptrdiff_t", "cdk", "sizeof(CCTK_REAL) * dk"] + }, + ""], + DeclareAssignVariable[DataType[], "dx", "ToReal(CCTK_DELTA_SPACE(0))"], + DeclareAssignVariable[DataType[], "dy", "ToReal(CCTK_DELTA_SPACE(1))"], + DeclareAssignVariable[DataType[], "dz", "ToReal(CCTK_DELTA_SPACE(2))"], + DeclareAssignVariable[DataType[], "dt", "ToReal(CCTK_DELTA_TIME)"] }; InitialiseFDSpacingVariablesFortran[] := { + AssignVariable["dt", "CCTK_DELTA_TIME"], AssignVariable["dx", "CCTK_DELTA_SPACE(1)"], AssignVariable["dy", "CCTK_DELTA_SPACE(2)"], AssignVariable["dz", "CCTK_DELTA_SPACE(3)"] } -InitialiseFDVariables[] := +InitialiseFDVariables[vectorise_] := CommentedBlock["Initialise finite differencing variables", { If[SOURCELANGUAGE == "Fortran", InitialiseFDSpacingVariablesFortran[], - InitialiseFDSpacingVariablesC[]], + InitialiseFDSpacingVariablesC[vectorise]], - DeclareAssignVariable["CCTK_REAL", "dxi", "1.0 / dx"], - DeclareAssignVariable["CCTK_REAL", "dyi", "1.0 / dy"], - DeclareAssignVariable["CCTK_REAL", "dzi", "1.0 / dz"], - DeclareAssignVariable["CCTK_REAL", "khalf", "0.5"], - DeclareAssignVariable["CCTK_REAL", "kthird", "1/3.0"], - DeclareAssignVariable["CCTK_REAL", "ktwothird", "2.0/3.0"], - DeclareAssignVariable["CCTK_REAL", "kfourthird", "4.0/3.0"], - DeclareAssignVariable["CCTK_REAL", "keightthird", "8.0/3.0"], - DeclareAssignVariable["CCTK_REAL", "hdxi", "0.5 * dxi"], - DeclareAssignVariable["CCTK_REAL", "hdyi", "0.5 * dyi"], - DeclareAssignVariable["CCTK_REAL", "hdzi", "0.5 * dzi"]}]; + DeclareAssignVariable[DataType[], "dxi", "INV(dx)"], + DeclareAssignVariable[DataType[], "dyi", "INV(dy)"], + DeclareAssignVariable[DataType[], "dzi", "INV(dz)"], + If[vectorise, + {DeclareAssignVariable[DataType[], "khalf", "ToReal(0.5)"], + DeclareAssignVariable[DataType[], "kthird", "ToReal(1.0/3.0)"], + DeclareAssignVariable[DataType[], "ktwothird", "ToReal(2.0/3.0)"], + DeclareAssignVariable[DataType[], "kfourthird", "ToReal(4.0/3.0)"], + DeclareAssignVariable[DataType[], "keightthird", "ToReal(8.0/3.0)"], + DeclareAssignVariable[DataType[], "hdxi", "kmul(ToReal(0.5), dxi)"], + DeclareAssignVariable[DataType[], "hdyi", "kmul(ToReal(0.5), dyi)"], + DeclareAssignVariable[DataType[], "hdzi", "kmul(ToReal(0.5), dzi)"]}, + {DeclareAssignVariable[DataType[], "khalf", "0.5"], + DeclareAssignVariable[DataType[], "kthird", "1/3.0"], + DeclareAssignVariable[DataType[], "ktwothird", "2.0/3.0"], + DeclareAssignVariable[DataType[], "kfourthird", "4.0/3.0"], + DeclareAssignVariable[DataType[], "keightthird", "8.0/3.0"], + DeclareAssignVariable[DataType[], "hdxi", "0.5 * dxi"], + DeclareAssignVariable[DataType[], "hdyi", "0.5 * dyi"], + DeclareAssignVariable[DataType[], "hdzi", "0.5 * dzi"]}]}]; GridName[x_] := If[SOURCELANGUAGE == "C", ToExpression[ToString[x] <> "[index]"], @@ -543,9 +558,9 @@ InitialiseGridLoopVariables[derivativesUsedSwitch_, addToStencilWidth_] := AssignVariable["kstart", arrayIndex["index_offset_z"]], "\n", - AssignVariable["iend", {arrayElement["cctk_lssh", "CCTK_LSSH_IDX(0,0)"], " - index_offset_x"}], - AssignVariable["jend", {arrayElement["cctk_lssh", "CCTK_LSSH_IDX(0,1)"], " - index_offset_y"}], - AssignVariable["kend", {arrayElement["cctk_lssh", "CCTK_LSSH_IDX(0,2)"], " - index_offset_z"}] + AssignVariable["iend", {"CCTK_LSSH(0,0) - index_offset_x"}], + AssignVariable["jend", {"CCTK_LSSH(0,1) - index_offset_y"}], + AssignVariable["kend", {"CCTK_LSSH(0,2) - index_offset_z"}] }, { @@ -554,9 +569,9 @@ InitialiseGridLoopVariables[derivativesUsedSwitch_, addToStencilWidth_] := AssignVariable["kstart", arrayIndex[0]], "\n", - AssignVariable["iend", arrayElement["cctk_lssh", "CCTK_LSSH_IDX(0,0)"]], - AssignVariable["jend", arrayElement["cctk_lssh", "CCTK_LSSH_IDX(0,1)"]], - AssignVariable["kend", arrayElement["cctk_lssh", "CCTK_LSSH_IDX(0,2)"]] + AssignVariable["iend", "CCTK_LSSH(0,0)"], + AssignVariable["jend", "CCTK_LSSH(0,1)"], + AssignVariable["kend", "CCTK_LSSH(0,2)"] }] ]; @@ -629,7 +644,7 @@ Options[GenericGridLoop] = ThornOptions; GenericGridLoop[functionName_, block_, opts:OptionsPattern[]] := If[OptionValue[UseLoopControl], - GenericGridLoopUsingLoopControl[functionName, block], + GenericGridLoopUsingLoopControl[functionName, block, OptionValue[UseVectors]], GenericGridLoopTraditional[block]]; GenericGridLoopTraditional[block_] := @@ -647,24 +662,26 @@ GenericGridLoopTraditional[block_] := } ]]]]; -GenericGridLoopUsingLoopControl[functionName_, block_] := +GenericGridLoopUsingLoopControl[functionName_, block_, vectorise_] := If[SOURCELANGUAGE == "C", CommentedBlock["Loop over the grid points", { "#pragma omp parallel\n", - "LC_LOOP3 (", functionName, ",\n", - " i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],\n", - " cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])\n", + If[vectorise, "LC_LOOP3VEC", "LC_LOOP3"] <> " (", functionName, ",\n", + " i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],\n", + " cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]", If[vectorise, {",\n", + " CCTK_REAL_VEC_SIZE"},""] <> ")\n", "{\n", indentBlock[ { - DeclareVariable["index", "// int"], - DeclareAssignVariable["int", "index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"], + (* DeclareVariable["index", "// int"], *) + (* DeclareAssignVariable["int", "index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"], *) + DeclareAssignVariable["ptrdiff_t", "index", "di*i + dj*j + dk*k"], block } ], "}\n", - "LC_ENDLOOP3 (", functionName, ");\n" + If[vectorise, "LC_ENDLOOP3VEC", "LC_ENDLOOP3"] <> " (", functionName, ");\n" } ], "" @@ -691,9 +708,9 @@ BoundaryLoop[block_] := AssignVariable[arrayElement["bmin", 0], "is_physbnd[0*2+0] ? 0 : imin[0]"], AssignVariable[arrayElement["bmin", 1], "is_physbnd[1*2+0] ? 0 : imin[1]"], AssignVariable[arrayElement["bmin", 2], "is_physbnd[2*2+0] ? 0 : imin[2]"], - AssignVariable[arrayElement["bmax", 0], "is_physbnd[0*2+1] ? cctk_from[CCTK_LSSH_IDX(0,0)] : imax[0]"], - AssignVariable[arrayElement["bmax", 1], "is_physbnd[1*2+1] ? cctk_from[CCTK_LSSH_IDX(0,1)] : imax[1]"], - AssignVariable[arrayElement["bmax", 2], "is_physbnd[2*2+1] ? cctk_from[CCTK_LSSH_IDX(0,2)] : imax[2]"]}], + AssignVariable[arrayElement["bmax", 0], "is_physbnd[0*2+1] ? CCTK_LSSH(0,0) : imax[0]"], + AssignVariable[arrayElement["bmax", 1], "is_physbnd[1*2+1] ? CCTK_LSSH(0,1) : imax[1]"], + AssignVariable[arrayElement["bmax", 2], "is_physbnd[2*2+1] ? CCTK_LSSH(0,2) : imax[2]"]}], CommentedBlock["Loop over all faces", loopOverInteger["dir", "0", "3", @@ -704,7 +721,7 @@ BoundaryLoop[block_] := {0, {AssignVariable[arrayElement["bmax", "dir"], {arrayElement["imin", "dir"], ""}], AssignVariable[arrayElement["bmin", "dir"], {0, ""}]}}, {1, {AssignVariable[arrayElement["bmin", "dir"], {arrayElement["imax", "dir"], "" }], - AssignVariable[arrayElement["bmax", "dir"], {"cctk_lssh[CCTK_LSSH_IDX(0,dir)]", ""}]}}]], + AssignVariable[arrayElement["bmax", "dir"], {"CCTK_LSSH(0,dir)", ""}]}}]], conditional[arrayElement["is_physbnd", "dir * 2 + face"], loopOverInteger["k", arrayElement["bmin",2], arrayElement["bmax",2], loopOverInteger["j", arrayElement["bmin",1], arrayElement["bmax",1], @@ -728,9 +745,9 @@ BoundaryWithGhostsLoop[block_] := AssignVariable[arrayElement["bmin", 0], "0"], AssignVariable[arrayElement["bmin", 1], "0"], AssignVariable[arrayElement["bmin", 2], "0"], - AssignVariable[arrayElement["bmax", 0], "cctk_lssh[CCTK_LSSH_IDX(0,0)]"], - AssignVariable[arrayElement["bmax", 1], "cctk_lssh[CCTK_LSSH_IDX(0,1)]"], - AssignVariable[arrayElement["bmax", 2], "cctk_lssh[CCTK_LSSH_IDX(0,2)]"]}], + AssignVariable[arrayElement["bmax", 0], "CCTK_LSSH(0,0)"], + AssignVariable[arrayElement["bmax", 1], "CCTK_LSSH(0,1)"], + AssignVariable[arrayElement["bmax", 2], "CCTK_LSSH(0,2)"]}], CommentedBlock["Loop over all faces", loopOverInteger["dir", "0", "3", @@ -741,7 +758,7 @@ BoundaryWithGhostsLoop[block_] := {0, {AssignVariable[arrayElement["bmax", "dir"], {arrayElement["imin", "dir"], ""}], AssignVariable[arrayElement["bmin", "dir"], {0, ""}]}}, {1, {AssignVariable[arrayElement["bmin", "dir"], {arrayElement["imax", "dir"], "" }], - AssignVariable[arrayElement["bmax", "dir"], {"cctk_lssh[CCTK_LSSH_IDX(0,dir)]", ""}]}}]], + AssignVariable[arrayElement["bmax", "dir"], {"CCTK_LSSH(0,dir)]", ""}]}}]], conditional[arrayElement["is_physbnd", "dir * 2 + face"], loopOverInteger["k", arrayElement["bmin",2], arrayElement["bmax",2], loopOverInteger["j", arrayElement["bmin",1], arrayElement["bmax",1], @@ -757,12 +774,28 @@ BoundaryWithGhostsLoop[block_] := ]} ]]]}; -conditional[condition_, block_] := - {"if (", condition, ")\n", +(* Remove call to ToReal from condtion. This should really instead be + done much earlier. Sometimes Conditional is called with a + conditional that is an expression, sometimes it is a list of + strings, so we need to handle both cases. *) +(* One approach to remove calls to ToReal is to wrap the condition + into a function call, such as e.g. Cond[x], which is then handled by + the vectorisation code in the same way the IfThen[x,y,z] expression + is handled there. *) +removeToRealFunc[condition_] := Replace[condition, {ToReal[x_] -> x}] +removeToRealString[condition_] := If[StringQ[condition], StringReplace[condition, "ToReal(" ~~ x__ ~~ ")" :> x], condition] +removeToReal[condition_] := Map[removeToRealString, removeToRealFunc[condition]] + +Conditional[condition_, block_] := + {"if (", removeToReal[condition], ")\n", CBlock[block]}; +Conditional[condition_, block1_, block2_] := + {"if (", removeToReal[condition], ")\n", + CBlock[block1], "else\n", CBlock[block2]}; + onceInGridLoop[block_] := - conditional["i == 5 && j == 5 && k == 5", + Conditional["i == 5 && j == 5 && k == 5", block]; InfoVariable[name_] := @@ -789,47 +822,166 @@ insertFile[name_] := Close[istream]; contents]; +vectoriseExpression[exprp_] := + Module[{isNotMinusOneQ, isNotTimesMinusOneQ, fmaRules, isNotKneg, arithRules, undoRules, expr}, + expr = exprp; + + (* Constants *) + (* expr = expr /. xx_Integer/; xx!=-1 :> ToReal[xx]; *) + expr = expr /. xx_Integer -> ToReal[xx]; + expr = expr /. xx_Real -> ToReal[xx]; + removeToRealRules = { + - ToReal[xx_] -> ToReal[- xx], + ToReal[xx_] + ToReal[yy_] -> ToReal[xx + yy], + ToReal[xx_] * ToReal[yy_] -> ToReal[xx * yy], + ToReal[xx_] == ToReal[yy_] -> ToReal[xx == yy], + ToReal[xx_] != ToReal[yy_] -> ToReal[xx != yy], + pow[xx_, ToReal[power_]] -> pow[xx, power], + (* keep the conditional expression scalar *) + IfThen[ToReal[xx_], yy_, zz_] -> IfThen[xx, yy, zz] + }; + expr = expr //. removeToRealRules; + + (* Replace all operators and functions *) + (* kneg, kadd, ksub, kmul, kdiv *) + isNotKneg[n_] := ! MatchQ[n,kneg[_]]; + arithRules = { + - xx_ -> kneg[xx], + xx_ * yy_ -> kmul[xx,yy], + xx_ / yy_ -> kdiv[xx,yy], + xx_ + yy_ -> kadd[xx,yy], + xx_ - yy_ -> ksub[xx,yy], + kmul[-1,xx_] -> kneg[xx], + kmul[xx_,-1] -> kneg[xx], + kmul[ToReal[-1],xx_] -> kneg[xx], + kmul[xx_,ToReal[-1]] -> kneg[xx], + (* kmul[xx_,INV[yy_]] -> kdiv[xx,yy], *) + (* kmul[INV[xx_],yy_] -> kdiv[yy,xx], *) + kdiv[xx_,kdiv[yy_,zz_]] -> kdiv[kmul[xx,zz],yy], + kdiv[kdiv[xx_,yy_],zz_] -> kdiv[xx,kmul[yy,zz]], + kmul[kneg[xx_],yy_] -> kneg[kmul[xx,yy]], + kmul[xx_,kneg[yy_]] -> kneg[kmul[xx,yy]], + kdiv[kneg[xx_],yy_] -> kneg[kdiv[xx,yy]], + kdiv[xx_,kneg[yy_]] -> kneg[kdiv[xx,yy]], + kadd[kneg[xx_],yy_] -> ksub[yy,xx], + ksub[kneg[xx_],yy_] -> kneg[kadd[xx,yy]], + kadd[xx_,kneg[yy_]] -> ksub[xx,yy], + ksub[xx_,kneg[yy_]] -> kadd[xx,yy], + kneg[ksub[xx_,yy_]] -> ksub[yy,xx], + Abs[xx_] -> kfabs[xx], + Log[xx_] -> klog[xx], + fabs[xx_] -> kfabs[xx], + fmax[xx_,yy_] -> kfmax[xx,yy], + fmin[xx_,yy_] -> kfmin[xx,yy], + sqrt[xx_] -> ksqrt[xx], + exp[xx_] -> kexp[xx], + log[xx_] -> klog[xx], + pow[xx_,yy_] -> kpow[xx,yy], + kneg[kneg[xx_]] -> xx, + kfabs[kneg[xx_]] -> kfabs[xx], + kfnabs[kneg[xx_]] -> kfnabs[xx], + kneg[kfabs[xx_]] -> kfnabs[xx], + kneg[kfnabs[xx_]] -> kfabs[xx] + }; + expr = expr //. arithRules; + + (* Undo some transformations *) + undoRules = { + IfThen[kmul[xx_,yy_], aa_, bb_] -> IfThen[xx*yy, aa, bb], + IfThen[kmul[xx_,yy_] != zz_, aa_, bb_] -> IfThen[xx*yy!=zz, aa, bb], + ToReal[kneg[xx_]] -> ToReal[-xx], + ToReal[kadd[xx_,yy_]] -> ToReal[xx+yy], + ToReal[ksub[xx_,yy_]] -> ToReal[xx-yy], + ToReal[kmul[xx_,yy_]] -> ToReal[xx*yy], + ToReal[xx_*kadd[yy_,zz_]] -> ToReal[xx*(yy+zz)], + kpow[xx_, kneg[power_]] -> kpow[xx, -power] + }; + expr = expr //. undoRules; + + (* FMA (fused multiply-add) instructions *) + (* kmadd (x,y,z) = xy+z + kmsub (x,y,z) = xy-z + knmadd(x,y,z) = -(xy+z) + knmsub(x,y,z) = -(xy-z) *) + fmaRules = { + kadd[kmul[xx_,yy_],zz_] -> kmadd[xx,yy,zz], + kadd[zz_,kmul[xx_,yy_]] -> kmadd[xx,yy,zz], + ksub[kmul[xx_,yy_],zz_] -> kmsub[xx,yy,zz], + ksub[zz_,kmul[xx_,yy_]] -> knmsub[xx,yy,zz], + kneg[kmadd [xx_,yy_,zz_]] -> knmadd[xx,yy,zz], + kneg[kmsub [xx_,yy_,zz_]] -> knmsub[xx,yy,zz], + kneg[knmadd[xx_,yy_,zz_]] -> kmadd [xx,yy,zz], + kneg[knmsub[xx_,yy_,zz_]] -> kmsub [xx,yy,zz] + (* we could match this and similar patterns + kmul[xx_, kadd[yy_, ToReal[+1]]] -> kmadd[xx, yy, xx], + kmul[xx_, kadd[yy_, ToReal[-1]]] -> kmsub[xx, yy, xx], + *) + }; + expr = expr //. fmaRules; + + Return[expr]]; + (* Take an expression x and replace occurrences of Powers with the C macros SQR, CUB, QAD *) -ReplacePowers[x_] := +ReplacePowers[expr_, vectorise_] := Module[{rhs}, - rhs = x /. Power[xx_, -1] -> INV[xx]; + rhs = expr /. Power[xx_, -1] -> INV[xx]; If[SOURCELANGUAGE == "C", Module[{}, - rhs = rhs //. Power[xx_, 2] -> SQR[xx]; - rhs = rhs //. Power[xx_, 3] -> CUB[xx]; - rhs = rhs //. Power[xx_, 4] -> QAD[xx]; + rhs = rhs /. Power[xx_, 2 ] -> SQR[xx]; + rhs = rhs /. Power[xx_, 3 ] -> CUB[xx]; + rhs = rhs /. Power[xx_, 4 ] -> QAD[xx]; + rhs = rhs /. Power[xx_, -2 ] -> INV[SQR[xx]]; + rhs = rhs /. Power[xx_, 1/2] -> sqrt[xx]; + rhs = rhs /. Power[xx_, -1/2] -> INV[sqrt[xx]]; + rhs = rhs /. Power[xx_, 0.5] -> sqrt[xx]; + rhs = rhs /. Power[xx_, -0.5] -> INV[sqrt[xx]]; + + (* + rhs = rhs /. 1/2 -> khalf + rhs = rhs /. -1/2 -> -khalf; - rhs = rhs //. xx_/2 -> khalf xx; - rhs = rhs //. (-1/2) -> -khalf; + rhs = rhs /. 1/3 -> kthird; + rhs = rhs /. -1/3 -> -kthird; - rhs = rhs //. xx_/3 -> kthird xx; - rhs = rhs //. (-1/3) -> -kthird; + rhs = rhs /. 2/3 -> ktwothird; + rhs = rhs /. -2/3 -> -ktwothird; - rhs = rhs //. 2/3 -> ktwothird; - rhs = rhs //. (-2/3) -> -ktwothird; + rhs = rhs /. 4/3 -> kfourthird; + rhs = rhs /. -4/3 -> -kfourthird; - rhs = rhs //. 4/3 -> kfourthird; - rhs = rhs //. (-4/3) -> -kfourthird; + rhs = rhs /. 8/3 -> keightthird; + rhs = rhs /. -8/3 -> -keightthird; + *) - rhs = rhs //. 8/3 -> keightthird; - rhs = rhs //. (-8/3) -> -keightthird; + (* Avoid rational numbers *) + rhs = rhs /. Rational[xx_,yy_] :> N[xx/yy, 30]; - rhs = rhs //. xx_ y_ + xx_ z_ -> xx(y+z); + rhs = rhs //. IfThen[cond1_,xx1_,yy1_] + IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :> IfThen[cond1, Simplify[ xx1 + xx2], Simplify[ yy1 + yy2]]; + rhs = rhs //. ff1_ IfThen[cond1_,xx1_,yy1_] + IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :> IfThen[cond1, Simplify[ff1 xx1 + xx2], Simplify[ff1 yy1 + yy2]]; + rhs = rhs //. IfThen[cond1_,xx1_,yy1_] + ff2_ IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :> IfThen[cond1, Simplify[ xx1 + ff2 xx2], Simplify[ yy1 + ff2 yy2]]; + rhs = rhs //. ff1_ IfThen[cond1_,xx1_,yy1_] + ff2_ IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :> IfThen[cond1, Simplify[ff1 xx1 + ff2 xx2], Simplify[ff1 yy1 + ff2 yy2]]; - rhs = rhs //. Power[E, power_] -> exp[power]; - rhs = rhs //. Power[xx_, 0.5] -> sqrt[xx]; + (* Is this still a good idea when FMA instructions are used? *) + rhs = rhs //. xx_ yy_ + xx_ zz_ -> xx (yy+zz); + rhs = rhs //. xx_ yy_ - xx_ zz_ -> xx (yy-zz); + + rhs = rhs /. Power[E, power_] -> exp[power]; (* there have been some problems doing the Max/Min replacement via the preprocessor for C, so we do it here *) - rhs = rhs //. Max[xx_, yy_] -> fmax[xx, yy]; - rhs = rhs //. Min[xx_, yy_] -> fmin[xx, yy]; + rhs = rhs /. Max[xx_, yy_] -> fmax[xx, yy]; + rhs = rhs /. Min[xx_, yy_] -> fmin[xx, yy]; + + rhs = rhs /. Power[xx_, power_] -> pow[xx, power]; - rhs = rhs //. Power[xx_, power_] -> pow[xx, power]], + If[vectorise === True, + rhs = vectoriseExpression[rhs]]; + ], - rhs = rhs //. Power[xx_, power_] -> xx^power + rhs = rhs /. Power[xx_, power_] -> xx^power ]; (* Print[rhs//FullForm];*) rhs @@ -841,180 +993,6 @@ CFormHideStrings[x_, opts___] := StringReplace[ToString[CForm[x,opts]], "\"" -> -(* Debug output *) -DebugCSE = True; - -(* Eliminate common subexpressions in a code sequence *) -CSE[code_] := Module[ - {expr, optexpr, - decomposed, locals, block, - block1, block2, temps1, stmts1, stmts2, stmts3, - replacevar, - stmts4, - stmts5, stmts6, stmts7}, - if [DebugCSE, Print["code\n", code, "\nendcode\n"]]; - - (* The code is passed in as list of {lhs,rhs} tuples. Turn this - list into a single expression, so that it can be optimised. *) - expr = code //. {a_, b__} -> CSequence[a, {b}] - //. {a_} -> a - //. (a_ -> b_) -> CAssign[a, b]; - If [DebugCSE, Print["expr\n", expr, "\nendexpr\n"]]; - - (* Optimise this expression *) - optexpr = Experimental`OptimizeExpression[expr]; - If [DebugCSE, Print["optexpr\n", optexpr, "\nendoptexpr\n"]]; - - (* This expression is a Mathematica expression. Decompose it into - the set of newly introduced local variables and the optimised - expression itself. *) - decomposed = - ReleaseHold[(Hold @@ optexpr) - /. Verbatim[Block][vars_, seq_] :> {vars, Hold[seq]}]; - - If[decomposed[[0]] =!= List, - (* If the optimiser didn't create a Block expression, we assume it - didn't do anything useful and return the original. *) - code, - - {locals, block} = decomposed; - If [DebugCSE, Print["locals\n", locals, "\nendlocals\n"]]; - If [DebugCSE, Print["block\n", block, "\nendblock\n"]]; - - block1 = block /. Hold[CompoundExpression[seq__]] :> Hold[{seq}]; - If [DebugCSE, Print["block1\n", block1, "\nendblock1\n"]]; - block2 = First[block1 //. Hold[{a___Hold, b_, c___}] - /; Head[Unevaluated[b]] =!= Hold - :> Hold[{a, Hold[b], c}]]; - If [DebugCSE, Print["block2\n", block2, "\nendblock2\n"]]; - - (* Temporaries, including a fake declaration for them *) - temps1 = Most[block2] //. Hold[lhs_ = rhs_] -> CAssign[CDeclare[lhs], rhs]; - If [DebugCSE, Print["temps1\n", temps1, "\nendtemps1\n"]]; - - (* Expression *) - stmts1 = ReleaseHold[Last[block2]]; - If [DebugCSE, Print["stmts1\n", stmts1, "\nendstmts1\n"]]; - - (* Turn CSequence back into a list *) - stmts2 = Flatten[{stmts1} //. CSequence[a_,b_] -> {a,b}]; - If [DebugCSE, Print["stmts2\n", stmts2, "\nendstmts2\n"]]; - - (* Combine temporaries and expression *) - stmts3 = Join[temps1, stmts2]; - If [DebugCSE, Print["stmts3\n", stmts3, "\nendstmts3\n"]]; - - (* Replace the internal names of the newly generated temporaries - with legal C names *) - replacevar = - Rule @@@ Transpose[{(*ToString[CForm[#]] & /@*) locals, - Symbol[ - StringReplace[StringReplace[ToString[#], {__ ~~ "`" ~~ a_ :> a}], - "$" -> "T"]] & /@ locals}]; - If [DebugCSE, Print["replacevar\n", replacevar, "\nendreplacevar\n"]]; - - stmts4 = stmts3 //. replacevar; - If [DebugCSE, Print["stmts4\n", stmts4, "\nendstmts4\n"]]; - - (* Sort statements topologically *) -(* - stmts5 = stmts4; -*) - If [DebugCSE, Print["A\n"]]; - stmts5 = - Module[{debug, - tmpVars, newVars, i, - stmtsLeft, stmtsDone, - lhs, rhs, any, contains, containsAny, - canDoStmts, cannotDoStmts, - selfStmts, selfVars, allVars, nonSelfVars}, - debug = False; - stmtsLeft = stmts4; - If [DebugCSE, Print["B\n"]]; - stmtsDone = {}; - If [DebugCSE, Print["C\n"]]; - (* lhs[x_] := x[[1]]; *) - lhs[x_] := x /. (CAssign[lhs_, rhs_] -> lhs); - If [DebugCSE, Print["D\n"]]; - (* rhs[x_] := x[[2]]; *) - rhs[x_] := x /. (CAssign[lhs_, rhs_] -> rhs); - If [DebugCSE, Print["E\n"]]; - (* any[xs_] := Fold[Or, False, xs]; *) - any[xs_] := MemberQ[xs, True]; - If [DebugCSE, Print["F\n"]]; - (* contains[e_, x_] := (e /. x -> {}) =!= e; *) - (* contains[e_, x_] := Count[{e}, x, Infinity] > 0; *) - contains[e_, x_] := MemberQ[{e}, x, Infinity]; - If [DebugCSE, Print["G\n"]]; - containsAny[e_, xs_] := any[Map[contains[e,#]&, xs]]; - If [DebugCSE, Print["H\n"]]; - getVars[stmts_] := Map[lhs, stmts] //. (CDeclare[lhs_] -> lhs); - - (* Rename temporary variables deterministically *) - tmpVars = Select[getVars[stmtsLeft], - StringMatchQ[ToString[#], "TT"~~__]&]; - newVars = Table[Symbol["T"<>ToString[1000000+i]], - {i, 1, Length[tmpVars]}]; - stmtsLeft = stmtsLeft /. MapThread[(#1->#2)&, {tmpVars, newVars}]; - - allVars = getVars[stmtsLeft]; - While[stmtsLeft =!= {}, - If[debug, Print["stmtsLeft = \n", stmtsLeft]]; - If[debug, Print["stmtsDone = \n", stmtsDone]]; - allVars = getVars[stmtsLeft]; - If[debug, Print["allVars = \n", allVars]]; - canDoStmts = - Select[stmtsLeft, Not[containsAny[rhs[#], allVars]] &]; - cannotDoStmts = - Select[stmtsLeft, containsAny[rhs[#], allVars] &]; - If[debug, Print["canDoStmts = \n", canDoStmts]]; - If[debug, Print["cannotDoStmts = \n", cannotDoStmts]]; - If[False && canDoStmts == {}, - (* Handle assignment where LHS and RHS access the same variables - (hopefully without taking derivatives!) *) - selfStmts = Select[stmtsLeft, contains[rhs[#], lhs[#]]]; - selfVars = getVars[selfStmts]; - nonSelfVars = Select[allVars, Not[contains[selfVars, #]] &]; - canDoStmts = - Select[stmtsLeft, Not[containsAny[rhs[#], nonSelfVars]] &]; - cannotDoStmts = - Select[stmtsLeft, containsAny[rhs[#], nonSelfVars] &]; - If[debug, Print["nonself/canDoStmts = \n", canDoStmts]]; - If[debug, Print["nonself/cannotDoStmts = \n", cannotDoStmts]]; - ]; - If[canDoStmts == {}, - (* Accept the first statement *) - canDoStmts = {First[stmtsLeft]}; - cannotDoStmts = Rest[stmtsLeft]; - If[debug, Print["takeone/canDoStmts = \n", canDoStmts]]; - If[debug, Print["takeone/cannotDoStmts = \n", cannotDoStmts]]; - ]; - If[canDoStmts == {}, ThrowError["canDoStmts == {}"]]; - stmtsDone = Join[stmtsDone, canDoStmts]; - If [DebugCSE, Print["I\n"]]; - stmtsLeft = cannotDoStmts; - If [DebugCSE, Print["J\n"]]; - ]; - If[debug, Print["stmtsLeft\n", stmtsLeft]]; - If[debug, Print["stmtsDone\n", stmtsDone]]; - stmtsDone]; - If [DebugCSE, Print["Z\n"]]; - - (* Turn CAssign statements back into (->) tuples *) - stmts6 = stmts5 //. CAssign[lhs_,rhs_] -> (lhs -> rhs); - If [DebugCSE, Print["stmts6\n", stmts6, "\nendstmts6\n"]]; - - (* Turn CDeclare statements into "faked" declarations *) - stmts7 = stmts6 - //. CDeclare[var_] - :> "CCTK_REAL const " <> - StringReplace[ToString[var], __ ~~ "`" -> ""]; - If [DebugCSE, Print["stmts7\n", stmts7, "\nendstmts7\n"]]; - - stmts7 - ] -]; - Quote[x_] := {"\"", x, "\""}; End[]; diff --git a/Tools/CodeGen/Differencing.m b/Tools/CodeGen/Differencing.m index ff6db9c..ab3d478 100644 --- a/Tools/CodeGen/Differencing.m +++ b/Tools/CodeGen/Differencing.m @@ -142,6 +142,7 @@ DZero::usage = ""; shift::usage = ""; spacing::usage = ""; ComponentDerivativeOperatorStencilWidth::usage = ""; +CheckStencil::usage = ""; Begin["`Private`"]; @@ -155,7 +156,7 @@ DZero[n_] := (DPlus[n] + DMinus[n])/2; (* User API *) (*************************************************************) -CreateDifferencingHeader[derivOps_, zeroDims_] := +CreateDifferencingHeader[derivOps_, zeroDims_, vectorise_] := Module[{componentDerivOps, dupsRemoved, expressions, componentDerivOps2, zeroDimRules, derivOps2, pDefs}, Map[DerivativeOperatorVerify, derivOps]; @@ -167,14 +168,12 @@ CreateDifferencingHeader[derivOps_, zeroDims_] := dupsRemoved = RemoveDuplicateRules[componentDerivOps2]; - mDefPairs = Map[ComponentDerivativeOperatorMacroDefinition, dupsRemoved]; + mDefPairs = Map[ComponentDerivativeOperatorMacroDefinition[#, vectorise] &, dupsRemoved]; pDefs = Union[Flatten[Map[First, mDefPairs]]]; expressions = Flatten[Map[#[[2]]&, mDefPairs]]; -(* expressions = Flatten[Map[ComponentDerivativeOperatorInlineDefinition, dupsRemoved]];*) - - {pDefs,Map[{#, "\n"} &, expressions]}]; + {pDefs, Map[{#, "\n"} &, expressions]}]; ordergfds[_[v1_,___], _[v2_,___]] := Order[v1,v2] != -1; @@ -205,6 +204,16 @@ ReplaceDerivatives[derivOps_, expr_, precompute_] := rules = Map[# :> evaluateDerivative[#] &, gfds]]; expr /. rules]; +(* Generate code to ensure that there are sufficient ghost and + boundary points for the passed derivative operators used in eqs *) +CheckStencil[derivOps_, eqs_, name_] := + Module[{gfds, rgzList, rgz}, + gfds = Map[GridFunctionDerivativesInExpression[{#}, eqs] &, derivOps]; + rgzList = MapThread[If[Length[#2] > 0, DerivativeOperatorStencilWidth[#1], {0,0,0}] &, {derivOps, gfds}]; + If[Length[rgzList] === 0, Return[{}]]; + rgz = Map[Max, Transpose[rgzList]]; + If[Max[rgz] == 0, {}, + {"GenericFD_EnsureStencilFits(cctkGH, ", Quote@name, ", ", Riffle[rgz,", "], ");\n"}]]; (*************************************************************) (* Misc *) @@ -212,15 +221,18 @@ ReplaceDerivatives[derivOps_, expr_, precompute_] := PrecomputeDerivative[d:pd_[gf_, inds___]] := Module[{}, - DeclareAssignVariable["CCTK_REAL", GridFunctionDerivativeName[d], evaluateDerivative[d]]]; + DeclareAssignVariable[DataType[], GridFunctionDerivativeName[d], evaluateDerivative[d]]]; evaluateDerivative[d:pd_[gf_, inds___]] := Module[{macroname}, macroName = ComponentDerivativeOperatorMacroName[pd[inds] -> expr]; - Return[ToString[macroName] <> "(" <> ToString[gf] <> ", i, j, k)"]]; + (* Return[ToString[macroName] <> "(" <> ToString[gf] <> ", i, j, k)"] *) + (* Return[ToString[macroName] <> "(" <> ToString[gf] <> ")"] *) + Return[ToString[macroName] <> "(&" <> ToString[gf] <> "[index])"] + ]; DeclareDerivative[d:pd_[gf_, inds___]] := - DeclareVariable[GridFunctionDerivativeName[d], "// CCTK_REAL"]; + DeclareVariable[GridFunctionDerivativeName[d], "// CCTK_REAL_VEC"]; (*************************************************************) @@ -251,8 +263,8 @@ sbpMacroDefinition[macroName_, d_] := FlattenBlock[{"#define ", macroName, "(u,i,j,k) (sbp_deriv_" <> ds <> "(i,j,k,sbp_" <> l <> "min,sbp_" <> l <> "max,d" <> ds <> ",u,q" <> ds <> ",cctkGH))"}] ]; -ComponentDerivativeOperatorMacroDefinition[componentDerivOp:(name_[inds___] -> expr_)] := - Module[{macroName, rhs, rhs2, i = "i", j = "j", k = "k", spacings, spacings2, pat, ss, num, den, newnum, signModifier, quotient, liName, rhs3, rhs4}, +ComponentDerivativeOperatorMacroDefinition[componentDerivOp:(name_[inds___] -> expr_), vectorise_] := + Module[{macroName, rhs, i = "i", j = "j", k = "k", spacings, spacings2, pat, ss, num, den, newnum, signModifier, quotient, liName, finalDef}, macroName = ComponentDerivativeOperatorMacroName[componentDerivOp]; @@ -265,23 +277,24 @@ ComponentDerivativeOperatorMacroDefinition[componentDerivOp:(name_[inds___] -> e If[expr === SBPDerivative[3], Return[sbpMacroDefinition[macroName, 3]]]; - rhs = DifferenceGF[expr, i, j, k]; + rhs = DifferenceGF[expr, i, j, k, vectorise]; +(* Print["rhs1 == ", FullForm[rhs]];*) spacings = {spacing[1] -> 1/"dxi", spacing[2] -> 1/"dyi", spacing[3] -> 1/"dzi"}; spacings2 = {spacing[1] -> "dx", spacing[2] -> "dy", spacing[3] -> "dz"}; - rhs2 = FullSimplify[rhs]; + rhs = FullSimplify[rhs]; -(* Print["rhs2 == ", FullForm[rhs2]];*) +(* Print["rhs2 == ", FullForm[rhs]];*) pat = (Times[spInExpr:(Power[spacing[_],_]...), (Rational[x_,y_])..., rest__]) | (rest__); (* Print["pat == ", pat//FullForm]; *) - If[MatchQ[rhs2, pat], + If[MatchQ[rhs, pat], (* Print["matches!"];*) - ss = Times[rhs2 /. pat -> spInExpr]; + ss = Times[rhs /. pat -> spInExpr]; (* Print["ss == ", ss];*) - num = rhs2 /. pat -> x; - den = rhs2 /. pat -> y; + num = rhs /. pat -> x; + den = rhs /. pat -> y; (* Print["num == ", num]; Print["den == ", den];*) If[{num, 1, 2} === {1, 2},(* Print["SEQ!"]; *) newnum = 1; den=1; signModifier = "", @@ -307,84 +320,117 @@ ComponentDerivativeOperatorMacroDefinition[componentDerivOp:(name_[inds___] -> e liName = "p" <> signModifier <> quotient <> ToString[Apply[SequenceForm,Simplify[1/(ss /. spacings2)],{0,Infinity}]]; (* Print["liName == ", liName];*) - rhs3 = rhs2 /. pat -> Times[liName, rest], - ThrowError["Partial derivative operator definition ", rhs2, - " does not match pattern ", pat]; - rhs3 = rhs2]; - -(* Print["rhs3 == ", rhs3];*) + rhs = rhs /. pat -> Times[liName, rest], +(* + rhs = (rhs /. pat -> Times[liName, rest]) / If[vectorise, liName, 1], (* horrible *) +*) +(* Print["!!!!!!!!DOES NOT MATCH!!!!!!!!!"];*) + rhs = rhs]; - pDefs = {{liName -> CFormHideStrings[ReplacePowers[num / den ss /. spacings2]]}}; +(* Print["rhs3 == ", FullForm[rhs]];*) -(* rhs4 = Factor[rhs3];*) + pDefs = {{liName -> CFormHideStrings[ReplacePowers[num / den ss /. spacings2, vectorise]]}}; - rhs4 = rhs3 //. (x_ a_ + x_ b_) -> x(a+b); - rhs5 = rhs4 //. (x_ a_ - x_ b_) -> x(a-b); +(* rhs = Factor[rhs];*) + rhs = rhs //. (x_ a_ + x_ b_) -> x (a+b); + rhs = rhs //. (x_ a_ - x_ b_) -> x (a-b); (* Print[componentDerivOp, ": "]; - Print[FullForm[rhs5]]; + Print[FullForm[rhs]]; Print[""];*) - rhs6 = CFormHideStrings[ReplacePowers[rhs5 /. spacings]]; - {pDefs, FlattenBlock[{"#define ", macroName, "(u,i,j,k) ", "(", rhs6, ")"}]}]; - -ComponentDerivativeOperatorInlineDefinition[componentDerivOp:(name_[inds___] -> expr_)] := - Module[{inlineName, rhs, rhs2, i = "i", j = "j", k = "k", spacings}, - - inlineName = ComponentDerivativeOperatorMacroName[componentDerivOp]; - - rhs = DifferenceGF[expr, i, j, k]; -(* rhs = DifferenceGFInline[expr, i, j, k];*) - spacings = {spacing[1] -> 1/"dxi", spacing[2] -> 1/"dyi", spacing[3] -> 1/"dzi"}; - rhs2 = CFormHideStrings[FullSimplify[ReplacePowers[rhs /. spacings]]]; - - DefineFunction[inlineName, "static inline CCTK_REAL", - "CCTK_REAL *u, int i, int j, int k", - {"return ", rhs2, ";\n"}]]; + rhs = CFormHideStrings[ReplacePowers[rhs /. spacings, vectorise]]; + (* Print["rhs=",FullForm[rhs]]; *) + (* {pDefs, FlattenBlock[{"#define ", macroName, "(u,i,j,k) ", "(", rhs, ")"}]} *) + finalDef = + If[vectorise, + {pDefs, FlattenBlock[{ + "#ifndef KRANC_DIFF_FUNCTIONS\n", + (* default, differencing operators are macros *) +(* + "# define ", macroName, "(u) ", "(kmul(", liName, ",", rhs, "))\n", +*) + "# define ", macroName, "(u) ", "(", rhs, ")\n", + "#else\n", + (* new, differencing operators are static functions *) +(* + "# define ", macroName, "(u) ", "(kmul(", liName, ",", macroName, "_impl((u),dj,dk)))\n", + "static CCTK_REAL_VEC ", macroName, "_impl(CCTK_REAL const* restrict const u, ptrdiff_t const dj, ptrdiff_t const dk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;\n", + "static CCTK_REAL_VEC ", macroName, "_impl(CCTK_REAL const* restrict const u, ptrdiff_t const dj, ptrdiff_t const dk)\n", + If[StringMatchQ[rhs, RegularExpression[".*\\bdir\\d\\b.*"]], + { "{ return ToReal(1e30); /* ERROR */ }\n" }, + { "{ return ", rhs, "; }\n" }], +*) +(* + "# define ", macroName, "(u) ", "(kmul(", liName, ",", macroName, "_impl(u,cdj,cdk)))\n", + "static CCTK_REAL_VEC ", macroName, "_impl(CCTK_REAL const* restrict const u, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;\n", + "static CCTK_REAL_VEC ", macroName, "_impl(CCTK_REAL const* restrict const u, ptrdiff_t const cdj, ptrdiff_t const cdk)\n", + (* We cannot handle dirN, + so we punt on all expressions that contain dirN *) + If[StringMatchQ[rhs, RegularExpression[".*\\bdir\\d\\b.*"]], + { "{ return ToReal(1e30); /* ERROR */ }\n" }, + { "{\n", + " ptrdiff_t const cdi=sizeof(CCTK_REAL);\n", + " return ", rhs, ";\n", + "}\n" }], +*) + "# define ", macroName, "(u) ", "(", macroName, "_impl(u,", liName, ",cdj,cdk))\n", + "static CCTK_REAL_VEC ", macroName, "_impl(CCTK_REAL const* restrict const u, CCTK_REAL_VEC const ", liName, ", ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;\n", + "static CCTK_REAL_VEC ", macroName, "_impl(CCTK_REAL const* restrict const u, CCTK_REAL_VEC const ", liName, ", ptrdiff_t const cdj, ptrdiff_t const cdk)\n", + (* We cannot handle dirN, + so we punt on all expressions that contain dirN *) + If[StringMatchQ[rhs, RegularExpression[".*\\bdir\\d\\b.*"]], + { "{ return ToReal(1e30); /* ERROR */ }\n" }, + { "{\n", + " ptrdiff_t const cdi=sizeof(CCTK_REAL);\n", + " return ", rhs, ";\n", + "}\n" }], + "#endif\n" + }]}, + + {pDefs, FlattenBlock[{"#define ", macroName, "(u) ", "(", rhs, ")"}]}]; + finalDef +]; ComponentDerivativeOperatorMacroName[componentDerivOp:(name_[inds___] -> expr_)] := Module[{stringName}, stringName = StringJoin[Map[ToString, Join[{name}, {inds}]]]; stringName]; - - +DerivativeOperatorStencilWidth[derivOp_] := + Map[Max, Transpose[Map[ComponentDerivativeOperatorStencilWidth, DerivativeOperatorToComponents[derivOp]]]]; ComponentDerivativeOperatorStencilWidth[componentDerivOp:(name_[inds___] -> expr_)] := - Module[{cases, nx, ny, nz}, - cases = Union[Flatten[Cases[{expr}, shift[_] | Power[shift[_],_], Infinity]]]; - Print[cases]; - - nx = Exponent[op, shift[1]]; - ny = Exponent[op, shift[2]]; - nz = Exponent[op, shift[3]]; - - ]; - - - + Module[{cases, nx, ny, nz, result}, + result = Table[ + cases = Union[Flatten[Cases[{expr}, shift[d] | Power[shift[d],_], Infinity]]]; + ns = Map[Exponent[#, shift[d]] &, cases]; + If[Length[ns] == 0, 0, Max[Abs[ns]]], {d, 1, 3}]; + + (* We do not know the run-time value of any shorthands used in + operator definitions. In all the current known cases, this + will be a "direction" which is +/- 1. In future, the + differencing mechanism will support shorthand arguments to + operators and this hack can be removed. *) + result = Replace[result, _Symbol -> 1, {-1}]; + + If[!And@@Map[NumericQ, result], + Throw["Stencil width is not numeric in "<>ToString[componentDerivOp]]]; + result]; (* Farm out each term of a difference operator *) -DifferenceGF[op_, i_, j_, k_] := - Module[{expanded}, - expanded = Expand[op]; - - If[Head[expanded] === Plus, - Apply[Plus, Map[DifferenceGFTerm[#, i, j, k] &, expanded]], - DifferenceGFTerm[expanded, i, j, k]]]; - -DifferenceGFInline[op_, i_, j_, k_] := +DifferenceGF[op_, i_, j_, k_, vectorise_] := Module[{expanded}, expanded = Expand[op]; If[Head[expanded] === Plus, - Apply[Plus, Map[DifferenceGFTermInline[#, i, j, k] &, expanded]], + Apply[Plus, Map[DifferenceGFTerm[#, i, j, k, vectorise] &, expanded]], DifferenceGFTerm[expanded, i, j, k]]]; (* Return the fragment of a macro definition for defining a derivative operator *) -DifferenceGFTerm[op_, i_, j_, k_] := +DifferenceGFTerm[op_, i_, j_, k_, vectorise_] := Module[{nx, ny, nz, remaining}, If[op === 0, @@ -409,10 +455,52 @@ DifferenceGFTerm[op_, i_, j_, k_] := "(int)(" <> ToString[CFormHideStrings[j+ny]] <> ")," <> "(int)(" <> ToString[CFormHideStrings[k+nz]] <> "))]", *) - remaining "(u)[index" <> - "+di*(" <> ToString[CFormHideStrings[nx]] <> ")" <> +(* + remaining "vec_loadu_maybe" <> + "(" <> ToString[CFormHideStrings[nx]] <> "," <> + "(u)[index" <> + "+di*(" <> ToString[CFormHideStrings[nx]] <> ")" <> + "+dj*(" <> ToString[CFormHideStrings[ny]] <> ")" <> + "+dk*(" <> ToString[CFormHideStrings[nz]] <> ")])", +*) +(* + remaining "vec_loadu_maybe(" <> ToString[CFormHideStrings[nx]] <> "," <> + "(u)[(" <> ToString[CFormHideStrings[nx]] <> ")" <> + "+dj*(" <> ToString[CFormHideStrings[ny]] <> ")" <> + "+dk*(" <> ToString[CFormHideStrings[nz]] <> ")])", +*) + + If[vectorise, +(* + remaining "vec_loadu_maybe3" <> + "(" <> ToString[CFormHideStrings[nx /. {dir1->1, dir2->1, dir3->1}]] <> "," <> + ToString[CFormHideStrings[ny /. {dir1->1, dir2->1, dir3->1}]] <> "," <> + ToString[CFormHideStrings[nz /. {dir1->1, dir2->1, dir3->1}]] <> "," <> + "(u)[(" <> ToString[CFormHideStrings[nx]] <> ")" <> + "+dj*(" <> ToString[CFormHideStrings[ny]] <> ")" <> + "+dk*(" <> ToString[CFormHideStrings[nz]] <> ")])", +*) + remaining "vec_loadu_maybe3" <> + "(" <> ToString[CFormHideStrings[nx /. {dir1->1, dir2->1, dir3->1}]] <> "," <> + ToString[CFormHideStrings[ny /. {dir1->1, dir2->1, dir3->1}]] <> "," <> + ToString[CFormHideStrings[nz /. {dir1->1, dir2->1, dir3->1}]] <> "," <> + "*(CCTK_REAL const*)&((char const*)(u))" <> + "[cdi*(" <> ToString[CFormHideStrings[nx]] <> ")" <> + "+cdj*(" <> ToString[CFormHideStrings[ny]] <> ")" <> + "+cdk*(" <> ToString[CFormHideStrings[nz]] <> ")])", + + remaining "(u)[" <> + "di*(" <> ToString[CFormHideStrings[nx]] <> ")" <> "+dj*(" <> ToString[CFormHideStrings[ny]] <> ")" <> - "+dk*(" <> ToString[CFormHideStrings[nz]] <> ")]", + "+dk*(" <> ToString[CFormHideStrings[nz]] <> ")]"], + + +(* + remaining "vec_loadu" <> + "(u[(" <> ToString[CFormHideStrings[nx]] <> ")" <> + "+dj*(" <> ToString[CFormHideStrings[ny]] <> ")" <> + "+dk*(" <> ToString[CFormHideStrings[nz]] <> ")])", +*) (* remaining "(u)[CCTK_GFINDEX3D(cctkGH,floor((" <> ToString[CFormHideStrings[i+nx]] <> ")+0.5),floor((" <> @@ -422,27 +510,6 @@ DifferenceGFTerm[op_, i_, j_, k_] := remaining "u(" <> ToString[FortranForm[i+nx]] <> "," <> ToString[FortranForm[j+ny]] <> "," <> ToString[FortranForm[k+nz]] <> ")"] ]; -(* Return the fragment of a function definition for defining a derivative - operator *) -DifferenceGFTermInline[op_, i_, j_, k_] := - Module[{nx, ny, nz, remaining}, - - If[op === 0, - Return[0]]; - - nx = Exponent[op, shift[1]]; - ny = Exponent[op, shift[2]]; - nz = Exponent[op, shift[3]]; - - remaining = op / (shift[1]^nx) / (shift[2]^ny) / (shift[3]^nz); - - If[Cases[{remaining}, shift[_], Infinity] != {}, - ThrowError["Could not parse difference operator:", op]]; - - remaining "(u)[CCTK_GFINDEX3D(cctkGH," <> ToString[CFormHideStrings[i+nx]] <> "," <> - ToString[CFormHideStrings[j+ny]] <> "," <> ToString[CFormHideStrings[k+nz]] <> ")]" - ]; - DerivativeOperatorGFDs[gf_]; diff --git a/Tools/CodeGen/Interface.m b/Tools/CodeGen/Interface.m index 29c198f..5f8f613 100644 --- a/Tools/CodeGen/Interface.m +++ b/Tools/CodeGen/Interface.m @@ -73,7 +73,7 @@ CreateKrancInterface[nonevolvedGroups_, evolvedGroups_, rhsGroups_, groups_, Module[{registerEvolved, (*registerConstrained,*) nonevolvedGroupStructures, evolvedGroupStructures, rhsGroupStructures, - groupStructures, interface}, + groupStructures, interface, getMap}, VerifyGroupNames[nonevolvedGroups]; VerifyGroupNames[evolvedGroups]; VerifyGroupNames[rhsGroups]; @@ -105,6 +105,13 @@ CreateKrancInterface[nonevolvedGroups_, evolvedGroups_, rhsGroups_, groups_, ArgString -> "CCTK_POINTER_TO_CONST IN cctkGH, CCTK_INT IN dir, CCTK_INT IN nsize, CCTK_INT OUT ARRAY imin, CCTK_INT OUT ARRAY imax, CCTK_REAL OUT ARRAY q, CCTK_INT IN table_handle" }; + getMap = + { + Name -> "MultiPatch_GetMap", + Type -> "CCTK_INT", + ArgString -> "CCTK_POINTER_TO_CONST IN cctkGH" + }; + (* For each group declared in this thorn, we need an entry in the interface file. Each evolved group needs an associated rhs @@ -127,10 +134,11 @@ CreateKrancInterface[nonevolvedGroups_, evolvedGroups_, rhsGroups_, groups_, interface = CreateInterface[implementation, inheritedImplementations, Join[includeFiles, {CactusBoundary`GetIncludeFiles[]}, - If[OptionValue[UseLoopControl], {"loopcontrol.h"}, {}]], + If[OptionValue[UseLoopControl], {"loopcontrol.h"}, {}], + If[OptionValue[UseVectors], {"vectors.h"}, {}]], groupStructures, UsesFunctions -> - Join[{registerEvolved, (*registerConstrained,*) diffCoeff}, + Join[{registerEvolved, (*registerConstrained,*) diffCoeff, getMap}, CactusBoundary`GetUsedFunctions[]]]; Return[interface]]; diff --git a/Tools/CodeGen/Jacobian.m b/Tools/CodeGen/Jacobian.m new file mode 100644 index 0000000..31f6b67 --- /dev/null +++ b/Tools/CodeGen/Jacobian.m @@ -0,0 +1,140 @@ + +(* Copyright 2011 Ian Hinder + + This file is part of Kranc. + + Kranc is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + Kranc is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Kranc; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*) + +BeginPackage["Jacobian`", {"Errors`", "Helpers`", "Kranc`", "Differencing`", "MapLookup`", "CodeGen`", "KrancGroups`"}]; + +JacobianQ; +InsertJacobian; +CreateJacobianVariables; +JacobianGenericFDParameters; +JacobianSymbols; +JacobianGroups; +JacobianCheckGroups; +JacobianConditionalGridFunctions; + +Begin["`Private`"]; + +Options[JacobianQ] = ThornOptions; +JacobianQ[opts:OptionsPattern[]] := + Length[OptionValue[Jacobian]] > 0; + +(* Assign a shorthand containing the Jacobian multiplied by the passed + 1st derivative *) +jacobianShorthand[d:(deriv_[var_, i_])] := + Module[{}, + derivToJacDeriv[d] -> + IfThen["use_jacobian", Sum[Symbol["J"<>ToString[j]<>ToString[i]] deriv[var, j], {j, 1 3}], deriv[var, i]] + ]; + +(* Assign a shorthand containing the Jacobian multiplied by the passed + 2nd derivative *) +jacobianShorthand[d:(deriv_[var_, i_,j_])] := + Module[{ip,jp}, + {ip,jp} = Sort[{i,j}]; (* dJ is symmetric in the last two indices *) + derivToJacDeriv[d] -> + IfThen["use_jacobian", Sum[Symbol["dJ"<>ToString[a]<>ToString[ip]<>ToString[jp]] deriv[var, a], {a, 1 3}] + + Sum[Symbol["J"<>ToString[a]<>ToString[i]] Symbol["J"<>ToString[b]<>ToString[j]] deriv[var, a, b], {a, 1 3}, {b, 1, 3}], + deriv[var, i, j]] + ]; + +(* Convert a 1st derivative to a Jacobian-multiplied derivative *) +derivToJacDeriv[deriv_[var_, i_]] := + Symbol["Global`Jac"<>ToString[deriv]<>ToString[i]<>ToString[var]]; + +(* Convert a 2nd derivative to a Jacobian-multiplied derivative *) +derivToJacDeriv[deriv_[var_, i_, j_]] := + Symbol["Global`Jac"<>ToString[deriv]<>ToString[i]<>ToString[j]<>ToString[var]]; + +(* Given a calculation containing partial derivatives, return a + version of the calculation with all the partial derivatives multiplied + by the Jacobian *) +Options[InsertJacobian] = ThornOptions; +InsertJacobian[calc_List, opts:OptionsPattern[]] := + Module[{pdDefs, derivs, newShortDefs, newShorts, combinedShorts, combinedEqs, combinedCalc, eqs, newEqs}, + pdDefs = OptionValue[PartialDerivatives]; + derivs = GridFunctionDerivativesInExpression[pdDefs, lookup[calc, Equations]]; + If[Length[derivs] === 0, Return[calc]]; + newShortDefs = Map[jacobianShorthand, derivs]; + newShorts = Map[First, newShortDefs]; + combinedShorts = Join[lookupDefault[calc, Shorthands, {}], newShorts]; + eqs = lookup[calc, Equations]; + newEqs = eqs /. (x_?(MemberQ[derivs, #] &) :> derivToJacDeriv[x]); + combinedEqs = Join[newShortDefs, newEqs]; + combinedCalc = mapReplace[mapReplace[mapEnsureKey[calc, Shorthands, {}], Shorthands, combinedShorts], Equations, combinedEqs]; + combinedCalc]; + +(* Define local pointers to the members of the Jacobian and Jacobian + derivatives groups *) +CreateJacobianVariables[] := +CommentedBlock["Jacobian variable pointers", + {"bool const use_jacobian = (!CCTK_IsFunctionAliased(\"MultiPatch_GetMap\") || MultiPatch_GetMap(cctkGH) != jacobian_identity_map)\n && strlen(jacobian_group) > 0;\n", + "if (use_jacobian && strlen(jacobian_derivative_group) == 0)\n", + "{\n", + " CCTK_WARN (1, \"GenericFD::jacobian_group and GenericFD::jacobian_derivative_group must both be set to valid group names\");\n", + "}\n\n", + "CCTK_REAL const *restrict jacobian_ptrs[9];\n", + "if (use_jacobian) GenericFD_GroupDataPointers(cctkGH, jacobian_group,\n", + " 9, jacobian_ptrs);\n", + "\n", + Table[{"CCTK_REAL const *restrict const J",i,j," = use_jacobian ? jacobian_ptrs[",(i-1)*3+j-1,"] : 0;\n"},{i,1,3},{j,1,3}], + "\n", + "CCTK_REAL const *restrict jacobian_derivative_ptrs[18];\n", + "if (use_jacobian) GenericFD_GroupDataPointers(cctkGH, jacobian_derivative_group,\n", + " 18, jacobian_derivative_ptrs);\n", + "\n", + Module[{syms = Flatten[Table[{"dJ",i,j,k},{i,1,3},{j,1,3},{k,j,3}],2]}, + MapIndexed[{"CCTK_REAL const *restrict const ", #1, " = use_jacobian ? jacobian_derivative_ptrs[", #2-1, "] : 0;\n"} &, syms]]}]; + +(* List of symbols which should be allowed in a calculation *) +JacobianSymbols[] := + Map[Symbol, Join[Flatten[Table[FlattenBlock[{"dJ",i,j,k}],{i,1,3},{j,1,3},{k,j,3}],2], + Flatten[Table[FlattenBlock[{"J",i,j}],{i,1,3},{j,1,3}],1]]]; + +(* Parameters to inherit from GenericFD *) +JacobianGenericFDParameters[] := + {{Name -> "jacobian_group", Type -> "CCTK_STRING"}, + {Name -> "jacobian_derivative_group", Type -> "CCTK_STRING"}, + {Name -> "jacobian_identity_map", Type -> "CCTK_INT"}}; + +(* The symbols which are used for the Jacobian variables in the + generated source code. These do not have to coincide with the + actual variable names, as the variable pointers are read using + CCTK_VarDataPtr. *) +JacobianGroups[] := + {{"unknown::unknown", {Global`J11, Global`J12, Global`J13, Global`J21, Global`J22, Global`J23, Global`J31, Global`J32, Global`J33}}, + {"unknown::unknown", {Global`dJ111, Global`dJ112, Global`dJ113, Global`dJ122, Global`dJ123, Global`dJ133, + Global`dJ211, Global`dJ212, Global`dJ213, Global`dJ222, Global`dJ223, Global`dJ233, + Global`dJ311, Global`dJ312, Global`dJ313, Global`dJ322, Global`dJ323, Global`dJ333}}}; + +JacobianCheckGroups[groups_] := + Module[{int}, + int = Intersection[allGroupVariables[groups], allGroupVariables[JacobianGroups[]]]; + If[int =!= {}, + Throw["Error: Some group variables conflict with reserved Jacobian variable names: " <> ToString[int]]]]; + +(* These gridfunctions are only given local variable copies if the use_jacobian variable is true *) +JacobianConditionalGridFunctions[] := + {("dJ" ~~ DigitCharacter ~~ DigitCharacter ~~ DigitCharacter) | ("J" ~~ DigitCharacter ~~ DigitCharacter), + "use_jacobian", + None}; + +End[]; + +EndPackage[]; diff --git a/Tools/CodeGen/Kranc.m b/Tools/CodeGen/Kranc.m index 2224d18..ffa98e7 100644 --- a/Tools/CodeGen/Kranc.m +++ b/Tools/CodeGen/Kranc.m @@ -22,7 +22,11 @@ BeginPackage["Kranc`"]; (* CodeGen.m *) -{INV, SQR, CUB, QAD, exp, pow, fmax, fmin, dx, dy, dz, khalf, kthird, ktwothird, kfourthird, keightthird}; +{INV, SQR, CUB, QAD, IfThen, ToReal, sqrt, exp, pow, fmax, fmin, + kmadd, kmsub, knmadd, knmsub, kpos, kneg, kadd, ksub, kmul, kdiv, + kfabs, kfmax, kfmin, ksqrt, kexp, klog, kpow, + dir1, dir2, dir3, dt, dx, dy, dz, + khalf, kthird, ktwothird, kfourthird, keightthird}; (* Helpers.m *) @@ -34,7 +38,7 @@ dummy; LoopPreIncludes, GroupImplementations, PartialDerivatives, NoSimplify, Boundary, Interior, InteriorNoSync, Where, AddToStencilWidth, Everywhere, normal1, normal2, normal3, INV, SQR, CUB, QAD, dot, pow, -exp, dx, dy, dz, idx, idy, idz, MinMod, VanLeer} +exp, dt, dx, dy, dz, idx, idy, idz, MinMod, VanLeer} {ConditionalOnKeyword, ConditionalOnKeywords, CollectList, Interior, InteriorNoSync, Boundary, BoundaryWithGhosts, Where, PreDefinitions, @@ -67,9 +71,11 @@ ThornOptions = ReflectionSymmetries -> {}, ZeroDimensions -> {}, UseLoopControl -> False, - UseCSE -> False, + UseVectors -> False, ProhibitAssignmentToGridFunctionsRead -> False, - IncludeFiles -> {}}; + IncludeFiles -> {}, + CSE -> False, + UseJacobian -> False}; (* Thorn.m *) diff --git a/Tools/CodeGen/KrancGroups.m b/Tools/CodeGen/KrancGroups.m index 05ebdea..13d53d2 100644 --- a/Tools/CodeGen/KrancGroups.m +++ b/Tools/CodeGen/KrancGroups.m @@ -52,6 +52,7 @@ AddGroupExtra; GroupTimelevels; allGroupVariables; NonevolvedTimelevels; +CheckGroups; Begin["`Private`"]; @@ -245,6 +246,18 @@ qualifyGFName[gfname_, allgroups_, defaultImp_] := allGroupVariables[groups_] := Flatten[Map[groupVariables, groups], 1]; +CheckGroups[groups_] := + Module[{vs, names}, +(* If[!MatchQ[{_String, {_Symbol ...}, {_Symbol -> _} ...}], + ThrowError["Groups structure should be of the form {name, {vars, ...}, extras}"]]; *) + + vs = Map[ToString, Union[Flatten[Map[groupVariables, groups]]]]; + names = Map[groupName, groups]; + + If[(int = Intersection[vs,names]) =!= {}, + ThrowError["Variable names and group names must be distinct. Group names which are also variable names:", int]]; + ]; + End[]; EndPackage[]; diff --git a/Tools/CodeGen/KrancTensor.m b/Tools/CodeGen/KrancTensor.m new file mode 100644 index 0000000..204bc60 --- /dev/null +++ b/Tools/CodeGen/KrancTensor.m @@ -0,0 +1,90 @@ +(* ::Package:: *) + +(* Copyright 2004-2010 + Sascha Husa, Ian Hinder, Christiane Lechner, Barry Wardell + + This file is part of Kranc. + + Kranc is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + Kranc is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Kranc; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*) + +(****************************************************************************) +(* Wrapper providing tensor support to Kranc (from TensorTools or xTensor) *) +(****************************************************************************) +If[!ValueQ[$KrancTensorPackage], $KrancTensorPackage = "TensorToolsKranc`"]; + +BeginPackage["KrancTensor`", {"Errors`", "KrancThorn`", "MapLookup`", "KrancGroups`", "Kranc`", $KrancTensorPackage}]; + +CreateKrancThornTT::usage = "Construct a Kranc thorn using tensor expressions."; + +(* FIXME: Move CreateGroupFromTensor here *) + +Begin["`Private`"]; + +(* -------------------------------------------------------------------------- + Tensor Tools + -------------------------------------------------------------------------- *) + +CreateKrancThornTT[groups_, parentDirectory_, thornName_, opts___] := + Module[{calcs, expCalcs, expGroups, options, derivs, expDerivs, reflectionSymmetries, declaredGroups}, + InfoMessage[Terse, "Processing tensorial arguments"]; + calcs = lookup[{opts}, Calculations]; + derivs = lookupDefault[{opts}, PartialDerivatives, {}]; + Map[CheckCalculationTensors, calcs]; + expCalcs = Map[makeCalculationExplicit, calcs]; + + InfoMessage[Info, "Group definitions:", groups]; + VerifyGroups[groups]; + + expDerivs = Flatten[Map[ExpandComponents,derivs],1]; + expGroups = Map[makeGroupExplicit, groups]; + options = Join[DeleteCases[{opts}, Calculations -> _], {Calculations -> expCalcs}]; + options = Join[DeleteCases[options, PartialDerivatives -> _], {PartialDerivatives -> expDerivs}]; + + declaredGroups = lookupDefault[{opts}, DeclaredGroups, {}]; + evolutionTimelevels = lookupDefault[{opts}, EvolutionTimelevels, 3]; + defaultEvolutionTimelevels = lookupDefault[{opts}, DefaultEvolutionTimelevels, evolutionTimelevels]; + InfoMessage[Info, "Declared groups: " <> ToString[declaredGroups]]; + InfoMessage[Terse, "Computing reflection symmetries"]; + reflectionSymmetries = computeReflectionSymmetries[declaredGroups, groups]; + InfoMessage[Info, "Reflection symmetries: ", reflectionSymmetries]; + + InfoMessage[Terse, "Creating (component-based) Kranc thorn"]; + + (* It is necessary to include the KrancThorn context here due to some annoying Needs[] dependency issue *) + KrancThorn`CreateKrancThorn[expGroups, parentDirectory, thornName, + Apply[Sequence, options], ReflectionSymmetries -> reflectionSymmetries]]; + +computeReflectionSymmetries[declaredGroups_, groups_] := + Module[{variables, syms}, + variables = variablesFromGroups[declaredGroups, groups]; + syms = Flatten[Map[ReflectionSymmetries, variables], 1]; + syms]; + +makeCalculationExplicit[calc_] := + mapValueMapMultiple[calc, + {Shorthands -> ExpandComponents, + CollectList -> ExpandComponents, + Equations -> ExpandComponents}]; + +makeGroupExplicit[g_] := + Module[{variables, newVariables, newGroup}, + variables = groupVariables[g]; + newVariables = DeleteDuplicates[ExpandComponents[variables]]; + newGroup = SetGroupVariables[g, newVariables]; + newGroup]; + +End[]; +EndPackage[]; diff --git a/Tools/CodeGen/KrancThorn.m b/Tools/CodeGen/KrancThorn.m index d41fb51..4d53ccb 100644 --- a/Tools/CodeGen/KrancThorn.m +++ b/Tools/CodeGen/KrancThorn.m @@ -1,3 +1,4 @@ +(* ::Package:: *) (* $Id$ *) @@ -27,12 +28,10 @@ BeginPackage["KrancThorn`", {"CodeGen`", "Thorn`", "MapLookup`", "KrancGroups`", "Differencing`", "CalculationFunction`", "Errors`", "Helpers`", "CactusBoundary`", - "TensorTools`", "Param`", "Schedule`", "Interface`", "Kranc`", + "KrancTensor`", "Param`", "Schedule`", "Interface`", "Kranc`", "Jacobian`", "ConservationCalculation`"}]; CreateKrancThorn::usage = "Construct a Kranc thorn"; -CreateKrancThornTT::usage = "Construct a Kranc thorn using TensorTools"; -CreateGroupFromTensor::usage = ""; Begin["`Private`"]; @@ -94,7 +93,7 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[ interface, evolvedGroupDefinitions, rhsGroupDefinitions, thornspec, allParams, boundarySources, reflectionSymmetries, realParamDefs, intParamDefs, - pDefs, useCSE, consCalcs, consCalcsIn, consGroups}, + pDefs, consCalcs, consCalcsIn, consGroups}, (* Parse named arguments *) @@ -125,9 +124,11 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[ partialDerivs = OptionValue[PartialDerivatives] ~Join~ ConservationDifferencingOperators[]; reflectionSymmetries = OptionValue[ReflectionSymmetries]; - useCSE = OptionValue[UseCSE]; coordGroup = {"grid::coordinates", {Kranc`x,Kranc`y,Kranc`z,Kranc`r}}; + + CheckGroups[groupsOrig]; + groups = Join[groupsOrig, {coordGroup}]; includeFiles = Join[includeFiles, {"GenericFD.h", "Symmetry.h", "sbp_calc_coeffs.h"}]; @@ -143,6 +144,8 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[ VerifyString[implementation]; VerifyGroupNames[declaredGroups]; + If[OptionValue[UseJacobian], JacobianCheckGroups[groups]]; + InfoMessage[Terse, "Creating startup file"]; startup = CreateStartupFile[thornName, thornName]; @@ -196,7 +199,7 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[ inheritedRealParams, inheritedIntParams, inheritedKeywordParams, extendedRealParams, extendedIntParams, extendedKeywordParams, evolutionTimelevels, defaultEvolutionTimelevels, - calcs]; + calcs, opts]; (* Construct the schedule file *) InfoMessage[Terse, "Creating schedule file"]; @@ -221,7 +224,10 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[ (* Write the differencing header file *) InfoMessage[Terse, "Creating differencing header file"]; - {pDefs, diffHeader} = CreateDifferencingHeader[partialDerivs, OptionValue[ZeroDimensions]]; + {pDefs, diffHeader} = CreateDifferencingHeader[partialDerivs, OptionValue[ZeroDimensions], OptionValue[UseVectors]]; + diffHeader = Join[ + If[OptionValue[UseVectors], {"#include \"vectors.h\"\n", "\n"}, {}], + diffHeader]; (* Add the predefinitions into the calcs *) calcs = Map[Join[#, {PreDefinitions -> pDefs}] &, calcs]; @@ -238,12 +244,12 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[ InfoMessage[Terse, "Creating calculation source files"]; calcSources = Map[CreateSetterSource[ {Join[#, {Parameters -> allParams, PartialDerivatives -> partialDerivs}]}, - False, useCSE, {}, implementation, opts] &, calcs]; + False, {}, implementation, opts] &, calcs]; calcFilenames = Map[lookup[#, Name] <> ext &, calcs]; (* Makefile *) InfoMessage[Terse, "Creating make file"]; - make = CreateMakefile[Join[{"Startup.c", "RegisterMoL.c", "RegisterSymmetries.c"}, calcFilenames, + make = CreateMakefile[Join[{"Startup.cc", "RegisterMoL.cc", "RegisterSymmetries.cc"}, calcFilenames, Map[lookup[#, Filename] &, boundarySources]]]; (* Put all the above together and generate the Cactus thorn *) @@ -255,10 +261,10 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[ Param -> param, Makefile -> make, Sources -> Join[{ - {Filename -> "Startup.c", Contents -> startup}, - {Filename -> "RegisterMoL.c", Contents -> molregister}, - {Filename -> "RegisterSymmetries.c", Contents -> symregister}, - {Filename -> "Differencing.h", Contents -> diffHeader}}, + {Filename -> "Startup.cc", Contents -> startup}, + {Filename -> "RegisterMoL.cc", Contents -> molregister}, + {Filename -> "RegisterSymmetries.cc", Contents -> symregister}, + {Filename -> "Differencing.h", Contents -> diffHeader}}, MapThread[{Filename -> #1, Contents -> #2} &, {calcFilenames, calcSources}], boundarySources]}; InfoMessage[Terse, "Creating thorn"]; @@ -327,6 +333,7 @@ createKrancMoLRegister[evolvedGroupNames_, nonevolvedGroupNames_, groups_, imple molregister = CreateMoLRegistrationSource[molspec, False]; Return[molregister]]; +<<<<<<< HEAD (* -------------------------------------------------------------------------- Tensors -------------------------------------------------------------------------- *) @@ -457,5 +464,7 @@ CheckCalculationTensors[calc_] := eqs = lookup[calc, Equations]; Map[CheckEquationTensors, eqs]]; +======= +>>>>>>> origin/master End[]; EndPackage[]; diff --git a/Tools/CodeGen/Optimize.m b/Tools/CodeGen/Optimize.m new file mode 100644 index 0000000..2f791da --- /dev/null +++ b/Tools/CodeGen/Optimize.m @@ -0,0 +1,190 @@ +(* ::Package:: *) + +(* Copyright 2011 Barry Wardell + + This file is part of Kranc. + + Kranc is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + Kranc is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Kranc; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*) + +BeginPackage["Optimize`", {"Kranc`", "Errors`"}]; + +EliminateCommonSubexpressions::usage = "EliminateCommonSubexpressions[calc] identifies common subexpressions in calc and introduces new shorthands for them."; +TopologicallySortEquations::usage = "TopologicallySortEquations[eqs, v] sorts eqs topologically, including Head v as a possible vertex." + +Begin["`Private`"]; + +CSEPrint[___] = null; +(* CSEPrint = Print; *) + +Options[EliminateCommonSubexpressions] = ThornOptions; +EliminateCommonSubexpressions[calc_List, OptionsPattern[]] := + Module[{eqs, shorts, name, pdDefs, derivs, newShorts, newEqs, allShorts, newCalc}, + name = (Name /. calc); + + InfoMessage[InfoFull, "Doing common subexpression elimination for "<>name]; + + eqs = (Equations /. calc) /. Equations -> {}; + shorts = (Shorthands /. calc) /. Shorthands -> {}; + + (* Get a list of symbols used for derivatives. We will not eliminate these as subexpressions. *) + pdDefs = OptionValue[PartialDerivatives]; + derivs = DeleteDuplicates[Head/@(First/@pdDefs)]; + + (* Generate new equations with subexpressions eliminated. *) + {newShorts, newEqs} = cse[eqs, Symbol["csetemp"], derivs]; + + If[Length[newShorts]>0, + InfoMessage[Info, "Extracted "<>ToString[Length[newShorts]]<>" common subexpressions from "<>name]; + ]; + + allShorts = Join[shorts, newShorts]; + newCalc = Join[calc /. {(Shorthands->_) -> Sequence[], (Equations->_) -> Sequence[]}, + {Shorthands->allShorts}, {Equations->newEqs}]; + + newCalc +]; + +cse[eqs_, v_, exceptions_, minSaving_:0] := + Module[{subexprs, replacements, replace, newEqs, defs, newDefs, i, relabelVars, allEqs, sortedEqs, newVars}, + (* Find all possible subexpressions and how many times they occur *) + CSEPrint["CSE"]; + CSEPrint["CSE: eqs=", eqs]; + subexprs = Reap[Scan[If[! AtomQ[#], Sow[#]] &, eqs[[All,2]], Infinity]]; + CSEPrint["CSE: subexprs=", subexprs]; + If[subexprs[[2]]=={}, Return[{{}, eqs}]]; + subexprs = Tally[subexprs[[2, 1]]]; + + (* Discard subexpressions which only appear once *) + subexprs = Select[subexprs, #[[2]] >= 2 &]; + + (* Sort subexpressions in ascending order by size (LeafCount) *) + subexprs = Sort[subexprs, LeafCount[#1] < LeafCount[#2] &]; + + (* Ony keep subexpressions larger than minSaving=(numoccurances-1)*size *) + subexprs = Select[subexprs, (#[[2]]-1) LeafCount[#[[1]]] >= minSaving &][[All,1]]; + + (* Discard some specific cases *) + subexprs = Cases[subexprs, Except[_?AtomQ]]; + subexprs = Cases[subexprs, Except[Times[-1, _?AtomQ]]]; (* -x *) + subexprs = Cases[subexprs, Except[Alternatives@@(Blank/@exceptions)]]; (* specified exceptions *) + subexprs = Cases[subexprs, Except[Times[-1, Alternatives@@(Blank/@exceptions)]]]; (* -exceptions *) + + (* Get the list of replacements for our original expression *) + replacements = Thread[subexprs -> Table[v[i], {i, Length[subexprs]}]]; + + (* Replace common subexpressions with new variables *) + (* Do not replace certain terms, e.g. the first argument of IfThen. *) + (* newEqs = eqs //. replacements; *) + CSEPrint["CSE: eqs=", eqs]; + CSEPrint["CSE: replacements=", replacements]; + replace[expr_] := Replace[Switch[expr, + IfThen[_,_,_], IfThen[expr[[1]], replace[expr[[2]]], replace[expr[[3]]]], + (* ToReal[_], ToReal[expr[[1]]], *) + _?AtomQ, expr, + _, Map[replace, expr]], + replacements]; + newEqs = FixedPoint[replace, eqs]; + CSEPrint["CSE: newEqs=", newEqs]; + + (* Build up definitions for the new variables *) + defs = Reverse/@replacements; + CSEPrint["CSE: defs=", defs]; + For[i = 2, i <= Length[subexprs], i++, + defs[[i,2]] = defs[[i,2]] /. replacements[[1;;i-1]]; + ]; + CSEPrint["CSE: defs=", defs]; + + (* Select only the definitions which are needed for the new expressions. + This accounts for cases where a subexpression appears multiple times, + but always as part of the same larger subexpression. For example, in + expr = Sqrt[(a+b)(a-b)c]+(a+b)(a-b)c+(a+b)d+Sqrt[(a+b)d+(a+b)c]; + we would identify the subexpressions + {v[1]->a+b,v[2]->d v[1],v[3]->a-b,v[4]->c v[1] v[3]}; + whereas all we really want it to identify is + {v[1]->a+b,v[2]->d v[1],v[4]->(a-b) c v[1]}; + and the introduction of v[3] is unnecessary. To achieve this, we only + keep temporary variables which appear in the expression after substitution + or which appear more than once in the definition of the temporary variables. + *) + newDefs = Select[defs, (Count[newEqs, #[[1]], Infinity] > 0) || + (Count[defs[[All,2]], #[[1]], Infinity] > 1) &]; + CSEPrint["CSE: newDefs=", newDefs]; + + (* Replace any temporaries eliminated by the previous procedure with their definition *) + newDefs = newDefs //. Complement[defs, newDefs]; + CSEPrint["CSE: newDefs2=", newDefs]; + + (* Check we actually have subexpressions to eliminate. Otherwise just return the original expression *) + If[Length[newDefs]==0, Return[{{}, eqs}]]; + + (* This is our new system of equations *) + allEqs = Join[newDefs, newEqs]; + CSEPrint["CSE: allEqs=", allEqs]; + + sortedEqs = Fold[InsertNewEquation, newEqs, Reverse[newDefs]]; + CSEPrint["CSE: sortedEqs=", sortedEqs]; + + (* Relabel new temporary variables so that they are sequential and C friendly *) + newVars = Select[sortedEqs[[All,1]], MemberQ[newDefs[[All, 1]],#]&]; + CSEPrint["CSE: newVars=", newVars]; + i = 0; + relabelVars = (# -> Symbol[ToString[v] <> ToString[i++]]) & /@ newVars; + + (* Return the list of new variables and the new equations *) + {newVars, sortedEqs} /. relabelVars +]; + +TopologicallySortEquations[eqs_] := Module[{lhs, rhs, lhsInrhs, dag, sortedVars, indVars, allVars, sortedEqs}, + lhs = eqs[[All,1]]; + rhs = eqs[[All,2]]; + + (* Generate an directed acyclic graph for the system of equations *) + lhsInrhs = DeleteDuplicates[Cases[{#}, _?(MemberQ[lhs, #] &), Infinity]] & /@ rhs; + dag = Graph[ Flatten[MapThread[Thread[Rule[#1, #2]] &, {lhsInrhs, lhs}]] ]; + + (* Topologically sort the DAG *) + sortedVars = Quiet[TopologicalSort[dag], TopologicalSort::argx]; + + (* Check if the topological sorting failed. This can happen if the graph for + the equations is cyclic. For example, we could have a->a+1 or {b->a*a, a->b} *) + If[SameQ[Head[sortedVars], TopologicalSort], + InfoMessage[Info, "Failed to topologically sort equations."]; + Return[$Failed] + ]; + + (* Some variables might be independent. Add them back in. *) + indVars = Complement[lhs, sortedVars]; + allVars = Join[indVars, sortedVars]; + + sortedEqs = Thread[allVars -> (allVars/.eqs)]; + sortedEqs +]; + +InsertNewEquation[oldEqs_, newEq_] := Module[{before}, + CSEPrint["InsertNewEquation oldEqs=", oldEqs, " newEq=", newEq]; + (* For some reason, we can be asked to insert an equation that is + not actually needed. This should not be the case. However, handle + it gracefully for now. *) + (* before = Position[oldEqs[[All,2]], newEq[[1]]][[1,1]]; *) + before = Position[oldEqs[[All,2]], newEq[[1]]]; + If[before=={}, + oldEqs, + Insert[oldEqs, newEq, before[[1,1]]]] +]; + +End[]; + +EndPackage[]; diff --git a/Tools/CodeGen/Param.m b/Tools/CodeGen/Param.m index 7821f40..1b4ad37 100644 --- a/Tools/CodeGen/Param.m +++ b/Tools/CodeGen/Param.m @@ -18,7 +18,7 @@ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA *) -BeginPackage["Param`", {"Thorn`", "Errors`", "Helpers`", "MapLookup`", "KrancGroups`", "Kranc`"}]; +BeginPackage["Param`", {"Thorn`", "Errors`", "Helpers`", "MapLookup`", "KrancGroups`", "Kranc`", "Jacobian`"}]; CreateKrancParam; MakeFullParamDefs; @@ -72,14 +72,16 @@ krancParamStructExtended[definition_, type_] := AllowedValues -> Map[{Value -> #, Description -> ""} &, allowedValues]}]; krancKeywordParamStruct[struct_] := -{ - Name -> lookup[struct, Name], - Type -> "KEYWORD", - Default -> lookup[struct, Default], - Description -> lookupDefault[struct, Description, lookup[struct, Name]], - Visibility -> lookupDefault[struct, Visibility, "private"], - AllowedValues -> Map[{Value -> #, Description -> #} &, lookup[struct, AllowedValues]] -}; + Join[ + {Name -> lookup[struct, Name], + Type -> "KEYWORD", + Default -> lookup[struct, Default], + Description -> lookupDefault[struct, Description, lookup[struct, Name]], + Visibility -> lookupDefault[struct, Visibility, "private"]}, + If[mapContains[struct, Steerable], + {Steerable -> lookup[struct, Steerable]}, + {}], + {AllowedValues -> Map[{Value -> #, Description -> #} &, lookup[struct, AllowedValues]]}]; MakeFullParamDefs[params_] := Module[{p}, @@ -116,12 +118,13 @@ extendParameters[imp_, reals_, ints_, keywords_] := Return[{Name -> imp, ExtendedParameters -> Join[realStructs, intStructs, keywordStructs]}], Return[{}]]]; +Options[CreateKrancParam] = ThornOptions; CreateKrancParam[evolvedGroups_, nonevolvedGroups_, groups_, thornName_, reals_, ints_, keywords_, inheritedReals_, inheritedInts_, inheritedKeywords_, extendedReals_, extendedInts_, extendedKeywords_, evolutionTimelevels_, defaultEvolutionTimelevels_, - calcs_] := + calcs_, opts:OptionsPattern[]] := Module[{nEvolved, evolvedMoLParam, evolvedGFs, (*constrainedMoLParam,*) genericfdStruct, realStructs, intStructs, allInherited, allExtended, implementationNames, molImplementation, @@ -146,7 +149,8 @@ CreateKrancParam[evolvedGroups_, nonevolvedGroups_, groups_, thornName_, Visibility -> "restricted", AccumulatorBase -> "MethodofLines::MoL_Num_Evolved_Vars", AllowedValues -> {{Value -> ToString[nEvolved] <> ":" <> ToString[nEvolved] , - Description -> "Number of evolved variables used by this thorn"}} + Description -> "Number of evolved variables used by this thorn"}}, + Steerable -> Recover }; (* @@ -171,7 +175,8 @@ CreateKrancParam[evolvedGroups_, nonevolvedGroups_, groups_, thornName_, Description -> "Number of active timelevels", Visibility -> "restricted", AllowedValues -> {{Value -> ToString[0] <> ":" <> ToString[evolutionTimelevels], - Description -> ""}} + Description -> ""}}, + Steerable -> Recover }; rhsTimelevelsParam = @@ -182,25 +187,22 @@ CreateKrancParam[evolvedGroups_, nonevolvedGroups_, groups_, thornName_, Description -> "Number of active RHS timelevels", Visibility -> "restricted", AllowedValues -> {{Value -> ToString[0] <> ":" <> ToString[evolutionTimelevels], - Description -> ""}} + Description -> ""}}, + Steerable -> Recover }; genericfdStruct = { Name -> "GenericFD", UsedParameters -> - {{Name -> "stencil_width", Type -> "CCTK_INT"}, - {Name -> "stencil_width_x", Type -> "CCTK_INT"}, - {Name -> "stencil_width_y", Type -> "CCTK_INT"}, - {Name -> "stencil_width_z", Type -> "CCTK_INT"}, - {Name -> "boundary_width", Type -> "CCTK_INT"}} + If[OptionValue[UseJacobian], JacobianGenericFDParameters[], {}] }; realStructs = Map[krancParamStruct[#, "CCTK_REAL", False] &, reals]; - verboseStruct = krancParamStruct[{Name -> "verbose", Default -> 0}, "CCTK_INT", False]; + verboseStruct = krancParamStruct[{Name -> "verbose", Default -> 0, Steerable -> Always}, "CCTK_INT", False]; intStructs = Map[krancParamStruct[#, "CCTK_INT", False] &, ints]; - calcEveryStructs = Map[krancParamStruct[{Name -> lookup[#, Name] <> "_calc_every", Default -> 1}, "CCTK_INT", False] &, calcs]; - calcOffsetStructs = Map[krancParamStruct[{Name -> lookup[#, Name] <> "_calc_offset", Default -> 0}, "CCTK_INT", False] &, calcs]; + calcEveryStructs = Map[krancParamStruct[{Name -> lookup[#, Name] <> "_calc_every", Default -> 1, Steerable -> Always}, "CCTK_INT", False] &, calcs]; + calcOffsetStructs = Map[krancParamStruct[{Name -> lookup[#, Name] <> "_calc_offset", Default -> 0, Steerable -> Always}, "CCTK_INT", False] &, calcs]; keywordStructs = Map[krancKeywordParamStruct, keywords]; allInherited = Join[inheritedReals, inheritedInts, inheritedKeywords]; diff --git a/Tools/CodeGen/TensorTools.m b/Tools/CodeGen/TensorTools.m index c6f4cfd..6579c63 100644 --- a/Tools/CodeGen/TensorTools.m +++ b/Tools/CodeGen/TensorTools.m @@ -217,10 +217,10 @@ IndexIsLower[TensorIndex[_, "u"]] := False; -------------------------------------------------------------------------- *) Format[TensorIndex[label_, "u"], OutputForm] := - Superscript[null,label]; + "u"<>ToString[label]; Format[TensorIndex[label_, "l"], OutputForm] := - Subscript[null,label]; + "l"<>ToString[label]; Format[TensorIndex[label_, "u"], StandardForm] := Superscript[null,label]; @@ -251,12 +251,12 @@ DefineTensor[T_] := Format[Tensor[T, is:((TensorIndex[_,_] | _Integer) ..) ], StandardForm] := PrecedenceForm[ - SequenceForm[T,is], + SequenceForm[T,"[",Sequence@@Riffle[{is},","],"]"], 10000]; Format[Tensor[T, is:((TensorIndex[_,_] | _Integer) ..) ], OutputForm] := PrecedenceForm[ - SequenceForm[T,is], + SequenceForm[T,"[",Sequence@@Riffle[{is},","],"]"], 10000]; (* Cannot get InputForm to work *) @@ -265,7 +265,7 @@ DefineTensor[T_] := HoldForm[T[is]];*) T[is:((TensorIndex[_,_] | _Integer) ..)] := Tensor[T, is]; - TensorAttributes[T] = {TensorWeight -> 1, Symmetries -> {}}; + TensorAttributes[T] = {TensorWeight -> 0, Symmetries -> {}}; T]; (* -------------------------------------------------------------------------- @@ -1235,10 +1235,11 @@ CheckTensors[x_, y_] := ]; CheckTensors[t:Tensor[k_, is__]] := - Module[{}, + Module[{is2}, (* Print["CheckTensors: Tensor: ", t];*) - If[!(Union[{is}] === Sort[{is}]), - ThrowError["Tensor has repeated indices: ", t, {is}]]; + is2 = Select[{is}, !NumericQ[#]&]; + If[!(Union[is2] === Sort[is2]), + ThrowError["Tensor has repeated indices: ", t, is2]]; True]; CheckTensors[t:f_[TensorIndex[__]..]] := diff --git a/Tools/CodeGen/TensorToolsKranc.m b/Tools/CodeGen/TensorToolsKranc.m new file mode 100644 index 0000000..1abda8b --- /dev/null +++ b/Tools/CodeGen/TensorToolsKranc.m @@ -0,0 +1,106 @@ +(* ::Package:: *) + +(* Copyright 2004-2010 + Sascha Husa, Ian Hinder, Christiane Lechner, Barry Wardell + + This file is part of Kranc. + + Kranc is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + Kranc is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Kranc; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*) + +(* This package provides a TensorTools wrapper for Kranc *) + +BeginPackage["TensorToolsKranc`", {"Kranc`", "KrancGroups`", "TensorTools`"}]; + +CreateGroupFromTensor::usage = ""; +ReflectionSymmetries::usage = "Produce a list of reflection symmetries of a tensor."; +ExpandComponents::usage = "ExpandComponents[expr] converts an expression containing abstract indices into one containing components instead." +IncludeCharacter::usage = "IncludeCharacter is an option for makeExplicit which specifies whether the character should also be included in the generated variable names." +TensorCharacterString::usage = "TensorCharacterString[tensor[inds]] returns a string consisting of a sequence of U's and D's representing the character of tensor." + +Begin["`Private`"]; + +ReflectionSymmetries[x___] := ReflectionSymmetriesOfTensor[x]; +ExpandComponents[expr_] := MakeExplicit[expr]; + +TensorCharacterString[k_, inds_] := + Module[{}, + InfoMessage[InfoFull, "Tensor attributes of " <> ToString[k], TensorAttributes[k]]; + If[HasTensorAttribute[k, TensorManualCartesianParities], + "ManualCartesian", + If[Length[inds] == 0, + "Scalar", + Apply[StringJoin, Map[If[IndexIsUpper[#], "U", "D"] &, inds]]]]]; + +CreateGroupFromTensor[T:Tensor[k_, is__]] := + CreateGroupFromTensor[k, {is}]; + +reflectionParityString[l_] := + Module[{chars}, + If[!ListQ[l] || !Length[l] == 3, + ThrowError["Expecting a list of three parities for TensorManualCartesianParities, must be 1 or -1"]]; + + chars= Map[Switch[#, -1, "-", +1, "+", _, + ThrowError["Expecting a list of three parities for TensorManualCartesianParities, must be 1 or -1"]] &, + l]; + + Apply[StringJoin, chars]]; + +CreateGroupFromTensor[k_, inds_] := + Module[{ttypeString, nInds, tags, group, vars}, + InfoMessage[InfoFull, "Creating group from tensor with kernel " <> ToString[k] <> " and indices " <> ToString[inds]]; + ttypeString = TensorCharacterString[k, inds]; + InfoMessage[InfoFull, "Tensor type string: ", ttypeString]; + nInds = Length[inds]; + If[nInds == 2 && GetTensorAttribute[k, Symmetries] == {{2,1},1}, + ttypeString = ttypeString <> "_sym"]; + If[nInds == 3 && GetTensorAttribute[k, Symmetries] == {{1,3,2},1}, + ttypeString = ttypeString <> "_sym"]; + tags = {"tensortypealias" -> ttypeString, "tensorweight" -> GetTensorAttribute[k, TensorWeight]}; + If[HasTensorAttribute[k, TensorSpecial], + tags = Append[tags, "tensorspecial" -> GetTensorAttribute[k, TensorSpecial]]]; + If[HasTensorAttribute[k, TensorManualCartesianParities], + tags = Append[tags, "cartesianreflectionparities" -> + reflectionParityString[GetTensorAttribute[k, TensorManualCartesianParities]]]]; + If[HasTensorAttribute[k, TensorParity], + tags = Append[tags, "tensorparity" -> GetTensorAttribute[k, TensorParity]]]; + + If[HasTensorAttribute[k, Checkpoint], + tags = Append[tags, "checkpoint" -> GetTensorAttribute[k, Checkpoint]]]; + + vars = If[nInds == 0, {k}, {Apply[Tensor, {k, Apply[Sequence,inds]}]}]; + group = CreateGroup[ToString[k] <> "_group", vars, {Tags -> tags}]; + Return[group]]; + +CreateGroupFromTensor[x_] := + If[IsTensor[x], + CreateGroupFromTensor[x, {}], + ThrowError["CreateGroupFromTensor: Not a tensor", x]]; + +CheckEquationTensors[eq_] := + Module[{}, + CheckTensors[eq]]; + +CheckCalculationTensors[calc_] := + Module[{eqs}, + + If[mapContains[calc, Shorthands], + CheckTensors[lookup[calc, Shorthands]]]; + + eqs = lookup[calc, Equations]; + Map[CheckEquationTensors, eqs]]; + +End[]; +EndPackage[]; diff --git a/Tools/CodeGen/Thorn.m b/Tools/CodeGen/Thorn.m index 06b3e8d..9381948 100644 --- a/Tools/CodeGen/Thorn.m +++ b/Tools/CodeGen/Thorn.m @@ -220,7 +220,8 @@ Options[CreateConfiguration] = ThornOptions; CreateConfiguration[opts:OptionsPattern[]] := {whoWhen["CCL"], "REQUIRES GenericFD\n", - If[OptionValue[UseLoopControl], "REQUIRES LoopControl\n", {}] + If[OptionValue[UseLoopControl], "REQUIRES LoopControl\n", {}], + If[OptionValue[UseVectors], "REQUIRES Vectors\n", {}] }; (* ------------------------------------------------------------------------ @@ -473,14 +474,19 @@ CreateSchedule[globalStorageGroups_, scheduledGroups_, scheduledFunctions_] := optional LoopPreIncludes -> {include file list}, Equations -> {{K11_rhs -> 2 A K11, ...}...}} *) -calculationMacros[] := +calculationMacros[vectorise_] := CommentedBlock["Define macros used in calculations", - Map[{"#define ", #, "\n"} &, - {"INITVALUE (42)", - "INV(x) ((1.0) / (x))" , - "SQR(x) ((x) * (x))" , - "CUB(x) ((x) * (x) * (x))" , - "QAD(x) ((x) * (x) * (x) * (x))"}]]; + Map[{"#define ", #, "\n"} &, + {"INITVALUE (42)", + "QAD(x) (SQR(SQR(x)))"} ~Join~ + If[vectorise, + {"INV(x) (kdiv(ToReal(1.0),x))", + "SQR(x) (kmul(x,x))", + "CUB(x) (kmul(x,SQR(x)))"}, + {"INV(x) ((1.0) / (x))", + "SQR(x) ((x) * (x))", + "CUB(x) ((x) * (x) * (x))"}] + ]]; (* Given a list of Calculation structures as defined above, create a CodeGen representation of a source file that defines a function for @@ -488,7 +494,7 @@ calculationMacros[] := Options[CreateSetterSource] = ThornOptions; -CreateSetterSource[calcs_, debug_, useCSE_, include_, imp_, +CreateSetterSource[calcs_, debug_, include_, imp_, opts:OptionsPattern[]] := Module[{}, @@ -503,21 +509,24 @@ CreateSetterSource[calcs_, debug_, useCSE_, include_, imp_, {IncludeSystemFile["assert.h"], IncludeSystemFile["math.h"], IncludeSystemFile["stdio.h"], - IncludeSystemFile["stdlib.h"]}, + IncludeSystemFile["stdlib.h"], + IncludeSystemFile["string.h"]}, {"\n"} ], Map[IncludeFile, Join[{"cctk.h", "cctk_Arguments.h", "cctk_Parameters.h", - (*"precomputations.h",*) "GenericFD.h", "Differencing.h"}, include, - If[OptionValue[UseLoopControl], {"loopcontrol.h"}, {}]]], - calculationMacros[], + (*"precomputations.h",*) "GenericFD.h", "Differencing.h"}, + include, + If[OptionValue[UseLoopControl], {"loopcontrol.h"}, {}], + If[OptionValue[UseVectors], {"vectors.h"}, {}]]], + calculationMacros[OptionValue[UseVectors]], (* For each function structure passed, create the function and insert it *) CalculationBoundariesFunction[First[calcs], imp], - Map[CreateCalculationFunction[# , debug, useCSE, opts] &, + Map[CreateCalculationFunction[# , debug, imp, opts] &, calcs]}]; @@ -749,10 +758,10 @@ CreateMoLBoundariesSource[spec_] := "if (CCTK_EQUALS(" <> boundpar <> ", \"none\" ) ||\n", " CCTK_EQUALS(" <> boundpar <> ", \"static\") ||\n", " CCTK_EQUALS(" <> boundpar <> ", \"flat\" ) ||\n", - " CCTK_EQUALS(" <> boundpar <> ", \"zero\" ) ) \n", + " CCTK_EQUALS(" <> boundpar <> ", \"zero\" ) )\n", "{\n", - " ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, \n", + " ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1,\n", " \"" <> fullgroupname <> "\", " <> boundpar <> ");\n", " if (ierr < 0)\n", @@ -771,10 +780,10 @@ CreateMoLBoundariesSource[spec_] := "if (CCTK_EQUALS(" <> boundpar <> ", \"none\" ) ||\n", " CCTK_EQUALS(" <> boundpar <> ", \"static\") ||\n", " CCTK_EQUALS(" <> boundpar <> ", \"flat\" ) ||\n", - " CCTK_EQUALS(" <> boundpar <> ", \"zero\" ) ) \n", + " CCTK_EQUALS(" <> boundpar <> ", \"zero\" ) )\n", "{\n", - " ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, \n", + " ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,\n", " \"" <> fullgfname <> "\", " <> boundpar <> ");\n", " if (ierr < 0)\n", @@ -1178,7 +1187,7 @@ headerComment2, " assert (len);\n\n", " for (d = 0; d < 3; ++d) {\n", -" assert (off[d] >= 0 && len[d] >= 0 && off[d] + len[d] <= cctk_lssh[CCTK_LSSH_IDX(0,d)]);\n", +" assert (off[d] >= 0 && len[d] >= 0 && off[d] + len[d] <= CCTK_LSSH(0,d));\n", " }\n\n", " assert (modes);\n", @@ -1340,7 +1349,7 @@ CreateStartupFile[thornName_, bannerText_] := tmp = {whoWhen["C"], IncludeFile["cctk.h"], - DefineFunction[thornName <> "_Startup", "int", "void", + DefineFunction[thornName <> "_Startup", "extern \"C\" int", "void", {DefineVariable["banner", "const char *", Quote[bannerText]], "CCTK_RegisterBanner(banner);\n", "return 0;\n"}]}; @@ -1354,7 +1363,7 @@ CreateStartupFile[thornName_, bannerText_] := Thorn creation ------------------------------------------------------------------------ *) -(* source = {Filename -> "MoLRegister.c", Contents -> "#include ..."} *) +(* source = {Filename -> "MoLRegister.cc", Contents -> "#include ..."} *) (* thorn = {Name -> "ClassicADMMolEvolve", Directory -> "ClassicADM", Interface -> i, Schedule -> s, Param -> p, Makefile -> m, diff --git a/Tools/CodeGen/xTensorKranc.m b/Tools/CodeGen/xTensorKranc.m new file mode 100644 index 0000000..f92072f --- /dev/null +++ b/Tools/CodeGen/xTensorKranc.m @@ -0,0 +1,136 @@ +(* ::Package:: *) + +(* Copyright 2010 Barry Wardell + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 +*) + +BeginPackage["xTensorKranc`", {"Differencing`", "Kranc`", "KrancGroups`", "xAct`xTensor`", "xAct`xCore`", "xAct`xCoba`"}]; + +CreateGroupFromTensor::usage = ""; +ReflectionSymmetries::usage = "Produce a list of reflection symmetries of a tensor."; +ExpandComponents::usage = "ExpandComponents[expr] converts an expression x containing abstract indices into one containing components instead." +IncludeCharacter::usage = "IncludeCharacter is an option for makeExplicit which specifies whether the character should also be included in the generated variable names." +TensorCharacterString::usage = "TensorCharacterString[tensor[inds]] returns a string consisting of a sequence of U's and D's representing the character of tensor." + +Begin["`Private`"]; + +(* FIXME: Add support for ManualCartesian attribute *) +TensorCharacterString[t_Symbol?xTensorQ[]] := "Scalar"; +TensorCharacterString[t_Symbol?xTensorQ[inds___]] := StringJoin[If[UpIndexQ[#],"U","D"]&/@{inds}]; + +Options[makeExplicit] = {IncludeCharacter -> False}; +SetAttributes[makeExplicit, Listable]; +e : makeExplicit[_Plus, opts:OptionsPattern[]] := Distribute[Unevaluated[e]]; +makeExplicit[x_Times, opts:OptionsPattern[]] := Map[makeExplicit[#, opts]&, x]; +makeExplicit[Power[x_,p_], opts:OptionsPattern[]] := Power[makeExplicit[x, opts],p]; +makeExplicit[x_?NumericQ, opts:OptionsPattern[]] := x; +makeExplicit[x_, {}, opts:OptionsPattern[]] := makeExplicit[x]; +makeExplicit[t_Symbol?xTensorQ[inds___], opts:OptionsPattern[]] := Module[{indexNumbers,character,indexString}, + indexNumbers=First/@{inds}; + + If[OptionValue[IncludeCharacter], + character = TensorCharacterString[t[inds]]; + indexString = StringJoin[character,ToString/@indexNumbers], + indexString = StringJoin[ToString/@indexNumbers] + ]; + SymbolJoin[PrintAs[t],Sequence@@indexString] +]; + +makeExplicit[x_, opts:OptionsPattern[]] := x; + +makeExplicit[(cd_?CovDQ)[ind_][expr_], opts:OptionsPattern[]] := Module[{indexNumbers}, + indexNumbers=First/@{ind}; + + Global`PDstandard2nd[makeExplicit[expr, opts], Sequence@@indexNumbers] +]; + +Options[ExpandComponents] = Options[ExpandComponents]; + +ExpandComponents[x_Rule, opts:OptionsPattern[makeExplicit]] := Thread[ExpandComponents[x[[1]], opts] -> ExpandComponents[x[[2]], opts]]; +ExpandComponents[dot[x_], opts:OptionsPattern[makeExplicit]] := dot/@ExpandComponents[x, opts]; +ExpandComponents[x_List, opts:OptionsPattern[makeExplicit]] := Flatten[Map[ExpandComponents[#, opts]&, x], 1]; +ExpandComponents[x_, opts:OptionsPattern[makeExplicit]] := + Module[{eqs, options}, + + eqs = ComponentArray[TraceBasisDummy[x]]; + options = Evaluate[FilterRules[{opts}, Options[makeExplicit]]]; + If[Length[options]==0, + makeExplicit[eqs], + makeExplicit[eqs, options] + ] +]; + + +(* Compute the reflection symmetries of a tensor *) +ReflectionSymmetries[t_Symbol?xTensorQ[inds__]] := + Module[{b=Global`Euclidean, cnums, components, componentIndices, counts}, + (* Get the compoent indices of the basis *) + cnums = CNumbersOf[b, VBundleOfBasis[b]]; + + (* Get a list of components of the tensor t in the basis b *) + components = Flatten[ComponentArray[ToBasis[b][t[inds]]]]; + + (* Get the indices of each component *) + componentIndices = Map[IndicesOf[b], components]; + + (* Count the number of instances of each basis index. *) + countInds[expr_, basis_, cinds_] := Map[(Count[expr,{#,basis}]+Count[expr,{#,-basis}])&, cinds]; + counts = Map[countInds[#, b, cnums]&, componentIndices]; + + (* For each instance, multiply by -1 *) + Thread[ExpandComponents[t[inds]] -> (-1)^counts] +]; + +ReflectionSymmetries[t_Symbol?xTensorQ[]] := t -> {1,1,1}; +ReflectionSymmetries[t_] := t -> {1, 1, 1}; + +(* FIXME: Implement this fully *) +GetTensorAttribute[t_Symbol?xTensorQ, TensorWeight] := WeightOfTensor[t]; + +CreateGroupFromTensor[t_Symbol?xTensorQ[inds__]] := Module[{tCharString, nInds, tags, vars, group}, + InfoMessage[InfoFull, "Creating group from tensor with kernel " <> SymbolName[t] <> " and indices " <> ToString[{inds}]]; + + (* Get a string representing the character of the tensor *) + tCharString = TensorCharacterString[t[inds]]; + InfoMessage[InfoFull, "Tensor character string: ", tCharString]; + + (* Check if the tensor is symmetric *) + nInds = Length[SlotsOfTensor[t]]; + If[SymmetryGroupOfTensor[t] == StrongGenSet[Range[nInds],GenSet[Cycles[Range[nInds]]]], + tCharString = tCharString <> "_sym"]; + + (* FIXME: Add tensorspecial, cartesianreflectionparities and tensorparity *) + tags = {"tensortypealias" -> tCharString, "tensorweight" -> GetTensorAttribute[t, TensorWeight]}; + + vars = If[nInds == 0, {t}, {t[inds]}]; + group = CreateGroup[SymbolName[t] <> "_group", vars, {Tags -> tags}]; + Return[group] +]; + +ReflectionSymmetries[x___]:= Throw["ReflectionSymmetries error: "<>ToString[x]]; +CreateGroupFromTensor[x___]:= Throw["CreateGroupFromTensor error: "<>ToString[x]]; + +CheckTensors[expr_] := Validate[expr]; + +End[]; +EndPackage[]; + + + + + + + diff --git a/Tools/External/Optimize.m b/Tools/External/Optimize.m deleted file mode 100644 index 9f51225..0000000 --- a/Tools/External/Optimize.m +++ /dev/null @@ -1,741 +0,0 @@ -(* :Title: Expression Optimization. *) - -(* :Author: Mark Sofroniou *) - -(* :Summary: - This package adds the procedure Optimize for performing common - sub-expression optimization on Mathematica expressions. - An optimized Module and an optimized compiled Module is also possible. - The function Horner factors uni-variate and multi-variate polynomials - in efficient computational form. - The function Cost measures the computational expense of numerical - evaluation of an expression. - The additional package Format.m (MathSource) enables the production of - efficient C and Fortran code. *) - -(* :Context: Optimize` *) - -(* :Package Version: 1.2 *) - -(* :Copyright: Copyright 1993-4, Mark Sofroniou. - Permission is hereby granted to modify and/or make copies of - this file for any purpose other than direct profit, or as part - of a commercial product, provided this copyright notice is left - intact. Sale, other than for the cost of media, is prohibited. - - Permission is hereby granted to reproduce part or all of - this file, provided that the source is acknowledged. *) - -(* :History: - Modified Optimize syntax and improved performance, August 1994. - Significantly revised and publicly released May, 1994. - Original Version by Mark Sofroniou, September, 1993. *) - -(* :Keywords: - Assign, C, CAssign, Common, Compile, Cost, CForm, FORTRAN, - FortranAssign, FortranForm, Horner, Optimize, Optimization, - Polynomial, Sub-Expression, Syntactic. *) - -(* :Source: - Mark Sofroniou, Ph.D. Thesis, Loughborough University, Loughborough, - Leicestershire LE11 3TU, England. *) - -(* :Mathematica Version: 2.2 *) - -(* :Limitations: - This package enables syntactic optimization in linear time. - Optimization is used in the sense of improving the arithmetic - operation count and compactness of the code, rather than the - best possible. - Expressions containing non-binary arithmetic operators (Plus - and Times) are optimized by matching only entire sub-expressions - possibly after extracting numeric coefficients. - These operations rarely dominate the computational cost. - Options enable control over the optimization process. - The issue of numerical stability has not been addressed. *) - -BeginPackage["Optimize`"] - -Cost::usage = "Cost[expr,options] returns a list of Mathematica operators -present in expr, providing a count of the basic arithmetic operations and -function calls. Cost has attribute HoldAll so that arguments are maintained -in unevaluated form. Some examples of the behaviour (using options):\n -1+2*3 is counted as an addition and a multiplication. -a+b+c is counted as two additions and a*b*c as two multiplications.\n -Power[E,x] is counted as the exponential function Exp[x]. -Power[x,-1] is counted as a division.\n -Power[x,y] is calculated as y-1 multiplications for integer y." - -Horner::usage = "Horner[poly] puts the polynomial poly in Horner or nested form -with respect to the default variables Variables[poly]. This is an efficient form -for numerical evaluation. Horner[poly,vars] specifies the variable ordering -explicitly as the List vars. Multinomial conversion is applied recursively. -Horner factorisation of rational polynomials is possible as Horner[p1/p2] or -Horner[p1/p2,v1,v2]. Pre-conditioning of the coefficients is sometimes necessary -to attain numerical stability and more efficient methods exist in this case. -Horner's rule is optimal if no pre-conditioning is assumed, requiring n -multiplications and n additions to evaluate a polynomial of degree n." - -Optimize::usage = "Optimize[expr,opts]\n -Optimize transforms expr into a sequence of optimization statements and an -optimized expression. Optimization is performed in linear time (O(n) operations) -providing an efficient means of reducing the arithmetic operation count - not the -best possible. The optimization performed is mainly syntactic (only exact common -sub-expressions are matched) with a few heuristics. Options opts control the -type of optimizations performed and the format of the output." - - -(* Options: *) - -CostDivide::usage = "CostDivide is an option of Cost specifying whether -to consider negative exponentiations as divisions. The setting All -corresponds to the standard evaluation procedure and the compiler. -The setting Share counts only one division in an expression such as -(x^-2)*(y^-2) in correspondence with OuputForm and C and FORTRAN code -translation. CostDivide may evaluate to All, Share or False. -See also CostPower" - -CostExp::usage = "CostExp is an option of Cost specifying whether to -consider Power[E,x] as a call to the exponential function Exp[x]. -CostExp may evaluate to True or False." - -CostNull::usage = "CostNull is an option of Cost specifying a -list of symbols to be omitted for the purposes of costing. This is -useful, for example, for removing named arrays from consideration -(which have the same syntax as functions). -CostNull may evaluate to a (possibly empty) list of symbols." - -CostPower::usage="CostPower is an option of Cost specifying whether to -consider Power[x,y] as calculated using two function calls, namely -Exp[y*Log[x]]. The setting Integer assumes in addition that integer -powers are calculated using y-1 multiplications. This is particularly -useful when costing expressions optimized using the option -OptimizePower->Binary. CostPower may evaluate to Integer, True or False. -See also CostDivide" - -Binary::usage = "Binary is an argument of the option OptimizePower -specifying whether to perform repeated binary factoring of exponents." - -OptimizeCoefficients::usage = "OptimizeCoefficients is an option of -Optimize specifying whether to extract numerical coefficients in -expressions with head Plus and Times. OptimizeCoefficients may -evaluate to True or False." - -OptimizeFunction::usage = "OptimizeFunction is an option of Optimize -specifying whether to optimize common sub-expressions which are -neither atoms (?AtomQ) nor expressions with head Plus, Power, -or Times. OptimizeFunction may evaluate to True or False." - -OptimizeNull::usage = "OptimizeNull is an option of Optimize specifying -a list of symbols to be disregarded during the optimization process. -This is useful, for example, for specifying a named array (which has the -same syntax as a function). OptimizeNull may evaluate to a (possibly empty) -list of symbols." - -OptimizePlus::usage = "OptimizePlus is an option of Optimize specifying -whether to extract common sub-expressions with head Plus. Only exact -sub-expressions are matched. Numerical coefficients may be extracted -using OptimizeCoefficients. OptimizePlus may evaluate to True or False." - -OptimizePower::usage = "OptimizePower is an option of Optimize -specifying whether to optimize sub-expressions with head Power. -Integer and rational powers can be factored efficiently by repeated -binary exponentiation. OptimizePower may evaluate to True, False, -or Binary." - -OptimizeProcedure::usage = "OptimizeProcedure is an option of Optimize -specifying the form of the output returned. The setting True returns a -Module wrapped in Hold. Function, Compile and SetDelayed all know -about Optimize and make use of the option OptimizeProcedure (automatically -removing Hold). Applying ReleaseHold recovers the original (unoptimized) -expression. With the setting False a list of transformation -rules is returned together with the optimized expression. A definition -{optrules,optexpr} = Optimize[expr] enables the original -(unoptimized) expression to be recovered as optexpr //. Dispatch[optrules]." - -OptimizeTimes::usage = "OptimizeTimes is an option of Optimize specifying -whether to extract common sub-expressions with head Times. Only exact -sub-expressions are matched. Numerical coefficients may be extracted -using OptimizeCoefficients. OptimizeTimes may evaluate to True or False." - -OptimizeVariable::usage = "OptimizeVariable is an option of Optimize -specifying a pair {symb,form} where symb is the optimization variable -to introduce and form is the type of format to use. -The format can be either Sequence (o1,o2,...) or Array (o[1],o[2],...). -Consecutively numbered variables are returned." - -(* Default optimization symbol. *) - -(* -SetAttributes[o,NProtectedAll]; -*) - -o::usage="o is the default optimization symbol introduced by -Optimize when performing common sub-expression optimization." - -Unprotect[Binary, Cost, CostDivide, CostExp, CostNull, CostPower, Horner, o, -Optimize, OptimizeCoefficients, OptimizeFunction, OptimizeNull, OptimizePlus, -OptimizePower, OptimizeProcedure, OptimizeTimes, OptimizeVariable]; - -Begin["`Private`"] - -(* Function to check the data types of Cost options with error - messages. *) - -Options[Cost] = {CostDivide->All,CostExp->True, -CostNull->{Block,Module,With,CompoundExpression,Hold}, -CostPower->Integer}; - -Cost::args = "The `1` did not evaluate to `2`."; - -costerrmssgs = {{"option CostDivide","All, Share or False"}, -{"option CostExp","True or False"}, -{"option CostNull","a (possibly empty) list of symbols"}, -{"option CostAsAtom","Integer, True or False"}}; - -CostOptTest[opts___]:= - Module[{costdiv,costexp,costnull,costpower,types,defaults=Options[Cost]}, - - optlist = {costdiv,costexp,costnull,costpower} = - Map[First,defaults] /. {opts} /. defaults; - - types = {MatchQ[costdiv,All|Share|False], - MatchQ[costexp,True|False], - MatchQ[costnull,{___Symbol}], - MatchQ[costpower,Integer|True|False]}; - - Check[ - MapThread[ - If[#1,#1,Message[Cost::args,Apply[Sequence,#2]]]&, - {types,costerrmssgs} - ]; optlist, (* Return list of option values. *) - $Failed, (* Option of wrong type. *) - Cost::args (* Check only for these messages. *) - ] - ]; - - -SetAttributes[{Cost,MainCost},HoldAll]; - -Cost[expr_,opts___?OptionQ]:= - Module[{optvals}, - optvals /; - And[ - (optvals = CostOptTest[opts])=!=$Failed, - optvals = MainCost[Unevaluated[expr],Evaluate[optvals]]; True - ] - ]; - - - -(* Function for counting cost of basic operations. *) - -MainCost[expr_,{costdiv_,costexp_,costnull_,costpower_}]:= - Block[{CostFunction}, - -(* Prevent evaluation during costing. *) - - SetAttributes[{CostFunction},HoldAll]; - -(* Rules to cost expressions. *) - - If[costexp, CostFunction[Power[E,_]]:= Exp ]; - -(* Integer powers as multiplications and divisions. *) - - If[costpower===Integer, - If[MatchQ[costdiv,All|Share], - CostFunction[Power[_,y_Integer?Positive]]:= (y-1) Times; - CostFunction[Power[_,y_Integer?Negative]]:= Divide + (-y-1) Times, - CostFunction[Power[_,y_Integer]]:= (Abs[y]-1) Times (* No divisions. *) - ] - ]; - -(* x^y calculated as Exp[y Log[x]]. *) - - If[MatchQ[costpower,Integer|True], - If[MatchQ[costdiv,All|Share], - CostFunction[Power[_,_?Negative]]:= Divide+Exp+Log ]; - CostFunction[_Power]:= Exp+Log - ]; - -(* Add back in occurrances of shared divisions. *) - - If[costdiv===Share, - -(* Pure reciprocal case. *) - - Literal[CostFunction[e:Times[Power[_,_?Negative]..]]]:= - (1-Length[Unevaluated[e]]) Divide + (Length[Unevaluated[e]]-1) Times; - -(* Mixed case: (a*b..)/(c*d..). Similar to above, but don't count first - division as a multiplication. Hence: Length[num]-1+Length[denom]-1 = Length[e]-2. *) - - Literal[CostFunction[e:(Times[___,Power[_,_?Negative],___]..)]]:= - (1-Count[Unevaluated[e],Power[_,_?Negative]]) Divide + - (Length[Unevaluated[e]]-2) Times; - ]; - -(* Arithmetic operators. *) - - CostFunction[x_Plus]:= (Length[Unevaluated[x]]-1) Plus; - CostFunction[x_Times]:= (Length[Unevaluated[x]]-1) Times; - -(* Ignore Lists. *) - - CostFunction[x_List]:= 0; - -(* Cost remaining functions. *) - - CostFunction[x_]:= Head[Unevaluated[x]]; - -(* Cost required operators. *) - - DeleteCases[ - If[Head[#]===Plus,Apply[List,#],{#}]&[ - Apply[Plus, - Map[CostFunction, Level[Unevaluated[{expr}],-2,Unevaluated] ] - ] - ], - Apply[Alternatives,Times[costnull,_.]] (* Remove specified operands. *) - ] - - ]; (* End of MainCost *) - - - -(* Horner's rule. *) - -(* Default variables. *) - -Horner[p1_/p2_]:= - Block[{$RecursionLimit=Infinity}, - HornerRule[ Expand[p1], Variables[p1] ]/ - HornerRule[ Expand[p2], Variables[p2] ] - ]; - -Horner[ser_SeriesData]:= - Block[{$RecursionLimit=Infinity}, - HornerRule[Expand[#],Variables[#]]& @ Normal[ser] - ]; - -Horner[poly_]:= - Block[{$RecursionLimit=Infinity}, - HornerRule[Expand[poly],Variables[poly]] - ]; - - -(* Specified variables. *) - -Horner[p1_/p2_,varsp1_?VectorQ,varsp2_?VectorQ]:= - Block[{$RecursionLimit=Infinity}, - Times[ - HornerRule[ Expand[p1], varsp1 ], - Power[ HornerRule[ Expand[p2], varsp2 ], -1] - ] - ]; - -Horner[ser_SeriesData,vars_?VectorQ]:= - Block[{$RecursionLimit=Infinity}, - HornerRule[ Expand[Normal[ser]], vars] - ]; - -Horner[poly_,vars_?VectorQ]:= - Block[{$RecursionLimit=Infinity}, - HornerRule[ Expand[poly], vars ] - ]; - -Horner[poly_,var_]:= - Block[{$RecursionLimit=Infinity}, - HornerRule[ Expand[poly], var ] - ]; - -(* No variables found by Variables. *) - -HornerRule[poly_,{}]:= poly; - -HornerRule[0,_]:= 0; - -(* Horner's rule for multi-variate polynomials as recursive univariate - decomposition. *) - -HornerRule[poly_,{v_,rem__}]:= - Fold[ - SumCoeffs, - 0, - -(* Pair off variable with coefficients as {var^exp,hcoeff} where hcoeff - is the coefficient in Horner form. *) - - Thread[{ - v^Offsets[#], - Map[ HornerRule[#,{rem}]&, GetCoeffs[poly,v,#] ] - }]& @ Reverse[ Exponent[poly,v,Union[{##}]&] ] (* Powers sorted in descending order. *) - ]; - -(* Horner's rule for uni-variate polynomials. *) - -HornerRule[poly_,{v_}]:= - Fold[ - SumCoeffs, - 0, - -(* Pair off variable with coefficients as {var^exp,coeff}. *) - - Thread[{ - v^Offsets[#], GetCoeffs[poly,v,#] - }]& @ Reverse[ Exponent[poly,v,Union[{##}]&] ] (* Powers sorted in descending order. *) - ]; - -(* Calculate offset powers as successive differences (not necessarily - incremental). *) - -Offsets[e:{_}]:= e; -Offsets[e_]:= Join[-Drop[e,1],{0}] + e; - -(* Accumulate Horner form of polynomial. *) - -SumCoeffs[sum_,{v_,c_}]:= v (c + sum); - -(* Can eventually be replaced by CoefficientList when bugs are fixed. *) - -SetAttributes[GetCoeffs,Listable]; -GetCoeffs[c_,v_,0]:= Coefficient[c,v,0]; -GetCoeffs[c_,v_,p_]:= Coefficient[c,v^p]; - - -(* Optimization of expressions. *) - -(* Function to check the data types of Optimize options with - error messages. *) - -Options[Optimize] = {OptimizeCoefficients->False,OptimizeFunction->True, -OptimizeNull->{List},OptimizePlus->True,OptimizePower->True, -OptimizeProcedure->False,OptimizeTimes->True,OptimizeVariable->{o,Sequence}}; - -Optimize::args = "The `1` did not evaluate to `2`."; - -opterrmsgs = { - {"option OptimizeCoefficients","True or False"}, - {"option OptimizeFunction","True or False"}, - {"option OptimizeNull","a (possibly empty) list of symbols"}, - {"option OptimizePlus","True or False"}, - {"option OptimizePower","True, False or Binary"}, - {"option OptimizeProcedure","True or False"}, - {"option OptimizeTimes","True or False"}, - {"option OptimizeVariable","a pair of the form {Symbol,Sequence|Array}"}}; - -OptimizeOptionTest[opts___]:= - Module[{defaults,types,optcoeff,optfunc,optlist,optnull,optplus, - optpower,optproc,opttimes,optoutv}, - - defaults = Options[Optimize]; - - optlist = {optcoeff,optfunc,optnull,optplus,optpower,optproc,opttimes, - optoutv} = Map[First,defaults] /. {opts} /. defaults; - - types = { - MatchQ[optcoeff,True|False], - MatchQ[optfunc,True|False], - MatchQ[optnull,{___Symbol}], - MatchQ[optplus,True|False], - MatchQ[optpower,Binary|True|False], - MatchQ[optproc,True|False], - MatchQ[opttimes,True|False], - MatchQ[optoutv,{_Symbol,Sequence|Array}]}; - - Check[ - MapThread[ - If[#1,#1,Message[Optimize::args,Apply[Sequence,#2]]]&, - {types,opterrmsgs} - ]; optlist, (* Return list of option values. *) - $Failed, (* Option of wrong type. *) - Optimize::args (* Check only for these messages. *) - ] - ]; (* End of OptimizeOptionTest. *) - - -(* Default optimization symbol is o. *) - -Optimize[expr_,opts___?OptionQ]:= - Module[{optvals}, - optvals /; - And[ - (optvals = OptimizeOptionTest[opts])=!=$Failed, - optvals = MainOptimize[expr,optvals]; True - ] - ]; - - - -(* Define optimization function. All expression heads are wrapped - in Hold to prevent re-evaluation during the optimization process. *) - -MainOptimize[expr_,{optcoeff_,optfunc_,optnull_,optplus_,optpower_,optproc_, - opttimes_,{outvar_,outform_}}]:= -Block[{$RecursionLimit=Infinity,BinPowDecomp,BinPowRule,downvals,HoldHead, - index=0,keeprules,MakeBins,MakeRule,outvar,OptRule,optexpr,optvar}, - -(* Store sub-expression by making a new DownValue for OptRule and replace - expression by a new optimization variable. *) - - SetAttributes[MakeRule,HoldAll]; - - MakeRule[e_]:= (OptRule[e]:=OptRule[e]=#; #)& @ optvar[++index]; - -(* Rules for reciprocals. *) - - If[optpower=!=False, - OptRule[e:Hold[Power][_,-1]]:= MakeRule[e]; - - OptRule[Hold[Power][x_,y_?(NumberQ[#]&&Negative[#]&)]]:= - OptRule[Hold[Power][OptRule[Hold[Power][x,-y]],-1]]; - ]; - -(* Rules to binary factor integer and rational powers. *) - - Switch[optpower, - -(* Rules for binary decomposition of positive powers. *) - - Binary, - -(* Save all binary power optimizations (even if they only occur once). *) - - MakeBins[e_]:= OptRule[e]=optvar[++index]; - -(* Calculate and store binary power decomposition. *) - - BinPowDecomp[p_]:= BinPowDecomp[p]= - 2^(-1+Flatten[Position[Reverse[IntegerDigits[p,2]],1]]); - -(* Decompose as products of binary powers. *) - - binprod[{p_}]:= p; - binprod[{p__}]:= OptRule[ Hold[Times][p] ]; - - OptRule[Hold[Power][x_,p_Integer]]:= - binprod[ BinPowRule[x,BinPowDecomp[p]] ]; - - SetAttributes[BinPowRule,{Listable}]; - -(* Binary power stopping criterion. *) - - BinPowRule[e_,1]:= e; - -(* Store all binary powers. *) - - BinPowRule[x_,p_]:= - BinPowRule[x,p]= - MakeBins[ Hold[Power][BinPowRule[x,Quotient[p,2]],2] ]; - -(* Rule for positive rational power with unit numerator. *) - - OptRule[e:Hold[Power][_,Rational[1,_]]]:= MakeBins[e]; - -(* Rules to binary decompose fractional exponents. *) - - OptRule[Hold[Power][x_,Rational[y_,z_]]]:= - If[y<z, - -(* Rule to decompose rational power as: x^(n/m) -> (x^(1/m))^n. *) - - OptRule[ Hold[Power][ OptRule[Hold[Power][x,Rational[1,z]]] ,y] ], - -(* Rule to decompose rational power as: x^(c/d) -> x^(q+r/d) -> x^(q)*x^(r/d). *) - - OptRule[ - Hold[Times][ - OptRule[ Hold[Power][ x, Quotient[y,z] ] ], - OptRule[ Hold[Power][ OptRule[Hold[Power][x,Rational[1,z]]], Mod[y,z]] ] - ] - ] - ]; - -(* Rule for general (non-numeric) and non-binary powers. *) - - OptRule[e:Blank[Hold[Power]]]:= MakeRule[e], - -(* Rule for literal powers. *) - - True, - - OptRule[e:Blank[Hold[Power]]]:= MakeRule[e], - - False, - - OptRule[e:Blank[Hold[Power]]]:= e (* Else disregard head Power. *) - ]; - - -(* Plus rules. *) - - If[optplus, - -(* Extract numerical coefficient. *) - - If[optcoeff, - OptRule[ Hold[Plus][n_?NumberQ,x_,y__] ]:= - OptRule[ Hold[Plus][ n, OptRule[Hold[Plus][x,y]] ] ] - ]; - -(* Store literal sub-expression. *) - - OptRule[e:Blank[Hold[Plus]]]:= MakeRule[e], - - OptRule[e:Blank[Hold[Plus]]]:= e (* Else disregard head Plus. *) - - ]; (* End of rules for expressions with head Plus. *) - - -(* Times rules. *) - - If[opttimes, - -(* Extract numerical coefficient. *) - - If[optcoeff, - OptRule[ Hold[Times][n_?NumberQ,x_,y__] ]:= - OptRule[ Hold[Times][ n, OptRule[Hold[Times][x,y]] ] ] - ]; - -(* Store literal sub-expression. *) - - OptRule[e:Blank[Hold[Times]]]:= MakeRule[e], - - OptRule[e:Blank[Hold[Times]]]:= e (* Else disregard head Times. *) - - ]; (* End of rules for expressions with head Times. *) - - -(* Rule for functions. *) - - If[optfunc, - OptRule[e_]:= MakeRule[e], - - OptRule[e_]:= e (* Else disregard remaining expressions. *) - ]; - - -(* Don't make rules for specified symbols, but enable argument evaluation. *) - - Map[(HoldHead[e:Blank[#]]:= Operate[Hold,e])&, optnull ]; - -(* Wrap up function heads to prevent re-evaluation. Enables held - functions arguments to be evaluated. *) - - HoldHead[e_]:= OptRule[Operate[Hold,e]]; - -(* Map function for creating optimization rules at non-atomic levels. - Include special case of binary power factorisation. *) - - optexpr = If[ Head[expr]===Power&&optpower===Binary, HoldHead[#], # ]& @ - Map[ HoldHead, expr, -2 ]; - - -(* Find occurances of temporary variables in optimized expression and - in optimization rules. Decide which temporary variables to keep and - which to discard by pattern look-up. *) - - downvals = DownValues[OptRule]; - -(* Discard rules by assigning optimization variable to subexpression. - Head replacement makes this efficient by avoiding re-evaluation. *) - - Cases[downvals,Literal[_:>(OptRule[rhs_]=lhs_)]:>(lhs=rhs)]; - -(* Substitute discarded rules into optimized expression. *) - - optexpr = optexpr; - -(* Save ordered list of repeated rules. *) - - keeprules = - Sort[ - Cases[downvals,Literal[_[OptRule[rhs:_]]:>lhs:_optvar]:>lhs->rhs] - ]; - - -(* Format optimization variable. Reset index for consecutively numbered - output variables. *) - - index = 0; - - If[outform===Sequence, - -(* Write optimization variable as consecutive symbols o1, o2, ... - Store indexing string function (not localised). *) - - istr[i_]:= istr[i]=ToString[i]; - (optvar[i_]:= optvar[i]=ToExpression[#<>istr[++index]])& @ ToString[outvar], - -(* Write optimization variable as a consecutive array o[1], o[2], ... *) - - optvar[i_]:= optvar[i]=outvar[++index] - ]; - - -(* Create a Module or leave as a list of replacement rules and - apply formatting rules for optimization variable. *) - - If[keeprules==={}, - If[optproc, expr, {{},expr} ], (* No optimization performed *) - If[optproc, - MakeProc[outvar,outform,ReleaseHold[keeprules],ReleaseHold[optexpr]], - ReleaseHold[ {keeprules, optexpr} ] - ] - ] - -]; (* End of MainOptimize.*) - - - -(* Convert output to a held module. *) - -MakeProc[optvar_,optform_,{optseq__Rule},optexpr_]:= - (Hold[ Module[#,optseq; optexpr] ] /. Rule->Set)& @ - If[optform===Sequence, First[Thread[{optseq},Rule]], {optvar}]; - -(* Create an optimized procedure. *) - -RemoveHold[expr_Hold]:= Apply[Unevaluated,expr]; -RemoveHold[expr_]:= Unevaluated[expr]; - -(* SetDelayed definition. *) - -Optimize /: - (lhs_:=Optimize[expr_,opts___?OptionQ]):= - (Unevaluated[lhs]:=#)& @ - RemoveHold[ Optimize[expr,OptimizeProcedure->True,opts] ]; - -(* Compile definition. *) - -Optimize /: - Compile[args_,Optimize[expr_,opts___?OptionQ],info___]:= - Compile[Unevaluated[args],#,info]& @ - RemoveHold[ Optimize[expr,OptimizeProcedure->True,opts] ]; - -(* Function definition. *) - -Optimize /: - Function[args_,Optimize[expr_,opts___?OptionQ],attr___]:= - ReleaseHold[ - Function[ - args, - Evaluate[Optimize[expr,OptimizeProcedure->True,opts]], - attr - ] - ]; - - -End[]; (* End `Private` Context. *) - -(* Protect exported symbols. *) - -SetAttributes[{Cost,Horner,Optimize},ReadProtected]; - -Protect[Binary, Cost, CostDivide, CostExp, CostNull, CostPower, Horner, o, -Optimize, OptimizeCoefficients, OptimizeFunction, OptimizeNull, OptimizePlus, -OptimizePower, OptimizeProcedure, OptimizeTimes, OptimizeVariable]; - -EndPackage[]; (* End package Context. *) |