​density-vz-and-Temperature-comparison

New Work

1. Continuing high res AGB run

AGB Paper

New notes

​Notes on AGB high resolution run 177, comparison with low resolution run 164.

Radiation Pressure Paper

Resubmitted. Hopefully not too long until we hear back.

Ray tracing

Almost working perfectly. Getting a few unexpected zero-length rays that I need to investigate, and a few strange rays toward the end.

Charge Exchange

Appears to be running well. Trying to make a movie of first six frames.

Magnetic Fields

Still having both pressure and divergence issues in the initialization. May need a brainstorming session for this one. See Mathematica notebook​ with equations to satisfy.

AMR profiling

Still waiting on Stampede.

New Work

Interpolated the OPAL EOS with the provided interpolation function. The result agrees with MESA in the region applicable to OPAL EOS.

Pressure
Internal energy

Next step

1. Now I am using a bash script that coordinates a Python script and two fortran scripts to do the interpolation, which is very non-efficient (computation time ~5 hours). Modifying the fortran script to get a one-step function may save a lot of time (It looks like MESA only use less then 1 second to do the calculations).
2. Inverting the tables to get other quantities as a function of density and internal energy.

AAS

Booked hotel for 3rd through 9th. No longer sure where I'm flying out of, but should hopefully know soon.

Radiation Pressure

Submitting Friday. Haven't gotten Ruth's comments yet.

Charge Exchange

See ​last week's blog post for test results.

Test comparison

Production run with John's "moderate" stellar wind is queued on Stampede2. Also queued non-rotating parameter space runs on BlueHive.

Magnetic fields

Added option for initial stellar B field that excludes planet (as if it were a perfect conductor). Getting divergence warning immediately after initialization. I believe the field should be divergence-free, but will this be a problem?

AMR profiling

New test still queued on Stampede.

Ray tracing

Tested patch output first. That's working well. Now just need to figure out best way to read mixed data types in Matlab, and should be able to plot rays without too much issue. Using 8-core, 3-level, 2-D 16x256 base resolution test mesh, which should leave enough room for each ray to be distinct even on the finest level.

Model Output vs Global (True) Values

Population and pCO2 Data From: Frank, Adam, and Woodruff Sullivan. "Sustainability and the astrobiological perspective: Framing human futures in a planetary context." Anthropocene 5 (2014): 32-41.
Temperature Data From: ​https://data.giss.nasa.gov/gistemp/graphs/graph_data/Global_Mean_Estimates_based_on_Land_and_Ocean_Data/graph.txt
How I Converted:

1. Got temperature anomaly data in Celsius from NASA Goddard's website (both smoothed and raw)
2. Converted the smoothed anomaly data into smoothed global temperature dat (in celsius) using Tony Del Genio's advice: People only report global temperature anomalies because surface weather stations only cover a fraction of the globe, and thus various assumptions are made by different groups to fill in all the places in the world with no data, so the actual temperature of Earth is only known to within a degree or so. People usually define the anomaly relative to the 1951-1980 reference period, whose mean anomaly should be close to zero. If you add a number like 14.5 deg. C to that, that should be as close to a real temperature as you can get. (will post his talk at the bottom of this post)
3. Then converted from celsius to kelvin (+273.15 deg. C), to compare that data to our model's output

Visualization of the Anthropocene With Actual Global Data

The above plot has been made using the data from Frank and Sullivan, 2014. From this plot, and with the assumption that the Anthropocene started at around 1945, it seems safe to define the Anthropocene physically as the point at which the relationship between human population and their carbon footprint switched from linear to exponential.

Tony's Talk: Albedos, Temperatures, and Habitability of Rocky Exoplanets

​Link to Tony's Talk

