diff options
Diffstat (limited to 'Tools/CodeGen/CalculationFunction.m')
-rw-r--r-- | Tools/CodeGen/CalculationFunction.m | 340 |
1 files changed, 112 insertions, 228 deletions
diff --git a/Tools/CodeGen/CalculationFunction.m b/Tools/CodeGen/CalculationFunction.m index 89a5fb3..1029925 100644 --- a/Tools/CodeGen/CalculationFunction.m +++ b/Tools/CodeGen/CalculationFunction.m @@ -1,6 +1,4 @@ -(* $Id$ *) - (* Copyright 2004 Sascha Husa, Ian Hinder, Christiane Lechner This file is part of Kranc. @@ -23,8 +21,9 @@ BeginPackage["sym`"]; {GridFunctions, Shorthands, Equations, t, DeclarationIncludes, -LoopPreIncludes, GroupImplementations, PartialDerivatives, Dplus1, NoSimplify, -Dplus2, Dplus3, Boundary, Interior, InteriorNoSync, Where, AddToStencilWidth, Everywhere, normal1, normal2, normal3} +LoopPreIncludes, GroupImplementations, PartialDerivatives, Dplus1, +NoSimplify, Dplus2, Dplus3, Boundary, Interior, InteriorNoSync, Where, +AddToStencilWidth, Everywhere, normal1, normal2, normal3} {INV, SQR, CUB, QAD, dot, pow, exp,dx,dy,dz, idx, idy, idz} @@ -37,9 +36,12 @@ BeginPackage["CalculationFunction`", {"CodeGen`", "sym`", "MapLookup`", "KrancGr CreateCalculationFunction::usage = ""; calculationUsedGroups::usage = ""; -UncommentSourceSync::usage = "UncommentSourceSync[source_struct] uncomments calls to CCTK_SyncGroup."; -GrepSyncGroups::usage = "GrepSyncGroups[source_struct] 'greps' a list of groups to be SYNC'ed from" <> - "code produced by CreateCalculationFunction." +UncommentSourceSync::usage = + "UncommentSourceSync[source_struct] uncomments calls to CCTK_SyncGroup."; + +GrepSyncGroups::usage = + "GrepSyncGroups[source_struct] 'greps' a list of groups to be SYNC'ed from code " <> + "produced by CreateCalculationFunction." VerifyCalculation::usage = ""; allVariables::usage = ""; @@ -50,16 +52,13 @@ Begin["`Private`"]; General Utility Functions -------------------------------------------------------------------------- *) +removeRHS[x_] := + Module[{string = ToString[x]}, + If[StringMatchQ[string, "*rhs"], + ToExpression@StringDrop[string, -3], + x]]; -(* remove the 'rhs' at the end of a symbol or string *) -removeRHS[x_] := Module[{string = ToString[x]}, - - If[StringMatchQ[string, "*rhs"], - ToExpression@StringDrop[string, -3], - x ] - ]; - -(* collect and simplify terms *) +(* Collect and simplify terms *) simpCollect[collectList_, eqrhs_, localvar_, debug_] := Module[{rhs, collectCoeff, all, localCollectList}, InfoMessage[InfoFull, localvar]; @@ -77,9 +76,7 @@ simpCollect[collectList_, eqrhs_, localvar_, debug_] := all = Collect[rhs, localCollectList, Simplify]; InfoMessage[InfoFull, "ByteCount[simplified rhs]: ", ByteCount@all]; - all - ]; - + all]; (* Take a grid function name and return a name suitable for use in a local computation *) @@ -91,10 +88,8 @@ localNameVectorised[x_] := ToExpression[ToString[x] <> "V"]; invertMap[m_] := Map[#[[2]] -> #[[1]] &, m]; (* A convenient list of the derivative symbols *) -derivativeHeads = {D1, D2, D3, - D11, - D21, D22, - D31, D32, D33, Dplus1, Dplus2, Dplus3}; +derivativeHeads = + {D1, D2, D3, D11, D21, D22, D31, D32, D33, Dplus1, Dplus2, Dplus3}; (* Take an expression containing derivatives (i.e. D1[<stuff> etc) and "hide" all the derivatives and their contents by replacing them @@ -107,7 +102,6 @@ hideDerivatives[x_] := derivatives = Union[Cases[x, _ ? (MemberQ[derivativeHeads, #] &)[__],Infinity]]; unhide = Map[Unique[] -> # &, derivatives]; hide = invertMap[unhide]; - {x /. hide, unhide}]; (* Apply the map (list of rules) to the expression x, but avoid replacing @@ -125,7 +119,6 @@ replaceDerivatives[x_, derivRules_] := replaceCustom = Flatten[derivRules,1]; x /. replaceStandard]; - (* Take a string s and break it into separate lines using l as a guide to the line length. If a word (sequence of non-whitespace characters) is longer than l, do not break it. Add two spaces of @@ -147,50 +140,42 @@ lineBreak[s_, l_] := Return[StringJoin[Riffle[words, " "]]]]; (* Return a CodeGen block which assigns dest by evaluating expr *) -assignVariableFromExpression[dest_, expr_, declare_] := Module[{tSym, type, cleanExpr, code}, - - tSym = Unique[]; - - type = If[StringMatchQ[ToString[dest], "dir*"], "int", "CCTK_REAL"]; - - cleanExpr = ReplacePowers[expr] /. sym`t -> tSym; +assignVariableFromExpression[dest_, expr_, declare_] := + Module[{tSym, type, cleanExpr, code}, + tSym = Unique[]; + type = If[StringMatchQ[ToString[dest], "dir*"], "int", "CCTK_REAL"]; + cleanExpr = ReplacePowers[expr] /. sym`t -> tSym; - If[SOURCELANGUAGE == "C", - code = If[declare, type <> " ", ""] <> ToString[dest] <> " = " <> ToString[cleanExpr, CForm, PageWidth -> Infinity] <> ";\n", - code = ToString@dest <> ".eq." <> ToString[cleanExpr, FortranForm, PageWidth -> 120] <> "\n" - ]; + If[SOURCELANGUAGE == "C", + code = If[declare, type <> " ", ""] <> ToString[dest] <> " = " <> + ToString[cleanExpr, CForm, PageWidth -> Infinity] <> ";\n", + code = ToString@dest <> ".eq." <> ToString[cleanExpr, FortranForm, PageWidth -> 120] + <> "\n"]; - If[SOURCELANGUAGE != "C", - Module[{}, - 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 = lineBreak[code, 70] <> "\n"; - -(* code = StringReplace[code, "Rule" -> " = "]; *) - 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}]; + 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 = lineBreak[code, 70] <> "\n"; + 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}]; (* This flag determines whether you want to generate debugging code to do CCTK_INFO on the variables as they are translated. This can help find problems in the construction of the translator maps. *) - debugInLoop = False; (* Derivative precomputation *) - oldDerivativesUsed[x_] := Union[Cases[x, _ ? (MemberQ[derivativeHeads, #] &)[_],Infinity]]; @@ -207,18 +192,13 @@ precomputeDerivative[d_] := printEq[eq_] := Module[{lhs, rhs, rhsSplit, split, rhsString}, - lhs = First@eq; rhs = Last@eq; - split[ x_ + y___] := { x, " + ..."}; split[-x_ + y___] := {-x, " + ..."}; split[ x_ ] := { x, ""}; - rhsSplit = split[Expand@ReplacePowers@rhs]; - rhsString = ToString@CForm[rhsSplit[[1]]] <> rhsSplit[[2]]; - InfoMessage[InfoFull, " " <> ToString@lhs <> " -> " <> rhsString]]; (* Return the names of any gridfunctions used in the calculation *) @@ -234,14 +214,12 @@ calculationUsedGFsLHS[calc_] := allGFs = allVariables[lookup[calc, Groups]]; Intersection[calcSymbols, allGFs]]; - calculationUsedGFsRHS[calc_] := Module[{calcSymbols, allGFs}, calcSymbols = calculationSymbolsRHS[calc]; allGFs = allVariables[lookup[calc, Groups]]; Intersection[calcSymbols, allGFs]]; - calculationUsedGroups[calc_] := Module[{gfs}, gfs = calculationUsedGFs[calc]; @@ -270,8 +248,7 @@ calculationAllUsedShorthands[calc_] := calculationSymbols[calc_] := Module[{allAtoms}, - allAtoms = Union[Level[lookup[calc, Equations], {-1}] (*, - Level[Map[Last, lookup[calc,PartialDerivatives]], {-1}] *) ]; + allAtoms = Union[Level[lookup[calc, Equations], {-1}]]; Cases[allAtoms, x_Symbol]]; calculationSymbolsLHS[calc_] := @@ -281,7 +258,8 @@ calculationSymbolsLHS[calc_] := calculationSymbolsRHS[calc_] := Module[{allAtoms}, - allAtoms = Union[Map[Last, Flatten@{lookup[calc, Equations], lookup[calc,PartialDerivatives]} ]]; + allAtoms = Union[Map[Last, Flatten@{lookup[calc, Equations], + lookup[calc,PartialDerivatives]} ]]; allAtoms = Union[Level[allAtoms, {-1}]]; Cases[allAtoms, x_Symbol]]; @@ -293,8 +271,6 @@ functionsInCalculation[calc_] := y = Union[Select[x, Context[#] === "Global`" &]]; y]; - - simplifyEquationList[eqs_] := Map[simplifyEquation, eqs]; @@ -304,10 +280,13 @@ simplifyEquation[lhs_ -> rhs_] := VerifyListContent[l_, type_, while_] := Module[{types}, If[!(Head[l] === List), - ThrowError["Expecting a list of ", type, " objects, but found the following, which is not a List object: ", l, while]]; + ThrowError["Expecting a list of ", type, + " objects, but found the following, which is not a List object: ", l, while]]; types = Union[Map[Head, l]]; If [!(types === {type}) && !(types === {}), - ThrowError["Expecting a list of ", type , " objects, but found the following types of object: ", ToString[types], " in ", l, while]]]; + ThrowError["Expecting a list of ", type , + " objects, but found the following types of object: ", + ToString[types], " in ", l, while]]]; VerifyCalculation[calc_] := Module[{calcName}, @@ -326,9 +305,8 @@ VerifyCalculation[calc_] := If[mapContains[calc, Equations], VerifyListContent[First[lookup[calc, Equations]], Rule, " while checking the equation" <> ToString[calcName]], - ThrowError["Invalid Calculation structure. Must contain Equations element: ", ToString[calc], " while checking the calculation called ", ToString[calcName]]]; - ]; - + ThrowError["Invalid Calculation structure. Must contain Equations element: ", + ToString[calc], " while checking the calculation called ", ToString[calcName]]]]; (* Remove equations in the calculation which assign to shorthands which are never used. Do not modify the Shorthands entry. This @@ -342,7 +320,8 @@ VerifyCalculation[calc_] := removed). *) removeUnusedShorthands[calc_] := - Module[{rhsShorthands, lhsShorthands, unusedButAssignedShorthands, removeShorts, eqs, neweqs, newCalc}, + Module[{rhsShorthands, lhsShorthands, unusedButAssignedShorthands, removeShorts, + eqs, neweqs, newCalc}, rhsShorthands = calculationRHSUsedShorthands[calc]; lhsShorthands = calculationLHSUsedShorthands[calc]; unusedButAssignedShorthands = Complement[lhsShorthands, rhsShorthands]; @@ -364,7 +343,6 @@ syncGroup[name_] := {"/* SYNC: "<>name<>" */\n", "/* sync via schedule instead of call CCTK_SyncGroup(i, cctkGH, \"", name, "\") -cut- */\n"} ]; - UncommentSourceSync[setRHS_]:= Module[{cleanRHS}, cleanRHS = Map[SafeStringReplace[#, "/* sync via schedule instead of ", ""]&, setRHS, Infinity]; cleanRHS = Map[SafeStringReplace[#, ") -cut- */", ")"]&, cleanRHS, Infinity]; (* for Fortran *) @@ -376,24 +354,23 @@ UncommentSourceSync[setRHS_]:= Module[{cleanRHS}, InsertSyncFuncName[setRHS_, funcName_]:= Map[SafeStringReplace[#, "/* SYNC: ", "/* SYNC: " <> funcName <> "//"]&, setRHS, Infinity]; -GrepSyncGroups[setRHS_]:= Module[{grepSYNC}, - grepSYNC = Map[PickMatch[#, "*SYNC*"] &, Flatten[setRHS, Infinity]]; - grepSYNC = Select[grepSYNC, IsNotEmptyString]; - - grepSYNC = Map[StringReplace[#, "/* SYNC: " -> ""]&, grepSYNC]; - grepSYNC = Map[StringReplace[#, "*/" -> ""]&, grepSYNC]; +GrepSyncGroups[setRHS_] := + Module[{grepSYNC}, + grepSYNC = Map[PickMatch[#, "*SYNC*"] &, Flatten[setRHS, Infinity]]; + grepSYNC = Select[grepSYNC, IsNotEmptyString]; - grepSYNC -]; + grepSYNC = Map[StringReplace[#, "/* SYNC: " -> ""]&, grepSYNC]; + grepSYNC = Map[StringReplace[#, "*/" -> ""]&, grepSYNC]; -GrepSyncGroups[x_, func_] := Module[{pick}, + grepSYNC]; +GrepSyncGroups[x_, func_] := + Module[{pick}, pick = GrepSyncGroups[x]; pick = Map[PickMatch[#, func <> "//*"] & , pick]; pick = Select[pick, IsNotEmptyString]; - pick = Map[StringReplace[#, func <> "//" -> ""] &, pick] - ]; + pick = Map[StringReplace[#, func <> "//" -> ""] &, pick]]; declareSubblockGFs[sbgfs_, counter_] := If[Length[sbgfs] == 0, @@ -401,7 +378,6 @@ declareSubblockGFs[sbgfs_, counter_] := Join[{DefineVariable[sbgfs[[1]], "CCTK_REAL", "*subblock_gfs[" <> ToString[counter] <> "]"]}, declareSubblockGFs[Rest[sbgfs], counter + 1]]]; - declarePreDefinitions[pDefs_] := CommentedBlock["Declare predefined quantities", Map[DeclareVariable[#, "// CCTK_REAL"] &, Map[First, pDefs]]]; @@ -410,7 +386,6 @@ definePreDefinitions[pDefs_] := CommentedBlock["Initialize predefined quantities", Map[DeclareAssignVariable["CCTK_REAL", #[[1]], #[[2]]] &, pDefs]]; - (* Calculation function generation *) Options[CreateCalculationFunction] = ThornOptions; @@ -456,7 +431,6 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] := InfoMessage[InfoFull, "Groups: ", Map[groupName, groups]]; If[Length@lookupDefault[cleancalc, CollectList, {}] > 0, - eqs = Table[Map[First[#] -> simpCollect[lookup[cleancalc, CollectList], Last[ #], First[#], debug]&, @@ -464,8 +438,7 @@ 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:"]; @@ -475,21 +448,10 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] := functionsPresent = functionsInCalculation[cleancalc]; -(* FIXME: Sascha does not understand the next lines and commented - it out in order to avod problems with using the exp function in BSSN *) -(* allowedFunctions = Map[Head[First[#]] &, pddefs]; *) -(* Print["allowedFunctions == ", allowedFunctions];*) - -(* unknownFunctions = Complement[functionsPresent, allowedFunctions]; - - If[Length[unknownFunctions] != 0, - ThrowError["The following functions are used in the calculation but are not defined:", - unknownFunctions, cleancalc]]; - -*) (* Check that there are no shorthands defined with the same name as a grid function *) If[!(Intersection[shorts, gfs] === {}), - ThrowError["The following shorthands are already declared as grid functions:", Intersection[shorts, gfs]]]; + ThrowError["The following shorthands are already declared as grid functions:", + Intersection[shorts, gfs]]]; (* Check that there are no unknown symbols in the calculation *) allSymbols = calculationSymbols[cleancalc]; @@ -501,8 +463,8 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] := If[unknownSymbols != {}, Module[{}, - ThrowError["Unknown symbols in calculation. Symbols are:", unknownSymbols, "Calculation is:", cleancalc]]]; - + ThrowError["Unknown symbols in calculation. Symbols are:", unknownSymbols, + "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[]", { @@ -510,17 +472,14 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] := "DECLARE_CCTK_PARAMETERS;\n\n", If[!OptionValue[UseLoopControl], DeclareGridLoopVariables[], {}], DeclareFDVariables[], - (* declarePrecomputedDerivatives[dsUsed], *) - (* DeclareDerivatives[pddefs, eqs], *) declareSubblockGFs[subblockGFs, 0], - declarePreDefinitions[pDefs], ConditionalOnParameterTextual["verbose > 1", "CCTK_VInfo(CCTK_THORNSTRING,\"Entering " <> bodyFunctionName <> "\");\n"], - ConditionalOnParameterTextual["cctk_iteration % " <> functionName <> "_calc_every != " <> functionName <> "_calc_offset", - "return;\n"], + ConditionalOnParameterTextual["cctk_iteration % " <> functionName <> "_calc_every != " <> + functionName <> "_calc_offset", "return;\n"], CommentedBlock["Include user-supplied include files", Map[IncludeFile, lookupDefault[cleancalc, DeclarationIncludes, {}]]], @@ -537,25 +496,16 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] := If[gfs != {}, { - (* This has been moved into equationLoop, in order *) - (* to be able to discriminate between the offsets necessary *) - (* depending on whether or not we take derivatives *) - - (* InitialiseGridLoopVariables[], *) - - - (* Have removed ability to include external header files here. - Can be put back when we need it. *) - - eqLoop = Map[equationLoop[#, cleancalc, dsUsed, gfs, shorts, subblockGFs, {}, groups, syncGroups, pddefs, where, addToStencilWidth, useCSE, opts] &, eqs]}, + eqLoop = Map[equationLoop[#, cleancalc, dsUsed, gfs, shorts, subblockGFs, {}, groups, + syncGroups, pddefs, where, addToStencilWidth, useCSE, opts] &, eqs]}, (* search for SYNCs *) If[numeq <= 1, GrepSYNC = GrepSyncGroups[eqLoop], GrepSYNC = {}; eqLoop = UncommentSourceSync[eqLoop]; - InfoMessage[Warning, "> 1 loop in thorn -> scheduling in source code, incompatible with Multipatch!"]; - ]; + InfoMessage[Warning, + "> 1 loop in thorn -> scheduling in source code, incompatible with Multipatch!"]]; InfoMessage[InfoFull, "grepSync from eqLoop: ",GrepSyncGroups[eqLoop]]; @@ -565,21 +515,28 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] := {}] }], Switch[where, - Everywhere, DefineCCTKSubroutine[functionName, "GenericFD_LoopOverEverything(cctkGH, &" <> bodyFunctionName <> ");\n"], - Interior, DefineCCTKSubroutine[functionName, "GenericFD_LoopOverInterior(cctkGH, &" <> bodyFunctionName <> ");\n"], - InteriorNoSync, DefineCCTKSubroutine[functionName, "GenericFD_LoopOverInterior(cctkGH, &" <> bodyFunctionName <> ");\n"], - Boundary, DefineCCTKSubroutine[functionName, "GenericFD_LoopOverBoundary(cctkGH, &" <> bodyFunctionName <> ");\n"], - BoundaryWithGhosts, DefineCCTKSubroutine[functionName, "GenericFD_LoopOverBoundaryWithGhosts(cctkGH, &" <> bodyFunctionName <> ");\n"], - PenaltyPrim2Char, DefineFunction[functionName, "CCTK_INT", + Everywhere, + DefineCCTKSubroutine[functionName, + "GenericFD_LoopOverEverything(cctkGH, &" <> bodyFunctionName <> ");\n"], + Interior, + DefineCCTKSubroutine[functionName, + "GenericFD_LoopOverInterior(cctkGH, &" <> bodyFunctionName <> ");\n"], + InteriorNoSync, + DefineCCTKSubroutine[functionName, + "GenericFD_LoopOverInterior(cctkGH, &" <> bodyFunctionName <> ");\n"], + Boundary, + DefineCCTKSubroutine[functionName, + "GenericFD_LoopOverBoundary(cctkGH, &" <> bodyFunctionName <> ");\n"], + BoundaryWithGhosts, + DefineCCTKSubroutine[functionName, + "GenericFD_LoopOverBoundaryWithGhosts(cctkGH, &" <> bodyFunctionName <> ");\n"], + PenaltyPrim2Char, + DefineFunction[functionName, "CCTK_INT", "CCTK_POINTER_TO_CONST const cctkGH_, CCTK_INT const dir, CCTK_INT const face, CCTK_REAL const * restrict const base_, CCTK_INT const * restrict const lbnd, CCTK_INT const * restrict const lsh, CCTK_INT const * restrict const from, CCTK_INT const * restrict const to, CCTK_INT const rhs_flag, CCTK_INT const num_modes, CCTK_POINTER const * restrict const modes, CCTK_POINTER const * restrict const speeds", - "GenericFD_PenaltyPrim2Char(cctkGH, &" <> bodyFunctionName <> ");\n"], - _, ThrowError["Unknown 'Where' entry in calculation " <> functionName <> ": " <> ToString[where]] - ] - }]; - - - - + "GenericFD_PenaltyPrim2Char(cctkGH, &" <> bodyFunctionName <> ");\n"], + _, + ThrowError["Unknown 'Where' entry in calculation " <> + functionName <> ": " <> ToString[where]]]}]; allVariables[groups_] := Flatten[Map[variablesInGroup, groups], 1]; @@ -607,16 +564,15 @@ checkShorthandAssignmentOrder[eqs_, shorthand_] := (* Make a list of booleans describing, for each equation, whether it uses the given shorthand *) useBooleans = Map[equationUsesShorthand[#, shorthand] &, eqs]; - uses = Position[useBooleans, True]; - lhss = Map[First, eqs]; (* The equation numbers that define this shorthand *) assignments = Position[lhss, shorthand]; If[Length[uses] == 0 && Length[assignments] >= 1, - InfoMessage[Warning, "WARNING: Shorthand ", shorthand, " is defined but not used in this equation list."]]; + InfoMessage[Warning, "WARNING: Shorthand ", shorthand, + " is defined but not used in this equation list."]]; If[Length[uses] == 0, Return[]]; @@ -630,7 +586,8 @@ checkShorthandAssignmentOrder[eqs_, shorthand_] := ThrowError["Shorthand", shorthand, "is not defined in this equation list", eqs]]; If[assignments[[1]] >= firstUse, - ThrowError["Shorthand", shorthand, "is used before it is defined in this equation list", eqs]]]; + ThrowError["Shorthand", shorthand, + "is used before it is defined in this equation list", eqs]]]; defContainsShorthand[def_, shorthands_] := Module[{allAtoms, c}, @@ -638,7 +595,6 @@ Module[{allAtoms, c}, c = Intersection[shorthands, allAtoms]; c != {}]; - (* Split the list of partial derivative definitions into those containing shorthands, and those that do not. *) @@ -658,8 +614,8 @@ pdCanonicalOrdering[name_[inds___] -> x_] := {}]]; (* Given an input list l, return a list L such that L[[i]] is True -only if l[[i]] is the first occurrence of l[[i]] in l and is not in -the list "already" *) + only if l[[i]] is the first occurrence of l[[i]] in l and is not in + the list "already" *) markFirst[l_List, already_List] := If[l =!= {}, {!MemberQ[already, First[l]]} ~Join~ markFirst[Rest[l], already ~Join~ {First[l]}], @@ -681,10 +637,6 @@ equationLoop[eqs_, lhss = Map[#[[1]] &, eqs]; orderings = Flatten[Map[pdCanonicalOrdering, pddefs], 1]; -(* Print["Canonical orderings: ", orderings]; - - Print["Applying orderings to:", eqs];*) - eqsOrdered = (eqs //. orderings); gfsInRHS = Union[Cases[rhss, _ ? (MemberQ[gfs,#] &), Infinity]]; @@ -704,7 +656,8 @@ equationLoop[eqs_, localGFs = Map[localName, gfs]; localMap = Map[# -> localName[#] &, gfs]; - derivSwitch = Join[oldDerivativesUsed[eqsOrdered], GridFunctionDerivativesInExpression[pddefs, eqsOrdered]] != {}; + derivSwitch = Join[oldDerivativesUsed[eqsOrdered], + GridFunctionDerivativesInExpression[pddefs, eqsOrdered]] != {}; gfsDifferentiated = Map[First, GridFunctionDerivativesInExpression[pddefs, eqsOrdered]]; @@ -719,25 +672,14 @@ equationLoop[eqs_, " on the ordering of the loop over grid points."]]; (* Replace the partial derivatives *) - {defsWithoutShorts, defsWithShorts} = splitPDDefsWithShorthands[pddefs, shorts]; -(* Map[ComponentDerivativeOperatorStencilWidth, pddefs];*) - (* This is for the custom derivative operators pddefs *) eqs2 = ReplaceDerivatives[defsWithoutShorts, eqsOrdered, True]; eqs2 = ReplaceDerivatives[defsWithShorts, eqs2, False]; checkEquationAssignmentOrder[eqs2, shorts]; - code = {(*InitialiseGridLoopVariables[derivSwitch, addToStencilWidth], *) - functionName = ToString@lookup[cleancalc, Name]; - -(* - calcCode = - Map[{assignVariableFromExpression[#[[1]], #[[2]]], "\n"} &, - replaceDerivatives[replaceWithDerivativesHidden[eqs2, localMap], {}] - ]; -*) + functionName = ToString@lookup[cleancalc, Name]; (* Replace grid functions with their local forms, and replace partial dervatives with their precomputed values *) @@ -758,25 +700,13 @@ equationLoop[eqs_, MapThread[{assignVariableFromExpression[#1[[1]], #1[[2]], #2], "\n"} &, {eqsReplaced, declare}]; + code = { Join[ - (* - CommentedBlock["Declare temporary variables for vectorisation", - DeclareVariablesInLoopVectorised[ - gfsInLHS, - Map[localNameVectorised, gfsInLHS], - Map[localName, gfsInLHS]]], - *) GenericGridLoop[functionName, {declarePrecomputedDerivatives[dsUsed], - (* DeclareDerivatives[pddefs, eqs], *) DeclareDerivatives[defsWithoutShorts, eqsOrdered], -(* - CommentedBlock["Assign local copies of grid functions", - Map[DeclareAssignVariableInLoop["CCTK_REAL", localName[#], GridName[#]] &, - gfsInRHS]], -*) CommentedBlock["Assign local copies of grid functions", Map[DeclareMaybeAssignVariableInLoop[ "CCTK_REAL", localName[#], GridName[#], @@ -784,26 +714,10 @@ equationLoop[eqs_, "*stress_energy_state"] &, gfsInRHS]], -(* - CommentedBlock["Check for nans", - Join[{"if (\n"}, - Map[{" isnan (", localName[#], ") ||\n"} &, - gfsInRHS], - {" 0)\n", - "{\n", - " CCTK_VInfo(CCTK_THORNSTRING, \"NaN found in RHS:\");\n", - " 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, \"from: %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"}, - Map[{" CCTK_VInfo(CCTK_THORNSTRING, \"", localName[#], ": %.17g\", (double)", localName[#], ");\n"} &, - gfsInRHS], - {"}\n"}], -*) - CommentedBlock["Assign local copies of subblock grid functions", - Map[DeclareAssignVariableInLoop["CCTK_REAL", localName[#], SubblockGridName[#]] &, - subblockGFs]], + Map[DeclareAssignVariableInLoop["CCTK_REAL", localName[#], + SubblockGridName[#]] &, + subblockGFs]], CommentedBlock["Include user supplied include files", Map[IncludeFile, incs]], @@ -821,47 +735,17 @@ equationLoop[eqs_, CommentedBlock["Copy local copies back to grid functions", Map[AssignVariableInLoop[GridName[#], localName[#]] &, gfsInLHS]], - (* - CommentedBlock["Copy local copies back to grid functions", - AssignVariablesInLoopVectorised[ - gfsInLHS, - Map[localNameVectorised, gfsInLHS], - Map[localName, gfsInLHS]]], - *) CommentedBlock["Copy local copies back to subblock grid functions", Map[AssignVariableInLoop[SubblockGridName[#], localName[#]] &, subblockGFs]], -(* - CommentedBlock["Check for nans", - Join[{"if (\n"}, - Map[{" isnan (", localName[#], ") ||\n"} &, - gfsInLHS], - {" 0)\n", - "{\n", - " CCTK_VInfo(CCTK_THORNSTRING, \"NaN found in LHS:\");\n", - " 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, \"from: %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"}, - Map[{" CCTK_VInfo(CCTK_THORNSTRING, \"", localName[#], ": %.17g\", (double)", localName[#], ");\n"} &, - gfsInLHS], - {"}\n"}], -*) - If[debugInLoop, Map[InfoVariable[GridName[#]] &, gfsInLHS], ""]}, opts]] }; lhsGroupNames = containingGroups[gfsInLHS, groups]; actualSyncGroups = Intersection[lhsGroupNames, syncGroups]; - (* This is nonsense. You need to synchronize only if the NEXT - loops contain derivatives. This cannot be determined here, so I am - changing the code so that synchronization is always performed *) - -(* If[Not@derivSwitch, actualSyncGroups = {}]; only sync when derivs are taken *) - If[Length@actualSyncGroups > 0, InfoMessage[InfoFull, "Synchronizing groups: ", actualSyncGroups]; syncCode = Map[syncGroup, actualSyncGroups]; |