From 1a0ca0077e424fa0029c48a260aa49231c7ad3c4 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Tue, 22 Apr 2025 14:24:26 -0700 Subject: [PATCH 1/9] Add pressure gradient design document --- components/omega/doc/design/PGrad.md | 197 +++++++++++++++++++++++++++ components/omega/doc/index.md | 1 + 2 files changed, 198 insertions(+) create mode 100644 components/omega/doc/design/PGrad.md diff --git a/components/omega/doc/design/PGrad.md b/components/omega/doc/design/PGrad.md new file mode 100644 index 000000000000..23d692bae3dc --- /dev/null +++ b/components/omega/doc/design/PGrad.md @@ -0,0 +1,197 @@ +(omega-design-pressure-grad)= +# 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 primative equations implemented in Omega. +In the non-Boussinesq model, the vertical coordinate will be pressure as opposed to height. +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. + +### 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 ore more high-order pressure gradients will be implemented and will 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. + +### 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 changes to the geoid + +### 2.5 Requirement: Pressure gradient for barotropic mode + +For split barotropic-baroclinic timestepping, the pressure gradient should provide the bottom pressure gradient tendency in the barotropic mode. + +### Desired: + +## 3 Algorithmic Formulation +The non-Boussinesq momentum equation is +$$ \frac{D \mathbf{u}_h}{D t } + f\boldsymbol{k}\times \mathbf{u}_h + \left(v\nabla_A p + \nabla_A \phi \right) = \boldsymbol{\mathcal{F}}. $$ + +where $\mathbf{u}_h$ is the horizontal velocity, $f$ is the Coriolis parameter, $v = \frac{1}{\rho}$ is the specific volume, $\rho = \rho(T,S,p)$ is the density, $p$ is the hydrostatic pressure, $\phi$ is the geopotential, and $\boldsymbol{\mathcal{F}}$ are the dissipative terms. +The operator $\nabla_A$ is the gradient along a constant surface, $A$, and the total derivative is + +$$ \frac{D \mathbf{u}_h}{D t } = \left( \frac{\partial}{\partial t} \right)_A + \mathbf{u}_h\cdot \nabla_A + \omega\frac{\partial}{\partial A}, $$ + +where $\omega$ is the cross coordinate flow. +In the layered non-Boussinesq equations, the prognostic variable is the pressure thickness $h_k$, so that the geometric thickness (in meters) is a diagnostic variable defined as: + +$$ \Delta z_k = v_k h_k. $$ + +The pressure at vertical cell interfaces is the found by summing the pressure thicknesses: + +$$ p_{K+1/2} = p_{surf} + g\sum_{k=1}^K h_k. $$ + +The geopotential at vertical cell integrates is found by summing the pressure thicknesses multiplied by the specific volume: + +$$ \phi_{K+1/2} = g\left(z_b + \sum_{k=K}^{N} v_k h_k\right), $$ + +where $z_b$ is the (positive-up) bottom depth. + +The discrete gradient operator at an edge is: + +$$ \nabla {(\cdot)} = \frac{1}{d_e} \sum_{i\in CE(e)} -n_{e,i} (\cdot)_i $$ + +where $d_e$ is the distance between cell centers, $CE(e)$ are the cells on edge $e$, and $n_{e,i}$ is the sign of the edge normal with respect to cell $i$. +Therefore the centered pressure gradient will be calculated as: + +$$ T^p_k = \frac{1}{d_e} \left( \widehat{v}_{k,e} \sum_{i \in CE(e)} n_{e,i} \overline{p}_{k,i} + \sum_{i\in CE(e)} n_{e,i} \overline{\phi}_{k,i}\right), $$ + +$$ = \frac{1}{d_e} \left( \sum_{i \in CE(e)} n_{e,i} \left( \widehat{v}_{k,e}\overline{p}_{k,i} + \overline{\phi}_{k,i} \right) \right), $$ + +with the vertical averaging operator defined as: + +$$ \overline{(\cdot)} = \frac{1}{2} \left((\cdot)_{k+1/2} + (\cdot)_{k-1/2} \right)$$ + +and the horizontal averaging operator: + +$$\widehat{(\cdot)} = \frac{1}{2} \left( (\cdot)_{i=1} + (\cdot)_{i=2}\right)$$ + +## 4 Design + +### 4.1 Data types and parameters + +The `PressureGradient` class will be used to perform the horizontal gradients of the pressure and geopotential +```c++ +class PressureGrad{ + public: + private: + PressureGradCentered CenteredPGrad; + PressureGradHighOrder HighOrderPGrad; + PressureGradType PressureGradChoice; + I4 NVertLevels; + I4 NChuncks; + Array2DI4 CellsOnEdge; + Array1DReal DvEdge; + Array1DReal EdgeSignOnCell; +} +``` +The user will select a pressure gradient option at runtime in the input yaml file under the pressure gradient section: +```yaml + PressureGrad: + PressureGradType: 'centered' +``` +An `enum class` will be used to specify options for the pressure gradient used for an Omega simulation: +```c++ +enum class PressureGradType{ + Centered, + HighOrder +} +``` +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. +### 4.2 Methods + +```c++ +class PressureGrad{ + public: + static PressureGrad *create(); + void computePressureGrad(); + void computePressureGradBtr(); + private: + +} +``` +#### 4.2.1 Creation + +The constructor will be responsible for storing any static mesh information as private variables and handling the selection of the user-specified pressure gradient option. +```c++ +PressureGrad::PressureGrad(const HorzMesh *Mesh, int NVertLevels, Config *Options); +``` + +The create method will take the same arguments as the constructor plus a name. +It calls the constructor to create a new pressure gradient instance, and put it in the static map of all pressure gradients. +It will return a pointer to the newly created object. +```c++ +PressureGrad *PressureGrad::create(const std::string &Name, const HorzMesh *Mesh, int NVertLevels, Config *Options); +``` + +#### 4.2.2 Initialization + +The init method will create the default pressure gradient and return an error code: +```c++ +int PressureGrad::init(); +``` + +#### 4.2.3 Retrieval + +There will be methods for getting the default and non-default pressure gradient instances: +```c++ +PressureGrad *PressureGrad::getDefault(); +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). +```c++ +void PressureGrad::computePressureGrad(const Array2DReal &Tend, + const Array2DReal &Pressure, + const Array2DReal &Geopotential, + const Array2DReal &SpecVol) { +OMEGA_SCOPE(LocCenteredPGrad, CenteredPGrad) +OMEGA_SCOPE(LocHighOrderPGrad, HighOrderPGrad) + if (PressureGradChoice == PressureGradType::Centered){ + parallelFor("pgrad-centered", {NCellsAll, NChunks}, + KOKKOS_LAMBDA(int ICell, int KChunk) { + LocCenteredPGrad(Tend, Pressure, Geopotential, SpecVol); + }); + } + else if (PressureGradChoice == PressureGradType::HighOrder){ + parallelFor("pgrad-highorder", {NCellsAll, NChunks}, + KOKKOS_LAMBDA(int ICell, int KChunk) { + LocHighOrderPGrad(Tend, Pressure, Geopotential, SpecVol); + }); + } +``` + +#### 4.2.5 Destruction and removal + +No operations are needed in the destructor. +The erase method will remove a named pressure gradient instance, whereas the clear method will remove all of +them. +Both will call the destructor in the process. +```c++ +void PressureGrad::erase(const std::string &Name); +void PressureGrad::clear(); +``` + +## Verification and Testing + +### Test: Spatial convergence to exact solution +For a given analytical $v$, $h$, and $b$, the spatial convergence of the pressure gradient can be assessed by computing errors on progressively finer meshes. + +### Test: Baroclinic gyre +The baroclinic gyre test case will test the pressure gradient term in the full non-Boussinesq equations. + diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md index 1dd11e400a43..ab948c4ea962 100644 --- a/components/omega/doc/index.md +++ b/components/omega/doc/index.md @@ -116,6 +116,7 @@ design/HorzMeshClass design/Logging design/MachEnv design/Metadata +design/PGrad design/IO design/IOStreams design/Reductions From 1e61c0f5cd60e981e32dabfcb18fe9c2c9c3d644 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 23 Apr 2025 09:57:35 -0500 Subject: [PATCH 2/9] Fix some lint --- components/omega/doc/design/PGrad.md | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/components/omega/doc/design/PGrad.md b/components/omega/doc/design/PGrad.md index 23d692bae3dc..c75649c5a83e 100644 --- a/components/omega/doc/design/PGrad.md +++ b/components/omega/doc/design/PGrad.md @@ -13,12 +13,12 @@ This means that Omega will need to compute both the pressure and geopotential gr ### 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. +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. ### 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 ore more high-order pressure gradients will be implemented and will replace the centered approach in production runs. +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 ore more high-order pressure gradients will be implemented and will 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 @@ -28,7 +28,7 @@ The pressure gradient framework should be flexible enough to support a high-orde ### 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 changes to the geoid +Additionally, other changes to the geoid ### 2.5 Requirement: Pressure gradient for barotropic mode @@ -106,8 +106,8 @@ The user will select a pressure gradient option at runtime in the input yaml fil An `enum class` will be used to specify options for the pressure gradient used for an Omega simulation: ```c++ enum class PressureGradType{ - Centered, - HighOrder + Centered, + HighOrder } ``` 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. @@ -156,7 +156,7 @@ PressureGrad *PressureGrad::get(const std::string &Name); The public `computePressureGrad` method will rely on private methods for each specific pressure gradient option (centered and high order). ```c++ -void PressureGrad::computePressureGrad(const Array2DReal &Tend, +void PressureGrad::computePressureGrad(const Array2DReal &Tend, const Array2DReal &Pressure, const Array2DReal &Geopotential, const Array2DReal &SpecVol) { @@ -167,7 +167,7 @@ OMEGA_SCOPE(LocHighOrderPGrad, HighOrderPGrad) KOKKOS_LAMBDA(int ICell, int KChunk) { LocCenteredPGrad(Tend, Pressure, Geopotential, SpecVol); }); - } + } else if (PressureGradChoice == PressureGradType::HighOrder){ parallelFor("pgrad-highorder", {NCellsAll, NChunks}, KOKKOS_LAMBDA(int ICell, int KChunk) { @@ -180,7 +180,7 @@ OMEGA_SCOPE(LocHighOrderPGrad, HighOrderPGrad) No operations are needed in the destructor. The erase method will remove a named pressure gradient instance, whereas the clear method will remove all of -them. +them. Both will call the destructor in the process. ```c++ void PressureGrad::erase(const std::string &Name); @@ -194,4 +194,3 @@ For a given analytical $v$, $h$, and $b$, the spatial convergence of the pressu ### Test: Baroclinic gyre The baroclinic gyre test case will test the pressure gradient term in the full non-Boussinesq equations. - From dc19bdfa55afcafe151ca78929e4f0f3e133673d Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Thu, 19 Jun 2025 13:10:01 -0700 Subject: [PATCH 3/9] Revise based on review comments --- components/omega/doc/Makefile | 2 +- components/omega/doc/design/PGrad.md | 120 +++++++++++++++++---------- 2 files changed, 75 insertions(+), 47 deletions(-) diff --git a/components/omega/doc/Makefile b/components/omega/doc/Makefile index 8539aa59fb0b..ffd39dbd1638 100644 --- a/components/omega/doc/Makefile +++ b/components/omega/doc/Makefile @@ -6,7 +6,7 @@ SPHINXOPTS ?= SPHINXBUILD ?= sphinx-build SOURCEDIR = . -BUILDDIR = _build +BUILDDIR = /global/cfs/cdirs/e3sm/www/sbrus/pgrad # Put it first so that "make" without argument is like "make help". help: diff --git a/components/omega/doc/design/PGrad.md b/components/omega/doc/design/PGrad.md index c75649c5a83e..742d85456212 100644 --- a/components/omega/doc/design/PGrad.md +++ b/components/omega/doc/design/PGrad.md @@ -3,8 +3,11 @@ ## 1 Overview -The pressure gradient will be responsible for computing the horizontal gradients of both the pressure and geopotential terms for the non-Boussinesq primative equations implemented in Omega. -In the non-Boussinesq model, the vertical coordinate will be pressure as opposed to height. +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 in m as in a Boussinesq model. + Some non-Boussinesq models are written in pressure coordinates (e.g. [de Szoeke and Samelson 2002](https://journals.ametsoc.org/view/journals/phoc/32/7/1520-0485_2002_032_2194_tdbtba_2.0.co_2.xml). + 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. @@ -18,66 +21,99 @@ This will allow for the use of a $p^\star$ coordinate, which functions similarly ### 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 ore more high-order pressure gradients will be implemented and will replace the centered approach in production runs. +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 will be implemented and will 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](https://doi.org/10.1016/j.ocemod.2008.02.001). ### 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 changes to the geoid +Additionally, other long-term changes to the geoid can be included in the geopotential. -### 2.5 Requirement: Pressure gradient for barotropic mode +### 2.5 Disired: Pressure gradient for barotropic mode For split barotropic-baroclinic timestepping, the pressure gradient should provide the bottom pressure gradient tendency in the barotropic mode. +This will be added in a future version when split time stepping is implemented. ### Desired: ## 3 Algorithmic Formulation -The non-Boussinesq momentum equation is -$$ \frac{D \mathbf{u}_h}{D t } + f\boldsymbol{k}\times \mathbf{u}_h + \left(v\nabla_A p + \nabla_A \phi \right) = \boldsymbol{\mathcal{F}}. $$ +### 3.1 Centered Pressure Gradient +In the layered non-Boussinesq [momentum equation](OmegaV1GoverningEqns.md#discrete-momentum) solved in Omega, the pressure gradient tendency term for edge $e$ and level $k$, $T^p_{e,k}$, includes the gradient of the pressure and the gradient of the geopotential, -where $\mathbf{u}_h$ is the horizontal velocity, $f$ is the Coriolis parameter, $v = \frac{1}{\rho}$ is the specific volume, $\rho = \rho(T,S,p)$ is the density, $p$ is the hydrostatic pressure, $\phi$ is the geopotential, and $\boldsymbol{\mathcal{F}}$ are the dissipative terms. -The operator $\nabla_A$ is the gradient along a constant surface, $A$, and the total derivative is +$$ +T^p_{e,k} = -\left[ \alpha_{i,k} \right]_e \nabla p_{i,k} - \nabla \Phi_{i,k}, +$$ -$$ \frac{D \mathbf{u}_h}{D t } = \left( \frac{\partial}{\partial t} \right)_A + \mathbf{u}_h\cdot \nabla_A + \omega\frac{\partial}{\partial A}, $$ +where the second term is necessary to account for tilted layers that occur when using a general vertical coordinate. +In this equation, $\alpha_{i,k}$ is the specific volume for cell $i$ at the mid-point of level $k$, $p_{i,k}$ is the pressure, and $\Phi_{i,k}$ is the geopotential. +The discrete gradient operator at an edge is: -where $\omega$ is the cross coordinate flow. -In the layered non-Boussinesq equations, the prognostic variable is the pressure thickness $h_k$, so that the geometric thickness (in meters) is a diagnostic variable defined as: +$$ + \nabla {(\cdot)} = \frac{1}{d_e} \sum_{i\in CE(e)} -n_{e,i} (\cdot)_i +$$ -$$ \Delta z_k = v_k h_k. $$ +where $d_e$ is the distance between cell centers, $CE(e)$ are the cells on edge $e$, and $n_{e,i}$ is the sign of the edge normal with respect to cell $i$. +The horizontal averaging operator is: -The pressure at vertical cell interfaces is the found by summing the pressure thicknesses: +$$ + [\cdot]_e = \frac{1}{2}\sum_{i\in CE(e)} (\cdot)_i +$$ -$$ p_{K+1/2} = p_{surf} + g\sum_{k=1}^K h_k. $$ +Therefore, the centered pressure gradient will be calculated as: -The geopotential at vertical cell integrates is found by summing the pressure thicknesses multiplied by the specific volume: +$$ + T^p_{e,k} = \frac{1}{d_e} \left( [\alpha_{i,k}]_e \sum_{i \in CE(e)} n_{e,i} p_{k,i} + \sum_{i\in CE(e)} n_{e,i} \Phi_{k,i}\right), +$$ -$$ \phi_{K+1/2} = g\left(z_b + \sum_{k=K}^{N} v_k h_k\right), $$ +$$ + = \frac{1}{d_e} \left( \sum_{i \in CE(e)} n_{e,i} \left( [\alpha_{i,k}]_ep_{k,i} + \Phi_{k,i} \right) \right), +$$ -where $z_b$ is the (positive-up) bottom depth. +### 3.2 Barotropic Pressure Gradient -The discrete gradient operator at an edge is: +When split baroclinic-barotropic time stepping is implemented in the future, the barotropic pressure gradient will be calculated by the pressure gradient class. +The barotropic pressure gradient is found by depth integrating the pressure gradient. +The pressure is -$$ \nabla {(\cdot)} = \frac{1}{d_e} \sum_{i\in CE(e)} -n_{e,i} (\cdot)_i $$ +$$ +p(z) = p_b - g \int^z_{-h} \rho dz^\prime, +$$ -where $d_e$ is the distance between cell centers, $CE(e)$ are the cells on edge $e$, and $n_{e,i}$ is the sign of the edge normal with respect to cell $i$. -Therefore the centered pressure gradient will be calculated as: +where $p_b$ is the bottom pressure. +The bottom pressure is the sum of the atmospheric surface pressure, $p_s$, and the pressure contribution of the water column: + +$$ +p(z) &= p_s + g\int_{-h}^\eta \rho dz - g \int^z_{-h} \rho dz^\prime, \\ + &= p_s + g\rho_0\widetilde{H} - g \int^z_{-h} \rho dz^\prime, +$$ -$$ T^p_k = \frac{1}{d_e} \left( \widehat{v}_{k,e} \sum_{i \in CE(e)} n_{e,i} \overline{p}_{k,i} + \sum_{i\in CE(e)} n_{e,i} \overline{\phi}_{k,i}\right), $$ +where the total water column pseudo height is expressed by -$$ = \frac{1}{d_e} \left( \sum_{i \in CE(e)} n_{e,i} \left( \widehat{v}_{k,e}\overline{p}_{k,i} + \overline{\phi}_{k,i} \right) \right), $$ +$$ +\widetilde{H} = \int_{-h}^\eta \frac{\rho}{\rho_0} dz. +$$ -with the vertical averaging operator defined as: +$\widetilde{H}$ is the prognositc variable in the barotropic continuity equation. +The vertical integral of the pressure gradient is -$$ \overline{(\cdot)} = \frac{1}{2} \left((\cdot)_{k+1/2} + (\cdot)_{k-1/2} \right)$$ +$$ +\frac{1}{\rho_0\widetilde{H}}\int^\eta_{-h} \nabla p dz &= \frac{1}{\rho_0\widetilde{H}}\int^\eta_{-h} \nabla \left( p_s + g\rho_0 \widetilde{H} - g \int^z_{-h} \rho dz^\prime \right) dz, \\ + &= \frac{H}{\rho_0\widetilde{H}}\nabla p_s + \frac{gH}{\widetilde{H}}\nabla \widetilde{H} - \frac{g}{\rho_0\widetilde{H}} \int_{-h}^\eta \left( \nabla \int_{-h}^z \rho dz^\prime\right) dz, +$$ -and the horizontal averaging operator: +where the height of the water column is represented by $H$. +The $1/\rho_0\widetilde{H}$ factor comes vertically integrating the material derivative and expressing the resulting barotropic momentum equation in non-conservative form. -$$\widehat{(\cdot)} = \frac{1}{2} \left( (\cdot)_{i=1} + (\cdot)_{i=2}\right)$$ +Therefore, the barotorpic pressure gradient term is discretized as: + +$$ +\overline{T}_e^p = g\left[ \frac{H_i}{\widetilde{H}_i} \right]_e\sum_{i \in CE(e)} n_{e,i}\widetilde{H}_e +$$ ## 4 Design @@ -88,8 +124,9 @@ The `PressureGradient` class will be used to perform the horizontal gradients of class PressureGrad{ public: private: + std::unique_ptr OmegaPressureGrad; PressureGradCentered CenteredPGrad; - PressureGradHighOrder HighOrderPGrad; + PressureGradHighOrder HighOrderPGrad; // To be implemented later PressureGradType PressureGradChoice; I4 NVertLevels; I4 NChuncks; @@ -116,9 +153,9 @@ The functions to compute the centered and high order pressure gradient terms wil ```c++ class PressureGrad{ public: - static PressureGrad *create(); + static PressureGrad *init(); + static PressureGrad *get(); void computePressureGrad(); - void computePressureGradBtr(); private: } @@ -130,31 +167,25 @@ The constructor will be responsible for storing any static mesh information as p PressureGrad::PressureGrad(const HorzMesh *Mesh, int NVertLevels, Config *Options); ``` -The create method will take the same arguments as the constructor plus a name. -It calls the constructor to create a new pressure gradient instance, and put it in the static map of all pressure gradients. -It will return a pointer to the newly created object. -```c++ -PressureGrad *PressureGrad::create(const std::string &Name, const HorzMesh *Mesh, int NVertLevels, Config *Options); -``` - #### 4.2.2 Initialization -The init method will create the default pressure gradient and return an error code: +The init method will create the pressure gradient and return an error code: ```c++ int 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 methods for getting the default and non-default pressure gradient instances: +There will be a method for getting a pointer to the pressure gradient instance: ```c++ -PressureGrad *PressureGrad::getDefault(); -PressureGrad *PressureGrad::get(const std::string &Name); +PressureGrad *PressureGrad::get(); ``` #### 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. ```c++ void PressureGrad::computePressureGrad(const Array2DReal &Tend, const Array2DReal &Pressure, @@ -179,11 +210,8 @@ OMEGA_SCOPE(LocHighOrderPGrad, HighOrderPGrad) #### 4.2.5 Destruction and removal No operations are needed in the destructor. -The erase method will remove a named pressure gradient instance, whereas the clear method will remove all of -them. -Both will call the destructor in the process. +The clear method will remove the pressure gradient instance, calling the destructor in the process. ```c++ -void PressureGrad::erase(const std::string &Name); void PressureGrad::clear(); ``` From 1de4d88eadca0f911021cc62ed5b66250a65bcfe Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Tue, 16 Sep 2025 14:48:41 -0700 Subject: [PATCH 4/9] Revisions based on V1 equations --- components/omega/doc/Makefile | 2 +- .../omega/doc/design/OmegaV1GoverningEqns.md | 2 + components/omega/doc/design/PGrad.md | 167 ++++++++++-------- 3 files changed, 100 insertions(+), 71 deletions(-) diff --git a/components/omega/doc/Makefile b/components/omega/doc/Makefile index ffd39dbd1638..8539aa59fb0b 100644 --- a/components/omega/doc/Makefile +++ b/components/omega/doc/Makefile @@ -6,7 +6,7 @@ SPHINXOPTS ?= SPHINXBUILD ?= sphinx-build SOURCEDIR = . -BUILDDIR = /global/cfs/cdirs/e3sm/www/sbrus/pgrad +BUILDDIR = _build # Put it first so that "make" without argument is like "make help". help: diff --git a/components/omega/doc/design/OmegaV1GoverningEqns.md b/components/omega/doc/design/OmegaV1GoverningEqns.md index 2262a437fe6c..fedd0e3fa2a5 100644 --- a/components/omega/doc/design/OmegaV1GoverningEqns.md +++ b/components/omega/doc/design/OmegaV1GoverningEqns.md @@ -605,6 +605,7 @@ $$ (vh-momentum-reynolds1) Here we have also moved the resolved component projection operator through the spatial integrals given the properties of the operator. Next we decompose the nonlinear products into resolved and unresolved components, this yields +(omega-v1-vh-momentum-reynolds2)= $$ \frac{d}{dt} \int_A \int_{\tilde{z}_k^{\text{bot}}}^{\tilde{z}_k^{\text{top}}} \left< {\bf u} \right> \, d\tilde{z} \, dA & + \int_{\partial A} \left( \int_{\tilde{z}_k^{\text{bot}}}^{\tilde{z}_k^{\text{top}}} \left( \left< {\bf u} \right> \otimes \left< {\bf u} \right> + \left< {\bf u}^\prime \otimes {\bf u}^\prime \right> \right) \, d\tilde{z} \right) \cdot \mathbf{n}_\perp dl \, dl \\ @@ -780,6 +781,7 @@ In the tracer equation, we note that surface fluxes (e.g. latent heat fluxes) wi Omega will only predict the layer average normal velocity, so we drop the bold face on the $u$ terms except for the product of primes, which is specified in the next section. +(omega-v1-momentum-eq)= $$ \frac{\partial u_{e,k}}{\partial t} & + \left[ {\bf k} \cdot \nabla \times u_{e,k} +f_v\right]_e\left(u_{e,k}^{\perp}\right) + \left[\nabla K\right]_e \\ diff --git a/components/omega/doc/design/PGrad.md b/components/omega/doc/design/PGrad.md index 742d85456212..72fc1190b81a 100644 --- a/components/omega/doc/design/PGrad.md +++ b/components/omega/doc/design/PGrad.md @@ -5,7 +5,7 @@ ## 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 in m as in a Boussinesq model. +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](https://journals.ametsoc.org/view/journals/phoc/32/7/1520-0485_2002_032_2194_tdbtba_2.0.co_2.xml). 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. @@ -18,6 +18,7 @@ This means that Omega will need to compute both the pressure and geopotential gr 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 @@ -33,24 +34,24 @@ The high order pressure gradient will be similar to [Adcroft et al. 2008](https: 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 Disired: Pressure gradient for barotropic mode +### 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. -This will be added in a future version when split time stepping is implemented. - -### Desired: +The details of the barotropic pressure gradient will be added in a future design document for split time stepping. ## 3 Algorithmic Formulation ### 3.1 Centered Pressure Gradient -In the layered non-Boussinesq [momentum equation](OmegaV1GoverningEqns.md#discrete-momentum) solved in Omega, the pressure gradient tendency term for edge $e$ and level $k$, $T^p_{e,k}$, includes the gradient of the pressure and the gradient of the geopotential, +In the layered non-Boussinesq {ref}`momentum equation ` solved in Omega, the pressure gradient tendency term for edge $e$ and layer $k$, $T^p_{e,k}$, includes the gradient of the geopotential, the gradient of a term involving pressure, and two terms evaluated at the cell interface: $$ -T^p_{e,k} = -\left[ \alpha_{i,k} \right]_e \nabla p_{i,k} - \nabla \Phi_{i,k}, +T^p_{e,k} = - \left(\nabla \Phi \right)_{e,k} - \frac{1}{\left[\tilde{h}_k\right]_e} \nabla \left( \tilde{h}_k \alpha_k p_k \right) + \frac{1}{\left[\tilde{h}_k\right]_e} \left\{ \left[ \alpha p \nabla \tilde{z}\right]_{e,k}^\text{top} - \left[ \alpha p \nabla \tilde{z}\right]_{e,k+1}^\text{top} \right\}. $$ -where the second term is necessary to account for tilted layers that occur when using a general vertical coordinate. -In this equation, $\alpha_{i,k}$ is the specific volume for cell $i$ at the mid-point of level $k$, $p_{i,k}$ is the pressure, and $\Phi_{i,k}$ is the geopotential. +The geopotential and interface terms are necessary to account for tilted layers that occur when using a general vertical coordinate, where the gradient operator is taken along layers. +In this equation, $\alpha_{i,k}$ specific volume, $p_{i,k}$ is the pressure, and $\Phi_{i,k}$ is the geopotential. +These three quantities are evaluated at the mid-point of layer $k$ of cell $i$ in the first two terms and at the cell interfaces in the third term along with the interface psudo-height, \tilde{z}. The discrete gradient operator at an edge is: $$ @@ -58,7 +59,7 @@ $$ $$ where $d_e$ is the distance between cell centers, $CE(e)$ are the cells on edge $e$, and $n_{e,i}$ is the sign of the edge normal with respect to cell $i$. -The horizontal averaging operator is: +The (cell-to-edge) horizontal averaging operator is: $$ [\cdot]_e = \frac{1}{2}\sum_{i\in CE(e)} (\cdot)_i @@ -67,53 +68,69 @@ $$ Therefore, the centered pressure gradient will be calculated as: $$ - T^p_{e,k} = \frac{1}{d_e} \left( [\alpha_{i,k}]_e \sum_{i \in CE(e)} n_{e,i} p_{k,i} + \sum_{i\in CE(e)} n_{e,i} \Phi_{k,i}\right), -$$ - -$$ - = \frac{1}{d_e} \left( \sum_{i \in CE(e)} n_{e,i} \left( [\alpha_{i,k}]_ep_{k,i} + \Phi_{k,i} \right) \right), -$$ - -### 3.2 Barotropic Pressure Gradient - -When split baroclinic-barotropic time stepping is implemented in the future, the barotropic pressure gradient will be calculated by the pressure gradient class. -The barotropic pressure gradient is found by depth integrating the pressure gradient. -The pressure is - -$$ -p(z) = p_b - g \int^z_{-h} \rho dz^\prime, -$$ - -where $p_b$ is the bottom pressure. -The bottom pressure is the sum of the atmospheric surface pressure, $p_s$, and the pressure contribution of the water column: - -$$ -p(z) &= p_s + g\int_{-h}^\eta \rho dz - g \int^z_{-h} \rho dz^\prime, \\ - &= p_s + g\rho_0\widetilde{H} - g \int^z_{-h} \rho dz^\prime, -$$ - -where the total water column pseudo height is expressed by - -$$ -\widetilde{H} = \int_{-h}^\eta \frac{\rho}{\rho_0} dz. -$$ - -$\widetilde{H}$ is the prognositc variable in the barotropic continuity equation. -The vertical integral of the pressure gradient is - -$$ -\frac{1}{\rho_0\widetilde{H}}\int^\eta_{-h} \nabla p dz &= \frac{1}{\rho_0\widetilde{H}}\int^\eta_{-h} \nabla \left( p_s + g\rho_0 \widetilde{H} - g \int^z_{-h} \rho dz^\prime \right) dz, \\ - &= \frac{H}{\rho_0\widetilde{H}}\nabla p_s + \frac{gH}{\widetilde{H}}\nabla \widetilde{H} - \frac{g}{\rho_0\widetilde{H}} \int_{-h}^\eta \left( \nabla \int_{-h}^z \rho dz^\prime\right) dz, -$$ - -where the height of the water column is represented by $H$. -The $1/\rho_0\widetilde{H}$ factor comes vertically integrating the material derivative and expressing the resulting barotropic momentum equation in non-conservative form. - -Therefore, the barotorpic pressure gradient term is discretized as: - -$$ -\overline{T}_e^p = g\left[ \frac{H_i}{\widetilde{H}_i} \right]_e\sum_{i \in CE(e)} n_{e,i}\widetilde{H}_e -$$ + T^p_{e,k} = \frac{1}{d_e}\sum_{i\in CE(e)} n_{e,i} \Phi_{i,k} + \frac{2}{d_e\sum_{i \in CE(e)}\tilde{h}_{i,k}}\left(\sum_{i \in CE(e)} n_{e,i} \tilde{h}_{i,k}\alpha_{i,k}p_{i,k}\right. \\ \left. - \frac{1}{2} \sum_{i\in CE(e)} \alpha_{i,k-1/2} p_{i,k-1/2}\sum_{i\in CE(e)} n_{e,i}\tilde{z}_{i,k-1/2} \right.\\ \left.+ \frac{1}{2} \sum_{i\in CE(e)} \alpha_{i,k+1/2} p_{i,k+1/2}\sum_{i\in CE(e)} n_{e,i}\tilde{z}_{i,k+1/2}\right), +$$ +where $k-1/2$ and $k+1/2$ refer to the top and bottom of layer $k$, respectively. + + +### 3.2 High-order Pressure Gradient +The high order pressure gradient will be based on the {ref}`full volume integral form ` of the geopotential and pressure terms: +$$ +T^p &= - \int_A \int_{\tilde{z}_k^{\text{bot}}}^{\tilde{z}_k^{\text{top}}} \rho_0 \, \left( \nabla \left<\Phi\right> \right) \, d\tilde{z} \, dA \\ +& - \int_{\partial A} \left( \int_{\tilde{z}_k^{\text{bot}}}^{\tilde{z}_k^{\text{top}}} \rho_0 \left(\left< \alpha \right> \left

