SPINS Science Files

SPINS comes with a variety of Science files which help with simulation data analysis. Below is a list of currently available functions with a brief description of their purpose and the assumptions behind their implementation. These Science files are in spins/src/Science.

Unless otherwise mentioned, the functions below work for both unmapped and mapped domains where the mapping is along the bottom of the domain (where the bottom is at the 0-index in the z-dimension, as opposed to at the $N_{z}-1$ index). Another way to say this is that z increases with the index.

Details and Notes

Notation

Warning! Multiple different notations are used below. Subscripts of $u$, $v$, and $w$ are always derivatives with respect to the subscript variable. However, subscripts of $\omega$ and $t$ are components. That is, $u_{x}$ and $\omega _{x}$ are the x-derivative of $u$ and the x-component of vorticity.

Units

The units of each of the functions have been built under the assumption that physical units are being used in the simulation. These computations will work for nondimensional simulations, but more care must be taken to ensure the scalings are appropriate. Physical units are typically:

Length	m
Time	s
Density	kg/m^3

The reference density should be 1000 kg/m^3 for the stress and dissipation calculations to be physically correct. The density field can be (and often is) non-dimensionalized to $\rho '=\rho /\rho _{0}-1$. If nondimensionalized, compute_Background_PE.cpp and compute_BPE_from_internal.cpp assume this form of scaling. These two functions also need to know if the density field is dimensional or not.

Boundary Stresses

These functions compute the surface stresses on the top and bottom boundaries. It is assumed that $z=z_{\text{max}}$ is the top and $z=z_{\text{min}}$ is the bottom.

Further calculations and saving to disk of various quantities of the stresses are made in the src/BaseCase directory.

bottom_slope.cpp

Compute the vector of the bottom slope at each x-grid point. This is used for calculating the stresses along the bottom boundary (cf. bottom_stress_x.cpp and bottom_stress_y.cpp).

bottom_stress_x.cpp

Compute the 2D vector of the x-component of stress on the bottom boundary. The z boundary condition must be no slip.

If the case is unmapped:

$t_{x}=\mu u_{z}$

where $\mu =\nu \rho _{0}$ is the dynamic viscosity.

If the case is mapped:

$t_{x}={\frac {\mu }{\alpha ^{2}(x)}}\left(2h'(x)(w_{z}-u_{x})+(u_{z}+w_{x})(1-h'^{2}(x))\right)$

where $\alpha (x)={\sqrt {1+h'^{2}(x)}}$ and where $h'(x)$ is the bottom slope.

bottom_stress_y.cpp

Compute the 2D vector of the y-component of stress on the bottom boundary. The z boundary condition must be no slip.

If the case is unmapped:

$t_{y}=\mu v_{z}$

where $\mu =\nu \rho _{0}$ is the dynamic viscosity.

If the case is mapped:

$t_{y}={\frac {\mu }{\alpha (x)}}\left(v_{z}-h'(x)v_{x}\right)$

where $\alpha (x)={\sqrt {1+h'^{2}(x)}}$ and where $h'(x)$ is the bottom slope.

top_stress_x.cpp

Compute the 2D vector of the x-component of stress on the top boundary. The z boundary condition must be no slip and have no mapping along the top boundary. Implementation for mapping is missing as of 2020 (it would entail copying what is done for the bottom boundary).

$t_{x}=-\mu u_{z}$

where $\mu =\nu \rho _{0}$ is the dynamic viscosity.

top_stress_y.cpp

Compute the 2D vector of the y-component of stress on the top boundary. The z boundary condition must be no slip and have no mapping along the top boundary. Implementation for mapping is missing as of 2020 (it's would entail copying what is done for the bottom boundary).

$t_{y}=-\mu v_{z}$

where $\mu =\nu \rho _{0}$ is the dynamic viscosity.

Energy Budget

See Energy Budget for more details.

compute_Background_PE.cpp

Compute the value of the background potential energy,

$E_{b}=g\int _{V}\rho _{*}z_{*}dV$

where $z_{*}$ is the depth of an element of the sorted density field which minimized the potential energy (ie. $\rho _{*}$). For more information see Winters et al. (1995). It is assumed that gravity points in the negative z direction.

compute_BPE_from_internal.cpp

Compute the rate of energy transfer from internal energy into background potential energy.

$\phi _{i}=-g\kappa \int _{A}\rho (z=z_{u})-\rho (z=z_{l})dA$

where $z_{u}$ and $z_{l}$ are the depths of the upper and lower boundaries, respectively. The integral is taken over the horizontal extent. For more information see Winters et al. (1995). It is assumed that gravity points in the negative z direction.

dissipation.cpp

Compute the viscous dissipation $\epsilon =2\mu e_{ij}e_{ij}$

where summation is implied over duplicate indices and $e_{ij}={\frac {1}{2}}({\frac {\partial u_{i}}{\partial x_{j}}}+{\frac {\partial u_{j}}{\partial x_{i}}})$ is the strain rate tensor.

Vorticity and Vorticity Equation

See the Vorticity equation

compute_vorticity.cpp

Compute all vorticity components.

compute_vort_x.cpp

Compute the x-component of vorticity: $\omega _{x}=w_{y}-v_{z}$

compute_vort_y.cpp

Compute the y-component of vorticity: $\omega _{y}=u_{z}-w_{x}$

compute_vort_z.cpp

Compute the z-component of vorticity: $\omega _{z}=v_{x}-u_{y}$

vortex_stretch_x.cpp

Compute the x-component of the vorticity production due to stretching/tilting

$\omega _{x}u_{x}+\omega _{y}u_{y}+\omega _{z}u_{z}$

vortex_stretch_y.cpp

Compute the y-component of the vorticity production due to stretching/tilting

$\omega _{x}v_{x}+\omega _{y}v_{y}+\omega _{z}v_{z}$

vortex_stretch_z.cpp

Compute the z-component of the vorticity production due to stretching/tilting

$\omega _{x}w_{x}+\omega _{y}w_{y}+\omega _{z}w_{z}$

Enstrophy equation

See the Enstrophy equation.

enstrophy_density.cpp

Compute the enstrophy density $\Omega ={\frac {1}{2}}\omega _{i}\omega _{i}={\frac {1}{2}}(\omega _{x}^{2}+\omega _{y}^{2}+\omega _{z}^{2})$

enstrophy_stretch_production.cpp

Compute the enstrophy production due to vortex tilting/stretching $\omega _{i}\omega _{j}e_{ij}=\omega _{i}\omega _{j}{\frac {\partial u_{i}}{\partial x_{j}}}$

where summation is implied over duplicate indices and $e_{ij}={\frac {1}{2}}({\frac {\partial u_{i}}{\partial x_{j}}}+{\frac {\partial u_{j}}{\partial x_{i}}})$ is the strain rate tensor.

Coherent Vortex Identification for 3D Flows

The characteristic polynomial of $u_{i,j}$ can be written as

$\lambda ^{3}-P\lambda ^{2}+Q\lambda -R=0,$

assuming that the flow is 3D. These coefficients are invariants—i.e., scalars—of the flow, and can be used to identify coherent structures that exist in all reference frames. They are given by

$P=\lambda _{1}+\lambda _{2}+\lambda _{3}=u_{i,i},$

$Q={\frac {1}{2}}(\lambda _{1}^{2}+\lambda _{2}^{2}+\lambda _{3}^{2})={\frac {1}{2}}(P^{2}-u_{i,j}u_{j,i})={\frac {1}{2}}(P^{2}+\Omega _{ij}\Omega _{ij}-S_{ij}S_{ij}),$

$R={\frac {1}{3}}\lambda _{1}\lambda _{2}\lambda _{3}=det(u_{i,j}).$

SPINS assumes the flow is incompressible, so P=0. The more useful variables are Q and R. As can be seen from the definition of Q, large positive values of Q imply that rotation dominates while large negative values of Q imply that strain dominates; this property therefore clarifies whether a region of high vorticity is due to a vortex, or due to a shear stress, respectively. While not immediately obvious, regions of large and positive R correspond to axial stretching in the flow while large and negative R correspond to planar stretching. As a consequence, intersecting iso-surfaces of large Q and large R likely contain coherent vortices.

This discussion of Q and R is adapted from Anton Baglaenko’s PhD thesis.

Q_invt.cpp

This calculates Q using

$Q=-{\frac {1}{2}}u_{i,j}u_{j,i}$

R_invt.cpp

This calculates R using

$R=det(u_{i,j})={\frac {1}{3}}u_{i,j}u_{j,k}u_{k,i}$

Other

compute_quadweights.cpp

Compute the quadrature weights for each dimension.

If the boundary condition is no-slip, use the Clenshaw-Curtis quadrature weights:

$w=...$

If the boundary condition is free-slip, use the trapezoid rule:

$w=L_{x}/N_{x}=\Delta x$

find_expansion.cpp

Returns the expansion types for each field.

Boundary Condition	Flow Variable
	u	v	w	otherwise
Fourier	Periodic	Periodic	Periodic	Periodic
No slip	Chebyshev	Chebyshev	Chebyshev	Chebyshev
Free slip in x	sine	cosine	cosine	cosine
Free slip in y	cosine	sine	cosine	cosine
Free slip in z	cosine	cosine	sine	cosine

Spins assumes that a variable ending in _x, _y, or _z is the x, y, or z derivative of that particular variable. This is needed to get the proper expansion type for the derivative of the given variable.

get_quad.cpp

Fairly certain that this is a deprecated function

overturning_2d.cpp

Details to be added. It's an old function and probably deprecated.

read_2d_restart.cpp

For reading 2D spins files and (likely) extending to 3D.

Unsure why this is in /Science. IO should also be cleaned up since there's a bunch of overlapping functions.

read_2d_slice.cpp

For reading 2D matlab files or 2D slices.

Unsure why this is in /Science. IO should also be cleaned up since there's a bunch of overlapping functions.

swap_trig.cpp

Change cosine expansion to sine expansion and vice versa. To be used in find_expansion for variables of derivatives.

Nonlinear Equations of State

Several equations of state are included as inline functions in Science.hpp. They will take temperature, salinity, or both as arguments and there are 3 available to choose from. The first is the full nonlinear equation of state from MacDougall et. al (2003) (JAOT) which is valid for 0-45 PSU and 0-40 $^{\circ }$C. Note that the current implementation assumes that the temperature is taken to be at surface pressure (the temperature is in situ). This is not a bad simplification for SPINS, as the model is incompressible.

nleos.cpp

nleos.cpp calls a full nonlinear equation of state from MacDougall et. al (2003) (JAOT) which is valid for 0-45 PSU and 0-40 $^{\circ }$C. The user is required to supply temperature and salinity.

$\rho =\rho (T,S)$

quadeos.cpp

quadeos.cpp calls an approximation of the EOS from MacDougall et. al. (2003) (JAOT) for use with temperature only. This EOS is intended for use between 0-10 $^{\circ }$C. It takes the form

$\rho (T)=\rho (T_{md})+C(T-T_{md})^{2}$

where

$T_{md}=3.98^{\circ }C$

lineos.cpp

lineos.cpp calls a linear, two-constituent (heat/salt) equation of state. The user selects the reference temperature and salinity within the case file, and using this information, the reference density and expansion coefficients are calculated numerically from MacDougall et. al. (2003). The EOS takes the form

$\rho (T)=\rho (T_{0},S_{0})[(1+\alpha (T-T_{0})+\beta (S-S_{0})]$

where

$\alpha ={\frac {1}{\rho (T_{0},S_{0})}}{\frac {\partial \rho }{\partial T}}$

and

$\beta ={\frac {1}{\rho (T_{0},S_{0})}}{\frac {\partial \rho }{\partial S}}$