To do
Here are things I will be working on over the next few weeks:
- I am currently running two simulations to see how pulsed jets evolve, one with Beta = 1 and another with Beta = 5.
- Search the literature to find the conditions required for pulses to form shocks. From there, I will run a simulation to test it.
- Learn how radial equilibrium is disrupted after passing through a shock for a cylindrically stratified jet. I've already worked out some aspects of this but theres more to be done here.
- Contact Jonathan to ask about warm starts, axisymmetric simulations, and auxiliary variables in astrobear.
Numerical Methods
I have extended my Riemann solver for use in a Godunov reconstructions. I then ran the results using Riemann-problem like setups (see table 6.2 and figures 6.8-6.12 of Riemann Solvers and Numerical Methods for Fluid Dynamics by Toro) and compared the Godunov solution (blue) to the exact Riemann solution (green). Plots are (top to bottom) density, velocity, pressure.
One caveat worth mentioning is that I didn't implement a method to compute timesteps and instead just use a sufficiently small fixed tilmestep. Not the most efficient way to do it but it allows me to verify that I understand the concept behind a Godunov solver, and it still ran fast enough where I couldn't tell
Since then I have implemented various approximate solvers. Examples include Two-Rarefaction Riemann Solver (TRRS), Primitive Variable Riemann Solver (PVRS), and the Roe solver. Pictured below is the same Godunov method using the Roe solver:
Approximate solvers allow us to bypass the iterative nature of the exact Riemann solver, which allows for a considerable speedup when simulating more complex problems. Most approximate solvers however suffer one or more shortcomings. For example, the roe solver tends to do well with discontinuities, but often fails with rarefactions as a result of entropy violation. One option for handling this is to use an adaptive solver, which selects a solver based the type of solution expected.
Whether exact or approximate solvers are used, the Godunov method used so far is accurate only to first order. Improvements can be achieved using various techniques, such as the Weighted Average Flux (WAF) and Total Variation Diminishing (TVD) schemes. Compared to the first order methods, discontinuities and rarefactions are much sharper when using higher order accuracy, at the expense of an increased chance of various numerical artifacts.
Below I have implemented a WAF-TVD scheme using an exact solver. A notable challenge with this scheme in particular (which gave me some trouble) is that for most solvers (including the exact solver) an extra step is required to handle rarefactions.
Lastly I have been reading about splitting schemes, which allow for solutions to equations with source terms. I have not yet implemented these, but these will be particularly important for me to understand since my research focuses on the effects of a source term (technically it's a sink term) in the energy equation.
Initial Runs with unequal flows and realistic cooling
The first run here is a reference case with equal flows and a realistic cooling curve.
For the next six runs, the top jet remains unchanged, with only the bottom jet being changed. First, we lower the density of the bottom jet while keeping the radius fixed.
Next, we have two runs where the radius is varied and the density is varied proportional to the inverse square of the radius; this keeps the mass of any 'gaussian pillbox' fixed (provided the pillbox contains the entire radius of the jet.
Here the radius is again varied, but this time density remains unaltered from the reference case. We see that a change in the radius of the jet does not cause the interaction region to move after collision, but does have a noticeable impact on the spray.
Finally, we change the jet velocities from 70 and 70 km/s to 80 and 60. In the first of these two runs density remains unaltered from the reference case, while in the second run density is increased to keep the mass flux fixed. Note that the top jet still has a higher ram pressure and thus the interaction region still moves.
Next step is to rerun a few of these at a larger scale. Some of the runs are being run in a 128x128x128 box (original is 64) at lower refinement to better study the spray, while another set of runs (not mutually exclusive) are being rerun in 32x32x16 at higher refinement to better resolve the interaction region. I began by rerunning the reference run in both cases, while I work on selecting 2-3 additional runs for each type.
Even with the restricted region, the higher level of refinement results in the interaction region run only about 10-15 frames per day. Currently it is on frame 117, shown below:
Meanwhile rerunning to focus on the spray produced nearly 900 frames within a period of about a week
Movies of 1D shocks that move with the shock
I've generated a new set of 1D movies. They are the same as the previous one except that the spatial coordinates are shifted so that 0.0 moves with the shock. This allows us to watch instabilities evolve as the shock moves. Values of Beta start at -1.0 in the top left, increasing by 0.5 moving to the right and 2.0 by moving down
Runs with larger Jets
From left to right: slice of density, slice of temperature, projection of density along the axis, projection of density along the side,
Beta = -1.0, Alpha increased x10:
(projections will be done when the run completes)
Beta = +1.0, Alpha increased x10:
Beta = +3.0, Alpha increased x10:
Courtesy of Visit, two additional movies of the last run sliced in a different location than usual:
Update 8 February: Added two additional values of Beta. Original post only had Beta=+1.0
Update 16 February: Added additional run at Beta=3.0 with a larger Alpha
Update 2 March: Added runs at Beta= ± 1.0 with larger Alpha
Integer Values of Beta
The following runs were conducted using values of Beta ranging from -1.0 (far left) to +3.0 (far right). For each type of plot, the integer values are in the top row while the half-integer values are in the second row.
Progress as of 19 November 2020
These first plots are taken from one-dimensional runs. These plots trace the position of the shock front as it moves through time. The red curve is taken as the point of maximum temperature, while the blue curve is taken as the point of steepest increase in the logarithm of density. On the right of each of these is the Fourier transform of the temperature plot. In cases where the temperature is constant throughout (which occurs if the shock isn't fully resolved; this only seems to occur for negative values of beta) the peak temperature location is taken to be the location of the previous peak.
It should be noted that there are usually two regions of high density increase; the forward of which usually coincides with temperature peak while the rear usually coincides with the point where temperature returns to the ambient value. The steeper of the two is usually, but not necessarily, the rear. This explains some of the sharp oscillations observed in the higher values of beta as the front gradient becomes as steep as the rear.
Beta | Position | Fourier Transform |
---|---|---|
-1.0 | ![]() | ![]() |
-0.5 | ![]() | ![]() |
0.0 | ![]() | ![]() |
+0.5 | ![]() | ![]() |
+1.0 | ![]() | ![]() |
+1.5 | ![]() | ![]() |
+2.0 | ![]() | ![]() |
+2.5 | ![]() | ![]() |
+3.0 | ![]() | ![]() |
These next sets switch to three dimensions, all using Beta=0 (constant cooling function). From left to right, the (estimated) cooling length is set to be approximately equal to five jet radii, one jet radius, half the jet radius, one fourth, and one eighth.
And another projection along the jet axis
Currently I am conducting runs with varying Beta, using the same value of alpha that produced the cooling length equal to five jet radii and all of the values of beta that were used in the one-dimensional runs. These runs are also conducted over a larger area. So far I have preliminary results for Beta=1 (Left) and Beta=2 (Right)
Varying Cooling Length
All of these runs use Beta=0, which corresponds a constant cooling function. From left to right, the (estimated) cooling length is set to be approximately equal to five jet radii, one jet radius, half the jet radius, one fourth, and one eighth.
And another projection along the jet axis
UPDATE 11/06: added cooling length = 5 Jet Radii
Update 5 October
3D runs are finished, hopefully they will be analysed this week.
For the one dimensional runs, I have now plotted location of the shock front using three different methods: the highest temperature (red), the highest velocity gradient (green), and the highest density gradient (blue). I found however that density gradient tends to be highest behind the shock front where the cold slab builds up, rather than at the initial density jump that occurs at the shock. These plots are from Beta=-1 to Beta=3.5 in half integer increments
Lastly I have built a reading list for the next few weeks from various papers which cite Strickland+Blondin (1995) and appear relevant based on their abstracts.
Long Run in 3D
The following run uses a larger radius than my previous 3D runs, and runs for a significantly longer amount of time
For those wishing to create large gifs like this (particularly if you want to do so on your own computer), I have a few tips to share:
- The 'convert' command is extremely resource intensive if you want to combine thousands of images. To save space and time, consider making smaller movies first and then combining (ignore the underscores, they are their to get the wiki to accept the following code)
_convert img/img00*.gif movie_00.gif _convert img/img01*.gif movie_01.gif ... _convert movie_*.gif movie.gif _rm movie_*.gif
- I have found other gif manipulators (eg gifsicle) perform better than convert
A Series of runs in one dimension
I have done a sequence of runs varying beta (the power of temperature dependence) for the one dimensional problem. These runs are set with the following parameters in order to achieve a cooling length of approximately 1 cm:
1E17 particles per cubic cm
Flow temperature is 720K
Mach 15
This gives us T_0 (shock temperature) of 5.06250E4 K. To achieve the desired cooling length at this temperature we use alpha = 8.27393E-23 ergs/s per cubic cm.
Beta varies between the runs.
Beta = -1.0
Beta = -0.5
Beta = 0.0
Beta = 0.5
Beta = 1.0
Beta = 1.5
Beta = 2.0
Beta = 2.5
Beta = 3.0
Update 9/5
Since my last post, I have done the following
I have done a sequence of runs varying beta (the power of temperature dependence) for the one dimensional problem. The values of beta (in order) are -1.0, -0.5, 0.0, 0.5, 1.0, 2.0. I also ran -2.0 but decided not to analyse it yet since -1.0 and -0.5 were similar.
I have done another reading of "Numerical Analysis of the Dynamic Stability of Radiative Shocks" (Strickland and Blondin 1995) I also did my first reading of "The structure of knots in variable stellar jets - I. Symmetric knots" (Falle and Raga 1993)
Lastly, I did a few more runs in the 3D Jet collision problem. First, I did a sequence of runs in which the Jets are (anti)parallel but not aligned, with the offset decreasing with each successive run. As I expected, a portion of each jet continues to propagate forward relatively unaltered by the collision. If the separation is not to great, another portion moves diagonally having been shaped by the collision; The angle at which this portion moves however did not depend as significantly as I had expected on the separation between the beams
A final run returns the jets to axial alignment, but makes one jet significantly hotter than the other
Update 8/13/20
This run is a repeat of one of my previous runs (two identical heavy jets), except with more controlled refinement:
most of the run only goes to level 2, but there is a sphere of radius 0.5 in the center that goes to level 4 and a larger sphere of radius 1.0 that goes to level 3. This allows us to see a higher resolution in the region of collision without wasting computational resources where they are not needed.
This is a run for a 1D Jet with cooling
Labels on the Axes will be added soon. For now, red is log of temperature, blue is log of density, and green is velocity divided by 100
I also read "Numerical Analysis of the Dynamic Stability of Radiative Shocks" (Strickland and Blondin 1995)
Update 8/16:
Labels have been added for axes, but still need to add colour key
My first blog post
To familiarise myself with AstroBear and various parameters, I have been running a problem which involves a collision between two jets
Trial 0: ran the problem as provided
Trial 1: Changed density of top jet from 6e16 to 3e16
Trial 2: Same as Trial 0, but increased density of ambient medium to 1e17 Took several attempts to get Visit to cooperate
Trial 3: Combine changes of both Trial 1 and Trial 2
Trial 4: Repeat Trial 3 with 160 frames, final time changed from 0.2 to 0.4
Trial 5: Same as Trial 4, but use ambient density 4e16 (i.e. between the densities of the jets). I accidentally overwrote the chombo files for this run so I will have to rerun it if I want anything other than this gif
Trial 6: Reran trial 2 at 4 levels of refinement. Since I'm running in debug this limited me to 22 frames (which I have since lost).