aboutsummaryrefslogtreecommitdiff
path: root/Tools/CodeGen/Differencing.m
diff options
context:
space:
mode:
Diffstat (limited to 'Tools/CodeGen/Differencing.m')
-rw-r--r--Tools/CodeGen/Differencing.m253
1 files changed, 160 insertions, 93 deletions
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_];
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-30 18:16:50 +0000