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

[felker]

Just opening this draft PR to start parsing the massive amount of changes and create an area for discussion.

• [x] Would definitely need some Wiki documentation before merging. At a minimum:
  • [x] New User Guide page "Newtonian radiation"
  • [x] https://github.com/PrincetonUniversity/athena/wiki/The-Input-File
  • [x] https://github.com/PrincetonUniversity/athena/wiki/Configuring
• [x] Consider excising separating particles, BlockFFT, thermal conduction from this PR? https://github.com/PrincetonUniversity/athena/tree/rad_newtonian/src/thermal_conduction
• [x] tst/multi_group/ needs to be reorganized and cleaned up.
  • [x] This would be the first time we are considering adding .ipynb notebooks to version control; do we want to do that?
  • [x] Remove the many duplicates of athena_read.py

[yanfeij]

@felker How to add equations in the documentation easily.

[munan]

I am working on merging the chemistry_scalar branch, which should be done in the next ~2 weeks before the Athena++ conference. There is also a radiation class there, although it's mostly just a container. See src/radiation on the branch.

[c-white]

@felker How to add equations in the documentation easily.

@yanfeij I just learned that as of very recently, github brought back mathjax support! So you can enclose basic latex commands in $...$ or $$...$$ and it should work.

[yanfeij]

I am working on merging the chemistry_scalar branch, which should be done in the next ~2 weeks before the Athena++ conference. There is also a radiation class there, although it's mostly just a container. See src/radiation on the branch.

I am afraid it will have exactly the same name as what I am using. Looks like you are just following six fixed rays, which can be a very special case of general radiation transport.

[yanfeij]

@felker How to add equations in the documentation easily.

@yanfeij I just learned that as of very recently, github brought back mathjax support! So you can enclose basic latex commands in $...$ or $$...$$ and it should work.

It works! Thank you.

[munan]

@yanfeij Yes, the six-ray integrator is a very special case. I am hoping the container Radiation class is general enough to be shared among different radiation modules. Perhaps we should discuss this. See src/radiation/radiation.hpp (or.cpp) in chemistry_scalar branch.

[yanfeij]

@yanfeij Yes, the six-ray integrator is a very special case. I am hoping the container Radiation class is general enough to be shared among different radiation modules. Perhaps we should discuss this. See src/radiation/radiation.hpp (or.cpp) in chemistry_scalar branch.

There are many more steps to do for general radiation transport. The implicit radiation transport is also using a completely different integrator. The data array structure and boundary values I need are also different. You can also take a look at the rad_newtonian branch to get an idea.

[c-white]

Just as a heads up, after this is merged I will work on merging the gr_rad branch as well. This will probably be more difficult, as there is a lot of overlap in files and functions but not details of what the code is doing.

Maybe we should decide among the following cases soon:

