Start by loading some boiler plate: matplotlib, numpy, scipy, json, functools, and a convenience class.
In [1]:
%matplotlib inline
import matplotlib
matplotlib.rcParams['figure.figsize'] = (10.0, 8.0)
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline
from scipy.interpolate import UnivariateSpline
import json
import pandas as pd
from functools import partial
class Foo: pass
And some more specialized dependencies:
Slict provides a convenient slice-able dictionary interfaceChest is an out-of-core dictionary that we'll hook directly to a globus remote using...glopen is an open-like context manager for remote globus files
In [2]:
from chest import Chest
from slict import CachedSlict
from glopen import glopen, glopen_many
Helper routines
In [3]:
def load_from_archive(names, arch):
cs = []
for name in names:
cs.append(Chest(path = "{:s}-results".format(name),
open = partial(glopen, endpoint=arch),
open_many = partial(glopen_many, endpoint=arch)))
scs = [CachedSlict(c) for c in cs]
ps = []
for name in names:
with glopen(
"{:s}.json".format(name), mode='r',
endpoint = arch,
) as f:
ps.append(json.load(f))
return cs, scs, ps
Configuration for this figure.
In [4]:
config = Foo()
config.names = [
# "Wilk/Wilk_kmin_2.5/Wilk_kmin_2.5",
# "Wilk/Wilk_kmin_3.5/Wilk_kmin_3.5",
# "Wilk/Wilk_kmin_4.5/Wilk_kmin_4.5",
"Wilk/Wilk_long/Wilk_long",
]
#config.arch_end = "maxhutch#alpha-admin/~/pub/"
#config.arch_end = "alcf#dtn_mira/projects/alpha-nek/experiments/"
config.arch_end = "alcf#dtn_mira/projects/PetaCESAR/maxhutch/"
height = 'H_exp'
Open a chest located on a remote globus endpoint and load a remote json configuration file.
In [5]:
cs, scs, ps = load_from_archive(config.names, config.arch_end);
We want to plot the spike depth, which is the 'H' field in the chest. Chests can prefetch lists of keys more quickly than individual ones, so we'll prefetch the keys we want.
In [6]:
for c,sc in zip(cs, scs):
c.prefetch(sc[:,height].full_keys())
Use a spline to compute the derivative of 'H' vs time: the Froude number.
In [7]:
spls = []
for sc, p in zip(scs, ps):
T = np.array(sc[:,height].keys())
H = np.array(sc[:,height].values()) #- 2 * np.sqrt(p['conductivity']* (T + p['delta']**2 / p['conductivity'] / 4))
spls.append(UnivariateSpline(T,
H,
k = 5,
s = 1.e-12))
Frs = [spl.derivative() for spl in spls]
Tss = [np.linspace(sc[:,height].keys()[0], sc[:,height].keys()[-1], 1000) for sc in scs]
In [8]:
Run37 = pd.DataFrame.from_csv('WRun37 4.49.56 PM 7_3_07.txt', sep='\t')
Run58 = pd.DataFrame.from_csv('WRun058 4.32.52 PM 7_3_07.txt', sep='\t')
Run78 = pd.DataFrame.from_csv('WRun078 4.49.56 PM 7_3_07.txt', sep='\t')
def plot_exp(data, n, fmt):
norm = .5*( np.sqrt(data["Atwood"]/(1-data["Atwood"])*data["Accel. [mm/sec^2]"]* 76 / n)
+ np.sqrt(data["Atwood"]/(1+data["Atwood"])*data["Accel. [mm/sec^2]"]* 76 / n))
axs.plot(
data["AvgAmp (mm)"] * n / 76,
data["Average Velocity"]/norm, fmt);
#data["Froude Average"], fmt);
return
Plot the Froude number, non-dimensionalized by the theoretical dependence on Atwood, acceleration, and wave number, vs the spike depth, normalized by wave-length.
The dotted line is the theoretical prediction of Goncharaov. The solid black line is the farthest that Wilkinson and Jacobs were able to get.
In [9]:
fig, axs = plt.subplots(1,1)
for p, spl, Fr, T in zip(ps, spls, Frs, Tss):
axs.plot(
spl(T) * p["kmin"],
Fr(T)/ np.sqrt(p["atwood"]*p["g"] / p["kmin"]),
label="{:3.1f} modes".format(p["kmin"]));
#axs.plot(Run37["AvgAmp (mm)"] * 2.5 / 76, Run37["Froude Average"], "bx");
#plot_exp(Run37, 2.5, "bx")
#plot_exp(Run78, 3.5, "gx")
plot_exp(Run58, 4.5, "bx")
axs.plot([0,10], [np.sqrt(1/np.pi), np.sqrt(1/np.pi)], 'k--')
axs.axvline(x=1.4, color='k');
axs.set_ylabel(r'Fr')
axs.set_xlabel(r'$h/\lambda$');
axs.legend(loc=4);
axs.set_xbound(0,3);
axs.set_ybound(0,1.5);
plt.savefig('Figure17_long.png')
In [10]:
fig, axs = plt.subplots(1,1)
for sc, p, spl, Fr, T in zip(scs, ps, spls, Frs, Tss):
axs.plot(
T,
spl(T) * p["kmin"],
label="{:3.1f} modes".format(p["kmin"]));
axs.plot(
sc[:,height].keys(),
np.array(sc[:,height].values())*p['kmin'],
'bo');
#axs.plot(Run37["Time (sec)"]-.5, Run37["AvgAmp (mm)"] * 2.5 / 76, "bx");
axs.plot(Run58["Time (sec)"]-.515, Run58["AvgAmp (mm)"] * 4.5 / 78, "bx");
#axs.plot(Run78["Time (sec)"]-.5, Run78["AvgAmp (mm)"] * 3.5 / 76, "gx");
axs.set_ylabel(r'$h/\lambda$')
axs.set_xlabel(r'T (s)');
axs.set_xbound(0.0,1.5);
axs.set_ybound(-0.0,4);
axs.legend(loc=4);
In [11]:
%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
%load_ext version_information
%version_information numpy, matplotlib, slict, chest, glopen, globussh
Out[11]: