Add a conservative limiter for the 'cons2' interpolation

[sbrdar]

Description

This addresses the issue #347 to add mass-conserving limiter for over-/under-shoots after the 2nd order conservative interpolation 'cons2'.

This implementation, albeit MPI-parallel, is limited to cell-centred data and assumes that spherical polygons are indexed as source cells (not true when masked cells are present).

Contributor Declaration

By opening this pull request, I affirm the following:

• All authors agree to the Contributor License Agreement.
• The code follows the project's coding standards.
• I have performed self-review and added comments where needed.
• I have added or updated tests to verify that my changes are effective and functional.
• I have run all existing tests and confirmed they pass.

💣💥☠️ Static Analyzer Report ☠️💥💣
https://sites.ecmwf.int/docs/atlas/static-analyzer/PR-350

[Copilot]

Pull request overview

This pull request adds a conservative mass-conserving limiter for the 2nd order conservative interpolation ('cons2') to address over-/under-shoots near strong gradients, as described in issue #347. The implementation includes:

Changes:

• Implements a new ConservativeSphericalPolygonInterpolationLimiter class that applies monotone limiting to prevent unphysical over/undershoots in 2nd order conservative interpolation
• Re-enables and implements the matrix-free 2nd order interpolation method that was previously disabled
• Adds two new analytical test functions: XStep and SlottedCylinder for testing the limiter
• Updates metadata key naming for consistency (e.g., 'errors.conservation' instead of 'errors.conservation_error')
• Adds limiter configuration options with support for 'none', 'zeroslope', and 'clip' strategies

Reviewed changes

Copilot reviewed 12 out of 12 changed files in this pull request and generated 5 comments.

Show a summary per file

File	Description
src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h	Adds limiter support, friend declarations, src_iparam_ for limiter, mass statistics, and gradient computation methods
src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc	Implements matrix-free 2nd order interpolation for all data type combinations, gradient computation, limiter integration, and updates metadata keys
src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolationLimiter.h	Defines the limiter class with detection and correction methods for over/undershoots
src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolationLimiter.cc	Implements the limiter with MPI-parallel support for detecting violations and applying conservative corrections
src/atlas/util/function/XStep_function.h/cc	Adds XStep analytical function from MDPI benchmarking paper for testing
src/atlas/util/function/SlottedCylinder_function.h/cc	Adds SlottedCylinder function from Zalesak 1979 paper for testing
src/tests/interpolation/test_interpolation_conservative.cc	Re-enables 2nd order matrix-free test and updates metadata key names
src/sandbox/interpolation/atlas-conservative-interpolation.cc	Adds command-line options for limiter configuration and new test functions
src/atlas/CMakeLists.txt	Adds new source files and removes SparseMatrixTriplet.h
src/atlas/interpolation/AssembleGlobalMatrix.cc	Removes blank lines (formatting only)

Comments suppressed due to low confidence (11)

src/tests/interpolation/test_interpolation_conservative.cc:165

• The variable 'err' is initialized to -1. on line 165, but this initialization value has no clear semantic meaning. If the metadata.get() calls on lines 167, 183, or 186 fail to retrieve values, 'err' will remain at -1., which could lead to confusing test failures. Consider using a more descriptive initialization or checking if the get operations succeed before using the value.

    double err = -1.;

src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc:1883

• The variable dual_area_inv is declared as volatile on line 1883, which prevents compiler optimizations and forces memory access on every read/write. The use of volatile is typically incorrect for multithreading synchronization (should use atomics or mutexes instead) and appears unnecessary here since this is not shared state or hardware-mapped memory. The volatile qualifier should be removed unless there's a specific documented reason for its use.

    volatile double dual_area_inv     = 0.;

src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc:1907

• Similar to line 1883, the variable dual_area_inv is declared as volatile on line 1907, which is likely unnecessary and can hurt performance. Consider removing the volatile qualifier unless there's a specific documented reason for its use.

    volatile double dual_area_inv     = 0.;

src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolationLimiter.h:50

• The function signature on line 50 takes a parameter named 'scell' but based on the usage in ConservativeSphericalPolygonInterpolationLimiter.cc line 216, this function is called with 'scsp_id'. This naming inconsistency between the declaration and the calling code is confusing. The parameter should be renamed to 'scsp_id' to match how it's being used, or the calling code should be clarified.

    void limit_contrib_from_source(idx_t scell, const Field& src_field, array::ArrayView<double,1>& tgt_lim_vals);

