From 04ec178d1ffcd3beeb8eb28067f8aaae17e5ccd5 Mon Sep 17 00:00:00 2001 From: danieldouglas92 Date: Sat, 24 Feb 2024 14:43:48 -0500 Subject: [PATCH] Address comments, add Assets --- include/world_builder/utilities.h | 4 +-- .../temperature/half_space_model.cc | 8 ++++-- .../temperature/plate_model.cc | 9 ++++-- .../temperature/mass_conserving.cc | 7 ++--- source/world_builder/parameters.cc | 22 +++++++-------- source/world_builder/utilities.cc | 28 +++++++++++++------ 6 files changed, 46 insertions(+), 32 deletions(-) diff --git a/include/world_builder/utilities.h b/include/world_builder/utilities.h index 8ac9ce2a5..dd80e4cce 100644 --- a/include/world_builder/utilities.h +++ b/include/world_builder/utilities.h @@ -465,8 +465,8 @@ namespace WorldBuilder std::vector> mid_oceanic_spreading_velocities, const std::unique_ptr &coordinate_system, const Objects::NaturalCoordinate &position_in_natural_coordinates_at_min_depth, - std::vector> subducting_plate_velocities = {{0.0}}, - std::vector ridge_migration_times = {0.0}); + const std::vector> &subducting_plate_velocities, + const std::vector &ridge_migration_times); // todo_effective /** diff --git a/source/world_builder/features/oceanic_plate_models/temperature/half_space_model.cc b/source/world_builder/features/oceanic_plate_models/temperature/half_space_model.cc index 06d26a266..b84b82bbf 100644 --- a/source/world_builder/features/oceanic_plate_models/temperature/half_space_model.cc +++ b/source/world_builder/features/oceanic_plate_models/temperature/half_space_model.cc @@ -84,7 +84,7 @@ namespace WorldBuilder "in degree Kelvin for this feature. If the model has an adiabatic gradient" "this should be the mantle potential temperature, and T = Tad + Thalf. "); - prm.declare_entry("spreading velocity", Types::OneOf(Types::Double(0.05),Types::Array(Types::ValueAtPoints(0.05, std::numeric_limits::max()))), + prm.declare_entry("spreading velocity", Types::OneOf(Types::Double(0.05),Types::Array(Types::ValueAtPoints(0.05, std::numeric_limits::max()))), "The spreading velocity of the plate in meter per year. " "This is the velocity with which one side moves away from the ridge."); @@ -154,6 +154,8 @@ namespace WorldBuilder Objects::NaturalCoordinate position_in_natural_coordinates_at_min_depth = Objects::NaturalCoordinate(position, *(world->parameters.coordinate_system)); position_in_natural_coordinates_at_min_depth.get_ref_depth_coordinate() += depth-min_depth; + std::vector> subducting_plate_velocities = {{0}}; + std::vector ridge_migration_times = {0.0}; double bottom_temperature_local = bottom_temperature; @@ -167,7 +169,9 @@ namespace WorldBuilder std::vector ridge_parameters = Utilities::calculate_ridge_distance_and_spreading(mid_oceanic_ridges, spreading_velocities_at_each_ridge_point, world->parameters.coordinate_system, - position_in_natural_coordinates_at_min_depth); + position_in_natural_coordinates_at_min_depth, + subducting_plate_velocities, + ridge_migration_times); diff --git a/source/world_builder/features/oceanic_plate_models/temperature/plate_model.cc b/source/world_builder/features/oceanic_plate_models/temperature/plate_model.cc index 9e63ecaf3..a0bbf28b4 100644 --- a/source/world_builder/features/oceanic_plate_models/temperature/plate_model.cc +++ b/source/world_builder/features/oceanic_plate_models/temperature/plate_model.cc @@ -78,7 +78,7 @@ namespace WorldBuilder prm.declare_entry("bottom temperature", Types::Double(-1), "The temperature in degree Kelvin which this feature should have"); - prm.declare_entry("spreading velocity", Types::OneOf(Types::Double(0.05),Types::Array(Types::ValueAtPoints(0.05, std::numeric_limits::max()))), + prm.declare_entry("spreading velocity", Types::OneOf(Types::Double(0.05),Types::Array(Types::ValueAtPoints(0.05, std::numeric_limits::max()))), "The spreading velocity of the plate in meter per year. " "This is the velocity with which one side moves away from the ridge."); @@ -150,7 +150,8 @@ namespace WorldBuilder Objects::NaturalCoordinate position_in_natural_coordinates_at_min_depth = Objects::NaturalCoordinate(position, *(world->parameters.coordinate_system)); position_in_natural_coordinates_at_min_depth.get_ref_depth_coordinate() += depth-min_depth; - + std::vector> subducting_plate_velocities = {{0}}; + std::vector ridge_migration_times = {0.0}; double bottom_temperature_local = bottom_temperature; if (bottom_temperature_local < 0) @@ -165,7 +166,9 @@ namespace WorldBuilder std::vector ridge_parameters = Utilities::calculate_ridge_distance_and_spreading(mid_oceanic_ridges, spreading_velocities_at_each_ridge_point, world->parameters.coordinate_system, - position_in_natural_coordinates_at_min_depth); + position_in_natural_coordinates_at_min_depth, + subducting_plate_velocities, + ridge_migration_times); const double thermal_diffusivity = this->world->thermal_diffusivity; const double age = ridge_parameters[1] / ridge_parameters[0]; diff --git a/source/world_builder/features/subducting_plate_models/temperature/mass_conserving.cc b/source/world_builder/features/subducting_plate_models/temperature/mass_conserving.cc index 9ecd231b1..e7f14142e 100644 --- a/source/world_builder/features/subducting_plate_models/temperature/mass_conserving.cc +++ b/source/world_builder/features/subducting_plate_models/temperature/mass_conserving.cc @@ -117,10 +117,10 @@ namespace WorldBuilder prm.declare_entry("density", Types::Double(3300), "The reference density of the subducting plate in $kg/m^3$"); - prm.declare_entry("plate velocity", Types::OneOf(Types::Double(0.05),Types::Array(Types::ValueAtPoints(0.05, std::numeric_limits::max()))), + prm.declare_entry("plate velocity", Types::OneOf(Types::Double(0.05),Types::Array(Types::ValueAtPoints(0.05, std::numeric_limits::max()))), "The velocity with which the plate subducts in meters per year. Default is 5 cm/yr"); - prm.declare_entry("subducting velocity", Types::OneOf(Types::Double(0.0), Types::Array(Types::Array(Types::Double(0.0), 1), 1)), + prm.declare_entry("subducting velocity", Types::OneOf(Types::Double(-1), Types::Array(Types::Array(Types::Double(-1), 1), 1)), "The velocity with which the ridge is moving through time, and how long the ridge " "has been moving. First value is the velocity, second is the time. Default is [0 cm/yr, 0 yr]"); @@ -188,7 +188,6 @@ namespace WorldBuilder ridge_spreading_velocities = prm.get_value_at_array("plate velocity"); subducting_velocities = prm.get_vector_or_double("subducting velocity"); - mantle_coupling_depth = prm.get("coupling depth"); forearc_cooling_factor = prm.get("forearc cooling factor"); @@ -286,7 +285,7 @@ namespace WorldBuilder const double depth_to_reference_surface = distance_from_planes.depth_reference_surface; const double total_segment_length = additional_parameters.total_local_segment_length; const double average_angle = distance_from_planes.average_angle; - + std::vector slab_ages = calculate_effective_trench_and_plate_ages(ridge_parameters, distance_along_plane); const double seconds_in_year = 60.0 * 60.0 * 24.0 * 365.25; // sec/y diff --git a/source/world_builder/parameters.cc b/source/world_builder/parameters.cc index eeea8f6f8..ae2d974dc 100644 --- a/source/world_builder/parameters.cc +++ b/source/world_builder/parameters.cc @@ -764,7 +764,7 @@ namespace WorldBuilder // now get the array of points. Value *array_of_doubles = Pointer((base + "/1").c_str()).Get(parameters); - double testing_value; + double value_in_array; if (array_of_doubles != nullptr) { for (size_t coordinate_i = 0; coordinate_i < array_of_doubles->Size(); ++coordinate_i ) @@ -776,7 +776,7 @@ namespace WorldBuilder // That means that there are exactly two values per item try { - testing_value = Pointer((base + "/1/" + std::to_string(coordinate_i) + "/" + std::to_string(coordinate_j)).c_str()).Get(parameters)->GetDouble(); + value_in_array = Pointer((base + "/1/" + std::to_string(coordinate_i) + "/" + std::to_string(coordinate_j)).c_str()).Get(parameters)->GetDouble(); } catch (...) { @@ -787,11 +787,10 @@ namespace WorldBuilder << "\"."); } - result.first.emplace_back(value); - result.second.emplace_back(testing_value); + result.second.emplace_back(value_in_array); } } - + result.first.emplace_back(value); } } } @@ -903,7 +902,6 @@ namespace WorldBuilder // start with adding the additional points with the default value // to do this we need the default value double default_value = 0; - std::vector please_work; bool is_array = true; if (Pointer((strict_base + "/" + name).c_str()).Get(parameters) != nullptr && Pointer((strict_base + "/" + name).c_str()).Get(parameters)->IsArray()) { @@ -949,8 +947,8 @@ namespace WorldBuilder // now get the array of points. // Value *coordinates_array = Pointer((base + "/0").c_str()).Get(parameters); Value *coordinates_array = Pointer((base).c_str()).Get(parameters); - double testing_value; - std::vector testing_array; + double value_in_array; + std::vector array_of_values; if (coordinates_array != nullptr) { for (size_t coordinate_i = 0; coordinate_i < coordinates_array->Size(); ++coordinate_i ) @@ -959,7 +957,7 @@ namespace WorldBuilder // That means that there are exactly two values per item try { - testing_value = Pointer(((base + "/") + std::to_string(coordinate_i)).c_str()).Get(parameters)->GetDouble(); + value_in_array = Pointer(((base + "/") + std::to_string(coordinate_i)).c_str()).Get(parameters)->GetDouble(); } catch (...) { @@ -969,16 +967,16 @@ namespace WorldBuilder << Pointer((base + "/" + std::to_string(i) + "/" + std::to_string(coordinate_i)).c_str()).Get(parameters)->GetString() << "\"."); } - testing_array.emplace_back(testing_value); + array_of_values.emplace_back(value_in_array); } - result.emplace_back(testing_array); + result.emplace_back(array_of_values); } } } else { - // case 1: there one value, not an array + // case 1: there's one value, not an array double value = 0; try { diff --git a/source/world_builder/utilities.cc b/source/world_builder/utilities.cc index 22053705a..e581e4b33 100644 --- a/source/world_builder/utilities.cc +++ b/source/world_builder/utilities.cc @@ -1277,8 +1277,8 @@ namespace WorldBuilder std::vector> mid_oceanic_spreading_velocities, const std::unique_ptr &coordinate_system, const Objects::NaturalCoordinate &position_in_natural_coordinates_at_min_depth, - std::vector> subducting_plate_velocities, - std::vector ridge_migration_times) + const std::vector> &subducting_plate_velocities, + const std::vector &ridge_migration_times) { double distance_ridge = std::numeric_limits::max(); @@ -1333,16 +1333,26 @@ namespace WorldBuilder const double spreading_velocity_point0 = mid_oceanic_spreading_velocities[relevant_ridge][i_coordinate]; const double spreading_velocity_point1 = mid_oceanic_spreading_velocities[relevant_ridge][i_coordinate + 1]; + // When subducting_velocities is not input by the user, default value is 0, which + // results in subducting velocity == spreading_velocity. When a single value is + // input by the user, subducting velocity != spreading_velocity, but + // subducting velocity is spatially constant. double subducting_velocity_point0 = subducting_plate_velocities[0][0]; double subducting_velocity_point1 = subducting_plate_velocities[0][0]; - ridge_migration_time = ridge_migration_times[relevant_ridge]; - - if (subducting_plate_velocities != std::vector> {{0}}) - { - subducting_velocity_point0 = subducting_plate_velocities[relevant_ridge][i_coordinate]; - subducting_velocity_point1 = subducting_plate_velocities[relevant_ridge][i_coordinate + 1]; - } + // When subducting_velocities is input as an array, spatial variation + if (subducting_plate_velocities[0].size() > 1) + { + WBAssertThrow(subducting_plate_velocities.size() == mid_oceanic_ridges.size() && \ + subducting_plate_velocities[relevant_ridge].size() == mid_oceanic_ridges[relevant_ridge].size(), + "subducting velocity and ridge coordinates must be the same dimension"); + WBAssertThrow(ridge_migration_times.size() == mid_oceanic_ridges.size(), + "the times for ridge migration specified in 'spreading velocity' must be the same dimension " + "as ridge coordinates."); + subducting_velocity_point0 = subducting_plate_velocities[relevant_ridge][i_coordinate]; + subducting_velocity_point1 = subducting_plate_velocities[relevant_ridge][i_coordinate + 1]; + ridge_migration_time = ridge_migration_times[relevant_ridge]; + } { // based on http://geomalgorithms.com/a02-_lines.html