Merge branch 'release/0.28.0'

* release/0.28.0: (52 commits) Update Changelog Version 0.28.0 ATLAS-344 Fix ConvexSphericalPolygon::compute_radius() ATLAS-344 Remove part of test_convexsphericalpolygon ATLAS-344 Use EPS and TOL in ConvexSphericalPolygon where possible Apply clang-format Make Docker-Ubuntu1604 happy (GNU 5.4) feature: bilinear-remapping interpolation (#85) Introduce regional_variable_resolution grid type via a new GridBuilder (#92) ATLAS-347 Fix build_nodes_remote_idx for Healpix by fixing the partition of periodic BC|EAST points Improve build_nodes_remote_idx error reporting ATLAS-347 remove not used periodic nodes from the partitions - fails in building remote indices Add default constructor for Projection::Jacobian (#93) Added methods to projection::Jacobian (#93) Add macro to check eckit version ATLAS-344 Fix PolygonXY::contains, requiring tolerances at edges ATLAS-346 Fix CellColumns::haloExchange for multiple element types Add += and -= operators to IndexView elements ATLAS-322 Add more tests to atlas_test_convexsphericalpolygon ATLAS-322 Refactor ConvexSphericalPolygon ...
ecmwf · Mar 2, 2022 · 3901c05 · 3901c05
2 parents d825fad + 4da9a9e
commit 3901c05
Show file tree

Hide file tree

Showing 114 changed files with 6,486 additions and 826 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -7,6 +7,27 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html
 
 ## [Unreleased]
 
+## [0.28.0] - 2021-03-02
+### Added
+- Assignment of ArrayView from ArrayView
+- Grid "regional_variable_resolution" via a new VariableResolutionProjection
+- CubedSphereDualMeshGenerator
+- VortexRollup function as analytical field for initialising data
+- ConvexSphericalPolygon utility class
+- Improve Projection::Jacobian
+- Initial implementation for bilinear interpolation for unstructured meshes
+
+### Changed
+- Use new eckit (1.19.0) Sparse and Dense linear algebra API
+- General robustness improvements to CubedSphere to using functionspaces with various halos
+
+### Fixed
+- Workarounds to fix compilation with Fujitsu compiler
+- Workarounds to avoid Cray compiler problems with certain flag combinations
+- CellColumns::haloExchange for meshes with multiple element types
+- Computation of HEALPix mesh remote indices.
+
+
 ## [0.27.0] - 2021-12-03
 ### Added
 - Adjoint interpolation with some restrictions
@@ -320,6 +341,7 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html
 ## 0.13.0 - 2018-02-16
 
 [Unreleased]: https://github.com/ecmwf/atlas/compare/master...develop
+[0.28.0]: https://github.com/ecmwf/atlas/compare/0.27.0...0.28.0
 [0.27.0]: https://github.com/ecmwf/atlas/compare/0.26.0...0.27.0
 [0.26.0]: https://github.com/ecmwf/atlas/compare/0.25.0...0.26.0
 [0.25.0]: https://github.com/ecmwf/atlas/compare/0.24.1...0.25.0

diff --git a/VERSION b/VERSION
@@ -1 +1 @@
-0.27.0
+0.28.0
diff --git a/doc/example-grids/regional_variable_resolution.yml b/doc/example-grids/regional_variable_resolution.yml
@@ -0,0 +1,29 @@
+# (C) Crown copyright 2022, Met Office
+#
+# 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.
+# regional projection using variable resolution
+
+type : "regional_variable_resolution"
+progression : 1.13
+inner :
+      xmin : 351.386944827586319
+      ymin : -5.008754172413662
+      xend : 366.363355172413776
+      yend : 8.665359620690706
+      dx : 0.6511482758621128
+outer :
+      xmin : 348.13120344827576
+      xend : 369.6190965517242
+      ymin : -8.264495551724226
+      yend : 11.92110100000127
+      dx : 1
+      width: 4.
+
+check :
+  size : 1088
+  lonlat(first) : [346.9838795150958504,-9.41181948490414122]
+  lonlat(last)  : [370.7664204849036764,13.06842493318118414]
+  bounding_box(n,w,s,e) : [13.06842493318118414,346.9838795150958504,-9.41181948490414122,370.766420484903676]
+  uid : aca7861d7ddeb4f2448c57484946d44e
+
diff --git a/doc/example-grids/update_uid.py b/doc/example-grids/update_uid.py
@@ -16,7 +16,7 @@ def run_program_get_error(file):
 
     from subprocess import Popen, PIPE
 
-    p = Popen(command, stdout=PIPE, stderr=PIPE)
+    p = Popen(command, stdout=PIPE, stderr=PIPE, shell=True)
     output, error = p.communicate()
     if p.returncode != 0:
        return str(error.strip(),'ascii')

diff --git a/src/apps/atlas-meshgen.cc b/src/apps/atlas-meshgen.cc
@@ -125,6 +125,7 @@ class Meshgen2Gmsh : public AtlasTool {
     std::string key;
     long halo;
     bool edges;
+    bool cells;
     bool brick;
     bool stats;
     bool info;
@@ -165,6 +166,7 @@ Meshgen2Gmsh::Meshgen2Gmsh(int argc, char** argv): AtlasTool(argc, argv) {
     add_option(new Separator("Advanced"));
     add_option(new SimpleOption<long>("halo", "Halo size"));
     add_option(new SimpleOption<bool>("edges", "Build edge datastructure"));
+    add_option(new SimpleOption<bool>("cells", "Build cells datastructure"));
     add_option(new SimpleOption<bool>("brick", "Build brick dual mesh"));
     add_option(new SimpleOption<bool>("stats", "Write statistics file"));
     add_option(new SimpleOption<bool>("info", "Write Info"));
@@ -203,6 +205,8 @@ int Meshgen2Gmsh::execute(const Args& args) {
 
     edges = false;
     args.get("edges", edges);
+    cells = false;
+    args.get("cells", cells);
     stats = false;
     args.get("stats", stats);
     info = false;
@@ -280,6 +284,7 @@ int Meshgen2Gmsh::execute(const Args& args) {
         throw;
     }
 
+
     if (grid.projection().units() == "degrees") {
         functionspace::NodeColumns nodes_fs(mesh, option::halo(halo));
     }
@@ -299,6 +304,10 @@ int Meshgen2Gmsh::execute(const Args& args) {
         }
     }
 
+    if (cells) {
+        functionspace::CellColumns cells_fs(mesh, option::halo(halo));
+    }
+
     if (stats) {
         build_statistics(mesh);
     }

diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt
@@ -93,6 +93,8 @@ projection/detail/ProjectionImpl.h
 projection/detail/ProjectionUtilities.h
 projection/detail/SchmidtProjection.cc
 projection/detail/SchmidtProjection.h
+projection/detail/VariableResolutionProjection.cc
+projection/detail/VariableResolutionProjection.h
 
 domain.h
 domain/Domain.cc
@@ -158,6 +160,8 @@ grid/detail/grid/LonLat.h
 grid/detail/grid/LonLat.cc
 grid/detail/grid/Regional.h
 grid/detail/grid/Regional.cc
+grid/detail/grid/RegionalVariableResolution.h
+grid/detail/grid/RegionalVariableResolution.cc
 grid/detail/grid/Healpix.h
 grid/detail/grid/Healpix.cc
 
@@ -332,6 +336,8 @@ mesh/actions/ExtendNodesGlobal.cc
 mesh/actions/BuildDualMesh.h
 mesh/actions/BuildCellCentres.cc
 mesh/actions/BuildCellCentres.h
+mesh/actions/Build2DCellCentres.cc
+mesh/actions/Build2DCellCentres.h
 mesh/actions/BuildConvexHull3D.cc
 mesh/actions/BuildConvexHull3D.h
 mesh/actions/BuildDualMesh.cc
@@ -365,6 +371,8 @@ meshgenerator/MeshGenerator.cc
 meshgenerator/MeshGenerator.h
 meshgenerator/detail/CubedSphereMeshGenerator.h
 meshgenerator/detail/CubedSphereMeshGenerator.cc
+meshgenerator/detail/CubedSphereDualMeshGenerator.h
+meshgenerator/detail/CubedSphereDualMeshGenerator.cc
 meshgenerator/detail/NodalCubedSphereMeshGenerator.h
 meshgenerator/detail/NodalCubedSphereMeshGenerator.cc
 meshgenerator/detail/DelaunayMeshGenerator.h
@@ -525,8 +533,12 @@ interpolation/Vector2D.cc
 interpolation/Vector2D.h
 interpolation/Vector3D.cc
 interpolation/Vector3D.h
+interpolation/element/Quad2D.h
+interpolation/element/Quad2D.cc
 interpolation/element/Quad3D.cc
 interpolation/element/Quad3D.h
+interpolation/element/Triag2D.cc
+interpolation/element/Triag2D.h
 interpolation/element/Triag3D.cc
 interpolation/element/Triag3D.h
 interpolation/method/Intersect.cc
@@ -537,12 +549,16 @@ interpolation/method/MethodFactory.cc
 interpolation/method/MethodFactory.h
 interpolation/method/PointIndex3.cc
 interpolation/method/PointIndex3.h
+interpolation/method/PointIndex2.cc
+interpolation/method/PointIndex2.h
 interpolation/method/PointSet.cc
 interpolation/method/PointSet.h
 interpolation/method/Ray.cc
 interpolation/method/Ray.h
 interpolation/method/fe/FiniteElement.cc
 interpolation/method/fe/FiniteElement.h
+interpolation/method/bil/BilinearRemapping.cc
+interpolation/method/bil/BilinearRemapping.h
 interpolation/method/knn/GridBox.cc
 interpolation/method/knn/GridBox.h
 interpolation/method/knn/GridBoxAverage.cc
@@ -650,6 +666,7 @@ array/helpers/ArrayInitializer.h
 array/helpers/ArrayAssigner.h
 array/helpers/ArrayWriter.h
 array/helpers/ArraySlicer.h
+array/helpers/ArrayCopier.h
 #array/Table.h
 #array/Table.cc
 #array/TableView.h
@@ -696,6 +713,8 @@ runtime/Exception.h
 util/Config.cc
 util/Config.h
 util/Constants.h
+util/ConvexSphericalPolygon.cc
+util/ConvexSphericalPolygon.h
 util/Earth.h
 util/GaussianLatitudes.cc
 util/GaussianLatitudes.h
@@ -723,6 +742,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

diff --git a/src/atlas/array/ArrayView.h b/src/atlas/array/ArrayView.h
@@ -31,6 +31,10 @@ namespace array {
     EXPLICIT_TEMPLATE_DECLARATION_TYPE_RANK(float, RANK); \
     EXPLICIT_TEMPLATE_DECLARATION_TYPE_RANK(double, RANK);
 
+// The Fujitsu compiler complains about missing references with the below
+// so we just skip it in that case
+#if !defined(__FUJITSU)
+
 // For each RANK in [1..9]
 EXPLICIT_TEMPLATE_DECLARATION(1)
 EXPLICIT_TEMPLATE_DECLARATION(2)
@@ -42,6 +46,8 @@ EXPLICIT_TEMPLATE_DECLARATION(7)
 EXPLICIT_TEMPLATE_DECLARATION(8)
 EXPLICIT_TEMPLATE_DECLARATION(9)
 
+#endif
+
 #undef EXPLICIT_TEMPLATE_DECLARATION_TYPE_RANK
 #undef EXPLICIT_TEMPLATE_DECLARATION
 

diff --git a/src/atlas/array/gridtools/GridToolsArray.cc b/src/atlas/array/gridtools/GridToolsArray.cc
@@ -162,49 +162,63 @@ class ArrayT_impl {
             case 1:
                 return construct(shape[0]);
             case 2: {
-                if (layout[0] == 0 && layout[1] == 1)
+                if (layout[0] == 0 && layout[1] == 1) {
                     return construct_with_layout<::gridtools::layout_map<0, 1>>(shape[0], shape[1]);
-                if (layout[0] == 1 && layout[1] == 0)
+                }
+                if (layout[0] == 1 && layout[1] == 0) {
                     return construct_with_layout<::gridtools::layout_map<1, 0>>(shape[0], shape[1]);
+                }
             }
             case 3: {
-                if (layout[0] == 0 && layout[1] == 1 && layout[2] == 2)
+                if (layout[0] == 0 && layout[1] == 1 && layout[2] == 2) {
                     return construct_with_layout<::gridtools::layout_map<0, 1, 2>>(shape[0], shape[1], shape[2]);
-                if (layout[0] == 0 && layout[1] == 2 && layout[2] == 1)
+                }
+                if (layout[0] == 0 && layout[1] == 2 && layout[2] == 1) {
                     return construct_with_layout<::gridtools::layout_map<0, 2, 1>>(shape[0], shape[1], shape[2]);
-                if (layout[0] == 1 && layout[1] == 0 && layout[2] == 2)
+                }
+                if (layout[0] == 1 && layout[1] == 0 && layout[2] == 2) {
                     return construct_with_layout<::gridtools::layout_map<1, 0, 2>>(shape[0], shape[1], shape[2]);
-                if (layout[0] == 1 && layout[1] == 2 && layout[2] == 0)
+                }
+                if (layout[0] == 1 && layout[1] == 2 && layout[2] == 0) {
                     return construct_with_layout<::gridtools::layout_map<1, 2, 0>>(shape[0], shape[1], shape[2]);
-                if (layout[0] == 2 && layout[1] == 0 && layout[2] == 1)
+                }
+                if (layout[0] == 2 && layout[1] == 0 && layout[2] == 1) {
                     return construct_with_layout<::gridtools::layout_map<2, 0, 1>>(shape[0], shape[1], shape[2]);
-                if (layout[0] == 2 && layout[1] == 1 && layout[2] == 0)
+                }
+                if (layout[0] == 2 && layout[1] == 1 && layout[2] == 0) {
                     return construct_with_layout<::gridtools::layout_map<2, 1, 0>>(shape[0], shape[1], shape[2]);
+                }
             }
             case 4: {
-                if (layout[0] == 0 && layout[1] == 1 && layout[2] == 2 && layout[3] == 3)
+                if (layout[0] == 0 && layout[1] == 1 && layout[2] == 2 && layout[3] == 3) {
                     return construct_with_layout<::gridtools::layout_map<0, 1, 2, 3>>(shape[0], shape[1], shape[2],
                                                                                       shape[3]);
-                if (layout[0] == 3 && layout[1] == 2 && layout[2] == 1 && layout[3] == 0)
+                }
+                if (layout[0] == 3 && layout[1] == 2 && layout[2] == 1 && layout[3] == 0) {
                     return construct_with_layout<::gridtools::layout_map<3, 2, 1, 0>>(shape[0], shape[1], shape[2],
                                                                                       shape[3]);
+                }
             }
             case 5: {
-                if (layout[0] == 0 && layout[1] == 1 && layout[2] == 2 && layout[3] == 3 && layout[4] == 4)
+                if (layout[0] == 0 && layout[1] == 1 && layout[2] == 2 && layout[3] == 3 && layout[4] == 4) {
                     return construct_with_layout<::gridtools::layout_map<0, 1, 2, 3, 4>>(shape[0], shape[1], shape[2],
                                                                                          shape[3], shape[4]);
-                if (layout[0] == 4 && layout[1] == 3 && layout[2] == 2 && layout[3] == 1 && layout[4] == 0)
+                }
+                if (layout[0] == 4 && layout[1] == 3 && layout[2] == 2 && layout[3] == 1 && layout[4] == 0) {
                     return construct_with_layout<::gridtools::layout_map<4, 3, 2, 1, 0>>(shape[0], shape[1], shape[2],
                                                                                          shape[3], shape[4]);
+                }
             }
             default: {
                 std::stringstream err;
-                if (shape.size() > 5)
+                if (shape.size() > 5) {
                     err << "shape not recognized";
+                }
                 else {
                     err << "Layout < ";
-                    for (size_t j = 0; j < layout.size(); ++j)
+                    for (size_t j = 0; j < layout.size(); ++j) {
                         err << layout[j] << " ";
+                    }
                     err << "> not implemented in Atlas.";
                 }
                 throw_Exception(err.str(), Here());
@@ -229,8 +243,9 @@ class ArrayT_impl {
             throw_Exception(err.str(), Here());
         }
 
-        if (array_.valid())
+        if (array_.valid()) {
             array_.syncHostDevice();
+        }
 
         Array* resized = Array::create<Value>(ArrayShape{(idx_t)c...});
 
@@ -397,7 +412,7 @@ Array* Array::create(DataType datatype, ArraySpec&& spec) {
 
 //------------------------------------------------------------------------------
 
-Array::~Array() {}
+Array::~Array() = default;
 
 //------------------------------------------------------------------------------
 
@@ -429,8 +444,9 @@ void ArrayT<Value>::insert(idx_t idx1, idx_t size1) {
     // if( hostNeedsUpdate() ) {
     //    updateHost();
     //}
-    if (not hasDefaultLayout())
+    if (not hasDefaultLayout()) {
         ATLAS_NOTIMPLEMENTED;
+    }
 
     ArrayShape nshape = shape();
     if (idx1 > nshape[0]) {