Merge pull request #78 from MChabanov/entropy

Entropy

AsterSeeds/param.ccl

@@ -2,6 +2,7 @@
 
 KEYWORD test_type "Type of test to set up" STEERABLE=never
 {
+  "None"   :: "Initial data should be written by another thorn"
   "1DTest" :: ""
   "2DTest" :: ""
   "3DTest" :: "" 
@@ -149,7 +150,14 @@ REAL dipole_z[2] "z-coordinate of the dipole center" STEERABLE=ALWAYS
   *:*     :: "Anything"
 } 0.0
 
+# entropy
+
+BOOLEAN set_entropy_postinitial "Set the entropy consistently in HydroBaseX_PostInitial" STEERABLE=always
+{
+} "yes"
+
 SHARES: EOSX
+USES KEYWORD evolution_eos
 USES CCTK_REAL poly_gamma
 USES CCTK_REAL poly_k
 USES CCTK_REAL gl_gamma

AsterSeeds/schedule.ccl

@@ -91,3 +91,14 @@ if (CCTK_Equals(test_type, "3DTest")) {
   }
 
 }
+
+#Initial conditions for entropy
+
+schedule SetEntropy IN HydroBaseX_PostInitial
+{
+  LANG: C
+
+  READS: HydroBaseX::rho(everywhere) HydroBaseX::eps(everywhere)
+  WRITES: HydroBaseX::entropy(everywhere)
+
+} "Set initial entropy"

AsterSeeds/src/SetEntropy.cxx

@@ -0,0 +1,40 @@
+#include <loop.hxx>
+#include <loop_device.hxx>
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include <cstdio>
+#include <cstdbool>
+#include <cmath> 
+
+#include "eos.hxx"
+#include "eos_idealgas.hxx"
+
+extern "C" void SetEntropy(CCTK_ARGUMENTS)
+{
+
+  using namespace EOSX;
+
+  DECLARE_CCTK_ARGUMENTSX_SetEntropy;
+  DECLARE_CCTK_PARAMETERS;
+
+  if (set_entropy_postinitial) {
+
+    eos::range rgeps(eps_min, eps_max), rgrho(rho_min, rho_max),
+         rgye(ye_min, ye_max);
+
+    const eos_idealgas eos_th(gl_gamma, particle_mass, rgeps, rgrho, rgye);
+
+    grid.loop_all_device<1, 1, 1>(
+          grid.nghostzones,
+          [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+
+            entropy(p.I) = eos_th.kappa_from_valid_rho_eps_ye(rho(p.I),eps(p.I),1.0);
+
+    });
+  }
+}
+
+

AsterSeeds/src/make.code.defn

@@ -7,8 +7,8 @@ SRCS =	\
   3D_tests.cxx \
   atmosphere.cxx \
   magTOV.cxx \
-  magBNS.cxx  
-	
+  magBNS.cxx \
+  SetEntropy.cxx
 
 # Subdirectories containing source files
 SUBDIRS = 

AsterX/configuration.ccl

@@ -1,3 +1,7 @@
 # Configuration definitions for thorn AsterX
 
 REQUIRES Loop EOSX Con2PrimFactory ReconX AsterUtils
+
+PROVIDES AsterX
+{
+}

AsterX/interface.ccl

@@ -8,7 +8,7 @@ USES INCLUDE HEADER: fixmath.hxx
 USES INCLUDE HEADER: loop_device.hxx
 USES INCLUDE HEADER: eos_1p.hxx eos_polytropic.hxx
 USES INCLUDE HEADER: eos.hxx eos_idealgas.hxx
-USES INCLUDE HEADER: c2p.hxx c2p_2DNoble.hxx c2p_1DPalenzuela.hxx
+USES INCLUDE HEADER: c2p.hxx c2p_2DNoble.hxx c2p_1DPalenzuela.hxx c2p_1DEntropy.hxx
 USES INCLUDE HEADER: reconstruct.hxx
 USES INCLUDE HEADER: aster_fd.hxx aster_interp.hxx aster_utils.hxx
 
@@ -23,11 +23,24 @@ CCTK_REAL mom TYPE=gf CENTERING={ccc} TAGS='rhs="momrhs" dependents="HydroBaseX:
 
 CCTK_REAL tau TYPE=gf CENTERING={ccc} TAGS='rhs="taurhs" dependents="HydroBaseX::eps HydroBaseX::press HydroBaseX::temperature HydroBaseX::Ye HydroBaseX::Bvec TmunuBaseX::eTij"' "Conserved internal energy density"
 
