Add function SolidBodyRotation

src/atlas/CMakeLists.txt

@@ -757,6 +757,8 @@ util/detail/BlackMagic.h
 util/detail/Cache.h
 util/detail/Debug.h
 util/detail/KDTree.h
+util/function/SolidBodyRotation.h
+util/function/SolidBodyRotation.cc
 util/function/SphericalHarmonic.h
 util/function/SphericalHarmonic.cc
 util/function/VortexRollup.h

src/atlas/util/function/SolidBodyRotation.cc

@@ -0,0 +1,149 @@
+/*
+ * (C) Copyright 2023 ECMWF.
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ * In applying this licence, ECMWF does not waive the privileges and immunities
+ * granted to it by virtue of its status as an intergovernmental organisation
+ * nor does it submit to any jurisdiction.
+ */
+
+#include <cmath>
+
+#include "atlas/util/Constants.h"
+
+#include "atlas/util/function/SolidBodyRotation.h"
+
+namespace atlas {
+
+namespace util {
+
+namespace function {
+
+SolidBodyRotation::SolidBodyRotation(const double beta, const double radius) {
+    sin_beta_ = std::sin(beta * Constants::degreesToRadians());
+    cos_beta_ = std::cos(beta * Constants::degreesToRadians());
+    radius_ = radius;
+}
+
+void SolidBodyRotation::wind(const double lon, const double lat, double& u, double& v) const {
+    double x = lon * Constants::degreesToRadians();
+    double y = lat * Constants::degreesToRadians();
+    double cos_x = std::cos(x);
+    double cos_y = std::cos(y);
+    double sin_x = std::sin(x);
+    double sin_y = std::sin(y);
+    u            = cos_y * cos_beta_ + cos_x * sin_y * sin_beta_;
+    v            = -sin_x * sin_beta_;
+}
+
+double SolidBodyRotation::u(const double lon, const double lat) const {
+    double x = lon * Constants::degreesToRadians();
+    double y = lat * Constants::degreesToRadians();
+    double cos_x = std::cos(x);
+    double cos_y = std::cos(y);
+    double sin_y = std::sin(y);
+    return cos_y * cos_beta_ + cos_x * sin_y * sin_beta_;
+}
+
+double SolidBodyRotation::v(const double lon, const double lat) const {
+    double x = lon * Constants::degreesToRadians();
+    double sin_x = std::sin(x);
+    return -sin_x * sin_beta_;
+}
+
+void SolidBodyRotation::vordiv(const double lon, const double lat, double& vor, double& div) const {
+    double x = lon * Constants::degreesToRadians();
+    double y = lat * Constants::degreesToRadians();
+
+    double cos_x = std::cos(x);
+    double cos_y = std::cos(y);
+    double sin_x = std::sin(x);
+    double sin_y = std::sin(y);
+
+    // Divergence = 1./(R*cos(y)) * ( d/dx( u ) + d/dy( v * cos(y) ) )
+    // Vorticity  = 1./(R*cos(y)) * ( d/dx( v ) - d/dy( u * cos(y) ) )
+    double ddx_u      = -sin_x * sin_y * sin_beta_;
+    double ddy_cosy_v = (-sin_x * sin_beta_) * (-sin_y);
+    double ddx_v      = -cos_x * sin_beta_;
+    double ddy_cosy_u =
+        2 * cos_y * (-sin_y) * cos_beta_ + (-sin_y) * cos_x * sin_y * sin_beta_ + cos_y * cos_x * cos_y * sin_beta_;
+
+    double metric = 1. / (radius_ * cos_y);
+
+    div = metric * (ddx_u + ddy_cosy_v);
+    vor = metric * (ddx_v - ddy_cosy_u);
+}
+
+double SolidBodyRotation::vorticity(const double lon, const double lat) const {
+    double x = lon * Constants::degreesToRadians();
+    double y = lat * Constants::degreesToRadians();
+
+    double cos_x = std::cos(x);
+    double cos_y = std::cos(y);
+    double sin_y = std::sin(y);
+
+    // Vorticity  = 1./(R*cos(y)) * ( d/dx( v ) - d/dy( u * cos(y) ) )
+    double ddx_v      = -cos_x * sin_beta_;
+    double ddy_cosy_u =
+        2 * cos_y * (-sin_y) * cos_beta_ + (-sin_y) * cos_x * sin_y * sin_beta_ + cos_y * cos_x * cos_y * sin_beta_;
+
+    double metric = 1. / (radius_ * cos_y);
+
+    return metric * (ddx_v - ddy_cosy_u);
+}
+
+double SolidBodyRotation::divergence(const double lon, const double lat) const {
+    double x = lon * Constants::degreesToRadians();
+    double y = lat * Constants::degreesToRadians();
+
+    double cos_y = std::cos(y);
+    double sin_x = std::sin(x);
+    double sin_y = std::sin(y);
+
+    // Divergence = 1./(R*cos(y)) * ( d/dx( u ) + d/dy( v * cos(y) ) )
+    double ddx_u      = -sin_x * sin_y * sin_beta_;
+    double ddy_cosy_v = (-sin_x * sin_beta_) * (-sin_y);
+
+    double metric = 1. / (radius_ * cos_y);
+
+    return metric * (ddx_u + ddy_cosy_v);
+}
+
+double SolidBodyRotation::windMagnitude(const double lon, const double lat) const {
+    return std::sqrt(windMagnitudeSquared(lon, lat));
+}
+
+double SolidBodyRotation::windMagnitudeSquared(const double lon, const double lat) const {
+    double u, v;
+    wind(lon, lat, u, v);
+    return u * u + v * v;
+}
+
+void SolidBodyRotation::windMagnitudeSquaredGradient(const double lon, const double lat, double& dfdx, double& dfdy) const {
+    double x = lon * Constants::degreesToRadians();
+    double y = lat * Constants::degreesToRadians();
+
+    double cos_x = std::cos(x);
+    double cos_y = std::cos(y);
+    double sin_x = std::sin(x);
+    double sin_y = std::sin(y);
+
+    double metric_y = 1. / radius_;
+    double metric_x = metric_y / cos_y;
+
+    double u    = cos_y * cos_beta_ + cos_x * sin_y * sin_beta_;
+    double v    = -sin_x * sin_beta_;
+    double dudx = metric_x * (-sin_x * sin_y * sin_beta_);
+    double dudy = metric_y * (-sin_y * cos_beta_ + cos_x * cos_y * sin_beta_);
+    double dvdx = metric_x * (-cos_x * sin_beta_);
+    double dvdy = metric_y * (0.);
+    dfdx        = 2 * u * dudx + 2 * v * dvdx;
+    dfdy        = 2 * u * dudy + 2 * v * dvdy;
+}
+
+}  // namespace function
+
+}  // namespace util
+
+}  // namespace atlas

