Merge pull request #104 from lanl/brryan/rad_torus

Radiative torus

cmake/Fluid.cmake

@@ -18,14 +18,6 @@ else()
   set(USE_VALENCIA 0 CACHE BOOL "Covariant formulation of GRMHD")
 endif()
 
-option(PHOEBUS_C2P_ROBUST_FLOORS "Use full floors for c2p robust" ON)
-
-if(PHOEBUS_C2P_ROBUST_FLOORS)
-  set(USE_C2P_ROBUST_FLOORS 1 CACHE BOOL "Use regular floors for c2p robust")
-else()
-  set(USE_C2P_ROBUST_FLOORS 0 CACHE BOOL "Use only sanity floors for c2p robust")
-endif()
-
 option(PHOEBUS_FLUX_AND_SRC_DIAGS "Stash flux divergence and sources for possible output" OFF)
 if(PHOEBUS_FLUX_AND_SRC_DIAGS)
   set(SET_FLUX_SRC_DIAGS 1 CACHE BOOL "Allocate and set these diagnostics")

inputs/linear_modes.pin

@@ -112,7 +112,6 @@ riemann = hll
 c2p_max_iter = 100
 recon = linear
 c2p_tol = 1.e-12
-#c2p_method = classic
 c2p_method = robust
 
 <hydro_modes>

inputs/radiation_advection.pin

@@ -105,6 +105,7 @@ absorption = true
 do_nu_electron = true
 do_nu_electron_anti = false
 do_nu_heavy = false
+src_solver = oned
 
 scattering_fraction = 1.0
 

inputs/torus.pin

@@ -125,6 +125,7 @@ sie_exp_floor = -1.0
 ceiling_type = ConstantGamSie
 sie0_ceiling = 1.0
 gam0_ceiling = 10.0
+enable_ox1_fmks_inflow_check = true
 
 <torus>
 target_beta = 100.

inputs/torus_rad.pin