src/atlas/util/function/XStep_function.h:26

• The documentation states "The longitude (lon) and latitude (lat) are assumed to be in degrees" but should clarify this is referring to the function parameters, not the internal calculations. Consider rewording to "The longitude (lon) and latitude (lat) parameters are in degrees" for clarity.

///           The longitude (lon) and latitude (lat) are assumed to be in degrees,

src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h:196

• The Statistics struct has been extended with a new 'mass' array (line 196) to track MASS_SRC, MASS_TGT, and MASS_LIMITER values. However, the implementation doesn't show initialization of this array in the constructor. While I cannot see the constructor in the changed regions, please ensure that the Statistics constructor initializes the mass array, similar to how counts, errors, memory, and time arrays are initialized.

        std::array<double, MASS_ENUM_SIZE> mass;

src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolationLimiter.cc:155

• The TODO comment "convert scell to scsp_id" indicates incomplete implementation. The code assumes scell equals scsp_id, but based on the PR description stating "assumes that spherical polygons are indexed as source cells (not true when masked cells are present)", this assumption may not hold. This is a known limitation that should be resolved before merging, or at minimum should include a runtime check to detect when this assumption is violated and fail gracefully.

                idx_t scsp_id = scell;  // TODO: convert scell to scsp_id

src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc:2088

• In the "CELL TO NODE" case (line 2067), the outer loop iterates over tcsp_id (line 2070), but the inner loop also uses a variable named tcsp_id (line 2073) that shadows the outer loop variable. This shadowing makes the code confusing and error-prone. The inner tcsp_id from tgt_node2csp[tnode] is different from the outer tcsp_id. Consider renaming the inner loop variable to something like tcsp_id_for_node to avoid confusion.

                for (idx_t tcsp_id = 0; tcsp_id < data_->tgt_.csp_size; ++tcsp_id) {
                    idx_t tnode = tgt_csp2node[tcsp_id];
                    double tgt_val = 0.;
                    for (const auto& tcsp_id: tgt_node2csp[tnode]) {
                        const auto& iparam  = tgt_iparam[tcsp_id];
                        for (idx_t i_scsp = 0; i_scsp < iparam.csp_ids.size(); ++i_scsp) {
                            idx_t scsp_id = iparam.csp_ids[i_scsp];
                            idx_t scell   = csp_to_cell(scsp_id, data_->src_);
                            const PointXYZ& src_barycentre = src_points[scell]; // TODO: this is a bad barycentre numerically
                            PointXYZ grad  = src_grads[scell];
                            grad           = grad - PointXYZ::mul(src_barycentre, PointXYZ::dot(grad, src_barycentre));
                            tgt_val += iparam.weights[i_scsp] * (src_vals(scell) + PointXYZ::dot(grad, iparam.centroids[i_scsp] - src_barycentre));
                        }
                    }
                    if (tgt_areas[tnode] > 0.) {
                        tgt_val /= tgt_areas[tnode];
                    }
                    tgt_vals(tnode) = tgt_val;
                }

src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc:2136

• Similar to the CELL TO NODE case, this NODE TO NODE case has variable shadowing on line 2121. The outer loop uses tcsp_id (line 2118), and the inner loop also uses tcsp_id (line 2121). The inner tcsp_id from tgt_node2csp[tnode] shadows the outer loop variable. Consider renaming the inner loop variable to avoid confusion and potential bugs.

                for (idx_t tcsp_id = 0; tcsp_id < data_->tgt_.csp_size; ++tcsp_id) {
                    idx_t tnode = tgt_csp2node[tcsp_id];
                    double tgt_val = 0.;
                    for (const auto& tcsp_id: tgt_node2csp[tnode]) {
                        const auto& iparam = tgt_iparam[tcsp_id];
                        for (idx_t i_scsp = 0; i_scsp < iparam.csp_ids.size(); ++i_scsp) {
                            idx_t scsp_id = iparam.csp_ids[i_scsp];
                            idx_t snode   = src_csp2node[scsp_id];
                            const PointXYZ& src_barycentre = src_points[snode]; // TODO: this is a bad barycentre numerically
                            PointXYZ grad  = src_grads[snode];
                            grad           = grad - PointXYZ::mul(src_barycentre, PointXYZ::dot(grad, src_barycentre));
                            tgt_val += iparam.weights[i_scsp] * (src_vals(snode) + PointXYZ::dot(grad, iparam.centroids[i_scsp] - src_barycentre));
                        }
                    }
                    if (tgt_areas[tnode] > 0.) {
                        tgt_val /= tgt_areas[tnode];
                    }
                    tgt_vals(tnode) = tgt_val;
                }

