Posts for the month of October 2019

CE Jet status

model A, run002, at 46 frames

Starting with P2 at the edge of the RG envelope.

Gets to the green dot below in the binary separation plot. The next job is in queue on stampede.

Evolution movies: rho-edge-on (with-v-vec) rho-face-on (with-v-vec) Temperature z-speed

Run002-model-A-Movies-and-plots

http://www.pas.rochester.edu/~yzou5/CE/Blogpost_20191029/compare_3_stp_runs_to_143.png

model C, run 003, at 11 frames

Starting with P2 at 1.5 times RG radius.

I might pause the run and fix the following problems first:

  • RG gets unstable around the edge of its envelope.
  • Jets' temperature has a sharp change aligned with the boundary between level 3 to level 4 refinement area. Seems like pressure protection, but the actual temperature is 1000 times higher (~1e-6 K) than the minTemp set in physics.data (1e-10 K).
  • May be a bug in the particle buffer code I haven't found out.

Evolution movies: rho-edge-on (with-v-vec) rho-face-on (with-v-vec) Temperature z-speed

Run002-model-B-Movies-and-plots

http://www.pas.rochester.edu/~yzou5/CE/Blogpost_20191029/separation_run003_model_C.png

Update 10/28

MHD

We've found the bug, in upwinded_emf. I'll be checking the implementation from Gardiner & Stone 2005, but for now it runs.

We now have a more physical issue, with infall from around the AMR boundaries which I believe is due to (magnetic) pressure gradients in the ambient material. Setting the total (thermal+magnetic) pressure to a constant equal to the far-field value creates more problems (pressure protections in the initialization), but I have one more idea to try (in addition to trying a different interpolation method, to see if the boundaries will smooth out that way). Jonathan pointed out there's also magnetic tension - I'm hoping that countering the pressure will be sufficient to make things stable long enough to start the wind.

http://www.pas.rochester.edu/~adebrech/MHD/infall.png

Charge Exchange

http://www.pas.rochester.edu/~adebrech/ChargeExchange/HD209458b_charge_exchange_species0018.png

http://www.pas.rochester.edu/~adebrech/ChargeExchange/HD209458b_charge_exchange0018.png

Also started a low-mass-loss-rate stellar wind run. Will check and see how that one's running soon.

COMMON ENVELOPE SIMULATIONS

New Work

  • Advanced Run 183 (restarted at frame 194 = 44.9 days of Run 177 now with extra refinement around particle 2)
  • With JC coded in new routines that conserve total gas energy including PE due to self-gravity but not including particles. These new routines have not yet been tested.

Results

  • Rate of energy gain after first periastron passage is much smaller in Run 183 compared to 177, as expected.

Notes

Notes on AGB high resolution run 183.

Next steps

  • Continue Run 183
  • Plot results from Run 185 (convergence test)
  • Further modify new energy conserving routines to take into account particle-gas and particle-particle energy terms
  • Further automate and improve simulation:
    • predefine refinement radius as a function of time (should vary smoothly rather than in sudden steps)
    • center of refinement region should transition smoothly in time from particle 1 to particle CM (instead of suddenly)
    • experiment with using true AMR, e.g. refine on density or pressure gradients within refinement shape

Update 10/21 - CEJet

Two production runs on Stampede

Run 002, model A

Jet-mass-loss-rate = 2 Msun/year, higher ambient profile (rho0 = 6.69e-9, p0 = 1e5), a0 = 49 Rsun.

Stopped at frame 28, waiting for the next queue.

Run 003, model C

Jet-mass-loss-rate = 2 Msun/year, lower ambient profile (rho0 = 1e-10, p0 = 2e3), a0 = 73.5 Rsun.

