Add chemistry module built on passive scalars module

[munan]

Merge chemistry module to the master branch

Prerequisite checklist

• My code follows the Athena++ Style Guide
• My change requires a change to the documentation.
• I have updated the documentation in the Wiki accordingly.
• I have added tests to cover my changes.
• All new and existing tests passed.

Please review the CONTRIBUTING.md file for detailed guidelines.

Description

The detailed description of the chemistry module can be found at the wiki page: https://github.com/PrincetonUniversity/athena/wiki/chemistry
An accompanying code paper will be submitted shortly and linked to the wiki.

Testing and validation

Added regression tests under "chemistry"

[felker]

Yes, CR and explicit RT are in the main integrator, while the implicit RT is completely separate. That is because implicit RT needs to work on the whole mesh and iterative for several times. So we just finish all the tasks in the main integrator and then perform the implicit RT task list

In this case, is the only coupling to the fluid via RadIntegrator::AddSourceTerms() called within pimradhylist->DoTaskListOneStage(wght); in rad_iteration.cpp? Is it not a numerical issue that it occurs long after TimeIntegratorTaskList::AddSourceTerms() from the main explicit fluid integrator?

Sorry if these questions are answered in your latest implicit rad paper, I dont have time to read it right now :)

[yanfeij]

Yes, CR and explicit RT are in the main integrator, while the implicit RT is completely separate. That is because implicit RT needs to work on the whole mesh and iterative for several times. So we just finish all the tasks in the main integrator and then perform the implicit RT task list

In this case, is the only coupling to the fluid via RadIntegrator::AddSourceTerms() called within pimradhylist->DoTaskListOneStage(wght); in rad_iteration.cpp? Is it not a numerical issue that it occurs long after TimeIntegratorTaskList::AddSourceTerms() from the main explicit fluid integrator?

Sorry if these questions are answered in your latest implicit rad paper, I dont have time to read it right now :)

What are the issues you have in mind? It is still VL integrator for example. At each stage, you add all the changes to U. It is not necessary to add all the source terms at the same location. It also shouldn't matter the order you add the source terms. The radiation source terms are added during both the predict and correct steps.

[felker]

I dont have any particular issue in mind. Just trying to mentally map out the extension of @c-white's diagram with implicit radiation: https://github.com/PrincetonUniversity/athena/wiki/presentations/athena_structure_mhd_refinement.pdf
(unless he wants to volunteer to make a newer complete version from scratch 😄 )

I also dont clearly understand the separation between:

• RadIntegrator::CalSourceTerms
• RadIntegrator::AddSourceTerms
• RadIntegrator::GetHydroSourceTerms

and which are used when in explicit vs. implicit radiation schemes.

I suppose it is fine to add any radiation fluid source terms in any order (and in any separate task list after the main TimeIntegratorTaskList), since HydroSourceTerms() is called after HydroIntegrate() in each explicit fluid integrator stage--- but then dont you need to re-call the CONS2PRIM after changing the conserved variables?

Any separate task list also wouldnt work correctly if the implicit radiation task list were called before TimeIntegratorTaskList, since CalculateFluxes() would be based on the CONS2PRIM/Primitives()/W(U) from the previous timestepper stage?

Is the scheme fully compatible with the other explicit time-integrators for the fluid, i.e. RK3, RK4?

---

Also, we should probably change ADD_RAD_SCR to ADD_RAD_SRC and similarly for FLX_AND_SCR in im_radhydro_task_list.cpp, im_radit_task_list.cpp, im_rad_task_list.hpp, assuming it is supposed to be "source" and not "source cosmic rays" or something.

Finally, I just noticed that TaskStatus AddSourceTermsIMRad(MeshBlock *pmb, int stage); appears in task_list.hpp and isnt defined or used anywhere. Are there other such cases?

[changgoo]

retest this

[yanfeij]

I dont have any particular issue in mind. Just trying to mentally map out the extension of @c-white's diagram with implicit radiation: https://github.com/PrincetonUniversity/athena/wiki/presentations/athena_structure_mhd_refinement.pdf (unless he wants to volunteer to make a newer complete version from scratch 😄 )

I also dont clearly understand the separation between:

• RadIntegrator::CalSourceTerms
• RadIntegrator::AddSourceTerms
• RadIntegrator::GetHydroSourceTerms

and which are used when in explicit vs. implicit radiation schemes.

I suppose it is fine to add any radiation fluid source terms in any order (and in any separate task list after the main TimeIntegratorTaskList), since HydroSourceTerms() is called after HydroIntegrate() in each explicit fluid integrator stage--- but then dont you need to re-call the CONS2PRIM after changing the conserved variables?

