Skip to content

Commit

Permalink
Merge pull request #80 from lanl/brryan/add_mocmc
Browse files Browse the repository at this point in the history
Add MOCMC
  • Loading branch information
brryan authored Jun 8, 2022
2 parents b8a9a4b + 6b01f91 commit c570f42
Show file tree
Hide file tree
Showing 62 changed files with 4,958 additions and 1,888 deletions.
44 changes: 35 additions & 9 deletions .github/workflows/tests.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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: |
Expand Down
2 changes: 1 addition & 1 deletion external/parthenon
42 changes: 24 additions & 18 deletions inputs/radiation_advection.pin
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -74,7 +79,7 @@ max_level = 3
<eos>
type = IdealGas
Gamma = 1.6666666666666666667
Cv = 1.0
Cv = 1.0

<physics>
hydro = false
Expand All @@ -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
Expand All @@ -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
38 changes: 19 additions & 19 deletions inputs/radiation_equilibration.pin
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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

Expand All @@ -73,7 +77,7 @@ max_level = 3
<eos>
type = IdealGas
Gamma = 1.6666666666666666667
Cv = 1.0
Cv = 1.0

<physics>
hydro = true
Expand All @@ -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
11 changes: 6 additions & 5 deletions inputs/thincooling.pin
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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
Expand Down
16 changes: 7 additions & 9 deletions scripts/python/neutrino_cooling.py
Original file line number Diff line number Diff line change
Expand Up @@ -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']
Expand Down Expand Up @@ -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
Expand Down
79 changes: 79 additions & 0 deletions scripts/python/phoebus_constants.py
Original file line number Diff line number Diff line change
@@ -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))
Loading

0 comments on commit c570f42

Please sign in to comment.