Currently running, has got to frame 5.


  • Each run is on 16 nodes (768 cores). Recent queue time on Stampede is about 3.5 days, so I submit jobs in a 3-day separation.
  • lSkipProfile = False for production quality. This takes time at every restart, (dt~13 hours between the first two chombos files). Later on each frames take about 2 - 1.5 hours, and is getting faster throughout the run time (so I didn't use the restarting at each 5 frames scheme).
  • The last run on Stampede was killed with the following message:

any ideas?

===================================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
=   PID 334217 RUNNING AT c498-002
=   EXIT CODE: 7
=   CLEANING UP REMAINING PROCESSES
=   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
===================================================================================
   Intel(R) MPI Library troubleshooting guide:
      https://software.intel.com/node/561764
===================================================================================
TACC:  MPI job exited with code: 7
TACC:  Shutdown complete. Exiting.
Starting replacement routine
/tmp/slurmd/job4506883/slurm_script: line 69: .true.: syntax error: operand expected (error token is ".true.")

Summary of the four models planned

  • A: 2 Msun/yr jet, higher ambient, starting at 49 Rsun.
  • B: 2 Msun/yr jet, lower ambient, starting at 49 Rsun.
  • C: 2 Msun/yr jet, lower ambient, starting at 73.5 Rsun.
  • D: no jet, lower ambient, starting at 73.5 Rsun.

Update 10/21

MHD

I think a CFL issue has been successfully ruled out (Jonathan, if you think a direct wave speed comparison would be useful, I can do that - but based on the tests already done, don't think they're necessary).

Is the prolongation of the aux fields still the most likely place for a bug? It may be time to sit down for a bug hunt.

Charge exchange

Still moving right along. Hoping we reach a steady state soon. Here's the latest frame. (At some point I'll get a movie, but waiting to hear back about the paper so we can hopefully archive the rad pressure run and make room….)

http://www.pas.rochester.edu/~adebrech/ChargeExchange/charge_exchange_species0017.png

http://www.pas.rochester.edu/~adebrech/ChargeExchange/charge_exchange0017.png

I meant to turn off the lower y boundary, turned off the upper one instead - been corrected, and shouldn't affect the simulation as it progresses (but I've queued the restart from an earlier frame anyway).

AMR line transfer

Set up a stack, still need to change ray tracing to use it. Think it should be a stack of rays rather than ray positions, which means the main ray tracing loop needs to be rejiggered again.

Once the MHD bug is worked out I should have time to go a little deeper into this, see if there's a better way to approximate the graph that represents the dependencies.

Dust in AstroBEAR - Update 2019/10/21

Current Objective

  • Debugging the dust drag advection routines (collisional drag)
  • Finishing the dust sputtering routine
  • Finishing the total mass calculation routine

Current Status

Drag Advection