src/atlas/util/function/SolidBodyRotation.h

@@ -0,0 +1,54 @@
+/*
+ * (C) Copyright 2023 ECMWF.
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ * In applying this licence, ECMWF does not waive the privileges and immunities
+ * granted to it by virtue of its status as an intergovernmental organisation nor
+ * does it submit to any jurisdiction.
+ */
+
+
+#pragma once
+
+namespace atlas {
+
+namespace util {
+
+namespace function {
+
+/// \brief An analytic function that provides solid body rotation winds on a sphere.
+///
+/// All angles must be provided in degrees.
+///
+class SolidBodyRotation {
+public:
+
+    SolidBodyRotation() : SolidBodyRotation(0., 1.) {}
+
+
+    SolidBodyRotation(const double beta) : SolidBodyRotation(beta, 1.) {}
+    SolidBodyRotation(const double beta, const double radius);
+
+    void wind(const double lon, const double lat, double& u, double& v) const;
+    void vordiv(const double lon, const double lat, double& vor, double& div) const;
+    double windMagnitude(const double lon, const double lat) const;
+    double u(const double lon, const double lat) const;
+    double v(const double lon, const double lat) const;
+    double vorticity(const double lon, const double lat) const;
+    double divergence(const double lon, const double lat) const;
+
+    double windMagnitudeSquared(const double lon, const double lat) const;
+    void windMagnitudeSquaredGradient(const double lon, const double lat, double& dfdx, double& dfdy) const;
+
+private:
+    double sin_beta_;
+    double cos_beta_;
+    double radius_;
+};
+
+}  // namespace function
+
+}  // namespace util
+
+}  // namespace atlas

