From dc486d8eca26196cd6f1891fa85c142cd980600c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 18 Mar 2026 06:58:15 +0100 Subject: [PATCH 01/17] New checkIfInside Backporting of fixes in 64c7d9786e80cb4eb83f2d14f45f4ee7342ddbc1 and bc3d8be3abf338c4e98973bd73af01272b7b736f. Here, we use the cross product to determine if a point is inside or not. If there is uncertainty about this, the caller must decide (i.e. if the point is too close to the boundary to reliably tell). We assume an uncertain point is inside if the caller is checking if a grid corner is inside a triangle, but outside if the caller is checking if a triangle vertex is inside a grid cell. --- .../bamg/include/ConservativeRemapping.hpp | 2 +- contrib/bamg/src/ConservativeRemapping.cpp | 77 ++++++++++++------- 2 files changed, 49 insertions(+), 30 deletions(-) diff --git a/contrib/bamg/include/ConservativeRemapping.hpp b/contrib/bamg/include/ConservativeRemapping.hpp index ec158dbd8..6f47b1731 100644 --- a/contrib/bamg/include/ConservativeRemapping.hpp +++ b/contrib/bamg/include/ConservativeRemapping.hpp @@ -39,7 +39,7 @@ inline void checkTriangle(BamgMesh* bamgmesh, std::vector const &gridCor inline bool visited(int current_triangle, std::vector const &triangles); // Check if points are inside polygon -inline bool checkIfInside(std::vector const &vertx, std::vector const &verty, double testx, double testy); +inline bool checkIfInside(const std::vector& vertx, const std::vector& verty, double testx, double testy, bool inclusive); // Check for intersection and add intersecting points to the list inline bool checkIfIntersecting(double X, double Y, double Xnext, double Ynext, std::vector const &gridCornerX, std::vector const &gridCornerY, // inputs diff --git a/contrib/bamg/src/ConservativeRemapping.cpp b/contrib/bamg/src/ConservativeRemapping.cpp index edf4d701e..e7bc33cea 100644 --- a/contrib/bamg/src/ConservativeRemapping.cpp +++ b/contrib/bamg/src/ConservativeRemapping.cpp @@ -362,7 +362,7 @@ inline void checkTriangle(BamgMesh* bamgmesh, std::vector const &gridCor Y[i] = bamgmesh->Vertices[3*nodeID[i]+1]; // If we're inside we note the point and call self for the surrounding triangles - inCell[i] = checkIfInside(gridCornerX, gridCornerY, X[i], Y[i]); + inCell[i] = checkIfInside(gridCornerX, gridCornerY, X[i], Y[i], false); if ( inCell[i] ) { points.push_back(std::make_pair(X[i],Y[i])); @@ -391,7 +391,7 @@ inline void checkTriangle(BamgMesh* bamgmesh, std::vector const &gridCor int counter = 0; for (int i=0; i const &triangles) return false; } -/* - * Check if points are inside polygon - * Short and works for both triangles and quadrangles (or any higher order polygon). - * - * from https://wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html - * - * Copyright (c) 1970-2003, Wm. Randolph Franklin - * - * Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: - * - * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimers. - * 2. Redistributions in binary form must reproduce the above copyright notice in the documentation and/or other materials provided with the distribution. - * 3. The name of W. Randolph Franklin may not be used to endorse or promote products derived from this Software without specific prior written permission. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ -inline bool checkIfInside(std::vector const &vertx, std::vector const &verty, double testx, double testy) + +// Check if points are inside polygon using a cross product +bool checkIfInside(const std::vector& vertx, const std::vector& verty, double testx, double testy, bool inclusive) { // Initilisation and sanity check - bool inside = false; - int nvert = vertx.size(); + const int nvert = vertx.size(); assert(nvert==verty.size()); - int i; - int j=nvert-1; - for (i=0; itesty) != (verty[j]>testy)) && - (testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) ) - inside = !inside; + // Check if the point is on a vertex of the triangle + const double eps = 1.e-3; + for (int i = 0; i < nvert; ++i) + if (std::abs(testx - vertx[i]) < eps && std::abs(testy - verty[i]) < eps) + return inclusive; - return inside; -} + // The cross product + auto cross = [&](double ax, double ay, double bx, double by, double px, double py) { + return (bx - ax)*(py - ay) - (by - ay)*(px - ax); + }; + + // Check the cross product for all sides + bool hasPos = false; + bool hasNeg = false; + bool hasMaybe = false; + const double epsx = 1e-8; + for (int i=0; i epsx) + hasPos = true; + else if (cp < -epsx) + hasNeg = true; + else + hasMaybe = true; + + // If we found both positive and negative cross products, the point is outside. + if (hasPos && hasNeg ) return false; + } + + // If we're uncertain, then we take our cue from the caller + if ( hasMaybe ) return inclusive; + + return true; + +} //checkIfInside // Check for intersection and add intersecting points to the list // https://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect From 4d6d593ad8032c536a535e38f7abe08aba36c607 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 18 Mar 2026 07:07:15 +0100 Subject: [PATCH 02/17] Add uncertainty to checkIfIntersecting A backport of e494fbd150dbf9bed65c3b566ad3876634f26c34, which says: A comparison with zero is not very useful on a finite-precision machine! Using 1e-6 seems to work well. --- contrib/bamg/src/ConservativeRemapping.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/contrib/bamg/src/ConservativeRemapping.cpp b/contrib/bamg/src/ConservativeRemapping.cpp index e7bc33cea..a5ab8ca2d 100644 --- a/contrib/bamg/src/ConservativeRemapping.cpp +++ b/contrib/bamg/src/ConservativeRemapping.cpp @@ -533,8 +533,8 @@ inline bool checkIfIntersecting(double X, double Y, double Xprev, double Yprev, double s2_x = gridCornerX[i] - gridCornerX[prev]; double s2_y = gridCornerY[i] - gridCornerY[prev]; - double det = -s2_x * s1_y + s1_x * s2_y; - if ( det == 0 ) + const double det = -s2_x * s1_y + s1_x * s2_y; + if ( std::abs(det) < 1e-6 ) continue; // The lines are parallel! double rdet = 1./det; From 5721238b14982e01e0a7e6ac5f6a633a2839bec2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 18 Mar 2026 07:08:59 +0100 Subject: [PATCH 03/17] Limit the reach of the intersection algorithm A backport of d391780bda81c13ebc9983a4fb93d6d349a97252, which says: checkIfIntersecting now ignores intersections that are exactly on the vertex. It never should have included those, and doing it like this is very marginally better in terms of total error. --- contrib/bamg/src/ConservativeRemapping.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/contrib/bamg/src/ConservativeRemapping.cpp b/contrib/bamg/src/ConservativeRemapping.cpp index a5ab8ca2d..8febb8b4f 100644 --- a/contrib/bamg/src/ConservativeRemapping.cpp +++ b/contrib/bamg/src/ConservativeRemapping.cpp @@ -542,12 +542,10 @@ inline bool checkIfIntersecting(double X, double Y, double Xprev, double Yprev, double t = ( s2_x * (Yprev - gridCornerY[prev]) - s2_y * (Xprev - gridCornerX[prev])) * rdet; /* - * Here we assume that the case of overlaping points is an - * intersection. It will result in double counting in some cases, but - * not doing it would result in us not catching all the points all the - * time. + * Here we assume that the case of overlaping points is not an + * intersection. The corner point is cought by checkIfInside anyway. */ - if (s >= 0 && s <= 1 && t >= 0 && t <= 1) + if (s > 0. && s < 1. && t > 0. && t < 1.) { // Intersection detected points.push_back(std::make_pair(Xprev + (t * s1_x), Yprev + (t * s1_y))); From d6f6c176543cd09f3bc3d2cb1c1f62ef12414b1e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 18 Mar 2026 07:10:22 +0100 Subject: [PATCH 04/17] Replace sub-optimal push-back in ConservativeRemappingWeights A backport of 90855a133b150566c293f7058eec5d038d9df88a, which says: I was fiddling with it ans saw that I had already done something similar in ConservativeRemappingFromMeshToMesh. It probably won't have much of an effect. Putting it here, since these functions have been massively rewritten anyway. --- contrib/bamg/src/ConservativeRemapping.cpp | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/contrib/bamg/src/ConservativeRemapping.cpp b/contrib/bamg/src/ConservativeRemapping.cpp index 8febb8b4f..da748e953 100644 --- a/contrib/bamg/src/ConservativeRemapping.cpp +++ b/contrib/bamg/src/ConservativeRemapping.cpp @@ -42,15 +42,13 @@ void ConservativeRemappingWeights(BamgMesh* bamgmesh, std::vector &gridX } // Copy the grid information - int grid_size = gridX.size(); + const int grid_size = gridX.size(); assert(grid_size==gridY.size()); - // Reset the size of gridP, triangles, and weights - // Doing this and using push_back is probably sub-optimal, but that's an - // optimisation for another day - gridP.resize(0); - triangles.resize(0); - weights.resize(0); + // Initialise gridP, triangles, and weights + gridP.resize(grid_size); + triangles.resize(grid_size); + weights.resize(grid_size); // Find which element each P-point hits - use -1 for land points double* elnum_out; @@ -72,7 +70,7 @@ void ConservativeRemappingWeights(BamgMesh* bamgmesh, std::vector &gridX if ( i_elnum_out >= 0 ) { // Save the ppoints - gridP.push_back(ppoint); + gridP[ppoint] = ppoint; // coordinates and size of the grid cell std::vector cornerX(4); @@ -94,8 +92,8 @@ void ConservativeRemappingWeights(BamgMesh* bamgmesh, std::vector &gridX checkTriangle(bamgmesh, cornerX, cornerY, i_elnum_out, local_triangles, local_weights); // Save the weights and triangle numbers - triangles.push_back(local_triangles); - weights.push_back(local_weights); + triangles[ppoint] = local_triangles; + weights[ppoint] = local_weights; } } xDelete(elnum_out); From 9442f5b0c51169bbef5f4f693d5ef8e43acd362d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 18 Mar 2026 07:13:54 +0100 Subject: [PATCH 05/17] Throw in a lot of 'const's It's a good thing, and helped catching two very minor bugs. --- contrib/bamg/src/ConservativeRemapping.cpp | 66 +++++++++++----------- 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/contrib/bamg/src/ConservativeRemapping.cpp b/contrib/bamg/src/ConservativeRemapping.cpp index da748e953..f7796ccd5 100644 --- a/contrib/bamg/src/ConservativeRemapping.cpp +++ b/contrib/bamg/src/ConservativeRemapping.cpp @@ -17,7 +17,7 @@ void ConservativeRemappingWeights(BamgMesh* bamgmesh, std::vector &gridX { // ---------- Initialisation ---------- // // Copy the triangle information - int numTriangles = bamgmesh->TrianglesSize[0]; + const int numTriangles = bamgmesh->TrianglesSize[0]; std::vector indexTr(3*numTriangles); std::vector elnum(numTriangles); for (int tr=0; tr &gridX } // Copy the node information - int numNodes = bamgmesh->VerticesSize[0]; + const int numNodes = bamgmesh->VerticesSize[0]; std::vector coordX(numNodes); std::vector coordY(numNodes); for (int id=0; id &gridX for (int ppoint=0; ppoint= 0 ) @@ -114,12 +114,12 @@ void ConservativeRemappingMeshToGrid(double* &interp_out, std::vector &i for (int i=0; i> points(num_corners); for (int corner=0; corner &i // ---------- Initialisation ---------- // // Copy the triangle information of the _old mesh - int numTriangles = bamgmesh_old->TrianglesSize[0]; + const int numTriangles = bamgmesh_old->TrianglesSize[0]; std::vector indexTr(3*numTriangles); std::vector elnum(numTriangles); for (int tr=0; tr &i } // Copy the node information - int numNodes = bamgmesh_old->VerticesSize[0]; + const int numNodes = bamgmesh_old->VerticesSize[0]; std::vector coordX(numNodes); std::vector coordY(numNodes); for (int id=0; id &i // Copy the triangle information of the _new mesh and calculate the barycentre // Keep the nomenclature for a grid (even if it's a mesh) - int grid_size = bamgmesh_new->TrianglesSize[0]; + const int grid_size = bamgmesh_new->TrianglesSize[0]; std::vector gridX(grid_size); std::vector gridY(grid_size); std::vector gridCornerX(3*grid_size); @@ -218,7 +218,7 @@ void ConservativeRemappingMeshToMesh(double* &interp_out, std::vector &i for (int i=0; i<3; ++i) { // grid corner == vertice - int id = bamgmesh_new->Triangles[4*tr+i] - 1; // Here we need C/C++ numbering + const int id = bamgmesh_new->Triangles[4*tr+i] - 1; // Here we need C/C++ numbering gridCornerX[3*tr+i] = bamgmesh_new->Vertices[3*id]; gridCornerY[3*tr+i] = bamgmesh_new->Vertices[3*id+1]; @@ -257,7 +257,7 @@ void ConservativeRemappingMeshToMesh(double* &interp_out, std::vector &i assert( elnum_out[ppoint] >= 0. ); // Carefully take the right integer value for element number - int i_elnum_out = std::round(elnum_out[ppoint]); + const int i_elnum_out = std::round(elnum_out[ppoint]); // Save the ppoint - for compatibility with ConservativeRemappingMeshToGrid gridP[ppoint] = ppoint; @@ -271,7 +271,7 @@ void ConservativeRemappingMeshToMesh(double* &interp_out, std::vector &i std::vector nodes_new(3); for ( int i=0; i<3; ++i ) { - int id = bamgmesh_new->Triangles[4*ppoint+i] - 1; // Here we need C/C++ numbering + const int id = bamgmesh_new->Triangles[4*ppoint+i] - 1; // Here we need C/C++ numbering /* The previous numbering of the nodes on the boundaries given by bamg is wrong. * Bamg adds these nodes at the beginning of the list of nodes. * In our case, these nodes are always the same as the boundary is neither @@ -333,17 +333,17 @@ inline void checkTriangle(BamgMesh* bamgmesh, std::vector const &gridCor /* * We must plant the flag here, before any recursion can happen. We must * also push_back on weights and remember where to write the final weight - * because some of the recursive calles may do a push_back of their own. + * because some of the recursive calls may do a push_back of their own. */ triangles.push_back(current_triangle); weights.push_back(0.); - int loc = weights.size()-1; + const int loc = weights.size()-1; // This is a list of the points we use to calculate the area std::vector> points; - int num_corners = gridCornerX.size(); - assert(num_corners = gridCornerY.size()); + const int num_corners = gridCornerX.size(); + assert(num_corners == gridCornerY.size()); // ---------- Now we do the three checks ---------- // @@ -365,10 +365,10 @@ inline void checkTriangle(BamgMesh* bamgmesh, std::vector const &gridCor { points.push_back(std::make_pair(X[i],Y[i])); - int num_elements = bamgmesh->NodalElementConnectivitySize[1]; - for (int j=0; jNodalElementConnectivitySize[1]; + for (int j = 0; j < num_elements; j++) { - int elt_num = bamgmesh->NodalElementConnectivity[num_elements*nodeID[i]+j] - 1; // Here we need C/C++ numbering + const int elt_num = bamgmesh->NodalElementConnectivity[num_elements*nodeID[i]+j] - 1; // Here we need C/C++ numbering // Negative elt_num means there are no more elements belonging to this node if ( elt_num < 0 ) continue; @@ -417,7 +417,7 @@ inline void checkTriangle(BamgMesh* bamgmesh, std::vector const &gridCor // We don't know which of the (up to) three connected triangles is the right one so we must search for (int j=0; j<3; j++) { - int elt_num = bamgmesh->ElementConnectivity[3*current_triangle+j] - 1; // Here we need C/C++ numbering + const int elt_num = bamgmesh->ElementConnectivity[3*current_triangle+j] - 1; // Here we need C/C++ numbering // Negative elt_num means there are no more elements belonging to this node if ( elt_num < 0 ) continue; @@ -429,7 +429,7 @@ inline void checkTriangle(BamgMesh* bamgmesh, std::vector const &gridCor int counter = 0; for (int k=0; k<3; ++k) { - int myID = bamgmesh->Triangles[4*elt_num+k] - 1; // Here we need C/C++ numbering + const int myID = bamgmesh->Triangles[4*elt_num+k] - 1; // Here we need C/C++ numbering if ( myID==nodeID[i] || myID==nodeID[prev] ) ++counter; } @@ -517,27 +517,27 @@ inline bool checkIfIntersecting(double X, double Y, double Xprev, double Yprev, std::vector> &points) // side-effect { // Initialise - int num_corners = gridCornerX.size(); - assert(num_corners = gridCornerY.size()); + const int num_corners = gridCornerX.size(); + assert(num_corners == gridCornerY.size()); bool ret_val = false; - double s1_x = X - Xprev; - double s1_y = Y - Yprev; + const double s1_x = X - Xprev; + const double s1_y = Y - Yprev; // Loop over the grid int prev=num_corners-1; for (int i=0; i> &points) // Calculate value of shoelace formula double area = 0.; - int n = points.size(); + const int n = points.size(); int j = n-1; for (int i=0; i < n; j=i++) area += (points[j].first + points[i].first) * (points[j].second - points[i].second); // Return half the absolute value return std::abs(area) * 0.5; -} +} // area // Sort the points in a clock wise manner around the centre of the cloud of points inline void sortClockwise(std::vector> &points) { // Calculate the centre point - int n = points.size(); + const int n = points.size(); double centreX = 0.; double centreY = 0.; for (auto it=points.begin(); it!=points.end(); ++it) @@ -586,7 +586,7 @@ inline void sortClockwise(std::vector> &points) centreX += it->first; centreY += it->second; } - double rn = 1./(double)n; + const double rn = 1./double(n); centreX *= rn; centreY *= rn; From afa7b4d8fb2786b2800676f8e6c3d822b5bcbf26 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 18 Mar 2026 07:15:58 +0100 Subject: [PATCH 06/17] Formatting changes This is just to make comparison with the work in PR #684 easier. --- contrib/bamg/src/ConservativeRemapping.cpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/contrib/bamg/src/ConservativeRemapping.cpp b/contrib/bamg/src/ConservativeRemapping.cpp index f7796ccd5..2d84b0ca2 100644 --- a/contrib/bamg/src/ConservativeRemapping.cpp +++ b/contrib/bamg/src/ConservativeRemapping.cpp @@ -20,7 +20,7 @@ void ConservativeRemappingWeights(BamgMesh* bamgmesh, std::vector &gridX const int numTriangles = bamgmesh->TrianglesSize[0]; std::vector indexTr(3*numTriangles); std::vector elnum(numTriangles); - for (int tr=0; trTriangles[4*tr]; @@ -132,7 +132,7 @@ void ConservativeRemappingMeshToGrid(double* &interp_out, std::vector &i interp_out[gridP[i]*nb_var+var] *= r_cell_area; } } -} +} // ConservativeRemappingMeshToGrid // Apply weights going from grid to mesh void ConservativeRemappingGridToMesh(double* &interp_out, std::vector &interp_in, int const nb_var, int const numElements, std::vector &gridP, std::vector> &triangles, std::vector> &weights) @@ -183,7 +183,7 @@ void ConservativeRemappingMeshToMesh(double* &interp_out, std::vector &i const int numTriangles = bamgmesh_old->TrianglesSize[0]; std::vector indexTr(3*numTriangles); std::vector elnum(numTriangles); - for (int tr=0; trTriangles[4*tr]; @@ -211,7 +211,7 @@ void ConservativeRemappingMeshToMesh(double* &interp_out, std::vector &i std::vector gridY(grid_size); std::vector gridCornerX(3*grid_size); std::vector gridCornerY(3*grid_size); - for (int tr=0; tr const &gridCor assert(counter>=1); assert(counter<=2); - if ( counter==2 && !visited(elt_num, triangles) ) + if ( counter == 2 && !visited(elt_num, triangles) ) checkTriangle(bamgmesh, gridCornerX, gridCornerY, elt_num, triangles, weights); } } @@ -447,17 +447,17 @@ inline void checkTriangle(BamgMesh* bamgmesh, std::vector const &gridCor // We've checked everything. Finish by calculating the area weights[loc] = area(points); -} +}//checkTriangle // Check if we've already visited this triangle inline bool visited(int current_triangle, std::vector const &triangles) { - for (auto it=triangles.begin(); it!=triangles.end(); ++it) + for (auto it = triangles.begin(); it != triangles.end(); ++it) if ( current_triangle == *it ) return true; return false; -} +} // visited // Check if points are inside polygon using a cross product @@ -551,7 +551,7 @@ inline bool checkIfIntersecting(double X, double Y, double Xprev, double Yprev, } } return ret_val; -} +} // checkIfIntersecting // Calculate the area inline double area(std::vector> &points) @@ -609,7 +609,7 @@ inline void sortClockwise(std::vector> &points) it->first += centreX; it->second += centreY; } -} +} // sortClockwise // The custom sorting function // https://stackoverflow.com/questions/6989100/sort-points-in-clockwise-order <- the second code exmple @@ -642,4 +642,4 @@ inline bool CWSort(std::pair p1, std::pair p2) return p2.first*p1.second < p2.second*p1.first; else return qa < qb; -} +} // CWSort From 927a4721b7fc2a6caf691353ae26d44a964ef822 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 18 Mar 2026 07:16:39 +0100 Subject: [PATCH 07/17] Add . to doubles Propably not very useful, but helps with comparing with the work in PR --- contrib/bamg/src/ConservativeRemapping.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/contrib/bamg/src/ConservativeRemapping.cpp b/contrib/bamg/src/ConservativeRemapping.cpp index 2d84b0ca2..1c198a566 100644 --- a/contrib/bamg/src/ConservativeRemapping.cpp +++ b/contrib/bamg/src/ConservativeRemapping.cpp @@ -125,7 +125,7 @@ void ConservativeRemappingMeshToGrid(double* &interp_out, std::vector &i for (int var=0; var &i std::vector gridCornerY(3*grid_size); for (int tr = 0; tr < grid_size; ++tr) { - gridX[tr] = 0; // barycentre - gridY[tr] = 0; + gridX[tr] = 0.; // barycentre + gridY[tr] = 0.; for (int i=0; i<3; ++i) { // grid corner == vertice @@ -225,8 +225,8 @@ void ConservativeRemappingMeshToMesh(double* &interp_out, std::vector &i gridX[tr] += gridCornerX[3*tr+i]; gridY[tr] += gridCornerY[3*tr+i]; } - gridX[tr] /= 3; - gridY[tr] /= 3; + gridX[tr] /= 3.; + gridY[tr] /= 3.; } // Initialise gridP, triangles, and weights From ddc13c87bb9e1935a4c6ca12cad442dbb26c4065 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Wed, 18 Mar 2026 07:17:34 +0100 Subject: [PATCH 08/17] Replace two continues with breaks It doesn't make very much difference, but using continue is correct. --- contrib/bamg/src/ConservativeRemapping.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/contrib/bamg/src/ConservativeRemapping.cpp b/contrib/bamg/src/ConservativeRemapping.cpp index 1c198a566..75a5a9cf6 100644 --- a/contrib/bamg/src/ConservativeRemapping.cpp +++ b/contrib/bamg/src/ConservativeRemapping.cpp @@ -370,7 +370,7 @@ inline void checkTriangle(BamgMesh* bamgmesh, std::vector const &gridCor { const int elt_num = bamgmesh->NodalElementConnectivity[num_elements*nodeID[i]+j] - 1; // Here we need C/C++ numbering // Negative elt_num means there are no more elements belonging to this node - if ( elt_num < 0 ) continue; + if ( elt_num < 0 ) break; if ( ! visited(elt_num, triangles) ) checkTriangle(bamgmesh, gridCornerX, gridCornerY, elt_num, triangles, weights); @@ -419,7 +419,7 @@ inline void checkTriangle(BamgMesh* bamgmesh, std::vector const &gridCor { const int elt_num = bamgmesh->ElementConnectivity[3*current_triangle+j] - 1; // Here we need C/C++ numbering // Negative elt_num means there are no more elements belonging to this node - if ( elt_num < 0 ) continue; + if ( elt_num < 0 ) break; /* * Find the triangle adjacent to the current one by comparing From 595eb16f3683a49df014971660232ff656a38577 Mon Sep 17 00:00:00 2001 From: Timothy Williams Date: Thu, 19 Mar 2026 11:40:14 +0100 Subject: [PATCH 09/17] fix typo in contrib/bamg/src/ConservativeRemapping.cpp --- contrib/bamg/src/ConservativeRemapping.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/contrib/bamg/src/ConservativeRemapping.cpp b/contrib/bamg/src/ConservativeRemapping.cpp index 75a5a9cf6..f777ca5ed 100644 --- a/contrib/bamg/src/ConservativeRemapping.cpp +++ b/contrib/bamg/src/ConservativeRemapping.cpp @@ -540,8 +540,8 @@ inline bool checkIfIntersecting(double X, double Y, double Xprev, double Yprev, const double t = ( s2_x * (Yprev - gridCornerY[prev]) - s2_y * (Xprev - gridCornerX[prev])) * rdet; /* - * Here we assume that the case of overlaping points is not an - * intersection. The corner point is cought by checkIfInside anyway. + * Here we assume that the case of overlapping points is not an + * intersection. The corner point is caught by checkIfInside anyway. */ if (s > 0. && s < 1. && t > 0. && t < 1.) { From 2aaee55a41829e22d9fccf406cbe6ee27decb63f Mon Sep 17 00:00:00 2001 From: tdcwilliams Date: Tue, 17 Mar 2026 13:39:05 +0100 Subject: [PATCH 10/17] Get helpful debugging flags from develop: change NEXTSIM_BUILD_TYPE: now can be VALGRIND, DEBUG, PROFILE or not set. VALGRIND has no optimisation; DEBUG has some optimisation and uses flags to give line numbers for memory errors/segmentation faults. PROFILE is the same as before. It can't be lower case or mixed case anymore. --- contrib/bamg/src/Makefile | 10 +++++++++- contrib/mapx/src/Makefile | 10 +++++++++- core/src/Makefile | 10 +++++++++- model/Makefile | 8 +++++++- 4 files changed, 34 insertions(+), 4 deletions(-) diff --git a/contrib/bamg/src/Makefile b/contrib/bamg/src/Makefile index 3cf4c7928..0904b23c3 100644 --- a/contrib/bamg/src/Makefile +++ b/contrib/bamg/src/Makefile @@ -17,7 +17,15 @@ CXXFLAGS += -std=c++14 CXXFLAGS += -pedantic -ftemplate-depth-256 -Wno-inline -D_MULTITHREADING_ -ifneq (,$(strip $(filter DEBUG Debug debug,$(NEXTSIM_BUILD_TYPE)))) +ifneq (,$(strip $(filter DEBUG,$(NEXTSIM_BUILD_TYPE)))) + CXXFLAGS += -g -O1 -fsanitize=address -fno-omit-frame-pointer -DNDEBUG + LDFLAGS += -fsanitize=address -fno-omit-frame-pointer +ifneq ($(KERNEL),Linux) + CXXFLAGS += -Wl,-no_pie +endif +endif + +ifneq (,$(strip $(filter VALGRIND,$(NEXTSIM_BUILD_TYPE)))) CXXFLAGS += -g -O0 -DNDEBUG ifneq ($(KERNEL),Linux) CXXFLAGS += -Wl,-no_pie diff --git a/contrib/mapx/src/Makefile b/contrib/mapx/src/Makefile index 14c2cf5ba..f5f9d7bce 100644 --- a/contrib/mapx/src/Makefile +++ b/contrib/mapx/src/Makefile @@ -18,7 +18,15 @@ ifdef NEXTSIM_COMPILE_VERBOSE CXXFLAGS += -v endif -ifneq (,$(strip $(filter DEBUG Debug debug,$(NEXTSIM_BUILD_TYPE)))) +ifneq (,$(strip $(filter DEBUG,$(NEXTSIM_BUILD_TYPE)))) + CCFLAGS += -g -O1 -fsanitize=address -fno-omit-frame-pointer -DNDEBUG + LDFLAGS += -fsanitize=address -fno-omit-frame-pointer +ifneq ($(KERNEL),Linux) + CCFLAGS += -Wl,-no_pie +endif +endif + +ifneq (,$(strip $(filter VALGRIND,$(NEXTSIM_BUILD_TYPE)))) CCFLAGS += -g -O0 -DNDEBUG ifneq ($(KERNEL),Linux) CCFLAGS += -Wl,-no_pie diff --git a/core/src/Makefile b/core/src/Makefile index 4ad3f889f..ed77ee034 100644 --- a/core/src/Makefile +++ b/core/src/Makefile @@ -25,7 +25,15 @@ ifdef USE_OASIS CHAN = MPI1 endif -ifneq (,$(strip $(filter DEBUG Debug debug,$(NEXTSIM_BUILD_TYPE)))) +ifneq (,$(strip $(filter DEBUG,$(NEXTSIM_BUILD_TYPE)))) + CXXFLAGS += -g -O1 -fsanitize=address -fno-omit-frame-pointer -DNDEBUG + LDFLAGS += -fsanitize=address -fno-omit-frame-pointer +ifneq ($(KERNEL),Linux) + CXXFLAGS += -Wl,-no_pie +endif +endif + +ifneq (,$(strip $(filter VALGRIND,$(NEXTSIM_BUILD_TYPE)))) CXXFLAGS += -g -O0 -DNDEBUG ifneq ($(KERNEL),Linux) CXXFLAGS += -Wl,-no_pie diff --git a/model/Makefile b/model/Makefile index e5dc418bd..5fdc209e7 100644 --- a/model/Makefile +++ b/model/Makefile @@ -19,7 +19,13 @@ ifdef USE_OASIS LDFLAGS += -L$(OASIS_DIR)/lib -lpsmile.$(CHAN) -lmct -lmpeu -lscrip $(LD_EXTRA_OASIS) endif -ifneq (,$(strip $(filter DEBUG Debug debug PROFILE Profile profile,$(NEXTSIM_BUILD_TYPE)))) +ifneq (,$(strip $(filter DEBUG,$(NEXTSIM_BUILD_TYPE)))) + CXXFLAGS := $(filter-out -O3 -pthread,$(CXXFLAGS)) + CXXFLAGS += -g -O1 -fsanitize=address -fno-omit-frame-pointer -DNDEBUG + LDFLAGS += -fsanitize=address -fno-omit-frame-pointer +endif + +ifneq (,$(strip $(filter VALGRIND PROFILE,$(NEXTSIM_BUILD_TYPE)))) CXXFLAGS := $(filter-out -O3 -pthread,$(CXXFLAGS)) CXXFLAGS += -g -O0 -DNDEBUG ifneq (,$(strip $(filter PROFILE Profile profile,$(NEXTSIM_BUILD_TYPE)))) From 2d88559f932e4275cb3fa05713b54673373dfc17 Mon Sep 17 00:00:00 2001 From: tdcwilliams Date: Thu, 19 Mar 2026 12:19:22 +0100 Subject: [PATCH 11/17] compiler flags to enable tracking of nan errors. Fix use of -DNDEBUG (gcc thing meaning "not debugging", which disables assert's) --- contrib/bamg/src/Makefile | 10 ++++++---- contrib/mapx/src/Makefile | 12 +++++++----- core/src/Makefile | 10 ++++++---- model/Makefile | 15 +++++++-------- 4 files changed, 26 insertions(+), 21 deletions(-) diff --git a/contrib/bamg/src/Makefile b/contrib/bamg/src/Makefile index 0904b23c3..5e1520452 100644 --- a/contrib/bamg/src/Makefile +++ b/contrib/bamg/src/Makefile @@ -15,18 +15,20 @@ CXXFLAGS := $(filter-out -fopenmp -pthread,$(CXXFLAGS)) CXXFLAGS += -std=c++14 -CXXFLAGS += -pedantic -ftemplate-depth-256 -Wno-inline -D_MULTITHREADING_ +CXXFLAGS += -pedantic -ftemplate-depth-256 -Wno-inline -D_MULTITHREADING_ -DNDEBUG ifneq (,$(strip $(filter DEBUG,$(NEXTSIM_BUILD_TYPE)))) - CXXFLAGS += -g -O1 -fsanitize=address -fno-omit-frame-pointer -DNDEBUG - LDFLAGS += -fsanitize=address -fno-omit-frame-pointer + CXXFLAGS := $(filter-out -O3 -pthread -DNDEBUG,$(CXXFLAGS)) + CXXFLAGS += -g -O1 + CXXFLAGS += -fsanitize=address,float-divide-by-zero -fno-omit-frame-pointer -fno-finite-math-only ifneq ($(KERNEL),Linux) CXXFLAGS += -Wl,-no_pie endif endif ifneq (,$(strip $(filter VALGRIND,$(NEXTSIM_BUILD_TYPE)))) - CXXFLAGS += -g -O0 -DNDEBUG + CXXFLAGS := $(filter-out -O3 -pthread -DNDEBUG,$(CXXFLAGS)) + CXXFLAGS += -g -O0 ifneq ($(KERNEL),Linux) CXXFLAGS += -Wl,-no_pie endif diff --git a/contrib/mapx/src/Makefile b/contrib/mapx/src/Makefile index f5f9d7bce..3e3840cf4 100644 --- a/contrib/mapx/src/Makefile +++ b/contrib/mapx/src/Makefile @@ -13,21 +13,23 @@ LIBRARYDIR=$(NEXTSIMDIR)/lib/ # add gcc option flags -CCFLAGS += -pedantic -Wno-inline +CCFLAGS += -pedantic -Wno-inline -DNDEBUG ifdef NEXTSIM_COMPILE_VERBOSE - CXXFLAGS += -v + CCFLAGS += -v endif ifneq (,$(strip $(filter DEBUG,$(NEXTSIM_BUILD_TYPE)))) - CCFLAGS += -g -O1 -fsanitize=address -fno-omit-frame-pointer -DNDEBUG - LDFLAGS += -fsanitize=address -fno-omit-frame-pointer + CCFLAGS := $(filter-out -O3 -pthread -DNDEBUG,$(CCFLAGS)) + CCFLAGS += -g -O1 -DDEBUGGING + CCFLAGS += -fsanitize=address,float-divide-by-zero -fno-omit-frame-pointer -fno-finite-math-only ifneq ($(KERNEL),Linux) CCFLAGS += -Wl,-no_pie endif endif ifneq (,$(strip $(filter VALGRIND,$(NEXTSIM_BUILD_TYPE)))) - CCFLAGS += -g -O0 -DNDEBUG + CCFLAGS := $(filter-out -O3 -pthread -DNDEBUG,$(CCFLAGS)) + CCFLAGS += -g -O0 -DDEBUGGING ifneq ($(KERNEL),Linux) CCFLAGS += -Wl,-no_pie endif diff --git a/core/src/Makefile b/core/src/Makefile index ed77ee034..4897ed911 100644 --- a/core/src/Makefile +++ b/core/src/Makefile @@ -16,7 +16,7 @@ LIBRARYDIR=$(NEXTSIMDIR)/lib/ CXXFLAGS += -std=c++14 # add g++ option flags -CXXFLAGS += -pedantic -ftemplate-depth-256 -Wno-inline -fpermissive +CXXFLAGS += -pedantic -ftemplate-depth-256 -Wno-inline -fpermissive -DNDEBUG ifdef USE_OASIS CXXFLAGS += -DOASIS @@ -26,15 +26,17 @@ ifdef USE_OASIS endif ifneq (,$(strip $(filter DEBUG,$(NEXTSIM_BUILD_TYPE)))) - CXXFLAGS += -g -O1 -fsanitize=address -fno-omit-frame-pointer -DNDEBUG - LDFLAGS += -fsanitize=address -fno-omit-frame-pointer + CXXFLAGS := $(filter-out -O3 -pthread -DNDEBUG,$(CXXFLAGS)) + CXXFLAGS += -g -O1 + CXXFLAGS += -fsanitize=address,float-divide-by-zero -fno-omit-frame-pointer -fno-finite-math-only ifneq ($(KERNEL),Linux) CXXFLAGS += -Wl,-no_pie endif endif ifneq (,$(strip $(filter VALGRIND,$(NEXTSIM_BUILD_TYPE)))) - CXXFLAGS += -g -O0 -DNDEBUG + CXXFLAGS := $(filter-out -O3 -pthread -DNDEBUG,$(CXXFLAGS)) + CXXFLAGS += -g -O0 ifneq ($(KERNEL),Linux) CXXFLAGS += -Wl,-no_pie endif diff --git a/model/Makefile b/model/Makefile index 5fdc209e7..3f0d44376 100644 --- a/model/Makefile +++ b/model/Makefile @@ -5,8 +5,7 @@ PROGNAME=nextsim.exec CXXFLAGS += -std=c++14 # add g++ option flags -CXXFLAGS += -ftemplate-depth-256 -Wno-inline \ - -DHAVE_CONFIG_H -D_MULTITHREADING_ +CXXFLAGS += -ftemplate-depth-256 -Wno-inline -DHAVE_CONFIG_H -D_MULTITHREADING_ -DNDEBUG ifdef NEXTSIM_COMPILE_VERBOSE CXXFLAGS += -v endif @@ -20,15 +19,15 @@ ifdef USE_OASIS endif ifneq (,$(strip $(filter DEBUG,$(NEXTSIM_BUILD_TYPE)))) - CXXFLAGS := $(filter-out -O3 -pthread,$(CXXFLAGS)) - CXXFLAGS += -g -O1 -fsanitize=address -fno-omit-frame-pointer -DNDEBUG - LDFLAGS += -fsanitize=address -fno-omit-frame-pointer + CXXFLAGS := $(filter-out -O3 -pthread -DNDEBUG,$(CXXFLAGS)) + CXXFLAGS += -g -O1 + CXXFLAGS += -fsanitize=address,float-divide-by-zero -fno-omit-frame-pointer -fno-finite-math-only endif ifneq (,$(strip $(filter VALGRIND PROFILE,$(NEXTSIM_BUILD_TYPE)))) - CXXFLAGS := $(filter-out -O3 -pthread,$(CXXFLAGS)) - CXXFLAGS += -g -O0 -DNDEBUG -ifneq (,$(strip $(filter PROFILE Profile profile,$(NEXTSIM_BUILD_TYPE)))) + CXXFLAGS := $(filter-out -O3 -pthread -DNDEBUG,$(CXXFLAGS)) + CXXFLAGS += -g -O0 +ifneq (,$(strip $(filter PROFILE,$(NEXTSIM_BUILD_TYPE)))) CXXFLAGS += -DWITHGPERFTOOLS endif endif From 24b8557e953c05e77bc5042e2f68cd10fe0dfaf0 Mon Sep 17 00:00:00 2001 From: tdcwilliams Date: Fri, 20 Mar 2026 09:30:28 +0100 Subject: [PATCH 12/17] Modify Lemieux basal stress calculation to avoid division by zero. Remove Bouillon one since it is never used. --- model/enums.hpp | 1 - model/finiteelement.cpp | 17 ++++------------- 2 files changed, 4 insertions(+), 14 deletions(-) diff --git a/model/enums.hpp b/model/enums.hpp index a979cf5ca..efe709d25 100644 --- a/model/enums.hpp +++ b/model/enums.hpp @@ -87,7 +87,6 @@ namespace setup { NONE = 0, LEMIEUX = 1, - BOUILLON = 2 }; enum class IceCategoryType diff --git a/model/finiteelement.cpp b/model/finiteelement.cpp index 6a9284665..228bcaef6 100644 --- a/model/finiteelement.cpp +++ b/model/finiteelement.cpp @@ -1369,7 +1369,6 @@ FiniteElement::initOptAndParam() const boost::unordered_map str2basal_stress= boost::assign::map_list_of ("none", setup::BasalStressType::NONE) ("lemieux", setup::BasalStressType::LEMIEUX) - ("bouillon", setup::BasalStressType::BOUILLON); M_basal_stress_type = this->getOptionFromMap("setup.basal_stress-type", str2basal_stress); //! \param M_basal_stress_type (string) Option on the type of basal stress (none, from Lemieux et al., 2016 or from Bouillon) LOG(DEBUG) <<"BASALSTRESTYPE= "<< (int) M_basal_stress_type <<"\n"; @@ -10276,7 +10275,7 @@ FiniteElement::explicitSolve() double max_keel_depth=28; // [m] from "A comprehensive analysis of the morphology of first-year sea ice ridges" double ice_to_keel_factor=19.28; // from "A comprehensive analysis of the morphology of first-year sea ice ridges" - double keel_depth; + double mean_keel_depth; // use (M_conc[cpt] * keel_depth) to avoid division by zero double critical_h; double critical_h_mod; double const min_water_depth = 2.; //m @@ -10290,21 +10289,13 @@ FiniteElement::explicitSolve() critical_h = 0.; critical_h_mod = 0.; break; - case setup::BasalStressType::BOUILLON: - // Sylvain's grounding scheme - // TODO: Remove this one - we've never used it - keel_depth = ice_to_keel_factor * std::sqrt(M_thick[cpt]/M_conc[cpt]); - keel_depth = std::min( keel_depth, max_keel_depth ); - critical_h = M_conc[cpt] * std::pow(depth_eff / ice_to_keel_factor, 2.); - critical_h_mod = M_conc[cpt] * std::pow(keel_depth / ice_to_keel_factor, 2.); - break; case setup::BasalStressType::LEMIEUX: // JF Lemieux's grounding (critical_h = h_c, critical_h_mod = h) // Limit keel depth (JF doesn't do that). - keel_depth = k1 * M_thick[cpt] / M_conc[cpt]; - keel_depth = std::min( keel_depth, max_keel_depth ); + mean_keel_depth = k1 * M_thick[cpt]; + mean_keel_depth = std::min( mean_keel_depth, M_conc[cpt] * max_keel_depth ); critical_h = M_conc[cpt] * depth_eff / k1; - critical_h_mod = M_conc[cpt] * keel_depth / k1; + critical_h_mod = mean_keel_depth / k1; break; } From 79fc27cb372571dab499d6475b54f9c8217625af Mon Sep 17 00:00:00 2001 From: tdcwilliams Date: Fri, 20 Mar 2026 10:19:07 +0100 Subject: [PATCH 13/17] modify some myi code to avoid division by zero --- model/finiteelement.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/model/finiteelement.cpp b/model/finiteelement.cpp index 228bcaef6..eade19854 100644 --- a/model/finiteelement.cpp +++ b/model/finiteelement.cpp @@ -6098,13 +6098,13 @@ FiniteElement::thermo(int dt) } else // on a non-reset day, myi is only modified by melting, not freezing { - // We ignore the young ice for now - double del_c_ratio = std::min(M_conc[i]/old_conc,1.); - double del_v_ratio = std::min(M_thick[i]/old_vol,1.); - if (del_v_ratio < 1.) // if there is some melt of old ice + if ((M_thick[i] < old_vol) && (old_conc > 0) && (old_vol > 0)) // if there is some melt of old ice { if (equal_melting) { + // We ignore the young ice for now + double const del_c_ratio = std::min(M_conc[i]/old_conc,1.); + double const del_v_ratio = std::min(M_thick[i]/old_vol,1.); del_ci_mlt_myi = std::min(0.,M_conc_myi[i]*(del_c_ratio-1.)); // <0 del_vi_mlt_myi = std::min(0.,M_thick_myi[i]*(del_v_ratio-1.)); // <0 } From 936d2bdab52a31742e7884fbb6d2328b8a9fb63f Mon Sep 17 00:00:00 2001 From: tdcwilliams Date: Fri, 20 Mar 2026 10:20:03 +0100 Subject: [PATCH 14/17] modify thermoWinton code to avoid division by zero --- model/finiteelement.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/model/finiteelement.cpp b/model/finiteelement.cpp index eade19854..e114e9a19 100644 --- a/model/finiteelement.cpp +++ b/model/finiteelement.cpp @@ -6806,7 +6806,7 @@ FiniteElement::thermoWinton(const double dt, const double conc, const double vol double f1 = h1/hi*2.; // Fraction of layer 1 ice found in the new layer 1 double Tbar = f1*( T1 + qi*Tfr_ice/(Crho*T1) ) + (1-f1)*T2; // (39) T1 = ( Tbar - std::sqrt(Tbar*Tbar - 4*Tfr_ice*qi/Crho) )/2.; // (38) - } else { + } else if (hi > 0.) { // Upper layer ice is added to the lower layer // T2 changes, but T1 not double f1 = (2.*h1-hi)/hi; // Fraction of layer 1 ice found in new layer 2 From c192c1a036f7dde61465b82d72ec2c621472298e Mon Sep 17 00:00:00 2001 From: tdcwilliams Date: Fri, 20 Mar 2026 10:20:35 +0100 Subject: [PATCH 15/17] mark some sources of memory leaks --- model/finiteelement.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/model/finiteelement.cpp b/model/finiteelement.cpp index e114e9a19..2dbeeca87 100644 --- a/model/finiteelement.cpp +++ b/model/finiteelement.cpp @@ -991,7 +991,7 @@ FiniteElement::checkReloadMainDatasets(double const CRtime) void FiniteElement::initBamg() { - bamgopt = new BamgOpts(); + bamgopt = new BamgOpts();//TODO memory leak bamgopt->Crack = 0; bamgopt->anisomax = 1e30; @@ -9706,6 +9706,7 @@ FiniteElement::readRestart(std::string const& name_str) std::vector misc_int; std::vector time_vec; + //TODO memory leak in this block - probably from bamg pointers if (M_rank == 0) { From a1d01078fc10d9d9120d3aee85288f5f9ccd12c3 Mon Sep 17 00:00:00 2001 From: tdcwilliams Date: Fri, 20 Mar 2026 10:35:23 +0100 Subject: [PATCH 16/17] fix typo --- model/finiteelement.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/model/finiteelement.cpp b/model/finiteelement.cpp index 2dbeeca87..1bf5024ca 100644 --- a/model/finiteelement.cpp +++ b/model/finiteelement.cpp @@ -1368,7 +1368,7 @@ FiniteElement::initOptAndParam() const boost::unordered_map str2basal_stress= boost::assign::map_list_of ("none", setup::BasalStressType::NONE) - ("lemieux", setup::BasalStressType::LEMIEUX) + ("lemieux", setup::BasalStressType::LEMIEUX); M_basal_stress_type = this->getOptionFromMap("setup.basal_stress-type", str2basal_stress); //! \param M_basal_stress_type (string) Option on the type of basal stress (none, from Lemieux et al., 2016 or from Bouillon) LOG(DEBUG) <<"BASALSTRESTYPE= "<< (int) M_basal_stress_type <<"\n"; From da5351c358a02b28fe1ab256d29887de454c1aa8 Mon Sep 17 00:00:00 2001 From: tdcwilliams Date: Fri, 20 Mar 2026 11:45:08 +0100 Subject: [PATCH 17/17] avoid division by zero in meltPonds function --- model/finiteelement.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/model/finiteelement.cpp b/model/finiteelement.cpp index 1bf5024ca..d1f233b26 100644 --- a/model/finiteelement.cpp +++ b/model/finiteelement.cpp @@ -6543,6 +6543,7 @@ FiniteElement::meltPonds(const int cpt, const double dt, const double hi, const double hIceMin = 0.1; // minimum ice thickness with ponds (m) const double concMin = 0.1; // minimum ice concentration with ponds const double max_lid_thickness = 0.3; // maximum lid thickness + const double min_lid_thickness = 1e-3; // minimum lid thickness const double ice_to_water = physical::rhoi/physical::rhow; const double snow_to_water = physical::rhos/physical::rhow; @@ -6592,13 +6593,13 @@ FiniteElement::meltPonds(const int cpt, const double dt, const double hi, (M_lid_volume[cpt]+M_pond_volume[cpt])/pond_depth); double delLidVolume = 0; // Volume increase is always positive! - if ( M_lid_volume[cpt] > 0. ) // a lid exits + if ( M_lid_volume[cpt] > 0. && D_pond_fraction[cpt] > 1e-11 ) // a lid exists { // Grow or melt the lid - lid volume is in water-equivalent meters /* Assume the pond water has the same salinity as sea ice and is at the * freezing point */ const double TPond = -M_freezingpoint_mu*physical::si; - const double lidThickness = M_lid_volume[cpt]*water_to_ice/D_pond_fraction[cpt]; + const double lidThickness = std::max(min_lid_thickness, std::min(max_lid_thickness, M_lid_volume[cpt]*water_to_ice/D_pond_fraction[cpt])); const double Qic = (TPond - M_tice[0][cpt]) / lidThickness * physical::ki; const double delLidThickness = ( std::min(Qia-Qic,0.) + Qic ) // surface + bottom *dt/(physical::rhoi*physical::Lf);