+ \left<\alpha^\prime p^\prime\right> \right) \, d\tilde{z} \right) dl \\ +& - \int_A \rho_0 \left[ \left< \alpha \right> \left

+ \left<\alpha^\prime \left(p \nabla \tilde{z}_k^{\text{top}}\right)^\prime\right> \right]_{\tilde{z} = \tilde{z}_k^{\text{top}}} \, dA \\ +& + \int_A \rho_0 \left[ \left< \alpha \right> \left

+ \left<\alpha^\prime \left(p \nabla \tilde{z}_k^{\text{bot}}\right)^\prime\right> \right]_{\tilde{z} = \tilde{z}_k^{\text{bot}}} \, dA. +$$ +To obtain the expression that will be used, we neglect the turblent correlations and drop the Reynold's average notation for single variables: +$$ +T^p &= - \int_A \int_{\tilde{z}_k^{\text{bot}}}^{\tilde{z}_k^{\text{top}}} \rho_0 \, \left( \nabla \Phi \right) \, d\tilde{z} \, dA \\ +& - \int_{\partial A} \left( \int_{\tilde{z}_k^{\text{bot}}}^{\tilde{z}_k^{\text{top}}} \rho_0 \left(\alpha p \right) \, d\tilde{z} \right) dl \\ +& - \int_A \rho_0 \left[ \alpha \left

