Pressure Gradient
1 Overview
The pressure gradient will be responsible for computing the horizontal gradients of both the pressure and geopotential terms for the non-Boussinesq primitive equations implemented in Omega. In the non-Boussinesq model, the conserved quantity is mass rather than volume. In Omega the prognostic variable \(\tilde{h}\) is a pseudo thickness, rather than geometric thickness as in a Boussinesq model. Some non-Boussinesq models are written in pressure coordinates (e.g. de Szoeke and Samelson 2002. However, Omega is written in general vertical coordinates and can reference either pressure \(p\) or distance \(z\) in the vertical. In a pure pressure coordinate the pressure gradient term disappears (since the pressure does not vary along lines of constant pressure), just as how the geopotential term disappears in a pure z coordinate model. However, similar to MPAS-Ocean’s support for tilted height coordinates, Omega will allow for tilted pressure coordinates. This means that Omega will need to compute both the pressure and geopotential gradients.
2 Requirements
2.1 Requirement: Support for tilted pressure coordinates for the non-Boussinesq primitive equations
The pressure gradient will compute the horizontal gradients of both the pressure and geopotential to support tilted pressure coordinates in the non-Boussinesq model. This will allow for the use of a \(p^\star\) coordinate, which functions similarly to the \(z^\star\) in the Boussinesq MPAS-Ocean model. Note that we use the name “\(p^\star\)” to refer to the vertical coordinate, in Omega it will be expressed in terms of the pseudo-height, \(\tilde{z}\), as opposed to pressure directly.
2.2 Requirement: Initial support for a simple centered pressure gradient
For initial global cases without ice shelf cavities, the pressure and geopotential gradients will be computed with a simple centered difference approximation. In later versions of Omega, one or more high-order pressure gradients may be implemented to replace the centered approach in production runs. However, the centered pressure gradient will remain an option for use in idealized testing.
2.3 Requirement: Flexibility to support a high-order pressure gradient
The centered pressure gradient will be insufficient for future versions of Omega that include ice shelf cavities and high resolution shelf breaks. The pressure gradient framework should be flexible enough to support a high-order pressure gradient in the future. The high order pressure gradient will be similar to Adcroft et al. 2008.
2.4 Requirement: Flexibility to support tidal forcing and sea level change
In later versions of Omega, the pressure gradient will need to be able to include tidal forcing in the geopotential term.
These tidal forcings include both the tidal potential and the self attraction and loading terms.
Additionally, other long-term changes to the geoid can be included in the geopotential.
This requirement is satisfied by tidal forcing terms being included in the geopotential calculation in the VertCoord class.
2.5 Desired: Pressure gradient for barotropic mode
For split barotropic-baroclinic timestepping, the pressure gradient should provide the bottom pressure gradient tendency in the barotropic mode. The details of the barotropic pressure gradient will be added in a future design document for split time stepping.
3 Algorithmic Formulation
In the layered non-Boussinesq equation solved in Omega, the pressure and geopotential gradient terms are taken along a generalized coordinate \(r\):
3.1 Centered Pressure Gradient
The tendency term for edge \(e\) and layer \(k\), \(T^p_{e,k}\) is discretized as:
3.2 High-order Pressure Gradient
The high order pressure gradient will be based on the full volume integral form of the geopotential and pressure terms:
4 Design
4.1 Data types and parameters
The PressureGradient class will be used to perform the horizontal gradients of the pressure and geopotential
class PressureGrad{
public:
private:
// Map of all pressure gradients
static std::map<std::string, std::unique_ptr<PressureGrad>> AllPGrads;
// Instances of functors
PressureGradCentered CenteredPGrad;
PressureGradHighOrder HighOrderPGrad1; // To be implemented later
PressureGradHighOrder HighOrderPGrad2; // Multiple high order options are likely in the future
// Pressure gradient choice from config
PressureGradType PressureGradChoice;
// Data required for computation (stored copies of HorzMesh/VCoord arrays)
Array2DI4 CellsOnEdge;
Array1DReal DcEdge;
Array1DReal EdgeSignOnCell;
Array1DI4 MinLayerEdgeBot;
Array1DI4 MaxLaterEdgeTop;
}
The functions to compute the centered and high order pressure gradient terms will be implemented as functors and the pressure gradient class will have private instances of these classes.
The PressureGrad class will store pointers to the static mesh and vertical coordinate needed in the calculation of the pressure gradient, which will be initialized on construction.
The class will support multiple named instances to allow for the possibility of a barotropic pressure gradient calculation on a different mesh in the future.
The user will select a pressure gradient option at runtime in the input yaml file under the pressure gradient section:
PressureGrad:
PressureGradType: 'centered'
An enum class will be used to specify options for the pressure gradient used for an Omega simulation:
enum class PressureGradType{
Centered,
HighOrder1,
HighOrder2
}
4.2 Methods
The PressureGrad class will have methods for creating, retrieving, and removing instances similar to other Omega classes such as Dcomp, HorzMesh, etc.
Computation of the pressure gradient will be performed by computePressureGrad, which will select the user’s chosen implementation from the available options.
It will also have a helper method to pre-compute the interface term, \(\gamma_{e,k-1/2}\), needed for the centered pressure gradient to separate layer interface and midpoint calculations.
class PressureGrad{
public:
// Creation methods
static PressureGrad *init();
static PressureGrad *create(const std::string &Name, const HorzMesh *Mesh,
const VertCoord *VCoord, Config *Options);
// Retrival methods
static PressureGrad *get(const std::string &Name);
static PressureGrad *getDefault();
// Removal methods
static void clear();
static void erase(const std::string &Name);
~PressureGrad();
// Main compute method
void computePressureGrad(Array2DReal &Tend, const OceanState *State,
const VertCoord *VCoord, const Eos *EqState,
const int TimeLevel) const;
}
4.2.1 Creation
The create method will be responsible for calling the constructor which will, store any static mesh information as private variables and handle the retrieval of the user-specified pressure gradient option.
PressureGrad *PressureGrad::create(const std::string &Name,
const HorzMesh *Mesh,
const VertCoord *VCoord,
Config *Options);
4.2.2 Initialization
The init method will create the default pressure gradient and return an error code:
static PressureGrad::init();
It calls the constructor to create a new pressure gradient instance, producing a static managed (unique) pointer to the single instance.
4.2.3 Retrieval
There will be a method for getting a pointer to the default pressure gradient instance:
PressureGrad *PressureGrad::getDefault();
or a named instance:
PressureGrad *PressureGrad::get(const std::string &Name);
4.2.4 Computation
The public computePressureGrad method will rely on private methods for each specific pressure gradient option (centered and high order).
Note that the functors called by computePressureGrad are responsible for computing the sum of the pressure gradient and geopotential gradient accumulated in the Tend output array.
void PressureGrad::computePressureGrad(const Array2DReal &Tend,
const OceanState *State,
const VertCoord *VCoord,
const Eos *EqState,
const I4 TimeLevel) {
OMEGA_SCOPE(LocCenteredPGrad, CenteredPGrad);
OMEGA_SCOPE(LocHighOrderPGrad, HighOrderPGrad);
OMEGA_SCOPE(LocMinLayerEdgeBot, MinLayerEdgeBot);
OMEGA_SCOPE(LocMaxLayerEdgeTop, MaxLayerEdgeTop);
const Array2DReal &PressureMid = VCoord->PressureMid;
const Array2DReal &PressureInterface = VCoord->PressureInterface;
const Array2DReal &Geopotential = VCoord->GeopotentialMid;
const Array2DReal &SpecVol = EqState->SpecVol;
const Array2DReal &ZInterface = VCoord->ZInterface;
Array2DReal LayerThick;
State->getLayerThickness(LayerThick, TimeLevel);
if (PressureGradChoice == PressureGradType::Centered) {
// computes centered geopotential and pressure gradient tendency
parallelForOuter(
"pgrad-centered", {NEdgesAll},
KOKKOS_LAMBDA(I4 IEdge, const TeamMember &Team) {
const int KMin = LocMinLayerEdgeBot(IEdge);
const int KMax = LocMaxLayerEdgeTop(IEdge);
const int KRange = vertRange(KMin, KMax);
parallelForInner(
Team, KRange,
INNER_LAMBDA(int KChunk) {
LocCenteredPGrad(Tend, IEdge, KChunk, PressureMid,
PressureInterface, ZInterface,
LocTidalPotential, LocSelfAttractionLoading,
SpecVol);
});
});
} else {
// computes high-order geopotential and pressure gradient tendency
parallelForOuter(
"pgrad-highorder", {NEdgesAll},
KOKKOS_LAMBDA(I4 IEdge, const TeamMember &Team) {
const int KMin = MinLayerEdgeBot(IEdge);
const int KMax = MaxLayerEdgeTop(IEdge);
const int KRange = vertRange(KMin, KMax);
parallelForInner(
Team, KRange,
INNER_LAMBDA(int KChunk) {
LocHighOrderPGrad(Tend, IEdge, KChunk, PressureMid,
Geopotential, SpecVol);
});
});
}
} // end compute pressure gradient
4.2.5 Destruction and removal
No operations are needed in the destructor. The there will be a method to remove all pressure gradient instances,
static void PressureGrad::clear();
or a named instance
static void erase(const std::string &Name);
Verification and Testing
Unit Test: Hydrostatic test with tilted layers
The pressure gradient term will be computed on an edge with two adjacent cells that have the same sea surface height, but internal layers that are tilted in terms of geometric height. The same linear temperature and salinity profiles will be used for both cells, thus the pressure gradient should be zero on the edge. The pseudo-height (and specific volume) of the layers will be computed via iteration, where the pressure is updated in the EOS evaluation until convergence is reached. The pressure gradient term will be compared to the zero exact solution at progressively finer resolutions to evaluate convergence.
Test: Spatial convergence to exact solution
For given analytical functions of \(T\), and \(S\) and with specified \(z\) layer spacing, the spatial convergence of the pressure gradient can be assessed by computing errors using a high-fidelity reference solution on progressively finer meshes.
Test: Seamount test
The seamount test in Polaris will be used to verify the pressure gradient’s ability to preserve the resting state of fluid in a case with tilted layers.
Test: Baroclinic gyre
The baroclinic gyre test case will test the pressure gradient term in the full non-Boussinesq equations.