src/sandbox/interpolation/atlas-conservative-interpolation.cc

@@ -45,6 +45,7 @@
 #include "atlas/output/Gmsh.h"
 #include "atlas/runtime/AtlasTool.h"
 #include "atlas/util/Config.h"
+#include "atlas/util/function/SolidBodyRotation.h"
 #include "atlas/util/function/SphericalHarmonic.h"
 #include "atlas/util/function/VortexRollup.h"
 
@@ -107,7 +108,8 @@ class AtlasParallelInterpolation : public AtlasTool {
         add_option(new eckit::option::Separator("Initial condition options"));
 
         add_option(new SimpleOption<std::string>(
-            "init", "Setup initial source field [ constant, spherical_harmonic, vortex_rollup (default) ]"));
+            "init", "Setup initial source field [ constant, spherical_harmonic, vortex_rollup (default), solid_body_rotation_wind_magnitude ]"));
+        add_option(new SimpleOption<double>("solid_body_rotation.angle", "Angle of solid body rotation (default = 0.)"));
         add_option(new SimpleOption<double>("vortex_rollup.t", "Value that controls vortex rollup (default = 0.5)"));
         add_option(new SimpleOption<double>("constant.value", "Value that is assigned in case init==constant)"));
         add_option(new SimpleOption<long>("spherical_harmonic.n", "total wave number 'n' of a spherical harmonic"));
@@ -147,6 +149,12 @@ std::function<double(const PointLonLat&)> get_init(const AtlasTool::Args& args)
         args.get("constant.value", value = 1.);
         return [value](const PointLonLat&) { return value; };
     }