1. All radiation -- whether written by @munan or @yanfeij or me or anyone else) will eventually live in the same radiation module (i.e., the same class Radiation and the same directory src/radiation/). This might make for a confusingly large module.
2. Each radiation implementation will be separate, unless they already coexist (maybe some of @yanfeij 's do). There might be some duplicated code, but at least those interested in one use case would not have to sift through code pertaining to vastly different cases.
3. Something in between.

[yanfeij]

Just as a heads up, after this is merged I will work on merging the gr_rad branch as well. This will probably be more difficult, as there is a lot of overlap in files and functions but not details of what the code is doing.

Maybe we should decide among the following cases soon:

1. All radiation -- whether written by @munan or @yanfeij or me or anyone else) will eventually live in the same radiation module (i.e., the same class Radiation and the same directory src/radiation/). This might make for a confusingly large module.
2. Each radiation implementation will be separate, unless they already coexist (maybe some of @yanfeij 's do). There might be some duplicated code, but at least those interested in one use case would not have to sift through code pertaining to vastly different cases.
3. Something in between.

The two modules I have share some common data structures. But I think it will be a lot of work to make all these "radiation" to use the same integrator or data arrays. Also we do not expect to turn on any two of these all at once. I think it is easier to keep each one separate. But do give them some specific names, not just radiation. For example, six_ray_tracing, gr_radiation or something like that.

[felker]

@jmstone should we schedule a Zoom call in advance of meeting next month? Or just discuss here asynchronously? Should we absolutely wait to merge until the Athena dev meeting?

[ecostriker]

@yanfeij I agree that it's a good idea to decide on a strategy. FYI, @jeonggyukim is also starting work on implementing his ray-tracing radiation, so it will be good to include him in the discussion.

[jmstone]

Glad to see radiation is being prepared for merging! I certainly don't think it needs to wait until the dev meeting (although that would be a convenient venue for discussions). Scheduling a zoom between interested parties before then sounds like a fine idea -- however I'll leave it to others to take the lead on that as needed.

[yanfeij]

@felker I added wiki pages, replaced all the unnecessary files to the master version. Now "only" 135 files changed. Most of the files are new ones that I need. Can you take a look again, particularly the files also used by other modules?

[jeonggyukim]

Just as a heads up, after this is merged I will work on merging the gr_rad branch as well. This will probably be more difficult, as there is a lot of overlap in files and functions but not details of what the code is doing. Maybe we should decide among the following cases soon:

1. All radiation -- whether written by @munan or @yanfeij or me or anyone else) will eventually live in the same radiation module (i.e., the same class Radiation and the same directory src/radiation/). This might make for a confusingly large module.
2. Each radiation implementation will be separate, unless they already coexist (maybe some of @yanfeij 's do). There might be some duplicated code, but at least those interested in one use case would not have to sift through code pertaining to vastly different cases.
3. Something in between.

The two modules I have share some common data structures. But I think it will be a lot of work to make all these "radiation" to use the same integrator or data arrays. Also we do not expect to turn on any two of these all at once. I think it is easier to keep each one separate. But do give them some specific names, not just radiation. For example, six_ray_tracing, gr_radiation or something like that.

I agree with @yanfeij that option 2 is the way to go. Both six ray tracing and adaptive ray tracing are mainly used for UV photochemistry and share some common features with each other. I just started working on adaptive ray tracing and I will need to talk with @munan about a strategy. Until we need simulations that require a hybrid radiation scheme, however, it would be fine to have @yanfeij 's radiation transport module implemented separately.

[munan]

@yanfeij It seems that we share the energy density of radiation (variable "ir") but not much else. I agree having separate names for different radiation modules is an easier option now. Maybe I can coordinate with @jeonggyukim to do "chem_rad" focused on the interaction between chemistry and radiation.

Let's try to set up a meeting here to discuss radiation in Athena++: https://www.when2meet.com/?19752096-Ml59s Everyone interested is welcome to join @jmstone @c-white @felker @ecostriker

[yanfeij]

energy density of radiation

I will not call "ir" as "energy density of radiation" :)

I just filled up the when2meeting.

[ecostriker]

@yanfeij It seems that we share the energy density of radiation (variable "ir") but not much else. I agree having separate names for different radiation modules is an easier option now. Maybe I can coordinate with @jeonggyukim to do "chem_rad" focused on the interaction between chemistry and radiation.

Let's try to set up a meeting here to discuss radiation in Athena++: https://www.when2meet.com/?19752096-Ml59s Everyone interested is welcome to join @jmstone @c-white @felker @ecostriker

@jeonggyukim It would be good if you can join -- do we need to readjust times earlier?

[jeonggyukim]

Sure, I can join. No need to adjust times.

[munan]

Let's do the zoom meeting for radiation on next Monday, April 24, 9am Princeton/3pm Germany/10pm Japan. Zoom link here: https://zoom.us/j/3476901962?pwd=STBQVHNIcDlPdzdyUzdyQnNwM1Qzdz09 I will prepare an agenda and send it to you on Monday. Looking forward to it! @yanfeij @jeonggyukim @jmstone @felker @c-white @ecostriker

[munan]

Reminder for meeting today in ~3h :) Agenda here: https://comfortable-lace-bd0.notion.site/Radiation-in-Athena-meeting-f2f8ab6e774a444da7b813a4c697453e @yanfeij @jeonggyukim @jmstone @felker @c-white @ecostriker

[munan]

And zoom link again: https://zoom.us/j/3476901962?pwd=STBQVHNIcDlPdzdyUzdyQnNwM1Qzdz09 @yanfeij @jeonggyukim @jmstone @felker @c-white @ecostriker

[felker]

Sorry I missed this! I got my timezones confused. Did anyone take notes?

EDIT: just saw the Notion link, thanks!

[munan]

Sorry I missed this! I got my timezones confused. Did anyone take notes?

EDIT: just saw the Notion link, thanks!

Yes everything is there. We are going with option 1 :)

[yanfeij]

I have changed my radiation module to nr_radiation, which represents non-relativistic radiation transport. All the files and wiki have been updated. Regression tests passed. It will be useful if the admin team can take a look. @felker @c-white @tomidakn

[felker]

I will probably get around to this next weekend FYI

[c-white]

I too will have more time to look this over once a certain workshop is over...