(See also report from last week https://www.pas.rochester.edu/astrobear/blog/fschmidt10142019)

  • dqdt for the dust velocity is way too large (in the order of 1e42) which (apart from being quite obviously incorrect) causes an issue with the RK integrators
  • I found a couple of bugs already (e.g. there was an issue with the transformation into the dust grain rest frame) but it's still giving me values that are too big. The frame of reference is changed multiple times throughout the calculation so I suspect more errors in there somewhere. Currently still debugging.

Sputtering Routine

  • Started implementing the sputtering routines based on Dwek & Arendt (1992), Tielens et al. (1994), and Nozawa et al. (2006). This is a work in progress.

Total Mass Calculation

  • Did some more work on this internal processing routine which calculates the total mass of dust in each bin and across all the bins for each timestep (not just each output timestep) and stores them in a file. This additional dust-related output can be turned on or off and will mainly be used for 3D simulations (in case we don't get a new HPC award next year; our local facilities can't handle many AMR 3D Chombo output files. So using this output option would allow us to cut down on the Chombo output while still getting dust mass information with high temporal resolution). Still a work in progress but should be done soon-ish.

Questions

Feedback on Gas

Accelerating the dust should remove energy from the system and growing or destroying grains should either generate or remove gas mass from the system. I've been thinking about how best to implement that:

  • For the grain destruction or grain growth, the most obvious choice would be to simply add or remove material to/from the gas during the source step. So

  • For the gas drag the defining quantity that's calculated is the force the dust grains experience. is then used to calculate the acceleration that each grain experiences. My current approach is to set , where is the total force experienced by the gas particles in one cell. The force which a single gas particles experiences is then and that then results in a declaration. All of that is currently calculated in the source term routine alongside the dust acceleration. Questions: Will this result in energy conservation problems?

Stiffness of Equations

I'm very certain that the current issue with the RK integrator error is down to a bug and not down to an intrinsic issue with the stiffness of the drag equations (since by definition the acceleration of the dust grains cannot exceed the acceleration of the gas particles). However, it is, in theory, possible that the destruction routines are very stiff under certain conditions (very high densities, very high velocities etc.). Question: I was wondering what could be done in that case?

Planet-Civ Evolution Project Update 10/21

Update on Magnetization of Asteroids by Solar Winds

Ran the Oran resistivity structure using our solar wind parameter. Both start with uniform B = 100 nT in y direction, constant density and field velocity ramping up.

Magnetic field pile up appears to nearly double for our solar wind conditions.

Our condition (300 particles/cc, 500 km/s, 1e6 K) has 5.7x Final B, 2 different planes are shown:

http://www.pas.rochester.edu/~aanand/movies/xy_a3cl.gif

http://www.pas.rochester.edu/~aanand/movies/xz_a3cl.gif

Oran condition (35 particles/cc, 700 km/s, 50,000 K) has 3.4x Final B:

http://www.pas.rochester.edu/~aanand/movies/or_1_1.7_a3cl_r_151.gif

http://www.pas.rochester.edu/~aanand/movies/oran_xz_a3cl.gif

Common parameters are listed here: http://www.pas.rochester.edu/~aanand/movies/problem.data

Currently need ~6 days of computation time to reach this state.

Moon Case

Ran Moon case using new boundary condition, and uniform field and particle density. Core does not diffuse out field till the last frame; need to start with 0 field inside the moon (work in progress, condition similar to Alex's case).

http://www.pas.rochester.edu/~aanand/movies/moon_0.2_1.7_150r.png

3D Viscosity

Extrapolated Code for 3D Viscosity. Yet to be tested, but should be able to add soon.

Parameters to Change

  1. Box height and width to be made double (y and z directions).
  1. Resolution studies.
  1. Artificially add relative magnetic permeability.
  1. Add Lineout plots through the center
  1. Run 6 initial cases, with both wind conditions (except for moon case), making 10 cases:
    1. Moon with core
    2. Moon without core
    3. Oran resistivity
    4. Conducting layer on top of case ©. Additional parameter needed.
    5. Uniform resistivity of value close to max resistivity of ©
    6. Case (e) with an outer conducting shell
    7. Case (f) with Front side conducting, back side resistive (like a comet).

Base parameter values

Critical Resistivity (Comp. Units): 6.12e6

Resistivity (fraction of critical):

  1. Moon with core: 1
  2. Moon without core:
  3. Oran resistivity:

a3cl: 1.63e-4 to 1.62e-2

a1ht: 2.45e-4 to 1.96e-3

  1. Conducting layer on top of case ©: Shell of 1.96e-4
  2. Uniform resistivity of value close to max resistivity of ©: 2e-3
  3. Case (e) with an outer conducting shell: Body - 2e-3 and Shell - 2e-4
  4. Case (f) with Front side conducting, back side resistive (like a comet).

Plots Needed

  1. Peak Magnification vs Effective Resistivity (with M*√γβ)

COMMON ENVELOPE SIMULATIONS

New Work

  1. Advanced Run 177 (AGB run with increased resolution at AMR level 5 around particle 1 compared to the old AGB run, Run 164) and found that after the first periastron passage, the rate of energy increase becomes dramatically higher (going from 10% to 60% of that of Run 164)
  2. The cause of this is likely that the resolution around particle 2 is not high enough to resolve the accreted mass around particle 2. Therefore, I have submitted two new runs which put the same extra resolution around particle 2 as well (each run has a different restart time, both restart times are before first periastron passage). I've also submitted a third run which does NOT change the center of the AMR level 3 spherical refinement zone from particle 1 to the particle CM just after periastron passage, to make sure that this was not the cause of the sudden increase in the rate of energy gain in Run 177. Waiting for Stampede queue.
  3. I read this paper and discussed with Jonathan. It gives a way to conserve energy explicitly including potential energy from self-gravity. We think that we could take it a step further and include the PE due to particle-gas interaction as well. This should not be too hard to implement.

Updated Notes

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

Update 10/14

MHD

Believe current problems are related to prolongation of the aux fields. Because we lose information when going from face-centered to cell-centered B fields, we can't just reconstruct the "correct" (uncorrected) aux fields from q. I believe this is the relevant piece of code, in ProlongateParentData:

DO i=1,nDim
  DO j=0,r**nDim-1
    IF (MOD(j/r**(i-1),r)==1) CYCLE
    DO n=1,nDim
      l(n)=MOD(j/r**(n-1),r)  
    END DO
    ic(1:nDim,1)=1+l(1:nDim)-rmbc; ic(1:nDim,2)=Info%mX(1:nDim)-r+1+l(1:nDim)+rmbc;ic(i,2)=ic(i,2)+2
    ip(1:nDim,1)=mB(1:nDim,1)-mbc; ip(1:nDim,2)=mB(1:nDim,2)+mbc;ip(i,2)=ip(i,2)+1
    Info%aux(ic(1,1):ic(1,2):r,ic(2,1):ic(2,2):r,ic(3,1):ic(3,2):r,i) =     &
    Parent%auxChild(ip(1,1):ip(1,2)  ,ip(2,1):ip(2,2)  ,ip(3,1):ip(3,2)  ,i)
  END DO
END DO

Haven't looked to see if there are any obvious errors yet. Outer ghost zones were the source of the error in last test.

AMR line transfer

Thinking some more about prioritization. Setting up a stack isn't a bad idea - may be the next thing I try, since it should be easy.

In the truly ideal case, we'd know which processors each ray depended on, and we'd prioritize based on number of dependencies. There may be a better way to approximate this than a stack. I'll continue thinking.

Charge Exchange

Still progressing. Here's the latest frame:

http://www.pas.rochester.edu/~adebrech/ChargeExchange/charge_exchange_species0014.png

http://www.pas.rochester.edu/~adebrech/ChargeExchange/charge_exchange0014.png

Dust in AstroBEAR - Update 2019/10/14

(Can't make it to today's meeting due to a time conflict with a seminar I am required to attend. Back next week!)

Current Objective

  • Developing the dust drag advection routines (collisional drag)

Current Status

  • Collisional drag subroutine is implemented though simulations crash with the following error message: "Src Error::RKOrder45 — timestep insignificant. Halting." I'm currently checking whether this is a bug in the source routine I've set up resulting in an unrealistically rapid dust acceleration or whether the system is too stiff.

Supporting Material

Gas Drag Equations (for Collisional Drag)

The new dust velocity after one time-setp is given by:

where

The drag force is given by

The sum over covers all particles species in the gas (we are using an average value for an oxygen dominated SNR). is given by:

where is

The remaining variables are (all in cgs):

  • : mass of a dust grain
  • : number density of gas species
  • : gas temperatire
  • : grain radius
  • : mass of a gas particle of species
  • : absolute value of the relative velocity between gas and dust
  • : timestep

The resulting drag force is then multiplied with a scaled unit vector to extract the drag forces in each dimension.

(Refernces: Kirchschlager et al. 2019, Baines et al. 1965, Draine & Salpeter 1979)

Coupled EBM Project Update 10/13

Goal

Make plots for the equilibrium temperatures as a function of pCO2 and orbital distance.
(Equilibrium Temperature = temperature at which model balances incoming solar radiation with outgoing long-wave terrestrial radiation, without coupling)
http://www.pas.rochester.edu/~esavitch/plots/pco2Vslogtemp.png http://www.pas.rochester.edu/~esavitch/plots/distanceVStemp.png

Visualizing How our Model Behaves at Different Orbital Radii

http://www.pas.rochester.edu/~esavitch/plots/diffDist.gif http://www.pas.rochester.edu/~esavitch/plots/diffDistSLO.gif

Dust in AstroBEAR - Update 2019/10/07

Passive Advection of Dust

Last Status

  • passive advection was working if the gas-to-dust density ratio was 1. For any other ratios, the flux calculations were way off resulting in a faulty advection and eventually in a termination of the simulation run.

Current Status

  • passive advection is now working for all gas-to-dust density ratios and bin numbers. The bug was in the conserved ↔ primitive conversion routines of EOS.f90. For passive advection the conversion routine should divide by or multiply with the gas density not the dust density.

Test Simulation: https://youtu.be/KJW2fLnDyG0

Next Steps

  • Processing subroutine to store total dust mass in domain for each bin → Put in Sweep After Step in sweep_scheme.f90?
  • Add other advection schemes (gas drag & plasma drag, starting with gas drag) → Put call to advection routines in Sweep Before Step in sweep_scheme.f90 but the actual calculations will be in dust.f90

Equqations

  • Conservation of Gas Momentum:

where

  • Conservation of dust momentum:

where

Update 10/7

AAS

Current version of the abstract:

The photoionization-driven evaporation of planetary atmospheres has emerged as a potentially fundamental process for planets on short period orbits. Using AstroBEAR, an AMR multiphysics code, in a co-rotating frame centered on the planet, we model in three dimensions the transfer of ionizing photons into the atmosphere, the subsequent launch of the wind and the wind's large scale evolution subject to tidal and non-inertial forces. We run simulations for planets of 0.263 and 0.07 Jupiter masses and stellar fluxes of 2x1013 and 2x1014 photons/cm2/s. These simulations reveal new, potentially observable planetary wind flow patterns, including the development, in some cases, of an extended neutral tail lagging behind the planet in its orbit. In addition, the role of radiation pressure in shaping exoplanet photoevaporation remains a topic of contention. Radiation pressure from the exoplanet's host star has been proposed as a mechanism to drive the escaping atmosphere into a "cometary" tail and explain the high velocities observed in systems where mass loss is occurring. Using simulations of a planet modeled after HD 209458b and subject to both ionizing and Lyman-α radiation, we demonstrate that for the Lyman-α flux expected for HD 209458, radiation pressure is unlikely to significantly affect photoevaporative winds or to explain the high velocities at which wind material is observed. Charge exchange between the stellar and planetary winds has also been suggested as a method for creating the observed Lyman-α absorption signature. We present results of simulations that explore the effect of charge exchange on our synthetic observations. In addition, we present simulations of the effect of stellar and planetary magnetic fields on our self-consistently launched wind.

MHD

Pressure protections around the planet. Some very small B fields in the planet - may need to fix the field to zero inside the planet at every step?

Also realized that I set the field to zero inside Rp only, not inside Rob. That may help some. Perhaps it's even worth fixing the field to zero out to within a small but significant distance from Rob?

http://www.pas.rochester.edu/~adebrech/MHD/pressure_protections.png

Also plan to try fixed-boundary global sim when I get a chance.

Charge Exchange

Still running fairly slowly. I expect it should reach steady state soon, though (assuming we don't enter a cycle with the stellar wind - I am using John's intermediate case).

http://www.pas.rochester.edu/~adebrech/ChargeExchange/charge_exchange_species0010.png

The wind is turning inward as it approaches the outside boundary - I think this is expected due to Coriolis/centrifugal force/differential rotation? But haven't really had time to investigate yet.

http://www.pas.rochester.edu/~adebrech/ChargeExchange/charge_exchange_species0010.png

AMR profiling

Actual profiling should finally run successfully on Stampede. Also on my to-do list is to make sure we're not waiting longer than necessary for local/MPI rays.

Coupled EBM Project Update 10/4

Equations As of Now




  • T = Temperature
  • P = pCO2
    • = per-capita carbon footprint
  • N = Population
    • = per-capita growth rate, dependent on temperature
    • =temperature resulting in optimum growth rate
    • = exponential ramping distance
    • = proxy for technology, kept constant here

http://www.pas.rochester.edu/~esavitch/plots/habitableZone/exInput.png

Methodology

  1. Run energy balance model until equilibrium is reached, with an average global temperature of
  2. Center the optimum growth rate for the civilization to that temperature, ie: set
  3. Turn on coupling functionality, which allows population to rise, causing more CO2 to be in the atmosphere, which inevitably raises global temperatures past a habitable point

Example outputs for Varying Orbital Distances

http://www.pas.rochester.edu/~esavitch/plots/habitableZone/0-94AU.png http://www.pas.rochester.edu/~esavitch/plots/habitableZone/0-95AU.png http://www.pas.rochester.edu/~esavitch/plots/habitableZone/0-98AU.png http://www.pas.rochester.edu/~esavitch/plots/habitableZone/1AU.png http://www.pas.rochester.edu/~esavitch/plots/habitableZone/1.1AU.png http://www.pas.rochester.edu/~esavitch/plots/habitableZone/1.3AU.png http://www.pas.rochester.edu/~esavitch/plots/habitableZone/1.524AU.png http://www.pas.rochester.edu/~esavitch/plots/habitableZone/1.7AU.png http://www.pas.rochester.edu/~esavitch/plots/habitableZone/1.8AU.png

notes

  1. 1.524AU is the orbital distance of mars
  2. 273.15K is the freezing temperature of water
  3. 373.15K is the boiling temperature of water