+    else if (init == "solid_body_rotation_wind_magnitude") {
+        double beta;
+        args.get("solid_body_rotation.angle", beta = 0.);
+        util::function::SolidBodyRotation sbr(beta);
+        return [sbr](const PointLonLat& p) { return sbr.windMagnitude(p.lon(), p.lat()); };
+    }
     else {
         if (args.has("init")) {
             Log::error() << "Bad value for \"init\": \"" << init << "\" not recognised." << std::endl;

src/tests/numerics/test_fvm_nabla_validation.cc

@@ -32,6 +32,7 @@
 #include "atlas/util/Constants.h"
 #include "atlas/util/CoordinateEnums.h"
 #include "atlas/util/Earth.h"
+#include "atlas/util/function/SolidBodyRotation.h"
 
 #include "tests/AtlasTestEnvironment.h"
 
@@ -79,77 +80,10 @@ static int metric_approach() {
     return 0;
 }
 
-class SolidBodyRotation {
-    double radius;
-    double beta;
-    double sin_beta;
-    double cos_beta;
-
-public:
-    SolidBodyRotation(const double _radius, const double _beta): radius{_radius}, beta{_beta} {
-        sin_beta = std::sin(beta);
-        cos_beta = std::cos(beta);
-    }
-    void wind(const double x, const double y, double& u, double& v) {
-        double cos_x = std::cos(x);
-        double cos_y = std::cos(y);
-        double sin_x = std::sin(x);
-        double sin_y = std::sin(y);
-        u            = cos_y * cos_beta + cos_x * sin_y * sin_beta;
-        v            = -sin_x * sin_beta;
-    }
-
-    void vordiv(const double x, const double y, double& vor, double& div) {
-        double cos_x = std::cos(x);
-        double cos_y = std::cos(y);
-        double sin_x = std::sin(x);
-        double sin_y = std::sin(y);
-
-        // Divergence = 1./(R*cos(y)) * ( d/dx( u ) + d/dy( v * cos(y) ) )
-        // Vorticity  = 1./(R*cos(y)) * ( d/dx( v ) - d/dy( u * cos(y) ) )
-        double ddx_u      = -sin_x * sin_y * sin_beta;
-        double ddy_cosy_v = (-sin_x * sin_beta) * (-sin_y);
-        double ddx_v      = -cos_x * sin_beta;
-        double ddy_cosy_u =
-            2 * cos_y * (-sin_y) * cos_beta + (-sin_y) * cos_x * sin_y * sin_beta + cos_y * cos_x * cos_y * sin_beta;
-
-        double metric = 1. / (radius * cos_y);
-
-        div = metric * (ddx_u + ddy_cosy_v);
-        vor = metric * (ddx_v - ddy_cosy_u);
-    }
-
-    void wind_magnitude_squared(const double x, const double y, double& f) {
-        double u, v;
-        wind(x, y, u, v);
-        f = u * u + v * v;
-    }
-
-    void wind_magnitude_squared_gradient(const double x, const double y, double& dfdx, double& dfdy) {
-        double cos_x = std::cos(x);
-        double cos_y = std::cos(y);
-        double sin_x = std::sin(x);
-        double sin_y = std::sin(y);
-
-        double metric_y = 1. / radius;
-        double metric_x = metric_y / cos_y;
-
-        double u    = cos_y * cos_beta + cos_x * sin_y * sin_beta;
-        double v    = -sin_x * sin_beta;
-        double dudx = metric_x * (-sin_x * sin_y * sin_beta);
-        double dudy = metric_y * (-sin_y * cos_beta + cos_x * cos_y * sin_beta);
-        double dvdx = metric_x * (-cos_x * sin_beta);
-        double dvdy = metric_y * (0.);
-        dfdx        = 2 * u * dudx + 2 * v * dvdx;
-        dfdy        = 2 * u * dudy + 2 * v * dvdy;
-    }
-};
-
 FieldSet analytical_fields(const fvm::Method& fvm) {
-    constexpr double deg2rad = M_PI / 180.;
     const double radius      = fvm.radius();
 
-    auto lonlat_deg = array::make_view<double, 2>(fvm.mesh().nodes().lonlat());
+    auto lonlat = array::make_view<double, 2>(fvm.mesh().nodes().lonlat());
 
     FieldSet fields;
     auto add_scalar_field = [&](const std::string& name) {
@@ -169,20 +103,20 @@ FieldSet analytical_fields(const fvm::Method& fvm) {
     auto div    = add_scalar_field("ref_div");
     auto vor    = add_scalar_field("ref_vor");
 
-    auto flow          = SolidBodyRotation{radius, beta_in_degrees() * deg2rad};
+    auto flow          = atlas::util::function::SolidBodyRotation{beta_in_degrees(), radius};
     auto is_ghost      = array::make_view<int, 1>(fvm.mesh().nodes().ghost());
     const idx_t nnodes = fvm.mesh().nodes().size();
     for (idx_t jnode = 0; jnode < nnodes; ++jnode) {
         if (is_ghost(jnode)) {
             continue;
         }
-        double x = lonlat_deg(jnode, LON) * deg2rad;
-        double y = lonlat_deg(jnode, LAT) * deg2rad;
+        double x = lonlat(jnode, LON);
+        double y = lonlat(jnode, LAT);
 
         flow.wind(x, y, u(jnode), v(jnode));
         flow.vordiv(x, y, vor(jnode), div(jnode));
-        flow.wind_magnitude_squared(x, y, f(jnode));
-        flow.wind_magnitude_squared_gradient(x, y, dfdx(jnode), dfdy(jnode));
+        f(jnode) = flow.windMagnitudeSquared(x, y);
+        flow.windMagnitudeSquaredGradient(x, y, dfdx(jnode), dfdy(jnode));
 
         uv(jnode, XX)     = u(jnode);
         uv(jnode, YY)     = v(jnode);