YAC 3.20.0
Yet Another Coupler
Loading...
Searching...
No Matches
interpolation_utils.h
Go to the documentation of this file.
1// Copyright (c) 2024 The YAC Authors
2//
3// SPDX-License-Identifier: BSD-3-Clause
4
5#ifndef INTERPOLATION_UTILS_H
6#define INTERPOLATION_UTILS_H
7
8#include <math.h>
9
10#include <yaxt.h>
11
12#include "utils_common.h"
13#include "ppm/ppm_xfuncs.h"
14#include "interpolation/interpolation_internal.h"
15#include "interpolation/interp_grid_internal.h"
16#include "instrument_internal.h"
17
118//------------------------------------------------------------------------------
119// Post-processing scale macros applied to the final result.
120//------------------------------------------------------------------------------
121
123#define NO_SCALING(RESULT) (RESULT)
125#define MULT(RESULT) ((RESULT) * scale_factor)
127#define ADD(RESULT) ((RESULT) + scale_summand)
129#define MULT_ADD(RESULT) ((RESULT) * scale_factor + scale_summand)
130
135#define FRAC_MASK_TOL (1e-12)
136
137//------------------------------------------------------------------------------
138// Inner kernels for target field computation
139//------------------------------------------------------------------------------
140
172#define COMPUTE_FIELD_STENCIL_FRAC_WGT(TGT_POS, WEIGHT, SCALE) \
173 { \
174 for (size_t l = 0; l < collection_size; ++l) { \
175 double const * restrict * curr_local_field_data = \
176 src_fields?src_fields[l]:NULL; \
177 double const * restrict * curr_local_frac_mask_data = \
178 src_frac_masks?src_frac_masks[l]:NULL; \
179 double const * restrict curr_remote_field_data = \
180 remote_src_fields[l * num_src_fields]; \
181 double const * restrict curr_remote_frac_mask_data = \
182 remote_src_frac_masks[l * num_src_fields]; \
183 double * restrict curr_tgt_field = tgt_field[l]; \
184 YAC_OMP_PARALLEL \
185 { \
186 YAC_OMP_FOR \
187 for (size_t i = 0; i < tgt_count; ++i) { \
188 double result = 0.0; \
189 double frac_weight_sum = 0.0; \
190 double weight_sum = 0.0; \
191 size_t const k_bound = prefix_num_src_per_tgt[i+1]; \
192 for (size_t k = prefix_num_src_per_tgt[i]; k < k_bound; ++k) { \
193 double const * restrict frac_mask_data; \
194 double const * restrict src_field_data; \
195 if (src_field_idx[k] == SIZE_MAX) { \
196 frac_mask_data = curr_remote_frac_mask_data; \
197 src_field_data = curr_remote_field_data; \
198 } else { \
199 frac_mask_data = \
200 curr_local_frac_mask_data[src_field_idx[k]]; \
201 src_field_data = curr_local_field_data[src_field_idx[k]]; \
202 } \
203 result += src_field_data[src_idx[k]] * (WEIGHT); \
204 frac_weight_sum += frac_mask_data[src_idx[k]] * (WEIGHT); \
205 weight_sum += (WEIGHT); \
206 } \
207 curr_tgt_field[(TGT_POS)] = \
208 (fabs(frac_weight_sum) > FRAC_MASK_TOL)? \
209 (SCALE(result / frac_weight_sum * weight_sum)): \
210 frac_mask_fallback_value; \
211 } \
212 } \
213 } \
214 }
215
240#define COMPUTE_FIELD_STENCIL_FRAC_NOWGT(TGT_POS, SCALE) \
241 { \
242 for (size_t l = 0; l < collection_size; ++l) { \
243 double const * restrict * curr_local_field_data = \
244 src_fields?src_fields[l]:NULL; \
245 double const * restrict * curr_local_frac_mask_data = \
246 src_frac_masks?src_frac_masks[l]:NULL; \
247 double const * restrict curr_remote_field_data = \
248 remote_src_fields[l * num_src_fields]; \
249 double const * restrict curr_remote_frac_mask_data = \
250 remote_src_frac_masks[l * num_src_fields]; \
251 double * restrict curr_tgt_field = tgt_field[l]; \
252 YAC_OMP_PARALLEL \
253 { \
254 YAC_OMP_FOR \
255 for (size_t i = 0; i < tgt_count; ++i) { \
256 double result = 0.0; \
257 double frac_weight_sum = 0.0; \
258 size_t const k_start = prefix_num_src_per_tgt[i]; \
259 size_t const k_bound = prefix_num_src_per_tgt[i+1]; \
260 for (size_t k = k_start; k < k_bound; ++k) { \
261 double const * restrict frac_mask_data; \
262 double const * restrict src_field_data; \
263 if (src_field_idx[k] == SIZE_MAX) { \
264 frac_mask_data = curr_remote_frac_mask_data; \
265 src_field_data = curr_remote_field_data; \
266 } else { \
267 frac_mask_data = \
268 curr_local_frac_mask_data[src_field_idx[k]]; \
269 src_field_data = curr_local_field_data[src_field_idx[k]]; \
270 } \
271 result += src_field_data[src_idx[k]]; \
272 frac_weight_sum += frac_mask_data[src_idx[k]]; \
273 } \
274 double const weight_sum = (double)(k_bound - k_start); \
275 curr_tgt_field[(TGT_POS)] = \
276 (fabs(frac_weight_sum) > FRAC_MASK_TOL)? \
277 (SCALE(result / frac_weight_sum * weight_sum)): \
278 frac_mask_fallback_value; \
279 } \
280 } \
281 } \
282 }
283
301#define COMPUTE_FIELD_STENCIL_NOFRAC_(TGT_POS, WEIGHT, SCALE) \
302 { \
303 for (size_t l = 0; l < collection_size; ++l) { \
304 double const * restrict * curr_local_field_data = \
305 src_fields?src_fields[l]:NULL; \
306 double const * restrict curr_remote_field_data = \
307 remote_src_fields[l * num_src_fields]; \
308 double * restrict curr_tgt_field = tgt_field[l]; \
309 YAC_OMP_PARALLEL \
310 { \
311 YAC_OMP_FOR \
312 for (size_t i = 0; i < tgt_count; ++i) { \
313 double result = 0.0; \
314 size_t const k_bound = prefix_num_src_per_tgt[i+1]; \
315 for (size_t k = prefix_num_src_per_tgt[i]; k < k_bound; ++k) { \
316 double const * restrict src_field_data; \
317 if (src_field_idx[k] == SIZE_MAX) { \
318 src_field_data = curr_remote_field_data; \
319 } else { \
320 src_field_data = curr_local_field_data[src_field_idx[k]]; \
321 } \
322 result += src_field_data[src_idx[k]] * (WEIGHT); \
323 } \
324 curr_tgt_field[(TGT_POS)] = SCALE(result); \
325 } \
326 } \
327 } \
328 }
329
350#define COMPUTE_FIELD_DIRECT_FRAC(SCALE) \
351 { \
352 for (size_t i = 0; i < collection_size; ++i) { \
353 for (size_t j = 0; j < num_src_fields; ++j) { \
354 memcpy(tgt_field[i * num_src_fields + j], \
355 src_fields[i][j], tgt_buffer_sizes[j]); \
356 size_t const tgt_count = \
357 tgt_buffer_sizes[j] / sizeof(***src_fields); \
358 YAC_OMP_PARALLEL \
359 { \
360 YAC_OMP_FOR \
361 for (size_t k = 0; k < tgt_count; ++k) { \
362 if (fabs(src_frac_masks[i][j][k]) > FRAC_MASK_TOL) \
363 tgt_field[i * num_src_fields + j][k] = \
364 SCALE( \
365 tgt_field[i * num_src_fields + j][k] / \
366 src_frac_masks[i][j][k]); \
367 else \
368 tgt_field[i * num_src_fields + j][k] = \
369 frac_mask_fallback_value; \
370 } \
371 } \
372 } \
373 } \
374 }
375
388#define COMPUTE_FIELD_DIRECT_NOFRAC(SCALE) \
389 { \
390 for (size_t i = 0; i < collection_size; ++i) { \
391 for (size_t j = 0; j < num_src_fields; ++j) { \
392 memcpy(tgt_field[i * num_src_fields + j], \
393 src_fields[i][j], tgt_buffer_sizes[j]); \
394 size_t const tgt_count = \
395 tgt_buffer_sizes[j] / sizeof(***src_fields); \
396 YAC_OMP_PARALLEL \
397 { \
398 YAC_OMP_FOR \
399 for (size_t k = 0; k < tgt_count; ++k) \
400 tgt_field[i * num_src_fields + j][k] = \
401 SCALE(tgt_field[i * num_src_fields + j][k]); \
402 } \
403 } \
404 } \
405 }
406
407//------------------------------------------------------------------------------
408// Master dispatch macros for target field computation
409//------------------------------------------------------------------------------
410
417#define COMPUTE_FIELD_STENCIL_FRAC(SCALE) \
418 { \
419 if (weights != NULL) { \
420 if (tgt_pos != NULL) \
421 COMPUTE_FIELD_STENCIL_FRAC_WGT(tgt_pos[i], weights[k], SCALE) \
422 else \
423 COMPUTE_FIELD_STENCIL_FRAC_WGT(i, weights[k], SCALE) \
424 } else { \
425 if (tgt_pos != NULL) \
426 COMPUTE_FIELD_STENCIL_FRAC_NOWGT(tgt_pos[i], SCALE) \
427 else \
428 COMPUTE_FIELD_STENCIL_FRAC_NOWGT(i, SCALE) \
429 } \
430 }
431
438#define COMPUTE_FIELD_STENCIL_NOFRAC(SCALE) \
439 { \
440 if (weights != NULL) { \
441 if (tgt_pos != NULL) \
442 COMPUTE_FIELD_STENCIL_NOFRAC_(tgt_pos[i], weights[k], SCALE) \
443 else \
444 COMPUTE_FIELD_STENCIL_NOFRAC_(i, weights[k], SCALE) \
445 } else { \
446 if (tgt_pos != NULL) \
447 COMPUTE_FIELD_STENCIL_NOFRAC_(tgt_pos[i], 1.0, SCALE) \
448 else \
449 COMPUTE_FIELD_STENCIL_NOFRAC_(i, 1.0, SCALE) \
450 } \
451 }
452
461#define COMPUTE_FIELD(COMPUTE_FIELD_FRAC, COMPUTE_FIELD_NOFRAC) \
462 { \
463 if (YAC_FRAC_MASK_VALUE_IS_VALID(frac_mask_fallback_value)) { \
464 if (scale_factor == 1.0) { \
465 if (scale_summand == 0.0) COMPUTE_FIELD_FRAC(NO_SCALING) \
466 else COMPUTE_FIELD_FRAC(ADD) \
467 } else { \
468 if (scale_summand == 0.0) COMPUTE_FIELD_FRAC(MULT) \
469 else COMPUTE_FIELD_FRAC(MULT_ADD) \
470 } \
471 } else { \
472 if (scale_factor == 1.0) { \
473 if (scale_summand == 0.0) COMPUTE_FIELD_NOFRAC(NO_SCALING) \
474 else COMPUTE_FIELD_NOFRAC(ADD) \
475 } else { \
476 if (scale_summand == 0.0) COMPUTE_FIELD_NOFRAC(MULT) \
477 else COMPUTE_FIELD_NOFRAC(MULT_ADD) \
478 } \
479 } \
480 }
481
482 // end of compute_tgt_field_macros group
483
582static inline void compute_tgt_field_stencil(
583 double const * restrict ** src_fields,
584 double const * restrict ** src_frac_masks,
585 double const * restrict * remote_src_fields,
586 double const * restrict * remote_src_frac_masks,
587 double * restrict * tgt_field,
588 size_t const * restrict tgt_pos,
589 size_t tgt_count, size_t const * restrict prefix_num_src_per_tgt,
590 double const * restrict weights,
591 size_t const * restrict src_field_idx,
592 size_t const * restrict src_idx,
593 size_t num_src_fields, size_t collection_size,
594 double frac_mask_fallback_value,
595 double scale_factor, double scale_summand) {
596
597 YAC_INSTRUMENT_START(YAC_INSTR_COMPUTE_TGT_FIELD);
598
599 // return if there is work to do
600 if (tgt_count != 0) {
601
602 // Generate code for target field computation
603 COMPUTE_FIELD(COMPUTE_FIELD_STENCIL_FRAC, COMPUTE_FIELD_STENCIL_NOFRAC)
604 }
605
606 YAC_INSTRUMENT_STOP(YAC_INSTR_COMPUTE_TGT_FIELD);
607}
608
657static inline void compute_tgt_field_direct(
658 double const * restrict ** src_fields,
659 double const * restrict ** src_frac_masks,
660 double * restrict * tgt_field,
661 size_t * restrict tgt_buffer_sizes,
662 size_t num_src_fields, size_t collection_size,
663 double frac_mask_fallback_value,
664 double scale_factor, double scale_summand) {
665
666 YAC_INSTRUMENT_START(YAC_INSTR_COMPUTE_TGT_FIELD_DIRECT);
667
668 // Generate code for target field computation
669 COMPUTE_FIELD(COMPUTE_FIELD_DIRECT_FRAC, COMPUTE_FIELD_DIRECT_NOFRAC)
670
671 YAC_INSTRUMENT_STOP(YAC_INSTR_COMPUTE_TGT_FIELD_DIRECT);
672}
673
674// Cleanup target field computation macros to avoid pollution of the
675// global namespace.
676#undef COMPUTE_FIELD_STENCIL_FRAC
677#undef COMPUTE_FIELD_STENCIL_NOFRAC
678#undef COMPUTE_FIELD_STENCIL_FRAC_WGT
679#undef COMPUTE_FIELD_STENCIL_FRAC_NOWGT
680#undef COMPUTE_FIELD_STENCIL_NOFRAC_
681#undef COMPUTE_FIELD_DIRECT_NOFRAC
682#undef COMPUTE_FIELD_DIRECT_FRAC
683#undef MULT_ADD
684#undef ADD
685#undef MULT
686#undef NO_SCALING
687#undef COMPUTE_FIELD
688#undef FRAC_MASK_TOL
689
699#define CHECK_WITH_FRAC_MASK() \
700 YAC_ASSERT_F( \
701 with_frac_mask == \
702 YAC_FRAC_MASK_VALUE_IS_VALID(frac_mask_fallback_value), \
703 "with_frac_mask does not match value provided to constructor\n" \
704 "(frac_mask_fallback_value = %lf with_frac_mask = %d " \
705 "with_frac_mask(constructor) %d)", frac_mask_fallback_value, \
706 frac_mask_fallback_value != YAC_FRAC_MASK_NO_VALUE, with_frac_mask)
707
708// -----------------------------------------------------------------------
709// Base assertion macros for interp_grid operations
710// -----------------------------------------------------------------------
711
712#define YAC_INTERP_GRID_ASSERT_F(INTERP_GRID, EXP, FMT, ...) \
713 YAC_ASSERT_F( \
714 (EXP), FMT "(source grid: %s; target grid: %s)", \
715 __VA_ARGS__, \
716 yac_interp_grid_get_src_grid_name(INTERP_GRID), \
717 yac_interp_grid_get_tgt_grid_name(INTERP_GRID))
718
719#define YAC_INTERP_GRID_ASSERT(INTERP_GRID, EXP, MSG) \
720 YAC_INTERP_GRID_ASSERT_F(INTERP_GRID, EXP, MSG "%s", "")
721
722#define YAC_INTERP_GRID_FIELD_ASSERT_F(INTERP_GRID, EXP, FMT, ...) \
723 YAC_ASSERT_F( \
724 (EXP), FMT \
725 "(source grid: %s; source field location: %s; " \
726 "target grid: %s; target field location: %s)", \
727 __VA_ARGS__, \
728 yac_interp_grid_get_src_grid_name(INTERP_GRID), \
729 yac_loc2str(yac_interp_grid_get_src_field_location(INTERP_GRID, 0)), \
730 yac_interp_grid_get_tgt_grid_name(INTERP_GRID), \
731 yac_loc2str(yac_interp_grid_get_tgt_field_location(INTERP_GRID)))
732
733#define YAC_INTERP_GRID_FIELD_ASSERT(INTERP_GRID, EXP, MSG) \
734 YAC_INTERP_GRID_FIELD_ASSERT_F(INTERP_GRID, EXP, MSG "%s", "")
735
736// -----------------------------------------------------------------------
737// Validation macros for interp_grid operations
738// -----------------------------------------------------------------------
739
740// Helper macro to check if source fields exist
741#define CHECK_INTERP_GRID_HAS_SRC_FIELDS(INTERP_GRID) \
742 YAC_INTERP_GRID_ASSERT( \
743 INTERP_GRID, yac_interp_grid_get_num_src_fields(INTERP_GRID) != 0, \
744 "Number of source fields is 0.")
745
746#define CHECK_INTERP_GRID_SRC_FIELD_IDX_VALID(INTERP_GRID, SRC_FIELD_IDX) \
747 CHECK_INTERP_GRID_HAS_SRC_FIELDS(INTERP_GRID) \
748 YAC_INTERP_GRID_ASSERT_F( \
749 INTERP_GRID, (SRC_FIELD_IDX) < yac_interp_grid_get_num_src_fields(INTERP_GRID), \
750 "invalid src_field_idx. Got %zu, expected < %zu", \
751 (size_t)(SRC_FIELD_IDX), yac_interp_grid_get_num_src_fields(INTERP_GRID))
752
753#define CHECK_INTERP_GRID_SRC_FIELD_LOCATION_CORNER_CELL_EDGE(INTERP_GRID, LOCATION) \
754 CHECK_INTERP_GRID_HAS_SRC_FIELDS(INTERP_GRID) \
755 YAC_INTERP_GRID_FIELD_ASSERT( \
756 INTERP_GRID, \
757 ((LOCATION) == YAC_LOC_CORNER) || \
758 ((LOCATION) == YAC_LOC_CELL) || \
759 ((LOCATION) == YAC_LOC_EDGE), \
760 "Expected source fields defined as CORNER, CELL, or EDGE")
761
762#define CHECK_INTERP_GRID_NUM_SRC_FIELDS_SINGLE(INTERP_GRID) \
763 CHECK_INTERP_GRID_HAS_SRC_FIELDS(INTERP_GRID) \
764 YAC_INTERP_GRID_FIELD_ASSERT_F( \
765 INTERP_GRID, \
766 yac_interp_grid_get_num_src_fields(INTERP_GRID) == 1, \
767 "Invalid number of source fields. Got %zu, expected 1.", \
768 yac_interp_grid_get_num_src_fields(INTERP_GRID))
769
770#define CHECK_INTERP_GRID_SRC_FIELD_IS_CELL(INTERP_GRID, SRC_FIELD_IDX) \
771 CHECK_INTERP_GRID_SRC_FIELD_IDX_VALID(INTERP_GRID, SRC_FIELD_IDX) \
772 YAC_INTERP_GRID_FIELD_ASSERT( \
773 INTERP_GRID, \
774 yac_interp_grid_get_src_field_location(INTERP_GRID, SRC_FIELD_IDX) == YAC_LOC_CELL, \
775 "Expected source fields defined as CELL")
776
777#define CHECK_INTERP_GRID_TGT_FIELD_IS_CELL(INTERP_GRID) \
778 YAC_INTERP_GRID_FIELD_ASSERT( \
779 INTERP_GRID, \
780 yac_interp_grid_get_tgt_field_location(INTERP_GRID) == YAC_LOC_CELL, \
781 "Expected target fields defined as CELL")
782
783// -----------------------------------------------------------------------
784// Interpolation buffer helpers
785// -----------------------------------------------------------------------
786
791
792struct yac_interpolation_buffer {
793
794 double ** buffer;
795 size_t * buffer_sizes;
796};
797
798struct yac_interpolation_buffer yac_interpolation_buffer_init(
799 Xt_redist * redists, size_t num_fields, size_t collection_size,
800 enum yac_interpolation_buffer_type type);
801
802struct yac_interpolation_buffer yac_interpolation_buffer_init_2(
803 Xt_redist * redists, size_t * min_buffer_sizes, size_t num_fields,
804 size_t collection_size, enum yac_interpolation_buffer_type type);
805
806struct yac_interpolation_buffer yac_interpolation_buffer_copy(
807 struct yac_interpolation_buffer buffer, size_t num_fields,
808 size_t collection_size);
809
810void yac_interpolation_buffer_free(struct yac_interpolation_buffer * buffer);
811
812#endif // INTERPOLATION_UTILS_H
num_fields
unsigned num_fields
Definition OASIS3_MCT_example.c:34
COMPUTE_FIELD
#define COMPUTE_FIELD(COMPUTE_FIELD_FRAC, COMPUTE_FIELD_NOFRAC)
Definition interpolation_utils.h:461
COMPUTE_FIELD_STENCIL_NOFRAC
#define COMPUTE_FIELD_STENCIL_NOFRAC(SCALE)
Definition interpolation_utils.h:438
COMPUTE_FIELD_DIRECT_FRAC
#define COMPUTE_FIELD_DIRECT_FRAC(SCALE)
Definition interpolation_utils.h:350
COMPUTE_FIELD_DIRECT_NOFRAC
#define COMPUTE_FIELD_DIRECT_NOFRAC(SCALE)
Definition interpolation_utils.h:388
COMPUTE_FIELD_STENCIL_FRAC
#define COMPUTE_FIELD_STENCIL_FRAC(SCALE)
Definition interpolation_utils.h:417
instrument_internal.h
YAC_INSTRUMENT_START
#define YAC_INSTRUMENT_START(ID)
Enter an instrumented region.
Definition instrument_internal.h:60
YAC_INSTRUMENT_STOP
#define YAC_INSTRUMENT_STOP(ID)
Exit an instrumented region entered with YAC_INSTRUMENT_START(ID).
Definition instrument_internal.h:66
interp_grid_internal.h
type
enum callback_type type
Definition interp_method_check.c:137
Xt_redist
struct Xt_redist_ * Xt_redist
Definition interpolation.h:29
interpolation_internal.h
Defines internal basic interpolation definitions.
yac_interpolation_buffer_copy
struct yac_interpolation_buffer yac_interpolation_buffer_copy(struct yac_interpolation_buffer buffer, size_t num_fields, size_t collection_size)
Definition interpolation_utils.c:116
compute_tgt_field_stencil
static void compute_tgt_field_stencil(double const *restrict **src_fields, double const *restrict **src_frac_masks, double const *restrict *remote_src_fields, double const *restrict *remote_src_frac_masks, double *restrict *tgt_field, size_t const *restrict tgt_pos, size_t tgt_count, size_t const *restrict prefix_num_src_per_tgt, double const *restrict weights, size_t const *restrict src_field_idx, size_t const *restrict src_idx, size_t num_src_fields, size_t collection_size, double frac_mask_fallback_value, double scale_factor, double scale_summand)
Compute target field values using stencil-based summation of source data (sparse-matrix-vector multip...
Definition interpolation_utils.h:582
yac_interpolation_buffer_init_2
struct yac_interpolation_buffer yac_interpolation_buffer_init_2(Xt_redist *redists, size_t *min_buffer_sizes, size_t num_fields, size_t collection_size, enum yac_interpolation_buffer_type type)
Definition interpolation_utils.c:100
yac_interpolation_buffer_free
void yac_interpolation_buffer_free(struct yac_interpolation_buffer *buffer)
Definition interpolation_utils.c:127
compute_tgt_field_direct
static void compute_tgt_field_direct(double const *restrict **src_fields, double const *restrict **src_frac_masks, double *restrict *tgt_field, size_t *restrict tgt_buffer_sizes, size_t num_src_fields, size_t collection_size, double frac_mask_fallback_value, double scale_factor, double scale_summand)
Copy source field values to target buffers, optionally applying fractional masking and scaling.
Definition interpolation_utils.h:657
yac_interpolation_buffer_type
yac_interpolation_buffer_type
Definition interpolation_utils.h:787
SEND_BUFFER
@ SEND_BUFFER
Definition interpolation_utils.h:788
RECV_BUFFER
@ RECV_BUFFER
Definition interpolation_utils.h:789
yac_interpolation_buffer_init
struct yac_interpolation_buffer yac_interpolation_buffer_init(Xt_redist *redists, size_t num_fields, size_t collection_size, enum yac_interpolation_buffer_type type)
Definition interpolation_utils.c:88
ppm_xfuncs.h
add versions of standard API functions not returning on error
yac_interpolation_buffer
Definition interpolation_utils.h:792
yac_interpolation_buffer::buffer
double ** buffer
Definition interpolation_utils.h:794
yac_interpolation_buffer::buffer_sizes
size_t * buffer_sizes
Definition interpolation_utils.h:795
collection_size
int collection_size
Definition test_dummy_coupling_c.c:24
buffer
double * buffer
Definition toy_coupling.c:39
utils_common.h