Here are the first three steps of charge exchange for a simulation where the densities of stellar and planetary material differ by five orders of magnitude.

 q(iHIIhot) =    55.0433622655540     
 q(iHcold) =    5504281.18319314     
 q(iHIIcold) =   0.000000000000000E+000
 q(iHhot) =   0.000000000000000E+000
 After step:
 q(iHIIhot) =    54.7826014982936     
 q(iHcold) =    5478205.36722786     
 q(iHIIcold) =   0.000000000000000E+000
 q(iHhot) =   0.000000000000000E+000
 Charge exchange coefficient is   3.455999999999997E-003
 Before iHIIhot protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 After iHIIhot protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 Before iHhot protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 After iHhot protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 Before iHIIcold protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 After iHIIcold protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 Before iHcold protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 After iHcold protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 Advanced level  0 to tnext= 0.6000E-08 with dt= 0.6000E-08 CFL= 0.5764E-07 max speed= 0.1921E+02 radiation sub-cycles=       1
 Info allocations    =   148.7 kb  148.7 kb
 message allocations =   ------    ------  
 sweep allocations   =   ------    122.5 kb
 filling fractions   = 
 Current efficiency  =  75%  16%  91% 
 Cell updates/second =        211       209 100%
 Wall Time Remaining =     9.1 day at frame    0.0 of     10
 AMR Speed-Up Factor =       0.2004E+03
 Before step:
 q(iHIIhot) =    54.7752475699149     
 q(iHcold) =    5478092.25506669     
 q(iHIIcold) =   6.222857684820517E-003
 q(iHhot) =   6.222857684820517E-003
 After step:
 q(iHIIhot) =    54.7751651392209     
 q(iHcold) =    5478084.01114326     
 q(iHIIcold) =   6.222848320107427E-003
 q(iHhot) =   6.222848320107427E-003
 Charge exchange coefficient is   3.455999999999997E-003
 Before iHIIhot protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 After iHIIhot protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 Before iHhot protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 After iHhot protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 Before iHIIcold protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 After iHIIcold protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 Before iHcold protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 After iHcold protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 Advanced level  0 to tnext= 0.3124E-04 with dt= 0.3124E-04 CFL= 0.3001E-03 max speed= 0.1921E+02 radiation sub-cycles=       1
 Info allocations    =   148.7 kb  148.7 kb
 message allocations =   ------    ------  
 sweep allocations   =   ------    122.5 kb
 filling fractions   = 
 Current efficiency  =  77%  14%  90% 
 Cell updates/second =        303       302 100%
 Wall Time Remaining =     1.8 min at frame    0.0 of     10
 AMR Speed-Up Factor =       0.2879E+03
 Before step:
 q(iHIIhot) =    22.3818513598937     
 q(iHcold) =    5478051.58190027     
 q(iHIIcold) =    32.3995362683480     
 q(iHhot) =    32.3995362683480     

----

 After step:
 q(iHIIhot) =    22.3818513491939     
 q(iHcold) =    5478051.57928144     
 q(iHIIcold) =    32.3995362528591     
 q(iHhot) =    32.3995362528591     
 Charge exchange coefficient is   3.455999999999997E-003
 Before iHIIhot protection:
 q(iHIIhot) =   -4.07465558395285     
 q(iHcold) =    5478025.12276545     
 q(iHIIcold) =    58.8560431859153     
 q(iHhot) =    58.8560431859153</span>

----

 After iHIIhot protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     

----

 Before iHhot protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     
 After iHhot protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     
 Before iHIIcold protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     
 After iHIIcold protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     
 Before iHcold protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     
 After iHcold protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624
 Advanced level  0 to tnext= 0.9368E-04 with dt= 0.6244E-04...

Jonathan was right on the money with the timestepping being too long. I've broken out the relevant blocks above. For the charge exchange step that follows the first line, we get an overall \Delta of

\Delta = \beta(HIIhot \cdot Hcold - Hhot \cdot HIIcold)dt = 3.456\times10^{-3}(22.382*5478052-32.4^2)6.244\times10^{-5} \approx 26

