Binder

This notebook teaches how to:
update groups and fields and use heartbeats.

In addition to the simframe requirements, the following packages are needed for this notebook: astropy, matplotlib.

3. Updating Groups and Fields

In this example we revisit the orbital integration from the previous tutorial. But this time we also want to calculate the equilibrium temperature on Earth. That is not a quantity that we need to integrate. We can simply calculate it from the other quantities. But we have to tell simframe how to do it.

Setting up the problem

The problem setup follows the previous example. Please have a look there for more details.

[1]:
from simframe import Frame
[2]:
sim = Frame(description="Earth-Sun system")
[3]:
sim.addgroup("Sun")
sim.addgroup("Earth")
[4]:
AU       = 1.495978707e11         # Astronomical unit [m]
day      = 8.64e4                 # Day [s]
G        = 6.6743e-11             # Gravitational constant [m³/kg/s²]
sigma_sb = 5.6703744191844314e-08 # Stefan-Boltzmann constant [kg/s³/K^4]
year     = 3.15576e7              # Year [s]

M_earth  = 5.972167867791379e24   # Mass of the Earth [kg]
M_sun    = 1.988409870698051e30   # Mass of the Sun [kg]
R_earth  = 6378100.               # Radius of the Earth [m]
R_sun    = 695700000.             # Radius of the Sun [m]

Earth

[5]:
import numpy as np
[6]:
sim.Earth.addfield("M", M_earth, description="Mass [kg]", constant=True)
sim.Earth.addfield("r", np.zeros(3), description="Position [m]")
sim.Earth.addfield("v", np.zeros(3), description="Velocity [m/s]")

Sun

[7]:
sim.Sun.addfield("M", M_sun, description="Mass [kg]", constant=True)
sim.Sun.addfield("r", np.zeros(3), description="Position [m]")
sim.Sun.addfield("v", np.zeros(3), description="Velocity [m/s]")

Integration variable

[8]:
sim.addintegrationvariable("t", 0., description="Time [s]")
[9]:
dt = 1.*day
[10]:
def f_dt(frame):
    return dt
[11]:
sim.t.updater = f_dt
[12]:
snapwidth = 5.*day
tmax = 2.*year
[13]:
sim.t.snapshots = np.arange(snapwidth, tmax+1, snapwidth)

Writer

[14]:
from simframe import writers
[15]:
sim.writer = writers.namespacewriter()
[16]:
sim.writer.verbosity = 0

Differential equations

[17]:
def dr_Earth(frame, x, Y):
    return frame.Earth.v

def dr_Sun(frame, x, Y):
    return frame.Sun.v
[18]:
# Gravitational acceleration
def ag(M, r, R):
    direction = r-R
    distance = np.linalg.norm(direction)
    return -G * M * direction / distance**3
[19]:
def dv_Earth(frame, x, Y):
    return ag(frame.Sun.M, frame.Earth.r, frame.Sun.r)

def dv_Sun(frame, x, Y):
    return ag(frame.Earth.M, frame.Sun.r, frame.Earth.r)
[20]:
sim.Earth.v.differentiator = dv_Earth
sim.Earth.r.differentiator = dr_Earth
sim.Sun.v.differentiator   = dv_Sun
sim.Sun.r.differentiator   = dr_Sun

Integrator

[21]:
from simframe import Integrator
[22]:
integrator = Integrator(sim.t, description="Leapfrog integratior")
[23]:
from simframe import schemes
from simframe.integration import Instruction

Leapfrog integrator

[24]:
instructions = [
    Instruction(schemes.expl_1_euler, sim.Sun.v,   fstep=0.5),
    Instruction(schemes.expl_1_euler, sim.Earth.v, fstep=0.5),
    Instruction(schemes.update,       sim.Sun.v             ),
    Instruction(schemes.update,       sim.Earth.v           ),
    Instruction(schemes.expl_1_euler, sim.Sun.r,   fstep=1.0),
    Instruction(schemes.expl_1_euler, sim.Earth.r, fstep=1.0),
    Instruction(schemes.update,       sim.Sun.r             ),
    Instruction(schemes.update,       sim.Earth.r           ),
    Instruction(schemes.expl_1_euler, sim.Sun.v,   fstep=0.5),
    Instruction(schemes.expl_1_euler, sim.Earth.v, fstep=0.5),
]
[25]:
integrator.instructions = instructions
[26]:
sim.integrator = integrator