\right]_{\tilde{z} = \tilde{z}_k^{\text{top}}} \, dA \\ +& + \int_A \rho_0 \left[ \alpha \left

\right]_{\tilde{z} = \tilde{z}_k^{\text{bot}}} \, dA. +$$ +These volume and area integrals will be computed using quadrature to account for the variablity of $\alpha$ with the recontructed values of temperature, salinity, and pressure atthe quadrature points. +The complete details for the high-order pressure gradient will be the subject of a future design document. + +%### 3.3 Barotropic Pressure Gradient +% +%When split baroclinic-barotropic time stepping is implemented in the future, the barotropic pressure gradient will be calculated by the pressure gradient class. +%The barotropic pressure gradient is found by depth integrating the pressure gradient. +%The pressure is +% +%$$ +%p(z) = p_b - g \int^z_{-h} \rho dz^\prime, +%$$ +% +%where $p_b$ is the bottom pressure. +%The bottom pressure is the sum of the atmospheric surface pressure, $p_s$, and the pressure contribution of the water column: +% +%$$ +%p(z) &= p_s + g\int_{-h}^\eta \rho dz - g \int^z_{-h} \rho dz^\prime, \\ +% &= p_s + g\rho_0\widetilde{H} - g \int^z_{-h} \rho dz^\prime, +%$$ +% +%where the total water column pseudo height is expressed by +% +%$$ +%\widetilde{H} = \int_{-h}^\eta \frac{\rho}{\rho_0} dz. +%$$ +% +%$\widetilde{H}$ is the prognositc variable in the barotropic continuity equation. +%The vertical integral of the pressure gradient is +% +%$$ +%\frac{1}{\rho_0\widetilde{H}}\int^\eta_{-h} \nabla p dz &= \frac{1}{\rho_0\widetilde{H}}\int^\eta_{-h} \nabla \left( p_s + g\rho_0 \widetilde{H} - g \int^z_{-h} \rho dz^\prime \right) dz, \\ +% &= \frac{H}{\rho_0\widetilde{H}}\nabla p_s + \frac{gH}{\widetilde{H}}\nabla \widetilde{H} - \frac{g}{\rho_0\widetilde{H}} \int_{-h}^\eta \left( \nabla \int_{-h}^z \rho dz^\prime\right) dz, +%$$ +% +%where the height of the water column is represented by $H$. +%The $1/\rho_0\widetilde{H}$ factor comes vertically integrating the material derivative and expressing the resulting barotropic momentum equation in non-conservative form. +% +%Therefore, the barotorpic pressure gradient term is discretized as: +% +%$$ +%\overline{T}_e^p = g\left[ \frac{H_i}{\widetilde{H}_i} \right]_e\sum_{i \in CE(e)} n_{e,i}\widetilde{H}_e +%$$ ## 4 Design @@ -126,10 +143,11 @@ class PressureGrad{ private: std::unique_ptr OmegaPressureGrad; PressureGradCentered CenteredPGrad; - PressureGradHighOrder HighOrderPGrad; // To be implemented later + PressureGradHighOrder HighOrderPGrad1; // To be implemented later + PressureGradHighOrder HighOrderPGrad2; // Multiple high order options are likely in the future PressureGradType PressureGradChoice; I4 NVertLevels; - I4 NChuncks; + I4 NChunks; Array2DI4 CellsOnEdge; Array1DReal DvEdge; Array1DReal EdgeSignOnCell; @@ -144,7 +162,8 @@ An `enum class` will be used to specify options for the pressure gradient used f ```c++ enum class PressureGradType{ Centered, - HighOrder + HighOrder1, + HighOrder2 } ``` 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. @@ -164,12 +183,12 @@ class PressureGrad{ The constructor will be responsible for storing any static mesh information as private variables and handling the selection of the user-specified pressure gradient option. ```c++ -PressureGrad::PressureGrad(const HorzMesh *Mesh, int NVertLevels, Config *Options); +PressureGrad::PressureGrad(const HorzMesh *Mesh, const VertCoord *VCoord, Config *Options); ``` #### 4.2.2 Initialization -The init method will create the pressure gradient and return an error code: +The init method will create the default pressure gradient and return an error code: ```c++ int PressureGrad::init(); ``` @@ -193,17 +212,25 @@ void PressureGrad::computePressureGrad(const Array2DReal &Tend, const Array2DReal &SpecVol) { OMEGA_SCOPE(LocCenteredPGrad, CenteredPGrad) OMEGA_SCOPE(LocHighOrderPGrad, HighOrderPGrad) +const Array1DI4 &MinLyrEdgeBot = VCoord->MinLayerEdgeBot; +const Array1DI4 &MaxLyrEdgeTop = VCoord->MaxLayerEdgeTop; if (PressureGradChoice == PressureGradType::Centered){ - parallelFor("pgrad-centered", {NCellsAll, NChunks}, - KOKKOS_LAMBDA(int ICell, int KChunk) { - LocCenteredPGrad(Tend, Pressure, Geopotential, SpecVol); + parallelForOuter( + {NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int NChunks = computeNChunks(MinLyrEdgeBottom, MaxLyrEdgeTop, IEdge); + parallelForInner(Team, NChunks, [=](const int KChunk) { + LocCenteredPGrad(Tend, IEdge, KChunk, Pressure, Geopotential, SpecVol); + }); }); } else if (PressureGradChoice == PressureGradType::HighOrder){ - parallelFor("pgrad-highorder", {NCellsAll, NChunks}, - KOKKOS_LAMBDA(int ICell, int KChunk) { - LocHighOrderPGrad(Tend, Pressure, Geopotential, SpecVol); - }); + parallelForOuter( + {NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { + const int NChunks = computeNChunks(MinLyrEdgeBottom, MaxLyrEdgeTop, IEdge); + parallelForInner(Team, NChunks, [=](const int KChunk) { + LocHighOrderPGrad(Tend, IEdge, KChunk, Pressure, Geopotential, SpecVol); + }); + }); } ``` @@ -218,7 +245,7 @@ void PressureGrad::clear(); ## Verification and Testing ### Test: Spatial convergence to exact solution -For a given analytical $v$, $h$, and $b$, the spatial convergence of the pressure gradient can be assessed by computing errors on progressively finer meshes. +For given analytical functions of $\alpha$, $h$, and $z$, the spatial convergence of the pressure gradient can be assessed by computing errors on progressively finer meshes. ### Test: Baroclinic gyre The baroclinic gyre test case will test the pressure gradient term in the full non-Boussinesq equations. From ea64c55184575741f4354a9cac933d9f91c08bc6 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Fri, 10 Oct 2025 09:01:09 -0500 Subject: [PATCH 5/9] Fix sign in pressure gradient design document Co-authored-by: Mark Petersen --- components/omega/doc/design/PGrad.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/omega/doc/design/PGrad.md b/components/omega/doc/design/PGrad.md index 72fc1190b81a..b8afb485d455 100644 --- a/components/omega/doc/design/PGrad.md +++ b/components/omega/doc/design/PGrad.md @@ -46,7 +46,7 @@ The details of the barotropic pressure gradient will be added in a future design In the layered non-Boussinesq {ref}`momentum equation ` solved in Omega, the pressure gradient tendency term for edge $e$ and layer $k$, $T^p_{e,k}$, includes the gradient of the geopotential, the gradient of a term involving pressure, and two terms evaluated at the cell interface: $$ -T^p_{e,k} = - \left(\nabla \Phi \right)_{e,k} - \frac{1}{\left[\tilde{h}_k\right]_e} \nabla \left( \tilde{h}_k \alpha_k p_k \right) + \frac{1}{\left[\tilde{h}_k\right]_e} \left\{ \left[ \alpha p \nabla \tilde{z}\right]_{e,k}^\text{top} - \left[ \alpha p \nabla \tilde{z}\right]_{e,k+1}^\text{top} \right\}. +T^p_{e,k} = - \left(\nabla \Phi \right)_{e,k} - \frac{1}{\left[\tilde{h}_k\right]_e} \nabla \left( \tilde{h}_k \alpha_k p_k \right) - \frac{1}{\left[\tilde{h}_k\right]_e} \left\{ \left[ \alpha p \nabla \tilde{z}\right]_{e,k}^\text{top} - \left[ \alpha p \nabla \tilde{z}\right]_{e,k+1}^\text{top} \right\}. $$ The geopotential and interface terms are necessary to account for tilted layers that occur when using a general vertical coordinate, where the gradient operator is taken along layers. From 086bc304fd71108d4d516a8e6a0026cb740cabb8 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Fri, 5 Dec 2025 12:25:32 -0800 Subject: [PATCH 6/9] Update implementation details and testing --- components/omega/doc/design/PGrad.md | 199 ++++++++++++++++++++------- 1 file changed, 151 insertions(+), 48 deletions(-) diff --git a/components/omega/doc/design/PGrad.md b/components/omega/doc/design/PGrad.md index b8afb485d455..00c09e2d68bd 100644 --- a/components/omega/doc/design/PGrad.md +++ b/components/omega/doc/design/PGrad.md @@ -22,7 +22,7 @@ Note that we use the name "$p^\star$" to refer to the vertical coordinate, in Om ### 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 will be implemented and will replace the centered approach in production runs. +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 @@ -50,27 +50,38 @@ T^p_{e,k} = - \left(\nabla \Phi \right)_{e,k} - \frac{1}{\left[\tilde{h}_k\right $$ The geopotential and interface terms are necessary to account for tilted layers that occur when using a general vertical coordinate, where the gradient operator is taken along layers. -In this equation, $\alpha_{i,k}$ specific volume, $p_{i,k}$ is the pressure, and $\Phi_{i,k}$ is the geopotential. -These three quantities are evaluated at the mid-point of layer $k$ of cell $i$ in the first two terms and at the cell interfaces in the third term along with the interface psudo-height, \tilde{z}. +Even for layers at constant pseudo-height (rather than geometric height), the geoponential gradient would still be non-zero as will the second term (unless density is constant along layers). +In this equation, $\alpha_{i,k}$ is the specific volume, $p_{i,k}$ is the pressure, and $\Phi_{i,k}$ is the geopotential. +These three quantities are evaluated at the mid-point of layer $k$ of cell $i$ in the first two terms and at the cell interfaces in the third term along with the interface pseudo-height, $\tilde{z}$. The discrete gradient operator at an edge is: $$ - \nabla {(\cdot)} = \frac{1}{d_e} \sum_{i\in CE(e)} -n_{e,i} (\cdot)_i + \nabla {(\cdot)} = \frac{1}{d_e} \sum_{i\in CE(e)} -n_{e,i} (\cdot)_i, $$ where $d_e$ is the distance between cell centers, $CE(e)$ are the cells on edge $e$, and $n_{e,i}$ is the sign of the edge normal with respect to cell $i$. The (cell-to-edge) horizontal averaging operator is: $$ - [\cdot]_e = \frac{1}{2}\sum_{i\in CE(e)} (\cdot)_i + [\cdot]_e = \frac{1}{2}\sum_{i\in CE(e)} (\cdot)_i. +$$ + +The $\gamma_{e,k-1/2} \equiv \left[ \alpha p \nabla \tilde{z}\right]_{e,k}^\text{top}/\left[\tilde{h}_k\right]_e$ term is computed using a four point average for all interior layer interfaces with a two-point averaged used for the top and bottom interfaces. +The interior interfaces are calculated as: +$$ +\gamma_{e,k-1/2} = \frac{\sum_{i\in CE(e)} \alpha_{i,k} p_{i,k} h_{i,k} + \alpha_{i,k-1} p_{i,k-1} h_{i,k-1}}{d_e\sum_{i\in CE(e)} h_{i,k} + h_{i,k-1}} \sum_{i\in CE(e)} n_{e,i}\tilde{z}_{i,k-1/2} +$$ +with the top and bottom computed as: +$$ +\gamma_{e,1/2} = \frac{\sum_{i\in CE(e)} \alpha_{i,k} p_{i,1} h_{i,1}}{d_e\sum_{i\in CE(e)} h_{i,1}} \sum_{i\in CE(e)} n_{e,i}\tilde{z}_{i,1/2}, \\ +\gamma_{e,K+1/2} = \frac{\sum_{i\in CE(e)} \alpha_{i,K} p_{i,k} h_{i,K}}{d_e\sum_{i\in CE(e)} h_{i,K}} \sum_{i\in CE(e)} n_{e,i}\tilde{z}_{i,K+1/2}. $$ Therefore, the centered pressure gradient will be calculated as: $$ - T^p_{e,k} = \frac{1}{d_e}\sum_{i\in CE(e)} n_{e,i} \Phi_{i,k} + \frac{2}{d_e\sum_{i \in CE(e)}\tilde{h}_{i,k}}\left(\sum_{i \in CE(e)} n_{e,i} \tilde{h}_{i,k}\alpha_{i,k}p_{i,k}\right. \\ \left. - \frac{1}{2} \sum_{i\in CE(e)} \alpha_{i,k-1/2} p_{i,k-1/2}\sum_{i\in CE(e)} n_{e,i}\tilde{z}_{i,k-1/2} \right.\\ \left.+ \frac{1}{2} \sum_{i\in CE(e)} \alpha_{i,k+1/2} p_{i,k+1/2}\sum_{i\in CE(e)} n_{e,i}\tilde{z}_{i,k+1/2}\right), + T^p_{e,k} = \frac{1}{d_e}\sum_{i\in CE(e)} n_{e,i} \Phi_{i,k} + \frac{2}{d_e\sum_{i \in CE(e)}\tilde{h}_{i,k}}\sum_{i \in CE(e)} n_{e,i} \tilde{h}_{i,k}\alpha_{i,k}p_{i,k} - \gamma_{e,k-1/2} + \gamma_{e,k+1/2}. $$ -where $k-1/2$ and $k+1/2$ refer to the top and bottom of layer $k$, respectively. ### 3.2 High-order Pressure Gradient @@ -81,14 +92,14 @@ T^p &= - \int_A \int_{\tilde{z}_k^{\text{bot}}}^{\tilde{z}_k^{\text{top}}} \rho_ & - \int_A \rho_0 \left[ \left< \alpha \right> \left

