Repository navigation
AR 13664
On 2024 May 10/11, the strongest geomagnetic storm since November 2003 has occurred, with a peak Dst index of -412 nT. The storm was caused by NOAA Active Region (AR) 13664, which was the source of a large number of fast coronal mass ejections (CMEs) and big flares, including 12 X-class flares.
We perform 3D magnetic field extrapolations with the NF2 nonlinear-force free code based on physics informed neural networks. This tutorial provides an overview on how to access the derived data products and how to load snapshots for further analysis.
Overview of the evolution of NOAA AR 13664 from 2024 May 5-00:00 to 11-04:36 UT. a) Evolution of the vertical
distribution of free magnetic energy. b)
GOES 1-8 Å soft X-ray flux. c) Observed SDO/HMI radial magnetic field maps, SDO/AIA 131 Å EUV maps, and the
modeled current density integrated along the vertical axis. We show the initial flux emergence (snapshot at 2024 May 7
00:00UT) and the major flare eruptions (indicated by dashed lines). The magenta line indicates the point where the
right edge of the SHARP reaches 60° longitude, beyond which the reliability of the vector magnetograms and the
modeled coronal field decreases. The circles in the integrated current density maps indicate points of flux domain
interactions. An animation of the full time series is available in the online journal. Left top: B_z of the modeled
photospheric magnetic field. Right top: EUV observations from SDO/AIA 131 Å EUV. Left bottom: maps of integrated
current density. Right bottom: maps of integrated free magnetic energy. The sequence shows the re-configuration of
magnetic domains (dark lines in current density maps) during strong solar eruptions.
All data products are publicly available through Globus website.
-
/nf2_series/: modeling result of NF2 snapshots. Contains all information for further evaluations and is compact (2 MB per file). -
/evaluation/pkl/: contains pickle files with evaluations -
/hdf5/: extrapolated magnetic field (B) of snapshots -
/vtk/: selected VTK snapshots for visualization with e.g. Paraview -
/sharps/: input SHARP vector magnetograms
The simulated volume is parametrized as a neural representation, where the weights of the neural network store all the information. This is a very efficient representation which makes the save states efficient for data transfer, however additional computing effort is required to transfer this representation to a grid for further evaluation.
Download sample nf2 file
import urllib.request
nf2_sample_path = 'sample.nf2'
urllib.request.urlretrieve("https://g-66affc.7a577b.6fbd.data.globus.org/nf2_series/20240507_172400_TAI.nf2", nf2_sample_path)Load NF2 file
from nf2.evaluation.output import CartesianOutput
nf2_model = CartesianOutput(nf2_sample_path)List all available metrics that can be loaded
from nf2.evaluation.output_metrics import metric_mapping
print(f"{'Metrics':<20}")
print(f"{'-'*20}")
for metric in metric_mapping:
print(f"{metric:<20}")- j
- b_nabla_bz
- alpha
- spherical_energy_gradient
- energy_gradient
- magnetic_helicity
- los_trv_azi
- free_energy
- squashing_factor
Load the extrapolation result from 0 to 20 Mm
nf2_out = nf2_model.load_cube(height_range=[0, 20], Mm_per_pixel=0.72, metrics=['j'], progress=True)Print overview of outputs
import numpy as np
# print all keys and subkeys
print(f"{'Keys':<20}{'Values':<20}")
print(f"{'-'*20}{'-'*20}")
for key in nf2_out.keys():
if isinstance(nf2_out[key], dict):
print(f"{key:<20}{list(nf2_out[key].keys())}")
elif isinstance(nf2_out[key], np.ndarray):
print(f"{key:<20}{nf2_out[key].shape}")
elif isinstance(nf2_out[key], float):
print(f"{key:<20}{nf2_out[key]:.2f}")
else:
print(f"{key:<20}")Keys Values
----------------------------------------
b (601, 257, 29, 3)
a (601, 257, 29, 3)
jac_matrix (601, 257, 29, 3, 3)
metrics ['j', 'alpha']
coords (601, 257, 29, 3)
Mm_per_pixel 0.72
Plot integrated current density
from matplotlib import pyplot as plt
from astropy import units as u
cm_per_pixel = nf2_out['Mm_per_pixel'] * 1e-8
current_density = np.linalg.norm(nf2_out['metrics']['j'], axis=-1)
integrated_current_density = current_density.sum(axis=2).to_value(u.G / u.s) * cm_per_pixel
extent = [0, integrated_current_density.shape[0] * nf2_out['Mm_per_pixel'], 0,
integrated_current_density.shape[1] * nf2_out['Mm_per_pixel']]
fig, ax = plt.subplots(1, 1, figsize=(5, 2), dpi=200)
im = ax.imshow(integrated_current_density.T, cmap='inferno', origin='lower', extent=extent, norm='log')
plt.colorbar(im, ax=ax, label=r'$\Vert \mathbf{J} \Vert$ [G cm s$^{-1}$]')
ax.set_xlabel("x [Mm]")
ax.set_ylabel("y [Mm]")
plt.tight_layout()
plt.show()
We do not provide all time steps as VTK and the files are typically inefficient for download. NF2 files can be directly converted from the command line.
nf2-to-vtk --nf2_path "<<path to .nf2 file>>" --out_path "<<path to save .vtk file>>"With an available GPU conversions are fast. When using CPU resources the conversion can be optimized by reducing the resolution or cropping the height range (most current carrying fields are <20Mm).
- Use
--Mm_per_pixel 0.72for bin2 of the SHARP resolution - Use
--height_range 0 20to only load the first 20 Mm.
Additional metrics to be loaded can be specified by --metrics. For example for the current density and free magnetic energy set --metrics j free_energy (see list above for more metrics).
For the computation of the squashing factor and twist an additional installation of FastQSL is required. Note that this installation requires an available GPU.
GPU-version:
pip install git+https://github.com/peijin94/FastQSL
pip install cupy
The 3D squashing factor and twist map can be loaded as additional metrics squashing_factor (see load VTK or load NF2 files).
HDF5 snapshots are the easiest way to obtain the extrapolated magnetic fields. They provide the magnetic field solution and the World-Coordinate-System (WCS) header information. They further provide the scaling in Mm per pixel and the timestamp. Further quantities like free magnetic energy, and current density can be loaded from the NF2 snapshots.
Download sample HDF5 file from Globus
import urllib.request
download_path = 'sample.hdf5'
urllib.request.urlretrieve("https://g-66affc.7a577b.6fbd.data.globus.org/hdf5/20240508_043600_TAI.hdf5", download_path)Read HDF5 file
# read hdf5 file
import h5py
import numpy as np
with h5py.File(download_path, "r") as f:
# primary keys --> B
print(f"Keys: {f.keys()}")
# meta information
print(f"Attributes:")
for k, v in f.attrs.items():
print(f"{k}: {v}")
# data
nf2_b = np.array(f['B'])
Mm_per_pixel = f.attrs['Mm_per_pixel']Plot 5 slices of the vertical magnetic field
from matplotlib import pyplot as plt
# select 5 slices across the height of the frame
height_idx = np.linspace(0, nf2_b.shape[2] - 1, 5, dtype=int)
nf2_heights = height_idx * Mm_per_pixel
# setup scale of frame
extent = [0, nf2_b.shape[0] * Mm_per_pixel, 0, nf2_b.shape[1] * Mm_per_pixel]
# plot slices
fig, axs = plt.subplots(1, 5, figsize=(30, 3), dpi=200)
for i, ax in enumerate(axs):
idx = height_idx[i]
h = nf2_heights[i]
b_slice = nf2_b[:, :, idx, 2]
# plot B_z
v_min_max = np.max(np.abs(b_slice))
im = ax.imshow(nf2_b[:, :, idx, 2].T, cmap='gray', vmin=-v_min_max, vmax=v_min_max, origin='lower', extent=extent)
# add colorbar
plt.colorbar(im, ax=ax, label='$B_z$ [G]')
ax.set_xlabel("x [Mm]")
ax.set_ylabel("y [Mm]")
ax.set_title(f"Height {h:.1f} Mm")
plt.tight_layout()
plt.show()
The derived quantities of the full time series are provided as pickle files (/evaluation/pkl/).
Download evaluation files
from tqdm import tqdm
import urllib.request
import os
import pandas as pd
pkl_path = 'evaluation_results'
t_start = '20240509T120000' # minimum date is: 20240505T000000
t_end = '20240510T000000' # maximum date is: 20240511T043600
cadence = '12min' # max cadence is 12min --> HMI cadence
os.makedirs(pkl_path, exist_ok=True)
for date in tqdm(pd.date_range(t_start, t_end, freq=cadence), desc='Downloading pkl files'):
date_str = date.strftime('%Y%m%d_%H%M%S_TAI')
pkl_url = f'https://g-66affc.7a577b.6fbd.data.globus.org/evaluation/pkl/{date_str}.pkl'
save_path = os.path.join(pkl_path, f'{date_str}.pkl')
if os.path.exists(save_path):
print(f"File {save_path} already exists, skipping")
continue
try:
urllib.request.urlretrieve(pkl_url, save_path)
except:
print(f"Failed to download {pkl_url}")Load result files
from nf2.evaluation.sharp.convert_series import load_results
import numpy as np
output = load_results(os.path.join(pkl_path, '*.pkl'))Print overview of result dictonary
def print_keys(output):
for key in output.keys():
if isinstance(output[key], dict):
print(f"{key:<20}{list(output[key].keys())}")
elif isinstance(output[key], np.ndarray):
print(f"{key:<20}{output[key].shape}")
elif isinstance(output[key], float):
print(f"{key:<20}{output[key]:.2f}")
else:
print(f"{key:<20}")
print(f"{'Keys':<20}{'Values':<20}")
print(f"{'-'*20}{'-'*20}")
print_keys(output)
print('\n\n')
print('Integrated quantities:')
print(f"{'Keys':<20}{'Values':<20}")
print(f"{'-'*20}{'-'*20}")
print_keys(output['integrated_quantities'])
print('\n\n')
print('maps:')
print(f"{'Keys':<20}{'Values':<20}")
print(f"{'-'*20}{'-'*20}")
print_keys(output['maps'])Keys Values
----------------------------------------
times (61,)
integrated_quantities['energy', 'free_energy']
metrics ['divergence', 'jxb', 'theta']
maps ['b_0', 'current_density_map', 'energy_map', 'free_energy_map', 'jxb_map']
height_distribution ['height_free_energy', 'height_current_density']
wcs (61, 1)
Mm_per_pixel 0.72
Integrated quantities:
Keys Values
----------------------------------------
energy (61,)
free_energy (61,)
maps:
Keys Values
----------------------------------------
b_0 (61, 601, 257)
current_density_map (61, 601, 257)
energy_map (61, 601, 257)
free_energy_map (61, 601, 257)
jxb_map (61, 601, 257)
Load current density maps as SunPy map
from sunpy.map import Map
idx = 0
wcs = output['wcs'][idx, 0]
data = output['maps']['current_density_map'][idx]
map = Map(data.T, wcs)
map.plot(cmap='inferno', norm='log')
Plot time series of free magnetic energy
import matplotlib.pyplot as plt
from astropy import units as u
fig, ax = plt.subplots(figsize=(10, 5))
ax.xaxis_date()
ax.plot(output['times'], output['integrated_quantities']['free_energy'].to_value(u.erg) * 1e-32)
ax.set_xlabel('Time')
ax.set_ylabel('Free Magnetic Energy [$10^{32}$ erg]')
plt.show()
publication in press (The Astrophysical Journal Letters)
@article{jarolim2024ar13664,
title={Magnetic Field Evolution of the Solar Active Region 13664},
author={Jarolim, Robert and Veronig, Astrid and Purkhart, Stefan and Zhang, Peijin and Rempel, Matthias},
journal={The Astrophysical Journal Letters},
year={2024},
publisher={IOP Publishing}
year={2024}
}