Fully-implicit cosmic-ray diffusion module using Multigrid

[tomidakn]

Here is my Christmas present.

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.

Description

This PR implements a fully-implicit cosmic-ray diffusion solver using the Multigrid method. For details, please see [[Cosmic‐ray diffusion with Multigrid]] in Wiki. This also includes some enhancements in the Multigrid solver. It is now generalized for complicated physics. Also, now users can set the Multigrid parameters at runtime and additional diagnosis information can be shown on the fly.

Note that this is quite different from Yan-Fei's [[Cosmic Ray Transport]] module. This module adopts the diffusion approximation solving the zeroth moment only, which is less accurate. On the other hand, as this is a fully-implicit method, it runs using the (long) timestep of the MHD part, without using the reduced speed of light technique.

Testing and validation

An anisotropic diffusion test problem is implemented, but more quantitative tests are necessary.

The capability is still limited, but I think it is already useful for some users, and I want users to help testing and improving it.

[tomidakn]

@changgoo , can you tell me why the regression tests failed?

[changgoo]

I'm going to forward the error message to you. Dr. Chang-Goo Kim Research Scholar Department of Astrophysical Sciences Princeton University http://changgoo.github.io/index.html

…

On Sun, Dec 24, 2023 at 8:37 AM Kengo TOMIDA ***@***.***> wrote: @changgoo <https://github.com/changgoo> , can you tell me why the regression tests failed? — Reply to this email directly, view it on GitHub <#550 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABNKZHZLZTQKTYVCU7M746DYLAVX5AVCNFSM6AAAAABBBJH6CSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNRYGUYTSMBVHA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[tomidakn]

Thanks! And it is still failing...

[felker]

Just merged chemistry in #496, so there are a bunch of conflicts with master now. I will review this PR after they are fixed.

And finally can mint a new release after this and #549 and #553 are merged?

[tomidakn]

I am happy to merge it, and I agree to tag it, 2023 or 2024.

[felker]

My biggest questions regarding this PR are about how the -crdiff and -cr capabilities interact or conflict?

Your edits to https://github.com/PrincetonUniversity/athena/wiki/Cosmic%E2%80%90ray-diffusion-with-Multigrid nicely contrast the two algorithms, but do we need to be worried about a user configuring Athena++ with both options? I cant imagine any reason why you would want to use both at the same time, so maybe we should have configure.py throw an error in such a case. On the other hand, you may want to run the same binary and switch between the two solvers at runtime for comparison--- but I don't think this is currently possible anyway.

Other comments:

• It is unfortunate that there is no consistency in variable names between the two CR algorithms. E.g. Dpara/Dperp vs. sigma_diff(0:2, ...), initial condition set via pcrdiff->ecr(k,j,i) vs. pcr->u_cr(4,k,j,i). Why does CR diffusion have Lambda, zeta_factor, but CR transport does not seem to have counterparts? Should they be added to CR transport?
• Some of the options aren't fully explained in https://github.com/PrincetonUniversity/athena/wiki/Cosmic%E2%80%90ray-diffusion-with-Multigrid. You might want to refer to the parent page and/or mention that MGI is the other mgmode possibility. Are the 4x parameters that are commented-out (npresmooth, ... , threshold) supposed to be commented-out? Is everyone required, or are there defaults?
• Add a small text subsection and link to CR Diffusion wiki page on https://github.com/PrincetonUniversity/athena/wiki/Cosmic-Ray-Transport
• Would be great if @yanfeij could briefly weigh-in, since he authored the other CR module.
• (not for this PR, necessarily) https://github.com/PrincetonUniversity/athena/wiki/Cosmic-Ray-Transport only describes cr/vmax and EnrollOpacityFunction(). Need to describe some other parameters and functions in the class, especially those related to CR streaming and source functions. Also would be good to include the counterpart of the "equations solved" copied from Jiang and Oh (2018), assuming that they are different from the CR diffusion equations with the inclusion of the velocity-related terms. This is related to making the names consistent between CR diffusion and transport modules. Some items from cr.hpp header:

  Real vlim;
  Real max_opacity;
  AthenaArray<Real> sigma_adv;  // KGF: show an example of how this is set in the enroll functions 
...
  // Function in problem generators to update opacity
...
  void EnrollStreamingFunction(CRStreamingFunc MyStreamingFunction);

  void EnrollUserCRSource(CRSrcTermFunc my_func);

