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 non-relativistic radiation transport
+#   -implicit_radiation  implicit radiation transport module
+#   -cr                  enable cosmic ray transport
 # ----------------------------------------------------------------------------------------
 
 # Modules
@@ -205,6 +208,24 @@
                     default='',
                     help='path to HDF5 libraries')
 
+# -nr_radiation argument
+parser.add_argument('-nr_radiation',
+                    action='store_true',
+                    default=False,
+                    help='enable non-relativistic radiative transfer')
+
+# -implicit_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
@@ -350,6 +371,14 @@ def c_to_cpp(arg):
         raise SystemExit('### CONFIGURE ERROR: '
                          + 'General EOS is incompatible with flux ' + args['flux'])
 
+if args['g'] and (args['nr_radiation'] or args['implicit_radiation']):
+    raise SystemExit('### CONFIGURE ERROR: '
+                     + ' GR is incompatible with nr_radiation or implicit_radiation')
+
+if args['nr_radiation'] and args['implicit_radiation']:
+    raise SystemExit('### CONFIGURE ERROR: '
+                     + ' nr_radiation and implicit_radiation cannot be used together')
+
 # --- Step 3. Set definitions and Makefile options based on above argument
 
 # Prepare dictionaries of substitutions to be made
@@ -444,6 +473,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)
@@ -811,30 +865,52 @@ def c_to_cpp(arg):
 elif args['grav'] == 'mg':
     self_grav_string = 'Multigrid'
 
-print('Your Athena++ distribution has now been configured with the following options:')
-print('  Problem generator:          ' + args['prob'])
-print('  Coordinate system:          ' + args['coord'])
-print('  Equation of state:          ' + args['eos'])
-print('  Riemann solver:             ' + args['flux'])
-print('  Magnetic fields:            ' + ('ON' if args['b'] else 'OFF'))
-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('  Frame transformations:      ' + ('ON' if args['t'] else 'OFF'))
-print('  Self-Gravity:               ' + self_grav_string)
-print('  Super-Time-Stepping:        ' + ('ON' if args['sts'] else 'OFF'))
-print('  Debug flags:                ' + ('ON' if args['debug'] else 'OFF'))
-print('  Code coverage flags:        ' + ('ON' if args['coverage'] else 'OFF'))
-print('  Linker flags:               ' + makefile_options['LINKER_FLAGS'] + ' '
-      + makefile_options['LIBRARY_FLAGS'])
-print('  Floating-point precision:   ' + ('single' if args['float'] else 'double'))
-print('  Number of ghost cells:      ' + args['nghost'])
-print('  MPI parallelism:            ' + ('ON' if args['mpi'] else 'OFF'))
-print('  OpenMP parallelism:         ' + ('ON' if args['omp'] else 'OFF'))
-print('  FFT:                        ' + ('ON' if args['fft'] else 'OFF'))
-print('  HDF5 output:                ' + ('ON' if args['hdf5'] else 'OFF'))
+
+def output_config(opt_descr, opt_choice, filehandle=None):
+    first_col_width = 32
+    first_col_indent = 2
+    descr_len = len(opt_descr)
+    right_pad_len = first_col_width - (descr_len + first_col_indent + 2)  # include colon
+    right_pad = right_pad_len*' ' if right_pad_len >= 0 else ''
+    line_str = first_col_indent*' ' + opt_descr + ': ' + right_pad + opt_choice
+    print(line_str)
+    if (filehandle is not None):
+        filehandle.write(line_str + '\n')
+
+
+# write the configuration optitions into a log file
+flog = open('./configure.log', 'w')
+
+output_config('Your Athena++ distribution has now been configured with the following options', '', flog)  # noqa
+output_config('Problem generator', args['prob'], flog)
+output_config('Coordinate system', args['coord'], flog)
+output_config('Equation of state', args['eos'], flog)
+output_config('Riemann solver', args['flux'], flog)
+output_config('Magnetic fields', ('ON' if args['b'] else 'OFF'), flog)
+output_config('Number of scalars', args['nscalars'], flog)
+output_config('Special relativity', ('ON' if args['s'] else 'OFF'), flog)
+output_config('General relativity', ('ON' if args['g'] else 'OFF'), flog)
+output_config('Radiative Transfer', ('ON' if args['nr_radiation'] else 'OFF'), flog)
+output_config('Implicit Radiation', ('ON' if args['implicit_radiation'] else 'OFF'), flog)
+output_config('Cosmic Ray Transport', ('ON' if args['cr'] else 'OFF'), flog)
+output_config('Frame transformations', ('ON' if args['t'] else 'OFF'), flog)
+output_config('Self-Gravity', self_grav_string, flog)
+output_config('Super-Time-Stepping', ('ON' if args['sts'] else 'OFF'), flog)
+output_config('Debug flags', ('ON' if args['debug'] else 'OFF'), flog)
+output_config('Code coverage flags', ('ON' if args['coverage'] else 'OFF'), flog)
+output_config('Linker flags', makefile_options['LINKER_FLAGS'] + ' '
+              + makefile_options['LIBRARY_FLAGS'], flog)
+output_config('Floating-point precision', ('single' if args['float'] else 'double'), flog)
+output_config('Number of ghost cells', args['nghost'], flog)
+output_config('MPI parallelism', ('ON' if args['mpi'] else 'OFF'), flog)
+output_config('OpenMP parallelism', ('ON' if args['omp'] else 'OFF'), flog)
+output_config('FFT', ('ON' if args['fft'] else 'OFF'), flog)
+output_config('HDF5 output', ('ON' if args['hdf5'] else 'OFF'), flog)
 if args['hdf5']:
