Add MOCMC

.github/workflows/tests.yml

@@ -30,22 +30,48 @@ jobs:
         - name: Run test scripts
           run: |
             cd tst/regression
-            mkdir build_minkowski
-            cd build_minkowski
+            mkdir build_minkowski_release
+            cd build_minkowski_release
             cmake -DCMAKE_BUILD_TYPE=Release -DPHOEBUS_ENABLE_UNIT_TESTS=ON -DMAX_NUMBER_CONSERVED_VARS=10 -DPHOEBUS_GEOMETRY=Minkowski -DPHOEBUS_CACHE_GEOMETRY=ON -DPARTHENON_DISABLE_HDF5_COMPRESSION=ON ../../../
             cmake --build . --parallel 4
             cd ..
-            mkdir build_fmks
-            cd build_fmks
+            mkdir build_fmks_release
+            cd build_fmks_release
             cmake -DCMAKE_BUILD_TYPE=Release -DPHOEBUS_ENABLE_UNIT_TESTS=ON -DMAX_NUMBER_CONSERVED_VARS=10 -DPHOEBUS_GEOMETRY=FMKS -DPHOEBUS_CACHE_GEOMETRY=ON -DPARTHENON_DISABLE_HDF5_COMPRESSION=ON ../../../
             cmake --build . --parallel 4
             cd ..
+            mkdir build_minkowski_debug
+            cd build_minkowski_debug
+            cmake -DCMAKE_BUILD_TYPE=Debug -DPHOEBUS_ENABLE_UNIT_TESTS=ON -DMAX_NUMBER_CONSERVED_VARS=10 -DPHOEBUS_GEOMETRY=Minkowski -DPHOEBUS_CACHE_GEOMETRY=ON -DPARTHENON_DISABLE_HDF5_COMPRESSION=ON ../../../
+            cmake --build . --parallel 4
+            cd ..
+            mkdir build_fmks_debug
+            cd build_fmks_debug
+            cmake -DCMAKE_BUILD_TYPE=Debug -DPHOEBUS_ENABLE_UNIT_TESTS=ON -DMAX_NUMBER_CONSERVED_VARS=10 -DPHOEBUS_GEOMETRY=FMKS -DPHOEBUS_CACHE_GEOMETRY=ON -DPARTHENON_DISABLE_HDF5_COMPRESSION=ON ../../../
+            cmake --build . --parallel 4
+            cd ..
             # TODO(BRR) add --use_gpu flag here if running on GPU
-            ./radiation_diffusion.py --executable build_minkowski/src/phoebus
-            ./bondi.py --executable build_fmks/src/phoebus
-            ./linear_modes.py --executable build_minkowski/src/phoebus
-            ./thincooling.py --executable build_minkowski/src/phoebus
-            ./thincooling_coolingfunction.py --executable build_minkowski/src/phoebus
+            # Release suite
+            ./radiation_diffusion.py --executable build_minkowski_release/src/phoebus
+            ./mocmc_diffusion.py --executable build_minkowski_release/src/phoebus
+            ./radiation_equilibration.py --executable build_minkowski_release/src/phoebus
+            ./mocmc_equilibration.py --executable build_minkowski_release/src/phoebus
+            ./bondi.py --executable build_fmks_release/src/phoebus
+            ./linear_modes.py --executable build_minkowski_release/src/phoebus
+            ./thincooling.py --executable build_minkowski_release/src/phoebus
+            ./thincooling_coolingfunction.py --executable build_minkowski_release/src/phoebus
+
+            # Debug suite
+            ./radiation_diffusion.py --executable build_minkowski_debug/src/phoebus
+            ./mocmc_diffusion.py --executable build_minkowski_debug/src/phoebus
+            ./radiation_equilibration.py --executable build_minkowski_debug/src/phoebus
+            ./mocmc_equilibration.py --executable build_minkowski_debug/src/phoebus
+            ./bondi.py --executable build_fmks_debug/src/phoebus
+            ./linear_modes.py --executable build_minkowski_debug/src/phoebus
+            ./thincooling.py --executable build_minkowski_debug/src/phoebus
+            ./thincooling_coolingfunction.py --executable build_minkowski_debug/src/phoebus
+
+            # Custom builds
             ./friedmann.py
         - name: Build code for unit testing
           run: |

