Merge eschnett/dgfe into master

Auxiliary/Cactus/KrancNumericalTools/GenericFD/param.ccl

@@ -32,12 +32,33 @@ CCTK_INT  boundary_width "width of boundary (fix later to use Cactus boundary ca
 } 1
 
 restricted:
+CCTK_INT assume_stress_energy_state "Assume stress_energy_state has a particular value" STEERABLE=RECOVER
+{
+  -1  :: "do not assume anything"
+  0:1 :: "assume off or on"
+} -1
+
+CCTK_INT assume_use_jacobian "Assume use_jacobian has a particular value" STEERABLE=RECOVER
+{
+  -1  :: "do not assume anything"
+  0:1 :: "assume off or on"
+} -1
+
 CCTK_STRING jacobian_group "Name of group containing Jacobian" STEERABLE=RECOVER
 {
   "" :: "String of the form <implementation>::<groupname>"
 } ""
 
-restricted:
+CCTK_STRING jacobian_determinant_group "Name of group containing Jacobian determinant" STEERABLE=RECOVER
+{
+  "" :: "String of the form <implementation>::<groupname>"
+} ""
+
+CCTK_STRING jacobian_inverse_group "Name of group containing Jacobian inverse" STEERABLE=RECOVER
+{
+  "" :: "String of the form <implementation>::<groupname>"
+} ""
+
 CCTK_STRING jacobian_derivative_group "Name of group containing Jacobian derivative" STEERABLE=RECOVER
 {
   "" :: "String of the form <implementation>::<groupname>"

Auxiliary/Cactus/KrancNumericalTools/GenericFD/src/GenericFD.h

@@ -32,6 +32,7 @@
 #  include <cmath>
 #endif
 #include <math.h>
+#include <sys/time.h>
 
 #include <cctk.h>
 

Tools/CodeGen/CalculationBoundaries.m

@@ -46,8 +46,9 @@
    implemented. *)
 
 CalculationBoundariesFunction[calc_List] :=