[changgoo]

Large python style errors. To check those errors, run this

pip install flake8
python -m flake8

Maybe autopep8 will fix most issues automatically.

[felker]

this PR has a ton of issues with style that would make the code impossible to maintain if we merged it as-is. I have been slowly fixing them by hand and script. @yanfeij any chance you can take a pass or two? A big help would be the following issues:

• Remove duplicated boiler plate code, like descriptions of "beam.cpp" in many problem generators where it shouldnt be, and entirely duplicated comments like:

  // The information stored in the array
  // b_angle is
  // b_angle[0]=sin_theta_b
  // b_angle[1]=cos_theta_b
  // b_angle[2]=sin_phi_b
  // b_angle[3]=cos_phi_b

• Replace Numerical Recipes code with open source license code

Currently:

Category 'build/header_guard' errors found: 14
Category 'build/include_alpha' errors found: 139
Category 'build/include_order' errors found: 8
Category 'build/include_what_you_use' errors found: 2
Category 'readability/braces' errors found: 22
Category 'readability/casting' errors found: 3
Category 'readability/namespace' errors found: 2
Category 'runtime/casting' errors found: 18
Category 'runtime/explicit' errors found: 3
Category 'whitespace/blank_line' errors found: 502
Category 'whitespace/braces' errors found: 1248
Category 'whitespace/end_of_line' errors found: 907
Category 'whitespace/indent' errors found: 100
Category 'whitespace/line_length' errors found: 37
Category 'whitespace/newline' errors found: 21
Category 'whitespace/semicolon' errors found: 15
Category 'whitespace/tab' errors found: 2
Total errors found: 3043

[felker]

@yanfeij I would suggest adding a setting to your editor to delete trailing whitespaces when saving a file, in the future

[felker]

I think we should remove and avoid any closing loop comments like

} end k,j,i
}// end radiation

Unless the nested loop is particularly deep, long, or confusing.

[yanfeij]

@yanfeij I would suggest adding a setting to your editor to delete trailing whitespaces when saving a file, in the future

@felker Thank you. I changed the setting in sublime I used.

[yanfeij]

I think we should remove and avoid any closing loop comments like

} end k,j,i
}// end radiation

Unless the nested loop is particularly deep, long, or confusing.

OK

[felker]

In outflow_rad.cpp, the docstrings for the functions are mixed; some say "VACUUM" BCs, some say "REFLECTING".

Also it says, for example:

\fn void RadBoundaryVariable::OutflowInnerX2(

for

void RadBoundaryVariable::OutflowOuterX3(

Can you delete the wrong docstrings, including the ones at the top of copy/pasted files? E.g.

 //! \file reflect.cpp
 //  \brief implementation of reflecting BCs in each dimension

[yanfeij]

In outflow_rad.cpp, the docstrings for the functions are mixed; some say "VACUUM" BCs, some say "REFLECTING".

Also it says, for example:

\fn void RadBoundaryVariable::OutflowInnerX2(

for

void RadBoundaryVariable::OutflowOuterX3(

Can you delete the wrong docstrings, including the ones at the top of copy/pasted files? E.g.

 //! \file reflect.cpp
 //  \brief implementation of reflecting BCs in each dimension

I just committed fixes of these docstrings.

[felker]

What is going on In plm_simple.cpp? Did you mean to duplicate the X1, X2, X3 function definitions?

❯ grep -nri "void R" plm_simple.cpp
plm_simple.cpp:35:void Reconstruction::PiecewiseLinearX1(
plm_simple.cpp:110:void Reconstruction::PiecewiseLinearX1(
plm_simple.cpp:209:void Reconstruction::PiecewiseLinearX2(
plm_simple.cpp:282:void Reconstruction::PiecewiseLinearX2(
plm_simple.cpp:383:void Reconstruction::PiecewiseLinearX3(
plm_simple.cpp:451:void Reconstruction::PiecewiseLinearX3(
plm_simple.cpp:545:void Reconstruction::PiecewiseLinearZeta(
plm_simple.cpp:596:void Reconstruction::PiecewiseLinearPsi(

Where is ALI_LEN defined?

#pragma omp simd simdlen(SIMD_WIDTH) aligned(qln,qrn,qcn,dqmn:ALI_LEN)

[yanfeij]

What is going on In plm_simple.cpp? Did you mean to duplicate the X1, X2, X3 function definitions?

❯ grep -nri "void R" plm_simple.cpp
plm_simple.cpp:35:void Reconstruction::PiecewiseLinearX1(
plm_simple.cpp:110:void Reconstruction::PiecewiseLinearX1(
plm_simple.cpp:209:void Reconstruction::PiecewiseLinearX2(
plm_simple.cpp:282:void Reconstruction::PiecewiseLinearX2(
plm_simple.cpp:383:void Reconstruction::PiecewiseLinearX3(
plm_simple.cpp:451:void Reconstruction::PiecewiseLinearX3(
plm_simple.cpp:545:void Reconstruction::PiecewiseLinearZeta(
plm_simple.cpp:596:void Reconstruction::PiecewiseLinearPsi(

Where is ALI_LEN defined?

#pragma omp simd simdlen(SIMD_WIDTH) aligned(qln,qrn,qcn,dqmn:ALI_LEN)

@felker Radiation variables (specific intensities) cannot use the default plm_simple functions, because ir takes the order[k,j,i,n] while the default the Athena order is [n,k,j,i]. I created overwritten functions to take array with different order so that I do not need to modify the default plm_simple functions.

ALI_LEN is removed. I will clean that up.

[yanfeij]

What is going on In plm_simple.cpp? Did you mean to duplicate the X1, X2, X3 function definitions?

❯ grep -nri "void R" plm_simple.cpp
plm_simple.cpp:35:void Reconstruction::PiecewiseLinearX1(
plm_simple.cpp:110:void Reconstruction::PiecewiseLinearX1(
plm_simple.cpp:209:void Reconstruction::PiecewiseLinearX2(
plm_simple.cpp:282:void Reconstruction::PiecewiseLinearX2(
plm_simple.cpp:383:void Reconstruction::PiecewiseLinearX3(
plm_simple.cpp:451:void Reconstruction::PiecewiseLinearX3(
plm_simple.cpp:545:void Reconstruction::PiecewiseLinearZeta(
plm_simple.cpp:596:void Reconstruction::PiecewiseLinearPsi(

Where is ALI_LEN defined?

#pragma omp simd simdlen(SIMD_WIDTH) aligned(qln,qrn,qcn,dqmn:ALI_LEN)

@felker Radiation variables (specific intensities) cannot use the default plm_simple functions, because ir takes the order[k,j,i,n] while the default the Athena order is [n,k,j,i]. I created overwritten functions to take array with different order so that I do not need to modify the default plm_simple functions.

ALI_LEN is removed. I will clean that up.

By the way, similar functions are also created for dc_simple

[felker]

Cool, thanks.

I dont understand these lines in reflect_rad.cpp

          Real *iri = &((*var_cc)(k,j,(iu-i+1),ifr*pmy_block_->pnrrad->nang));
          Real *iro = &((*var_cc)(k,j, iu+i, ifr*pmy_block_->pnrrad->nang));

They throw a cpplint warning:

Are you taking an address of a cast?  This is dangerous: could be a temp var.  Take the address before doing the \
cast, rather than after  [runtime/casting] [4]

[yanfeij]

Cool, thanks.

I dont understand these lines in reflect_rad.cpp

          Real *iri = &((*var_cc)(k,j,(iu-i+1),ifr*pmy_block_->pnrrad->nang));
          Real *iro = &((*var_cc)(k,j, iu+i, ifr*pmy_block_->pnrrad->nang));

They throw a cpplint warning:

Are you taking an address of a cast?  This is dangerous: could be a temp var.  Take the address before doing the \
cast, rather than after  [runtime/casting] [4]

So the last dimension has the size nfreq*nang. For the boundary condition, I need to distinguish all the angles in each frequency group. So I try to access the memory address at the beginning of each frequency group. I do not understand the recommendation. Shall we take the address first and move the pointer by some amount?

[felker]

Category 'runtime/casting' errors found: 18
Category 'whitespace/blank_line' errors found: 58
Category 'whitespace/braces' errors found: 98
Category 'whitespace/end_of_line' errors found: 1
Category 'whitespace/indent' errors found: 17
Category 'whitespace/line_length' errors found: 9
Total errors found: 201

progress ...

[yanfeij]

Cool, thanks.

I dont understand these lines in reflect_rad.cpp

          Real *iri = &((*var_cc)(k,j,(iu-i+1),ifr*pmy_block_->pnrrad->nang));
          Real *iro = &((*var_cc)(k,j, iu+i, ifr*pmy_block_->pnrrad->nang));

They throw a cpplint warning:

Are you taking an address of a cast?  This is dangerous: could be a temp var.  Take the address before doing the \
cast, rather than after  [runtime/casting] [4]

So the last dimension has the size nfreq*nang. For the boundary condition, I need to distinguish all the angles in each frequency group. So I try to access the memory address at the beginning of each frequency group. I do not understand the recommendation. Shall we take the address first and move the pointer by some amount?

[yanfeij]

Is there any automatic script to fix these?

[felker]

working on understanding that runtime/casting error.

Why did you close the issue?

[yanfeij]

working on understanding that runtime/casting error.

Why did you close the issue?

It is mistake. I tried to close a comment. Didn't realize it closed the whole thing.

[felker]

Another reminder: it is best to use std::pow, std::exp, std::log, std::abs instead of pow, exp, log, abs, fabs.

[felker]

First pass of style fixes complete. Still should fix a lot of comments, delete unnecessary blank lines.

Now let's see which regression tests fail.

runtime/casting was just getting confused by the syntax of Real *iri = &((*var_cc)(k,j,(iu-i+1),ifr*pmy_block_->pnrrad->nang)); since (*var_cc*) looks like a type cast. I avoided it by creating a temporary reference variable (hopefully is still correct, @yanfeij can you check the last commit?)

[yanfeij]

First pass of style fixes complete. Still should fix a lot of comments, delete unnecessary blank lines.

Now let's see which regression tests fail.

runtime/casting was just getting confused by the syntax of Real *iri = &((*var_cc)(k,j,(iu-i+1),ifr*pmy_block_->pnrrad->nang)); since (*var_cc*) looks like a type cast. I avoided it by creating a temporary reference variable (hopefully is still correct, @yanfeij can you check the last commit?)

It looks correct to me.

[felker]

Previously, these convenience macros are (always?) translated to double literals:

#define TINY_NUMBER 1.0e-20
#define HUGE_NUMBER 1.0e+36

TINY_NUMBER only appears 27 times in master, and within only a couple std:: math function calls, like:

twid_asq = std::max((gm1*(hp-0.5*vsq)-(gm1-1.0)*x), TINY_NUMBER);

These were fine in all cases, even though vsq, etc. are Real and therefore float when the code is configured with -float, because the other floating point literals like 1.0 are double so the first argument in the expression gets converted to double. This is the one exception, since cost is a double:

cost_ = std::min(cost, TINY_NUMBER);

(@tomidakn should cost always be double precision, even if the code was configured with -float ?)

But in this radiation branch, the macro appears in many std::max, min calls, with the first argument being a bare Real variable:

irn[n] = std::max(irn[n], TINY_NUMBER);

which causes a compile/link error when configuring the code with -float. The solution is to change:

#define TINY_NUMBER Real(1.0e-20)
#define HUGE_NUMBER Real(1.0e+36)

which is what @yanfeij did, but then we would need to change

cost_ = std::min(cost, static_cast<double>(TINY_NUMBER));

Or instead, change many radiation calls like:

irn[n] = std::max(irn[n], static_cast<Real>(TINY_NUMBER));

Thoughts?

[yanfeij]

Previously, these convenience macros are (always?) translated to double literals:

#define TINY_NUMBER 1.0e-20
#define HUGE_NUMBER 1.0e+36

TINY_NUMBER only appears 27 times in master, and within only a couple std:: math function calls, like:

twid_asq = std::max((gm1*(hp-0.5*vsq)-(gm1-1.0)*x), TINY_NUMBER);

These were fine in all cases, even though vsq, etc. are Real and therefore float when the code is configured with -float, because the other floating point literals like 1.0 are double so the first argument in the expression gets converted to double. This is the one exception, since cost is a double:

cost_ = std::min(cost, TINY_NUMBER);

(@tomidakn should cost always be double precision, even if the code was configured with -float ?)

But in this radiation branch, the macro appears in many std::max, min calls, with the first argument being a bare Real variable:

irn[n] = std::max(irn[n], TINY_NUMBER);

which causes a compile/link error when configuring the code with -float. The solution is to change:

#define TINY_NUMBER Real(1.0e-20)
#define HUGE_NUMBER Real(1.0e+36)

which is what @yanfeij did, but then we would need to change

cost_ = std::min(cost, static_cast<double>(TINY_NUMBER));

Or instead, change many radiation calls like:

irn[n] = std::max(irn[n], static_cast<Real>(TINY_NUMBER));

Thoughts?

I will be happy to change just the radiation files and keep others unchanged.

[yanfeij]

@felker MPI regression test seems broken in the rad_newtonian branch. When the code compiles, it now complains:

src/mesh/amr_loadbalance.cpp:448:16: error: use of undeclared identifier 'bnx1' int bssame = bnx1bnx2bnx3*nx4_tot;