inputs/radiation_advection.pin

@@ -12,8 +12,8 @@
 
 <phoebus>
 problem = radiation_advection
-ix1_bc = outflow
-ox1_bc = outflow
+#ix1_bc = outflow
+#ox1_bc = outflow
 
 <parthenon/job>
 problem_id  = rad_adv   # problem ID: basename of output filenames
@@ -22,7 +22,7 @@ problem_id  = rad_adv   # problem ID: basename of output filenames
 variables = r.c.E, &
             r.c.F, &
             r.p.J, &
-            r.p.H 
+            r.p.H
 
 file_type   = hdf5      # Tabular data dump
 dt          = 1.e-1     # time increment between outputs
@@ -42,8 +42,13 @@ nghost = 4
 nx1         = 512      # Number of zones in X1-direction
 x1min       = 0        # minimum value of X1
 x1max       = 1        # maximum value of X1
-ix1_bc      = user     # Inner-X1 boundary condition flag
-ox1_bc      = user     # Outer-X1 boundary condition flag
+#ix1_bc      = user     # Inner-X1 boundary condition flag
+#ox1_bc      = user     # Outer-X1 boundary condition flag
+# TODO(BRR) Currently an issue with asymptotic preserving transport in samples
+# in pure scattering problems that makes this problem misbehave with outflow
+# boundaries
+ix1_bc = periodic
+ox1_bc = periodic
 
 nx2         = 1         # Number of zones in X2-direction
 x2min       = -1        # minimum value of X2
@@ -59,7 +64,7 @@ ox3_bc      = periodic  # Outer-X3 boundary condition flfgag
 
 num_threads = 1         # maximum number of OMP threads
 
-bc_vars = primitive 
+bc_vars = primitive
 
 <parthenon/meshblock>
 nx1 = 512
@@ -74,7 +79,7 @@ max_level = 3
 <eos>
 type = IdealGas
 Gamma = 1.6666666666666666667
-Cv = 1.0 
+Cv = 1.0
 
 <physics>
 hydro = false
@@ -91,23 +96,24 @@ c2p_max_iter = 100
 Ye = true
 
 <radiation>
-cfl = 0.5 
+cfl = 0.5
 method = moment_m1
-nu_min = 1.e15
-nu_max = 1.e22
-nu_bins = 200
+nu_min = 1.e-2
+nu_max = 1.e2
+nu_bins = 10
 absorption = true
 do_nu_electron = true
 do_nu_electron_anti = false
 do_nu_heavy = false
 
-scattering_fraction = 1.0 
-B_fake = 0.5
-use_B_fake = true
+scattering_fraction = 1.0
 
-<opacity> 
-type = gray 
-gray_kappa = 1.e3 
+<radiation/mocmc>
+nsamp_per_zone = 32
+
+<opacity>
+type = gray
+gray_kappa = 1.e3
 
 <radiation_advection>
 J = 1.0
@@ -116,5 +122,5 @@ Hy = 0.0
 Hz = 0.0
 vx = 0.3
 width = 0.0333
-kappas_init = 1.e3
+tau = 1.e3
 boost_profile = true

inputs/radiation_equilibration.pin

@@ -19,9 +19,13 @@ bc_ox1 = outflow
 problem_id  = rad_eql   # problem ID: basename of output filenames
 
 <parthenon/output1>
-variables = r.c.E, &
+variables = p.density, &
+            r.p.J, &
+            r.p.H, &
+            r.c.E, &
             r.c.F, &
-            c.energy 
+            p.energy, &
+            c.energy
 
 file_type   = hdf5      # Tabular data dump
 dt          = 1.e-1     # time increment between outputs
@@ -38,11 +42,11 @@ nghost = 4
 #refinement = adaptive
 #numlevel = 3
 
-nx1         = 32      # Number of zones in X1-direction
+nx1         = 4        # Number of zones in X1-direction
 x1min       = 0        # minimum value of X1
 x1max       = 1        # maximum value of X1
-ix1_bc      = periodic     # Inner-X1 boundary condition flag
-ox1_bc      = periodic     # Outer-X1 boundary condition flag
+ix1_bc      = periodic # Inner-X1 boundary condition flag
+ox1_bc      = periodic # Outer-X1 boundary condition flag
 
 nx2         = 1         # Number of zones in X2-direction
 x2min       = -1        # minimum value of X2
