<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Stokes solver for asymptotic flow
Inport away.
Step1:  Load a frame from a real simulation.
Step2:  Load the governing properties from the frame.
Step3:  Load the last midplane slice of the scalar field and manipulate it into a periodic box.
Step4:  Make sure it looks OK.
Step5:  Stokes is linear and we have periodic boundaries, so we can solve it directly using Fourier transforms and the frequency-space Green's function (which is diagonal).
Step6:  Look ok?
Step7:  Now we want to turn this into something Darcy-Weisbach-esque.  We don't have uniform forcing, so we take an average.
Step8:  Rayleight-Taylor types like the Froude number, which doesn't really make sense
Step9:  Instead, we normalize "the right way" using the viscosity
Step10:  Print everything out.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib
matplotlib.rcParams['figure.figsize'] = (10.0, 16.0)
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
from scipy import fftpack
from numpy import fft
import json
from functools import partial
class Foo: pass
from chest import Chest
from slict import CachedSlict
from glopen import glopen, glopen_many
Explanation: Stokes solver for asymptotic flow
Inport away.
End of explanation
name = "HighAspect/HA_viscosity_4.0E-4/HA_viscosity_4.0E-4"
arch = "alcf#dtn_mira/projects/alpha-nek/experiments"
c = Chest(path="{:s}-results".format(name), 
                       open=partial(glopen, endpoint=arch),
                       open_many = partial(glopen_many, endpoint=arch))
sc = CachedSlict(c)
c.prefetch(sc[:,'t_xy'].full_keys())
with glopen(
            "{:s}.json".format(name), mode='r',
            endpoint = arch,
            ) as f:
    p = json.load(f)
Explanation: Load a frame from a real simulation.
End of explanation
L = 1./p["kmin"]
Atwood = p["atwood"]
g = p["g"]
viscosity = p["viscosity"]
Explanation: Load the governing properties from the frame.
End of explanation
T_end = sc[:,'t_xy'].keys()[-1]
phi_raw = sc[T_end, 't_xy']
phi_raw = np.concatenate((phi_raw, np.flipud(phi_raw)), axis=0)
phi_raw = np.concatenate((phi_raw, np.flipud(phi_raw)), axis=0)
phi_raw = np.concatenate((phi_raw, np.fliplr(phi_raw)), axis=1)
phi_raw = np.concatenate((phi_raw, np.fliplr(phi_raw)), axis=1)
raw_shape = phi_raw.shape
nx = raw_shape[0]
ny = raw_shape[0]
phi = phi_raw[nx/8:5*nx/8, ny/8:5*ny/8]
nx = phi.shape[0]
ny = phi.shape[1]
Explanation: Load the last midplane slice of the scalar field and manipulate it into a periodic box.
End of explanation
plt.figure()
plt.imshow(phi)
plt.colorbar();
Explanation: Make sure it looks OK.
End of explanation
# Setup the frequencies
dx = L / ny
X = np.tile(np.linspace(0, L, nx), (ny, 1))
Y = np.tile(np.linspace(0, L, ny), (nx, 1)).transpose()
rfreqs = fft.rfftfreq(nx, dx) * 2 * np.pi;
cfreqs = fft.fftfreq(nx, dx)* 2 * np.pi;
rones = np.ones(rfreqs.shape[0]);
cones = np.ones(cfreqs.shape[0]);
# RHS comes from the forcing
F = phi * Atwood * g / viscosity
# Transform forward
p1 = fft.rfftn(F)
# Green's function
p1 = p1 / (np.square(np.outer(cfreqs, rones)) + np.square(np.outer(cones, rfreqs)))
p1[0,0] = 0
# Transform back
w = fft.irfftn(p1)
Explanation: Stokes is linear and we have periodic boundaries, so we can solve it directly using Fourier transforms and the frequency-space Green's function (which is diagonal).
End of explanation
plt.figure()
plt.imshow(w)
plt.colorbar();
Explanation: Look ok?
End of explanation
A_tilde = np.sum(np.abs(phi))/ (nx * ny)
Explanation: Now we want to turn this into something Darcy-Weisbach-esque.  We don't have uniform forcing, so we take an average.
End of explanation
Froude = np.sum(np.abs(w)) / np.sqrt(g * Atwood * A_tilde * L) / (nx * ny)
Explanation: Rayleight-Taylor types like the Froude number, which doesn't really make sense:
$$ \text{Fr} = \frac{u}{\sqrt{A g L}} $$
End of explanation
Right = np.sum(np.abs(w)) * viscosity / (g * Atwood * A_tilde * L**2)/ (nx*ny)
Explanation: Instead, we normalize "the right way" using the viscosity:
End of explanation
dff = 64 / 16 * 14.227
print("L={:f}, A={:f}, A_til={:f}, g={:f}, nu={:f}.  D-W is {:f}".format(
       L, Atwood, A_tilde, g, viscosity, 1./(dff*2)))
print("  Froude: {:10f} |  Right: {:10f}".format(Froude, Right))
print("  C1 = {:f} * C0 ".format(1./Right))
Explanation: Print everything out.
End of explanation
<END_TASK>