Initial conditions

In this example we want to send the Earth on an eccentric orbit to have a seasonal change in the equilibrium temperature. For this we simply incline the initial velocities of Earth an Sun by an angle alpha.

[27]:
r_Earth_ini = np.array([AU, 0., 0.])
r_Sun_ini = np.zeros(3)
# Center of mass
COM_ini = (M_earth*r_Earth_ini + M_sun*r_Sun_ini) / (M_earth+M_sun)
# Offset both positions
r_Earth_ini -= COM_ini
r_Sun_ini -= COM_ini
[28]:
mu = M_earth*M_sun / (M_earth+M_sun)
[29]:
alpha = 30. * np.pi/180  # conversion from deg to rad
v_Earth_ini = np.sqrt(G*M_sun/M_earth*mu/AU)
v_Earth_ini = np.array([-np.sin(alpha)*v_Earth_ini, np.cos(alpha)*v_Earth_ini, 0.])
v_Sun_ini   = np.sqrt(G*M_earth/M_sun*mu/AU)
v_Sun_ini   = np.array([np.sin(alpha)*v_Sun_ini, -np.cos(alpha)*v_Sun_ini, 0.])
[30]:
sim.Earth.r = r_Earth_ini
sim.Earth.v = v_Earth_ini
sim.Sun.r   = r_Sun_ini
sim.Sun.v   = v_Sun_ini

Starting

[31]:
sim.run()
Execution time: 0:00:00

Reading and plotting

[32]:
data = sim.writer.read.all()
[33]:
import matplotlib.pyplot as plt

def plot_orbits(data):
    fig, ax = plt.subplots(dpi=150)
    ax.set_aspect(1)
    ax.axis("off")
    fig.set_facecolor("#000000")
    imax = data.t.shape[0]
    for i in range(imax):
        alpha = np.maximum(i/imax-0.1, 0.5)
        ax.plot(data.Sun.r[i, 0], data.Sun.r[i, 1], "o", c="#FFFF00", markersize=4, alpha=alpha)
        ax.plot(data.Earth.r[i, 0], data.Earth.r[i, 1], "o", c="#0000FF", markersize=1, alpha=alpha)
    ax.plot(data.Sun.r[-1, 0], data.Sun.r[-1, 1], "o", c="#FFFF00", markersize=16)
    ax.plot(data.Earth.r[-1, 0], data.Earth.r[-1, 1], "o", c="#0000FF", markersize=4)
    ax.set_xlim(-1.5*AU, 1.5*AU)
    ax.set_ylim(-1.5*AU, 1.5*AU)
    fig.tight_layout()
    plt.show()
[34]:
plot_orbits(data)
_images/3_updating_47_0.png

The Earth is now on an eccentric orbit around the sound. By the spacing of the snapshots, which are constant in time, you can even tell that the Earth is faster around the perihelion and slower around the aphelion.

Adding Temperature

To calculate the temperature Earth and Sun need additional fields for their temperatures and their radii. In additional to that, we give the Earth an additional field for the Earth-Sun distance.

We initialize the Earth’s temperature with zero and update it as soon as we have a function for it.

[35]:
sim.Earth.addfield("d", AU, description="Earth-Sun distance [m]")
sim.Earth.addfield("R", R_earth, description="Radius [m]")
sim.Earth.addfield("T", 0., description="Temperature [K]")
[36]:
sim.Earth
[36]:
Group
-----
    d            : Field (Earth-Sun distance [m])
    M            : Field (Mass [kg]), constant
    r            : Field (Position [m])
    R            : Field (Radius [m])
    T            : Field (Temperature [K])
    v            : Field (Velocity [m/s])
  -----
[37]:
sim.Sun.addfield("R", R_sun, description="Radius [m]", constant=True)
sim.Sun.addfield("T", 5778, description="Effective temperature [K]", constant=True)
[38]:
sim.Sun
[38]:
Group
-----
    M            : Field (Mass [kg]), constant
    r            : Field (Position [m])
    R            : Field (Radius [m]), constant
    T            : Field (Effective temperature [K]), constant
    v            : Field (Velocity [m/s])
  -----

Skipping Fields in Output Files

It is possible to not write certain fields into the output. This can be especially useful if you have large fields, that consume a significant amount of memory, that you don’t need for your data analysis. Simply the the save attribute to False.

[39]:
sim.Sun.T.save = False

Updating fields

The energy the sun is emitting per time, i.e., it’s luminosity \(L\), is given by

\(\Large L_\odot = 4\pi R_\odot^2 \sigma_\mathrm{SB} T_\mathrm{eff}^4\).

The energy flux arriving at Earth, i.e., the Solar constant \(S\) is then given by

\(\Large S = \frac{L}{4\pi d^2}\).

Only one hemisphere of the Earth is illuminated at any time, so the energy that the Earth is receiving per unit time, neglecting any albedo effects, is given by

\(\Large P_\mathrm{in} = \pi R_\oplus^2 S\).

We assume that the temperature on Earth is in euilibrium, i.e., there is no difference between day and night. The energy the Earth is emitting per unit time is then given by

\(\Large P_\mathrm{out} = 4\pi R_\oplus^2 \sigma_\mathrm{SB} T^4\).

In equilibrium both are equal \(\left( P_\mathrm{in} = P_\mathrm{out} \right)\) and we can solve for \(T\)

\(\Large T = T_\mathrm{eff} \sqrt[4]{\frac{R_\odot^2}{4d^2}}\).

So the equilibrium temperature on Earth depends on the Sun’s effective temperature \(T_\mathrm{eff}\), the Sun’s radius \(R_\odot\), and the Earth-Sun distance \(d\). We can now write a function that takes the frame object as argument and returns Earth’s temperature.

[40]:
def T(frame):
    """Function computes the equilibrium temperature of the Earth"""
    return frame.Sun.T * (frame.Sun.R**2 / (4.*frame.Earth.d**2))**0.25

We can now assign this function to the updater of the temperature field

[41]:
sim.Earth.T.updater = T

The function for calculating the Earth-Sun distance is pretty simple and can be simply assigned to the updater of the field.

[42]:
def d(sim):
    """Function computed the Sun-Earth distance"""
    return np.linalg.norm(sim.Earth.r - sim.Sun.r)
[43]:
sim.Earth.d.updater = d

Right now the temperature is zero as we initialized it.

[44]:
sim.Earth.T
[44]:
0.0

We can now use our the updater of the temperature field to calculate it’s initial value

[45]:
sim.Earth.T.update()
[46]:
sim.Earth.T
[46]:
278.6190681198289

Updating groups

But before we can rerun the simulation we have to tell simframe how to update groups. This is not done automatically, because for one, not all fields and groups need to be updated, and second, the order of update matters. In our case, we do not need to update the group of the Sun, because it’s temperature is not changing, and we need to update the Earth-Sun distance before we update the Earth’s temperature, because we need the distance for it.

The only update operation that is performed once per time step is the update of the frame object Frame.update(). From here we have to delegate tasks down the tree structure.

There are in principle two methods of updating groups. One is by writing a function and assigning it to the group’s updater. Here we write a function that is calling Earth’s updater and assign it to theupdater of the frame.

[47]:
def update_Earth(frame):
    frame.Earth.update()
[48]:
sim.updater = update_Earth

The second method is by assigning a list of group attributes to an updater. The updater is then calling the updaters of these attributes in exactly that order.

[49]:
sim.Earth.updater = ["d", "T"]    # Earth-Sun distance first, then temperature

If the updater of a group has been set with a list, the order can be displayed. The attribute will return None, if no updater is set or if the updater was set with a function.

[50]:
sim.Earth.updateorder
[50]:
['d', 'T']

Resetting to initial conditions

We can now reset to the initial conditions we have saved earlier, reset the namespace writer, and run the simulation again.

