From e18b7d02df6b1ee2343498af2a8fc1acf46a1c79 Mon Sep 17 00:00:00 2001 From: "clara.bayley" Date: Mon, 9 Sep 2024 23:03:47 +0200 Subject: [PATCH 1/4] fix: fix calculation of source index for winds using yac edges fields --- libs/coupldyn_yac/yac_cartesian_dynamics.cpp | 108 ++++++++++++++----- libs/coupldyn_yac/yac_cartesian_dynamics.hpp | 17 ++- 2 files changed, 91 insertions(+), 34 deletions(-) diff --git a/libs/coupldyn_yac/yac_cartesian_dynamics.cpp b/libs/coupldyn_yac/yac_cartesian_dynamics.cpp index 9a13e463..465e21e9 100644 --- a/libs/coupldyn_yac/yac_cartesian_dynamics.cpp +++ b/libs/coupldyn_yac/yac_cartesian_dynamics.cpp @@ -9,7 +9,7 @@ * Author: Wilton Loch (WL) * Additional Contributors: Clara Bayley (CB) * ----- - * Last Modified: Saturday 24th August 2024 + * Last Modified: Monday 9th September 2024 * Modified By: CB * ----- * License: BSD 3-Clause "New" or "Revised" License @@ -132,48 +132,98 @@ void create_grid_and_points_definitions(const Config &config, const std::array &target_array, const size_t ndims_north, const size_t ndims_east, const size_t vertical_levels, - double conversion_factor = 1.0) const { + const double conversion_factor = 1.0) const { int info, error; yac_cget(yac_field_id, vertical_levels, yac_raw_data, &info, &error); - for (size_t j = 0; j < ndims_north; j++) { - for (size_t i = 0; i < ndims_east; i++) { - for (size_t k = 0; k < vertical_levels; k++) { - auto vertical_idx = k; - auto source_idx = j * ndims_east + i; - auto ii = (ndims_east * j + i) * vertical_levels + k; - target_array[ii] = yac_raw_data[vertical_idx][source_idx] / conversion_factor; + switch (grid_points) { + case 0: + /* handles fields defined on centres of cells of horizontal + (2-D) grid (grid_points = CENTRES), e.g. temp or wvel + */ + for (size_t j = 0; j < ndims_north; j++) { + for (size_t i = 0; i < ndims_east; i++) { + for (size_t k = 0; k < vertical_levels; k++) { + auto vertical_idx = k; + auto source_idx = ndims_east * j + i; + auto ii = (ndims_east * j + i) * vertical_levels + k; + target_array[ii] = yac_raw_data[vertical_idx][source_idx] / conversion_factor; + } + } } - } + return; + case 1: + /* handles fields defined on longitude edges of cells of + horizontal grid (grid_points = LONGITUDE_EDGES), e.g uvel + */ + for (size_t j = 0; j < ndims_north; j++) { + for (size_t i = 0; i < ndims_east; i++) { + for (size_t k = 0; k < vertical_levels; k++) { + auto vertical_idx = k; + auto edges_ndims_east = ndims_east * 2 - 1; + auto source_idx = edges_ndims_east * j + ndims_east - 1 + i; + auto ii = (ndims_east * j + i) * vertical_levels + k; + target_array[ii] = yac_raw_data[vertical_idx][source_idx] / conversion_factor; + } + } + } + return; + + case 2: + /* handles fields defined on latitude edges of cells of + horizontal grid (grid_points = LATITUDE_EDGES), e.g. vvel + */ + for (size_t j = 0; j < ndims_north; j++) { + for (size_t i = 0; i < ndims_east; i++) { + for (size_t k = 0; k < vertical_levels; k++) { + auto vertical_idx = k; + auto edges_ndims_east = ndims_east * 2 + 1; + auto source_idx = edges_ndims_east * j + i; + auto ii = (ndims_east * j + i) * vertical_levels + k; + target_array[ii] = yac_raw_data[vertical_idx][source_idx] / conversion_factor; + } + } + } + return; } } /* This subroutine is the main entry point for receiving data from YAC. * It checks the dimensionality of the simulation based on the config data. */ void CartesianDynamics::receive_fields_from_yac() { - receive_yac_field(temp_yac_id, yac_raw_cell_data, temp, ndims[NORTHWARD], ndims[EASTWARD], - ndims[VERTICAL], dlc::TEMP0); - receive_yac_field(pressure_yac_id, yac_raw_cell_data, press, ndims[NORTHWARD], ndims[EASTWARD], - ndims[VERTICAL], dlc::P0); - receive_yac_field(qvap_yac_id, yac_raw_cell_data, qvap, ndims[NORTHWARD], ndims[EASTWARD], - ndims[VERTICAL]); - receive_yac_field(qcond_yac_id, yac_raw_cell_data, qcond, ndims[NORTHWARD], ndims[EASTWARD], - ndims[VERTICAL]); - - receive_yac_field(vertical_wind_yac_id, yac_raw_vertical_wind_data, wvel, ndims[NORTHWARD], - ndims[EASTWARD], ndims[VERTICAL] + 1, dlc::W0); - - receive_yac_field(eastward_wind_yac_id, yac_raw_edge_data, uvel, ndims[NORTHWARD], - ndims[EASTWARD] + 1, ndims[VERTICAL], dlc::W0); - - receive_yac_field(northward_wind_yac_id, yac_raw_edge_data, vvel, ndims[NORTHWARD] + 1, - ndims[EASTWARD], ndims[VERTICAL], dlc::W0); + enum grid_points { + CENTRES, + LONGITUDE_EDGES, + LATITUDE_EDGES, + }; + + receive_yac_field(CENTRES, temp_yac_id, yac_raw_cell_data, temp, ndims[NORTHWARD], + ndims[EASTWARD], ndims[VERTICAL], dlc::TEMP0); + receive_yac_field(CENTRES, pressure_yac_id, yac_raw_cell_data, press, ndims[NORTHWARD], + ndims[EASTWARD], ndims[VERTICAL], dlc::P0); + receive_yac_field(CENTRES, qvap_yac_id, yac_raw_cell_data, qvap, ndims[NORTHWARD], + ndims[EASTWARD], ndims[VERTICAL]); + receive_yac_field(CENTRES, qcond_yac_id, yac_raw_cell_data, qcond, ndims[NORTHWARD], + ndims[EASTWARD], ndims[VERTICAL]); + + receive_yac_field(CENTRES, vertical_wind_yac_id, yac_raw_vertical_wind_data, wvel, + ndims[NORTHWARD], ndims[EASTWARD], ndims[VERTICAL] + 1, dlc::W0); + + receive_yac_field(LONGITUDE_EDGES, eastward_wind_yac_id, yac_raw_edge_data, uvel, + ndims[NORTHWARD], ndims[EASTWARD] + 1, ndims[VERTICAL], dlc::W0); + + receive_yac_field(LATITUDE_EDGES, northward_wind_yac_id, yac_raw_edge_data, vvel, + ndims[NORTHWARD] + 1, ndims[EASTWARD], ndims[VERTICAL], dlc::W0); } CartesianDynamics::CartesianDynamics(const Config &config, const std::array i_ndims, diff --git a/libs/coupldyn_yac/yac_cartesian_dynamics.hpp b/libs/coupldyn_yac/yac_cartesian_dynamics.hpp index 63ee13c7..2b457bd9 100644 --- a/libs/coupldyn_yac/yac_cartesian_dynamics.hpp +++ b/libs/coupldyn_yac/yac_cartesian_dynamics.hpp @@ -9,7 +9,7 @@ * Author: Clara Bayley (CB) * Additional Contributors: * ----- - * Last Modified: Saturday 24th August 2024 + * Last Modified: Wednesday 4th September 2024 * Modified By: CB * ----- * License: BSD 3-Clause "New" or "Revised" License @@ -123,10 +123,17 @@ struct CartesianDynamics { /* Public call to receive data from YAC * If the problem is 2D turns into a wrapper for receive_hor_slice_from_yac */ void receive_fields_from_yac(); - void receive_yac_field(unsigned int yac_field_id, double **yac_raw_data, - std::vector &target_array, const size_t ndims_north, - const size_t ndims_east, const size_t vertical_levels, - double conversion_factor) const; + + /* + fill's target_array with values from yac_raw_data at multiplied by their + conversion factor. Special handling when grid_points = 1 or 2 such that only + edge grid points for every other layer are used (to correctly choose only edges + in that direction and not another when using YAC's edge point definition). + */ + void receive_yac_field(const unsigned int grid_points, const unsigned int yac_field_id, + double **yac_raw_data, std::vector &target_array, + const size_t ndims_north, const size_t ndims_east, + const size_t vertical_levels, const double conversion_factor) const; }; /* type satisfying CoupledDyanmics solver concept From a76614a06b6f0e86b80302a44cbec9fb4b62fda6 Mon Sep 17 00:00:00 2001 From: "clara.bayley" Date: Wed, 11 Sep 2024 21:00:13 +0200 Subject: [PATCH 2/4] fix: avoid parastici edge case when coord initially lies exactly on gridbox boundary --- libs/cartesiandomain/cartesianmotion.hpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/libs/cartesiandomain/cartesianmotion.hpp b/libs/cartesiandomain/cartesianmotion.hpp index 12aba5a7..90da8b4e 100644 --- a/libs/cartesiandomain/cartesianmotion.hpp +++ b/libs/cartesiandomain/cartesianmotion.hpp @@ -9,7 +9,7 @@ * Author: Clara Bayley (CB) * Additional Contributors: * ----- - * Last Modified: Tuesday 16th April 2024 + * Last Modified: Wednesday 11th September 2024 * Modified By: CB * ----- * License: BSD 3-Clause "New" or "Revised" License @@ -54,8 +54,11 @@ struct CartesianCheckBounds { KOKKOS_INLINE_FUNCTION void operator()(const unsigned int idx, const Kokkos::pair bounds, const double coord) const { + constexpr double atol = 1e-14; // tolerance to bounds check + const double lower_lim = bounds.first - atol; + const double upper_lim = bounds.second + atol; const bool bad_gbxindex((idx != outofbounds_gbxindex()) && - ((coord < bounds.first) | (coord >= bounds.second))); + ((coord < lower_lim) || (coord >= upper_lim))); assert((!bad_gbxindex) && "SD not in previous gbx nor a neighbour." From 49b8f0119d3b88b7c31d447dc64c92adff0c38be Mon Sep 17 00:00:00 2001 From: "clara.bayley" Date: Wed, 11 Sep 2024 21:43:06 +0200 Subject: [PATCH 3/4] refactor: try example for entire bubble domain --- examples/bubble3d/bubble3d_inputfiles.py | 6 +++--- examples/bubble3d/src/config/bubble3d_config.yaml | 14 +++++++------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/examples/bubble3d/bubble3d_inputfiles.py b/examples/bubble3d/bubble3d_inputfiles.py index f417200d..cbaaa6d1 100644 --- a/examples/bubble3d/bubble3d_inputfiles.py +++ b/examples/bubble3d/bubble3d_inputfiles.py @@ -73,13 +73,13 @@ def main( zgrid = get_zgrid(icon_grid_file, num_vertical_levels) # [m] xgrid = [ 0, - 30000, + 100000, 2500, ] # evenly spaced xhalf coords [m] # distance must match longitude in config file ygrid = [ 0, - 12000, - 6000, + 20000, + 5000, ] # evenly spaced xhalf coords [m] # distance must match latitudes in config file ### --- settings for initial superdroplets --- ### diff --git a/examples/bubble3d/src/config/bubble3d_config.yaml b/examples/bubble3d/src/config/bubble3d_config.yaml index b42da1c5..dd52da5c 100644 --- a/examples/bubble3d/src/config/bubble3d_config.yaml +++ b/examples/bubble3d/src/config/bubble3d_config.yaml @@ -25,13 +25,13 @@ ### SDM Runtime Parameters ### domain: nspacedims : 3 # no. of spatial dimensions to model - ngbxs : 576 # total number of Gbxs - maxnsupers: 576 # maximum number of SDs + ngbxs : 3840 # total number of Gbxs + maxnsupers: 3840 # maximum number of SDs timesteps: CONDTSTEP : 2 # time between SD condensation [s] COLLTSTEP : 2 # time between SD collision [s] - MOTIONTSTEP : 3 # time between SDM motion [s] + MOTIONTSTEP : 5 # time between SDM motion [s] COUPLTSTEP : 60 # time between dynamic couplings [s] OBSTSTEP : 60 # time between SDM observations [s] T_END : 7200 # time span of integration from 0s to T_END [s] @@ -54,7 +54,7 @@ outputdata: coupled_dynamics: type: yac - lower_longitude: -0.9424777965 - upper_longitude: 0.9424777965 - lower_latitude: -0.392699082 - upper_latitude: 0.392699082 + lower_longitude: -3.29867229 + upper_longitude: 2.98451302 + lower_latitude: -1.2575 + upper_latitude: 1.2575 From db3f114fe7907537506acaa9290e512211fb8f58 Mon Sep 17 00:00:00 2001 From: "clara.bayley" Date: Fri, 20 Sep 2024 23:26:51 +0200 Subject: [PATCH 4/4] refactor: alternative attempt to fix calculation of source index for winds using yac edges fields --- libs/coupldyn_yac/yac_cartesian_dynamics.cpp | 14 +++++++++----- libs/coupldyn_yac/yac_cartesian_dynamics.hpp | 3 ++- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/libs/coupldyn_yac/yac_cartesian_dynamics.cpp b/libs/coupldyn_yac/yac_cartesian_dynamics.cpp index 465e21e9..6857fe6f 100644 --- a/libs/coupldyn_yac/yac_cartesian_dynamics.cpp +++ b/libs/coupldyn_yac/yac_cartesian_dynamics.cpp @@ -9,7 +9,7 @@ * Author: Wilton Loch (WL) * Additional Contributors: Clara Bayley (CB) * ----- - * Last Modified: Monday 9th September 2024 + * Last Modified: Friday 20th September 2024 * Modified By: CB * ----- * License: BSD 3-Clause "New" or "Revised" License @@ -170,8 +170,8 @@ void CartesianDynamics::receive_yac_field(const unsigned int grid_points, for (size_t i = 0; i < ndims_east; i++) { for (size_t k = 0; k < vertical_levels; k++) { auto vertical_idx = k; - auto edges_ndims_east = ndims_east * 2 - 1; - auto source_idx = edges_ndims_east * j + ndims_east - 1 + i; + auto source_idx = (ndims_east - 1) * j + ndims_east * j; + source_idx += std::min(2 * i + 1, 2 * ndims_east - 2); auto ii = (ndims_east * j + i) * vertical_levels + k; target_array[ii] = yac_raw_data[vertical_idx][source_idx] / conversion_factor; } @@ -187,8 +187,12 @@ void CartesianDynamics::receive_yac_field(const unsigned int grid_points, for (size_t i = 0; i < ndims_east; i++) { for (size_t k = 0; k < vertical_levels; k++) { auto vertical_idx = k; - auto edges_ndims_east = ndims_east * 2 + 1; - auto source_idx = edges_ndims_east * j + i; + auto source_idx = ndims_east * j + (ndims_east + 1) * j; + if (j < ndims_north - 1) { + source_idx += 2 * i; + } else { + source_idx += i; + } auto ii = (ndims_east * j + i) * vertical_levels + k; target_array[ii] = yac_raw_data[vertical_idx][source_idx] / conversion_factor; } diff --git a/libs/coupldyn_yac/yac_cartesian_dynamics.hpp b/libs/coupldyn_yac/yac_cartesian_dynamics.hpp index 2b457bd9..f61bb306 100644 --- a/libs/coupldyn_yac/yac_cartesian_dynamics.hpp +++ b/libs/coupldyn_yac/yac_cartesian_dynamics.hpp @@ -9,7 +9,7 @@ * Author: Clara Bayley (CB) * Additional Contributors: * ----- - * Last Modified: Wednesday 4th September 2024 + * Last Modified: Friday 20th September 2024 * Modified By: CB * ----- * License: BSD 3-Clause "New" or "Revised" License @@ -23,6 +23,7 @@ #ifndef LIBS_COUPLDYN_YAC_YAC_CARTESIAN_DYNAMICS_HPP_ #define LIBS_COUPLDYN_YAC_YAC_CARTESIAN_DYNAMICS_HPP_ +#include #include #include #include