IllinoisGRMHD is incompatible with setting TmunuBase::stress_energy_at_RHS = "no"

[Gabriele Bozzola]

IllinoisGRMHD schedules the conservative-to-primitive methods in the group AddToTmunu. This group is not added to the scheduler unless TmunuBase::stress_energy_at_RHS = "yes". If this parameter is set to no, IllinoisGRHMD doesn’t do the conservative-to-primitive step, quickly resulting in NaNs.

‌

[Roland Haas]

Assigning to IllinoisGRMHD maintainer (@Zach Etienne ).

[Zach Etienne]

Thanks for this report @Gabriele Bozzola . My postdoc @Leonardo Werneck and I will add an error message to IllinoisGRMHD if both parameters aren’t set correctly. Hopefully this is fixed before the upcoming release.

[Gabriele Bozzola]

In the future, it would be nice to have a permanent fix. I am running some simulations where I need to set stress_energy_at_RHS to no, and what I am currently doing is to change TmunuBase to add always the AddToTmunu group. I think that a clean solution is to change how the primitive solver is scheduled in IllinoisGRMHD avoiding using the AddToTmunu group.

[Gabriele Bozzola]

A opened a PR (https://bitbucket.org/zach_etienne/wvuthorns/pull-requests/8/fix-2497) which should fix this problem.

[Gabriele Bozzola]

Upon further inspection, the PR does not solve the problem. The group SetTmunu is not added to the scheduler when TmunuBase::stress_energy_at_RHS = "no"

I don’t know how to fix this in a clean and easy way. The scheduler looks like this

            GROUP HydroBase_Con2Prim: Convert from conservative to primitive variables
          GROUP SetTmunu: Group for calculating the stress-energy tensor
            TmunuBase::TmunuBase_ZeroTmunu: Initialise the stress-energy tensor to zero
            GROUP AddToTmunu: Add to the stress-energy tensor here
               IllinoisGRMHD::IllinoisGRMHD_conserv_to_prims: Compute primitive variables from conservatives. This is non-trivial, requiring a Newton-Raphson root-finder.
               IllinoisGRMHD::IllinoisGRMHD_outer_boundaries_on_P_rho_b_vx_vy_vz: Apply outflow-only, flat BCs on {P,rho_b,vx,vy,vz}. Outflow only == velocities pointed inward from outer boundary are set to zero.            

IllinoisGRMHD computes Tmunu when it is doing primitive-recovery (in IllinoisGRMHD_convser_to_prims), so, it has to be scheduled after TmunuBase_ZeroTmunu. However, the group SetTmunu does not exist at all when `TmunuBase::stress_energy_at_RHS = "no"

[Roland Haas]

Part of the issue may be that ILGRMHD’s Tmunu code actually does non-Tmunu things namely con2prim as far as I remember. So not scheduling the Tmunu routines is not possible.

[Gabriele Bozzola]

It’s the opposite: it is the con2prim that does non-con2prim things. My workaround to achieve what I want is to change `Tmunubase`

[Roland Haas]

@Zach Etienne is there a chance of this being addressed? It has been wrong for a while? This would fall under maintaining code that one has contributed to the Einstein Toolkit.

[Leonardo Rosa Werneck]

I will fix this for the upcoming release.

[Leonardo Werneck]

After reviewing the issue more carefully, I think it is more prudent to fix it for the November release, as there won’t be enough time to fully test the code in time for this release. One possible solution is to modify IllinoisGRMHD’s schedule.ccl to add the con2prim function to the HydroBase_Con2Prim bin when TmunuBase::stress_energy_at_RHS="no", e.g.,

if(stress_energy_at_RHS)
{
  schedule IllinoisGRMHD_conserv_to_prims in AddToTmunu after compute_B_and_Bstagger_from_A
  {
    LANG: C
  } "Compute primitive variables from conservatives. This is non-trivial, requiring a Newton-Raphson root-finder."

  schedule IllinoisGRMHD_outer_boundaries_on_P_rho_b_vx_vy_vz in AddToTmunu after IllinoisGRMHD_conserv_to_prims
  {
    # We must sync {P,rho_b,vx,vy,vz} here.
    SYNC: grmhd_primitives_allbutBi
    LANG: C
  } "Apply outflow-only, flat BCs on {P,rho_b,vx,vy,vz}. Outflow only == velocities pointed inward from outer boundary are set to zero."
}
else
{
  schedule IllinoisGRMHD_conserv_to_prims in HydroBase_Con2Prim after compute_B_and_Bstagger_from_A
  {
    LANG: C
  } "Compute primitive variables from conservatives. This is non-trivial, requiring a Newton-Raphson root-finder."

  schedule IllinoisGRMHD_outer_boundaries_on_P_rho_b_vx_vy_vz in HydroBase_Con2Prim after IllinoisGRMHD_conserv_to_prims
  {
    # We must sync {P,rho_b,vx,vy,vz} here.
    SYNC: grmhd_primitives_allbutBi
    LANG: C
  } "Apply outflow-only, flat BCs on {P,rho_b,vx,vy,vz}. Outflow only == velocities pointed inward from outer boundary are set to zero."
}

This fixes the scheduling problem, but the code will still keep adding to Tmunu unless one also sets IllinoisGRMHD::update_tmunu = "no" in the parfile.

Could you please test this to see if it provides the functionality that you seek @Gabriele Bozzola ? If I misunderstood something, please let me know.

[Vikram Manikantan]

Hi Leonardo,

In your suggested modification above, do you mean to replace the current lines that have ‘SetTmunu’ or just append this IF statement to the end of the file? See below for existing code in schedule.ccl.

schedule IllinoisGRMHD_conserv_to_prims in SetTmunu after (compute_B_and_Bstagger_from_A, TmunuBase_ZeroTmunu)
{
  LANG: C
} "Compute primitive variables from conservatives. This is non-trivial, requiring a Newton-Raphson root-finder."

schedule IllinoisGRMHD_outer_boundaries_on_P_rho_b_vx_vy_vz in SetTmunu after IllinoisGRMHD_conserv_to_prims
{
# We must sync {P,rho_b,vx,vy,vz} here.
  SYNC: grmhd_primitives_allbutBi
  LANG: C
} "Apply outflow-only, flat BCs on {P,rho_b,vx,vy,vz}. Outflow only == velocities pointed inward from outer boundary are set to zero."

‌

[Vikram Manikantan]

Hello @Leonardo Werneck ,

I implemented the changes in schedule.ccl in IlllinoisGRMHD as suggested, but it is still producing NaNs. See the output log of the simulation below.

Please also see my above comment about replacing the SetTmunu group versus appending to the scheduler. Let me know what is the correct approach. Thanks.

‌

[Leonardo Werneck]

Hi @Vikram Manikantan ! Would it be possible for you to share a parameter file that results in this error?

[Vikram Manikantan]

Parameter file that produces this issue

[Vikram Manikantan]

Hi @Leonardo Werneck ! Yes, I am attaching it above.

I also want to follow up on a previous question. Currently there is a SetTmunu scheduled in schedule.ccl and in your solution there is an AddToTmunu in the if statement. I am wondering if both should be in the file, or just the AddToTmunu?

‌

[Samuel Cupp]

I’m not sure what your schedule.ccl currently looks like, but the master branch of `IllinoisGRMHD` only has those two functions (which I’ll refer to as con2prim) in AddToTmunu, not SetTmunu. Did you also add it to SetTmunu? All Leo did was change the scheduling from something that looks like

schedule con2prim in AddToTmunu

to

if(stress_energy_at_RHS)
{
  schedule con2prim in AddToTmunu
}
else
{
  schedule con2prim in HydroBase_Con2Prim
}

‌

[Samuel Cupp]

If you could also post the scheduling output at the start of the run, that might help me know what is wrong.

[Vikram Manikantan]

Hi Samuel,

I adjusted my schedule.ccl so it looks like the master branch - the SetTmunu was an error on my end. My file contains the if(stress_energy_at_RHS) as described by Leo and yourself.

I am attaching an out file from my most recent test run. Let me know if that is enough and if there is anything else I can provide. Thanks for your help!

[Vikram Manikantan]

[dot]out file from most recent test run

[Samuel Cupp]

If I make the change in IllinoisGRMHD’s schedule.ccl and run the test parfile magnetizedTOV.par, I don’t get any issues with stress_energy_at_RHS on or off. As such, I’m not sure if the issue is with IllinoisGRMHD itself. Particularly, the values that are nan when Con2Prim starts are the B field and the metric. The B field I could envision potentially having an issue if the setup for the scheduling were wrong, but the fact that all the spacetime quantities are nan seems to suggest there could be something wrong with the spacetime evolution or ID.

‌

[Samuel Cupp]

My previous comment was incorrect, as I forgot to also set update_Tmunu. I do get nans when I set that. I haven’t tried to debug this, but my guess is that moving this to a different bin leads to the A->B routine happening in a different place relative to con2prim, which causes this issue.

[Samuel Cupp]

To have more correct scheduling, change the in HydroBase_Con2Prim to in HydroBase_PostStep after HydroBase_Con2Prim.
This is because that con2prim group also runs in PostPostInitial where IllinoisGRMHD manually schedules functions. This should only put it in the places where AddToTmunu appears. I still get nans, but I don’t have time right now to track down where they come from. If you track down where they first appear, please let me know.

[Vikram Manikantan]

The scheduling is already in HydroBase_Con2Prim after IllinoisGRMHD_conserv_to_prims. Do you mean to make it in HydroBase_PostStep after HydroBase_Con2Prim after IllinoisGRMHD_conserv_to_prims or just in HydroBase_PostStep after HydroBase_Con2Prim?

[Vikram Manikantan]

Also, I do not have update_Tmunu in my scheduler anywhere. What do you mean when you say you set this?

[Samuel Cupp]

It’s been a while since I looked at it, but I think it should be in HydroBase_PostStep after HydroBase_Con2Prim. update_Tmunu is a parameter, not something in the scheduler. If you don’t want IllinoisGRMHD to ever set Tmunu, then you need to set that to no in the parfile. It defaults to yes. I didn’t get nans with update_Tmunu=yes, but I did with update_Tmunu=no.

[Roland Haas]

There was some discussion on this in the Einstein Toolkit call just now. The situation seems that one can enable Tmunu storage but disable Tmunu being actually set by setting TmunuBase::stress_energy_storage = yes and TmunuBase::stress_energy_at_RHS = no . TmunuBase relies on running in MoL_PostStep. Since Tmunu is a derived quantity like the primitives it may be most consistent with other codes to have Tmunu set at INITIAL at the same time say HydroBase_Con2PrimInitial is running by scheduling the SetTmunu group in there instead of only in MoL_PostStep. This would rely on checking on ADMSbase_SetADMVars and possible touch HydroBase as well in case setting Tmunu uses primitive variables.

[Vikram Manikantan]

I will set update_Tmunu=yes and see if that fixes it for me - at least we’ll be on the same page then. However, in the long term, this will make the simulation more expensive.

Shouldn’t Con2Prim work regardless of the update_Tmunu setting?

[Samuel Cupp]

Well, I don't really know what you're doing, but update_Tmunu=”no” means that IGM will never set Tmunu. I'm guessing you want it to be set at least at initialization.

The core issue is presumably that Tmunu is never set. All the parfiles I have will expect Tmunu to be set, and it would be a miracle if things didn't go to nan.