diff --git a/README.md b/README.md index 8f6e2a6..d8ff6c2 100644 --- a/README.md +++ b/README.md @@ -260,7 +260,7 @@ Options: --yac Use YAC interpolation with default method (avg_arith) --yac-method Use YAC interpolation with specific method; click the method button in the GUI to cycle methods at runtime - --yac-3d Per-depth masked interpolation for 3D variables + --yac-3d Fractional fill-value masking for 3D variables -h, --help Show help message ``` @@ -282,7 +282,7 @@ Options (uterm): --cutoff Cutoff latitude for polar view (default: 60) --yac Use YAC interpolation (default: avg_arith) --yac-method Use YAC interpolation method (requires WITH_YAC=1) - --yac-3d Per-depth masked interpolation for 3D variables + --yac-3d Fractional fill-value masking for 3D variables -h, --help Show help ``` @@ -341,7 +341,7 @@ YAC interpolation (requires `make WITH_YAC=1`): ./ushow temp.nc -m mesh.nc --yac-method nnn4dist # 4-NN distance-weighted ./ushow temp.nc -m mesh.nc --yac-method nnn4gauss # 4-NN Gaussian-weighted ./ushow temp.nc -m mesh.nc --yac-method avg_arith # Cell averaging -./ushow temp.nc -m mesh.nc --yac-3d # Per-depth masking for 3D ocean data +./ushow temp.nc -m mesh.nc --yac-3d # Fractional masking for 3D ocean data ./ushow temp.nc -m mesh.nc --yac-3d --yac-method nnn4dist # Combine with specific method ``` @@ -399,7 +399,7 @@ The test suite includes: - **test_file_netcdf**: NetCDF file I/O - **test_file_mitgcm**: MITgcm MDS binary file I/O - **test_file_zarr**: Zarr file I/O (when built with `WITH_ZARR=1`) -- **test_yac_3d**: Per-depth masked YAC interpolation (when built with `WITH_YAC=1`) +- **test_yac_3d**: Fractional fill-value masking for YAC interpolation (when built with `WITH_YAC=1`) - **test_yac_click**: Time series click-to-node lookup in YAC mode (when built with `WITH_YAC=1`) - **test_integration**: End-to-end workflow tests @@ -460,13 +460,11 @@ By default, ushow uses a fast KDTree nearest-neighbor lookup to map unstructured NNN methods work with any grid type. Averaging methods require element connectivity — for grids that lack it (reduced Gaussian, HEALPix, etc.), ushow auto-generates a triangulation from latitude bands. For MPAS grids, native dual-mesh connectivity (`cellsOnVertex`) is read directly from the file. -### Per-Depth Masking (`--yac-3d`) +### Fractional Fill-Value Masking (`--yac-3d`) For 3D ocean variables (e.g., temperature at 50 depth levels), deeper levels have more land masking — fewer valid source points. Without `--yac-3d`, a single interpolation is used for all depths, so fill values from land points can bleed into the result at deeper levels. -With `--yac-3d`, ushow builds a separate masked interpolation for each depth level lazily (on first visit). Each level only uses valid (non-fill) source points, producing clean coastlines at all depths. Surface levels that have no masking reuse the base interpolation with no overhead. - -The first visit to each new depth level builds a masked interpolation (timing depends on grid size); subsequent visits are instant (cached). The cache is cleared automatically when you switch variables or cycle interpolation methods. +With `--yac-3d`, ushow generates a dynamic fractional mask from the source field before each interpolation step: valid points receive weight 1.0, fill points 0.0. The mask is passed to YAC's fractional interpolation (`yac_interpolation_execute_frac`), which rescales partially-masked target cells and assigns the fill value to fully-masked cells. This produces clean coastlines at every depth level. ## Data Flow diff --git a/src/regrid_yac.c b/src/regrid_yac.c index 1088e3f..1ac34d0 100644 --- a/src/regrid_yac.c +++ b/src/regrid_yac.c @@ -23,6 +23,10 @@ /* YAXT must be initialized before YAC grid operations */ #include +/* Grid name constants used when registering source/target grids with YAC */ +#define YAC_SOURCE_GRID_NAME "source" +#define YAC_TARGET_GRID_NAME "target" + /* Store YAC interpolation object for execute-based apply */ struct USYacRegrid { size_t target_nx, target_ny; @@ -35,6 +39,9 @@ struct USYacRegrid { /* YAC interpolation object (handles redistribution + weight apply) */ struct yac_interpolation *interp; + /* YAC interpolation weights (retained for reuse / inspection) */ + struct yac_interp_weights *weights; + /* Temporary double buffers for float<->double conversion */ double *src_buf; /* [n_src_points] */ double *tgt_buf; /* [n_tgt_points] */ @@ -47,13 +54,58 @@ struct USYacRegrid { double laea_R; USTargetConfig config; /* stored copy for 3D rebuild path */ - /* Per-depth cache (--yac-3d) */ - USMesh *mesh_ref; /* non-owned, for lazy creation */ - struct yac_interpolation **depth_cache; /* [depth_cache_n], NULL = not yet built */ - size_t depth_cache_n; /* 0 = 3d mode not enabled */ + /* Fractional-mask interpolation (--yac-3d) */ + int frac_enabled; /* use frac-mask apply path */ + struct yac_interpolation *interp_frac; /* single interp with frac masking */ + double *frac_buf; /* [n_src_points] frac mask buffer */ }; static int mpi_initialized_by_us = 0; +static int yaxt_initialized_by_us = 0; + +/* ---- Method descriptor table ------------------------------------------- */ + +typedef enum { INTERP_KIND_NNN, INTERP_KIND_AVG, INTERP_KIND_CONSERV } InterpKind; + +typedef struct { + USYacMethod method; + const char *name; + int needs_connectivity; + InterpKind kind; + union { + struct { int type; size_t n; double max_dist; double scale; } nnn; + struct { int type; int partial; } avg; + struct { int order; int enforced; int partial; int norm; } conserv; + } p; +} USYacMethodInfo; + +#define NNN(t,n,md,sc) .kind = INTERP_KIND_NNN, .p.nnn = {(t),(n),(md),(sc)} +#define AVG(t,pa) .kind = INTERP_KIND_AVG, .p.avg = {(t),(pa)} +#define CONSERV(o,e,pa,n) .kind = INTERP_KIND_CONSERV, .p.conserv = {(o),(e),(pa),(n)} + +static const USYacMethodInfo method_table[] = { + { YAC_METHOD_NNN_1, "nnn1", 0, NNN(YAC_INTERP_NNN_DIST, 1, 0.0, 0.0) }, + { YAC_METHOD_NNN_4_DIST, "nnn4dist", 0, NNN(YAC_INTERP_NNN_DIST, 4, 0.0, 0.0) }, + { YAC_METHOD_NNN_4_GAUSS, "nnn4gauss", 0, NNN(YAC_INTERP_NNN_GAUSS, 4, 0.0, + YAC_INTERP_NNN_GAUSS_SCALE_DEFAULT) }, + { YAC_METHOD_AVERAGE_ARITH, "avg_arith", 1, AVG(YAC_INTERP_AVG_ARITHMETIC, 1) }, + { YAC_METHOD_AVERAGE_DIST, "avg_dist", 1, AVG(YAC_INTERP_AVG_DIST, 1) }, + { YAC_METHOD_AVERAGE_BARY, "avg_bary", 1, AVG(YAC_INTERP_AVG_BARY, 1) }, + { YAC_METHOD_CONSERVATIVE_1, "conserv1", 1, CONSERV(1, 0, 1, YAC_INTERP_CONSERV_DESTAREA) }, + { YAC_METHOD_CONSERVATIVE_2, "conserv2", 1, CONSERV(2, 0, 1, YAC_INTERP_CONSERV_DESTAREA) }, +}; + +#undef NNN +#undef AVG +#undef CONSERV + +#define METHOD_COUNT (sizeof(method_table) / sizeof(method_table[0])) + +static const USYacMethodInfo *method_lookup(USYacMethod m) { + for (size_t i = 0; i < METHOD_COUNT; i++) + if (method_table[i].method == m) return &method_table[i]; + return NULL; +} static double get_time_seconds(void) { struct timeval tv; @@ -74,13 +126,15 @@ int yac_regrid_init(void) { /* YAXT must be initialized before any YAC grid operations */ if (!xt_initialized()) { xt_initialize(MPI_COMM_WORLD); + yaxt_initialized_by_us = 1; } return 0; } void yac_regrid_finalize(void) { - if (xt_initialized()) { + if (yaxt_initialized_by_us && xt_initialized()) { xt_finalize(); + yaxt_initialized_by_us = 0; } if (mpi_initialized_by_us) { int finalized = 0; @@ -93,34 +147,15 @@ void yac_regrid_finalize(void) { } const char *yac_method_name(USYacMethod m) { - switch (m) { - case YAC_METHOD_NNN_1: return "nnn1"; - case YAC_METHOD_NNN_4_DIST: return "nnn4dist"; - case YAC_METHOD_NNN_4_GAUSS: return "nnn4gauss"; - case YAC_METHOD_AVERAGE_ARITH: return "avg_arith"; - case YAC_METHOD_AVERAGE_DIST: return "avg_dist"; - case YAC_METHOD_AVERAGE_BARY: return "avg_bary"; - case YAC_METHOD_CONSERVATIVE_1: return "conserv1"; - case YAC_METHOD_CONSERVATIVE_2: return "conserv2"; - default: return "unknown"; - } + const USYacMethodInfo *info = method_lookup(m); + return info ? info->name : "unknown"; } int yac_method_parse(const char *name, USYacMethod *out) { if (!name || !out) return -1; - struct { const char *str; USYacMethod m; } table[] = { - {"nnn1", YAC_METHOD_NNN_1}, - {"nnn4dist", YAC_METHOD_NNN_4_DIST}, - {"nnn4gauss", YAC_METHOD_NNN_4_GAUSS}, - {"avg_arith", YAC_METHOD_AVERAGE_ARITH}, - {"avg_dist", YAC_METHOD_AVERAGE_DIST}, - {"avg_bary", YAC_METHOD_AVERAGE_BARY}, - {"conserv1", YAC_METHOD_CONSERVATIVE_1}, - {"conserv2", YAC_METHOD_CONSERVATIVE_2}, - }; - for (size_t i = 0; i < sizeof(table)/sizeof(table[0]); i++) { - if (strcmp(name, table[i].str) == 0) { - *out = table[i].m; + for (size_t i = 0; i < METHOD_COUNT; i++) { + if (strcmp(name, method_table[i].name) == 0) { + *out = method_table[i].method; return 0; } } @@ -128,20 +163,12 @@ int yac_method_parse(const char *name, USYacMethod *out) { } int yac_method_needs_connectivity(USYacMethod m) { - switch (m) { - case YAC_METHOD_AVERAGE_ARITH: - case YAC_METHOD_AVERAGE_DIST: - case YAC_METHOD_AVERAGE_BARY: - case YAC_METHOD_CONSERVATIVE_1: - case YAC_METHOD_CONSERVATIVE_2: - return 1; - default: - return 0; - } + const USYacMethodInfo *info = method_lookup(m); + return info ? info->needs_connectivity : 0; } /* Build YAC source grid from USMesh */ -static struct yac_basic_grid *build_source_grid(USMesh *mesh, USYacMethod method) { +static struct yac_basic_grid *build_source_grid(USMesh *mesh) { /* For 2D curvilinear or 1D structured grids with known dimensions, * use YAC's native curve_2d support. This avoids the need for * auto-generated triangulation which fails on grids with folds @@ -152,11 +179,11 @@ static struct yac_basic_grid *build_source_grid(USMesh *mesh, USYacMethod method size_t nbr_vertices[2] = {mesh->orig_nx, mesh->orig_ny}; int cyclic[2] = {0, 0}; return yac_basic_grid_curve_2d_deg_new( - "source", nbr_vertices, cyclic, mesh->lon, mesh->lat); + YAC_SOURCE_GRID_NAME, nbr_vertices, cyclic, mesh->lon, mesh->lat); } - if (mesh->n_elements > 0 && mesh->elem_nodes != NULL && - yac_method_needs_connectivity(method)) { + /* If connectivity information is available */ + if (mesh->n_elements > 0 && mesh->elem_nodes != NULL) { /* Unstructured grid with cell connectivity (e.g. FESOM triangles) */ int *num_verts_per_cell = malloc(mesh->n_elements * sizeof(int)); if (!num_verts_per_cell) return NULL; @@ -165,26 +192,7 @@ static struct yac_basic_grid *build_source_grid(USMesh *mesh, USYacMethod method } struct yac_basic_grid *grid = yac_basic_grid_unstruct_deg_new( - "source", - mesh->n_points, - mesh->n_elements, - num_verts_per_cell, - mesh->lon, - mesh->lat, - mesh->elem_nodes); - - free(num_verts_per_cell); - return grid; - } else if (mesh->n_elements > 0 && mesh->elem_nodes != NULL) { - /* Have connectivity but using NNN method - still provide full grid */ - int *num_verts_per_cell = malloc(mesh->n_elements * sizeof(int)); - if (!num_verts_per_cell) return NULL; - for (size_t i = 0; i < mesh->n_elements; i++) { - num_verts_per_cell[i] = mesh->n_vertices; - } - - struct yac_basic_grid *grid = yac_basic_grid_unstruct_deg_new( - "source", + YAC_SOURCE_GRID_NAME, mesh->n_points, mesh->n_elements, num_verts_per_cell, @@ -197,94 +205,81 @@ static struct yac_basic_grid *build_source_grid(USMesh *mesh, USYacMethod method } else { /* Point cloud (e.g. HEALPix, or unstructured without connectivity) */ return yac_basic_grid_cloud_deg_new( - "source", + YAC_SOURCE_GRID_NAME, mesh->n_points, mesh->lon, mesh->lat); } } -/* Convert lon/lat (degrees) to 3D Cartesian on unit sphere */ -static void lonlat_deg_to_xyz(double lon_deg, double lat_deg, double xyz[3]) { - double lon_rad = lon_deg * M_PI / 180.0; - double lat_rad = lat_deg * M_PI / 180.0; - double cos_lat = cos(lat_rad); - xyz[0] = cos_lat * cos(lon_rad); - xyz[1] = cos_lat * sin(lon_rad); - xyz[2] = sin(lat_rad); -} - -/* Build target grid (equirectangular or LAEA polar) */ -static struct yac_basic_grid *build_target_grid( +/* Build equirectangular target grid (global or regional --box) */ +static struct yac_basic_grid *build_target_grid_equirect( double resolution, const USTargetConfig *config, size_t *out_nx, size_t *out_ny) { - /* Default to global equirectangular if no config */ - USTargetConfig default_cfg; - if (!config) { - target_config_init_default(&default_cfg); - config = &default_cfg; - } + double lon_min = config->lon_min, lon_max = config->lon_max; + double lat_min = config->lat_min, lat_max = config->lat_max; - if (config->projection == PROJ_EQUIRECTANGULAR) { - /* Equirectangular (global or regional --box) */ - double lon_min = config->lon_min, lon_max = config->lon_max; - double lat_min = config->lat_min, lat_max = config->lat_max; - - size_t nx = (size_t)((lon_max - lon_min) / resolution); - size_t ny = (size_t)((lat_max - lat_min) / resolution); - if (nx < 1) nx = 1; - if (ny < 1) ny = 1; - - size_t n_lon_verts = nx + 1; - size_t n_lat_verts = ny + 1; - - double *lon_verts = malloc(n_lon_verts * sizeof(double)); - double *lat_verts = malloc(n_lat_verts * sizeof(double)); - if (!lon_verts || !lat_verts) { - free(lon_verts); - free(lat_verts); - return NULL; - } + size_t nx = (size_t)((lon_max - lon_min) / resolution); + size_t ny = (size_t)((lat_max - lat_min) / resolution); + if (nx < 1) nx = 1; + if (ny < 1) ny = 1; - double dlon = (lon_max - lon_min) / (double)nx; - double dlat = (lat_max - lat_min) / (double)ny; + size_t n_lon_verts = nx + 1; + size_t n_lat_verts = ny + 1; - for (size_t i = 0; i <= nx; i++) - lon_verts[i] = lon_min + i * dlon; - for (size_t j = 0; j <= ny; j++) - lat_verts[j] = lat_min + j * dlat; + double *lon_verts = malloc(n_lon_verts * sizeof(double)); + double *lat_verts = malloc(n_lat_verts * sizeof(double)); + if (!lon_verts || !lat_verts) { + free(lon_verts); + free(lat_verts); + return NULL; + } - size_t nbr_vertices[2] = {n_lon_verts, n_lat_verts}; - int cyclic[2] = {0, 0}; + double dlon = (lon_max - lon_min) / (double)nx; + double dlat = (lat_max - lat_min) / (double)ny; - struct yac_basic_grid *grid = yac_basic_grid_reg_2d_deg_new( - "target", nbr_vertices, cyclic, lon_verts, lat_verts); + for (size_t i = 0; i <= nx; i++) + lon_verts[i] = lon_min + i * dlon; + for (size_t j = 0; j <= ny; j++) + lat_verts[j] = lat_min + j * dlat; - free(lon_verts); - free(lat_verts); + size_t nbr_vertices[2] = {n_lon_verts, n_lat_verts}; + int cyclic[2] = {0, 0}; - /* Add cell center coordinates (required for YAC_LOC_CELL interpolation) */ - size_t n_cells = nx * ny; - yac_coordinate_pointer cell_coords = malloc(n_cells * sizeof(*cell_coords)); - if (cell_coords) { - for (size_t j = 0; j < ny; j++) { - double lat_center = lat_min + (j + 0.5) * dlat; - for (size_t i = 0; i < nx; i++) { - double lon_center = lon_min + (i + 0.5) * dlon; - lonlat_deg_to_xyz(lon_center, lat_center, - cell_coords[j * nx + i]); - } + struct yac_basic_grid *grid = yac_basic_grid_reg_2d_deg_new( + YAC_TARGET_GRID_NAME, nbr_vertices, cyclic, lon_verts, lat_verts); + + free(lon_verts); + free(lat_verts); + + /* Add cell center coordinates (required for YAC_LOC_CELL interpolation) */ + size_t n_cells = nx * ny; + yac_coordinate_pointer cell_coords = malloc(n_cells * sizeof(*cell_coords)); + if (cell_coords) { + for (size_t j = 0; j < ny; j++) { + double lat_center = lat_min + (j + 0.5) * dlat; + for (size_t i = 0; i < nx; i++) { + double lon_center = lon_min + (i + 0.5) * dlon; + lonlat_to_cartesian(lon_center, lat_center, + &cell_coords[j * nx + i][0], + &cell_coords[j * nx + i][1], + &cell_coords[j * nx + i][2]); } - yac_basic_grid_add_coordinates_nocpy(grid, YAC_LOC_CELL, cell_coords); } - - *out_nx = nx; - *out_ny = ny; - return grid; + yac_basic_grid_add_coordinates_nocpy(grid, YAC_LOC_CELL, cell_coords); } - /* LAEA polar projection */ + *out_nx = nx; + *out_ny = ny; + return grid; +} + +/* Build LAEA polar target grid (north or south pole) */ +static struct yac_basic_grid *build_target_grid_laea( + double resolution, const USTargetConfig *config, + size_t *out_nx, size_t *out_ny) { + double R = EARTH_RADIUS_M; int pole = (config->projection == PROJ_LAEA_NORTH) ? 1 : -1; double extent = laea_extent_from_cutoff(config->cutoff_lat, pole, R); @@ -329,7 +324,7 @@ static struct yac_basic_grid *build_target_grid( size_t nbr_vertices[2] = {nvx, nvy}; int cyclic[2] = {0, 0}; struct yac_basic_grid *grid = yac_basic_grid_curve_2d_deg_new( - "target", nbr_vertices, cyclic, vert_lon, vert_lat); + YAC_TARGET_GRID_NAME, nbr_vertices, cyclic, vert_lon, vert_lat); free(vert_lon); free(vert_lat); @@ -344,10 +339,16 @@ static struct yac_basic_grid *build_target_grid( double px = -extent + (i + 0.5) * dx; double lo, la; if (laea_inverse(px, py, pole, R, &lo, &la) == 0) { - lonlat_deg_to_xyz(lo, la, cell_coords[j * n + i]); + lonlat_to_cartesian(lo, la, + &cell_coords[j * n + i][0], + &cell_coords[j * n + i][1], + &cell_coords[j * n + i][2]); } else { - lonlat_deg_to_xyz(0.0, (pole > 0) ? 90.0 : -90.0, - cell_coords[j * n + i]); + double pole_lat = (pole > 0) ? 90.0 : -90.0; + lonlat_to_cartesian(0.0, pole_lat, + &cell_coords[j * n + i][0], + &cell_coords[j * n + i][1], + &cell_coords[j * n + i][2]); } } } @@ -363,48 +364,48 @@ static struct yac_basic_grid *build_target_grid( return grid; } +/* Build target grid — dispatches to projection-specific builder */ +static struct yac_basic_grid *build_target_grid( + double resolution, const USTargetConfig *config, + size_t *out_nx, size_t *out_ny) { + + /* Default to global equirectangular if no config */ + USTargetConfig default_cfg; + if (!config) { + target_config_init_default(&default_cfg); + config = &default_cfg; + } + + if (config->projection == PROJ_EQUIRECTANGULAR) + return build_target_grid_equirect(resolution, config, out_nx, out_ny); + + /* PROJ_LAEA_NORTH or PROJ_LAEA_SOUTH */ + return build_target_grid_laea(resolution, config, out_nx, out_ny); +} + /* Configure interpolation stack based on method */ static struct yac_interp_stack_config *build_interp_stack(USYacMethod method) { + const USYacMethodInfo *info = method_lookup(method); + if (!info) return NULL; + struct yac_interp_stack_config *config = yac_interp_stack_config_new(); if (!config) return NULL; - switch (method) { - case YAC_METHOD_NNN_1: + switch (info->kind) { + case INTERP_KIND_NNN: yac_interp_stack_config_add_nnn( - config, YAC_INTERP_NNN_DIST, 1, 0.0, 0.0); - break; - case YAC_METHOD_NNN_4_DIST: - yac_interp_stack_config_add_nnn( - config, YAC_INTERP_NNN_DIST, 4, 0.0, 0.0); - break; - case YAC_METHOD_NNN_4_GAUSS: - yac_interp_stack_config_add_nnn( - config, YAC_INTERP_NNN_GAUSS, 4, 0.0, - YAC_INTERP_NNN_GAUSS_SCALE_DEFAULT); - break; - case YAC_METHOD_AVERAGE_ARITH: - yac_interp_stack_config_add_average( - config, YAC_INTERP_AVG_ARITHMETIC, 1); + config, info->p.nnn.type, info->p.nnn.n, + info->p.nnn.max_dist, info->p.nnn.scale); break; - case YAC_METHOD_AVERAGE_DIST: + case INTERP_KIND_AVG: yac_interp_stack_config_add_average( - config, YAC_INTERP_AVG_DIST, 1); + config, info->p.avg.type, info->p.avg.partial); break; - case YAC_METHOD_AVERAGE_BARY: - yac_interp_stack_config_add_average( - config, YAC_INTERP_AVG_BARY, 1); - break; - case YAC_METHOD_CONSERVATIVE_1: - yac_interp_stack_config_add_conservative( - config, 1, 0, 1, YAC_INTERP_CONSERV_DESTAREA); - break; - case YAC_METHOD_CONSERVATIVE_2: + case INTERP_KIND_CONSERV: yac_interp_stack_config_add_conservative( - config, 2, 0, 1, YAC_INTERP_CONSERV_DESTAREA); + config, info->p.conserv.order, info->p.conserv.enforced, + info->p.conserv.partial, info->p.conserv.norm); break; - default: - yac_interp_stack_config_delete(config); - return NULL; } return config; @@ -424,112 +425,6 @@ static enum yac_location get_tgt_location(USYacMethod method) { return YAC_LOC_CELL; } -/* - * Build a single YAC interpolation object. - * If src_mask is non-NULL, applies it to the source grid corners. - * Returns NULL on failure. - */ -static struct yac_interpolation * -build_interpolation(USMesh *mesh, USYacMethod method, - double resolution, - const int *src_mask, - const USTargetConfig *config) { - /* 1. Build grids */ - struct yac_basic_grid *src_grid = build_source_grid(mesh, method); - if (!src_grid) return NULL; - - /* Apply source mask if provided */ - size_t mask_idx = SIZE_MAX; - if (src_mask) { - mask_idx = yac_basic_grid_add_mask( - src_grid, YAC_LOC_CORNER, src_mask, mesh->n_points, "depth_mask"); - } - - size_t nx, ny; - struct yac_basic_grid *tgt_grid = build_target_grid(resolution, config, &nx, &ny); - if (!tgt_grid) { - yac_basic_grid_delete(src_grid); - return NULL; - } - - /* 2. Create distributed grid pair */ - struct yac_dist_grid_pair *grid_pair = - yac_dist_grid_pair_new(src_grid, tgt_grid, MPI_COMM_SELF); - if (!grid_pair) { - yac_basic_grid_delete(src_grid); - yac_basic_grid_delete(tgt_grid); - return NULL; - } - - /* 3. Set up interpolation fields */ - struct yac_interp_field src_field = { - .location = get_src_location(method), - .coordinates_idx = SIZE_MAX, - .masks_idx = mask_idx - }; - struct yac_interp_field tgt_field = { - .location = get_tgt_location(method), - .coordinates_idx = 0, - .masks_idx = SIZE_MAX - }; - - struct yac_interp_grid *interp_grid = yac_interp_grid_new( - grid_pair, "source", "target", 1, &src_field, tgt_field); - if (!interp_grid) { - yac_dist_grid_pair_delete(grid_pair); - yac_basic_grid_delete(src_grid); - yac_basic_grid_delete(tgt_grid); - return NULL; - } - - /* 4. Configure interpolation method */ - struct yac_interp_stack_config *stack_config = build_interp_stack(method); - if (!stack_config) { - yac_interp_grid_delete(interp_grid); - yac_dist_grid_pair_delete(grid_pair); - yac_basic_grid_delete(src_grid); - yac_basic_grid_delete(tgt_grid); - return NULL; - } - - struct interp_method **methods = yac_interp_stack_config_generate(stack_config); - yac_interp_stack_config_delete(stack_config); - if (!methods) { - yac_interp_grid_delete(interp_grid); - yac_dist_grid_pair_delete(grid_pair); - yac_basic_grid_delete(src_grid); - yac_basic_grid_delete(tgt_grid); - return NULL; - } - - /* 5. Compute weights */ - struct yac_interp_weights *weights = - yac_interp_method_do_search(methods, interp_grid); - yac_interp_method_delete(methods); - if (!weights) { - yac_interp_grid_delete(interp_grid); - yac_dist_grid_pair_delete(grid_pair); - yac_basic_grid_delete(src_grid); - yac_basic_grid_delete(tgt_grid); - return NULL; - } - - /* 6. Create interpolation object */ - struct yac_interpolation *interp = yac_interp_weights_get_interpolation( - weights, - YAC_MAPPING_ON_TGT, 1, YAC_FRAC_MASK_NO_VALUE, - 1.0, 0.0, NULL, 1, 1); - - /* 7. Cleanup intermediates */ - yac_interp_weights_delete(weights); - yac_interp_grid_delete(interp_grid); - yac_dist_grid_pair_delete(grid_pair); - yac_basic_grid_delete(src_grid); - yac_basic_grid_delete(tgt_grid); - - return interp; -} - USYacRegrid *yac_regrid_create(USMesh *mesh, double resolution, USYacMethod method, const USTargetConfig *config) { if (!mesh || mesh->n_points == 0) { @@ -573,7 +468,7 @@ USYacRegrid *yac_regrid_create(USMesh *mesh, double resolution, USYacMethod meth printf("YAC: Creating regrid with method '%s'...\n", yac_method_name(method)); /* 1. Build grids */ - struct yac_basic_grid *src_grid = build_source_grid(mesh, method); + struct yac_basic_grid *src_grid = build_source_grid(mesh); if (!src_grid) { fprintf(stderr, "YAC: Failed to create source grid\n"); return NULL; @@ -630,7 +525,8 @@ USYacRegrid *yac_regrid_create(USMesh *mesh, double resolution, USYacMethod meth }; struct yac_interp_grid *interp_grid = yac_interp_grid_new( - grid_pair, "source", "target", 1, &src_field, tgt_field); + grid_pair, YAC_SOURCE_GRID_NAME, YAC_TARGET_GRID_NAME, 1, + &src_field, tgt_field); if (!interp_grid) { fprintf(stderr, "YAC: Failed to create interp grid\n"); yac_dist_grid_pair_delete(grid_pair); @@ -678,35 +574,16 @@ USYacRegrid *yac_regrid_create(USMesh *mesh, double resolution, USYacMethod meth return NULL; } - /* 6. Create interpolation object (handles data redistribution + weight apply) */ - struct yac_interpolation *interp = yac_interp_weights_get_interpolation( - weights, - YAC_MAPPING_ON_TGT, /* apply weights on target side */ - 1, /* collection_size */ - YAC_FRAC_MASK_NO_VALUE, /* disable fractional masking */ - 1.0, /* scaling_factor */ - 0.0, /* scaling_summand */ - NULL, /* yaxt_exchanger_name */ - 1, /* is_source */ - 1 /* is_target */ - ); - - /* 7. Cleanup intermediates (interpolation object is self-contained) */ - yac_interp_weights_delete(weights); + /* 6. Cleanup intermediates (weights kept in struct for lazy interp creation) */ yac_interp_grid_delete(interp_grid); yac_dist_grid_pair_delete(grid_pair); yac_basic_grid_delete(src_grid); yac_basic_grid_delete(tgt_grid); - if (!interp) { - fprintf(stderr, "YAC: Failed to create interpolation object\n"); - return NULL; - } - - /* 8. Allocate result struct */ + /* 7. Allocate result struct (interp created lazily on first apply) */ USYacRegrid *r = calloc(1, sizeof(USYacRegrid)); if (!r) { - yac_interpolation_delete(interp); + yac_interp_weights_delete(weights); return NULL; } @@ -716,7 +593,8 @@ USYacRegrid *yac_regrid_create(USMesh *mesh, double resolution, USYacMethod meth r->resolution = resolution; r->n_src_points = mesh->n_points; r->n_tgt_points = nx * ny; - r->interp = interp; + r->interp = NULL; /* created lazily on first apply */ + r->weights = weights; r->projection = config->projection; r->laea_R = EARTH_RADIUS_M; r->config = *config; @@ -751,8 +629,20 @@ USYacRegrid *yac_regrid_create(USMesh *mesh, double resolution, USYacMethod meth return r; } -void yac_regrid_apply(const USYacRegrid *r, const float *src, float fill, float *dst) { - if (!r || !src || !dst || !r->interp) return; +void yac_regrid_apply(USYacRegrid *r, const float *src, float fill, float *dst) { + if (!r || !src || !dst || !r->weights) return; + + /* Lazy-create interpolation object on first call */ + if (!r->interp) { + r->interp = yac_interp_weights_get_interpolation( + r->weights, + YAC_MAPPING_ON_TGT, 1, YAC_FRAC_MASK_NO_VALUE, + 1.0, 0.0, NULL, 1, 1); + if (!r->interp) { + fprintf(stderr, "YAC: Failed to create interpolation object\n"); + return; + } + } size_t n_src = r->n_src_points; size_t n_tgt = r->n_tgt_points; @@ -823,96 +713,76 @@ USYacMethod yac_regrid_get_method(const USYacRegrid *r) { return r ? r->method : YAC_METHOD_NNN_1; } -void yac_regrid_enable_3d(USYacRegrid *r, USMesh *mesh) { - if (r) r->mesh_ref = mesh; -} - -void yac_regrid_clear_depth_cache(USYacRegrid *r) { - if (!r || !r->depth_cache) return; - for (size_t i = 0; i < r->depth_cache_n; i++) { - if (r->depth_cache[i] && r->depth_cache[i] != r->interp) { - yac_interpolation_delete(r->depth_cache[i]); - } - } - free(r->depth_cache); - r->depth_cache = NULL; - r->depth_cache_n = 0; +void yac_regrid_enable_frac(USYacRegrid *r) { + if (r) r->frac_enabled = 1; } -int yac_regrid_is_3d(const USYacRegrid *r) { - return r && r->mesh_ref != NULL; +int yac_regrid_frac_enabled(const USYacRegrid *r) { + return r && r->frac_enabled; } -void yac_regrid_apply_3d(USYacRegrid *r, size_t depth_idx, size_t n_depths, - const float *src, float fill, float *dst) { - if (!r || !src || !dst) return; - - /* Lazy-allocate depth cache array */ - if (!r->depth_cache) { - r->depth_cache = calloc(n_depths, sizeof(struct yac_interpolation *)); - if (!r->depth_cache) return; - r->depth_cache_n = n_depths; - } - - /* Build masked interpolation for this depth if not cached */ - if (!r->depth_cache[depth_idx]) { - /* Scan source data for fill values to build mask */ - size_t n = r->n_src_points; - int all_valid = 1; - int *mask = malloc(n * sizeof(int)); - if (!mask) { - /* Fallback to unmasked */ - yac_regrid_apply(r, src, fill, dst); +void yac_regrid_apply_frac(USYacRegrid *r, + const float *src, float fill, float *dst) { + if (!r || !src || !dst || !r->weights) return; + + /* Lazy-create fractional-mask interpolation from cached weights. + * Uses fill value as frac_mask_fallback_value so that fully-masked + * target points receive the fill value automatically. */ + if (!r->interp_frac) { + r->interp_frac = yac_interp_weights_get_interpolation( + r->weights, + YAC_MAPPING_ON_TGT, 1, (double)fill, + 1.0, 0.0, NULL, 1, 1); + if (!r->interp_frac) { + fprintf(stderr, "YAC frac: failed to create frac-mask interpolation\n"); return; } - for (size_t i = 0; i < n; i++) { - mask[i] = (src[i] != fill) ? 1 : 0; - if (!mask[i]) all_valid = 0; - } - - if (all_valid) { - /* No masking needed — reuse base interpolation (sentinel) */ - r->depth_cache[depth_idx] = r->interp; - free(mask); - } else { - double t0 = get_time_seconds(); - struct yac_interpolation *masked = - build_interpolation(r->mesh_ref, r->method, - r->resolution, mask, &r->config); - free(mask); - - if (masked) { - r->depth_cache[depth_idx] = masked; - } else { - /* Fallback to unmasked */ - r->depth_cache[depth_idx] = r->interp; + /* Allocate fractional mask buffer alongside */ + if (!r->frac_buf) { + r->frac_buf = malloc(r->n_src_points * sizeof(double)); + if (!r->frac_buf) { + fprintf(stderr, "YAC frac: failed to allocate frac mask buffer\n"); + return; } - double t1 = get_time_seconds(); - printf("YAC 3D: Built masked interpolation for depth %zu (%.2fs)\n", - depth_idx, t1 - t0); } } - /* Execute interpolation using the cached object */ - struct yac_interpolation *interp = r->depth_cache[depth_idx]; size_t n_src = r->n_src_points; size_t n_tgt = r->n_tgt_points; + /* Convert source data to double and build fractional mask: + * 1.0 = valid source point, 0.0 = masked (fill). + * The source field must be pre-multiplied by the frac mask before + * passing to yac_interpolation_execute_frac. */ double *src_d = r->src_buf; - for (size_t i = 0; i < n_src; i++) - src_d[i] = (double)src[i]; + double *frac = r->frac_buf; + for (size_t i = 0; i < n_src; i++) { + double valid = (src[i] != fill) ? 1.0 : 0.0; + frac[i] = valid; + src_d[i] = (double)src[i] * valid; + } + /* Initialize target to fill */ double *tgt_d = r->tgt_buf; for (size_t i = 0; i < n_tgt; i++) tgt_d[i] = (double)fill; - double *src_field_ptr = src_d; + /* Execute with fractional masking: + * src_fields[collection][field][local], src_frac_masks same layout */ + double *src_field_ptr = src_d; double **src_fields_ptr = &src_field_ptr; double ***src_collection = &src_fields_ptr; + + double *frac_field_ptr = frac; + double **frac_fields_ptr = &frac_field_ptr; + double ***frac_collection = &frac_fields_ptr; + double **tgt_collection = &tgt_d; - yac_interpolation_execute(interp, src_collection, tgt_collection); + yac_interpolation_execute_frac( + r->interp_frac, src_collection, frac_collection, tgt_collection); + /* Convert back to float */ for (size_t i = 0; i < n_tgt; i++) dst[i] = (float)tgt_d[i]; } @@ -928,8 +798,10 @@ int yac_regrid_is_regional(const USYacRegrid *r) { void yac_regrid_free(USYacRegrid *r) { if (!r) return; - yac_regrid_clear_depth_cache(r); + if (r->interp_frac) yac_interpolation_delete(r->interp_frac); if (r->interp) yac_interpolation_delete(r->interp); + if (r->weights) yac_interp_weights_delete(r->weights); + free(r->frac_buf); free(r->src_buf); free(r->tgt_buf); free(r); diff --git a/src/regrid_yac.h b/src/regrid_yac.h index 72ac6cb..5800ed3 100644 --- a/src/regrid_yac.h +++ b/src/regrid_yac.h @@ -61,7 +61,7 @@ USYacRegrid *yac_regrid_create(USMesh *mesh, double resolution, USYacMethod meth * fill: fill value for missing/invalid data * dst: output data [target_ny * target_nx], must be preallocated */ -void yac_regrid_apply(const USYacRegrid *r, const float *src, float fill, float *dst); +void yac_regrid_apply(USYacRegrid *r, const float *src, float fill, float *dst); /* * Get target grid dimensions. @@ -95,25 +95,26 @@ int yac_method_needs_connectivity(USYacMethod m); USYacMethod yac_regrid_get_method(const USYacRegrid *r); /* - * Enable per-depth masked interpolation (--yac-3d). + * Enable fractional-mask mode (--yac-3d). + * In contrast to yac_regrid_apply, the source field is checked for fill + * values, a dynamic fractional mask is generated accordingly, and applied. */ -void yac_regrid_enable_3d(USYacRegrid *r, USMesh *mesh); +void yac_regrid_enable_frac(USYacRegrid *r); /* - * Clear per-depth cache (call on variable change). + * Check if fractional-mask mode is enabled. */ -void yac_regrid_clear_depth_cache(USYacRegrid *r); +int yac_regrid_frac_enabled(const USYacRegrid *r); /* - * Check if 3d mode is enabled. + * Apply regridding with dynamic fractional masking. + * The source field is scanned for fill values; valid points receive a + * fractional mask of 1.0, fill points 0.0. The result is computed via + * yac_interpolation_execute_frac so that partially-masked target cells + * are weighted correctly and fully-masked cells receive the fill value. */ -int yac_regrid_is_3d(const USYacRegrid *r); - -/* - * Apply with per-depth masking. Builds/caches masked interpolation on first access. - */ -void yac_regrid_apply_3d(USYacRegrid *r, size_t depth_idx, size_t n_depths, - const float *src, float fill, float *dst); +void yac_regrid_apply_frac(USYacRegrid *r, + const float *src, float fill, float *dst); /* * Check if regrid uses a non-global region (box or polar projection). diff --git a/src/ushow.c b/src/ushow.c index 9b9943e..f301ec8 100644 --- a/src/ushow.c +++ b/src/ushow.c @@ -574,7 +574,7 @@ static void on_yac_method_cycle(void) { &options.target_config); } if (options.yac_3d) - yac_regrid_enable_3d(yac_regrid_ptr, mesh); + yac_regrid_enable_frac(yac_regrid_ptr); view_set_yac_regrid(view, yac_regrid_ptr); x_update_render_mode_label(yac_method_name( yac_regrid_get_method(yac_regrid_ptr))); @@ -871,7 +871,7 @@ static void print_usage(const char *prog) { fprintf(stderr, " nnn1, nnn4dist, nnn4gauss,\n"); fprintf(stderr, " avg_arith, avg_dist, avg_bary,\n"); fprintf(stderr, " conserv1, conserv2\n"); - fprintf(stderr, " --yac-3d Per-depth masked interpolation for 3D variables\n"); + fprintf(stderr, " --yac-3d Fractional fill-value masking for 3D variables\n"); #endif fprintf(stderr, " --light Use light theme (default: dark)\n"); fprintf(stderr, " -h, --help Show this help\n"); @@ -1252,7 +1252,7 @@ int main(int argc, char *argv[]) { return 1; } if (options.yac_3d) - yac_regrid_enable_3d(yac_regrid_ptr, mesh); + yac_regrid_enable_frac(yac_regrid_ptr); } else #endif { diff --git a/src/ushow.defines.h b/src/ushow.defines.h index 266682b..ed8afe5 100644 --- a/src/ushow.defines.h +++ b/src/ushow.defines.h @@ -280,7 +280,7 @@ typedef struct { USTargetConfig target_config; /* Target grid configuration */ #ifdef HAVE_YAC int yac_method; /* YAC interpolation method (-1 = disabled) */ - int yac_3d; /* Per-depth masked interpolation */ + int yac_3d; /* Fractional fill-value masking */ #endif } USOptions; diff --git a/src/uterm.c b/src/uterm.c index fce811f..702f040 100644 --- a/src/uterm.c +++ b/src/uterm.c @@ -254,7 +254,7 @@ static void print_usage(const char *prog) { fprintf(stderr, " nnn1, nnn4dist, nnn4gauss,\n"); fprintf(stderr, " avg_arith, avg_dist, avg_bary,\n"); fprintf(stderr, " conserv1, conserv2\n"); - fprintf(stderr, " --yac-3d Per-depth masked interpolation for 3D variables\n"); + fprintf(stderr, " --yac-3d Fractional fill-value masking for 3D variables\n"); #endif fprintf(stderr, " -h, --help Show this help\n\n"); @@ -1111,7 +1111,7 @@ int main(int argc, char *argv[]) { return 1; } if (options.yac_3d) - yac_regrid_enable_3d(yac_regrid_ptr, mesh); + yac_regrid_enable_frac(yac_regrid_ptr); } else #endif { diff --git a/src/view.c b/src/view.c index 1398b7a..4203481 100644 --- a/src/view.c +++ b/src/view.c @@ -51,12 +51,6 @@ int view_set_variable(USView *view, USVar *var, USMesh *mesh, USRegrid *regrid) view->mesh = mesh; view->regrid = regrid; -#ifdef HAVE_YAC - if (view->yac_regrid) { - yac_regrid_clear_depth_cache((USYacRegrid*)view->yac_regrid); - } -#endif - /* Get dimension info - use fileset total if available */ if (view->fileset) { #ifdef HAVE_ZARR @@ -510,10 +504,10 @@ int view_update(USView *view) { #ifdef HAVE_YAC if (view->yac_regrid) { USYacRegrid *yr = (USYacRegrid*)view->yac_regrid; - if (yac_regrid_is_3d(yr) && view->n_depths > 1) { - yac_regrid_apply_3d(yr, view->depth_index, view->n_depths, - view->raw_data, view->variable->fill_value, - view->regridded_data); + if (yac_regrid_frac_enabled(yr) && view->n_depths > 1) { + yac_regrid_apply_frac(yr, + view->raw_data, view->variable->fill_value, + view->regridded_data); } else { yac_regrid_apply(yr, view->raw_data, view->variable->fill_value, view->regridded_data); diff --git a/tests/test_yac_3d.c b/tests/test_yac_3d.c index 8eee656..a76b024 100644 --- a/tests/test_yac_3d.c +++ b/tests/test_yac_3d.c @@ -1,9 +1,9 @@ /* - * test_yac_3d.c - Unit tests for per-depth masked YAC interpolation (--yac-3d) + * test_yac_3d.c - Unit tests for fractional fill-value masking (--yac-3d) * * Requires: make WITH_YAC=1 test-yac-3d * - * Tests the depth-cache API: enable_3d, is_3d, clear_depth_cache, apply_3d. + * Tests the frac-mask API: enable_frac, frac_enabled, apply_frac. * Uses a small point cloud so YAC setup is fast (~0.1s per interpolation). */ @@ -69,44 +69,29 @@ static void teardown(void) { TEST(yac3d_not_enabled_by_default) { ASSERT_TRUE(setup()); - ASSERT_FALSE(yac_regrid_is_3d(g_regrid)); + ASSERT_FALSE(yac_regrid_frac_enabled(g_regrid)); teardown(); return 1; } TEST(yac3d_enable) { ASSERT_TRUE(setup()); - yac_regrid_enable_3d(g_regrid, g_mesh); - ASSERT_TRUE(yac_regrid_is_3d(g_regrid)); + yac_regrid_enable_frac(g_regrid); + ASSERT_TRUE(yac_regrid_frac_enabled(g_regrid)); teardown(); return 1; } -TEST(yac3d_is_3d_null_safe) { - ASSERT_FALSE(yac_regrid_is_3d(NULL)); - return 1; -} - -TEST(yac3d_clear_empty_cache) { - ASSERT_TRUE(setup()); - yac_regrid_enable_3d(g_regrid, g_mesh); - /* Clearing before any apply should be safe (no-op) */ - yac_regrid_clear_depth_cache(g_regrid); - ASSERT_TRUE(yac_regrid_is_3d(g_regrid)); - teardown(); - return 1; -} - -TEST(yac3d_clear_null_safe) { - yac_regrid_clear_depth_cache(NULL); +TEST(yac3d_frac_enabled_null_safe) { + ASSERT_FALSE(yac_regrid_frac_enabled(NULL)); return 1; } TEST(yac3d_apply_all_valid_reuses_base) { - /* When all source data is valid, depth cache should use the base - * interpolation (sentinel). Verify by checking output is reasonable. */ + /* When all source data is valid, frac mask is all-ones and the base + * interpolation weights apply unchanged. */ ASSERT_TRUE(setup()); - yac_regrid_enable_3d(g_regrid, g_mesh); + yac_regrid_enable_frac(g_regrid); size_t n_src = g_mesh->n_points; float *src = malloc(n_src * sizeof(float)); @@ -121,8 +106,7 @@ TEST(yac3d_apply_all_valid_reuses_base) { float fill = 1.0e20f; fill_all_valid(src, n_src, fill); - /* Apply for depth 0 of 5 depths */ - yac_regrid_apply_3d(g_regrid, 0, 5, src, fill, dst); + yac_regrid_apply_frac(g_regrid, src, fill, dst); /* Output should have some non-fill values */ int n_valid = 0; @@ -138,9 +122,9 @@ TEST(yac3d_apply_all_valid_reuses_base) { } TEST(yac3d_apply_with_mask_builds_new_interp) { - /* When some source points are fill, a new masked interpolation is built. */ + /* When some source points are fill, the frac mask excludes them. */ ASSERT_TRUE(setup()); - yac_regrid_enable_3d(g_regrid, g_mesh); + yac_regrid_enable_frac(g_regrid); size_t n_src = g_mesh->n_points; float *src = malloc(n_src * sizeof(float)); @@ -156,8 +140,7 @@ TEST(yac3d_apply_with_mask_builds_new_interp) { /* Every 3rd point is fill (masked) */ fill_with_mask(src, n_src, fill, 3); - /* Apply for depth 2 of 5 depths */ - yac_regrid_apply_3d(g_regrid, 2, 5, src, fill, dst); + yac_regrid_apply_frac(g_regrid, src, fill, dst); /* Output should have some non-fill values */ int n_valid = 0; @@ -173,11 +156,10 @@ TEST(yac3d_apply_with_mask_builds_new_interp) { } TEST(yac3d_cached_second_call_fast) { - /* Second call to same depth should use cached interpolation. - * We can't measure time precisely but verify it doesn't crash - * and produces same results. */ + /* Second call with same mask should reuse the frac interpolation + * object and produce identical results. */ ASSERT_TRUE(setup()); - yac_regrid_enable_3d(g_regrid, g_mesh); + yac_regrid_enable_frac(g_regrid); size_t n_src = g_mesh->n_points; float *src = malloc(n_src * sizeof(float)); @@ -194,10 +176,10 @@ TEST(yac3d_cached_second_call_fast) { float fill = 1.0e20f; fill_with_mask(src, n_src, fill, 4); - /* First call builds the cached interpolation */ - yac_regrid_apply_3d(g_regrid, 1, 3, src, fill, dst1); + /* First call builds the frac interpolation */ + yac_regrid_apply_frac(g_regrid, src, fill, dst1); /* Second call should reuse it */ - yac_regrid_apply_3d(g_regrid, 1, 3, src, fill, dst2); + yac_regrid_apply_frac(g_regrid, src, fill, dst2); /* Results should be identical */ for (size_t i = 0; i < n_tgt; i++) { @@ -214,7 +196,7 @@ TEST(yac3d_cached_second_call_fast) { TEST(yac3d_different_depths_different_masks) { /* Two depth levels with different masks should produce different results */ ASSERT_TRUE(setup()); - yac_regrid_enable_3d(g_regrid, g_mesh); + yac_regrid_enable_frac(g_regrid); size_t n_src = g_mesh->n_points; float *src_d0 = malloc(n_src * sizeof(float)); @@ -237,8 +219,8 @@ TEST(yac3d_different_depths_different_masks) { /* Depth 1: heavy masking (every 2nd point) */ fill_with_mask(src_d1, n_src, fill, 2); - yac_regrid_apply_3d(g_regrid, 0, 2, src_d0, fill, dst0); - yac_regrid_apply_3d(g_regrid, 1, 2, src_d1, fill, dst1); + yac_regrid_apply_frac(g_regrid, src_d0, fill, dst0); + yac_regrid_apply_frac(g_regrid, src_d1, fill, dst1); /* The results should differ (different masks, different source data) */ int n_differ = 0; @@ -255,10 +237,10 @@ TEST(yac3d_different_depths_different_masks) { return 1; } -TEST(yac3d_clear_then_rebuild) { - /* After clearing cache, next apply should rebuild */ +TEST(yac3d_reapply_consistent) { + /* Applying twice should produce the same result */ ASSERT_TRUE(setup()); - yac_regrid_enable_3d(g_regrid, g_mesh); + yac_regrid_enable_frac(g_regrid); size_t n_src = g_mesh->n_points; float *src = malloc(n_src * sizeof(float)); @@ -275,12 +257,10 @@ TEST(yac3d_clear_then_rebuild) { float fill = 1.0e20f; fill_with_mask(src, n_src, fill, 5); - /* Build cache */ - yac_regrid_apply_3d(g_regrid, 0, 3, src, fill, dst1); + yac_regrid_apply_frac(g_regrid, src, fill, dst1); - /* Clear and rebuild */ - yac_regrid_clear_depth_cache(g_regrid); - yac_regrid_apply_3d(g_regrid, 0, 3, src, fill, dst2); + /* Reapply */ + yac_regrid_apply_frac(g_regrid, src, fill, dst2); /* Results should be the same */ for (size_t i = 0; i < n_tgt; i++) { @@ -295,9 +275,9 @@ TEST(yac3d_clear_then_rebuild) { } TEST(yac3d_non_3d_apply_still_works) { - /* When 3d is not enabled, regular apply should work fine */ + /* When frac mode is not enabled, regular apply should work fine */ ASSERT_TRUE(setup()); - ASSERT_FALSE(yac_regrid_is_3d(g_regrid)); + ASSERT_FALSE(yac_regrid_frac_enabled(g_regrid)); size_t n_src = g_mesh->n_points; float *src = malloc(n_src * sizeof(float)); @@ -333,7 +313,7 @@ int main(void) { return 1; } - int result = run_all_tests("YAC 3D Per-Depth Masking"); + int result = run_all_tests("YAC Fractional Masking"); yac_regrid_finalize(); return result; diff --git a/tests/test_yac_regional.c b/tests/test_yac_regional.c index 9db6647..e921642 100644 --- a/tests/test_yac_regional.c +++ b/tests/test_yac_regional.c @@ -453,8 +453,8 @@ TEST(box_3d_rebuild_respects_config) { USYacRegrid *r = yac_regrid_create(mesh, 5.0, YAC_METHOD_NNN_1, &cfg); ASSERT_NOT_NULL(r); - yac_regrid_enable_3d(r, mesh); - ASSERT_TRUE(yac_regrid_is_3d(r)); + yac_regrid_enable_frac(r); + ASSERT_TRUE(yac_regrid_frac_enabled(r)); size_t n_src = mesh->n_points; float *src = malloc(n_src * sizeof(float)); @@ -468,7 +468,7 @@ TEST(box_3d_rebuild_respects_config) { float fill = 1.0e20f; - /* Create data with some fill values to force a 3D rebuild */ + /* Create data with some fill values to exercise frac masking */ for (size_t i = 0; i < n_src; i++) { if (i % 3 == 0) src[i] = fill; @@ -476,8 +476,8 @@ TEST(box_3d_rebuild_respects_config) { src[i] = 273.0f + (float)i * 0.01f; } - /* This triggers build_interpolation with stored config */ - yac_regrid_apply_3d(r, 0, 3, src, fill, dst); + /* This triggers frac-mask interpolation with stored config */ + yac_regrid_apply_frac(r, src, fill, dst); /* Grid dims should still match box (not global) */ size_t nx2, ny2;