+ \left<\alpha^\prime \left(p \nabla \tilde{z}_k^{\text{top}}\right)^\prime\right> \right]_{\tilde{z} = \tilde{z}_k^{\text{top}}} \, dA \\ & + \int_A \rho_0 \left[ \left< \alpha \right> \left

+ \left<\alpha^\prime \left(p \nabla \tilde{z}_k^{\text{bot}}\right)^\prime\right> \right]_{\tilde{z} = \tilde{z}_k^{\text{bot}}} \, dA. $$ -To obtain the expression that will be used, we neglect the turblent correlations and drop the Reynold's average notation for single variables: +To obtain the expression that will be used, we neglect the turbulent correlations and drop the Reynold's average notation for single variables: $$ T^p &= - \int_A \int_{\tilde{z}_k^{\text{bot}}}^{\tilde{z}_k^{\text{top}}} \rho_0 \, \left( \nabla \Phi \right) \, d\tilde{z} \, dA \\ & - \int_{\partial A} \left( \int_{\tilde{z}_k^{\text{bot}}}^{\tilde{z}_k^{\text{top}}} \rho_0 \left(\alpha p \right) \, d\tilde{z} \right) dl \\ & - \int_A \rho_0 \left[ \alpha \left

\right]_{\tilde{z} = \tilde{z}_k^{\text{top}}} \, dA \\ & + \int_A \rho_0 \left[ \alpha \left

\right]_{\tilde{z} = \tilde{z}_k^{\text{bot}}} \, dA. $$ -These volume and area integrals will be computed using quadrature to account for the variablity of $\alpha$ with the recontructed values of temperature, salinity, and pressure atthe quadrature points. +These volume and area integrals will be computed using quadrature to account for the variability of $\alpha$ with the reconstructed values of temperature, salinity, and pressure at the quadrature points. The complete details for the high-order pressure gradient will be the subject of a future design document. %### 3.3 Barotropic Pressure Gradient @@ -141,18 +152,31 @@ The `PressureGradient` class will be used to perform the horizontal gradients of class PressureGrad{ public: private: - std::unique_ptr OmegaPressureGrad; + // Map of all pressure gradients + static std::map> 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; - I4 NVertLevels; - I4 NChunks; + + // Data required for computation (stored copies of HorzMesh/VCoord arrays) Array2DI4 CellsOnEdge; Array1DReal DvEdge; Array1DReal EdgeSignOnCell; + Array1DI4 MinLayerEdgeBot; + Array1DI4 MaxLaterEdgeTop; + + // Array for interface produce needed in centered pressure gradient + Array2DReal InterfaceProduct; } ``` +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: ```yaml PressureGrad: @@ -166,39 +190,67 @@ enum class PressureGradType{ HighOrder2 } ``` -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. ### 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. ```c++ class PressureGrad{ public: + // Creation methods static PressureGrad *init(); - static PressureGrad *get(); - void computePressureGrad(); + 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); private: + // Helper compute method + void computeInterfaceProduct(const Array2DReal &PressureMid, + const Array2DReal &SpecVol, + const Array2DReal &LayerThick, + const Array2DReal &ZInterface); } ``` #### 4.2.1 Creation -The constructor will be responsible for storing any static mesh information as private variables and handling the selection of the user-specified pressure gradient option. +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. ```c++ -PressureGrad::PressureGrad(const HorzMesh *Mesh, const VertCoord *VCoord, Config *Options); +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: ```c++ -int PressureGrad::init(); +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 pressure gradient instance: +There will be a method for getting a pointer to the default pressure gradient instance: ```c++ -PressureGrad *PressureGrad::get(); +PressureGrad *PressureGrad::getDefault(); +``` +or a named instance: +```c++ +PressureGrad *PressureGrad::get(const std::string &Name); ``` #### 4.2.4 Computation @@ -207,45 +259,96 @@ The public `computePressureGrad` method will rely on private methods for each sp 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. ```c++ void PressureGrad::computePressureGrad(const Array2DReal &Tend, - const Array2DReal &Pressure, - const Array2DReal &Geopotential, - const Array2DReal &SpecVol) { -OMEGA_SCOPE(LocCenteredPGrad, CenteredPGrad) -OMEGA_SCOPE(LocHighOrderPGrad, HighOrderPGrad) -const Array1DI4 &MinLyrEdgeBot = VCoord->MinLayerEdgeBot; -const Array1DI4 &MaxLyrEdgeTop = VCoord->MaxLayerEdgeTop; - if (PressureGradChoice == PressureGradType::Centered){ - parallelForOuter( - {NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { - const int NChunks = computeNChunks(MinLyrEdgeBottom, MaxLyrEdgeTop, IEdge); - parallelForInner(Team, NChunks, [=](const int KChunk) { - LocCenteredPGrad(Tend, IEdge, KChunk, Pressure, Geopotential, SpecVol); - }); - }); - } - else if (PressureGradChoice == PressureGradType::HighOrder){ - parallelForOuter( - {NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) { - const int NChunks = computeNChunks(MinLyrEdgeBottom, MaxLyrEdgeTop, IEdge); - parallelForInner(Team, NChunks, [=](const int KChunk) { - LocHighOrderPGrad(Tend, IEdge, KChunk, Pressure, Geopotential, SpecVol); - }); - }); - } + 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 alpha*p*grad(z) term at edge interfaces + computeInterfaceProduct(PressureMid, SpecVol, + LayerThick, ZInterface); + + // 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, + Geopotential, LayerThick, + LocInterfaceProduct, 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 clear method will remove the pressure gradient instance, calling the destructor in the process. +The there will be a method to remove all pressure gradient instances, ```c++ -void PressureGrad::clear(); +static void PressureGrad::clear(); ``` +or a named instance +```c++ +static void erase(const std::string &Name); +``` + ## Verification and Testing +### 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 $\alpha$, $h$, and $z$, the spatial convergence of the pressure gradient can be assessed by computing errors on progressively finer meshes. +This will be compared to an exact solution computed using a high-order quadrature. + +### 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. From 16fdee7ee741334774e88f5b5f3dd5edb27359df Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Fri, 20 Feb 2026 13:24:04 -0800 Subject: [PATCH 7/9] Revise for grad(montgomery potential) - p grad(alpha) form --- components/omega/doc/design/PGrad.md | 89 +++++++++++++--------------- 1 file changed, 40 insertions(+), 49 deletions(-) diff --git a/components/omega/doc/design/PGrad.md b/components/omega/doc/design/PGrad.md index 00c09e2d68bd..886a77abc0e2 100644 --- a/components/omega/doc/design/PGrad.md +++ b/components/omega/doc/design/PGrad.md @@ -6,8 +6,8 @@ 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](https://journals.ametsoc.org/view/journals/phoc/32/7/1520-0485_2002_032_2194_tdbtba_2.0.co_2.xml). - However, Omega is written in general vertical coordinates and can reference either pressure $p$ or distance $z$ in the vertical. +Some non-Boussinesq models are written in pressure coordinates (e.g. [de Szoeke and Samelson 2002](https://journals.ametsoc.org/view/journals/phoc/32/7/1520-0485_2002_032_2194_tdbtba_2.0.co_2.xml). +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. @@ -42,48 +42,53 @@ For split barotropic-baroclinic timestepping, the pressure gradient should provi The details of the barotropic pressure gradient will be added in a future design document for split time stepping. ## 3 Algorithmic Formulation -### 3.1 Centered Pressure Gradient -In the layered non-Boussinesq {ref}`momentum equation ` solved in Omega, the pressure gradient tendency term for edge $e$ and layer $k$, $T^p_{e,k}$, includes the gradient of the geopotential, the gradient of a term involving pressure, and two terms evaluated at the cell interface: - +In the layered non-Boussinesq equation solved in Omega, the pressure and geopotential gradient terms are taken along a generalized coordinate $r$: +$$ +\alpha \nabla_z p + \nabla_z \Phi = \alpha\left( \nabla_r p - \frac{\partial r}{\partial z}\frac{\partial p}{\partial r}\nabla_r z\right) + \nabla_r \Phi - \frac{\partial r}{\partial z}\frac{\partial \Phi}{\partial r} \nabla_r z. +$$ +Using the hydrostatic relation +$$ +\frac{\partial p}{\partial z} &= - \rho g, \\ +\frac{\partial p}{\partial r}\frac{\partial r}{\partial z} &= - \rho g, \\ +\frac{\partial r}{\partial z} &= -\rho g \frac{\partial r}{\partial p}. +$$ +to subtitute for $\partial r/\partial z$, and substituting $\Phi = gz + \phi_{TP} + \phi_{SAL}$ yields: +$$ +\alpha \nabla_z p &= \alpha\left( \nabla_r p + \rho g \nabla_r z\right) + \nabla_r \Phi - \frac{\partial r}{\partial z}\frac{\partial gz}{\partial r} \nabla_r z, \\ +&= \alpha\nabla_r p + g \nabla_r z + \nabla_r \Phi - g \nabla_r z, \\ +&= \nabla_r \alpha p - p\nabla_r \alpha + g\nabla_r z + \nabla_r\left(\Phi - gz\right)\\ +&= \nabla_r\left(\alpha p + gz \right) - p\nabla_r\alpha + \nabla_r\left(\phi_{TP} + \phi_{SAL}\right). $$ -T^p_{e,k} = - \left(\nabla \Phi \right)_{e,k} - \frac{1}{\left[\tilde{h}_k\right]_e} \nabla \left( \tilde{h}_k \alpha_k p_k \right) - \frac{1}{\left[\tilde{h}_k\right]_e} \left\{ \left[ \alpha p \nabla \tilde{z}\right]_{e,k}^\text{top} - \left[ \alpha p \nabla \tilde{z}\right]_{e,k+1}^\text{top} \right\}. +The Montgomery potential is defined as: +$$ +M = \alpha p + gz. $$ -The geopotential and interface terms are necessary to account for tilted layers that occur when using a general vertical coordinate, where the gradient operator is taken along layers. -Even for layers at constant pseudo-height (rather than geometric height), the geoponential gradient would still be non-zero as will the second term (unless density is constant along layers). -In this equation, $\alpha_{i,k}$ is the specific volume, $p_{i,k}$ is the pressure, and $\Phi_{i,k}$ is the geopotential. -These three quantities are evaluated at the mid-point of layer $k$ of cell $i$ in the first two terms and at the cell interfaces in the third term along with the interface pseudo-height, $\tilde{z}$. +### 3.1 Centered Pressure Gradient +The tendency term for edge $e$ and layer $k$, $T^p_{e,k}$, is discretized as: +$$ +T^p_{e,k} &= \frac{1}{2}\left[\nabla M_{k+1/2} + \nabla M_{k-1/2}\right] - \left[p_k\right]_e \nabla\alpha_k + \nabla(\phi_{TP} + \phi_{SAL}), \\ +$$ +where the gradient of the Montgomery potential is evaluated at layer interfaces and then averaged to the cell mid point. +In caclulating $M$ at the layer interfaces, $p$ and $z$ are evaluated at their natural interface locations, and $\alpha$, is held constant throughout the cell, i.e., +$$ +M_{k+1/2} &= \alpha_k p_{k+1/2} + gz_{k+1/2}\\ +M_{k-1/2} &= \alpha_k p_{k-1/2} + gz_{k-1/2} +$$ The discrete gradient operator at an edge is: - $$ \nabla {(\cdot)} = \frac{1}{d_e} \sum_{i\in CE(e)} -n_{e,i} (\cdot)_i, $$ - where $d_e$ is the distance between cell centers, $CE(e)$ are the cells on edge $e$, and $n_{e,i}$ is the sign of the edge normal with respect to cell $i$. The (cell-to-edge) horizontal averaging operator is: - $$ [\cdot]_e = \frac{1}{2}\sum_{i\in CE(e)} (\cdot)_i. $$ - -The $\gamma_{e,k-1/2} \equiv \left[ \alpha p \nabla \tilde{z}\right]_{e,k}^\text{top}/\left[\tilde{h}_k\right]_e$ term is computed using a four point average for all interior layer interfaces with a two-point averaged used for the top and bottom interfaces. -The interior interfaces are calculated as: -$$ -\gamma_{e,k-1/2} = \frac{\sum_{i\in CE(e)} \alpha_{i,k} p_{i,k} h_{i,k} + \alpha_{i,k-1} p_{i,k-1} h_{i,k-1}}{d_e\sum_{i\in CE(e)} h_{i,k} + h_{i,k-1}} \sum_{i\in CE(e)} n_{e,i}\tilde{z}_{i,k-1/2} -$$ -with the top and bottom computed as: -$$ -\gamma_{e,1/2} = \frac{\sum_{i\in CE(e)} \alpha_{i,k} p_{i,1} h_{i,1}}{d_e\sum_{i\in CE(e)} h_{i,1}} \sum_{i\in CE(e)} n_{e,i}\tilde{z}_{i,1/2}, \\ -\gamma_{e,K+1/2} = \frac{\sum_{i\in CE(e)} \alpha_{i,K} p_{i,k} h_{i,K}}{d_e\sum_{i\in CE(e)} h_{i,K}} \sum_{i\in CE(e)} n_{e,i}\tilde{z}_{i,K+1/2}. -$$ - Therefore, the centered pressure gradient will be calculated as: - $$ - T^p_{e,k} = \frac{1}{d_e}\sum_{i\in CE(e)} n_{e,i} \Phi_{i,k} + \frac{2}{d_e\sum_{i \in CE(e)}\tilde{h}_{i,k}}\sum_{i \in CE(e)} n_{e,i} \tilde{h}_{i,k}\alpha_{i,k}p_{i,k} - \gamma_{e,k-1/2} + \gamma_{e,k+1/2}. + T^p_{e,k} = \frac{1}{2d_e}\left(\sum_{i\in CE(e)} -n_{e,i} M_{i,k+1/2} + \sum_{i\in CE(e)} -n_{e,i} M_{i,k-1/2}\right) + \frac{1}{2}\left(\sum_{i\in CE(e)} p_{i,k}\right) \frac{1}{d_e} \sum_{i\in CE(e)} -n_{e,i}\alpha_{i,k} + \frac{1}{d_e} \sum_{i\in CE(e)} -n_{e,i} (\phi_{TP,i} + \phi_{SAL,i}). $$ - ### 3.2 High-order Pressure Gradient The high order pressure gradient will be based on the {ref}`full volume integral form ` of the geopotential and pressure terms: $$ @@ -165,13 +170,10 @@ class PressureGrad{ // Data required for computation (stored copies of HorzMesh/VCoord arrays) Array2DI4 CellsOnEdge; - Array1DReal DvEdge; + Array1DReal DcEdge; Array1DReal EdgeSignOnCell; Array1DI4 MinLayerEdgeBot; Array1DI4 MaxLaterEdgeTop; - - // Array for interface produce needed in centered pressure gradient - Array2DReal InterfaceProduct; } ``` 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. @@ -214,14 +216,7 @@ class PressureGrad{ // Main compute method void computePressureGrad(Array2DReal &Tend, const OceanState *State, const VertCoord *VCoord, const Eos *EqState, - const int TimeLevel); - private: - // Helper compute method - void computeInterfaceProduct(const Array2DReal &PressureMid, - const Array2DReal &SpecVol, - const Array2DReal &LayerThick, - const Array2DReal &ZInterface); - + const int TimeLevel) const; } ``` #### 4.2.1 Creation @@ -279,10 +274,6 @@ void PressureGrad::computePressureGrad(const Array2DReal &Tend, if (PressureGradChoice == PressureGradType::Centered) { - // computes alpha*p*grad(z) term at edge interfaces - computeInterfaceProduct(PressureMid, SpecVol, - LayerThick, ZInterface); - // computes centered geopotential and pressure gradient tendency parallelForOuter( "pgrad-centered", {NEdgesAll}, @@ -295,8 +286,9 @@ void PressureGrad::computePressureGrad(const Array2DReal &Tend, Team, KRange, INNER_LAMBDA(int KChunk) { LocCenteredPGrad(Tend, IEdge, KChunk, PressureMid, - Geopotential, LayerThick, - LocInterfaceProduct, SpecVol); + PressureInterface, ZInterface, + LocTidalPotential, LocSelfAttractionLoading, + SpecVol); }); }); @@ -337,15 +329,14 @@ static void erase(const std::string &Name); ## Verification and Testing -### Test: Hydrostatic test with tilted layers +### 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 $\alpha$, $h$, and $z$, the spatial convergence of the pressure gradient can be assessed by computing errors on progressively finer meshes. -This will be compared to an exact solution computed using a high-order quadrature. +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. From c52994dd1bf84de4d6f5de2e596902436b24adec Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Tue, 17 Mar 2026 07:44:11 -0700 Subject: [PATCH 8/9] Fix linting issues --- components/omega/doc/design/PGrad.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/components/omega/doc/design/PGrad.md b/components/omega/doc/design/PGrad.md index 886a77abc0e2..fba21d2c2a88 100644 --- a/components/omega/doc/design/PGrad.md +++ b/components/omega/doc/design/PGrad.md @@ -57,7 +57,7 @@ $$ \alpha \nabla_z p &= \alpha\left( \nabla_r p + \rho g \nabla_r z\right) + \nabla_r \Phi - \frac{\partial r}{\partial z}\frac{\partial gz}{\partial r} \nabla_r z, \\ &= \alpha\nabla_r p + g \nabla_r z + \nabla_r \Phi - g \nabla_r z, \\ &= \nabla_r \alpha p - p\nabla_r \alpha + g\nabla_r z + \nabla_r\left(\Phi - gz\right)\\ -&= \nabla_r\left(\alpha p + gz \right) - p\nabla_r\alpha + \nabla_r\left(\phi_{TP} + \phi_{SAL}\right). +&= \nabla_r\left(\alpha p + gz \right) - p\nabla_r\alpha + \nabla_r\left(\phi_{TP} + \phi_{SAL}\right). $$ The Montgomery potential is defined as: $$ @@ -223,7 +223,7 @@ class PressureGrad{ 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. ```c++ -PressureGrad *PressureGrad::create(const std::string &Name, +PressureGrad *PressureGrad::create(const std::string &Name, const HorzMesh *Mesh, const VertCoord *VCoord, Config *Options); @@ -293,7 +293,7 @@ void PressureGrad::computePressureGrad(const Array2DReal &Tend, }); } else { - + // computes high-order geopotential and pressure gradient tendency parallelForOuter( "pgrad-highorder", {NEdgesAll}, @@ -331,7 +331,7 @@ static void erase(const std::string &Name); ### 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 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. From b82a4ae8755550d5061a7a3c6c46a52092c3d689 Mon Sep 17 00:00:00 2001 From: Steven Brus Date: Fri, 20 Mar 2026 14:42:34 -0700 Subject: [PATCH 9/9] Switch baroclinic gyre->overflow and add urls for polaris tests --- components/omega/doc/design/PGrad.md | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/components/omega/doc/design/PGrad.md b/components/omega/doc/design/PGrad.md index fba21d2c2a88..b263bb2e37a4 100644 --- a/components/omega/doc/design/PGrad.md +++ b/components/omega/doc/design/PGrad.md @@ -337,9 +337,10 @@ The pressure gradient term will be compared to the zero exact solution at progre ### 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. +This will be implemented as a two-column test in Polaris. ### 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. +The [seamount test](https://docs.e3sm.org/polaris/main/users_guide/ocean/tasks/seamount.html) 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. +### Test: Overflow test +The [overflow test](https://docs.e3sm.org/polaris/main/users_guide/ocean/tasks/overflow.html) in Polaris will test the pressure gradient term in the full non-Boussinesq equations.