+# DEnt is the conserved "evolved" entropy.
+# The corresponding primitive is stored in HydroBaseX::entropy.
+# However, note that this gf stores the "evolved" entropy which is 
+# not necessarily the "physical" entropy. This depends on the EOS, e.g.
+# for the ideal gas we have entropy = p rho^(-gamma).
+# The distinction between "evolved" and "physical" entropy is made
+# explicit in EOSX where functions with "entropy_..." and "kappa_..." 
+# refer to "physical" and "evolved" entropy, respectively. In this 
+# evolution thorn we always refer to DEnt and entropy to describe 
+# the "evolved" quantities. 
+CCTK_REAL DEnt TYPE=gf CENTERING={ccc} TAGS='rhs="DEntrhs" dependents="HydroBaseX::entropy AsterX::flux_x AsterX::flux_y AsterX::flux_z"' "Advected entropy density"
+
 CCTK_REAL flux_x TYPE=gf CENTERING={vcc} TAGS='checkpoint="no"'
 {
   fxdens
   fxmomx fxmomy fxmomz
   fxtau
+  fxDEnt
   fxBx fxBy fxBz
 } "Fluxes in x direction"
 
@@ -36,6 +49,7 @@ CCTK_REAL flux_y TYPE=gf CENTERING={cvc} TAGS='checkpoint="no"'
   fydens
   fymomx fymomy fymomz
   fytau
+  fyDEnt
   fyBx fyBy fyBz
 } "Fluxes in y direction"
 
@@ -44,6 +58,7 @@ CCTK_REAL flux_z TYPE=gf CENTERING={ccv} TAGS='checkpoint="no"'
   fzdens
   fzmomx fzmomy fzmomz
   fztau
+  fzDEnt
   fzBx fzBy fzBz
 } "Fluxes in z direction"
 
@@ -56,6 +71,7 @@ CCTK_REAL momrhs TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"'
 
 CCTK_REAL taurhs TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' "Conserved internal energy density RHS" 
 
+CCTK_REAL DEntrhs TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' "Advected entropy density RHS"
 
 CCTK_REAL Avec_x TYPE=gf CENTERING={cvv} TAGS='rhs="Avec_x_rhs"' "x-component of vector potential"
 CCTK_REAL Avec_y TYPE=gf CENTERING={vcv} TAGS='rhs="Avec_y_rhs"' "y-component of vector potential"

AsterX/par/Scaling_magTOV_Z4c_unigrid.par

@@ -110,7 +110,6 @@ Con2PrimFactory::unit_test = "yes"
 Con2PrimFactory::B_lim = 1e8
 Con2PrimFactory::vw_lim = 1e8
 Con2PrimFactory::Ye_lenient = "yes"
-Con2PrimFactory::rho_strict = 6.4e-05
 
 EOSX::evolution_eos = "IdealGas"
 EOSX::gl_gamma = 2.0

AsterX/par/magTOV_Cowling_unigrid.par

@@ -110,7 +110,6 @@ Con2PrimFactory::unit_test = "yes"
 Con2PrimFactory::B_lim = 1e8
 Con2PrimFactory::vw_lim = 1e8
 Con2PrimFactory::Ye_lenient = "yes"
-Con2PrimFactory::rho_strict = 6.4e-05
 
 EOSX::evolution_eos = "IdealGas"
 EOSX::gl_gamma = 2.0

AsterX/par/magTOV_Z4c_AMR.par

@@ -121,7 +121,6 @@ Con2PrimFactory::unit_test = "yes"
 Con2PrimFactory::B_lim = 1e8
 Con2PrimFactory::vw_lim = 1e8
 Con2PrimFactory::Ye_lenient = "yes"
-Con2PrimFactory::rho_strict = 6.4e-05
 
 EOSX::evolution_eos = "IdealGas"
 EOSX::gl_gamma = 2.0

AsterX/param.ccl

@@ -15,6 +15,10 @@ BOOLEAN use_uct "Shall we use the Upwind-CT method to compute the electric field
 {
 } "no"
 