Since HIIhot is only ~22, this gives us HIIhot of ~-4 after charge exchange is complete.

This run is from after I implemented protections for any species less than 0. The section that follows the first line shows that 4 is added back to HIIhot and Hcold, and removed from HIIcold and Hhot, reversing the excess charge exchange that brought HIIhot negative.

This method of protection can result in large oscillations arount equilibrium. Better would be to prevent any one species from changing more than e.g. 10%, as we do with the other quantities, but this is more computationally expensive thanks to the amount of subcycling that would be required in some places (namely those with large \Delta in comparison to any one species).

AAS

Registered. Have flight info ready, will get hotel once I have flight confirmed.

Abstract: Charge exchange? Should have results by January. Poster, iPoster, or iPoster-Plus?

Radiation Pressure

Ruth said minor comments are on the way. Haven't heard anything from Eric.

Charge Exchange

Last week's problem with no hot neutrals appearing was exactly as expected. That's been fixed. But in the second frame I'm still ending up with stellar material right by the planet, for some reason (which is what the changes that caused the previous issue were prompted by…). I've added some debug statements - if that doesn't work, may be time for another small test problem, though it seems odd that it happens only near the planet - suggests to me that it's a problem with the simulation, or possibly some other multispecies routine (species protections? These should preserve 0s, though), rather than the charge exchange routines. In particular, the fact that there's no problem for an entire frame suggests it's something about the changes caused by the stellar wind.

First frame

Second frame

Freezes on BlueHive when compiled with optimization level 3, before reading in from restart file. Runs fine when compiled with optimization level 2. Should probably figure out why at some point, but for now I'm letting them run.

AMR Profiling

Trying to run on Stampede. Minor problem with profiler, trying suggested workarounds now.

Ray Tracing

