Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Improve robustness of conservative spherical polygon interpolation for high resolution meshes #130

Merged
merged 3 commits into from
May 3, 2023
Merged

Improve robustness of conservative spherical polygon interpolation for high resolution meshes #130

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
20094ca
New algorithm for ConvexSphericalPolygon intersection; interpolation …
sbrdar Mar 8, 2023
56454e9
Fix and diagnose problems with ConservativeSphericalInterpolation
sbrdar Mar 8, 2023
5233ae5
Configurable eps in compare_pointxyz via ATLAS_COMPAREPOINTXYZ_EPS_FA…
wdeconinck Apr 26, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1,133 changes: 743 additions & 390 deletions src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc
View file
Open in desktop

Large diffs are not rendered by default.

42 changes: 21 additions & 21 deletions src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h
View file
Open in desktop
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