Skip to content

Commit

Permalink
Merge pull request #130 from ecmwf/feature/highres-conservative-spher…
Browse files Browse the repository at this point in the history 
…ical-polygon-interpolation

Improve robustness of conservative spherical polygon interpolation for high resolution meshes
  • Loading branch information
@wdeconinck
wdeconinck authored May 3, 2023
2 parents a6a049b + 5233ae5 commit 8440795
Show file tree
Hide file tree
Showing 7 changed files with 1,496 additions and 968 deletions.
1,133 changes: 743 additions & 390 deletions src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc
View file

Large diffs are not rendered by default.

42 changes: 21 additions & 21 deletions src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,10 +25,8 @@ class ConservativeSphericalPolygonInterpolation : public Method {
struct InterpolationParameters { // one polygon intersection
std::vector<idx_t> cell_idx; // target cells used for intersection
std::vector<PointXYZ> centroids; // intersection cell centroids
std::vector<double> src_weights; // intersection cell areas

// TODO: tgt_weights can be computed on the fly
std::vector<double> tgt_weights; // (intersection cell areas) / (target cell area)
std::vector<double> weights; // intersection cell areas
std::vector<double> tgt_weights; // (intersection cell areas) / (target cell area) // TODO: tgt_weights vector should be removed for the sake of highres
};

private:
Expand All @@ -55,7 +53,6 @@ class ConservativeSphericalPolygonInterpolation : public Method {
std::vector<std::vector<idx_t>> src_node2csp_;
std::vector<std::vector<idx_t>> tgt_node2csp_;


// Timings
struct Timings {
double source_polygons_assembly{0};
Expand Down Expand Up @@ -102,18 +99,20 @@ class ConservativeSphericalPolygonInterpolation : public Method {
std::array<int, 4> counts;
enum Errors
{
SRC_PLG_L1 = 0, // index, over/undershoot in source subpolygon creation
SRC_PLG_LINF,
TGT_PLG_L1, // index, over/untershoot in target subpolygon creation
TGT_PLG_LINF,
GEO_L1, // index, cumulative area mismatch in polygon intersections
GEO_LINF, // index, like GEO_L1 but in L_infinity norm
GEO_DIFF, // index, difference in earth area coverages
SRC_SUBPLG_L1 = 0, // index, \sum_{cell of mesh} {cell.area - \sum_{subpol of cell} subpol.area}
SRC_SUBPLG_LINF, // index, \max_-||-
TGT_SUBPLG_L1, // see above
TGT_SUBPLG_LINF, // see above
SRC_INTERSECTPLG_L1, // index, \sum_{cell of src-mesh} {cell.area - \sum_{intersect of cell} intersect.area}
SRC_INTERSECTPLG_LINF, // index, \max_-||-
TGT_INTERSECTPLG_L1, // see above
TGT_INTERSECTPLG_LINF, // see above
SRCTGT_INTERSECTPLG_DIFF, // index, 1/(unit_sphere.area) ( \sum_{scell} scell.area - \sum{tcell} tcell.area )
REMAP_CONS, // index, error in mass conservation
REMAP_L2, // index, error accuracy for given analytical function
REMAP_LINF // index, like REMAP_L2 but in L_infinity norm
};
std::array<double, 10> errors;
std::array<double, 11> errors;

double tgt_area_sum;
double src_area_sum;
Expand All @@ -123,8 +122,8 @@ class ConservativeSphericalPolygonInterpolation : public Method {
Statistics();
Statistics(const Metadata&);

void accuracy(const Interpolation& interpolation, const Field target,
std::function<double(const PointLonLat&)> func);
Metadata accuracy(const Interpolation& interpolation, const Field target,
std::function<double(const PointLonLat&)> func);


// compute difference field of source and target mass
Expand All @@ -139,6 +138,7 @@ class ConservativeSphericalPolygonInterpolation : public Method {
void do_setup(const FunctionSpace& src_fs, const FunctionSpace& tgt_fs) override;
void do_setup(const Grid& src_grid, const Grid& tgt_grid, const interpolation::Cache&) override;
void do_execute(const Field& src_field, Field& tgt_field, Metadata&) const override;
void do_execute(const FieldSet& src_fields, FieldSet& tgt_fields, Metadata&) const override;

void print(std::ostream& out) const override;

Expand All @@ -160,17 +160,17 @@ class ConservativeSphericalPolygonInterpolation : public Method {
void intersect_polygons(const CSPolygonArray& src_csp, const CSPolygonArray& tgt_scp);
Matrix compute_1st_order_matrix();
Matrix compute_2nd_order_matrix();
void dump_intersection(const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array,
void dump_intersection(const std::string, const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array,
const std::vector<idx_t>& plg_2_idx_array) const;
template <class TargetCellsIDs>
void dump_intersection(const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array,
void dump_intersection(const std::string, const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array,
const TargetCellsIDs& plg_2_idx_array) const;
std::vector<idx_t> sort_cell_edges(Mesh& mesh, idx_t cell_id) const;
std::vector<idx_t> sort_node_edges(Mesh& mesh, idx_t cell_id) const;
std::vector<idx_t> get_cell_neighbours(Mesh& mesh, idx_t jcell) const;
std::vector<idx_t> get_node_neighbours(Mesh& mesh, idx_t jcell) const;
CSPolygonArray get_polygons_celldata(Mesh& mesh) const;
CSPolygonArray get_polygons_nodedata(Mesh& mesh, std::vector<idx_t>& csp2node,
std::vector<idx_t> get_cell_neighbours(Mesh&, idx_t jcell) const;
std::vector<idx_t> get_node_neighbours(Mesh&, idx_t jcell) const;
CSPolygonArray get_polygons_celldata(FunctionSpace) const;
CSPolygonArray get_polygons_nodedata(FunctionSpace, std::vector<idx_t>& csp2node,
std::vector<std::vector<idx_t>>& node2csp,
std::array<double, 2>& errors) const;

Expand Down
Loading

0 comments on commit 8440795

Please sign in to comment.