src/atlas/util/function/SlottedCylinder_function.h:26

• The documentation states "The longitude (lon) and latitude (lat) are assumed to be in radians" but the implementation in SlottedCylinder_function.cc converts from degrees to radians using Constants::degreesToRadians(). The documentation should state that lon and lat are assumed to be in degrees, not radians.

///           The longitude (lon) and latitude (lat) are assumed to be in radians.

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

• The help text for the 'init' option has been updated to list available options in alphabetical order with square brackets notation. However, the list on line 112 shows "slotted_cylinder, solid_body_rotation_wind_magnitude, spherical_harmonic, vortex_rollup (default), xstep" which is not fully alphabetical (xstep should come before vortex_rollup if following strict alphabetical order, or the options should be kept in some other logical grouping). Consider either making it fully alphabetical or adding a comment explaining the ordering convention.

            "init", "Setup initial source field [constant, slotted_cylinder, solid_body_rotation_wind_magnitude, spherical_harmonic, vortex_rollup (default), xstep ]"));

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[codecov-commenter]

Codecov Report

❌ Patch coverage is 58.88158% with 125 lines in your changes missing coverage. Please review.
✅ Project coverage is 79.13%. Comparing base (57639cb) to head (8a6ea84).
⚠️ Report is 12 commits behind head on develop.

Files with missing lines	Patch %	Lines
...onservativeSphericalPolygonInterpolationLimiter.cc	12.28%	100 Missing ⚠️
src/atlas/util/function/XStep_function.cc	0.00%	10 Missing ⚠️
...rc/atlas/util/function/SlottedCylinder_function.cc	0.00%	9 Missing ⚠️
...tured/ConservativeSphericalPolygonInterpolation.cc	96.83%	5 Missing ⚠️
...ConservativeSphericalPolygonInterpolationLimiter.h	0.00%	1 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           develop     #350      +/-   ##
===========================================
- Coverage    79.69%   79.13%   -0.56%     
===========================================
  Files          912      918       +6     
  Lines        71413    62525    -8888     
===========================================
- Hits         56910    49482    -7428     
+ Misses       14503    13043    -1460     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[wdeconinck]

CI is failing. Can you replace the EXPECTS that fail with EXPECT_APPROX_EQ ? Then we see in a blink what the values are. It is not evident for me which values from the log are being used.

[Copilot]

Pull request overview

Copilot reviewed 13 out of 14 changed files in this pull request and generated 3 comments.

Comments suppressed due to low confidence (5)

src/tests/interpolation/test_interpolation_conservative.cc:258

• Two statements ended up on the same line: EXPECT(...);remap_stat[...].... This is easy to miss in reviews and makes future diffs harder to read; please split into separate lines.

    Log::info() << "\n2nd order conservation (new < ref)                   =  (" << std::abs(err_mv) << " < " << tol[5] << ")" << std::endl;
    EXPECT(std::abs(err_mv) < tol[5]);remap_stat[RemapStats::CONS2_MFREE].get("errors.conservation", err_mf);
    Log::info() << "2nd order matrix-free conservation (new < ref)       =  (" << std::abs(err_mf) << " < " << tol[5] << ")" << std::endl;

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

• The json output uses config.getBool("matrix-free", false) but the option defined/used elsewhere in this tool is matrix_free. Unless LocalConfiguration normalizes -/_ for lookups, this will always report false even when --matrix_free is set. Use the same key (matrix_free) consistently.

        output.set("setup.interpolation.limiter", config.getString("limiter", "none"));
        output.set("setup.interpolation.validate", config.getBool("validate", false));
        output.set("setup.interpolation.matrix_free", config.getBool("matrix-free", false));
        output.set("setup.init", config.getString("init", "vortex_rollup"));

src/atlas/util/function/SlottedCylinder.h:25

• Typo in the paper title: "Algorthim" → "Algorithm".

