Add special relativistic Euler code and parameter file

2 files changed, 198 insertions, 4 deletions

diff --git a/Examples/EulerSR.m b/Examples/EulerSR.m
new file mode 100644
index 0000000..2c58c96
--- /dev/null
+++ b/Examples/EulerSR.m
@@ -0,0 +1,165 @@
+
+Get["KrancThorn`"];
+
+SetEnhancedTimes[False];
+
+(**************************************************************************************)
+(* Tensors *)
+(**************************************************************************************)
+
+(* Register the tensor quantities with the TensorTools package *)
+Map[DefineTensor, {Den, S, tau, rho, v, epsi, W, h, p}];
+
+(**************************************************************************************)
+(* Groups *)
+(**************************************************************************************)
+
+evolvedGroups = Map[CreateGroupFromTensor, {Den, S[lj], tau}];
+nonevolvedGroups = Map[CreateGroupFromTensor, 
+{
+  rho, v[uj], epsi, W, h, p
+}];
+
+declaredGroups = Join[evolvedGroups, nonevolvedGroups];
+declaredGroupNames = Map[First, declaredGroups];
+
+groups = declaredGroups;
+
+(**************************************************************************************)
+(* Initial data *)
+(**************************************************************************************)
+
+initialShockCalc =
+{
+  Name -> "eulersr_initial_shock",
+  Schedule -> {"at CCTK_INITIAL as eulersr_initial"},
+  ConditionalOnKeyword -> {"initial_data", "shock"},
+  Equations ->
+  {
+    rho -> rhoR0 StepFunction[x-0.5] + rhoL0 (1-StepFunction[x-0.5]),
+    v[1] -> vR0 StepFunction[x-0.5] + vL0 (1-StepFunction[x-0.5]),
+    v[2] -> 0,
+    v[3] -> 0,
+    epsi -> epsiR0 StepFunction[x-0.5] + epsiL0 (1-StepFunction[x-0.5])
+  }
+};
+
+(**************************************************************************************)
+(* Evolution equations *)
+(**************************************************************************************)
+
+(* Euler's equation is dot[u] + PD[F[ui],li] = 0
+
+   with
+
+     u = {D, S, tau}
+
+   and
+
+     DF[ui] = D v[ui]
+     SF[ui,lj] = S[lj] v[ui] + p Euc[lj,ui]
+     tauF[ui] = v[ui](tau + p)
+
+*)
+
+eulerCons =
+{
+  Name -> "eulersr_cons_calc",
+  Shorthands -> {pBar, Z, Ssq, vsq, pEOS, f, cs, df, Wx},
+
+  Primitives -> {rho, v[ui], epsi},
+
+  Equations ->
+  {
+    flux[Den,ui] -> Den v[ui],
+    flux[S[lj],ui] -> S[lj] v[ui] + ((gamma-1) rho epsi (* This term is p *)) Euc[ui,lj],
+    flux[tau,ui] -> v[ui](tau + (gamma-1) rho epsi (* This term is p *))
+  },
+
+  ConservedEquations ->
+  {
+    Wx -> 1 - v[ui] v[uj] Euc[li,lj],
+    W -> Wx^(-1/2),
+    p -> (gamma-1) rho epsi,
+    h -> 1 + epsi + p/rho,
+
+    Den -> rho W,
+    S[li] -> rho h W^2 v[uj] Euc[li,lj],
+    tau -> rho h W^2 - p - Den
+  },
+
+  PrimitiveEquations ->
+  {
+    (* To compute p, given Den, S[ui], tau and a guess for p (pBar),
+         Z = tau + Den + pBar
+         S2 = S[li] S[lj] Euc[ui,uj]
+         v2 = S2/Z^2
+         W = (1-v2)^(-1/2)
+         rho = Den/W
+         h = Z/(rho W^2)
+         epsi = h-1-pBar/rho
+         pNew = (gamma - 1) rho epsi
+         f = pNew - pBar
+         cs = Sqrt[gamma (gamma-1) epsi/h]
+         df = v2 cs^2 - 1
+         
+         -> p (until f is sufficiently small) (also get rho, epsi)
+    *)
+
+    pBar -> p, (* from previous timestep *)
+    (* Start loop *)
+
+    Sequence@@Join@@Table[
+    {Z -> tau + Den + pBar,
+    Ssq -> S[li] S[lj] Euc[ui,uj],
+    vsq -> Ssq/Z^2,
+    W -> (1-vsq)^(-1/2),
+    rho -> Den/W,
+    h -> Z/(rho W^2),
+    epsi -> h-1-pBar/rho,
+    pEOS -> (gamma - 1) rho epsi,
+    f -> pEOS - pBar,
+    cs -> Sqrt[gamma (gamma-1) epsi/h],
+    df -> vsq cs^2 - 1,
+    pBar -> pBar - f/df},
+      {i, 1, 5}],
+
+    (* end of loop *)
+
+    p -> pBar,
+
+    v[ui] -> S[lj] Euc[ui,uj] / (rho h W^2)
+  }
+}
+
+(**************************************************************************************)
+(* Parameters *)
+(**************************************************************************************)
+
+realParameters = {sigma, v0, amp, rhoR0, rhoL0, vR0, vL0, epsiR0, epsiL0, gamma};
+
+keywordParameters = {
+  {
+    Name -> "initial_data",
+    Default -> "shock",
+    AllowedValues -> {"shock"}
+  }
+};
+
+(**************************************************************************************)
+(* Construct the thorn *)
+(**************************************************************************************)
+
+calculations = 
+{
+  initialShockCalc
+};
+
+consCalculations = {eulerCons};
+
+CreateKrancThornTT[groups, ".", "EulerSR", 
+  Calculations -> calculations,
+  ConservationCalculations -> consCalculations,
+  DeclaredGroups -> declaredGroupNames,
+  RealParameters -> realParameters,
+  KeywordParameters -> keywordParameters];
diff --git a/Examples/euler_shock.par b/Examples/euler_shock.par
index 4178f98..f4d873e 100644
--- a/Examples/euler_shock.par
+++ b/Examples/euler_shock.par
@@ -40,7 +40,7 @@ CoordBase::xmax                         = 2
 CoordBase::ymax                         = 0.1
 CoordBase::zmax                         = 0.1
 
