VL+CT support gravity source terms #186
Conversation
For testing purposes, I have introduced support for an inclined linear Jeans Wave (users can choose to make it stable or unstable based on the value of the Gravitational constant)
Also did some minor cleanup.
While I was here, I decided to refactor the circularly polarized alfven wave so that it used lambda functions to define its vector initializers. Doing this makes the code a lot easier to read - all of the code specific to initializing a Circularly-Polarized Alfven Wave is now located in a single spot (previously it was split into 2 different places). I would like to clean this all up a little more, but it's not a great use of my time
- describe the changes to the inclined_wave initializer - describe the actual Jeans wave test
|
@peeples this is ready for review |
There was a problem hiding this comment.
Generally looks good to me.
Would it be worth making an EnzoSource base class, and have EnzoSourceGravity (and potentially other types of sources) be a class derived from EnzoSource? I understand this might be too much work for now and you would prefer to leave that for a a future PR. Could feedback and accretion from / onto star / sink particles be handled in a similar way?
|
@stefanarridge this looks ready to merge, modulo that one conflict? |
Yes I'd say so. |
- fixed an oversight in EnzoMethodMHDVlct::compute where we never initialized the map of acceleration fields - fixed an oversight in EnzoSourceGravity where we were accidentally trying to read acceleration fields from the wrong map. - fixed an oversight in the default constructor for array_StringIndRdOnlyMap that was revealed during debugging. Previously, when we would default construct an array, we would not properly initialize the size and capacity members to zero. This test problem appears to break in a similar way to the PPM test problem. For that reason I have removed the Jeans Wave from the automated testing infrastructure
|
After merging with main, I discovered a massive oversight. It turns out, we weren't actually adding the gravity source terms in the VL+CT solver (I had forgotten to initialize the array map that holds the acceleration fields)... This revealed a few obvious minor bugs in my implementation that I have now fixed. Previously, the test problem must have been evolving a pure hydro wave (which explains the mysterious second order error convergence). After fixing this, an x-axis aligned Jeans wave with the VL+CT solver now looks something like this:
This looks a lot like the result for the PPM solver: All of this seems to suggest that there is probably an issue with my initial conditions (or possibly a problem with the Gravity solver, but that seems a lot less likely...). You can use these input files to reproduce these results: grav_input.tar.gz.
Yep, that's exactly what I had in mind, but I was leaving that for a future PR. (As an aside, there is also a source term for the dual energy formalism, but that has a different signature that it probably wouldn't be subclassed from the same base class). We should talk more about this next week, because it may make more sense to keep certain types of source terms in separate methods. |
|
Uh, actually I just figured it all out (based on a comment @jwise77 made a while back during his review of PR #89). Everything works out extremely nicely when I set Here's the aligned wave for VL+CT: Here's the aligned wave for PPM: This also cleared up a weird inconsistency I was getting about choosing the correct stopping time |
|
As I mentioned, I have done some isothermal collapse runs with the PPM and VL+CT methods and gravity. Here is a link to a Google drive folder with a couple of gifs. They show the evolution of the density field and the position of the sink particle. https://drive.google.com/drive/folders/105ha9f_DmbcWKA6zZxcl09K3UpuKEzdy?usp=sharing I am also attaching a couple of tarballs with more results: Unfortunately there are currently some missing files which I am unable to access at the moment due to technical difficulties. Both contain a "shu.in" which is the parameter file used to run the problem. In the "ppm" directory, there is a "radial_profiles.png", and a "mass.png". These are attempts to reproduce Figures 10 and 9 from this paper: The "vlct" directory contains "shu.out" is is the Enzo-E output file. It shows that the simulation crashed at some point. Unfortunately I don't have access to the python scripts used to make the figures at the moment. Edit: I have managed to get some of missing files now, namely, the python scripts and the "mass" and "radial profile" plots for the vlct run, which you can extract from this tarball: I will try and do these runs again with the "alpha" parameter set to zero. To run with the same version of Enzo-E, checkout commit 23b44da from this repository: https://github.com/stefanarridge/enzo-e. This branch is the subject of PR #200. Let me know if you have any questions. |
Co-authored-by: Greg Bryan <gbryan@astro.columbia.edu>
… test run the problem with the PPM solver.
|
I updated the answer-test and added tests for the PPM solver. Here is the L1 Error-Norm scaling. Both hydro solvers achieve first-order convergence.
If all tests pass, let's merge
Let me know if you reproduced this issue or if you don't have time to do that. If this remains an issue, I think we should address this outside of this PR. |
If you're referring to the issue of how the sphere collapses into a cross shape, I suspect this may have something to do with the initial conditions, and I am going to test this soon. In any case, this does not need to be addressed in this PR. |
Previously the drift time for the gas density was set to the value of `alpha` (taken from `Method:pm_deposit:alpha`). As @stefanarridge pointed out in PR enzo-project#89 this didn't make a lot of sense. I also ran into issues with introducing a Stable Jeans Wave regression test in PR enzo-project#186 that required me to set `Method:pm_deposit:alpha` to 0. For comparison, the drift time for the gravitating mass particles is set to `alpha * dt / cosmo_a`, in which: - `alpha` again comes from `Method:pm_deposit:alpha` - `dt` is the time-step during the current cycle - `cosmo_a` is the scale factor computed at `time + alpha * dt` (where `time` is the time at the start of the current cycle) To investigate this issue I checked the original code from enzo-dev in `Grid_DepositBaryons.C`. I concluded that enzo-dev sets the gas density's drift-timestep should get set in exactly the same way as the Particle's drift timestep. However, when using the PPM and Zeus Hydro-solvers, the drift timestep is set to 0 (this is consistent with a comment @johnwise made during his review of PR enzo-project#189). Since our 2 primary hydro-solvers in Enzo-E (Ppm and VL+CT) currently handle these source terms with the same temporal order as these solvers, we now force the gas-deposition drift timestep to be 0 in `EnzoMethodPmDeposit`. We will need to revisit this exact behavior when we eventually modify the VL+CT solver to handle self-gravity in a higher temporal-order manner.
…posit Previously the drift time for the gas density was set to the value of `alpha` (taken from `Method:pm_deposit:alpha`). As @stefanarridge pointed out in PR enzo-project#89 this didn't make a lot of sense. I also ran into issues with introducing a Stable Jeans Wave regression test in PR enzo-project#186 that required me to set `Method:pm_deposit:alpha` to 0. For comparison, the drift time for the gravitating mass particles is set to `alpha * dt / cosmo_a`, in which: - `alpha` again comes from `Method:pm_deposit:alpha` - `dt` is the time-step during the current cycle - `cosmo_a` is the scale factor computed at `time + alpha * dt` (where `time` is the time at the start of the current cycle) To investigate this issue I checked the original code from enzo-dev in `Grid_DepositBaryons.C`. I concluded that enzo-dev sets the gas density's drift-timestep should get set in exactly the same way as the Particle's drift timestep. However, when using the PPM and Zeus Hydro-solvers, the drift timestep is set to 0 (this is consistent with a comment @johnwise made during his review of PR enzo-project#189). Since our 2 primary hydro-solvers in Enzo-E (Ppm and VL+CT) currently handle these source terms with the same temporal order as these solvers, we now force the gas-deposition drift timestep to be 0 in `EnzoMethodPmDeposit`. We will need to revisit this exact behavior when we eventually modify the VL+CT solver to handle self-gravity in a higher temporal-order manner.
This PR introduces gravity source terms to VL+CT solver.
Overview
The implementation is very simplistic.
The implementation should have the same temporal order as the PPM solver. I adapted the source term formulas from Enzo’s Runge Kutta solver.
Tests
Implemented Test
At this time, I have only implemented a single test. Essentially we run an inclined stable Jeans wave over one period and compute the L1 error norm of the final snapshot compared to the initial snapshot (this is the same procedure that is used by our inclined HD and MHD tests). While introducing this test I made some changes to the infrastructure in
EnzoInitialInclinedWave(there's still a lot more that could be improved), but I took care to make sure I wasn't breaking other tests.The following shows the evolution of a high resolution simulation on a uniform grid with a resolution of (256,128,128):
The above figure was run with a much earlier version of this branch, but from everything I've seen the results shouldn't have changed much. The above figure includes a data point for every cell in the simulation plotted at its position along the (incline) axis that the wave is propagating along.
I have included a plot of the L1 error norm as function of the number of cells. This used a more recent version of the code (from just before I rebased the branch onto the cmake branch (so compiler settings were slightly different). Following the convention that Athena/Athena++ use for their inclined wave convergence study, a domain of
(256,128,128)corresponds toN=128:The dashed line shows 2nd order convergence. I'm not sure why we are getting higher order convergence (presumably this is just related to this particular test problem). I tried changing
Method:gravity:orderfrom 4 (the default) to 2 and saw no difference (beyond round-off error) in this plot.As I've mentioned on the slack, the PPM solver encounters issues with this test problem. I'll open an issue to document this
Other Tests
We could definitely use more tests (especially ones with analytic solutions).
An obvious choice is the unstable Jeans wave (since$\lambda/\lambda_J = 1.5$ .
EnzoInitialInclinedWavecan already generate the initial conditions) The following shows the early time evolution whenTo finish implementing this test we would need to compare against the analytic solution (presumably we could reuse/adapt the machinery from the Shock Tube tests). The analytic solution is the same as the one used in Athena.
Alternative implementation
As an aside: other choices could be made for integrating the solver with the gravity solver.
As a reminder: a van Leer Predictor-Corrector scheme (the the VL in VL+CT) employs the following procedure:
u_begin, to compute fluxes,f_begin.u_half(fromu_beginandf_begin)f_halfu_beginandf_halfAlternative 1
A simple way to try to improve the solver’s integration with gravity is to reuse the same acceleration source terms in both the partial and full time steps. It’s super easy to adapt the code to do this, but I think we should have some more test cases before we do this (it’s not completely obvious to me how well this should work).
This is probably what we want to be doing in the short-term.
Alternative 2
The optimal approach treatment with higher temporal order (what athena++ does) would:
u_half["density"](and the expected particle positions at the half timestep)As in athena++, it would be fairly trivial to use the momentum source terms described by Mullen+ 2021 under this approach - this ensures momentum conservative.
I’ve been thinking about this sort of solution for a while (it has applications for other source terms like radiative cooling) and have some ideas for how to approach it. Unfortunately, it opens a can of worms:
"mhd_vlct"and"ppm"to remain drop-in replacements for each other.So, I’ve been holding off on this for the short term