Castro
Castro copied to clipboard
AMReX-Astro
Provide a 3d rotating star test problem
A three dimensional rotating star test would provide a relevant astrophysical test problem for Castro as well as allow a direct comparison with Octo-Tiger.
it would be easiest with some simple initial conditions. And in particular, we need to understand how to do the self-consistent initialization to create the initial model in equilibrium if it is not spherical (because of the rotation), so we should use the same method there.
The SCF implementation we do have already can account for a rotating reference frame and can generate non-spherical stars, although I don't have much practice with it yet. I don't have much experience generating the equilibrium setup in the rotating frame and then converting it to the inertial frame, though it seems straightforward to try (I tried this once a while back in the WD merger problem but didn't get very far).
@dmarce1 I figured I'd try something out to get the conversation started. If you take a look at the scf_tests/single_star problem setup, you'll find a setup that invokes the Hachisu self-consistent field method to set up a single star. Although we haven't published on our implementation yet, you can find a brief writeup in the documentation; it is methodologically similar to the way the LSU group has done it based on the papers in the 1990s and 2000s that I am familiar with.
Right now the problem setup uses the electron degenerate Helmholtz EOS since that is what we mostly used for the nuclear astro science problems, but I spent some time this week validating it to make sure it would also work for a polytropic EOS. Our Microphysics repository has the EOS for a polytrope and this can be accessed by setting EOS_DIR=polytrope in the makefile for the problem. The user can select the constant K and the index n, but we also have shortcuts for a fully degenerate electron gas that is either relativistic (gamma = 4/3) or non-relativistic (gamma = 5/3). So for example you can get a relativistic electron-degenerate WD with eos.polytrope_type = 2 in the inputs file.
Using a setup with this EOS, I came up with a semi-arbitrary non-rotating spherical initial model using the SCF method, with radius 1.5445e9 cm and central density 1.0e7 g/cc. The mass of this model happens to be about a Chandrasekhar mass. Keeping the equatorial radius constant and arbitrarily selecting a polar radius of 1.4e9 cm, I then varied the central density until I got a model with a reasonably low virial error, which turned out to be at a central density of 1.245e7 g/cc. (The model has a comparable mass to the nonrotating case, still about 1.4 solar masses.) The rotation period of this model for these parameters is about 60 seconds. I've uploaded the inputs files for the rotating and nonrotating cases to a branch on my fork in case you want to look at them.
So my thinking is, we can generate the SCF initial rotating model and then evolve the star for some number of fixed rotation periods, and compare the density profile of the star to its initial value; the closer it is, the better the equilibrium is maintained. For comparison, we can evolve the analogous non-rotating model for the same amount of simulation time to get a baseline estimate of how well equilibrium is maintained at that resolution, and comparing the two will give us a sense of how well the code deals with rotation.
At this point I mainly want to know if you think I'm on the right track, ignoring the specific details. I know I have some homework to do if this is the path we want to go down. First, I've never really validated the rotating models obtained with my implementation of the SCF method, so I don't know whether they are correct and/or high quality. Second, I don't have any infrastructure to automate the process of obtaining a low-virial-error initial model with SCF (I've just been tinkering by hand via trial and error the few times I've needed to do this). Third, I assume you'll want to study binaries at some point, and right now the implementation only does a single star (though I generally know how to do the binary star case in SCF and have implemented it before). This would be a good opportunity for me to harden and improve the implementation, I just don't want to go too far down that path if you have something else in mind.
@maxpkatz Yes, it looks like you're on the right track. Octo-Tiger has a separate tool for producing the 2d rotating star. It uses a 2d SCF with a 2d Poisson solver. It's not highly optimized (the boundary values for gravity are computed by direct summation) but it gets the job done. We generally use a 100x100 grid with the larger axis of the star occupying 90% of the grid length and the smaller axis 60%. We get a virial error around 5x10^-5. The tool outputs a binary file then Octo-Tiger reads it an does a bicubic interpolation to map it to the grid. If you think it'd be useful I can point you to it.
It doesn't look like our SCF initial model is doing a good job yet establishing a polytropic star that's actually in equilibrium. Here is the case I mentioned above without any rotation (just attempting to use SCF to generate a star in stationary hydrostatic equilibrium). For a 128^3 grid I can carefully pick a radius that yields a virial error of ~1e-5, but the star collapses over the timescale of a few tens of seconds. This does not seem to improve with resolution. I am not sure yet what the problem is; the initial profile looks right to me by eye for an n = 3 polytrope.
OK, for now I'm going to stay away from gamma = 4/3 polytropes and their fun physical properties (I knew they were unstable to gravitational collapse but I didn't know they were unstable to auto-exploding!). Here is a gamma = 5/3 polytrope (non-relativistic degenerate electrons, eos.polytrope_type = 1), I will use this going forward.
Is the n=5/3 working? I would have expected the n=3 to do weird things if evolved with an ideal gas of gamma=5/3
The first plot was non-rotating n = 3, gamma = 4/3 polytrope. I think this is a thermodynamically consistent EOS but doesn't seem to result in a stable star. The second plot was non-rotating n = 1.5, gamma = 5/3 polytrope and this seems to work and is stable. Next step is to build a rotating version of that star and see how it does.
Yes that should be stable. I'm confused as to why it wasn't. Perhaps since it is exactly on the line between stability and instability, numerically for some reason the 4/3 ends up just barely unstable while the 5/3 just barely stable.
I was wondering if any progress had been made lately with the rotating star and/or if there is anything we can do to help.
@dmarce1 Sorry for the delay, I was on vacation for a little while. I hit a bit of a roadblock because I still don't feel that our SCF implementation (even for a single star) is robust. I have not lost interest in the problem but I don't have a lot of time to think creatively about it. If you were to offer help, and that would be gracious but certainly not expected, it would be to take a look at what we're doing there and offer thoughts about how we could improve the implementation.