-    print('  HDF5 precision:             ' + ('double' if args['h5double'] else 'single'))
-print('  Compiler:                   ' + args['cxx'])
-print('  Compilation command:        ' + makefile_options['COMPILER_COMMAND'] + ' '
-      + makefile_options['PREPROCESSOR_FLAGS'] + ' ' + makefile_options['COMPILER_FLAGS'])
+    output_config('HDF5 precision', ('double' if args['h5double'] else 'single'), flog)
+output_config('Compiler', args['cxx'], flog)
+output_config('Compilation command', makefile_options['COMPILER_COMMAND'] + ' '
+              + makefile_options['PREPROCESSOR_FLAGS'] + ' '
+              + makefile_options['COMPILER_FLAGS'], flog)
+
+flog.close()

inputs/cosmic_ray/athinput.cr_diffusion

@@ -0,0 +1,54 @@
+<comment>
+problem   = cosmic ray diffusion
+reference =
+configure = --prob=cr_diffusion
+
+<job>
+problem_id = cr  # problem ID: basename of output filenames
+
+
+<time>
+cfl_number = 0.3  # The Courant, Friedrichs, & Lewy (CFL) Number
+nlim       = -1   # cycle limit
+tlim       = 0.2  # time limit
+ncycle_out = 1    # interval for stdout summary info
+
+<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
+
+
+<meshblock>
+nx1 = 16
+nx2 = 256
+nx3 = 1
+
+
+<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.rad_linearwave

@@ -0,0 +1,58 @@
+<comment>
+problem   = Radiation linear wave test
+reference =
+configure = --prob=rad_linearwave -nr_radiation
+or
+configure = --prob=rad_linearwave -implicit_radiation
+
+<job>
+problem_id = radwave  # problem ID: basename of output filenames
+
+<output1>
+file_type = hst   # History data dump
+dt        = 0.01  # time increment between outputs
+
+
+<time>
+cfl_number = 0.4                 # The Courant, Friedrichs, & Lewy (CFL) Number
+nlim       = 100000              # cycle limit
+tlim       = 0.7745966144169111  # time limit
+ncycle_out = 1                   # interval for stdout summary info
+
+<mesh>
+nx1        = 256       # Number of zones in X1-direction
+x1min      = 0.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        = 8         # 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
+
+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
+
+<meshblock>
+nx1 = 64
+nx2 = 8
+nx3 = 1
+
+<hydro>
+gamma = 1.6666666666667  # gamma = C_p/C_v
+
+<radiation>
+nmu         = 1
+prat        = 0.01
+crat        = 10
+error_limit = 1.e-12
+taucell     = 5
+
+<problem>
+regime        = 1
+compute_error = true