Tutorial: Interpolating simulation fields in time
[1]:
try:
from google.colab import drive
%pip install -Uq git+https://github.com/seap-udea/fargopy
except ImportError:
print("Not running in Colab, skipping installation")
%load_ext autoreload
%autoreload 2
!mkdir -p ./gallery/
Not running in Colab, skipping installation
In this notebook, we illustrate how to interpolate between snapshots, ie. in the time direction, data from FARGO3D
With the interpolated data, we can generate enhanced plots and animations that provide deeper insights into the simulation results. This allows for better visualization and understanding of the physical phenomena represented in the data.
Before starting
If you are in Google Colab, install the latest version of the package:
For this tutorial you will need the following external modules and tools:
[2]:
import fargopy as fp
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from celluloid import Camera
from tqdm import tqdm
from IPython.display import HTML
Running FARGOpy version 1.1.1.
NOTE: Since alpha versions (<=0.X.X) a major refactor has been done in versions 1.1.X.
Please check the documentation for more information.
Let’s FARGOpy
Download, for instance, the 3D simulation of a disk with a Jovian planet:
[3]:
PATH = fp.Simulation.download_precomputed('p3disoj')
sim = fp.Simulation(output_dir=PATH)
Precomputed output directory '/tmp/p3disoj' already exist
Your simulation is now connected with '/Users/jzuluaga/fargo3d/'
Now you are connected with output directory '/tmp/p3disoj'
Found a variables.par file in '/tmp/p3disoj', loading properties
Loading variables
85 variables loaded
Simulation in 3 dimensions
Loading domain in spherical coordinates:
Variable phi: 128 [[0, np.float64(-3.117048960983623)], [-1, np.float64(3.117048960983623)]]
Variable r: 64 [[0, np.float64(0.5078125)], [-1, np.float64(1.4921875)]]
Variable theta: 32 [[0, np.float64(1.4231400767948967)], [-1, np.float64(1.5684525767948965)]]
Number of snapshots in output directory: 11
Planets found in summary.dat
Name: Jupiter, Initial pos: [1.0, 0.001, 0.0], Mass: 0.001
Now we can load multiples fields as a list from specific snapshots:
[4]:
fields_interp = sim.load_field(
fields='gasdens',
slice='phi=0',
snapshot=(0,10)
)
Xrot = fields_interp.var1_mesh[0] # X mesh
Yrot = fields_interp.var2_mesh[0] # Y mesh
Zrot = fields_interp.var3_mesh[0] # Z mesh
Time Interpolation
1D + Time
Lets see the time interpolatión in action:
[14]:
fields = sim.load_field(
fields='gasdens',
slice="phi=0,theta= 89 deg",
snapshot=[0,10])
xsim = fields.var1_mesh[0]
ysim = fields.var2_mesh[0]
zsim = fields.var3_mesh[0]
rsim = np.sqrt(xsim**2 + ysim**2 + zsim**2)
times = np.linspace(0, 1.0, 100) # 100 timesteps
[15]:
snap_table = fields.snapshot_time_table
snapshots = snap_table["Snapshot"].values
norm_times = snap_table["Normalized_time"].values
times = np.linspace(0, 1.0, 100) # 100 timesteps
tolerance = (norm_times[1] - norm_times[0])
plt.ioff()
fig, axs = plt.subplots(1, 1, figsize=(8, 6))
camera = Camera(fig)
for t in times:
fp.Plot.fargopy_mark(axs)
rmin = sim.domains.r.min()
rmax = sim.domains.r.max()
rhosim = fields.gasdens_mesh[0]
rs=np.linspace(rmin, rmax, 100)
dens = fields.evaluate(field='gasdens',time=t, var1=rs)
axs.plot(rs, dens, color='orangered', label='Space-time interpolated data')
idx = np.argmin(np.abs(norm_times - t))
if np.abs(norm_times[idx] - t) < tolerance:
snap = int(snapshots[idx])
fields = sim.load_field(
fields='gasdens',
slice="phi=0,theta=1.56",
snapshot=snap
)
r = np.sqrt(fields.var1_mesh[0]**2+fields.var2_mesh[0]**2+fields.var3_mesh[0]**2)
rhosim = fields.gasdens_mesh[0]
axs.scatter(rsim, rhosim, color='dodgerblue', label=f'Simulated data')
axs.set_xlabel('$x$ [au]')
axs.set_ylabel(r'$\rho$ [g/cm]')
axs.text(0.03, 0.98, f'Normalized time: {t:.2f}',
transform=axs.transAxes, ha='left', va='top',
fontsize=12, bbox=dict(facecolor='white', alpha=0.7, edgecolor='none'))
handles, labels = axs.get_legend_handles_labels()
by_label = dict(zip(labels, handles))
axs.legend(by_label.values(), by_label.keys())
axs.grid(0.3)
camera.snap()
plt.ion()
animation = camera.animate(interval=150, blit=False)
animation.save('gallery/fargopy-tutorial-interpolation_4.gif', writer='pillow') # Save animation
print(f"Generating animation (it may take a while, go for a coffee)...")
HTML(animation.to_jshtml())
Generating animation (it may take a while, go for a coffee)...
[15]:
2D + Time
Now we can obtain the density at any point in space and at any moment in time, allowing us to observe its temporal evolution in greater detail.
[5]:
steps=100 # Number of time steps for the interpolation
times= np.linspace(0, 1.0, steps)
Let’s create an animation to visualize this behavior dynamically.
[6]:
plt.ioff()
fig,axs = plt.subplots(1,1,figsize=(6,6))
cmap = 'Spectral_r'
camera = Camera(fig)
xmin,xmax = fields_interp.var1_mesh[0].min(), fields_interp.var1_mesh[0].max()
zmin,zmax = fields_interp.var3_mesh[0].min(), fields_interp.var3_mesh[0].max()
axs.set_xlim(xmin, xmax)
axs.set_ylim(zmin, zmax)
axs.set_xlabel('$x$ [au]')
axs.set_ylabel('$z$ [au]')
fp.Plot.fargopy_mark(axs)
times= np.linspace(0, 1.0, 100)
for t in tqdm(times):
dens = fields_interp.evaluate(
time=t,
var1=Xrot,
var3=Zrot
)
dens = np.log10(dens * sim.URHO)
pc = axs.contourf(Xrot, Zrot, dens, cmap=cmap, levels=100)
axs.contour(Xrot, Zrot, dens, levels=20, colors='k', linewidths=0.5)
axs.text(0.03, 0.98, f'Normalized time: {t:.2f}',transform=axs.transAxes,
ha='left', va='top', fontsize=12, bbox=dict(facecolor='white', alpha=0.7, edgecolor='none'))
camera.snap()
plt.colorbar(pc, ax=axs, label=f'Gas Density Log($g/cm^3$)')
plt.ion()
animation = camera.animate(interval=100,blit=True)
animation.save('gallery/fargopy-tutorial-interpolation.gif', writer='pillow')
print(f"Generating animation (it may takes a while, go for a coffee)...")
HTML(animation.to_jshtml())
100%|██████████| 100/100 [00:02<00:00, 42.92it/s]
Generating animation (it may takes a while, go for a coffee)...
[6]: