Issue538 seasurfaceheight (#662)

# Add sea-surface height Fixes #538 ### Task List - [x] Defined the tests that specify a complete and functioning change (*It may help to create a [design specification & test specification](../../../wiki/Specification-Template)*) - [x] Implemented the source code change that satisfies the tests - [ ] Documented the feature by providing worked example - [ ] Updated the README or other documentation - [x] Completed the pre-Request checklist below --- # Change Description This PR adds the contribution of the sea-surface-height gradient to both the mEVP and BBM codes (through VPCGDynamicsKernel and BrittleDynamicsKernel). SSH is read in from IOceanBoundary and passed to the dynamical core. It can also be output through ConfigOutput as "ssh". --- # Test Description No unit test was implemented for this feature, and no quantitative test exists. Qualitative testing shows that setting a fixed sea-surface-height gradient gives drift speed close to back-of-the-envelope estimates (13 cm/s modelled, 20 cm/s from back-of-the-envelope). --- # Documentation Impact None --- # Other Details None --- ### Pre-Request Checklist - [x] The requirements of this pull request are fully captured in an issue or design specification and are linked and summarised in the description of this PR - [x] No new warnings are generated - [x] The documentation has been updated (or an issue has been created to track the corresponding change) - [x] Methods and Tests are commented such that they can be understood without having to obtain additional context - [x] This PR/Issue is labelled as a bug/feature/enhancement/breaking change - [x] File dates have been updated to reflect modification date - [x] This change conforms to the conventions described in the README
nextsimhub · Oct 8, 2024 · ad86763 · ad86763
2 parents 334846d + 8115736
commit ad86763
Show file tree

Hide file tree

Showing 24 changed files with 340 additions and 51 deletions.
diff --git a/core/src/include/ModelComponent.hpp b/core/src/include/ModelComponent.hpp
@@ -50,6 +50,7 @@ namespace Protected {
         = "EXT_SSS"; // sea surface salinity from coupling or forcing, PSU
     inline constexpr TextTag EXT_SST
         = "EXT_SST"; // sea surface temperature from coupling or forcing, ˚C
+    inline constexpr TextTag SSH = "SSH"; // sea surface height, m
     inline constexpr TextTag EVAP_MINUS_PRECIP
         = "E-P"; // E-P atmospheric freshwater flux, kg s⁻¹ m⁻²
     // Derived fields, calculated once per timestep

diff --git a/core/src/include/gridNames.hpp b/core/src/include/gridNames.hpp
@@ -29,6 +29,7 @@ static const std::string uWindName = "uwind";
 static const std::string vWindName = "vwind";
 static const std::string uOceanName = "uocean";
 static const std::string vOceanName = "vocean";
+static const std::string sshName = "ssh";
 
 static const std::string uIOStressName = "uiostress";
 static const std::string vIOStressName = "viostress";

diff --git a/core/src/modules/DynamicsModule/BBMDynamics.cpp b/core/src/modules/DynamicsModule/BBMDynamics.cpp
@@ -80,6 +80,7 @@ void BBMDynamics::update(const TimestepTime& tst)
     kernel.setData(vWindName, vwind.data());
     kernel.setData(uOceanName, uocean.data());
     kernel.setData(vOceanName, vocean.data());
+    kernel.setData(sshName, ssh.data());
 
     /*
      * Ice velocity components are stored in the dynamics, and not changed by the model outside the

diff --git a/core/src/modules/DynamicsModule/MEVPDynamics.cpp b/core/src/modules/DynamicsModule/MEVPDynamics.cpp
@@ -69,6 +69,7 @@ void MEVPDynamics::update(const TimestepTime& tst)
     kernel.setData(vWindName, vwind.data());
     kernel.setData(uOceanName, uocean.data());
     kernel.setData(vOceanName, vocean.data());
+    kernel.setData(sshName, ssh.data());
 
     // kernel.setData(uName, uice);
     // kernel.setData(vName, vice);

diff --git a/core/src/modules/include/IDynamics.hpp b/core/src/modules/include/IDynamics.hpp
@@ -34,6 +34,7 @@ class IDynamics : public ModelComponent {
         , vwind(getStore())
         , uocean(getStore())
         , vocean(getStore())
+        , ssh(getStore())
         , m_usesDamage(usesDamageIn)
     {
         getStore().registerArray(Shared::DAMAGE, &damage, RW);
@@ -95,6 +96,7 @@ class IDynamics : public ModelComponent {
     ModelArrayRef<Protected::WIND_V> vwind;
     ModelArrayRef<Protected::OCEAN_U> uocean;
     ModelArrayRef<Protected::OCEAN_V> vocean;
+    ModelArrayRef<Protected::SSH> ssh;
 
     // Does this implementation of the dynamics use damage?
     bool m_usesDamage;

diff --git a/core/src/modules/include/ProtectedArrayNames.ipp b/core/src/modules/include/ProtectedArrayNames.ipp
@@ -1,44 +1,45 @@
 /*!
  * @file ProtectedArrayNames.ipp
  *
- * @date 1 Jul 2024
+ * @date 23 Aug 2024
  * @author Tim Spain <[email protected]>
  * @author Einar Ólason <[email protected]>
  */
 
 // External ProtectedArray names must be unique wrt to the external SharedArray names as well
 { "damage", "DAMAGE0" }, // Ice thickness, cell average, m
-    { "hice", "H_ICE_cell" }, // Ice thickness, cell average, m
-    { "cice", "C_ICE0" }, // Ice concentration
-    { "hsnow", "H_SNOW_cell" }, // Snow depth, cell average, m
-    { "tice", "T_ICE0" }, // Ice temperature, ˚C
-    { "tair", "T_AIR" }, // Air temperature, ˚C
-    { "dew2m", "DEW_2M" }, // Dew point at 2 m, ˚C
-    { "pair", "P_AIR" }, // sea level air pressure, Pa
-    { "mixrat", "MIXRAT" }, // water vapour mass mixing ratio
-    { "sw_in", "SW_IN" }, // incoming shortwave flux, W m⁻²
-    { "lw_in", "LW_IN" }, // incoming longwave flux, W m⁻²
-    { "mld", "MLD" }, // mixed layer depth, m
-    { "snowfall", "SNOWFALL" }, // snow fall, kg m⁻² s⁻¹
-    { "sss", "SSS" }, // sea surface salinity, PSU
-    { "sst", "SST" }, // sea surface temperature ˚C
-    { "sst_ext", "EXT_SST" }, // External sea surface temperature ˚C
-    { "sss_ext", "EXT_SSS" }, // External sea surface salinity PSU
-    { "eminusp", "E-P" }, // E-P atmospheric freshwater flux, kg s⁻¹ m⁻²
-    { "mlcp", "CPML" }, // Mixed layer bulk heat capacity J K⁻¹ m⁻²
-    { "tf", "TF" }, // Ocean freezing temperature, ˚C
-    { "wind_speed", "WIND_SPEED" }, // Wind speed, m s⁻¹
-    { "wind_u", "WIND_U" }, // wind velocity x component, m s⁻¹
-    { "wind_v", "WIND_V" }, // wind velocity y component, m s⁻¹
-    { "hice_true_pro", "HTRUE_ICE" }, // Ice thickness, ice average, m
-    { "hsnow_true_pro", "HTRUE_SNOW" }, // Snow thickness, ice average, m
-    { "ocean_u", "OCEAN_U" }, // x(east)-ward ocean current, m s⁻¹
-    { "ocean_v", "OCEAN_V" }, // y(north)-ward ocean current, m s⁻¹
-    { "u", "ICE_U" }, // x(east)-ward ice velocity, m s⁻¹
-    { "v", "ICE_V" }, // y(north)-ward ice velocity, m s⁻¹
-    { "sst_slab", "SLAB_SST" }, // Slab ocean surface temperature ˚C
-    { "sss_slab", "SLAB_SSS" }, // Slab ocean surface salinity PSU
-    { "qdw", "SLAB_QDW" }, // Slab ocean temperature nudging heat flux, W m⁻²
-    { "fdw", "SLAB_FDW" }, // Slab ocean salinity nudging water flux, kg s⁻¹ m⁻²
+{ "hice", "H_ICE_cell" }, // Ice thickness, cell average, m
+{ "cice", "C_ICE0" }, // Ice concentration
+{ "hsnow", "H_SNOW_cell" }, // Snow depth, cell average, m
+{ "tice", "T_ICE0" }, // Ice temperature, ˚C
+{ "tair", "T_AIR" }, // Air temperature, ˚C
+{ "dew2m", "DEW_2M" }, // Dew point at 2 m, ˚C
+{ "pair", "P_AIR" }, // sea level air pressure, Pa
+{ "mixrat", "MIXRAT" }, // water vapour mass mixing ratio
+{ "sw_in", "SW_IN" }, // incoming shortwave flux, W m⁻²
+{ "lw_in", "LW_IN" }, // incoming longwave flux, W m⁻²
+{ "mld", "MLD" }, // mixed layer depth, m
+{ "snowfall", "SNOWFALL" }, // snow fall, kg m⁻² s⁻¹
+{ "sss", "SSS" }, // sea surface salinity, PSU
+{ "sst", "SST" }, // sea surface temperature ˚C
+{ "sst_ext", "EXT_SST" }, // External sea surface temperature ˚C
+{ "sss_ext", "EXT_SSS" }, // External sea surface salinity PSU
+{ "eminusp", "E-P" }, // E-P atmospheric freshwater flux, kg s⁻¹ m⁻²
+{ "mlcp", "CPML" }, // Mixed layer bulk heat capacity J K⁻¹ m⁻²
+{ "tf", "TF" }, // Ocean freezing temperature, ˚C
+{ "wind_speed", "WIND_SPEED" }, // Wind speed, m s⁻¹
+{ "wind_u", "WIND_U" }, // wind velocity x component, m s⁻¹
+{ "wind_v", "WIND_V" }, // wind velocity y component, m s⁻¹
+{ "hice_true_pro", "HTRUE_ICE" }, // Ice thickness, ice average, m
+{ "hsnow_true_pro", "HTRUE_SNOW" }, // Snow thickness, ice average, m
+{ "ocean_u", "OCEAN_U" }, // x(east)-ward ocean current, m s⁻¹
+{ "ocean_v", "OCEAN_V" }, // y(north)-ward ocean current, m s⁻¹
+{ "u", "ICE_U" }, // x(east)-ward ice velocity, m s⁻¹
+{ "v", "ICE_V" }, // y(north)-ward ice velocity, m s⁻¹
+{ "sst_slab", "SLAB_SST" }, // Slab ocean surface temperature ˚C
+{ "sss_slab", "SLAB_SSS" }, // Slab ocean surface salinity PSU
+{ "qdw", "SLAB_QDW" }, // Slab ocean temperature nudging heat flux, W m⁻²
+{ "fdw", "SLAB_FDW" }, // Slab ocean salinity nudging water flux, kg s⁻¹ m⁻²
+{ "ssh", "SSH" }, // Slab ocean salinity nudging water flux, kg s⁻¹ m⁻²
 { "taux", "IO_STRESS_U" }, // Ice-ocean stress x(east) direction, Pa
 { "tauy", "IO_STRESS_V" }, // Ice-ocean stress x(east) direction, Pa
diff --git a/dynamics/src/CGDynamicsKernel.cpp b/dynamics/src/CGDynamicsKernel.cpp
@@ -37,6 +37,9 @@ void CGDynamicsKernel<DGadvection>::initialise(
     cgH.resize_by_mesh(*smesh);
     cgA.resize_by_mesh(*smesh);
 
+    uGradSeasurfaceHeight.resize_by_mesh(*smesh);
+    vGradSeasurfaceHeight.resize_by_mesh(*smesh);
+
     dStressX.resize_by_mesh(*smesh);
     dStressY.resize_by_mesh(*smesh);
 
@@ -115,6 +118,141 @@ template <int DGadvection> void CGDynamicsKernel<DGadvection>::prepareAdvection(
     dgtransport->prepareAdvection(u, v);
 }
 
+template <int DGadvection>
+void CGDynamicsKernel<DGadvection>::ComputeGradientOfSeaSurfaceHeight(
+    const DGVector<1>& seaSurfaceHeight)
+{
+    // First transform to CG1 Vector and do all computations in CG1
+    CGVector<1> cgSeasurfaceHeight(*smesh);
+    Interpolations::DG2CG(*smesh, cgSeasurfaceHeight, seaSurfaceHeight);
+
+    CGVector<1> uGrad(*smesh);
+    CGVector<1> vGrad(*smesh);
+    uGrad.setZero();
+    vGrad.setZero();
+
+    // parallelization in stripes
+    for (size_t p = 0; p < 2; ++p) {
+#pragma omp parallel for schedule(static)
+        for (size_t cy = 0; cy < smesh->ny; ++cy) {
+            //!< loop over all cells of the mesh
+            if (cy % 2 == p) {
+                size_t eid = smesh->nx * cy;
+                size_t cg1id = cy * (smesh->nx + 1);
+                for (size_t cx = 0; cx < smesh->nx; ++cx, ++eid, ++cg1id) {
+                    // get local CG nodes
+                    Eigen::Vector<Nextsim::FloatType, 4> loc_cgSSH = { cgSeasurfaceHeight(cg1id),
+                        cgSeasurfaceHeight(cg1id + 1), cgSeasurfaceHeight(cg1id + smesh->nx + 1),
+                        cgSeasurfaceHeight(cg1id + smesh->nx + 1 + 1) };
+
+                    // compute grad
+                    Eigen::Vector<Nextsim::FloatType, 4> tx = pmap->dX_SSH[eid] * loc_cgSSH;
+                    Eigen::Vector<Nextsim::FloatType, 4> ty = pmap->dY_SSH[eid] * loc_cgSSH;
+
+                    // add global vector
+                    uGrad(cg1id) -= tx(0);
+                    uGrad(cg1id + 1) -= tx(1);
+                    uGrad(cg1id + smesh->nx + 1) -= tx(2);
+                    uGrad(cg1id + smesh->nx + 1 + 1) -= tx(3);
+                    vGrad(cg1id) -= ty(0);
+                    vGrad(cg1id + 1) -= ty(1);
+                    vGrad(cg1id + smesh->nx + 1) -= ty(2);
+                    vGrad(cg1id + smesh->nx + 1 + 1) -= ty(3);
+                }
+            }
+        }
+    }
+
+    // scale with mass
+#pragma omp parallel for
+    for (size_t i = 0; i < uGrad.rows(); ++i) {
+        uGrad(i) /= pmap->lumpedcg1mass(i);
+        vGrad(i) /= pmap->lumpedcg1mass(i);
+    }
+
+    ///// correct boundary (just extend in last elements)
+    size_t cg1row = smesh->nx + 1;
+    size_t topleft = smesh->ny * cg1row;
+    for (size_t i = 1; i < smesh->nx; ++i) // bottom / top
+    {
+        uGrad(i) = uGrad(i + cg1row);
+        vGrad(i) = vGrad(i + cg1row);
+        uGrad(topleft + i) = uGrad(topleft + i - cg1row);
+        vGrad(topleft + i) = vGrad(topleft + i - cg1row);
+    }
+    for (size_t i = 1; i < smesh->ny; ++i) // left / right
+    {
+        uGrad(i * cg1row) = uGrad(i * cg1row + 1);
+        vGrad(i * cg1row) = vGrad(i * cg1row + 1);
+
+        uGrad(i * cg1row + cg1row - 1) = uGrad(i * cg1row + cg1row - 1 - 1);
+        vGrad(i * cg1row + cg1row - 1) = vGrad(i * cg1row + cg1row - 1 - 1);
+    }
+    // corners
+    uGrad(0) = uGrad(cg1row + 1);
+    vGrad(0) = vGrad(cg1row + 1);
+    uGrad(smesh->nx) = uGrad(smesh->nx + cg1row - 1);
+    vGrad(smesh->nx) = vGrad(smesh->nx + cg1row - 1);
+    uGrad(smesh->ny * cg1row) = uGrad((smesh->ny - 1) * cg1row + 1);
+    vGrad(smesh->ny * cg1row) = vGrad((smesh->ny - 1) * cg1row + 1);
+    uGrad((smesh->ny + 1) * cg1row - 1) = uGrad((smesh->ny) * cg1row - 1 - 1);
+    vGrad((smesh->ny + 1) * cg1row - 1) = vGrad((smesh->ny) * cg1row - 1 - 1);
+
+    // Interpolate to CG2 (maybe own function in interpolation?)
+    if (CGDEGREE == 1) {
+        uGradSeasurfaceHeight = uGrad;
+        vGradSeasurfaceHeight = vGrad;
+    } else {
+        // outer nodes
+        size_t icg1 = 0;
+#pragma omp parallel for
+        for (size_t iy = 0; iy <= smesh->ny; ++iy) {
+            size_t icg2 = (2 * smesh->nx + 1) * 2 * iy;
+            for (size_t ix = 0; ix <= smesh->nx; ++ix, ++icg1, icg2 += 2) {
+                uGradSeasurfaceHeight(icg2) = uGrad(icg1);
+                vGradSeasurfaceHeight(icg2) = vGrad(icg1);
+            }
+        }
+        // along lines
+#pragma omp parallel for
+        for (size_t iy = 0; iy <= smesh->ny; ++iy) // horizontal
+        {
+            size_t icg1 = (smesh->nx + 1) * iy;
+            size_t icg2 = (2 * smesh->nx + 1) * 2 * iy + 1;
+            for (size_t ix = 0; ix < smesh->nx; ++ix, ++icg1, icg2 += 2) {
+                uGradSeasurfaceHeight(icg2) = 0.5 * (uGrad(icg1) + uGrad(icg1 + 1));
+                vGradSeasurfaceHeight(icg2) = 0.5 * (vGrad(icg1) + vGrad(icg1 + 1));
+            }
+        }
+#pragma omp parallel for
+        for (size_t iy = 0; iy < smesh->ny; ++iy) // vertical
+        {
+            size_t icg1 = (smesh->nx + 1) * iy;
+            size_t icg2 = (2 * smesh->nx + 1) * (2 * iy + 1);
+            for (size_t ix = 0; ix <= smesh->nx; ++ix, ++icg1, icg2 += 2) {
+                uGradSeasurfaceHeight(icg2) = 0.5 * (uGrad(icg1) + uGrad(icg1 + cg1row));
+                vGradSeasurfaceHeight(icg2) = 0.5 * (vGrad(icg1) + vGrad(icg1 + cg1row));
+            }
+        }
+
+        // midpoints
+#pragma omp parallel for
+        for (size_t iy = 0; iy < smesh->ny; ++iy) // vertical
+        {
+            size_t icg1 = (smesh->nx + 1) * iy;
+            size_t icg2 = (2 * smesh->nx + 1) * (2 * iy + 1) + 1;
+            for (size_t ix = 0; ix < smesh->nx; ++ix, ++icg1, icg2 += 2) {
+                uGradSeasurfaceHeight(icg2) = 0.25
+                    * (uGrad(icg1) + uGrad(icg1 + 1) + uGrad(icg1 + cg1row)
+                        + uGrad(icg1 + cg1row + 1));
+                vGradSeasurfaceHeight(icg2) = 0.25
+                    * (vGrad(icg1) + vGrad(icg1 + 1) + vGrad(icg1 + cg1row)
+                        + vGrad(icg1 + cg1row + 1));
+            }
+        }
+    }
+}
+
 template <int DGadvection> void CGDynamicsKernel<DGadvection>::prepareIteration(const DataMap& data)
 {
     // interpolate ice height and concentration to local cg Variables
@@ -123,6 +261,10 @@ template <int DGadvection> void CGDynamicsKernel<DGadvection>::prepareIteration(
     Interpolations::DG2CG(*smesh, cgA, data.at(ciceName));
     VectorManipulations::CGAveragePeriodic(*smesh, cgA);
 
+    // Reinit the gradient of the sea surface height. Not done by
+    // DataMap as seaSurfaceHeight is always dG(0)
+    ComputeGradientOfSeaSurfaceHeight(DynamicsKernel<DGadvection, DGstressComp>::seaSurfaceHeight);
+
     // limit A to [0,1] and H to [0, ...)
     cgA = cgA.cwiseMin(1.0);
     cgA = cgA.cwiseMax(1.e-4);