• Should we rename pgen/cr_diffusion.cpp to pgen/cr_diffusion_transport.cpp to make it obvious and contrast it with the new pgen/cr_diffusion_mg.cpp?
• In general, it would be good to have more examples of CR transport pgens, and at least one that shows how to use EnrollUserCRSource(), etc.

[felker]

see above

[felker]

@tomidakn did you get a chance to look at this? With #560 merged, this is the last thing that would make it in before releasing 2024.1. Unless @munan wants us to wait for #553 ?

[tomidakn]

I am currently quite busy, but will work on it over this weekend (or next week).

[tomidakn]

My biggest questions regarding this PR are about how the -crdiff and -cr capabilities interact or conflict?

Currently, they are just independent modules to solve similar problems and do not share anything.

Your edits to https://github.com/PrincetonUniversity/athena/wiki/Cosmic%E2%80%90ray-diffusion-with-Multigrid nicely contrast the two algorithms, but do we need to be worried about a user configuring Athena++ with both options? I cant imagine any reason why you would want to use both at the same time, so maybe we should have configure.py throw an error in such a case. On the other hand, you may want to run the same binary and switch between the two solvers at runtime for comparison--- but I don't think this is currently possible anyway.

I think it is totally up to users' discretion. One can use CR diffusion for low-energy CR and Yanfei's CR for high-energy CR, for example, although I think such use cases are VERY limited. So I am happy to make them exclusive, but theoretically they are not and that's why I left so.

It is unfortunate that there is no consistency in variable names between the two CR algorithms. E.g. Dpara/Dperp vs. sigma_diff(0:2, ...), initial condition set via pcrdiff->ecr(k,j,i) vs. pcr->u_cr(4,k,j,i). Why does CR diffusion have Lambda, zeta_factor, but CR transport does not seem to have counterparts? Should they be added to CR transport?

First of all, my CR diffusion is designed for more specific problems with different approximation. For my purpose, this is to calculate the ionization rate in star forming clouds.

I can change ecr to u_cr, but I do not think they have to be the same as they have different array dimensions and are in different modules.

Lambda and Zeta are used to calculate the ionization rate which can be used for non-ideal MHD. In Yanfei's CR, one can implement such terms using the source function, but Lambda has to be specified for the full-implicit calculation. I am going to move zeta_factor outside the module.

Some of the options aren't fully explained in https://github.com/PrincetonUniversity/athena/wiki/Cosmic%E2%80%90ray-diffusion-with-Multigrid. You might want to refer to the parent page and/or mention that MGI is the other mgmode possibility. Are the 4x parameters that are commented-out (npresmooth, ... , threshold) supposed to be commented-out? Is everyone required, or are there defaults?
Add a small text subsection and link to CR Diffusion wiki page on https://github.com/PrincetonUniversity/athena/wiki/Cosmic-Ray-Transport

I am going to modify the documentation.

Should we rename pgen/cr_diffusion.cpp to pgen/cr_diffusion_transport.cpp to make it obvious and contrast it with the new pgen/cr_diffusion_mg.cpp?

cr_diffusion_transport sounds confusing to me.

Would be great if @yanfeij could briefly weigh-in, since he authored the other CR module.
(not for this PR, necessarily) https://github.com/PrincetonUniversity/athena/wiki/Cosmic-Ray-Transport only describes cr/vmax and EnrollOpacityFunction(). Need to describe some other parameters and functions in the class, especially those related to CR streaming and source functions. Also would be good to include the counterpart of the "equations solved" copied from Jiang and Oh (2018), assuming that they are different from the CR diffusion equations with the inclusion of the velocity-related terms. This is related to making the names consistent between CR diffusion and transport modules. Some items from cr.hpp header:
In general, it would be good to have more examples of CR transport pgens, and at least one that shows how to use EnrollUserCRSource(), etc.

@yanfeij ?

[felker]

Documentation changes look good, except I noticed that for
https://github.com/PrincetonUniversity/athena/wiki/Self-Gravity-with-Multigrid
you removed omega and both #npresmooth and #npostsmooth remain commented out.

And for https://github.com/PrincetonUniversity/athena/wiki/Cosmic%E2%80%90ray-diffusion-with-Multigrid
threshold remains commented-out

[tomidakn]

I've pushed slightly updated code, and I also updated the wiki.