@@ -58,10 +62,10 @@ ox3_bc      = periodic  # Outer-X3 boundary condition flfgag
 
 num_threads = 1         # maximum number of OMP threads
 
-bc_vars = primitive 
+bc_vars = primitive
 
 <parthenon/meshblock>
-nx1 = 32
+nx1 = 4
 nx2 = 1
 nx3 = 1
 
@@ -73,7 +77,7 @@ max_level = 3
 <eos>
 type = IdealGas
 Gamma = 1.6666666666666666667
-Cv = 1.0 
+Cv = 1.0
 
 <physics>
 hydro = true
@@ -90,23 +94,19 @@ c2p_max_iter = 100
 Ye = true
 
 <radiation>
-cfl = 0.8 
+cfl = 0.8
 method = moment_eddington
-nu_min = 1.e15
-nu_max = 1.e22
-nu_bins = 200
+nu_min = 1.e-2
+nu_max = 1.e2
+nu_bins = 10
 absorption = true
 do_nu_electron = true
 do_nu_electron_anti = false
 do_nu_heavy = false
 
-scattering_fraction = 0.0 
-B_fake = 0.5
-use_B_fake = false
-
-<opacity> 
-type = gray 
-gray_kappa = 1.0
+<opacity>
+type = gray
+gray_kappa = 0.2
 
 <radiation_equilibration>
 J = 0.2

inputs/thincooling.pin

@@ -13,7 +13,7 @@ variables = p.density,  &
             temperature
 
 file_type   = hdf5      # Tabular data dump
-dt          = 1.e-2     # time increment between outputs
+dt          = 1.e-3     # time increment between outputs
 
 <parthenon/time>
 nlim        = -1       # cycle limit
@@ -27,8 +27,7 @@ nghost = 4
 
 nx1         = 4         # Number of zones in X1-direction
 x1min       = 0         # minimum value of X1
-#x1max       = 1.e-7     # maximum value of X1 (for cooling_function)
-x1max       = 1.e-4     # maximum value of X1 (for monte_carlo)
+x1max       = 5.e-4     # maximum value of X1 (for monte_carlo)
 ix1_bc      = periodic  # Inner-X1 boundary condition flag
 ox1_bc      = periodic  # Outer-X1 boundary condition flag
 
@@ -81,8 +80,10 @@ Ye = true
 method = monte_carlo
 nu_min = 1.e15
 nu_max = 0.9999e17
-nu_bins = 200
-tune_emission = 1.e-39
+nu_bins = 100
+tune_emission = 1.e-43
+#tune_emission = 1.e-41
+#tune_emission = 1.e-34
 num_particles = 16
 opacity_model = tophat
 absorption = false

scripts/python/neutrino_cooling.py

@@ -20,25 +20,23 @@
 import os
 from subprocess import call, DEVNULL
 import glob
-from parthenon_tools import phdf
 import time
 from enum import Enum
+from phoedf import phoedf
 
 NeutrinoSpecies = Enum('NeutrinoSpecies', 'electron electronanti')
 
 dfnams = np.sort(glob.glob(DUMP_NAMES))
-dfile0 = phdf.phdf(dfnams[0])
+dfile0 = phoedf(dfnams[0])
 
 cl = 2.99792458e10
 mp = 1.672621777e-24
 h = 6.62606957e-27
 
-T_unit = dfile0.Params['eos/time_unit']
-L_unit = dfile0.Params['eos/length_unit']
-M_unit = dfile0.Params['eos/mass_unit']
-U_unit = M_unit / L_unit**3 * cl**2
-numin = dfile0.Params['radiation/nu_min']
-numax = dfile0.Params['radiation/nu_max']
+T_unit = dfile0.TimeCodeToCGS
+U_unit = dfile0.EnergyDensityCodeToCGS
+numin = dfile0.Params['opacity/numin']
+numax = dfile0.Params['opacity/numax']
 do_nu_electron = dfile0.Params['radiation/do_nu_electron']
 do_nu_electron_anti = dfile0.Params['radiation/do_nu_electron_anti']
 do_nu_heavy = dfile0.Params['radiation/do_nu_heavy']