@@ -0,0 +1,191 @@
+
+<phoebus>
+problem = torus
+ix1_bc = outflow
+ox1_bc = outflow
+ix2_bc = reflect
+ox2_bc = reflect
+
+<phoebus/mesh>
+bc_vars = primitive
+
+<parthenon/job>
+problem_id  = torus       # problem ID: basename of output filenames
+
+<parthenon/output0>
+file_type = rst
+dt = 200.0
+
+<parthenon/output1>
+variables = p.density,       &
+#            c.density,       &
+            p.velocity,      &
+#            c.momentum,      &
+            p.energy,        &
+#            c.energy,        &
+#            pressure,      	 &
+            r.p.J,           &
+            r.p.H,           &
+#            r.c.E,           &
+#            r.c.F,           &
+            g.c.coord,       &
+            g.n.coord,       &
+            g.c.detgam,      &
+#            g.c.bcon,        &
+            g.c.alpha,       &
+            g.c.gcov,        &
+#            src_terms,       &
+#            flux_divergence, &
+            r.i.kappaH,      &
+            p.bfield
+
+file_type   = hdf5      # Tabular data dump
+dt          = 5.0       # time increment between outputs
+
+<parthenon/time>
+nlim        = -1 # cycle limit
+tlim        = 2000.       # time limit
+integrator  = rk2       # time integration algorithm
+ncycle_out  = 1         # interval for stdout summary info
+dt_init_fact = 1.e-6
+
+<parthenon/mesh>
+nghost = 4
+#refinement = adaptive
+#numlevel = 3
+
+nx1         = 256           # Number of zones in X1-direction
+#nx1 = 128
+x1min       = 0.59          # minimum value of X1; overwritten for the torus problem
+x1max       = 3.69          # maximum value of X1; set by <coordinates>/r_outer for the torus problem
+ix1_bc      = user          # Inner-X1 boundary condition flag
+ox1_bc      = user          # Outer-X1 boundary condition flag
+
+nx2         = 256           # Number of zones in X2-direction
+#nx2 = 128
+x2min       = 0             # minimum value of X2
+x2max       = 1             # maximum value of X2. X2 = 1 -> theta = pi
+ix2_bc      = user          # Inner-X2 boundary condition flag
+ox2_bc      = user          # Outer-X2 boundary condition flag
+
+nx3         = 1             # Number of zones in X3-direction
+#nx3 = 64
+x3min       = 0             # minimum value of X3
+x3max       = 6.28318530718 # maximum value of X3. 2*pi
+ix3_bc      = periodic      # Inner-X3 boundary condition flag
+ox3_bc      = periodic      # Outer-X3 boundary condition flfgag
+
+num_threads = 1         # maximum number of OMP threads
+
+# Commenting out meshblock sizes makes the simulation 1 meshblock
+# <parthenon/meshblock>
+# nx1 = 128
+# nx2 = 128
+# nx3 = 1
+
+<parthenon/refinement0>
+field = c.c.bulk.rho
+method = derivative_order_1
+max_level = 3
+
+<eos>
+type = IdealGas
+Gamma = 1.66666666666666666666666667
+Cv = 1.0
+
+<physics>
+hydro = true
+rad = true
+
+<fluid>
+xorder = 2
+cfl = 0.8
+riemann = llf
+recon = weno5
+c2p_max_iter = 100
+c2p_tol = 1.e-8
+c2p_floor_scale_fac = 1.e-10
+c2p_fail_on_floors = false
+c2p_fail_on_ceilings = false
+c2p_method = robust
+Ye = false
+mhd = true
+
+<geometry>
+a = 0.9375
+
+<coordinates>
+derefine_poles = false
+r_outer = 250
+
+<radiation>
+cfl = 0.5
+recon = weno5
+method = moment_eddington
+absorption = true
+do_nu_electron = true
+do_nu_electron_anti = false
+do_nu_heavy = false
+src_solver = oned # fourd
+src_use_oned_backup = false # true
+src_rootfind_eps = 1.e-6
+src_rootfind_tol = 1.e-8
+src_rootfind_maxiter = 100
+oned_fixup_strategy = ignore_dJ
+recon_fixup_strategy = bounds
+closure_c2p_strategy = robust
+
+#scattering_fraction = 1.0
+scattering_fraction = 0.5
+
+<opacity>
+type = gray
+gray_kappa = 1.e-10 # midplane tau/zone >1e3
+
+<units>
+scale_free = false
+geom_mass_msun = 1.
+fluid_mass_g = 1.e26
+
+<fixup>
+enable_floors = true
+enable_mhd_floors = true
+enable_rad_floors = true
+enable_ceilings = true
+enable_rad_ceilings = true
+enable_flux_fixup = false # true
+enable_c2p_fixup = true
+enable_source_fixup = true
+report_c2p_fails = false
+report_source_fails = false
+c2p_failure_force_fixup_both = false
+enable_ox1_fmks_inflow_check = true
+
+fluid_c2p_failure_strategy = interpolate
+rad_c2p_failure_strategy = interpolate
+src_failure_strategy = floors
+
+floor_type = ExpX1RhoSie
+rho0_floor = 1.0e-5
+sie0_floor = 1.0e-2
+rho_exp_floor = -2.0
+sie_exp_floor = -1.0
+ceiling_type = ConstantGamSie
+sie0_ceiling = 1.0
+gam0_ceiling = 10.0
+rad_floor_type = ExpX1J
+J0_floor = 1.e-6
+J_exp_floor = -3.0
+rad_ceiling_type = ConstantXi0
+xi0_ceiling = 0.95
+
+<torus>
+target_beta = 100.
+harm_beta_normalization = true
+u_jitter = 0.05
+kappa = 0.01
+rin = 6.0
+rmax = 12.0
+n_inside_horizon = 5
+magnetized = true
+initial_radiation = thermal

scripts/python/phoebus_eos.py

@@ -26,6 +26,9 @@ def T_from_rho_u_Ye(self, rho, u, Ye):
   def u_from_rho_T_Ye(self, rho, T, Ye):
     raise NotImplementedError()
 
+  def P_from_rho_u_Ye(self, rho, u, Ye):
+    raise NotImplementedError()
+
 # -- Ideal gas EOS
 class IdealEOS(BaseEOS):
   def __init__(self, params):
@@ -40,6 +43,9 @@ def T_from_rho_u_Ye(self, rho, u, Ye):
   def u_from_rho_T_Ye(self, rho, T, Ye):
     return rho * self.Cv * T
 
+  def P_from_rho_u_Ye(self, rho, u, Ye):
+    return self.gm1 * u
+
 # ---------------------------------------------------------------------------- #
 # -- Dictionary of EOS types
 eos_type_dict = {'IdealGas' : IdealEOS}