Any separate task list also wouldnt work correctly if the implicit radiation task list were called before TimeIntegratorTaskList, since CalculateFluxes() would be based on the CONS2PRIM/Primitives()/W(U) from the previous timestepper stage?

Is the scheme fully compatible with the other explicit time-integrators for the fluid, i.e. RK3, RK4?

Also, we should probably change ADD_RAD_SCR to ADD_RAD_SRC and similarly for FLX_AND_SCR in im_radhydro_task_list.cpp, im_radit_task_list.cpp, im_rad_task_list.hpp, assuming it is supposed to be "source" and not "source cosmic rays" or something.

Finally, I just noticed that TaskStatus AddSourceTermsIMRad(MeshBlock *pmb, int stage); appears in task_list.hpp and isnt defined or used anywhere. Are there other such cases?

so, no, radiation is not used with RK3 and RK4. And there is no need to do it since implicit RT is first order anyway. I think anytime you update conservative variable, you just need to perform another CONS2PRIM, then it is fine.

[felker]

I want to revisit this and merge it soon, in the next few weeks. I think there was one or two outstanding design issues, but I am having trouble finding where I summarized them... #496 (comment) seems to be one of my more recent summary comments, but that was on July 7th, before we merged radiation on July 16th and then resolved conflicts + got all the tests passing--- so I think that is stale.

Might just be the 6x leftover TODO from @tomidakn's reviews on July 6th and 9th: https://github.com/PrincetonUniversity/athena/pull/496/files#r1254336650
But nothing pressing. We could then mint a new tag and release version, v23.0?

Any thoughts @changgoo @yanfeij @c-white @jmstone (@munan, only if you are back)?

[felker]

bump--- any thoughts?

[munan]

I, too, am getting extremely confused with the many task lists and how they interact with each other at this point. After this PR, we would have the following classes (and their parent-child relationships):

• TaskList

  • TimeIntegratorTaskList: allocated and called in main.cpp
  • ChemRadiationIntegratorTaskList: allocated and called in main.cpp (before main task list). Coupled to TimeIntegratorTaskList via IntegrateChemistry() task only on the last stage of the multistage time-integrator? I dont fully understand this, can you clarify @munan ?

The IntegrateChemistry() solves the source term of the chemistry ODE, and it is added in the TimeIntegratorTaskList after the advection of passive scalars. The ChemRadiationIntegratorTaskList solves the radiation that is related to chemistry, and yes it is called in main.cpp before the time integrator task list.

• SuperTimeStepTaskList: allocated and called in main.cpp, before and after the main TimeIntegratorTaskList task list, depending on the operator splitting scheme. Coupling through INT_HYD etc. tasks in SuperTimeStepTaskList (accurate @pdmullen ?)

• GravityBoundaryTaskList: allocated and called in either fft_gravity.cpp or mg_gravity.cpp. I assume this is basically coupled only through TimeIntegratorTaskList::AddSourceTerms()

• MultigridTaskList: allocated in mg_gravity.cpp called in multigrid_driver.cpp

• IMRadTaskList: all 3x derived classes are allocated in radiation_implicit.cpp and called in rad_iteration.cpp.

  • IMRadITTaskList
  • IMRadComptTaskList
  • IMRadHydroTaskList

Are both CR transport and explicit radiation just fully integrated in the main TimeIntegratorTaskList? Is my understanding correct @yanfeij? How does implicit radiation couple to the main task list?

The many TaskID tags that were previously "reserved" for chemistry (by being commented-out) were unused in this PR, so I removed the comments entirely. I believe this is because the original plan was to have all the chemistry tasks intertwined in the original task list, but that was abandoned @munan?

Yes, You are right, they should be removed. I also found another old CALC_CHMFLX TaskID that is unused. I will remove that too.

Related #502.

[munan]

I want to revisit this and merge it soon, in the next few weeks. I think there was one or two outstanding design issues, but I am having trouble finding where I summarized them... #496 (comment) seems to be one of my more recent summary comments, but that was on July 7th, before we merged radiation on July 16th and then resolved conflicts + got all the tests passing--- so I think that is stale.

Might just be the 6x leftover TODO from @tomidakn's reviews on July 6th and 9th: https://github.com/PrincetonUniversity/athena/pull/496/files#r1254336650 But nothing pressing. We could then mint a new tag and release version, v23.0?

Any thoughts @changgoo @yanfeij @c-white @jmstone (@munan, only if you are back)?

I am finally back part-time! I will put merging chemistry as my priority now. I went through your comments and here are what I can find that need to be done:

• Check all the TODOs
• change configure option —ode_solver to —chem_ode_solver
• delete obsolete task ID CALC_CHMFLX