[51]:
sim.Earth.r = r_Earth_ini
sim.Earth.v = v_Earth_ini
sim.Sun.r   = r_Sun_ini
sim.Sun.v   = v_Sun_ini
sim.t       = 0.
sim.writer.reset()
[52]:
sim.run()
Execution time: 0:00:01

Reading and plotting

[53]:
data_T = sim.writer.read.all()
[54]:
def to_celsius(K):
    return K - 273.15
[55]:
def plot_temperature(data):
    fig, ax = plt.subplots(dpi=150)
    ax.plot(data.t/year, to_celsius(data.Earth.T), c="C3", label="Temperature")
    ax.plot(-2., 0., label="Earth-Sun distance", c="C0")
    ax.set_xlabel("Time [years]")
    ax.set_ylabel("Temperature [°C]")
    ax.set_xlim(data.t[0]/year, data.t[-1]/year)
    ax.set_ylim(-100., 150.)
    axr = ax.twinx()
    axr.plot(data.t/year, data.Earth.d/AU, c="C0")
    axr.set_ylabel("Distance [AU]")
    axr.set_ylim(0., 2.0)
    ax.legend()
    fig.tight_layout()
    plt.show()
[56]:
plot_temperature(data_T)
_images/3_updating_84_0.png

The heartbeat concept

Every updater is actually performing three operations. The updater operation itself, a “systole”, that is executed before the update operation and a “diastole” that is executed after the update operation.

In this simple case we want to print Earth’s temperature before and after it’s calculation.

[57]:
def T_sys(frame):
    """This function prints the temperature before the update process."""
    msg = "{:6s}: T = {:5.1f} K".format("Before", sim.Earth.T)
    print(msg)
[58]:
def T_dia(frame):
    """This function prints the temperature after the update process."""
    msg = "{:6s}: T = {:5.1f} K".format("After", sim.Earth.T)
    print(msg)

We can assign these functions to systole and diastole of the temperature updater.

[59]:
sim.Earth.T.updater.systole = T_sys
sim.Earth.T.updater.diastole = T_dia
[60]:
sim.Earth.T = 0.
[61]:
sim.Earth.T.update()
Before: T =   0.0 K
After : T = 278.0 K

Note: The updater of an integration variable is not directly setting its new value. Only the integrator is setting the new value of the integration variable after all integration instructions have been executed successfully. This has the consequence, that the diastole of an integration variable DOES NOT have access to the updated value.

Retrieving updater information

It is possible to retrieve information about the currently assigned updater functions.

[62]:
sim.Earth.T.updater
[62]:
Heartbeat
---------

Systole:  <function T_sys at 0x7f6f8bb1d9e0>
Updater:  <function T at 0x7f6f8a5aaf20>
Diastole: <function T_dia at 0x7f6f8bb1da80>

Docstrings
----------

Systole:
This function prints the temperature before the update process.

Updater:
Function computes the equilibrium temperature of the Earth

Diastole:
This function prints the temperature after the update process.

The source code of the assigned functions can also be displayed by adressing the respective sub-level. For example the systole or the actual updater. This may not work if the assigned function is not a Python function.

[63]:
sim.Earth.T.updater.systole
[63]:
Updater
-------

Signature: T_sys(frame)
Source:
def T_sys(frame):
    """This function prints the temperature before the update process."""
    msg = "{:6s}: T = {:5.1f} K".format("Before", sim.Earth.T)
    print(msg)

File: /tmp/ipykernel_23644/2382827243.py
Type: function
[64]:
sim.Earth.T.updater.updater
[64]:
Updater
-------

Signature: T(frame)
Source:
def T(frame):
    """Function computes the equilibrium temperature of the Earth"""
    return frame.Sun.T * (frame.Sun.R**2 / (4.*frame.Earth.d**2))**0.25

File: /tmp/ipykernel_23644/4057753429.py
Type: function

If the updater of a group was defined by assigning a list as in this example, this will be shown as well.

[65]:
sim.Earth.updater
[65]:
Heartbeat
---------

Systole:  None
Updater:  list_updater
Diastole: None

Docstrings
----------

Systole:
None

Updater:
The attributes in this group are updated in the order:
['d', 'T'].

Diastole:
None