/// \detailed The formula is found in
///           "Fully Multidimensional Flux-Corrected Transport Algorthim for Fluids"
///           by Steven T. Zalesak, JCP 1979
///           as given in their Fig. 11

src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolationLimiter.cc:16

• This implementation relies on standard library types/functions (std::set, std::unordered_map, std::numeric_limits, std::accumulate, std::find, getenv/std::atof) without including the corresponding headers. This can break compilation depending on transitive includes. Add the missing standard headers (e.g. <set>, <unordered_map>, <limits>, <numeric>, <algorithm>, <cstdlib>).

#include <fstream>
#include <iomanip>
#include <unordered_set>
#include <vector>

#include "ConservativeSphericalPolygonInterpolationLimiter.h"

src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc:2533

• Data::footprint() adds memory_of(src_.csp_cell_index) twice (second time on line 2532), and never accounts for tgt_.csp_cell_index. This makes the reported footprint incorrect.

    mem_total += memory_of(src_.csp_index);
    mem_total += memory_of(src_.csp_cell_index);
    mem_total += memory_of(tgt_.csp_index);
    mem_total += memory_of(src_.csp_cell_index);
    return mem_total;

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

send_marked_scells is constructed with an initial size of send_marked_scells_set.size() and then populated via emplace_back, which leaves the first N entries default-initialized (likely 0) and appends N real entries. This will send incorrect global indices in allGatherv and can trigger invalid limiter corrections. Build the vector by reserving and push_back-ing (or fill by index) so that only the intended IDs are sent.

Suggested change

	std::vector<gidx_t> send_marked_scells(send_marked_scells_set.size());
	const auto src_global_index = array::make_view<gidx_t, 1>(src_fs_.global_index());
	for (auto idx : send_marked_scells_set) {
	gidx_t gidx = src_global_index(idx);
	send_marked_scells.emplace_back(gidx);
	std::vector<gidx_t> send_marked_scells;
	send_marked_scells.reserve(send_marked_scells_set.size());
	const auto src_global_index = array::make_view<gidx_t, 1>(src_fs_.global_index());
	for (auto idx : send_marked_scells_set) {
	gidx_t gidx = src_global_index(idx);
	send_marked_scells.push_back(gidx);

[Copilot]

recv_loc_scells[gid] uses operator[], which will insert a default entry when gid is not present on this rank. Since allGatherv collects global indices from all ranks, many gids will not be in the local map, and this will incorrectly map them to local index 0 and apply limiter corrections to the wrong source cell. Use find()/contains() and skip non-local IDs (or exchange local indices instead of global IDs).

Suggested change

	idx_t scell = recv_loc_scells[recv_marked_scells_buf.buffer[i_gid]];
	const gidx_t gid = recv_marked_scells_buf.buffer[i_gid];
	auto it = recv_loc_scells.find(gid);
	if (it == recv_loc_scells.end()) {
	// Global ID not present on this rank; skip non-local source cells
	continue;
	}
	idx_t scell = it->second;

[Copilot]

tgt_mass is computed from tgt_field after the limiter has already modified it, but then tgt_mass += limiter_mass_change applies the limiter correction a second time. This will distort the conservation error and the reported mass.tgt_after_limiter when the limiter is enabled. Also, inv_src_mass is computed before the allReduce, so the relative conservation percentage will be wrong in MPI runs; compute it after reducing src_mass.

Suggested change

	tgt_mass += limiter_mass_change;
	double inv_src_mass = 1.;
	if (src_mass > 0.) {
	inv_src_mass = 1. / src_mass;
	}
	ATLAS_TRACE_MPI(ALLREDUCE) { mpi::comm().allReduceInPlace(&src_mass, 1, eckit::mpi::sum()); }
	ATLAS_TRACE_MPI(ALLREDUCE) { mpi::comm().allReduceInPlace(&tgt_mass, 1, eckit::mpi::sum()); }
	ATLAS_TRACE_MPI(ALLREDUCE) { mpi::comm().allReduceInPlace(&src_mass, 1, eckit::mpi::sum()); }
	ATLAS_TRACE_MPI(ALLREDUCE) { mpi::comm().allReduceInPlace(&tgt_mass, 1, eckit::mpi::sum()); }
	double inv_src_mass = 1.;
	if (src_mass > 0.) {
	inv_src_mass = 1. / src_mass;
	}