Add non-relativistic radiation transport solvers (explicit and implicit) and cosmic ray transport

.gitignore

@@ -1,6 +1,7 @@
 # Configured outputs
 Makefile
 src/defs.hpp
+configure.log
 
 # Object and binary directories
 obj/

Makefile.in

@@ -30,6 +30,7 @@ SRC_FILES := $(wildcard src/*.cpp) \
 	     $(wildcard src/bvals/cc/fft_grav/*.cpp) \
 	     $(wildcard src/bvals/cc/hydro/*.cpp) \
 	     $(wildcard src/bvals/cc/mg/*.cpp) \
+             $(wildcard src/bvals/cc/nr_radiation/*.cpp) \
 	     $(wildcard src/bvals/fc/*.cpp) \
 	     $(wildcard src/bvals/orbital/*.cpp) \
 	     $(wildcard src/bvals/utils/*.cpp) \
@@ -45,6 +46,12 @@ SRC_FILES := $(wildcard src/*.cpp) \
 	     $(wildcard src/hydro/*.cpp) \
 	     $(wildcard src/hydro/srcterms/*.cpp) \
 	     $(wildcard src/hydro/hydro_diffusion/*.cpp) \
+             $(wildcard src/nr_radiation/*.cpp) \
+             $(wildcard src/nr_radiation/integrators/*.cpp) \
+             $(wildcard src/nr_radiation/integrators/srcterms/*.cpp) \
+             $(wildcard src/nr_radiation/implicit/*.cpp) \
+             $(wildcard src/cr/*.cpp) \
+             $(wildcard src/cr/integrators/*.cpp) \
 	     src/hydro/rsolvers/$(RSOLVER_DIR)$(RSOLVER_FILE) \
 	     $(wildcard src/inputs/*.cpp) \
 	     $(wildcard src/mesh/*.cpp) \

README.md

@@ -17,7 +17,7 @@ athena
 	  <img width="345" height="345" src="https://user-images.githubusercontent.com/1410981/115276281-759d8580-a108-11eb-9fc9-833480b97f95.png">
 </p>
 
-Athena++ GRMHD code and adaptive mesh refinement (AMR) framework
+Athena++ radiation GRMHD code and adaptive mesh refinement (AMR) framework
 
 Please read [our contributing guidelines](./CONTRIBUTING.md) for details on how to participate.
 

configure.py

@@ -7,36 +7,39 @@
 # Makefile.in and src/defs.hpp.in respectively.
 #
 # The following options are implememted:
-#   -h  --help        help message
-#   --prob=name       use src/pgen/name.cpp as the problem generator
-#   --coord=xxx       use xxx as the coordinate system
-#   --eos=xxx         use xxx as the equation of state
-#   --flux=xxx        use xxx as the Riemann solver
-#   --nghost=xxx      set NGHOST=xxx
-#   --nscalars=xxx    set NSCALARS=xxx
-#   -eos_table        enable EOS table
-#   -b                enable magnetic fields
-#   -s                enable special relativity
-#   -g                enable general relativity
-#   -t                enable interface frame transformations for GR
-#   -debug            enable debug flags (-g -O0); override other compiler options
-#   -coverage         enable compiler-dependent code coverage flags
-#   -float            enable single precision (default is double)
-#   -mpi              enable parallelization with MPI
-#   -omp              enable parallelization with OpenMP
-#   -hdf5             enable HDF5 output (requires the HDF5 library)
-#   --hdf5_path=path  path to HDF5 libraries (requires the HDF5 library)
-#   -fft              enable FFT (requires the FFTW library)
-#   --fftw_path=path  path to FFTW libraries (requires the FFTW library)
-#   --grav=xxx        use xxx as the self-gravity solver
-#   --cxx=xxx         use xxx as the C++ compiler (works w/ or w/o -mpi)
-#   --ccmd=name       use name as the command to call the (non-MPI) C++ compiler
-#   --mpiccmd=name    use name as the command to call the MPI C++ compiler
-#   --gcovcmd=name    use name as the command to call the gcov utility
-#   --cflag=string    append string whenever invoking compiler/linker
-#   --include=path    use -Ipath when compiling
-#   --lib_path=path   use -Lpath when linking
-#   --lib=xxx         use -lxxx when linking
+#   -h  --help           help message
+#   --prob=name          use src/pgen/name.cpp as the problem generator
+#   --coord=xxx          use xxx as the coordinate system
+#   --eos=xxx            use xxx as the equation of state
+#   --flux=xxx           use xxx as the Riemann solver
+#   --nghost=xxx         set NGHOST=xxx
+#   --nscalars=xxx       set NSCALARS=xxx
+#   -eos_table           enable EOS table
+#   -b                   enable magnetic fields
+#   -s                   enable special relativity
+#   -g                   enable general relativity
+#   -t                   enable interface frame transformations for GR
+#   -debug               enable debug flags (-g -O0); override other compiler options
+#   -coverage            enable compiler-dependent code coverage flags
+#   -float               enable single precision (default is double)
+#   -mpi                 enable parallelization with MPI
+#   -omp                 enable parallelization with OpenMP
+#   -hdf5                enable HDF5 output (requires the HDF5 library)
+#   --hdf5_path=path     path to HDF5 libraries (requires the HDF5 library)
+#   -fft                 enable FFT (requires the FFTW library)
+#   --fftw_path=path     path to FFTW libraries (requires the FFTW library)
+#   --grav=xxx           use xxx as the self-gravity solver
+#   --cxx=xxx            use xxx as the C++ compiler (works w/ or w/o -mpi)
+#   --ccmd=name          use name as the command to call the (non-MPI) C++ compiler
+#   --mpiccmd=name       use name as the command to call the MPI C++ compiler
+#   --gcovcmd=name       use name as the command to call the gcov utility
+#   --cflag=string       append string whenever invoking compiler/linker
+#   --include=path       use -Ipath when compiling
+#   --lib_path=path      use -Lpath when linking
+#   --lib=xxx            use -lxxx when linking
+#   -nr_radiation        turn on radiation transport
+#   -implicit_radiation  implicit radiation transport module
+#   -cr                  enable cosmic ray transport
 # ----------------------------------------------------------------------------------------
 
 # Modules
@@ -205,6 +208,24 @@
                     default='',
                     help='path to HDF5 libraries')
 
+# -radiation argument
+parser.add_argument('-nr_radiation',
+                    action='store_true',
+                    default=False,
+                    help='enable non-relativistic radiative transfer')
+
+# -radiation argument
+parser.add_argument('-implicit_radiation',
+                    action='store_true',
+                    default=False,
+                    help='enable radiative transfer')
+
+# -cosmic ray argument
+parser.add_argument('-cr',
+                    action='store_true',
+                    default=False,
+                    help='enable cosmic ray transport')
+
 # The main choices for --cxx flag, using "ctype[-suffix]" formatting, where "ctype" is the
 # major family/suite/group of compilers and "suffix" may represent variants of the
 # compiler version and/or predefined sets of compiler options. The C++ compiler front ends
@@ -444,6 +465,31 @@ def c_to_cpp(arg):
     if not args['t']:
         makefile_options['RSOLVER_FILE'] += '_no_transform'
 
+
+# -radiation argument
+definitions['NRAD_VARIABLES'] = '0'
+
+if args['nr_radiation']:
+    definitions['NR_RADIATION_ENABLED'] = '1'
+    definitions['NRAD_VARIABLES'] = '14'
+else:
+    definitions['NR_RADIATION_ENABLED'] = '0'
+
+if args['implicit_radiation']:
+    definitions['IM_RADIATION_ENABLED'] = '1'
+    definitions['NRAD_VARIABLES'] = '14'
+else:
+    definitions['IM_RADIATION_ENABLED'] = '0'
+
+# -cr argument
+definitions['NCR_VARIABLES'] = '0'
+if args['cr']:
+    definitions['CR_ENABLED'] = '1'
+    definitions['NCR_VARIABLES'] = '4'
+else:
+    definitions['CR_ENABLED'] = '0'
+
+
 # --cxx=[name] argument
 if args['cxx'] == 'g++':
     # GCC is C++11 feature-complete since v4.8.1 (2013-05-31)
@@ -817,6 +863,9 @@ def c_to_cpp(arg):
 print('  Number of scalars:          ' + args['nscalars'])
 print('  Special relativity:         ' + ('ON' if args['s'] else 'OFF'))
 print('  General relativity:         ' + ('ON' if args['g'] else 'OFF'))
+print('  Radiative Transfer:         ' + ('ON' if args['nr_radiation'] else 'OFF'))
+print('  Implicit Radiation:         ' + ('ON' if args['implicit_radiation'] else 'OFF'))
+print('  Cosmic Ray Transport:       ' + ('ON' if args['cr'] else 'OFF'))
 print('  Frame transformations:      ' + ('ON' if args['t'] else 'OFF'))
 print('  Self-Gravity:               ' + self_grav_string)
 print('  Super-Time-Stepping:        ' + ('ON' if args['sts'] else 'OFF'))
@@ -835,3 +884,44 @@ def c_to_cpp(arg):
 print('  Compiler:                   ' + args['cxx'])
 print('  Compilation command:        ' + makefile_options['COMPILER_COMMAND'] + ' '
       + makefile_options['PREPROCESSOR_FLAGS'] + ' ' + makefile_options['COMPILER_FLAGS'])
+
+# write the configuration optitions into a log file
+flog = open('./configure.log', 'w')
+flog.write('Your Athena++ distribution has now been configured with the following options:' + '\n')  # noqa
+flog.write('  Problem generator:          ' + args['prob'] + '\n')
+flog.write('  Coordinate system:          ' + args['coord'] + '\n')
+flog.write('  Equation of state:          ' + args['eos'] + '\n')
+flog.write('  Riemann solver:             ' + args['flux'] + '\n')
+flog.write('  Magnetic fields:            ' + ('ON' if args['b'] else 'OFF') + '\n')
+flog.write('  Number of scalars:          ' + args['nscalars'] + '\n')
+flog.write('  Special relativity:         ' + ('ON' if args['s'] else 'OFF') + '\n')
+flog.write('  General relativity:         ' + ('ON' if args['g'] else 'OFF') + '\n')
+flog.write('  Radiative Transfer:         ' + ('ON' if args['nr_radiation'] else 'OFF')
+           + '\n')
+flog.write('  Implicit Radiation:         ' + ('ON' if args['implicit_radiation']
+                                               else 'OFF') + '\n')
+flog.write('  Cosmic Ray Transport:       ' + ('ON' if args['cr'] else 'OFF') + '\n')
+flog.write('  Frame transformations:      ' + ('ON' if args['t'] else 'OFF') + '\n')
+flog.write('  Self-Gravity:               ' + self_grav_string + '\n')
+flog.write('  Super-Time-Stepping:        ' + ('ON' if args['sts'] else 'OFF') + '\n')
+flog.write('  Debug flags:                ' + ('ON' if args['debug'] else 'OFF') + '\n')
+flog.write('  Code coverage flags:        ' + ('ON' if args['coverage'] else 'OFF')
+           + '\n')
+flog.write('  Linker flags:               ' + makefile_options['LINKER_FLAGS'] + ' '
+           + makefile_options['LIBRARY_FLAGS'] + '\n')
+flog.write('  Floating-point precision:   ' + ('single' if args['float'] else 'double')
+           + '\n')
+flog.write('  Number of ghost cells:      ' + args['nghost'] + '\n')
+flog.write('  MPI parallelism:            ' + ('ON' if args['mpi'] else 'OFF') + '\n')
+flog.write('  OpenMP parallelism:         ' + ('ON' if args['omp'] else 'OFF') + '\n')
+flog.write('  FFT:                        ' + ('ON' if args['fft'] else 'OFF') + '\n')
+flog.write('  HDF5 output:                ' + ('ON' if args['hdf5'] else 'OFF') + '\n')
+if args['hdf5']:
+    flog.write('  HDF5 precision:             ' + ('double' if args['h5double']
+                                                   else 'single') + '\n')
+flog.write('  Compiler:                   ' + args['cxx'] + '\n')
+flog.write('  Compilation command:        ' + makefile_options['COMPILER_COMMAND'] + ' '
+           + makefile_options['PREPROCESSOR_FLAGS'] + ' '
+           + makefile_options['COMPILER_FLAGS'] + '\n')
+
+flog.close()

inputs/cosmic_ray/athinput.cr_diffusion

@@ -0,0 +1,74 @@
+<comment>
+problem   = cosmic ray diffusion
+reference = 
+configure = --prob=cr_diffusion
+
+<job>
+problem_id = cr  # problem ID: basename of output filenames
+
+#<output1>
+#file_type = hdf5  # Binary data dump
+#dt        = 0.2   # time increment between outputs
+#variable  = prim
+
+#<output2>
+#file_type = hst    # Binary data dump
+#dt        = 0.001  # time increment between outputs
+
+#<output3>
+#file_type = rst   # Binary data dump
+#dt        = 0.02  # time increment between outputs
+
+<time>
+cfl_number = 0.3  # The Courant, Friedrichs, & Lewy (CFL) Number
+nlim       = -1   # cycle limit
+tlim       = 0.2  # time limit
+
+<mesh>
+nx1    = 16        # Number of zones in X1-direction
+x1min  = -1.0      # minimum value of X1
+x1max  = 1.0       # maximum value of X1
+ix1_bc = periodic  # inner-X1 boundary flag
+ox1_bc = periodic  # inner-X1 boundary flag
+
+nx2    = 256      # Number of zones in X2-direction
+x2min  = -1.0     # minimum value of X2
+x2max  = 1.0      # maximum value of X2
+ix2_bc = outflow  # inner-X2 boundary flag
+ox2_bc = outflow  # inner-X2 boundary flag
+
+nx3    = 1         # Number of zones in X3-direction
+x3min  = -1.0      # minimum value of X3
+x3max  = 1.0       # maximum value of X3
+ix3_bc = periodic  # inner-X3 boundary flag
+ox3_bc = periodic  # inner-X3 boundary flag
+
+#refinement = static
+
+<meshblock>
+nx1 = 16 
+nx2 = 256
+nx3 = 1
+
+<refinement1>
+x1min = 1
+x1max = 23.4414
+x2min = 0.698132 
+x2max = 2.44346
+x3min = 0.0
+x3max = 6.283185307179586 
+level = 2
+
+<hydro>
+gamma  = 1.6666666666667  # gamma = C_p/C_v
+dfloor = 1.e-8
+pfloor = 1.e-7
+
+<cr>
+vmax     = 100
+src_flag = 0
+
+<problem>
+v0        = 0
+sigma     = 1.e3
+direction = 1

inputs/radiation/athinput.atmosphere

@@ -0,0 +1,74 @@
+<comment>
+problem   = Radiation beam test
+reference = 
+configure = --prob=beam
+
+<job>
+problem_id = atmosphere  # problem ID: basename of output filenames
+
+#<output1>
+#file_type = hst   # History data dump
+#dt        = 0.01  # time increment between outputs
+
+#<output2>
+#file_type = vtk   # Binary data dump
+#variable  = prim  # variables to be output
+#dt        = 1.0   # time increment between outputs
+
+#<output3>
+#file_type = rst      # Binary data dump
+#dt        = 0.00001  # time increment between outputs
+
+<output4>
+file_type = hdf5  # Tabular data dump
+variable  = prim  # variables to be output
+dt        = 1.0   # time increment between outputs
+
+<time>
+cfl_number = 0.4    # The Courant, Friedrichs, & Lewy (CFL) Number
+nlim       = -1     # cycle limit
+tlim       = 100.0  # time limit
+
+<mesh>
+nx1        = 1280  # Number of zones in X1-direction
+x1min      = -10   # minimum value of X1
+x1max      = 10    # maximum value of X1
+ix1_bc     = user  # inner-X1 boundary flag
+ox1_bc     = user  # inner-X1 boundary flag
+ix1_rad_bc = user
+ox1_rad_bc = user
+
+nx2        = 1         # Number of zones in X2-direction
+x2min      = -0.5      # minimum value of X2
+x2max      = 0.5       # maximum value of X2
+ix2_bc     = periodic  # inner-X2 boundary flag
+ox2_bc     = periodic  # inner-X2 boundary flag
+ix2_rad_bc = periodic
+ox2_rad_bc = periodic
+
+nx3        = 1         # Number of zones in X3-direction
+x3min      = -0.5      # minimum value of X3
+x3max      = 0.5       # maximum value of X3
+ix3_bc     = periodic  # inner-X3 boundary flag
+ox3_bc     = periodic  # inner-X3 boundary flag
+ix3_rad_bc = periodic
+ox3_rad_bc = periodic
+
+<meshblock>
+
+nx1 = 1280
+nx2 = 1
+nx3 = 1 
+
+<hydro>
+gamma = 1.6666666666667  # gamma = C_p/C_v
+
+<radiation>
+nmu         = 1
+Prat        = 0.0
+Crat        = 1000.0
+error_limit = 1.e-10
+iteration   = 2
+taucell     = 2
+
+<problem>