Skip to content
Snippets Groups Projects
Commit dfcc66f6 authored by Moritz Hanke's avatar Moritz Hanke
Browse files

fixes issue in HCSBB interpolation method

* HCSBB could cause an issue for regular cells that has two of its vertices at a pole. The triangulation would cae on of the two triangle to be a single line, which was not correctly handled.
parent 1dc77cb6
No related branches found
No related tags found
No related merge requests found
......@@ -198,7 +198,8 @@ struct bnd_triangle_reorder {
// compute the spherical barycentric coordinates of the given vertices with
// respect to the given triangle
static void compute_sb_coords(
double * sb_coords, size_t num_vertices, double triangle[3][3]) {
double * sb_coords, size_t num_vertices, double triangle[3][3],
char const * caller) {
double A[3][3];
lapack_int n = 3, nrhs = (lapack_int) num_vertices, lda = n, ldx = n, ipiv[3];
......@@ -210,10 +211,19 @@ static void compute_sb_coords(
// where: A is the matrix consisting of the vertex coordinates of the three
// corners of the triangle
// we compute b by solving this linear system using LAPACK
YAC_ASSERT(
YAC_ASSERT_F(
!LAPACKE_dgesv(
LAPACK_COL_MAJOR, n, nrhs, &A[0][0], lda, ipiv, sb_coords, ldx),
"ERROR: internal error (could not solve linear 3x3 system)")
"ERROR(%s::compute_sb_coords): "
"internal error (could not solve linear 3x3 system)\n"
"(vector: (% .6e;% .6e;% .6e)\n"
" triangle: ((% .6e;% .6e;% .6e),\n"
" (% .6e;% .6e;% .6e),\n"
" (% .6e;% .6e;% .6e))",
caller, sb_coords[0], sb_coords[1], sb_coords[2],
triangle[0][0], triangle[0][1], triangle[0][2],
triangle[1][0], triangle[1][1], triangle[1][2],
triangle[2][0], triangle[2][1], triangle[2][2])
}
static inline int compare_size_t(const void * a, const void * b) {
......@@ -488,8 +498,10 @@ static void compute_d_sbb_polynomials_3d(
// compute the spherical barycentric coordinates of the point and the
// direction
compute_sb_coords(&(sb_coord_p[0]), 1, bnd_triangle);
compute_sb_coords(&(sb_coord_d[0]), 1, bnd_triangle);
compute_sb_coords(
&(sb_coord_p[0]), 1, bnd_triangle, "compute_d_sbb_polynomials_3d");
compute_sb_coords(
&(sb_coord_d[0]), 1, bnd_triangle, "compute_d_sbb_polynomials_3d");
#define POW_0(x) (1.0)
#define POW_1(x) ((x))
......@@ -836,7 +848,8 @@ static void compute_triangle_coefficients(
double b_g[3][3] = {{g[0][0], g[0][1], g[0][2]},
{g[1][0], g[1][1], g[1][2]},
{g[2][0], g[2][1], g[2][2]}};
compute_sb_coords(&(b_g[0][0]), 3, corner_coordinate_xyz);
compute_sb_coords(
&(b_g[0][0]), 3, corner_coordinate_xyz, "compute_triangle_coefficients");
// computing quadratic spherical Bernstein-Bezier polynomials for the edge
// middle points (e.g. B_020(w0)) using w0 as an example
......@@ -1819,15 +1832,23 @@ static int check_polygon(
memcpy(
triangle_coords[j], point_coords[triangle[j]], 3 * sizeof(double));
sb_coords[0] = check_coord[0];
sb_coords[1] = check_coord[1];
sb_coords[2] = check_coord[2];
compute_sb_coords(sb_coords, 1, triangle_coords);
// if the point is in the current triangle
ret = (sb_coords[0] > -SB_COORD_TOL) &&
(sb_coords[1] > -SB_COORD_TOL) &&
(sb_coords[2] > -SB_COORD_TOL);
// check whether this is actually a triangle
if ((
ret =
!points_are_identically(triangle_coords[0], triangle_coords[1]) &&
!points_are_identically(triangle_coords[0], triangle_coords[2]) &&
!points_are_identically(triangle_coords[1], triangle_coords[2]))) {
sb_coords[0] = check_coord[0];
sb_coords[1] = check_coord[1];
sb_coords[2] = check_coord[2];
compute_sb_coords(sb_coords, 1, triangle_coords, "check_polygon");
// if the point is in the current triangle
ret = (sb_coords[0] > -SB_COORD_TOL) &&
(sb_coords[1] > -SB_COORD_TOL) &&
(sb_coords[2] > -SB_COORD_TOL);
}
}
return ret;
......
......@@ -3,6 +3,7 @@
// SPDX-License-Identifier: BSD-3-Clause
#include <stdlib.h>
#include <string.h>
#include "tests.h"
#include "test_common.h"
......@@ -764,6 +765,179 @@ int main(void) {
yac_basic_grid_delete(grids[0]);
}
{ // test in which the source grid goes up to the pole
// parallel interpolation process
// corner and cell ids for a 7 x 7 grid
// 56-----57-----58-----59-----60-----61-----62-----63
// | | | | | | | |
// | 42 | 43 | 44 | 45 | 46 | 47 | 48 |
// | | | | | | | |
// 48-----49-----50-----51-----52-----53-----54-----55
// | | | | | | | |
// | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
// | | | | | | | |
// 40-----41-----42-----43-----44-----45-----46-----47
// | | | | | | | |
// | 28 | 29 | 30 | 31 | 32 | 33 | 34 |
// | | | | | | | |
// 32-----33-----34-----35-----36-----37-----38-----39
// | | | | | | | |
// | 21 | 22 | 23 | 24 | 25 | 26 | 27 |
// | | | | | | | |
// 24-----25-----26-----27-----28-----29-----30-----31
// | | | | | | | |
// | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
// | | | | | | | |
// 16-----17-----18-----19-----20-----21-----22-----23
// | | | | | | | |
// | 07 | 08 | 09 | 10 | 11 | 12 | 13 |
// | | | | | | | |
// 08-----09-----10-----11-----12-----13-----14-----15
// | | | | | | | |
// | 00 | 01 | 02 | 03 | 04 | 05 | 06 |
// | | | | | | | |
// 00-----01-----02-----03-----04-----05-----06-----07
//
// the grid is distributed among the processes as follows:
// (index == process)
//
// 4---4---4---4---4---4---4---4
// | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
// 4---4---4---4---4---4---4---4
// | 4 | 4 | 4 | 4 | 4 | 4 | 4 |
// 3---3---3---3---4---4---4---4
// | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
// 2---2---2---2---3---3---3---3
// | 2 | 2 | 2 | 2 | 3 | 3 | 3 |
// 1---1---1---1---2---2---2---2
// | 1 | 1 | 1 | 1 | 3 | 3 | 3 |
// 0---0---0---0---1---1---1---1
// | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
// 0---0---0---0---0---0---0---0
// | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
// 0---0---0---0---0---0---0---0
//
//---------------
// setup
//---------------
int is_tgt = split_comm_size == 1;
double coordinates_x[2][10] =
{{72.00,72.25,72.50,72.75,73.00,73.25,73.5,73.75},{-1.0}};
double coordinates_y[2][10] =
{{90.0,89.75,89.50,89.25,89.00,88.75,88.50,88.25},{-1.0}};
size_t const num_cells[2][2] = {{7,7}, {9,9}};
size_t local_start[2][5][2] = {{{0,0},{0,2},{0,3},{4,2},{0,5}}, {{0,0}}};
size_t local_count[2][5][2] = {{{7,2},{4,1},{4,2},{3,3},{7,2}}, {{9,9}}};
int with_halo = 0;
for (size_t i = 0; i < 10; ++i) {
coordinates_x[1][i] = 72.50 + 0.90 * ((double)i / 9.0);
coordinates_y[1][i] = 89.06 + 0.90 * ((double)i / 9.0);
}
for (size_t i = 0; i <= num_cells[is_tgt][0]; ++i)
coordinates_x[is_tgt][i] *= YAC_RAD;
for (size_t i = 0; i <= num_cells[is_tgt][1]; ++i)
coordinates_y[is_tgt][i] *= YAC_RAD;
struct yac_basic_grid_data grid_data =
yac_generate_basic_grid_data_reg2d(
coordinates_x[is_tgt], coordinates_y[is_tgt], num_cells[is_tgt],
local_start[is_tgt][split_comm_rank],
local_count[is_tgt][split_comm_rank], with_halo);
struct yac_basic_grid * grids[2] =
{yac_basic_grid_new(grid_names[is_tgt], grid_data),
yac_basic_grid_empty_new(grid_names[is_tgt^1])};
if (!is_tgt) {
yac_coordinate_pointer vertex_coordinates =
xmalloc(grid_data.num_vertices * sizeof(*vertex_coordinates));
memcpy(
vertex_coordinates, grid_data.vertex_coordinates,
grid_data.num_vertices * sizeof(*vertex_coordinates));
// set x and y coordinate of all pole points to 0.0 in order to
// reproduce an error that occured in ICON
for (size_t i = 0; i < grid_data.num_vertices; ++i) {
if (fabs(1.0 - fabs(vertex_coordinates[i][2])) < 1e-9) {
vertex_coordinates[i][0] = 0.0;
vertex_coordinates[i][1] = 0.0;
}
}
yac_basic_grid_add_coordinates_nocpy(
grids[0], YAC_LOC_CORNER, vertex_coordinates);
}
struct yac_dist_grid_pair * grid_pair =
yac_dist_grid_pair_new(grids[0], grids[1], MPI_COMM_WORLD);
struct yac_interp_field src_fields[] =
{{.location = YAC_LOC_CORNER, .coordinates_idx = 0, .masks_idx = SIZE_MAX}};
size_t num_src_fields = sizeof(src_fields) / sizeof(src_fields[0]);
struct yac_interp_field tgt_field =
{.location = YAC_LOC_CORNER, .coordinates_idx = SIZE_MAX, .masks_idx = SIZE_MAX};
struct yac_interp_grid * interp_grid =
yac_interp_grid_new(grid_pair, grid_names[0], grid_names[1],
num_src_fields, src_fields, tgt_field);
struct interp_method * method_stack[] =
{yac_interp_method_hcsbb_new(),
yac_interp_method_fixed_new(-1.0), NULL};
struct yac_interp_weights * weights =
yac_interp_method_do_search(method_stack, interp_grid);
for (size_t i = 0; i < NUM_REORDER_TYPES; ++i) {
struct yac_interpolation * interpolation =
yac_interp_weights_get_interpolation(
weights, reorder_types[i], 1,
YAC_FRAC_MASK_NO_VALUE, 1.0, 0.0);
// check generated interpolation
{
double * src_field = NULL;
double ** src_fields = &src_field;
double * tgt_field = NULL;
double * ref_tgt_field = NULL;
if (is_tgt) {
tgt_field = xmalloc(grid_data.num_vertices * sizeof(*tgt_field));
ref_tgt_field =
xmalloc(grid_data.num_vertices * sizeof(*ref_tgt_field));
compute_field_data_XYZ(
grid_data.vertex_coordinates, 100, ref_tgt_field);
for (size_t i = 0; i < 100; ++i) tgt_field[i] = -1;
} else {
src_field = xmalloc(grid_data.num_vertices * sizeof(*src_field));
compute_field_data_XYZ(
grid_data.vertex_coordinates, grid_data.num_vertices, src_field);
}
yac_interpolation_execute(interpolation, &src_fields, &tgt_field);
if (is_tgt)
for (size_t i = 0; i < 100; ++i)
if (fabs(1.0 - ref_tgt_field[i] / tgt_field[i]) > 1e-2)
PUT_ERR("wrong interpolation result");
free(ref_tgt_field);
free(tgt_field);
free(src_field);
}
yac_interpolation_delete(interpolation);
}
yac_interp_weights_delete(weights);
yac_interp_method_delete(method_stack);
yac_interp_grid_delete(interp_grid);
yac_dist_grid_pair_delete(grid_pair);
yac_basic_grid_delete(grids[1]);
yac_basic_grid_delete(grids[0]);
}
MPI_Comm_free(&split_comm);
xt_finalize();
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment