Make VortexRollup function publicly available, and remove it from var…

…ious places (#87)

src/atlas/CMakeLists.txt

@@ -728,6 +728,8 @@ util/detail/BlackMagic.h
 util/detail/Cache.h
 util/detail/Debug.h
 util/detail/KDTree.h
+util/function/VortexRollup.h
+util/function/VortexRollup.cc
 )
 
 list( APPEND atlas_internals_srcs

src/atlas/util/function/VortexRollup.cc

@@ -0,0 +1,43 @@
+/*
+ * (C) Copyright 2013 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 "atlas/util/Constants.h"
+#include "atlas/util/Earth.h"
+
+#include "atlas/util/function/VortexRollup.h"
+
+namespace atlas {
+
+namespace util {
+
+namespace function {
+
+double vortex_rollup(double lon, double lat, double t) {
+    lon *= Constants::degreesToRadians();
+    lat *= Constants::degreesToRadians();
+
+    auto sqr                  = [](const double x) { return x * x; };
+    auto sech                 = [](const double x) { return 1. / std::cosh(x); };
+    constexpr double two_pi   = 2. * M_PI;
+    const double lambda_prime = std::atan2(-std::cos(lon - two_pi * t), std::tan(lat));
+    const double rho          = 3. * std::sqrt(1. - sqr(std::cos(lat)) * sqr(std::sin(lon - two_pi * t)));
+    double omega              = 0.;
+    double a                  = Earth::radius();
+    if (rho != 0.) {
+        omega = 0.5 * 3 * std::sqrt(3) * a * two_pi * sqr(sech(rho)) * std::tanh(rho) / rho;
+    }
+    return -std::tanh(0.2 * rho * std::sin(lambda_prime - omega / a * t));
+}
+
+}  // namespace function
+
+}  // namespace util
+
+}  // namespace atlas

src/atlas/util/function/VortexRollup.h

@@ -0,0 +1,39 @@
+/*
+ * (C) Copyright 2013 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 a vortex rollup on a 2D sphere
+///
+/// \detailed The formula is found in
+///           "A Lagrangian Particle Method with Remeshing for Tracer Transport on the Sphere"
+///           by Peter Bosler, James Kent, Robert Krasny, Christiane Jablonowski, JCP 2015
+///           as the tracer density in Test Case 1.
+///           The longitude (lon) and latitude (lat) are assumed to be in degrees,
+///           The time parameter associated with the vortex rollup is set by (t).
+///
+///           The time it takes for the counter-rotating vortices along
+///           the equator to return to its original position takes
+///           time t = 1.0;
+///
+double vortex_rollup(double lon, double lat, double t);
+
+}  // namespace function
+
+}  // namespace util
+
+}  // namespace atlas

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

@@ -14,39 +14,17 @@
 #include "eckit/config/Resource.h"
 #include "eckit/log/Plural.h"
 
+#include "PartitionedMesh.h"
 #include "atlas/field.h"
 #include "atlas/functionspace.h"
 #include "atlas/interpolation.h"
 #include "atlas/linalg/sparse/Backend.h"
 #include "atlas/runtime/AtlasTool.h"
 #include "atlas/runtime/Log.h"
-
-#include "PartitionedMesh.h"
+#include "atlas/util/function/VortexRollup.h"
 
 using namespace atlas;
 
-auto vortex_rollup = [](double lon, double lat, double t) {
-    // lon and lat in radians!
-
-    // Formula found in "A Lagrangian Particle Method with Remeshing for Tracer Transport on the Sphere"
-    // by Peter Bosler, James Kent, Robert Krasny, CHristiane Jablonowski, JCP 2015
-
-    auto sqr           = [](const double x) { return x * x; };
-    auto sech          = [](const double x) { return 1. / std::cosh(x); };
-    const double T     = 1.;
-    const double Omega = 2. * M_PI / T;
-    t *= T;
-    const double lambda_prime = std::atan2(-std::cos(lon - Omega * t), std::tan(lat));
-    const double rho          = 3. * std::sqrt(1. - sqr(std::cos(lat)) * sqr(std::sin(lon - Omega * t)));
-    double omega              = 0.;
-    double a                  = util::Earth::radius();
-    if (rho != 0.) {
-        omega = 0.5 * 3 * std::sqrt(3) * a * Omega * sqr(sech(rho)) * std::tanh(rho) / rho;
-    }
-    double q = 1. - std::tanh(0.2 * rho * std::sin(lambda_prime - omega / a * t));
-    return q;
-};
-
 class AtlasParallelInterpolation : public AtlasTool {
     int execute(const AtlasTool::Args& args) override;
     std::string briefDescription() override { return "Demonstration of parallel interpolation"; }
@@ -245,7 +223,7 @@ int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) {
             //            src_scalar_1( j ) = -std::tanh( y / 10. * std::cos( 50. / std::sqrt( x * x + y * y ) ) -
             //                                            x / 10. * std::sin( 50. / std::sqrt( x * x + y * y ) ) );
 
-            src_scalar_1(j) = vortex_rollup(lon, lat, 1.);
+            src_scalar_1(j) = util::function::vortex_rollup(lonlat(j, 0), lonlat(j, 1), 1.);
         }
     }
 

src/sandbox/interpolation/atlas-parallel-structured-interpolation.cc

@@ -19,6 +19,7 @@
 #include "atlas/runtime/AtlasTool.h"
 #include "atlas/runtime/Log.h"
 #include "atlas/util/CoordinateEnums.h"
+#include "atlas/util/function/VortexRollup.h"
 
 
 using namespace atlas;
@@ -87,32 +88,6 @@ static Config processed_config(const eckit::Configuration& _config) {
     return config;
 }
 
-
-double vortex_rollup(double lon, double lat, double t) {
-    // lon and lat in degrees!
-
-    // Formula found in "A Lagrangian Particle Method with Remeshing for Tracer Transport on the Sphere"
-    // by Peter Bosler, James Kent, Robert Krasny, Christiane Jablonowski, JCP 2015
-
-    lon *= M_PI / 180.;
-    lat *= M_PI / 180.;
-
-    auto sqr           = [](const double x) { return x * x; };
-    auto sech          = [](const double x) { return 1. / std::cosh(x); };
-    const double T     = 1.;
-    const double Omega = 2. * M_PI / T;
-    t *= T;
-    const double lambda_prime = std::atan2(-std::cos(lon - Omega * t), std::tan(lat));
-    const double rho          = 3. * std::sqrt(1. - sqr(std::cos(lat)) * sqr(std::sin(lon - Omega * t)));
-    double omega              = 0.;
-    double a                  = util::Earth::radius();
-    if (rho != 0.) {
-        omega = 0.5 * 3 * std::sqrt(3) * a * Omega * sqr(sech(rho)) * std::tanh(rho) / rho;
-    }
-    double q = 1. - std::tanh(0.2 * rho * std::sin(lambda_prime - omega / a * t));
-    return q;
-};
-
 int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) {
     ATLAS_TRACE("AtlasParallelInterpolation::execute");
     auto source_gridname = args.getString("source", "O32");
@@ -180,7 +155,7 @@ int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) {
         else {
             idx_t k = args.getInt("vortex-rollup", 0);
             atlas_omp_parallel_for(idx_t n = 0; n < size; ++n) {
-                source(n) = vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2);
+                source(n) = util::function::vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2);
             }
         }
     }

src/tests/acceptance_tests/atest_mgrids.cc

@@ -29,11 +29,10 @@
 #include "atlas/util/Config.h"
 #include "atlas/util/Constants.h"
 #include "atlas/util/CoordinateEnums.h"
+#include "atlas/util/function/VortexRollup.h"
 
 using namespace atlas;
 
-double vortex_rollup(double lon, double lat, double t, double mean);
-
 //------------------------------------------------------------------------------
 
 class Program : public AtlasTool {
@@ -100,7 +99,7 @@ int Program::execute(const Args& args) {
 
     double meanA = 1.;
     for (idx_t n = 0; n < fsA.size(); ++n) {
-        A(n) = vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 1., meanA);
+        A(n) = meanA + util::function::vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 1.);
     }
     fieldA.set_dirty();
     fieldA.haloExchange();
@@ -204,33 +203,6 @@ int Program::execute(const Args& args) {
 
 //------------------------------------------------------------------------------
 
-double vortex_rollup(double lon, double lat, double t, double mean) {
-    // lon and lat in degrees!
-
-    // Formula found in "A Lagrangian Particle Method with Remeshing for Tracer Transport on the Sphere"
-    // by Peter Bosler, James Kent, Robert Krasny, CHristiane Jablonowski, JCP 2015
-
-    lon *= util::Constants::degreesToRadians();
-    lat *= util::Constants::degreesToRadians();
-
-    auto sqr           = [](const double x) { return x * x; };
-    auto sech          = [](const double x) { return 1. / std::cosh(x); };
-    const double T     = 1.;
-    const double Omega = 2. * M_PI / T;
-    t *= T;
-    const double lambda_prime = std::atan2(-std::cos(lon - Omega * t), std::tan(lat));
-    const double rho          = 3. * std::sqrt(1. - sqr(std::cos(lat)) * sqr(std::sin(lon - Omega * t)));
-    double omega              = 0.;
-    double a                  = util::Earth::radius();
-    if (rho != 0.) {
-        omega = 0.5 * 3 * std::sqrt(3) * a * Omega * sqr(sech(rho)) * std::tanh(rho) / rho;
-    }
-    double q = mean - std::tanh(0.2 * rho * std::sin(lambda_prime - omega / a * t));
-    return q;
-};
-
-//------------------------------------------------------------------------------
-
 int main(int argc, char** argv) {
     Program tool(argc, argv);
     return tool.start();

src/tests/functionspace/test_cubedsphere_functionspace.cc

@@ -18,6 +18,7 @@
 #include "atlas/meshgenerator/detail/cubedsphere/CubedSphereUtility.h"
 #include "atlas/option.h"
 #include "atlas/util/CoordinateEnums.h"
+#include "atlas/util/function/VortexRollup.h"
 
 #include "tests/AtlasTestEnvironment.h"
 
@@ -27,11 +28,6 @@ namespace test {
 // Allow small differences in the last few bits of a double aprroximately equal to 1
 constexpr double epsilon = std::numeric_limits<double>::epsilon() * 16;
 
-
-double testFunction(double lon, double lat) {
-    return std::sin(3 * lon * M_PI / 180) * std::sin(2 * lat * M_PI / 180);
-}
-
 template <typename BaseFunctionSpace>
 void testFunctionSpace(const functionspace::CubedSphereColumns<BaseFunctionSpace>& functionspace) {
     // Make field.
@@ -63,7 +59,7 @@ void testFunctionSpace(const functionspace::CubedSphereColumns<BaseFunctionSpace
         EXPECT(!ghostView(index));
 
         // Set field values.
-        fieldView(index) = testFunction(lonLatView(index, LON), lonLatView(index, LAT));
+        fieldView(index) = util::function::vortex_rollup(lonLatView(index, LON), lonLatView(index, LAT), 1.0);
         ++testFuncCallCount;
     });
 
@@ -80,7 +76,8 @@ void testFunctionSpace(const functionspace::CubedSphereColumns<BaseFunctionSpace
         EXPECT(index == functionspace.index(t, i, j));
 
         // Set field values.
-        EXPECT_APPROX_EQ(fieldView(index), testFunction(lonLatView(index, LON), lonLatView(index, LAT)), epsilon);
+        EXPECT_APPROX_EQ(fieldView(index),
+                         util::function::vortex_rollup(lonLatView(index, LON), lonLatView(index, LAT), 1.0), epsilon);
 
         ++testFuncCallCount;
     });

src/tests/interpolation/test_interpolation_structured2D.cc

@@ -21,6 +21,7 @@
 #include "atlas/meshgenerator.h"
 #include "atlas/output/Gmsh.h"
 #include "atlas/util/CoordinateEnums.h"
+#include "atlas/util/function/VortexRollup.h"
 
 #include "tests/AtlasTestEnvironment.h"
 
@@ -116,31 +117,6 @@ FunctionSpace output_functionspace(const Grid& grid) {
     return NodeColumns{output_mesh};
 }
 
-double vortex_rollup(double lon, double lat, double t) {
-    // lon and lat in degrees!
-
-    // Formula found in "A Lagrangian Particle Method with Remeshing for Tracer Transport on the Sphere"
-    // by Peter Bosler, James Kent, Robert Krasny, CHristiane Jablonowski, JCP 2015
-
-    lon *= M_PI / 180.;
-    lat *= M_PI / 180.;
-
-    auto sqr           = [](const double x) { return x * x; };
-    auto sech          = [](const double x) { return 1. / std::cosh(x); };
-    const double T     = 1.;
-    const double Omega = 2. * M_PI / T;
-    t *= T;
-    const double lambda_prime = std::atan2(-std::cos(lon - Omega * t), std::tan(lat));
-    const double rho          = 3. * std::sqrt(1. - sqr(std::cos(lat)) * sqr(std::sin(lon - Omega * t)));
-    double omega              = 0.;
-    double a                  = util::Earth::radius();
-    if (rho != 0.) {
-        omega = 0.5 * 3 * std::sqrt(3) * a * Omega * sqr(sech(rho)) * std::tanh(rho) / rho;
-    }
-    double q = 1. - std::tanh(0.2 * rho * std::sin(lambda_prime - omega / a * t));
-    return q;
-};
-
 CASE("which scheme?") {
     Log::info() << scheme().getString("type") << std::endl;
 }
@@ -162,7 +138,7 @@ CASE("test_interpolation_structured using functionspace API") {
         auto lonlat = array::make_view<double, 2>(input_fs.xy());
         auto source = array::make_view<double, 1>(field_source);
         for (idx_t n = 0; n < input_fs.size(); ++n) {
-            source(n) = vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 1.);
+            source(n) = util::function::vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 1.);
         }
 
         EXPECT(field_source.dirty());
@@ -202,7 +178,7 @@ CASE("test_interpolation_structured using grid API") {
 
         idx_t n{0};
         for (auto p : input_grid.lonlat()) {
-            src_data[n++] = vortex_rollup(p.lon(), p.lat(), 1.);
+            src_data[n++] = util::function::vortex_rollup(p.lon(), p.lat(), 1.);
         }
 
         // Wrap memory in atlas Fields and interpolate
@@ -259,7 +235,7 @@ CASE("test_interpolation_structured using fs API multiple levels") {
     auto source = array::make_view<double, 2>(field_source);
     for (idx_t n = 0; n < input_fs.size(); ++n) {
         for (idx_t k = 0; k < 3; ++k) {
-            source(n, k) = vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2);
+            source(n, k) = util::function::vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2);
         }
     }
     interpolation.execute(field_source, field_target);
@@ -281,7 +257,7 @@ struct AdjointTolerance {
 template <>
 const double AdjointTolerance<double>::value = 2.e-14;
 template <>
-const float AdjointTolerance<float>::value = 4.e-6;
+const float AdjointTolerance<float>::value = 2.e-5;
 
 
 template <typename Value>
@@ -307,7 +283,7 @@ void test_interpolation_structured_using_fs_API_for_fieldset() {
         auto source = array::make_view<Value, 2>(field_source);
         for (idx_t n = 0; n < input_fs.size(); ++n) {
             for (idx_t k = 0; k < 3; ++k) {
-                source(n, k) = vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2);
+                source(n, k) = util::function::vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2);
             }
         }
     }
@@ -483,7 +459,7 @@ CASE("ATLAS-315: Target grid with domain West of 0 degrees Lon") {
     auto source        = array::make_view<double, 1>(field_src);
     constexpr double k = 1;
     for (idx_t n = 0; n < source.size(); ++n) {
-        source(n) = vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2);
+        source(n) = util::function::vortex_rollup(lonlat(n, LON), lonlat(n, LAT), 0.5 + double(k) / 2);
     }
 
     interpolation.execute(field_src, field_tgt);