@@ -84,7 +82,7 @@ def get_u(t):
 Ye_code = np.zeros(dfnams.size)
 u_code = np.zeros(dfnams.size)
 for n, dfnam in enumerate(dfnams):
-  dfile = phdf.phdf(dfnam)
+  dfile = phoedf(dfnam)
   t_code[n] = dfile.Time*T_unit
   Ye_code[n] = dfile.Get("p.ye").mean()
   u_code[n] = dfile.Get("p.energy").mean()*U_unit

scripts/python/phoebus_constants.py

@@ -0,0 +1,79 @@
+# © 2022. Triad National Security, LLC. All rights reserved.  This
+# program was produced under U.S. Government contract
+# 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
+# is operated by Triad National Security, LLC for the U.S.  Department
+# of Energy/National Nuclear Security Administration. All rights in
+# the program are reserved by Triad National Security, LLC, and the
+# U.S. Department of Energy/National Nuclear Security
+# Administration. The Government is granted for itself and others
+# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+# license in this material to reproduce, prepare derivative works,
+# distribute copies to the public, perform publicly and display
+# publicly, and to permit others to do so.
+
+from math import pi
+cgs = {}
+cgs['KM']      = 1e5
+cgs['CL']      = 2.99792458e10
+cgs['QE']      = 4.80320680e-10
+cgs['ME']      = 9.1093826e-28
+cgs['MP']      = 1.67262171e-24
+cgs['MN']      = 1.67492728e-24
+cgs['HPL']     = 6.6260693e-27
+cgs['HBAR']    = 1.0545717e-27
+cgs['KBOL']    = 1.3806505e-16
+cgs['GNEWT']   = 6.6742e-8
+cgs['SIG']     = 5.670400e-5
+cgs['AR']      = 7.5657e-15
+cgs['THOMSON'] = 0.665245873e-24
+cgs['JY']      = 1.e-23
+cgs['PC']      = 3.085678e18
+cgs['AU']      = 1.49597870691e13
+cgs['MSOLAR']  = 1.989e33
+cgs['RSOLAR']  = 6.96e10
+cgs['LSOLAR']  = 3.827e33
+cgs['EV']      = 1.60217653e-12
+cgs['MEV']     = 1.0e6*cgs['EV']
+cgs['GEV']     = 1.0e9*cgs['EV']
+cgs['K']       = 1.380648780669e-16
+cgs['GK']      = 1.0e9*cgs['K']
+cgs['GFERM']   = 1.435850814907447e-49
+cgs['GA']      = -1.272323
+cgs['S2THW']   = 0.222321
+cgs['alphafs'] = 1./137.
+cgs['NUSIGMA0'] = (4.0*(cgs['GFERM']**2)
+                   *((cgs['ME']*cgs['CL']**2)**2)
+                   /(pi*(cgs['HBAR']*cgs['CL'])**4))
+
+scalefree = {}
+scalefree['KM']      = 1.
+scalefree['CL']      = 1.
+scalefree['QE']      = 1.
+scalefree['ME']      = 1.
+scalefree['MP']      = 1.
+scalefree['MN']      = 1.
+scalefree['HPL']     = 1.
+scalefree['HBAR']    = scalefree['HPL']/(2.*pi)
+scalefree['KBOL']    = 1.
+scalefree['GNEWT']   = 1.
+scalefree['SIG']     = 1.
+scalefree['AR']      = 1.
+scalefree['THOMSON'] = 1.
+scalefree['JY']      = 1.
+scalefree['PC']      = 1.
+scalefree['AU']      = 1.
+scalefree['MSOLAR']  = 1.
+scalefree['RSOLAR']  = 1.
+scalefree['LSOLAR']  = 1.
+scalefree['EV']      = 1.
+scalefree['MEV']     = 1.0e6*scalefree['EV']
+scalefree['GEV']     = 1.0e9*scalefree['EV']
+scalefree['K']       = 1.
+scalefree['GK']      = 1.0e9*scalefree['K']
+scalefree['GFERM']   = 1.
+scalefree['GA']      = -1.
+scalefree['S2THW']   = 1.
+scalefree['alphafs'] = 1./137.
+scalefree['NUSIGMA0'] = (4.0*(scalefree['GFERM']**2)
+                   *((scalefree['ME']*scalefree['CL']**2)**2)
+                   /(pi*(scalefree['HBAR']*scalefree['CL'])**4))