I will start working on these. Anything else?

[munan]

OK, I just resolved the issues in the previous comments. The chemistry regression tests are passed but I haven't checked the other regression tests yet. @felker Are we ready to merge after checking all the regression tests?

[changgoo]

usage: configure.py [-h]
                    [--prob {hgb,field_loop,read_vtk,scalar_diff,jeans,fft,shu_osher,poisson,thermal_relaxation,collapse,chem_Pn_grid,chem_uniform_sixray,chem_pdr_moving,chem_G14Sod,rad_linearwave,resist,cr_diffusion,thermal_multigroup,twoibw,dmr,gr_blast,quirk,binary_gravity,turb,shock_tube,slotted_cylinder,from_array,gr_geodesic_infall,gr_linear_wave,hb3,blast,visc,disk,gr_mhd_inflow,chem_uniform,jet,rt,magnoh,gr_bondi,cpaw,shk_cloud,rotor,orszag_tang,lw_implode,linear_wave,chem_onecell_sixray,mignone_advection,ssheet,noh,gr_torus,strat,chem_turb,eos_test,chem_H2,field_loop_poles,kh,default_pgen,gr_shock_tube,jgg,beam}]
                    [--coord {cartesian,cylindrical,spherical_polar,minkowski,sinusoidal,tilted,schwarzschild,kerr-schild,gr_user}]
                    [--eos {adiabatic,isothermal,general/eos_table,general/hydrogen,general/ideal}]
                    [--flux {default,hlle,hllc,lhllc,hlld,lhlld,roe,llf}]
                    [--nghost NGHOST] [--nscalars NSCALARS]
                    [--nspecies NSPECIES] [-b] [-sts] [-s] [-g] [-t] [-debug]
                    [-coverage] [-float] [-mpi] [-omp] [--grav {none,fft,mg}]
                    [-fft] [--fftw_path FFTW_PATH]
                    [--chemistry {gow17,H2,kida,G14Sod}]
                    [--kida_rates {H2,gow17,nitrogen,nitrogen_gas_Sipila}]
                    [--chem_radiation {const,six_ray}]
                    [--chem_ode_solver {cvode,forward_euler}]
                    [--cvode_path CVODE_PATH] [-hdf5] [-h5double]
                    [--hdf5_path HDF5_PATH] [-nr_radiation]
                    [-implicit_radiation] [-cr]
                    [--cxx {g++,g++-simd,icpx,icpc,icpc-debug,icpc-phi,cray,bgxlc++,clang++,clang++-simd,clang++-apple}]
                    [--ccmd CCMD] [--mpiccmd MPICCMD] [--gcovcmd GCOVCMD]
                    [--cflag CFLAG] [--include INCLUDE] [--lib_path LIB_PATH]
                    [--lib LIB]
configure.py: error: unrecognized arguments: --ode_solver=forward_euler
Exception occurred
Traceback (most recent call last):
  File "/scratch/gpfs/changgoo/jenkins/workspace/tigris/athena-PR/tst/regression/scripts/utils/athena.py", line 41, in configure
    subprocess.check_call(configure_command + global_config_args,
  File "/usr/licensed/anaconda3/2023.3/lib/python3.10/subprocess.py", line 369, in check_call
    raise CalledProcessError(retcode, cmd)
subprocess.CalledProcessError: Command '['/usr/licensed/anaconda3/2023.3/bin/python', 'configure.py', '--chemistry=H2', '--ode_solver=forward_euler', '-float', '--prob=chem_pdr_moving', '--cxx=icpx', '--cflag=-Wall -Wextra -Werror -pedantic -pedantic-errors -Wno-unused-private-field -Wno-unused-variable -Wno-unused-parameter -Wunknown-pragmas -Wno-unused-function -Wshorten-64-to-32 -Wno-array-bounds -Wno-pass-failed -fp-model=precise']' returned non-zero exit status 2.

Test failed with this message.

[changgoo]

It looks like you've changed ode_solver option to chem_ode_solver

[munan]

It looks like you've changed ode_solver option to chem_ode_solver

I just fixed this and the pgen_compile regression test should pass now. Could you try again? Can I run regression test on jenkins or do you have to do it?

[changgoo]

Any new commit triggers new tests. This time, it failed with some unknown reason. I'm gonna retest this on stellar manually.

[changgoo]

retest this on stellar

[tomidakn]

I'm sorry for being late as I've been very busy these days. And (belated) welcome back, @munan !

Although I only quickly skimmed the code, it looks good to me as long as they pass the tests. I like the use of TaskList for Six-ray radiation transport. It is nice if we could merge it during the winter break.

[changgoo]

retest this on stellar