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

58 revise layer distribution #59

Merged
merged 14 commits into from
Jan 22, 2025
Merged

58 revise layer distribution #59

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
14d1dc8
Ignore profiles for salinity, simply add by volumes.
erwin-mulder Dec 5, 2024
3268d47
Rewrite of layer distribution, taking into account the sign of the qu…
erwin-mulder Dec 5, 2024
a1f18a5
Removed usage of profiles for concentration quantities (salinity).
erwin-mulder Dec 5, 2024
189dc3b
Improved logging of salinity layers.
erwin-mulder Dec 5, 2024
711a83e
Split volumes in 'to' and 'from' parts, we use proper volumes for the…
erwin-mulder Dec 10, 2024
e41cb84
Added masking for salinity layers based on discharge result.
erwin-mulder Dec 10, 2024
c7c7f36
Fixed typo from->to.
erwin-mulder Dec 10, 2024
826bac0
Fixed quantity names for get
erwin-mulder Dec 10, 2024
9a3e5c3
Renamed variables to adhere to new terminology.
erwin-mulder Dec 17, 2024
c476452
Fixed log message for phase wise update step.
erwin-mulder Dec 18, 2024
b603f91
Made open phasse a a bit more lenient towards inner and outer head di…
erwin-mulder Dec 19, 2024
65aba47
Fixed division by zero while correcting phase values in case lake or …
erwin-mulder Jan 7, 2025
5c4e34f
Addressed PR 58 review comments.
erwin-mulder Jan 21, 2025
79b2b1d
Typo in comment
erwin-mulder Jan 21, 2025
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
94 changes: 56 additions & 38 deletions src/dimr_bmi.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,20 +120,26 @@ int set_var(const char *key, void *src_ptr) {
// Set dest_ptr for intended lock and quantity.
if (match_key(quantity, "salinity_lake")) {
dest_ptr = config.locks[lock_index].parameters3d.salinity_lake;
dest_len = config.locks[lock_index].lake_volumes.num_volumes;
dest_len = config.locks[lock_index].from_lake_volumes.num_volumes;
} else if (match_key(quantity, "head_lake")) {
dest_ptr = &config.locks[lock_index].parameters.head_lake;
} else if (match_key(quantity, "salinity_sea")) {
dest_ptr = config.locks[lock_index].parameters3d.salinity_sea;
dest_len = config.locks[lock_index].sea_volumes.num_volumes;
dest_len = config.locks[lock_index].from_sea_volumes.num_volumes;
} else if (match_key(quantity, "head_sea")) {
dest_ptr = &config.locks[lock_index].parameters.head_sea;
} else if (match_key(quantity, "water_volume_lake")) {
dest_ptr = config.locks[lock_index].lake_volumes.volumes;
dest_len = config.locks[lock_index].lake_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_sea")) {
dest_ptr = config.locks[lock_index].sea_volumes.volumes;
dest_len = config.locks[lock_index].sea_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_from_lake")) {
dest_ptr = config.locks[lock_index].from_lake_volumes.volumes;
dest_len = config.locks[lock_index].from_lake_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_from_sea")) {
dest_ptr = config.locks[lock_index].from_sea_volumes.volumes;
dest_len = config.locks[lock_index].from_sea_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_to_lake")) {
dest_ptr = config.locks[lock_index].to_lake_volumes.volumes;
dest_len = config.locks[lock_index].to_lake_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_to_sea")) {
dest_ptr = config.locks[lock_index].to_sea_volumes.volumes;
dest_len = config.locks[lock_index].to_sea_volumes.num_volumes;
} else if (match_key(quantity, "temperature_lake")) {
dest_ptr = &config.locks[lock_index].parameters.temperature_lake;
} else if (match_key(quantity, "temperature_sea")) {
Expand Down Expand Up @@ -185,40 +191,46 @@ int get_var(const char *key, void **dst_ptr) {
// Set source based on key(s)...
if (match_key(quantity, "mass_transport_lake")) {
source_ptr = config.locks[lock_index].results3d.mass_transport_lake;
source_len = config.locks[lock_index].lake_volumes.num_volumes;
source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
} else if (match_key(quantity, "salt_load_lake")) {
source_ptr = config.locks[lock_index].results3d.salt_load_lake;
source_len = config.locks[lock_index].lake_volumes.num_volumes;
source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
} else if (match_key(quantity, "discharge_from_lake")) {
source_ptr = config.locks[lock_index].results3d.discharge_from_lake;
source_len = config.locks[lock_index].lake_volumes.num_volumes;
source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
} else if (match_key(quantity, "discharge_to_lake")) {
source_ptr = config.locks[lock_index].results3d.discharge_to_lake;
source_len = config.locks[lock_index].lake_volumes.num_volumes;
source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
} else if (match_key(quantity, "salinity_to_lake")) {
source_ptr = config.locks[lock_index].results3d.salinity_to_lake;
source_len = config.locks[lock_index].lake_volumes.num_volumes;
source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
} else if (match_key(quantity, "mass_transport_sea")) {
source_ptr = config.locks[lock_index].results3d.mass_transport_sea;
source_len = config.locks[lock_index].sea_volumes.num_volumes;
source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
} else if (match_key(quantity, "salt_load_sea")) {
source_ptr = config.locks[lock_index].results3d.salt_load_sea;
source_len = config.locks[lock_index].sea_volumes.num_volumes;
source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
} else if (match_key(quantity, "discharge_from_sea")) {
source_ptr = config.locks[lock_index].results3d.discharge_from_sea;
source_len = config.locks[lock_index].sea_volumes.num_volumes;
source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
} else if (match_key(quantity, "discharge_to_sea")) {
source_ptr = config.locks[lock_index].results3d.discharge_to_sea;
source_len = config.locks[lock_index].sea_volumes.num_volumes;
source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
} else if (match_key(quantity, "salinity_to_sea")) {
source_ptr = config.locks[lock_index].results3d.salinity_to_sea;
source_len = config.locks[lock_index].sea_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_lake")) {
source_ptr = config.locks[lock_index].lake_volumes.volumes;
source_len = config.locks[lock_index].lake_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_sea")) {
source_ptr = config.locks[lock_index].sea_volumes.volumes;
source_len = config.locks[lock_index].sea_volumes.num_volumes;
source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_from_lake")) {
source_ptr = config.locks[lock_index].from_lake_volumes.volumes;
source_len = config.locks[lock_index].from_lake_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_from_sea")) {
source_ptr = config.locks[lock_index].from_sea_volumes.volumes;
source_len = config.locks[lock_index].from_sea_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_to_lake")) {
source_ptr = config.locks[lock_index].to_lake_volumes.volumes;
source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_to_sea")) {
source_ptr = config.locks[lock_index].to_sea_volumes.volumes;
source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
} else if (match_key(quantity, "salinity_sea")) {
// NOTE: This is really a GET_VALUE_PTR(), called before the update.
source_ptr = config.locks[lock_index].parameters3d.salinity_sea;
Expand Down Expand Up @@ -267,6 +279,7 @@ int update(double dt) {
if (config.current_time == config.start_time && delta_time > 0) {
// Check if timestep is compatible with all loaded timeseries.
if (!sealock_delta_time_ok(&config.locks[lock_index], delta_time)) {
log_error("delta_time (%d) >= shortest phase in the timeseries!\n", delta_time);
return DIMR_BMI_FAILURE;
}
}
Expand Down Expand Up @@ -308,34 +321,39 @@ int get_var_shape(char *key, int dims[DIMR_BMI_MAXDIMS]) { // dims -> int[6]

// Set source based on key(s)...
if (match_key(quantity, "mass_transport_lake")) {
source_len = config.locks[lock_index].lake_volumes.num_volumes;
source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
} else if (match_key(quantity, "salt_load_lake")) {
source_len = config.locks[lock_index].lake_volumes.num_volumes;
source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
} else if (match_key(quantity, "discharge_from_lake")) {
source_len = config.locks[lock_index].lake_volumes.num_volumes;
source_len = config.locks[lock_index].from_lake_volumes.num_volumes;
} else if (match_key(quantity, "discharge_to_lake")) {
source_len = config.locks[lock_index].lake_volumes.num_volumes;
source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
} else if (match_key(quantity, "salinity_to_lake")) {
source_len = config.locks[lock_index].lake_volumes.num_volumes;
source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
} else if (match_key(quantity, "mass_transport_sea")) {
source_len = config.locks[lock_index].sea_volumes.num_volumes;
source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
} else if (match_key(quantity, "salt_load_sea")) {
source_len = config.locks[lock_index].sea_volumes.num_volumes;
source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
} else if (match_key(quantity, "discharge_from_sea")) {
source_len = config.locks[lock_index].sea_volumes.num_volumes;
source_len = config.locks[lock_index].from_sea_volumes.num_volumes;
} else if (match_key(quantity, "discharge_to_sea")) {
source_len = config.locks[lock_index].sea_volumes.num_volumes;
source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
} else if (match_key(quantity, "salinity_to_sea")) {
source_len = config.locks[lock_index].sea_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_lake")) {
source_len = config.locks[lock_index].lake_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_sea")) {
source_len = config.locks[lock_index].sea_volumes.num_volumes;
source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_from_lake")) {
source_len = config.locks[lock_index].from_lake_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_from_sea")) {
source_len = config.locks[lock_index].from_sea_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_to_lake")) {
source_len = config.locks[lock_index].to_lake_volumes.num_volumes;
} else if (match_key(quantity, "water_volume_to_sea")) {
source_len = config.locks[lock_index].to_sea_volumes.num_volumes;
} else {
log_debug("Unhandled get_var('%s', @%p)\n", key, dims);
source_len = 1;
}

memset(dims, 0, DIMR_BMI_MAXDIMS * sizeof(int));
dims[0] = source_len;
log_info("%s yielded %d for quantity '%s' of lock %d.\n", __func__, source_len,
quantity, lock_index);
Expand Down
73 changes: 47 additions & 26 deletions src/io_layer_distribution.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -139,6 +139,9 @@ int io_normalize_profile(profile_t *profile) {
for (int i = index_after_zero; i <= max_index; i++) {
profile->relative_discharge_from_lock[i] *= correction_factor_after;
}
// Keep zero information.
profile->relative_z_zero = z_zero;
profile->start_sign = start_sign;

// Integrate top.
const double check_integral_before =
Expand Down Expand Up @@ -203,16 +206,16 @@ double integrate_piecewise_linear_profile(const profile_t *profile, const double
return integrated_profile;
}

// Distribute the total_discharge over layers.
// The discharge that each layer receives is given by a relative profile
int distribute_discharge_over_layers(double total_discharge, const profile_t *profile,
// Distribute the total_quantity over layers.
// The quantity that each layer receives is given by a relative profile
int distribute_discharge_over_layers(double total_quantity, const profile_t *profile,
const layers_t *layers,
layered_discharge_t *layered_discharge_result) {
layered_discharge_t *layered_quantity_result) {
assert(layers != NULL);
assert(layers->number_of_layers > 0);
assert(layered_discharge_result != NULL);
assert(layered_discharge_result->number_of_layers == layers->number_of_layers);
assert(layered_discharge_result->discharge_per_layer != NULL);
assert(layered_quantity_result != NULL);
assert(layered_quantity_result->number_of_layers == layers->number_of_layers);
assert(layered_quantity_result->discharge_per_layer != NULL);
assert(profile->relative_z_position != NULL);
assert(profile->relative_discharge_from_lock != NULL);
assert(profile->number_of_positions > 0);
Expand All @@ -222,30 +225,48 @@ int distribute_discharge_over_layers(double total_discharge, const profile_t *pr
return -1;
}

double previous_volume = 0.0;
double next_volume = 0.0;
double profile_integral = 0;
double relative_z_start = 0.0;
double relative_z_end = 1.0;
double relative_z_prev = 0.0;
double relative_z = 0.0;

log_debug("quantity total = %g\n", total_discharge);
if (fsign(total_quantity) == profile->start_sign) {
relative_z_end = profile->relative_z_zero;
} else {
relative_z_start = profile->relative_z_zero;
}

double relative_quantity_layer = 0.0;

log_debug("quantity total = %g\n", total_quantity);
log_debug("num_layers = %d\n", layers->number_of_layers);
for (int layer = 0; layer < layers->number_of_layers; ++layer) {
log_debug("normalized_layer_volume[%d] = %g\n", layer,
layers->normalized_target_volumes[layer]);
next_volume += layers->normalized_target_volumes[layer];
const double relative_discharge_layer =
integrate_piecewise_linear_profile(profile, previous_volume, next_volume);
log_debug("profile_layer_volume [%d] = %g\n", layer, relative_discharge_layer);
double layer_discharge = 0;
if (relative_discharge_layer * total_discharge > 0 || layers->number_of_layers == 1) {
layer_discharge = fabs(relative_discharge_layer) * total_discharge;
profile_integral += relative_discharge_layer;
if (layers->number_of_layers > 1) {
for (int layer = 0; layer < layers->number_of_layers; ++layer) {
vrielin marked this conversation as resolved.
Show resolved Hide resolved
relative_z += layers->normalized_target_volumes[layer];
if (relative_z < relative_z_start || relative_z_prev > relative_z_end) {
relative_quantity_layer = 0.0;
} else if (relative_z >= relative_z_start && relative_z_prev < relative_z_start) {
relative_quantity_layer = integrate_piecewise_linear_profile(profile, relative_z_start, relative_z);
} else if (relative_z > relative_z_end && relative_z_prev >= relative_z_start) {
relative_quantity_layer = integrate_piecewise_linear_profile(profile, relative_z_prev, relative_z_end);
} else {
relative_quantity_layer = integrate_piecewise_linear_profile(profile, relative_z_prev, relative_z);
}
const double layer_quantity = fabs(relative_quantity_layer) * total_quantity;
profile_integral += relative_quantity_layer;
layered_quantity_result->discharge_per_layer[layer] = layer_quantity;
log_debug("layer quantity [%d] = %g (== %g * %g)\n", layer, layer_quantity,
fabs(relative_quantity_layer), total_quantity);
relative_z_prev = relative_z;
}
layered_discharge_result->discharge_per_layer[layer] = layer_discharge;
log_debug("layer quantity [%d] = %g (?= %g * %g)\n", layer, layer_discharge,
fabs(relative_discharge_layer), total_discharge);
previous_volume = next_volume;
} else {
layered_quantity_result->discharge_per_layer[0] = total_quantity;
profile_integral = 1.0;
log_debug("layer quantity [0] = %g (== 1.0 * %g)\n", total_quantity, total_quantity);
}

log_debug("profile integral = %g (should be +/-1)\n\n", profile_integral);

return 0;
}
9 changes: 7 additions & 2 deletions src/io_layer_distribution.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,13 @@ static const struct profile_struct {
relative_z_position; // Relative positions in the z-direction for which profile is provided (0 = bed level, 1 = water level)
double *
relative_discharge_from_lock; // Relative discharge at corresponding relative z-position (non-negative), need not be normalized
} profile_default = {
.number_of_positions = 0, .relative_z_position = NULL, .relative_discharge_from_lock = NULL};
double relative_z_zero;
int start_sign;
} profile_default = {.number_of_positions = 0,
.relative_z_position = NULL,
.relative_discharge_from_lock = NULL,
.relative_z_zero = 1.0,
.start_sign = 0};

typedef struct profile_struct profile_t;

Expand Down
Loading