-  Module[{selectGroup, groups, groupNamesSet, imp},
+  Module[{funcName, selectGroup, groups, groupNamesSet, imp},
 
+  funcName = lookup[calc, Name];
   imp = lookup[calc,Implementation];
 
   (* If the calculation sets every grid point, we don't need to apply
@@ -69,8 +70,10 @@
     groups = lookup[calc, Groups];
     groupNamesSet = qualifyGroupName[#, imp] & /@ groupsSetInCalc[calc, groups];
 
-    DefineCCTKFunction[lookup[calc, Name] <> "_SelectBCs", "void",
+    DefineCCTKFunction[funcName <> "_SelectBCs", "void",
     {
+      "if (cctk_iteration % " <> funcName <> "_calc_every != " <> funcName <> "_calc_offset)\n",
+      "  return;\n",
       DefineVariable["ierr",   "CCTK_INT", "0"],
 (*      DefineVariable["table",  "CCTK_INT", "-1"],*)
       Map[selectGroup, groupNamesSet],

Tools/CodeGen/CalculationFunction.m

@@ -193,7 +193,7 @@ pathalogical enough (e.g. {s1 -> s2, s2 -> s1} would not be
     ignoreGroups = {"TmunuBase::stress_energy_scalar", "TmunuBase::stress_energy_vector",
       "TmunuBase::stress_energy_tensor"};
     groupNames2 = Select[groupNames, !MemberQ[ignoreGroups, #] &];
-    {"\nconst char *const groups[] = {\n  ",
+    {"\nconst char* const groups[] = {\n  ",
     Riffle[Map[Quote,groupNames2], ",\n  "],
     "};\n",
     "GenericFD_AssertGroupStorage(cctkGH, ", Quote[calcName],", ", Length[groupNames2], ", groups);\n"}];
@@ -487,7 +487,7 @@ pathalogical enough (e.g. {s1 -> s2, s2 -> s1} would not be
   (* Check that there are no unknown symbols in the calculation *)
   allSymbols = calculationSymbols[cleancalc];
   knownSymbols = Join[lookupDefault[cleancalc, AllowedSymbols, {}], gfs, shorts, parameters,
-    {dx,dy,dz,dt,idx,idy,idz,t, Pi, E, Symbol["i"], Symbol["j"], Symbol["k"], normal1, normal2,
+    {dx,dy,dz,dt,idx,idy,idz,t, usejacobian, Pi, E, Symbol["i"], Symbol["j"], Symbol["k"], normal1, normal2,
     normal3, tangentA1, tangentA2, tangentA3, tangentB1, tangentB2, tangentB3},
     If[useJacobian, JacobianSymbols[], {}]];
 
@@ -563,7 +563,7 @@ pathalogical enough (e.g. {s1 -> s2, s2 -> s1} would not be
            "    enum {nequations = nvars};",
            "    enum {nexternal = 3*nvars};",
            "    enum {nbitmasks = 0};",
-           "    static bool const pure = false;",
+           "    static const bool pure = false;",
            "  };",
            "} // namespace",
            "",
@@ -574,7 +574,7 @@ pathalogical enough (e.g. {s1 -> s2, s2 -> s1} would not be
            "  typedef hrscc::CLaw<DGFE_"<>name<>"> claw;",
            "  typedef hrscc::traits<DGFE_"<>name<>">::state_t state_t;",
            "  typedef hrscc::traits<DGFE_"<>name<>"> variables_t;",
-           "  static int const nvars = variables_t::nvars;",
+           "  static const int nvars = variables_t::nvars;",
            "  ",
            "  DGFE_"<>name<>"();",
            "  ",
@@ -596,7 +596,7 @@ pathalogical enough (e.g. {s1 -> s2, s2 -> s1} would not be
              "    }"},
                  {dir, 1, 3}],
            "    default:",
-           "      assert(0);",
+           "      CCTK_BUILTIN_UNREACHABLE();",
            "    }",
            "    ",
            Map["    observer.flux[dir][variables_t::i"<>ToString[#]<>"] = flux"<>ToString[#]<>"L;" &, vars],
@@ -628,10 +628,10 @@ pathalogical enough (e.g. {s1 -> s2, s2 -> s1} would not be
            "",
            "",
            "namespace {",
-           "  int varindex(char const* const varname)",
+           "  int varindex(const char* const varname)",
            "  {",
-           "    int const vi = CCTK_VarIndex(varname);",
-           "    if (vi<0) CCTK_WARN(CCTK_WARN_ABORT, \"Internal error\");",
+           "    const int vi = CCTK_VarIndex(varname);",
+           "    if (vi<0) CCTK_ERROR(\"Internal error\");",
            "    return vi;",
            "  }",
            "}",
@@ -661,9 +661,9 @@ pathalogical enough (e.g. {s1 -> s2, s2 -> s1} would not be
            "#undef PDstandardNth1",
            "#undef PDstandardNth2",
            "#undef PDstandardNth3",
-           "#define PDstandardNth1(u) (solver->diff<hrscc::policy::x>(&(u)[-index], i,j,k))",
-           "#define PDstandardNth2(u) (solver->diff<hrscc::policy::y>(&(u)[-index], i,j,k))",
-           "#define PDstandardNth3(u) (solver->diff<hrscc::policy::z>(&(u)[-index], i,j,k))",
+           "#define PDstandardNth1(u) (solver->wdiff<hrscc::policy::x>(&(u)[-index], i,j,k))",
+           "#define PDstandardNth2(u) (solver->wdiff<hrscc::policy::y>(&(u)[-index], i,j,k))",
+           "#define PDstandardNth3(u) (solver->wdiff<hrscc::policy::z>(&(u)[-index], i,j,k))",
            "",
            "",
            ""
@@ -707,7 +707,7 @@ pathalogical enough (e.g. {s1 -> s2, s2 -> s1} would not be
     (* We could (or probably should) write this into a source file of its own *)
     If[OptionValue[UseOpenCL],
        {
-         "char const *const source =\n"
+         "const char* const source =\n"
        },
        {
        }],
@@ -750,13 +750,13 @@ pathalogical enough (e.g. {s1 -> s2, s2 -> s1} would not be
            groupNames = groupsInCalculation[cleancalc, imp];
            groupNames = Select[groupNames, !MemberQ[ignoreGroups, #] &];
            {
-             "char const *const groups[] = {\n  ",
+             "const char* const groups[] = {\n  ",
                Riffle[Join[Map[Quote, groupNames], {"NULL"}], ",\n  "],
                "};\n\n"
            }
          ],
          "static struct OpenCLKernel *kernel = NULL;\n",
-         "char const *const sources[] = {differencing, source, NULL};\n",
+         "const char* const sources[] = {differencing, source, NULL};\n",
          "OpenCLRunTime_CallKernel(cctkGH, CCTK_THORNSTRING, \"" <> functionName <> "\",\n",
          "                         sources, groups, NULL, NULL, NULL, -1,\n",
          "                         imin, imax, &kernel);\n\n"
@@ -937,7 +937,7 @@ it is only possible to do this outside all if(){} statements. *)
     assignLocalFunctions[gs_, useVectors_, useJacobian_, NameFunc_] :=
       Module[{conds, varPatterns, varsInConds, simpleVars, code},
         conds =
-          {{"eT" ~~ _ ~~ _, "*stress_energy_state", "ToReal(0.0)"}}; (* This should be passed as an option *)
+          {{"eT" ~~ _ ~~ _, "assume_stress_energy_state>=0 ? assume_stress_energy_state : *stress_energy_state", "ToReal(0.0)"}}; (* This should be passed as an option *)
         If[useJacobian,
           conds = Append[conds, JacobianConditionalGridFunctions[]]];
 
@@ -1015,7 +1015,7 @@ it is only possible to do this outside all if(){} statements. *)
                     gfsInLHS] }],
             OptionValue[UseVectors],
             CommentedBlock["Copy local copies back to grid functions",
-              { PrepareStorePartialVariableInLoop["i", "kimin", "kimax"],
+              { PrepareStorePartialVariableInLoop["i", "vecimin", "vecimax"],
                 Map[StorePartialVariableInLoop[gridName[#], localName[#]] &,
                     gfsInLHS] }],
             True,
@@ -1043,7 +1043,7 @@ it is only possible to do this outside all if(){} statements. *)
           "Calculate temporaries and grid functions", 
           If[OptionValue[UseVectors],
              {
-               PrepareStorePartialVariableInLoop["i", "kimin", "kimax"],
+               PrepareStorePartialVariableInLoop["i", "vecimin", "vecimax"],
                Map[StorePartialVariableInLoop[FlattenBlock@gridName[#[[1]]], #[[2]]] &, eqs2]
              },
              Map[

Tools/CodeGen/CodeGenC.m

@@ -103,7 +103,7 @@
 DefFn[
   DeclareVariableNoInit[name:(_String|_Symbol), type_String] :=
   If[SOURCELANGUAGE == "C",
-     {type, " ",    name, " CCTK_ATTRIBUTE_UNUSED ",EOL[]},
+     {type, " ",    name, " CCTK_ATTRIBUTE_UNUSED",EOL[]},
      {type, " :: ", name, EOL[]} (* no value init here to avoid implicit SAVE attribute *)]];
 
 DefFn[
@@ -140,15 +140,15 @@
 
 DefFn[
   DefineVariable[name:(_String|_Symbol), type_String, value:CodeGenBlock] :=
-  {type, " ", name, " CCTK_ATTRIBUTE_UNUSED ", " = ", value, EOL[]}];
+  {type, " ", name, " CCTK_ATTRIBUTE_UNUSED", " = ", value, EOL[]}];
 
 DefFn[
   AssignVariable[dest:(_String|_Symbol), src:CodeGenBlock] :=
   {dest, " = ", src, EOL[]}];
 
 DefFn[
   DeclareAssignVariable[type_String, dest:(_String|_Symbol), src:CodeGenBlock] :=
-  {type, " /*const*/ ", dest, " CCTK_ATTRIBUTE_UNUSED ", " = ", src, EOL[]}];
+  {"const ", type, " ", dest, " CCTK_ATTRIBUTE_UNUSED", " = ", src, EOL[]}];
 
 (* comments are always done C-style because they are killed by cpp anyway *) 
 DefFn[
@@ -215,13 +215,14 @@
 
 DefFn[
   switchOption[{value:(_String|_Symbol|_?NumberQ), block:CodeGenBlock}] :=
-  {"case ", value, ":\n", IndentBlock[{block,"break;\n"}]}];
+  {"case ", value, ":\n", CBlock[{block,"break;\n"}]}];
 (* Outer list unnecessary? *)
 
 DefFn[
   SwitchStatement[var:(_String|_Symbol), pairs__] :=
-  {"switch(", var, ")\n",
-   CBlock[{Riffle[Map[switchOption, {pairs}],"\n"]}]}];
+  {"switch (", var, ")\n",
+   CBlock[{Riffle[Map[switchOption, {pairs}],"\n"],
+           "default:\n", IndentBlock[{"CCTK_BUILTIN_UNREACHABLE();\n"}]}]}];
 
 DefFn[
   Conditional[condition:CodeGenBlock, block:CodeGenBlock] :=

Tools/CodeGen/CodeGenCactus.m

@@ -112,7 +112,7 @@
 
 DefFn[
   DeclareAssignVariableInLoop[type_String, dest:(_String|_Symbol), src:(_String|_Symbol)] :=
-  {type, " /*const*/ ", dest, " CCTK_ATTRIBUTE_UNUSED = vec_load(", src, ")", EOL[]}];
+  {"const ", type, dest, " CCTK_ATTRIBUTE_UNUSED = vec_load(", src, ")", EOL[]}];
 
 DefFn[
   MaybeAssignVariableInLoop[dest:(_String|_Symbol), src:(_String|_Symbol), cond:Boolean] :=
@@ -384,20 +384,50 @@
   If[SOURCELANGUAGE == "C",  
      CommentedBlock[
        "Loop over the grid points",
-       {"#pragma omp parallel\n",
-        If[vectorise, "CCTK_LOOP3STR", "CCTK_LOOP3"],
-        "(", functionName, ",\n",
-        "  i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],\n",
-        "  cctk_ash[0],cctk_ash[1],cctk_ash[2]",
-        If[vectorise, {",\n", "  kimin,kimax, CCTK_REAL_VEC_SIZE"}, ""],
-        ")\n",
-        "{\n",
-        IndentBlock[
-          {(* DeclareVariable["index", "// int"], *)
-           (* DeclareAssignVariable["int", "index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"], *)
-           DeclareAssignVariable["ptrdiff_t", "index", "di*i + dj*j + dk*k"],
-           block}], "}\n",
-        If[vectorise, "CCTK_ENDLOOP3STR", "CCTK_ENDLOOP3"] <> "(", functionName, ");\n"}],
+       { "/* Circumvent a compiler bug on Blue Gene/Q */\n",
+         "const int imin0=imin[0];\n",
+         "const int imin1=imin[1];\n",
+         "const int imin2=imin[2];\n",
+         "const int imax0=imax[0];\n",
+         "const int imax1=imax[1];\n",
+         "const int imax2=imax[2];\n",
+         "// #undef VEC_COUNT\n",
+         "// #define VEC_COUNT(x) x\n",
+         "// double vec_iter_timer;\n",
+         "// {\n",
+         "//   timeval tv;\n",
+         "//   gettimeofday(&tv, NULL);\n",
+         "//   vec_iter_timer = -(tv.tv_sec + 1.0e-6 * tv.tv_usec);\n",
+         "// }\n",
+         "// ptrdiff_t vec_iter_counter = 0;\n",
+         "// ptrdiff_t vec_op_counter = 0;\n",
+         "// ptrdiff_t vec_mem_counter = 0;\n",
+         "#pragma omp parallel // reduction(+: vec_iter_counter, vec_op_counter, vec_mem_counter)\n",
+         If[vectorise, "CCTK_LOOP3STR", "CCTK_LOOP3"],
+         "(", functionName, ",\n",
+         "  i,j,k, imin0,imin1,imin2, imax0,imax1,imax2,\n",
+         "  cctk_ash[0],cctk_ash[1],cctk_ash[2]",
+         If[vectorise, {",\n", "  vecimin,vecimax, CCTK_REAL_VEC_SIZE"}, ""],
+         ")\n",
+         "{\n",
+         IndentBlock[
+           {(* DeclareVariable["index", "// int"], *)
+            (* DeclareAssignVariable["int", "index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"], *)
+            DeclareAssignVariable["ptrdiff_t", "index", "di*i + dj*j + dk*k"],
+            If[vectorise,
+               "// vec_iter_counter+=CCTK_REAL_VEC_SIZE;\n",
+               "// ++vec_iter_counter;\n"],
+            block}],
+         "}\n",
+         If[vectorise, "CCTK_ENDLOOP3STR", "CCTK_ENDLOOP3"] <> "(", functionName, ");\n",
+         "// {\n",
+         "//   timeval tv;\n",
+         "//   gettimeofday(&tv, NULL);\n",
+         "//   vec_iter_timer += tv.tv_sec + 1.0e-6 * tv.tv_usec;\n",
+         "// }\n",
+         "// CCTK_VInfo(CCTK_THORNSTRING, \"function="<>functionName<>" time=%g points=%td fp_ops=%td mem_ops=%td\", vec_iter_timer, vec_iter_counter, vec_op_counter, vec_mem_counter);\n",
+         "// #undef VEC_COUNT\n",
+         "// #define VEC_COUNT(x)\n"}],
      (* else *)
      ""]];
 
@@ -576,15 +606,21 @@
 
     (* Optimise *)
     expr = expr //. {
-      kneg[ToReal[a_]]    -> ToReal[-a],
-      kmul[ToReal[-1],x_] -> kneg[x],
-      kmul[x_,ToReal[-1]] -> kneg[x],
-      kneg[kneg[x_]]      -> x,
+      kneg[ToReal[a_]]      -> ToReal[-a],
+      kmul[ToReal[-1],x_]   -> kneg[x],
+      kmul[ToReal[-1.0],x_] -> kneg[x],
+      kmul[x_,ToReal[-1]]   -> kneg[x],
+      kmul[x_,ToReal[-1.0]] -> kneg[x],
+      kneg[kneg[x_]]        -> x,
 
       kadd[ToReal[0],x_]             -> x,
+      kadd[ToReal[0.0],x_]           -> x,
       kadd[x_,ToReal[0]]             -> x,
+      kadd[x_,ToReal[0.0]]           -> x,
       ksub[ToReal[0],x_]             -> kneg[x],
+      ksub[ToReal[0.0],x_]           -> kneg[x],
       ksub[x_,ToReal[0]]             -> x,
+      ksub[x_,ToReal[0.0]]           -> x,
       kadd[kneg[x_],y_]              -> ksub[y,x],
       ksub[kneg[x_],y_]              -> kneg[kadd[x,y]],
       kadd[x_,kneg[y_]]              -> ksub[x,y],
@@ -602,11 +638,21 @@
            kadd[ToReal[a_],y_]]      -> kadd[ToReal[a],kadd[x,y]],
 
       kmul[ToReal[0],x_]             -> ToReal[0],
+      kmul[ToReal[0.0],x_]           -> ToReal[0],
       kmul[x_,ToReal[0]]             -> ToReal[0],
+      kmul[x_,ToReal[0.0]]           -> ToReal[0],
       kmul[ToReal[+1],x_]            -> x,
+      kmul[ToReal[+1.0],x_]          -> x,
       kmul[x_,ToReal[+1]]            -> x,
+      kmul[x_,ToReal[+1.0]]          -> x,
+      kmul[ToReal[-1],x_]            -> kneg[x],
+      kmul[ToReal[-1.0],x_]          -> kneg[x],
+      kmul[x_,ToReal[-1]]            -> kneg[x],
+      kmul[x_,ToReal[-1.0]]          -> kneg[x],
       kdiv[ToReal[0],x_]             -> ToReal[0],
+      kdiv[ToReal[0.0],x_]           -> ToReal[0],
       (* kdiv[x_,ToReal[0]]           -> ToReal[nan], *)
+      (* kdiv[x_,ToReal[0.0]]         -> ToReal[nan], *)
       kdiv[x_,ToReal[y_]]            -> kmul[x,ToReal[1/y]],
       kdiv[x_,kdiv[y_,z_]]           -> kdiv[kmul[x,z],y],
       kdiv[kdiv[x_,y_],z_]           -> kdiv[x,kmul[y,z]],
@@ -749,8 +795,6 @@
            *)
 
         (* Handle Piecewise function *)
-        (* TODO: This does not work with vectorisation, since IfThen
-           there expects a constant condition *)
         rhs = rhs /. Piecewise -> piecewise1
                   //. piecewise1[pairs_List, val_:0] :>
                          If[pairs==={}, val,
@@ -764,7 +808,7 @@
         (* Avoid rational numbers *)
         rhs = rhs /. xx_Rational :> N[xx, 30];
         (* Avoid integers *)
-        rhs = rhs /. xx_Integer :> 1.0*xx;
+        (* rhs = rhs /. xx_Integer :> 1.0*xx; *)
 
         (* Simple optimisations *)
         rhs = rhs /. IfThen[_, aa_, aa_] -> aa;