Add post-processing callback for initial condition

[EdwinChan]

This feature request references PrincetonUniversity/athena-public-version@20145fe.

According to the wiki, MeshBlock::ProblemGenerator() describes the initial condition mesh block by mesh block. However, mesh blocks (or their owning processes) can't communicate with one another at this stage. This is inconvenient if the initial condition must satisfy some global constraint. A fix is to introduce a post-processing stage for the initial condition.

For example, we may want the total mass in the simulation domain to match some target value. If the total mass isn't analytic, we must calculate the total mass from the output, and start a new run with different parameters. This approach is very error-prone.

The proposed diff is as follows; the naming of the callback is up to the code maintainers. The callback is an "epilogue" to the initial condition, so it is in some sense the opposite of Mesh::InitUserMeshData(), which is the prelude:

diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp
index 70f2ef1..a018d5b 100644
--- a/src/mesh/mesh.cpp
+++ b/src/mesh/mesh.cpp
@@ -1380,6 +1381,10 @@ void Mesh::Initialize(int res_flag, ParameterInput *pin) {
       for (int i=0; i<nblocal; ++i) {
         MeshBlock *pmb = my_blocks(i);
         pmb->ProblemGenerator(pin);
+      }
+      Mesh::ProblemGeneratorEpilogue();
+      for (int i=0; i<nblocal; ++i) {
+        MeshBlock *pmb = my_blocks(i);
         pmb->pbval->CheckUserBoundaries();
       }
     }

[felker]

It is an interesting idea. I would like to see a full implementation of this before offering my opinion one way or the other. And also some more use cases that would really motivate this new feature (in your example there are likely multiple less intrusive alternatives to getting the same initial condition).

Also, what does 20145fe have to do with this?

[EdwinChan]

@felker Thanks for the reply. I'm offering this up as an idea; I'd be interested to see what you guys think would be the best way to implement it.

For what it's worth, my use case is I want to scale the magnetic field so that the volume-integrated plasma beta matches some specified value. Because I have multiple magnetic-field topologies, changing topology would mean manually changing the magnetic field strength as well, which is a disaster waiting to happen.

The commit reference is to give a reference point to interpret the text and the diff. For example, issues at matplotlib have started requiring detailed version information, so I'm trying to do something similar.

[c-white]

I have the same issue with tuning plasma beta in every simulation I run, so I definitely understand the situation.

Fortunately, this particular problem converges very quickly with the workflow of dump initial conditions -> measure beta -> adjust field strength parameter.

One thing that occurs to me, though, is there might well be cases where the tuning amounts to some nonlinear root find requiring multiple iterations, with the initial conditions depending in a complicated way on the outcome of some analysis. If the Epilogue function determines convergence has not been attained, will it have to signal the code to retry ProblemGenerator with new parameters? This could rapidly get complicated.

I'll also add the one can add (and people have done this hack) MPI calls to ProblemGenerator, with care. If you have 1 block / rank, it will just work. With more than that, you have to make sure only 1 block per rank is directly involved in communication.

[EdwinChan]

Manual iteration is certainly doable, but it's always a good idea to take the human factor out of mechanical operations.

Regarding initial conditions that need iteration, one could make MeshBlock::ProblemGenerator() a no-op and have all the iterations done in the epilogue. The user could define, for example, setup_trial_ic(MeshBlock *mb, T params) and call it in the iteration. In the worst case, the epilogue could fail fast and exit, inviting user intervention.

[stale]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.