+BOOLEAN use_entropy_fix "Shall we use inversion based on entropy as a backup?" STEERABLE=always
+{
+} "no"
+
 KEYWORD flux_type "Flux solver" STEERABLE=always
 {
   "LxF" :: ""
@@ -108,7 +112,6 @@ USES CCTK_REAL p_atmo
 USES CCTK_REAL Ye_atmo
 USES CCTK_REAL B_lim
 USES CCTK_REAL vw_lim
-USES CCTK_REAL rho_strict
 USES BOOLEAN Ye_lenient
 USES CCTK_INT max_iter
 USES CCTK_REAL c2p_tol
@@ -120,6 +123,8 @@ USES CCTK_REAL alp_thresh
 USES CCTK_REAL rho_BH
 USES CCTK_REAL eps_BH
 USES CCTK_REAL vwlim_BH
+USES CCTK_REAL cons_error_limit
+USES BOOLEAN use_z
 
 SHARES: ReconX
 USES KEYWORD reconstruction_method

AsterX/schedule.ccl

@@ -1,8 +1,8 @@
 # Schedule definitions for thorn AsterX
 
-STORAGE: dens mom tau dB Psi HydroBaseX::Bvec dBx_stag dBy_stag dBz_stag
+STORAGE: dens mom tau DEnt dB Psi HydroBaseX::Bvec dBx_stag dBy_stag dBz_stag
 STORAGE: flux_x flux_y flux_z
-STORAGE: densrhs momrhs taurhs Avec_x_rhs Avec_y_rhs Avec_z_rhs Psi_rhs
+STORAGE: densrhs momrhs taurhs DEntrhs Avec_x_rhs Avec_y_rhs Avec_z_rhs Psi_rhs
 STORAGE: ADMBaseX::metric ADMBaseX::lapse ADMBaseX::shift ADMBaseX::curv
 STORAGE: Aux_in_RHSof_A_Psi
 STORAGE: TmunuBaseX::eTtt TmunuBaseX::eTti TmunuBaseX::eTij
@@ -69,12 +69,13 @@ SCHEDULE AsterX_Prim2Con_Initial IN AsterX_InitialGroup AFTER AsterX_ComputeBFro
   LANG: C
   READS: ADMBaseX::metric(interior)
   READS: HydroBaseX::rho(interior) HydroBaseX::vel(interior) HydroBaseX::eps(interior) HydroBaseX::press(interior) HydroBaseX::Bvec(interior)
-  WRITES: dens(interior) tau(interior) mom(interior) dB(interior)
+  READS: HydroBaseX::entropy(interior)
+  WRITES: dens(interior) tau(interior) DEnt(interior) mom(interior) dB(interior)
   WRITES: Psi(everywhere)
   WRITES: saved_prims
   WRITES: zvec
   WRITES: svec
-  SYNC: dens tau mom dB
+  SYNC: dens tau DEnt mom dB
   SYNC: saved_prims
   SYNC: zvec
   SYNC: svec
@@ -86,15 +87,15 @@ SCHEDULE AsterX_Sync AT postregrid
 {
   LANG: C
   OPTIONS: global
-  SYNC: dens tau mom Avec_x Avec_y Avec_z Psi
+  SYNC: dens tau mom DEnt Avec_x Avec_y Avec_z Psi
   SYNC: saved_prims
 } "Synchronize"
 
 SCHEDULE AsterX_Sync IN ODESolvers_PostStep
 {
   LANG: C
   OPTIONS: global
-  SYNC: dens tau mom Avec_x Avec_y Avec_z Psi
+  SYNC: dens tau mom DEnt Avec_x Avec_y Avec_z Psi
 } "Synchronize"
 
 
@@ -124,21 +125,23 @@ SCHEDULE AsterX_Con2Prim IN AsterX_Con2PrimGroup AFTER AsterX_ComputedBFromdBsta
   LANG: C
   READS: ADMBaseX::metric(interior)
   READS: ADMBaseX::lapse(everywhere)
-  READS: dens(interior) tau(interior) mom(interior) dB(interior)
+  READS: dens(interior) tau(interior) DEnt(interior) mom(interior) dB(interior)
   READS: saved_prims(interior)
   READS: Avec_x(interior) Avec_y(interior) Avec_z(interior)
   WRITES: con2prim_flag(interior)
-  WRITES: HydroBaseX::rho(interior) HydroBaseX::vel(interior) HydroBaseX::eps(interior) HydroBaseX::press(interior) HydroBaseX::Bvec(interior)
+  WRITES: HydroBaseX::rho(interior) HydroBaseX::vel(interior) HydroBaseX::eps(interior) HydroBaseX::press(interior) HydroBaseX::Bvec(interior) 
+  WRITES: HydroBaseX::entropy(interior)
   WRITES: saved_prims(interior)
   WRITES: zvec(interior)
   WRITES: svec(interior)
-  WRITES: dens(interior) tau(interior) mom(interior) dB(interior)
+  WRITES: dens(interior) tau(interior) DEnt(interior) mom(interior) dB(interior)
   SYNC: con2prim_flag
-  SYNC: HydroBaseX::rho HydroBaseX::vel HydroBaseX::eps HydroBaseX::press HydroBaseX::Bvec
+  SYNC: HydroBaseX::rho HydroBaseX::vel HydroBaseX::eps HydroBaseX::press HydroBaseX::Bvec 
+  SYNC: HydroBaseX::entropy
   SYNC: saved_prims
   SYNC: zvec
   SYNC: svec
-  SYNC: dens tau mom dB
+  SYNC: dens tau DEnt mom dB
 } "Calculate primitive variables from conservative variables"
 
 
