EBM Project Update
Project Description
Figure out how to run Jacob's Fortran Energy Balance Model through Mathematica. Put in inputs within Mathematica notebook, then receive/extract all data through the same notebook.
Idea 1: Use MathLink (status: not working)
- MathLink
- allows external programs both to call Mathematica, and to be called by Mathematica. By using MathLink, you can, for example, treat Mathematica essentially like a subroutine embedded inside an external program. You can also use MathLink to let Mathematica call individual functions inside an external program. (MathLink (as well as LibraryLink) are designed to be used in conjunction with programs written in C/C++. )
Where I'm at here is I have outlined the 4 steps/goals that if accomplished would allow me to successfully use MathLink. I have successfully accomplished the first two steps, yet the other two remain unworking…
- Successfully test MathLink using the example file they have built into Mathematica by compiling the example C++ program and template file into an executable file which is then run through Mathematica.
- SUCCESS
- Figure out how to call a Fortran 90 subroutine from a C++ program.
- SUCCESS
- Successfully test MathLink's ability to call an external FORTRAN subroutine from Mathematica by moving the meat of the example C++ program to a FORTRAN test subroutine, and then compile the example C++ program, the test FORTRAN subroutine, and the example template file into an executable file which is then run through Mathematica.
- UNSUCCESSFUL (I keep on running into various bugs)
- Do the same as step 3, yet instead of the example C++ program, create a wrapper which will call the EBM driver.f90 instead of the test FORTRAN subroutine.
- UNSUCCESSFUL (can only work if I figure out step 3)
Idea 2: Use LibraryLink (status: semi-working)
- LibrayLink
- provides a powerful way to connect external code to the Wolfram Language, enabling high-speed and memory-efficient execution. It does this by allowing dynamic libraries to be directly loaded into the Wolfram Language kernel so that functions in the libraries can be immediately called from the Wolfram Language.
I think I finally have LibraryLink working at least well enough for us to use it in conjunction with Jacob's EBM model. When I run it through Mathematica, I can now…
- Successfully extract:
- globally averaged temperature (double)
- minimum global temperatures by lat band (array size 18)
- maximum global temperatures by lat band (array size 18)
- Unsuccessfully extract:
- average global temperatures by lat band (array size 18)
I think this may not be an issue for our purposes, as it seems like we only need the globally averaged temperature. In addition, if we want an approximated average temp array, I can calculate that just by averaging the min and max temp arrays, as I did in the Mathematica notebook.
Gif Demonstrating that LibraryLink is Working by Showing how the Temperature Arrays change with Obliquity*…
*(the avg temp array used here was calculated by averaging the min and max temp arrays, it was NOT calculated directly from the EBM program)
Convection project update
Project description
Phase 1: Putting OPAL Equation of state (EOS) into AstroBEAR
Phase 2: Run simulations will new EOS
OPAL EOS
History: Developed mainly by Forrest J. Rogers in 1990s, revised in 1996 and 2002
Includes
- non-relativistic Fermi-Dirac electrons
- classical ions
- all stages of ionization and excitation
- molecular hydrogen
- degenerate Columb correction
- quantum electron diffraction
- electron exchange
- pressure ionization
- terms arising from the so-called ladder diagrams of full quantum theory
Excludes
- pseudopotential method for going to higher order in electron-electro and electron-ion interaction (as far as in Rogers 1996)
Accurate to the order of
- Quantum diagrammatic procedure are used to calculate terms to order
- In the case of hydrogen, it agrees with -order correction
Online OPAL EOS table
https://opalopacity.llnl.gov/opal.html
- 4 versions, the lastest is created in 2005 and updated in 2006.
- To determine which table we should use, we need Mentality and Hydrogen or Helium concentration or .
- Data ranges from to ; to
- Some interpolation code written in Fortran is also provided.
Ideal gas EOS
Assumes ideal, adiabatic and monoatomic gas
- adiabatic:
- Energy:
Comparision between Ideal gas EOS and OPAL EOS =
Red: OPAL
Green: Ideal
Pressure - density - temperature
Pressure - density - temperature - energy
X axis: density
Y axis: energy
By the way…
In run 143, frame 46 (reduSol)
min | max | |
rho | 1.483e-09 | 0.003275 |
T | 2721.12 | 8.695e06 |
P | 16007.5 | 1.0592e12 |
Next Step
- Understanding the difference between Ideal EOS and OPAL EOS
| Pressure comparison is much harder because pressure is a dependent variable in the table.
- Understand and use the Interpolation code
- Try to put them in AstroBEAR
Update 5/24
MHD proposal figs
Also, a test of the streamlines. I like the LICs better, since the streamlines have a habit of not ending.
Radiation Pressure Paper
Found issue with postprocessing. Reprocessing synthetic observations. Here's a first look at the difference. The difference is significant, but doesn't change our conclusions.
Charge Exchange
Current frame
We also need to decide on the high-med-low for the stellar wind. Here are the current parameters (taken from John's simulation, except the mass loss rate, which is the solar mass loss rate - which corresponds to a 4x higher wind density than their strongest wind):
Parameter | Value |
cm | |
K | |
cm/s | |
g/s | |
AMR line transfer
Timing run is next up (for comparison with non-AMR and new optimizations), along with work on optimizations.
Also working on a way to output propagation of rays across grid for a step (just for a nice visualization):
Need: ray position, synchronized time
Ray position is easy.
Synchronized time pseudocode:
call MPI_barrier ! To synchronize time - unless MPI_wtime_is_global, in which case all processors are synchronized already and this becomes much easier initTime = MPI_wtime call PropagateRay subroutine PropagateRay time = MPI_wtime - initTime MPI_file_write(time,position) ! Or might be easier to save an array of time, position and write it out after all rays are processed, in processor order (so we don't have to calculate the buffer)
COMMON ENVELOPE SIMULATIONS
New Work
- 2D movies for force paper for fiducial run 143 (q=0.5), using de-resolved data.
- 1D plots for force paper for runs 143 (q=0.5), 149 (q=0.25) and run 151 (q=0.125), using original data (not de-resolved data).
- Redoing force post-processing using reservation on bluehive (pp seems to be more stable with a reservation, but results should hardly change, if at all).
- Deleted most of AGB run from bluehive, keeping only 1/5 of frames, to save space (the other frames are still available for now on stampede scratch, but I may move them to backup drive).
- Analysis: begun comparison with Macleod+17b
Force Paper
New notes with 1D plots
Movies
all movies are in the corotating (with the particle orbit) reference frame of particle 2, and particle 2 is at center of the frame, while particle 1 is situated on the
Particularly interesting movies are highlighted in bold font
Force per unit volume between particle 2 and gas along phi-direction (magnitude), with contours also showing the same quantity plotted in color with contour values equal to values labeled on the color bar, and with vectors on particle 2 showing velocity of particle 2 relative to particle 1, and on particles 1 and 2 showing net force on each particle exerted by the gas , slice through orbital plane (view size ):
Force per unit volume along phi direction
Force per unit volume between particle 2 and gas along direction of particle 2 velocity relative to particle 1 (magnitude), with contours also showing the same quantity plotted in color with contour values equal to values labeled on the color bar, and with vectors on particle 2 showing velocity of particle 2 relative to particle 1, and on particles 1 and 2 showing net force on each particle exerted by the gas, slice through orbital plane (view size
Force per unit volume along velocity wrt particle 1
Force per unit volume between particle 2 and gas along phi-direction (magnitude), with contours also showing the same quantity plotted in color with contour values equal to values labeled on the color bar, and with vectors on particle 2 showing velocity of particle 2 relative to particle 1, and on particles 1 and 2 showing net force on each particle exerted by the gas , slice through orbital plane (view size
Force per unit volume along phi direction
Force per unit volume between particle 2 and gas along direction of particle 2 velocity relative to particle 1 (magnitude), with vectors on particle 2 showing velocity of particle 2 relative to particle 1, and on particles 1 and 2 showing net force on each particle exerted by the gas, slice through orbital plane (view size
Force per unit volume along velocity wrt particle 1
Density, slice through orbital plane, with velocity vectors (view size
Density
Density normalized to
Normalized density, slice through orbital plane
Density normalized to
Normalized density, perpendicular to orbital plane
Mach number, in frame of particle 2 corotating with particle orbit, slice through orbital plane (view size
Mach number, slice through orbital plane.
Mach number, in frame of particle 2 corotating with particle orbit, perpendicular to orbital plane (view from particle 1) (view size
Mach number, perpendicular to orbital plane.
Phi-component of velocity about particle 1, slice through orbital plane, in the frame corotating with the particle 1-particle 2 orbit, with vectors showing the same quantity. Units are km/s. In this movie, particle 1 is to the LEFT of center (view size Phi-component of velocity about particle 1, in frame corotating with particle 2
).Phi-component of velocity about particle 1, slice perpendicular to orbital plane (view from particle 1), in the frame corotating with the particle 1-particle 2 orbit, with vectors showing the same quantity. Units are km/s (view size Phi-component of velocity about particle 1, in frame corotating with particle 2
).Phi-component of velocity about particle 1, slice through orbital plane, in the frame corotating with the particle 1-particle 2 orbit, with vectors showing the same quantity, now NORMALIZED to the velocity Normalized phi-component of velocity about particle 1, in frame corotating with particle 2
computed using the initial RGB profile. Note the changing colorbar and zoom level.Next steps
- Understanding 1D plots in terms of 2D movies (in progress, for discussion at meeting)
- Comparison with MacLeod+17b (in progress, for discussion at meeting)
- Comparison with Rheichardt+19 (in progress, for discussion at meeting)
- Comparison with Ricker+Taam12 (in progress, for discussion at meeting)
- Understand better the BHL theory/Dodd+McCrea theory and adjust analytic formulae if necessary, also estimate factor.
- Redo movies using full resolution data (after deciding which ones are most important)
- Movies for runs 149 and 151
- Post-processing of simulations to extract forces: I'm re-doing it with a reservation on bluehive to ensure that the results are correct as there had been some issues when I had done this the first time without a reservation (in progress)
- Movies and analysis for run 132 (subgrid accretion run) (partly complete)
- Calculate force from loss of angular momentum
abc
abc
Update 05/13
PN paper
draft updated on 05/05: figures paper draft
CE Jet
Last time we got the jet turned on, and added tracer to it.
Should look at the initialization of temperature, and measure the boundary flows.
Eric suggested starting jet from RLOF.
Others
- Start working on the Toro book.
- Radiation transport in CE.
- MESA, and plan for summer travel.
Update 5/13
Radiation Pressure Paper
All changes received so far made. I'm ready to submit.
Also, we still need to upgrade the Overleaf account to enable sharing again.
Charge Exchange
Running successfully on Stampede. On 20 nodes (960 cores), predicted time to completion is currently 1.8 months (got 2 frames in 48 hours). I expect this will speed up significantly once the stellar wind is more completely introduced and we reach closer to steady state.
Cost to allocation would be 1.8 mo x 30 days/mo x 24 hr/day x 20 nodes = 26000 node-hours. Much higher than I'd like.
Current frame
We also need to decide on the high-med-low for the stellar wind. Here are the current parameters (taken from John's simulation, except the mass loss rate, which is the solar mass loss rate - which corresponds to a 4x higher wind density than their strongest wind):
Parameter | Value |
cm | |
K | |
cm/s | |
g/s | |
AMR line transfer
Final functional fixes made (except a couple relating to subcycling). Some optimizations remaining, but could start using it for runs. Still want to do timing comparisons between old and new routines, and would be good to get a baseline for comparison of various optimizations.