pdmullen
commented
4 years ago Greetings, @umgraftw!
I hope you are staying safe/healthy in these scary times... I have been thinking about energy conservation with self-gravity recently too, so maybe I can help tackle this question!
(1) A small update: in the Athena++ public repository, when history outputs are enabled, the gravitational potential energy is computed when self-gravity is enabled here: https://github.com/PrincetonUniversity/athena-public-version/blob/d535095157c9b6dd67829dab4ce605d21da6efbb/src/outputs/history.cpp#L101-L105 However, this is using the incorrect gravitational potential. The potential used in the above is the potential associated with the density distribution before executing the final stage of TimeIntegratorTaskList. The potential we truly desire here is the potential associated with the density distribution after executing the final stage. Therefore, in the development repository for Athena++ we have moved the potential solve to after executing ptlist->DoTaskListOneStage(). You can apply this hotfix by changing line 449-455 of main.cpp to
for (int stage=1; stage<=ptlist->nstages; ++stage) {
ptlist->DoTaskListOneStage(pmesh, stage);
if (SELF_GRAVITY_ENABLED == 1) // fft (flag 0 for discrete kernel, 1 for continuous)
pmesh->pfgrd->Solve(stage, 0);
else if (SELF_GRAVITY_ENABLED == 2) // multigrid
pmesh->pmgrd->Solve(stage);
}This fix will likely be implemented in the next public release. If you have been analyzing energy conservation with history outputs, this fix may produce some change (although I don't think this is the major issue here).
(2) I didn't check your pgen and athinput combination very closely (apologies), but one thing I would urge you to recheck and be absolutely confident of is that the mean density that you set with SetMeanDensity() is the mean density of the entire mesh. The mean density appears in the divergence of the gravitational stress tensor when applying Jeans' swindle. If you don't set this correctly, you will see conservation issues (see Issue #36).
(3) Even after you fix (1) and potentially address (2), total energy will not be conserved to machine precision in the public version of Athena++ for self-gravitating flows. To clear a potential point of confusion, this has nothing to do with the FFT Poisson solver, rather, the issue arises from how momenta and total energy updates from self-gravity are handled. Jiang et al. (2013) (https://ui.adsabs.harvard.edu/abs/2013NewA...19...48J/abstract) described a numerical scheme that can guarantee momentum and total energy conservation (to round-off error) by taking the divergence of momentum and energy "gravity fluxes". Athena++ updates the momentum (due to gravity) using such fluxes, but has not implemented the energy analogue (the energy gravity flux is a bit more involved, see Jiang et al. 2013 for details). Hanawa (2019) (https://ui.adsabs.harvard.edu/abs/2019JPhCS1225a2015H/abstract) describes a numerical scheme where one can actually guarantee momentum and total energy conservation if you apply gravity source terms in a special way. Again, this is not the scheme employed by the public version of Athena++.
(4) It is difficult to address when departures from total energy conservation "should be concerning". This is highly problem specific. Again, I would point you to Jiang et al. 2013 to see some of the consequences of violating total energy conservation.
(5) I have recently forked the private repo of Athena++ and developed a branch that applies the Hanawa (2019) numerical scheme for self-gravitating flows. The results are promising. For the linear wave Jean's instability problem (with lambda = 2 lambda_J), the Hanawa (2019) scheme conserves total energy to ~round-off error. This particular problem setup follows the collapse of an unstable wave into a sheet, and subsequently into a cylinder. The plots below show the different components of the energy during the collapse and the total energy conservation, respectively:
The same plot with the public version of Athena++ has total energy varying dramatically. These development efforts are still in early stages, but I am happy to chat about them if you are interested.
(6) Finally, this may be a bit off-the-wall (and may seem counterintuitive), but maybe you could try turning off gravity momentum fluxes in your Athena++ clone (comment out AddGravityFlux() in calculate_fluxes.cpp), and implement a momentum source term for self-gravity (i.e., in src/hydro/srcterms/self_gravity.cpp). For example, for X1-momentum add
cons(IM1,k,j,i) -= dtodx1*prim(IDN,k,j,i)*(phir-phil);
How does this change affect your energy conservation? (If you perform this test, I recommend that you use the rk2 integrator). Hanawa (2019) find that the gravity flux treatment can produce large errors in low density regions subject to strong gravity. This effect does not plague the source term treatment. Maybe this somehow connects to your energy conservation troubles (i.e., maybe these errors lead to hitting a density/pressure floor frequently?).
Yikes! Sorry for the water hose....hopefully something here helps.
Patrick
umgraftw
Greetings,
According to user @pdmullen in #15 one should not expect energy to be conserved to machine precision using the in-built FFT gravity solver, and pdmullen further provided a graph showing a comparison of E{tot,0} to E{tot} versus time.
I am noticing that my simulations are more extreme with their change in total energy versus time, which is concerning for me.
I was hoping that someone could help me to understand why energy isn't conserved in the FFT self gravity solver and under what circumstances that should be concerning to me (ie can I trust the dynamics of the system?)
Here is a simple test we did:
And here is the input file
We are hoping to study the dynamics of a cylindrical filament in a turbulent rarefied medium. Any help understanding the change in total energy with the FFT self gravity solver would be greatly appreciated.
Thanks!