@@ -153,9 +156,10 @@ SCHEDULE AsterX_Fluxes IN AsterX_RHSGroup
   READS: ADMBaseX::metric(everywhere)
   READS: ADMBaseX::lapse(everywhere)
   READS: ADMBaseX::shift(everywhere)
-  READS: dens(everywhere) tau(everywhere) mom(everywhere)
+  READS: dens(everywhere) tau(everywhere) mom(everywhere) DEnt(everywhere)
   READS: HydroBaseX::rho(everywhere) HydroBaseX::vel(everywhere) HydroBaseX::press(everywhere) HydroBaseX::eps(everywhere)
   READS: HydroBaseX::Bvec(everywhere)
+  READS: HydroBaseX::entropy(everywhere)
   READS: dBx_stag(everywhere) dBy_stag(everywhere) dBz_stag(everywhere)
   READS: zvec_x(everywhere) zvec_y(everywhere) zvec_z(everywhere)
   READS: svec_x(everywhere) svec_y(everywhere) svec_z(everywhere)
@@ -181,8 +185,8 @@ SCHEDULE AsterX_SourceTerms IN AsterX_RHSGroup AFTER AsterX_Fluxes
   READS: zvec_x(everywhere), zvec_y(everywhere), zvec_z(everywhere)
   READS: svec_x(everywhere), svec_y(everywhere), svec_z(everywhere)
   READS: HydroBaseX::Bvec(everywhere)
-  WRITES: densrhs(interior) taurhs(interior) momrhs(interior)
-  SYNC: densrhs taurhs momrhs
+  WRITES: densrhs(interior) taurhs(interior) momrhs(interior) DEntrhs(interior)
+  SYNC: densrhs taurhs momrhs DEntrhs
 } "Calculate the source terms and compute the RHS of the hydro equations"
 
 SCHEDULE AsterX_RHS IN AsterX_RHSGroup AFTER AsterX_SourceTerms
@@ -191,15 +195,15 @@ SCHEDULE AsterX_RHS IN AsterX_RHSGroup AFTER AsterX_SourceTerms
   READS: ADMBaseX::metric(everywhere) ADMBaseX::lapse(everywhere) ADMBaseX::shift(everywhere)
   READS: HydroBaseX::vel(everywhere) HydroBaseX::press(everywhere)
   READS: flux_x(everywhere) flux_y(everywhere) flux_z(everywhere)
-  READS: densrhs(everywhere) taurhs(everywhere) momrhs(everywhere)
+  READS: densrhs(everywhere) taurhs(everywhere) momrhs(everywhere) DEntrhs(everywhere)
   READS: Psi(everywhere)
   READS: Aux_in_RHSof_A_Psi(everywhere)
   READS: dBx_stag(everywhere) dBy_stag(everywhere) dBz_stag(everywhere)
   READS: vtilde_xface(everywhere) vtilde_yface(everywhere) vtilde_zface(everywhere)
   READS: a_xface(everywhere) a_yface(everywhere) a_zface(everywhere)
-  WRITES: densrhs(interior) taurhs(interior) momrhs(interior)
+  WRITES: densrhs(interior) taurhs(interior) momrhs(interior) DEntrhs(interior)
   WRITES: Avec_x_rhs(interior) Avec_y_rhs(interior) Avec_z_rhs(interior) Psi_rhs(interior)
-  SYNC: densrhs taurhs momrhs
+  SYNC: densrhs taurhs momrhs DEntrhs
   SYNC: Avec_x_rhs Avec_y_rhs Avec_z_rhs Psi_rhs
 } "Update the RHS of the hydro equations with the flux contributions"