Queued on BlueHive. One job should start on Friday; also requested job with shorter time allocation, to see if it might start earlier (don't need long).

Code Location on GitHub

​https://github.com/ethansavitch/coupled_ebm

Jacob's Answers To My Questions

1. Purpose of the Stochastic Parameter?

The stochastic parameter in the EBM equation was used in my glacial-interglacial cycles paper (cited on your poster), which was done to explore the concept of "stochastic resonance" in a climate system. The idea is that the Earth system has multiple stable climate states (ice free, large ice cap, small ice cap, ice covered), and that variability in the Earth system itself might contribute to periodic transitions between states (which for small random variability would be between the large and small ice cap solutions). The stochastic parameter thus is a hand-wavy way of approximating the effect of sub-model-scale processes that would exert a random noise upon equilibrium climate. The net effect of this is that the hysteresis loop decreases in width, since it becomes slightly easier to access different states. There's no need to include the stochastic parameter for accuracy. I included a longer explanation here in case it is useful at a later stage in the project.

2. What's the purpose of the soladj=solar forcing adjustment parameter?

The solar forcing adjustment is used to time step the model through a prescribed increase in stellar brightness as the star ages through the main sequence. If you look near line 265 of the code, you will see that when soladj is true, the model reads in the eccentricity, precession, obliquity, and solar constant from the file data/insolaout.dat. This data is from Laskar's astronomical solutions for Earth's past and future (​http://vo.imcce.fr/insola/earth/online/earth/earth.html). If soladj is false, then these values are all treated as fixed with time. You only need the solar forcing adjustment feature if you are interested in long time integrations with realistic representation of solar brightening and orbital evolution.

3. How best to incorporate weathering?

If you want CO2 weathering as described by Williams & Kasting (Icarus, 1997), then you could uncomment line 637, where I have added a comment "removed co2 adjusting.". However, the timescales we are talking about are on the order of thousands to tens of thousands of years, where shorter-term variation in atmospheric CO2 will dominate. The carbonate-silicate cycle operates with a residence time of about half a million years, and only matters for Earth climate when we consider timescales into the hundreds of thousands of years (see, e.g., ​Berner, Lasaga, & Garrels, Am. J. Sci., 1983). So I would say leave the code as-is, and don't worry about carbonate-silicate weathering. If you do want to experiment with the effects of weathering, then the solution in 3b of your blog post sounds like a good approach to me. (constant annual per-capita reduction in CO2, parameterized with \omega)

ToDo Next

1. Clean/Fix Jupyter Notebook. All parts are currently working except for the namelist functionality.
2. Compare output of coupled EBM with global average temperatures.
3. Increase accuracy of the growth rates temperature dependencies, which will hopefully make our output more in line with the global averages.

Attempts on OPAL EOS and MESA ​https://www.pas.rochester.edu/~yishengtu/research_files/CEE_gp_meeting_09042019/

CE Jet - Stampede run 1

Setups

• 73 frames done.
• Mass-loss-rate = 2 Msun/year
• Jet_radius = 2.25 Rsun (16 cells)
• working on restart on Stampede

Separation compares to fiducial run

Movies

​Density-edgeon

​Density-faceon

​Compare-Jet-fiducial-edgeon

​Compare-Jet-fiducial-faceon

Things to do

questions:

• how does the jet evolve?

○ quenched? ○ when, and why ○ look at tracer, zoom in of the jet, force balance, ram pressure (balance with gravity and thermal pressure)

• does the ambient matter?

○ compare run with lower ambient (another run)

• does the launch time matter?

○ start from farther out (another run) ○ and a fiducial (no jet) run to compare

• what happens at later times?

runs

1. continue this run, 2 Msun/yr jet
2. a lower ambient run
3. larger separation, with low ambient
4. no jet, larger separation, low ambient

• AGB?

PN paper

Paper is coming along. ​paper I spent some time fixing the calculations for section 4.3. ​notes Need to finish the conclusion.

Movies

​production runs

​tests_without_IW

New Work

1. Ran test runs + analysis + comparison with fiducial run, for AGB CE simulation

AGB Paper

New notes

​Notes on AGB test runs.

Movies

Pressure profile evolution, inner 0.6 Rsun (Note Rsoft = 2.4 Rsun)
​Run 164
​Run 164 frames 0-1
​Run 170
​Run 170 frames 0-1
​Run 173
​Run 173 frames 0-1
​Run 174

Next steps

• 2-3 more test runs to try to improve energy conservation further
• Tracer? Ambient, AGB core (e.g. to track energy)

Accumulation of magnetic field could have been observed at the grids when the AMR level changed. After 3 weeks, we discovered 2 indexing errors.

Images: ​http://www.pas.rochester.edu/~aanand/plots/?C=M;O=D

AAS

Need to register by 9/25, abstracts by 10/8 (if we want to submit something). I'm not currently an AAS member - not sure what the turn-around time is to become one, but it's a significant difference in registration fee.

Flights are ~$1200, similar for (AAS) hotel, though looks like about 2/3 price is available.

Radiation Pressure

Sent around what should be the nearly final version of the ​radiation pressure paper along with the response to the referee. I'll plan to submit as soon as I have a response from everyone (hopefully end of the week or beginning of next).

Charge Exchange

I've removed everything that I believe can eliminate the added hot neutrals from charge exchange (ionizing radiation, reverse charge exchange, collisional ionization [not turned on originally]), and checked that the charge exchange rate is being calculated correctly in the actual line transfer step (rather than in diagnostics). We have rates up to order 10 in some places, so with a dt of 0.1 over a frame, we should end up with some significant amount of hot neutrals after a frame now. Still getting none, however. Feel like I have to be missing something stupid now.

AMR profiling

Stampede has the Intel profiler available. Haven't had a chance to use it yet, but as soon as the charge exchange is working I'll check it out.

Non-rotating charge exchange

Feel like this should wait until we think we have charge exchange working, unless there's some reason it might give a clue to what's going wrong? I can't think of any, though. Ready and waiting, though.