-CoordBase::dx                           = 0.04
+CoordBase::dx                           = 0.0025
 CoordBase::dy                           = 0.1
 CoordBase::dz                           = 0.1
 
@@ -80,6 +80,7 @@ Carpet::num_integrator_substeps         = 4
 
 Cactus::terminate                     = "time"
 Cactus::cctk_final_time               = 0.2
+Cactus::cctk_itlast               = 2
 
 Time::dtfac                           = 0.25
 MethodOfLines::ode_method             = "RK4"
@@ -124,7 +125,7 @@ Carpet::check_for_poison                = no
 #############################################################
 
 IO::out_dir                   = $parfile
-IO::out_fileinfo              = "all"
+IO::out_fileinfo              = "none"
 
 CarpetIOBasic::outInfo_every        = 1
 CarpetIOBasic::outInfo_vars         = "Euler::rho Euler::DenF Euler::v1 Euler::p Euler::Den Euler::S1 Euler::En"
@@ -133,9 +134,37 @@ CarpetIOScalar::outScalar_every = 1
 CarpetIOScalar::outScalar_vars = ""
 CarpetIOScalar::outScalar_reductions = "minimum maximum norm2"
 
-IOASCII::out1D_every          = 1
+IOASCII::out1D_every          = 16
 IOASCII::out1D_x              = "yes"
-IOASCII::out1D_vars           = "Euler::rho Euler::v1 Euler::p Euler::Den Euler::S1 Euler::En"
+IOASCII::out1D_y              = "no"
+IOASCII::out1D_z              = "no"
+IOASCII::out1D_d              = "no"
+IOASCII::out1D_vars           = "
+Euler::Den
+Euler::S1
+Euler::En
+Euler::rho
+Euler::v1
+Euler::p
+
+Euler::DenLeft
+Euler::SLeft1
+Euler::EnLeft
+Euler::rhoLeft
+Euler::vLeft1
+Euler::pLeft
+
+Euler::DenRight
+Euler::SRight1
+Euler::EnRight
+Euler::rhoRight
+Euler::vRight1
+Euler::pRight
+
+Euler::DenF
+Euler::SF1
+Euler::EnF
+"
 
 CarpetIOASCII::out_precision = 19
 CarpetIOASCII::out3D_ghosts   = "yes"