diff options
| -rw-r--r-- | m/Makefile | 8 | ||||
| -rw-r--r-- | m/grhydro.m | 154 | ||||
| -rw-r--r-- | m/hydro.m | 241 |
3 files changed, 402 insertions, 1 deletions
@@ -1,6 +1,6 @@ # -*-Makefile-*- -all: McLachlan_ADM.out McLachlan_BSSN.out WaveToy.out +all: McLachlan_ADM.out McLachlan_BSSN.out WaveToy.out hydro.out @echo @echo "The Cactus thorns are up to date." @echo @@ -14,11 +14,17 @@ McLachlan_BSSN.out: McLachlan_BSSN.m rm -rf ML_BSSN* ./runmath.sh $^ for thorn in ML_BSSN*; do ./copy-if-changed.sh $$thorn ../$$thorn; done + WaveToy.out: WaveToy.m rm -rf ML_WaveToy ML_FOWaveToy ./runmath.sh $^ for thorn in ML_WaveToy ML_FOWaveToy; do ./copy-if-changed.sh $$thorn ../$$thorn; done +hydro.out: hydro.m + rm -rf ML_hydro + ./runmath.sh $^ + for thorn in ML_hydro; do ./copy-if-changed.sh $$thorn ../$$thorn; done + clean: rm -rf McLachlan_ADM.out McLachlan_ADM.err ML_ADM* rm -rf McLachlan_BSSN.out McLachlan_BSSN.err ML_BSSN* diff --git a/m/grhydro.m b/m/grhydro.m new file mode 100644 index 0000000..a71a28b --- /dev/null +++ b/m/grhydro.m @@ -0,0 +1,154 @@ +Get["KrancThorn`"]; + +(*SetDebugLevel[InfoFull];*) + +SetEnhancedTimes[False]; +SetSourceLanguage["C"]; + +(**************************************************************************** + 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], + + PDplus[i_] -> DPlus[i], + PDminus[i_] -> DMinus[i], + PDplus[i_,j_] -> DPlus[i] DPlus[j], + + PDonesided2nd[1] -> dir[1] (-shift[1]^(2 dir[1]) + 4 shift[1]^dir[1] - 3 )/ + (2 spacing[1]), + PDonesided2nd[2] -> dir[2] (-shift[2]^(2 dir[2]) + 4 shift[2]^dir[2] - 3 )/ + (2 spacing[2]), + PDonesided2nd[3] -> dir[3] (-shift[3]^(2 dir[3]) + 4 shift[3]^dir[3] - 3 )/ + (2 spacing[3]) +}; + +(*PD = PDstandard2nd;*) + +(**************************************************************************** + Tensors +****************************************************************************) + +(* Register all the tensors that will be used with TensorTools *) +Map[DefineTensor, +{ + rho, vel, eps, press, + dens, mom, tau +}]; + +(* Register the TensorTools symmetries (this is very simplistic) *) +Map[AssertSymmetricDecreasing, +{ +}]; + +(* Determinants of the metrics in terms of their components + (Mathematica symbolic expressions) *) +gDet = Det[MatrixOfComponents[g[la,lb]]]; + +(**************************************************************************** + Groups +****************************************************************************) + +SetGroupTimelevels[g_,tl_] = Join[g, {Timelevels -> tl}]; + +evolvedGroups = + {CreateGroupFromTensor [dens ], (* mass density *) + CreateGroupFromTensor [mom[la]], (* momentum density *) + CreateGroupFromTensor [tau ]}; (* energy density *) +evaluatedGroups = + {CreateGroupFromTensor [rho ], (* mass density *) + CreateGroupFromTensor [vel[ua]], (* velocity *) + CreateGroupFromTensor [eps ], (* specific internal energy *) + CreateGroupFromTensor [press ]}; (* pressure *) + +declaredGroups = Join [evolvedGroups, evaluatedGroups]; +declaredGroupNames = Map [First, declaredGroups]; + +groups = Join[declaredGroups]; + +(******************************************************************************) +(* Initial data *) +(******************************************************************************) + +initialCalc = +{ + Name -> BSSN <> "_vacuum", + Schedule -> {"IN ADMBase_InitialData"}, + ConditionalOnKeyword -> {"my_initial_data", "vacuum"}, + Equations -> + { + rho -> 0, + vel[ua] -> 0, + eps -> 0 + } +}; + +(******************************************************************************) +(* Convert from primitive to conserved variables *) +(******************************************************************************) + +prim2conCalc = +{ + Name -> BSSN <> "_con2prim", + Schedule -> {"AT initial AFTER ADMBase_PostInitial"}, + Equations -> + { + dens -> rho, + mom[la] -> rho vel[ua], + tau -> (1/2) rho vel[ua] vel[la] + rho eps + } +}; + +(******************************************************************************) +(* Convert from conserved to primitive variables *) +(******************************************************************************) + +con2primCalc = +{ + Name -> BSSN <> "_con2prim", + Schedule -> {"IN " <> BSSN <>"_con2primGroup"}, + Equations -> + { + rho -> dens, + vel[ua] -> mom[la] / dens, + eps -> tau / dens - (1/2) vel[ua] vel[la], + + press -> Gamma rho eps + + alpha PD[vel[ua],la] + } +}; + +(******************************************************************************) +(* Evolution equations *) +(******************************************************************************) + +evolCalc = +{ + Name -> BSSN <> "_RHS", + Schedule -> {"IN " <> BSSN <>"_evolCalcGroup"}, + Shorthands -> {rhov[ua], momv[la,ub], tauv[ua]}, + Equations -> + { + (* dt rho + div rho v = 0 *) + rhov[ua] -> mom[la], + dot[dens] -> - PD[rhov[ua],la], + + (* dt pi + div (pi v + P) = 0 *) + momv[la,ub] -> mom[la] vel[ub] + KD[la,ub] press, + dot[mom[la]] -> - PD[momv[la,ub],lb], + + (* dt tau + div (tauv + P) = 0 *) + tauv[ua] -> tau vel[ua] + press vel[ua], + dot[tau] -> - PD[tauv[ua],la] + } +}; diff --git a/m/hydro.m b/m/hydro.m new file mode 100644 index 0000000..1b90f22 --- /dev/null +++ b/m/hydro.m @@ -0,0 +1,241 @@ +$Path = Join[$Path, {"~/Calpha/kranc/Tools/CodeGen", + "~/Calpha/kranc/Tools/MathematicaMisc"}]; + +Get["KrancThorn`"]; + +(*SetDebugLevel[InfoFull];*) + +SetEnhancedTimes[False]; +SetSourceLanguage["C"]; + +(**************************************************************************** + Derivatives +****************************************************************************) + +derivOrder = 2; + +KD = KroneckerDelta; + +derivatives = +{ + PDstandardNth[i_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i], + PDstandardNth[i_, i_] -> StandardCenteredDifferenceOperator[2,derivOrder/2,i], + PDstandardNth[i_, j_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i] + StandardCenteredDifferenceOperator[1,derivOrder/2,j] +}; + +FD = PDstandardNth; +ResetJacobians; +DefineJacobian[PD, FD, KD, Zero3]; + +(**************************************************************************** + Tensors +****************************************************************************) + +(* Register all the tensors that will be used with TensorTools *) +Map[DefineTensor, +{ + rho, vel, eps, press, vol, + mass, mom, ene, + massflux, momflux, eneflux +}]; + +(* Determinants of the metrics in terms of their components + (Mathematica symbolic expressions) *) +gDet = Det[MatrixOfComponents[g[la,lb]]]; + +(**************************************************************************** + Groups +****************************************************************************) + +SetGroupTimelevels[g_,tl_] = Join[g, {Timelevels -> tl}]; + +evolvedGroups = { + CreateGroupFromTensor [mass ], (* mass *) + CreateGroupFromTensor [mom[la]], (* momentum *) + CreateGroupFromTensor [ene ] (* total energy *) +}; +evaluatedGroups = { + CreateGroupFromTensor [rho ], (* mass density *) + CreateGroupFromTensor [vel[ua]], (* velocity *) + CreateGroupFromTensor [eps ], (* specific internal energy *) + CreateGroupFromTensor [press ], (* pressure *) + CreateGroupFromTensor [vol ], (* volume *) + CreateGroupFromTensor [massflux[ua]], (* mass flux *) + CreateGroupFromTensor [momflux[la,ub]], (* momentum flux *) + CreateGroupFromTensor [eneflux[ua]] (* energy flux *) +}; + +declaredGroups = Join [evolvedGroups, evaluatedGroups]; +declaredGroupNames = Map [First, declaredGroups]; + +groups = Join[declaredGroups]; + +(******************************************************************************) +(* Initial data *) +(******************************************************************************) + +vacuumCalc = +{ + Name -> "hydro_vacuum", + Schedule -> {"IN ADMBase_InitialData"}, + ConditionalOnKeyword -> {"initial_data", "vacuum"}, + Equations -> + { + rho -> 0, + vel[ua] -> 0, + eps -> 0 + } +}; + +soundWaveCalc = +{ + Name -> "hydro_soundWave", + Schedule -> {"IN ADMBase_InitialData"}, + ConditionalOnKeyword -> {"initial_data", "sound wave"}, + Equations -> + { + rho -> 1.0, + vel[ua] -> A Sin[2 Pi x / L], + eps -> 1.0 + } +}; + +(******************************************************************************) +(* Convert from primitive to conserved variables *) +(******************************************************************************) + +prim2conCalc = +{ + Name -> "hydro_prim2con", + Schedule -> {"AT initial AFTER ADMBase_PostInitial"}, + Equations -> + { + vol -> h^3, + mass -> vol rho, + mom[la] -> mass vel[la], + ene -> (1/2) mass vel[ua] vel[la] + mass eps + } +}; + +(******************************************************************************) +(* Convert from conserved to primitive variables *) +(******************************************************************************) + +con2primCalc = +{ + Name -> "hydro_con2prim", + Schedule -> {"IN hydro_con2primGroup"}, + Equations -> + { + rho -> mass / vol, + vel[ua] -> mom[ua] / mass, + eps -> ene / mass - (1/2) vel[ua] vel[la], + + press -> Gamma rho eps + + alpha PD[vel[ua],la] + } +}; + +(******************************************************************************) +(* Evolution equations *) +(******************************************************************************) + +fluxCalc = +{ + Name -> "hydro_fluxes", + Schedule -> {"IN hydro_evolCalcGroup"}, + Equations -> + { + massflux[ua] -> mass vel[ua], + momflux [la,ub] -> mom[la] vel[ub] + KD[la,ub] press / vol, + eneflux [ua] -> (ene + press / vol) vel[ua] + } +}; + +evolCalc = +{ + Name -> "hydro_RHS", + Schedule -> {"IN hydro_evolCalcGroup AFTER hydro_fluxes"}, + Where -> Interior, + Equations -> + { + (* dt rho + div rho v = 0 *) + (* dt m + div m v = 0 *) + dot[mass] -> - PD[massflux[ua],la], + + (* dt pi + div (pi v + P) = 0 *) + (* dt p + div (p v + P/V) = 0 *) + dot[mom[la]] -> - PD[momflux[la,ub],lb], + + (* dt tau + div (tau v + P) = 0 *) + (* dt t + div (t v + P/V) = 0 *) + dot[ene] -> - PD[eneflux[ua],la] + } +}; + +(******************************************************************************) +(* Parameters *) +(******************************************************************************) + +keywordParameters = +{ + { + Name -> "initial_data", + (* Visibility -> "restricted", *) + (* Description -> "ddd", *) + AllowedValues -> {"vacuum", "sound wave"}, + Default -> "vacuum" + } +}; + +realParameters = +{ + { + Name -> h, + Description -> "grid spacing", + Default -> 0.01 + }, + { + Name -> A, + Description -> "sound wave amplitude", + Default -> 0.001 + }, + { + Name -> L, + Description -> "sound wave wavelength", + Default -> 1 + }, + { + Name -> alpha, + Description -> "artificial viscosity coefficient", + Default -> 0 + }, + { + Name -> Gamma, + Description -> "polytropic exponent", + Default -> 4/3 + } +}; + +(******************************************************************************) +(* Construct the thorns *) +(******************************************************************************) + +calculations = +{ + vacuumCalc, + soundWaveCalc, + prim2conCalc, + con2primCalc, + evolCalc +}; + +CreateKrancThornTT [groups, ".", "ML_hydro", + Calculations -> calculations, + DeclaredGroups -> declaredGroupNames, + PartialDerivatives -> derivatives, + UseLoopControl -> True, + KeywordParameters -> keywordParameters, + RealParameters -> realParameters +]; |