From 6c00e94ceddc63c277b3d9ca2416b2c4e566f899 Mon Sep 17 00:00:00 2001 From: danieldouglas92 Date: Fri, 16 Feb 2024 10:20:49 -0600 Subject: [PATCH 1/2] Add variable spreading for ridges --- .../temperature/half_space_model.h | 4 +- .../temperature/plate_model.h | 3 +- .../temperature/mass_conserving.h | 4 +- include/world_builder/parameters.h | 8 +- include/world_builder/types/value_at_points.h | 2 + include/world_builder/utilities.h | 3 +- .../features/continental_plate.cc | 4 +- .../composition/uniform.cc | 4 +- .../grains/random_uniform_distribution.cc | 4 +- .../grains/uniform.cc | 4 +- .../temperature/adiabatic.cc | 4 +- .../temperature/linear.cc | 4 +- .../temperature/uniform.cc | 4 +- source/world_builder/features/mantle_layer.cc | 4 +- .../composition/uniform.cc | 4 +- .../grains/random_uniform_distribution.cc | 4 +- .../mantle_layer_models/grains/uniform.cc | 4 +- .../temperature/adiabatic.cc | 4 +- .../mantle_layer_models/temperature/linear.cc | 4 +- .../temperature/uniform.cc | 4 +- .../world_builder/features/oceanic_plate.cc | 4 +- .../composition/uniform.cc | 4 +- .../grains/random_uniform_distribution.cc | 4 +- .../oceanic_plate_models/grains/uniform.cc | 4 +- .../temperature/adiabatic.cc | 4 +- .../temperature/half_space_model.cc | 34 +++-- .../temperature/linear.cc | 4 +- .../temperature/plate_model.cc | 38 +++-- .../temperature/plate_model_constant_age.cc | 4 +- .../temperature/uniform.cc | 4 +- .../temperature/mass_conserving.cc | 34 ++++- source/world_builder/parameters.cc | 133 ++++++++++++++++++ source/world_builder/types/value_at_points.cc | 8 +- source/world_builder/utilities.cc | 38 +++-- .../gwb-dat/cartesian_variable_spreading.dat | 14 ++ tests/gwb-dat/cartesian_variable_spreading.wb | 17 +++ .../screen-output.log | 9 ++ tests/gwb-dat/screen-output.log | 9 ++ .../gwb-dat/spherical_variable_spreading.dat | 15 ++ tests/gwb-dat/spherical_variable_spreading.wb | 17 +++ .../screen-output.log | 9 ++ tests/unit_tests/unit_test_world_builder.cc | 18 +-- 42 files changed, 409 insertions(+), 96 deletions(-) create mode 100644 tests/gwb-dat/cartesian_variable_spreading.dat create mode 100644 tests/gwb-dat/cartesian_variable_spreading.wb create mode 100644 tests/gwb-dat/cartesian_variable_spreading/screen-output.log create mode 100644 tests/gwb-dat/screen-output.log create mode 100644 tests/gwb-dat/spherical_variable_spreading.dat create mode 100644 tests/gwb-dat/spherical_variable_spreading.wb create mode 100644 tests/gwb-dat/spherical_variable_spreading/screen-output.log diff --git a/include/world_builder/features/oceanic_plate_models/temperature/half_space_model.h b/include/world_builder/features/oceanic_plate_models/temperature/half_space_model.h index 723653d68..519a87ca6 100644 --- a/include/world_builder/features/oceanic_plate_models/temperature/half_space_model.h +++ b/include/world_builder/features/oceanic_plate_models/temperature/half_space_model.h @@ -81,7 +81,6 @@ namespace WorldBuilder const double feature_min_depth, const double feature_max_depth) const override final; - private: // plate model temperature submodule parameters double min_depth; @@ -90,8 +89,9 @@ namespace WorldBuilder Objects::Surface max_depth_surface; double top_temperature; double bottom_temperature; - double spreading_velocity; + std::pair,std::vector> spreading_velocities; std::vector > > mid_oceanic_ridges; + std::vector> spreading_velocities_at_each_ridge_point; Operations operation; }; diff --git a/include/world_builder/features/oceanic_plate_models/temperature/plate_model.h b/include/world_builder/features/oceanic_plate_models/temperature/plate_model.h index 233a1b846..d36cc14aa 100644 --- a/include/world_builder/features/oceanic_plate_models/temperature/plate_model.h +++ b/include/world_builder/features/oceanic_plate_models/temperature/plate_model.h @@ -93,8 +93,9 @@ namespace WorldBuilder Objects::Surface max_depth_surface; double top_temperature; double bottom_temperature; - double spreading_velocity; + std::pair,std::vector> spreading_velocities; std::vector > > mid_oceanic_ridges; + std::vector> spreading_velocities_at_each_ridge_point; Operations operation; }; diff --git a/include/world_builder/features/subducting_plate_models/temperature/mass_conserving.h b/include/world_builder/features/subducting_plate_models/temperature/mass_conserving.h index 22247fd2c..c0f260971 100644 --- a/include/world_builder/features/subducting_plate_models/temperature/mass_conserving.h +++ b/include/world_builder/features/subducting_plate_models/temperature/mass_conserving.h @@ -100,6 +100,7 @@ namespace WorldBuilder const double min_temperature, const double background_temperature, const double temperature_, + const double plate_velocity, const double effective_plate_age, const double adjusted_distance) const; @@ -109,7 +110,8 @@ namespace WorldBuilder double min_depth; double max_depth; double density; - double plate_velocity; + std::pair,std::vector> plate_velocities; + std::vector> plate_velocities_at_each_ridge_point; double mantle_coupling_depth; double forearc_cooling_factor; double thermal_conductivity; diff --git a/include/world_builder/parameters.h b/include/world_builder/parameters.h index ef290f826..7c8df6e33 100644 --- a/include/world_builder/parameters.h +++ b/include/world_builder/parameters.h @@ -123,7 +123,13 @@ namespace WorldBuilder const std::vector > &addition_points = {}); /** - * A specialized version of get which can return vectors/arrays. + * A specialized version of get which can return a values at times type. + * \param name The name of the entry to retrieved + */ + std::pair,std::vector> get_value_at_array(const std::string &name); + + /** + * A specialized verions of get which can return vectors/arrays. * This version is designed for the plugin system. * \param name The name of the entry to retrieved */ diff --git a/include/world_builder/types/value_at_points.h b/include/world_builder/types/value_at_points.h index e55139ac7..3055c3802 100644 --- a/include/world_builder/types/value_at_points.h +++ b/include/world_builder/types/value_at_points.h @@ -40,6 +40,7 @@ namespace WorldBuilder * A constructor */ ValueAtPoints(const double default_value, + size_t max_values_in_array, std::vector> default_points_ = std::vector>()); /** @@ -61,6 +62,7 @@ namespace WorldBuilder const std::string &documentation) const override final; double default_value; + double max_values_in_array; std::vector > default_points; protected: diff --git a/include/world_builder/utilities.h b/include/world_builder/utilities.h index 77544ad59..fafafdac7 100644 --- a/include/world_builder/utilities.h +++ b/include/world_builder/utilities.h @@ -440,9 +440,10 @@ namespace WorldBuilder */ std::pair calculate_ridge_distance_and_spreading(std::vector>> mid_oceanic_ridges, - const double spreading_velocity, + std::vector> mid_oceanic_spreading_velocities, const std::unique_ptr &coordinate_system, const Objects::NaturalCoordinate &position_in_natural_coordinates_at_min_depth); + } // namespace Utilities } // namespace WorldBuilder diff --git a/source/world_builder/features/continental_plate.cc b/source/world_builder/features/continental_plate.cc index d5fb47a3d..09fa5d48f 100644 --- a/source/world_builder/features/continental_plate.cc +++ b/source/world_builder/features/continental_plate.cc @@ -64,9 +64,9 @@ namespace WorldBuilder { prm.declare_entry("", Types::Object(required_entries), "Continental plate object. Requires properties `model` and `coordinates`."); - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth from which this feature is present"); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth to which this feature is present"); prm.declare_entry("temperature models", Types::PluginSystem("", Features::ContinentalPlateModels::Temperature::Interface::declare_entries, {"model"}), diff --git a/source/world_builder/features/continental_plate_models/composition/uniform.cc b/source/world_builder/features/continental_plate_models/composition/uniform.cc index 93fbf6209..6913d5436 100644 --- a/source/world_builder/features/continental_plate_models/composition/uniform.cc +++ b/source/world_builder/features/continental_plate_models/composition/uniform.cc @@ -65,10 +65,10 @@ namespace WorldBuilder "Uniform compositional model. Sets constant compositional field."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the composition of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the composition of this feature is present."); prm.declare_entry("compositions", Types::Array(Types::UnsignedInt(),0), diff --git a/source/world_builder/features/continental_plate_models/grains/random_uniform_distribution.cc b/source/world_builder/features/continental_plate_models/grains/random_uniform_distribution.cc index 115b436f4..7db823110 100644 --- a/source/world_builder/features/continental_plate_models/grains/random_uniform_distribution.cc +++ b/source/world_builder/features/continental_plate_models/grains/random_uniform_distribution.cc @@ -67,9 +67,9 @@ namespace WorldBuilder "to a single value or to a random distribution."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the composition of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the composition of this feature is present."); prm.declare_entry("compositions", Types::Array(Types::UnsignedInt(),0), diff --git a/source/world_builder/features/continental_plate_models/grains/uniform.cc b/source/world_builder/features/continental_plate_models/grains/uniform.cc index 3ae4857c1..78c8b7310 100644 --- a/source/world_builder/features/continental_plate_models/grains/uniform.cc +++ b/source/world_builder/features/continental_plate_models/grains/uniform.cc @@ -64,10 +64,10 @@ namespace WorldBuilder "Uniform grains model. All grains start exactly the same."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the composition of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the composition of this feature is present."); prm.declare_entry("compositions", Types::Array(Types::UnsignedInt(),0), diff --git a/source/world_builder/features/continental_plate_models/temperature/adiabatic.cc b/source/world_builder/features/continental_plate_models/temperature/adiabatic.cc index db97ae853..70cd83344 100644 --- a/source/world_builder/features/continental_plate_models/temperature/adiabatic.cc +++ b/source/world_builder/features/continental_plate_models/temperature/adiabatic.cc @@ -66,10 +66,10 @@ namespace WorldBuilder "Adiabatic temperature model. Uses global values by default."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the composition of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the composition of this feature is present."); prm.declare_entry("potential mantle temperature", Types::Double(-1), diff --git a/source/world_builder/features/continental_plate_models/temperature/linear.cc b/source/world_builder/features/continental_plate_models/temperature/linear.cc index bedbef2c5..9244eb5ac 100644 --- a/source/world_builder/features/continental_plate_models/temperature/linear.cc +++ b/source/world_builder/features/continental_plate_models/temperature/linear.cc @@ -66,10 +66,10 @@ namespace WorldBuilder "Linear temperature model. Can be set to use an adiabatic temperature at the boundaries."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the composition of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the composition of this feature is present."); prm.declare_entry("top temperature", Types::Double(293.15), diff --git a/source/world_builder/features/continental_plate_models/temperature/uniform.cc b/source/world_builder/features/continental_plate_models/temperature/uniform.cc index cd65ac46f..d186a4dd3 100644 --- a/source/world_builder/features/continental_plate_models/temperature/uniform.cc +++ b/source/world_builder/features/continental_plate_models/temperature/uniform.cc @@ -64,10 +64,10 @@ namespace WorldBuilder "Uniform temperature model. Set the temperature to a constan value."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the composition of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the composition of this feature is present."); prm.declare_entry("temperature", Types::Double(293.15), diff --git a/source/world_builder/features/mantle_layer.cc b/source/world_builder/features/mantle_layer.cc index 4d237c330..7ac9b9520 100644 --- a/source/world_builder/features/mantle_layer.cc +++ b/source/world_builder/features/mantle_layer.cc @@ -58,9 +58,9 @@ namespace WorldBuilder { prm.declare_entry("", Types::Object(required_entries), "Mantle layer object. Requires properties `model` and `coordinates`."); - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth from which this feature is present"); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth to which this feature is present"); prm.declare_entry("temperature models", Types::PluginSystem("", Features::MantleLayerModels::Temperature::Interface::declare_entries, {"model"}), diff --git a/source/world_builder/features/mantle_layer_models/composition/uniform.cc b/source/world_builder/features/mantle_layer_models/composition/uniform.cc index d3828bb9c..dd137ec39 100644 --- a/source/world_builder/features/mantle_layer_models/composition/uniform.cc +++ b/source/world_builder/features/mantle_layer_models/composition/uniform.cc @@ -62,9 +62,9 @@ namespace WorldBuilder "Uniform compositional model. Sets constant compositional field."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the composition of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the composition of this feature is present."); prm.declare_entry("compositions", Types::Array(Types::UnsignedInt(),0), "A list with the labels of the composition which are present there."); diff --git a/source/world_builder/features/mantle_layer_models/grains/random_uniform_distribution.cc b/source/world_builder/features/mantle_layer_models/grains/random_uniform_distribution.cc index 01ed7a346..a5d139840 100644 --- a/source/world_builder/features/mantle_layer_models/grains/random_uniform_distribution.cc +++ b/source/world_builder/features/mantle_layer_models/grains/random_uniform_distribution.cc @@ -65,9 +65,9 @@ namespace WorldBuilder "to a single value or to a random distribution."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the composition of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the composition of this feature is present."); prm.declare_entry("compositions", Types::Array(Types::UnsignedInt(),0), diff --git a/source/world_builder/features/mantle_layer_models/grains/uniform.cc b/source/world_builder/features/mantle_layer_models/grains/uniform.cc index bd2252d85..6bb45713d 100644 --- a/source/world_builder/features/mantle_layer_models/grains/uniform.cc +++ b/source/world_builder/features/mantle_layer_models/grains/uniform.cc @@ -63,9 +63,9 @@ namespace WorldBuilder "Uniform grains model. All grains start exactly the same."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the composition of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the composition of this feature is present."); prm.declare_entry("compositions", Types::Array(Types::UnsignedInt(),0), diff --git a/source/world_builder/features/mantle_layer_models/temperature/adiabatic.cc b/source/world_builder/features/mantle_layer_models/temperature/adiabatic.cc index 026bda854..5b68f5d57 100644 --- a/source/world_builder/features/mantle_layer_models/temperature/adiabatic.cc +++ b/source/world_builder/features/mantle_layer_models/temperature/adiabatic.cc @@ -65,10 +65,10 @@ namespace WorldBuilder "Adiabatic temperature model. Uses global values by default."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the temperature of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the temperature of this feature is present."); prm.declare_entry("potential mantle temperature", Types::Double(-1), diff --git a/source/world_builder/features/mantle_layer_models/temperature/linear.cc b/source/world_builder/features/mantle_layer_models/temperature/linear.cc index db3ad32c6..7d9668ba9 100644 --- a/source/world_builder/features/mantle_layer_models/temperature/linear.cc +++ b/source/world_builder/features/mantle_layer_models/temperature/linear.cc @@ -64,10 +64,10 @@ namespace WorldBuilder "Linear temperature model. Can be set to use an adiabatic temperature at the boundaries."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the temperature of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the temperature of this feature is present."); prm.declare_entry("top temperature", Types::Double(293.15), diff --git a/source/world_builder/features/mantle_layer_models/temperature/uniform.cc b/source/world_builder/features/mantle_layer_models/temperature/uniform.cc index 2e47d23d3..3abf18e21 100644 --- a/source/world_builder/features/mantle_layer_models/temperature/uniform.cc +++ b/source/world_builder/features/mantle_layer_models/temperature/uniform.cc @@ -63,10 +63,10 @@ namespace WorldBuilder "Uniform temperature model. Set the temperature to a constan value."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the temperature of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the temperature of this feature is present."); prm.declare_entry("temperature", Types::Double(293.15), diff --git a/source/world_builder/features/oceanic_plate.cc b/source/world_builder/features/oceanic_plate.cc index a88aa13b6..9996efe35 100644 --- a/source/world_builder/features/oceanic_plate.cc +++ b/source/world_builder/features/oceanic_plate.cc @@ -60,9 +60,9 @@ namespace WorldBuilder { prm.declare_entry("", Types::Object(required_entries), "Oceanic plate object. Requires properties `model` and `coordinates`."); - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth from which this feature is present"); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth to which this feature is present"); prm.declare_entry("temperature models", Types::PluginSystem("", Features::OceanicPlateModels::Temperature::Interface::declare_entries, {"model"}), diff --git a/source/world_builder/features/oceanic_plate_models/composition/uniform.cc b/source/world_builder/features/oceanic_plate_models/composition/uniform.cc index 62665a6f4..722d056c2 100644 --- a/source/world_builder/features/oceanic_plate_models/composition/uniform.cc +++ b/source/world_builder/features/oceanic_plate_models/composition/uniform.cc @@ -63,9 +63,9 @@ namespace WorldBuilder "Uniform compositional model. Sets constant compositional field."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the composition of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the composition of this feature is present."); prm.declare_entry("compositions", Types::Array(Types::UnsignedInt(),0), "A list with the labels of the composition which are present there."); diff --git a/source/world_builder/features/oceanic_plate_models/grains/random_uniform_distribution.cc b/source/world_builder/features/oceanic_plate_models/grains/random_uniform_distribution.cc index 0302d4b28..ae3de62fc 100644 --- a/source/world_builder/features/oceanic_plate_models/grains/random_uniform_distribution.cc +++ b/source/world_builder/features/oceanic_plate_models/grains/random_uniform_distribution.cc @@ -66,9 +66,9 @@ namespace WorldBuilder "to a single value or to a random distribution."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the composition of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the composition of this feature is present."); prm.declare_entry("compositions", Types::Array(Types::UnsignedInt(),0), diff --git a/source/world_builder/features/oceanic_plate_models/grains/uniform.cc b/source/world_builder/features/oceanic_plate_models/grains/uniform.cc index f9c2a6d4b..c94783764 100644 --- a/source/world_builder/features/oceanic_plate_models/grains/uniform.cc +++ b/source/world_builder/features/oceanic_plate_models/grains/uniform.cc @@ -64,9 +64,9 @@ namespace WorldBuilder "Uniform grains model. All grains start exactly the same."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the composition of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the composition of this feature is present."); prm.declare_entry("compositions", Types::Array(Types::UnsignedInt(),0), diff --git a/source/world_builder/features/oceanic_plate_models/temperature/adiabatic.cc b/source/world_builder/features/oceanic_plate_models/temperature/adiabatic.cc index 1c6dcd73f..b54e5f71e 100644 --- a/source/world_builder/features/oceanic_plate_models/temperature/adiabatic.cc +++ b/source/world_builder/features/oceanic_plate_models/temperature/adiabatic.cc @@ -65,10 +65,10 @@ namespace WorldBuilder "Adiabatic temperature model. Uses global values by default."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the temperature of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the temperature of this feature is present."); prm.declare_entry("potential mantle temperature", Types::Double(-1), 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 60ee31597..ff13a7107 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 @@ -49,7 +49,6 @@ namespace WorldBuilder max_depth(NaN::DSNAN), top_temperature(NaN::DSNAN), bottom_temperature(NaN::DSNAN), - spreading_velocity(NaN::DSNAN), operation(Operations::REPLACE) { this->world = world_; @@ -68,10 +67,10 @@ namespace WorldBuilder "Half space cooling mode"); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the temperature of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the temperature of this feature is present." "Because half-space reaches background temperature asymptotically, " "this value should be ~2 times the nominal plate thickness of 100 km" ); @@ -84,7 +83,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::Double(-1), + prm.declare_entry("spreading velocity", Types::OneOf(Types::Double(0.01),Types::Array(Types::ValueAtPoints(0.01, 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."); @@ -106,7 +105,7 @@ namespace WorldBuilder operation = string_operations_to_enum(prm.get("operation")); top_temperature = prm.get("top temperature"); bottom_temperature = prm.get("bottom temperature"); - spreading_velocity = prm.get("spreading velocity")/31557600; // m/seconds + spreading_velocities = prm.get_value_at_array("spreading velocity"); mid_oceanic_ridges = prm.get_vector>>("ridge coordinates"); const double dtr = prm.coordinate_system->natural_coordinate_system() == spherical ? Consts::PI / 180.0 : 1.0; @@ -115,8 +114,24 @@ namespace WorldBuilder { ridge_coordinate *= dtr; } - } + unsigned int index_x = 0; + unsigned int index_y = 0; + unsigned int ridge_point_index = 0; + for (index_x = 0; index_x < mid_oceanic_ridges.size(); index_x++) + { + std::vector spreading_rates_for_ridge; + for (index_y = 0; index_y < mid_oceanic_ridges[index_x].size(); index_y++) + { + if (spreading_velocities.second.size() <= 1) + spreading_rates_for_ridge.push_back(spreading_velocities.first[0]); + else + spreading_rates_for_ridge.push_back(spreading_velocities.second[ridge_point_index]); + ridge_point_index += 1; + } + spreading_velocities_at_each_ridge_point.push_back(spreading_rates_for_ridge); + } + } double HalfSpaceModel::get_temperature(const Point<3> &position, @@ -137,7 +152,6 @@ namespace WorldBuilder *(world->parameters.coordinate_system)); position_in_natural_coordinates_at_min_depth.get_ref_depth_coordinate() += depth-min_depth; - double bottom_temperature_local = bottom_temperature; if (bottom_temperature_local < 0) @@ -148,7 +162,7 @@ namespace WorldBuilder } std::pair ridge_parameters = Utilities::calculate_ridge_distance_and_spreading(mid_oceanic_ridges, - spreading_velocity, + spreading_velocities_at_each_ridge_point, world->parameters.coordinate_system, position_in_natural_coordinates_at_min_depth); @@ -165,12 +179,12 @@ namespace WorldBuilder << ". Relevant variables: bottom_temperature_local = " << bottom_temperature_local << ", top_temperature = " << top_temperature << ", max_depth = " << max_depth - << ", spreading_velocity = " << spreading_velocity + << ", spreading_velocity = " << ridge_parameters.first << ", thermal_diffusivity = " << thermal_diffusivity << ", age = " << age << '.'); WBAssert(std::isfinite(temperature), "Temperature inside half-space cooling model is not a finite: " << temperature << ". Relevant variables: bottom_temperature_local = " << bottom_temperature_local << ", top_temperature = " << top_temperature - << ", spreading_velocity = " << spreading_velocity + << ", spreading_velocity = " << ridge_parameters.first << ", thermal_diffusivity = " << thermal_diffusivity << ", age = " << age << '.'); diff --git a/source/world_builder/features/oceanic_plate_models/temperature/linear.cc b/source/world_builder/features/oceanic_plate_models/temperature/linear.cc index f339e9785..51642cd6f 100644 --- a/source/world_builder/features/oceanic_plate_models/temperature/linear.cc +++ b/source/world_builder/features/oceanic_plate_models/temperature/linear.cc @@ -65,10 +65,10 @@ namespace WorldBuilder "Linear temperature model. Can be set to use an adiabatic temperature at the boundaries."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the temperature of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the temperature of this feature is present."); prm.declare_entry("top temperature", Types::Double(293.15), 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 c97b063cc..062ad0de2 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 @@ -48,7 +48,6 @@ namespace WorldBuilder max_depth(NaN::DSNAN), top_temperature(NaN::DSNAN), bottom_temperature(NaN::DSNAN), - spreading_velocity(NaN::DSNAN), operation(Operations::REPLACE) { this->world = world_; @@ -67,10 +66,10 @@ namespace WorldBuilder "Plate model."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the temperature of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the temperature of this feature is present."); prm.declare_entry("top temperature", Types::Double(293.15), @@ -79,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::Double(-1), + prm.declare_entry("spreading velocity", Types::OneOf(Types::Double(0.01),Types::Array(Types::ValueAtPoints(0.01, 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."); @@ -102,7 +101,7 @@ namespace WorldBuilder operation = string_operations_to_enum(prm.get("operation")); top_temperature = prm.get("top temperature"); bottom_temperature = prm.get("bottom temperature"); - spreading_velocity = prm.get("spreading velocity")/31557600; + spreading_velocities = prm.get_value_at_array("spreading velocity"); mid_oceanic_ridges = prm.get_vector>>("ridge coordinates"); const double dtr = prm.coordinate_system->natural_coordinate_system() == spherical ? Consts::PI / 180.0 : 1.0; @@ -111,6 +110,23 @@ namespace WorldBuilder { ridge_coordinate *= dtr; } + + unsigned int index_x = 0; + unsigned int index_y = 0; + unsigned int ridge_point_index = 0; + for (index_x = 0; index_x < mid_oceanic_ridges.size(); index_x++) + { + std::vector spreading_rates_for_ridge; + for (index_y = 0; index_y < mid_oceanic_ridges[index_x].size(); index_y++) + { + if (spreading_velocities.second.size() <= 1) + spreading_rates_for_ridge.push_back(spreading_velocities.first[0]); + else + spreading_rates_for_ridge.push_back(spreading_velocities.second[ridge_point_index]); + ridge_point_index += 1; + } + spreading_velocities_at_each_ridge_point.push_back(spreading_rates_for_ridge); + } } @@ -145,7 +161,7 @@ namespace WorldBuilder const int sommation_number = 100; std::pair ridge_parameters = Utilities::calculate_ridge_distance_and_spreading(mid_oceanic_ridges, - spreading_velocity, + spreading_velocities_at_each_ridge_point, world->parameters.coordinate_system, position_in_natural_coordinates_at_min_depth); @@ -159,10 +175,10 @@ namespace WorldBuilder { temperature = temperature + (bottom_temperature_local - top_temperature) * ((2 / (double(i) * Consts::PI)) * std::sin((double(i) * Consts::PI * depth) / max_depth) * - std::exp((((spreading_velocity * max_depth)/(2 * thermal_diffusivity)) - - std::sqrt(((spreading_velocity*spreading_velocity*max_depth*max_depth) / + std::exp((((ridge_parameters.first * max_depth)/(2 * thermal_diffusivity)) - + std::sqrt(((ridge_parameters.first*ridge_parameters.first*max_depth*max_depth) / (4*thermal_diffusivity*thermal_diffusivity)) + double(i) * double(i) * Consts::PI * Consts::PI)) * - ((spreading_velocity * age) / max_depth))); + ((ridge_parameters.first * age) / max_depth))); } @@ -170,12 +186,12 @@ namespace WorldBuilder << ". Relevant variables: bottom_temperature_local = " << bottom_temperature_local << ", top_temperature = " << top_temperature << ", max_depth = " << max_depth - << ", spreading_velocity = " << spreading_velocity + << ", spreading_velocity = " << ridge_parameters.first << ", thermal_diffusivity = " << thermal_diffusivity << ", age = " << age << '.'); WBAssert(std::isfinite(temperature), "Temparture inside plate model is not a finite: " << temperature << ". Relevant variables: bottom_temperature_local = " << bottom_temperature_local << ", top_temperature = " << top_temperature - << ", spreading_velocity = " << spreading_velocity + << ", spreading_velocity = " << ridge_parameters.first << ", thermal_diffusivity = " << thermal_diffusivity << ", age = " << age << '.'); diff --git a/source/world_builder/features/oceanic_plate_models/temperature/plate_model_constant_age.cc b/source/world_builder/features/oceanic_plate_models/temperature/plate_model_constant_age.cc index 845eae8cc..81040d2b7 100644 --- a/source/world_builder/features/oceanic_plate_models/temperature/plate_model_constant_age.cc +++ b/source/world_builder/features/oceanic_plate_models/temperature/plate_model_constant_age.cc @@ -66,10 +66,10 @@ namespace WorldBuilder "Plate model, but with a fixed age."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the temperature of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the temperature of this feature is present."); prm.declare_entry("top temperature", Types::Double(293.15), diff --git a/source/world_builder/features/oceanic_plate_models/temperature/uniform.cc b/source/world_builder/features/oceanic_plate_models/temperature/uniform.cc index 9ceace245..255828256 100644 --- a/source/world_builder/features/oceanic_plate_models/temperature/uniform.cc +++ b/source/world_builder/features/oceanic_plate_models/temperature/uniform.cc @@ -64,10 +64,10 @@ namespace WorldBuilder "Uniform temperature model. Set the temperature to a constan value."); // Declare entries of this plugin - prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0.))), + prm.declare_entry("min depth", Types::OneOf(Types::Double(0),Types::Array(Types::ValueAtPoints(0., 2.))), "The depth in meters from which the temperature of this feature is present."); - prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max()))), + prm.declare_entry("max depth", Types::OneOf(Types::Double(std::numeric_limits::max()),Types::Array(Types::ValueAtPoints(std::numeric_limits::max(), 2.))), "The depth in meters to which the temperature of this feature is present."); prm.declare_entry("temperature", Types::Double(293.15), 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 3b8cd0ac9..5667df28a 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 @@ -28,8 +28,10 @@ #include "world_builder/types/array.h" #include "world_builder/types/bool.h" #include "world_builder/types/double.h" +#include "world_builder/types/one_of.h" #include "world_builder/types/object.h" #include "world_builder/types/point.h" +#include "world_builder/types/value_at_points.h" #include "world_builder/utilities.h" #include "world_builder/world.h" @@ -48,7 +50,6 @@ namespace WorldBuilder min_depth(NaN::DSNAN), max_depth(NaN::DSNAN), density(NaN::DSNAN), - plate_velocity(NaN::DSNAN), mantle_coupling_depth(NaN::DSNAN), forearc_cooling_factor(NaN::DSNAN), thermal_conductivity(NaN::DSNAN), @@ -116,7 +117,7 @@ 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::Double(0.05), + prm.declare_entry("plate velocity", Types::OneOf(Types::Double(0.01),Types::Array(Types::ValueAtPoints(0.01, std::numeric_limits::max()))), "The velocity with which the plate subducts in meters per year. Default is 5 cm/yr"); prm.declare_entry("coupling depth", Types::Double(100e3), @@ -180,7 +181,7 @@ namespace WorldBuilder density = prm.get("density"); thermal_conductivity = prm.get("thermal conductivity"); - plate_velocity = prm.get("plate velocity"); + plate_velocities = prm.get_value_at_array("plate velocity"); mantle_coupling_depth = prm.get("coupling depth"); forearc_cooling_factor = prm.get("forearc cooling factor"); @@ -216,6 +217,24 @@ namespace WorldBuilder { ridge_coordinate *= dtr; } + + unsigned int index_x = 0; + unsigned int index_y = 0; + unsigned int ridge_point_index = 0; + for (index_x = 0; index_x < mid_oceanic_ridges.size(); index_x++) + { + std::vector plate_velocities_for_ridge; + for (index_y = 0; index_y < mid_oceanic_ridges[index_x].size(); index_y++) + { + if (plate_velocities.second.size() <= 1) + plate_velocities_for_ridge.push_back(plate_velocities.first[0]); + else + plate_velocities_for_ridge.push_back(plate_velocities.second[ridge_point_index]); + ridge_point_index += 1; + } + plate_velocities_at_each_ridge_point.push_back(plate_velocities_for_ridge); + } + std::string reference_model_name_str = prm.get("reference model name"); if (reference_model_name_str=="plate model") reference_model_name = plate_model; @@ -235,6 +254,8 @@ namespace WorldBuilder { const double distance_from_plane = distance_from_planes.distance_from_plane; + const int plate_model_summation_number = 100; // for the plate model + if (distance_from_plane <= max_depth && distance_from_plane >= min_depth) { @@ -242,8 +263,9 @@ namespace WorldBuilder const Point<3> trench_point = distance_from_planes.closest_trench_point; const Objects::NaturalCoordinate trench_point_natural = Objects::NaturalCoordinate(trench_point, *(world->parameters.coordinate_system)); + std::pair ridge_parameters = Utilities::calculate_ridge_distance_and_spreading(mid_oceanic_ridges, - plate_velocity, + plate_velocities_at_each_ridge_point, world->parameters.coordinate_system, trench_point_natural); @@ -251,6 +273,7 @@ namespace WorldBuilder const double cm2m = 100; // 100 cm/m const double my = 1.0e6; // 1e6 y/my const double seconds_in_year = 60.0 * 60.0 * 24.0 * 365.25; // sec/y + const double plate_velocity = ridge_parameters.first * seconds_in_year; // m/yr const double age_at_trench = ridge_parameters.second / plate_velocity; // m/(m/y) = yr const double plate_age_sec = age_at_trench * seconds_in_year; // y --> seconds @@ -455,7 +478,7 @@ namespace WorldBuilder } // Call the analytic solution to compute the temperature - temperature = get_temperature_analytic(top_heat_content, min_temperature, background_temperature, temperature_, effective_plate_age, adjusted_distance); + temperature = get_temperature_analytic(top_heat_content, min_temperature, background_temperature, temperature_, plate_velocity, effective_plate_age, adjusted_distance); } else { @@ -477,6 +500,7 @@ namespace WorldBuilder const double min_temperature, const double background_temperature, const double temperature_, + const double plate_velocity, const double effective_plate_age, const double adjusted_distance) const { diff --git a/source/world_builder/parameters.cc b/source/world_builder/parameters.cc index 8eaac431b..920cd5c21 100644 --- a/source/world_builder/parameters.cc +++ b/source/world_builder/parameters.cc @@ -651,6 +651,139 @@ namespace WorldBuilder return result; } + + std::pair,std::vector> Parameters::get_value_at_array(const std::string &name) + { + // There are two cases: + // 1. One double provided: use the default value everywhere. Return first with one value and second with size 0. + // 2. Other: fill the vectors with the default value and addition points and then add new point. + // If a value without points is encountered, the additional points are used. + std::pair,std::vector> result; + + const std::string strict_base = this->get_full_json_path(); + + // start with adding the additional points with the default value + // to do this we need the default value + double default_value = 0; + bool is_array = true; + if (Pointer((strict_base + "/" + name).c_str()).Get(parameters) != nullptr && Pointer((strict_base + "/" + name).c_str()).Get(parameters)->IsArray()) + { + const std::string value_def_path = get_full_json_schema_path() + "/" + name + "/oneOf/1/items/items/anyOf/0/default value"; + Value *value_def = Pointer(value_def_path.c_str()).Get(declarations); + WBAssertThrow(value_def != nullptr, + "internal error: could not retrieve the default value at: " + << value_def_path); + + // Since the default value is set in the code, if it fails it is an internal error, not a user error. + // So no try/catch needed. + default_value = value_def->GetDouble(); + } + else + { + is_array = false; + Value *value_def = Pointer((get_full_json_schema_path() + "/" + name + "/oneOf/0/default value").c_str()).Get(declarations); + WBAssertThrow(value_def != nullptr, + "internal error: could not retrieve the default value at: " + <GetDouble(); + } + + + // check if there is a user defined value + if (Pointer((strict_base + "/" + name).c_str()).Get(parameters) != nullptr) + { + // there is a user defined value, so either case 2, 3 or 4. + if (is_array) + { + Value *array = Pointer((strict_base + "/" + name).c_str()).Get(parameters); + + // go through all the points in order + for (size_t i = 0; i < array->Size(); ++i ) + { + // now parse a single value_at_point. + const std::string base = (strict_base + "/").append(name).append("/").append(std::to_string(i)); + // Let's assume that the file is correct, because it has been checked with the json schema. + // So there are exactly two values, the first value is a double, the second an array of 2d arrays (points). + + // Get the double + double value; + Value *value_pointer = Pointer((base + "/0").c_str()).Get(parameters); + + WBAssertThrow(value_pointer != nullptr, "internal error: this should not happen."); + + try + { + value = value_pointer->GetDouble(); + } + catch (...) + { + WBAssertThrow(false, "Could not convert values of " << base << "/0 into doubles. " + << "The provided value was \"" << Pointer((base + "/0").c_str()).Get(parameters)->GetString() << "\"."); + } + + // now get the array of points. + Value *coordinates_array = Pointer((base + "/1").c_str()).Get(parameters); + double testing_value; + if (coordinates_array != nullptr) + { + for (size_t coordinate_i = 0; coordinate_i < coordinates_array->Size(); ++coordinate_i ) + { + Value *coordinate_j_array = Pointer((base + "/1/" + std::to_string(coordinate_i)).c_str()).Get(parameters); + for (size_t coordinate_j = 0; coordinate_j < coordinate_j_array->Size(); ++coordinate_j) + { + // Let's assume that the file is correct, because it has been checked with the json schema. + // 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(); + } + catch (...) + { + WBAssertThrow(false, "Could not convert values of " << base + "/1/" + std::to_string(coordinate_i) + "/" + std::to_string(coordinate_j) + << " into a Point<2> array, because it could not convert the 1st sub-elements into doubles. " + << "The provided value was \"" + << Pointer((base + "/1/" + std::to_string(coordinate_i) + "/" + std::to_string(coordinate_j)).c_str()).Get(parameters)->GetString() + << "\"."); + } + + result.first.emplace_back(value); + result.second.emplace_back(testing_value); + } + } + + } + } + } + else + { + // case 1: there one value, not an array + double value = 0; + try + { + value = Pointer((strict_base + "/" + name).c_str()).Get(parameters)->GetDouble(); + } + catch (...) + { + WBAssertThrow(false, "Could not convert values of " << strict_base << "/" << name << " into a double. " + << "The provided value was \"" << Pointer((strict_base + "/" + name).c_str()).Get(parameters)->GetString() << "\"."); + } + result.first.emplace_back(value); + } + } + else + { + // there is no user defined value. Case one: return the default value and no points + result.first.emplace_back(default_value); + } + + return result; + } + + template<> std::vector > Parameters::get_vector(const std::string &name) diff --git a/source/world_builder/types/value_at_points.cc b/source/world_builder/types/value_at_points.cc index d36b5375a..d7eedd37d 100644 --- a/source/world_builder/types/value_at_points.cc +++ b/source/world_builder/types/value_at_points.cc @@ -29,9 +29,11 @@ namespace WorldBuilder { ValueAtPoints::ValueAtPoints(const double default_value_, + size_t max_values_in_array_, std::vector> default_points_) : default_value(default_value_), + max_values_in_array(max_values_in_array_), default_points(std::move(std::move(default_points_))) { this->type_name = Types::type::ValueAtPoints; @@ -42,6 +44,7 @@ namespace WorldBuilder ValueAtPoints::ValueAtPoints(ValueAtPoints const &other) : default_value(other.default_value), + max_values_in_array(other.max_values_in_array), default_points(other.default_points) { this->type_name = Types::type::ValueAtPoints; @@ -52,7 +55,6 @@ namespace WorldBuilder = default; - void ValueAtPoints::write_schema(Parameters &prm, const std::string &name, @@ -67,7 +69,7 @@ namespace WorldBuilder Pointer((base + "/type").c_str()).Set(declarations,"array"); Pointer((base + "/additionalProperties").c_str()).Set(declarations,false); Pointer((base + "/minItems").c_str()).Set(declarations,1); - Pointer((base + "/maxItems").c_str()).Set(declarations,2); + Pointer((base + "/maxItems").c_str()).Set(declarations, max_values_in_array); Pointer((base + "/description").c_str()).Set(declarations,documentation.c_str()); { @@ -80,7 +82,7 @@ namespace WorldBuilder Pointer((base + "/items/anyOf/1/items/type").c_str()).Set(declarations,"array"); Pointer((base + "/items/anyOf/1/items/minItems").c_str()).Set(declarations,1); - Pointer((base + "/items/anyOf/1/items/maxItems").c_str()).Set(declarations,2); + Pointer((base + "/items/anyOf/1/items/maxItems").c_str()).Set(declarations, max_values_in_array); Pointer((base + "/items/anyOf/1/items/items/type").c_str()).Set(declarations,"number"); } diff --git a/source/world_builder/utilities.cc b/source/world_builder/utilities.cc index 8b236cecf..d32be28ae 100644 --- a/source/world_builder/utilities.cc +++ b/source/world_builder/utilities.cc @@ -1223,12 +1223,13 @@ namespace WorldBuilder std::pair calculate_ridge_distance_and_spreading(std::vector>> mid_oceanic_ridges, - const double spreading_velocity, + std::vector> mid_oceanic_spreading_velocities, const std::unique_ptr &coordinate_system, const Objects::NaturalCoordinate &position_in_natural_coordinates_at_min_depth) { double distance_ridge = std::numeric_limits::max(); + double spreading_velocity_at_ridge = 0; // first find if the coordinate is on this side of a ridge unsigned int relevant_ridge = 0; @@ -1277,6 +1278,9 @@ namespace WorldBuilder const Point<2> segment_point0 = mid_oceanic_ridges[relevant_ridge][i_coordinate]; const Point<2> segment_point1 = mid_oceanic_ridges[relevant_ridge][i_coordinate + 1]; + 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]; + { // based on http://geomalgorithms.com/a02-_lines.html const Point<2> v = segment_point1 - segment_point0; @@ -1293,19 +1297,37 @@ namespace WorldBuilder // If you want to have infinite lines, use only the else statement. if (c1 <= 0) - Pb1=segment_point0; + { + Pb1=segment_point0; + spreading_velocity_at_ridge = spreading_velocity_point0; + } else if (c <= c1) - Pb1=segment_point1; + { + Pb1=segment_point1; + spreading_velocity_at_ridge = spreading_velocity_point1; + } else - Pb1=segment_point0 + (c1 / c) * v; + { + Pb1=segment_point0 + (c1 / c) * v; + spreading_velocity_at_ridge = spreading_velocity_point0 + (spreading_velocity_point1 - spreading_velocity_point0) * (c1 / c) * 1; + } Point<2> Pb2(coordinate_system->natural_coordinate_system()); if (c2 <= 0) - Pb2=segment_point0; + { + Pb2=segment_point0; + spreading_velocity_at_ridge = spreading_velocity_point0; + } else if (c <= c2) - Pb2=segment_point1; + { + Pb2=segment_point1; + spreading_velocity_at_ridge = spreading_velocity_point1; + } else - Pb2=segment_point0 + (c2 / c) * v; + { + Pb2=segment_point0 + (c1 / c) * v; + spreading_velocity_at_ridge = spreading_velocity_point0 + (spreading_velocity_point1 - spreading_velocity_point0) * (c1 / c) * 1; + } Point<3> compare_point1(coordinate_system->natural_coordinate_system()); Point<3> compare_point2(coordinate_system->natural_coordinate_system()); @@ -1330,7 +1352,7 @@ namespace WorldBuilder } } std::pair result; - result.first = spreading_velocity; + result.first = spreading_velocity_at_ridge / 31557600; // m/s; result.second = distance_ridge; return result; } diff --git a/tests/gwb-dat/cartesian_variable_spreading.dat b/tests/gwb-dat/cartesian_variable_spreading.dat new file mode 100644 index 000000000..c79f091aa --- /dev/null +++ b/tests/gwb-dat/cartesian_variable_spreading.dat @@ -0,0 +1,14 @@ +# This is a comment in the data +# file. +# Now define parameters: +# dim = 3 +# compositions = 0 +# x y z d T +250e3 -90e3 0 50e3 +250e3 -60e3 0 50e3 +250e3 -30e3 0 50e3 +250e3 -1e3 0 50e3 +250e3 1e3 0 50e3 +250e3 30e3 0 50e3 +250e3 60e3 0 50e3 +250e3 90e3 0 50e3 diff --git a/tests/gwb-dat/cartesian_variable_spreading.wb b/tests/gwb-dat/cartesian_variable_spreading.wb new file mode 100644 index 000000000..8d278f6e2 --- /dev/null +++ b/tests/gwb-dat/cartesian_variable_spreading.wb @@ -0,0 +1,17 @@ +{ + "version": "0.6", + "coordinate system":{"model":"cartesian"}, + "gravity model":{"model":"uniform", "magnitude":10}, + "surface temperature":273, "force surface temperature":true, + "potential mantle temperature":1673, "thermal expansion coefficient":3.1e-5, + "specific heat":1250, "thermal diffusivity":1.0e-6, + "features": + [ + { "model":"oceanic plate", "name":"Subducting", "max depth":250e3,"min depth":0, + "coordinates" :[[500e3, -100e3],[500e3, 100e3],[0, 100e3],[0, -100e3]], + "temperature models":[ + {"model":"half space model", "min depth":0, "max depth":100e3, "spreading velocity":[ [0,[[0.001, 0.0005, 0.0001]]], [1,[[0.000005, 0.0005]]] ], + "top temperature":273, + "ridge coordinates": [[[100e3,-110e3], [150e3, -60e3], [100e3, 0]], [[400e3, 0], [400e3, 110e3]]]}]} + ] +} diff --git a/tests/gwb-dat/cartesian_variable_spreading/screen-output.log b/tests/gwb-dat/cartesian_variable_spreading/screen-output.log new file mode 100644 index 000000000..8d7356469 --- /dev/null +++ b/tests/gwb-dat/cartesian_variable_spreading/screen-output.log @@ -0,0 +1,9 @@ +# x y z d g T tag +250e3 -90e3 0 50e3 751.596 0 +250e3 -60e3 0 50e3 761.356 0 +250e3 -30e3 0 50e3 751.596 0 +250e3 -1e3 0 50e3 728.246 0 +250e3 1e3 0 50e3 329.759 0 +250e3 30e3 0 50e3 489.645 0 +250e3 60e3 0 50e3 574.848 0 +250e3 90e3 0 50e3 639.408 0 diff --git a/tests/gwb-dat/screen-output.log b/tests/gwb-dat/screen-output.log new file mode 100644 index 000000000..4489ffda0 --- /dev/null +++ b/tests/gwb-dat/screen-output.log @@ -0,0 +1,9 @@ +# x y z d g T +250e3 -90e3 0 50e3 751.596 +250e3 -60e3 0 50e3 761.356 +250e3 -30e3 0 50e3 751.596 +250e3 -1e3 0 50e3 728.246 +250e3 1e3 0 50e3 329.759 +250e3 30e3 0 50e3 489.645 +250e3 60e3 0 50e3 574.848 +250e3 90e3 0 50e3 639.408 diff --git a/tests/gwb-dat/spherical_variable_spreading.dat b/tests/gwb-dat/spherical_variable_spreading.dat new file mode 100644 index 000000000..d2eb40124 --- /dev/null +++ b/tests/gwb-dat/spherical_variable_spreading.dat @@ -0,0 +1,15 @@ +# This is a comment in the data +# file. +# Now define parameters: +# dim = 3 +# compositions = 0 +# convert spherical = true +# R long lat d T +6321e3 22.5 16 50e3 +6321e3 22.5 18 50e3 +6321e3 22.5 20 50e3 +6321e3 22.5 22 50e3 +6321e3 22.5 23 50e3 +6321e3 22.5 25 50e3 +6321e3 22.5 27 50e3 +6321e3 22.5 29 50e3 diff --git a/tests/gwb-dat/spherical_variable_spreading.wb b/tests/gwb-dat/spherical_variable_spreading.wb new file mode 100644 index 000000000..04969e073 --- /dev/null +++ b/tests/gwb-dat/spherical_variable_spreading.wb @@ -0,0 +1,17 @@ +{ + "version": "0.6", + "coordinate system":{"model":"spherical", "depth method":"starting point"}, + "gravity model":{"model":"uniform", "magnitude":10}, + "surface temperature":273, "force surface temperature":true, + "potential mantle temperature":1673, "thermal expansion coefficient":3.1e-5, + "specific heat":1250, "thermal diffusivity":1.0e-6, + "features": + [ + { "model":"oceanic plate", "name":"Subducting", "max depth":250e3,"min depth":0, + "coordinates" :[[15, 15],[15, 30],[30, 30],[30, 15]], + "temperature models":[ + {"model":"half space model", "min depth":0, "max depth":100e3, "spreading velocity":[ [0,[[0.005, 0.0001, 0.005]]], [0, [[0.001, 0.00005]]] ], + "top temperature":273, + "ridge coordinates": [[[20,10], [21, 16], [20, 22.5]], [[25, 22.5], [25, 31]]]}]} + ] +} diff --git a/tests/gwb-dat/spherical_variable_spreading/screen-output.log b/tests/gwb-dat/spherical_variable_spreading/screen-output.log new file mode 100644 index 000000000..41fac7bf3 --- /dev/null +++ b/tests/gwb-dat/spherical_variable_spreading/screen-output.log @@ -0,0 +1,9 @@ +# x y z d g T tag +6321e3 22.5 16 50e3 1315.37 0 +6321e3 22.5 18 50e3 1258.52 0 +6321e3 22.5 20 50e3 1208.46 0 +6321e3 22.5 22 50e3 1166.44 0 +6321e3 22.5 23 50e3 696.071 0 +6321e3 22.5 25 50e3 647.504 0 +6321e3 22.5 27 50e3 588.45 0 +6321e3 22.5 29 50e3 510.556 0 diff --git a/tests/unit_tests/unit_test_world_builder.cc b/tests/unit_tests/unit_test_world_builder.cc index 390ded245..a7565be1c 100644 --- a/tests/unit_tests/unit_test_world_builder.cc +++ b/tests/unit_tests/unit_test_world_builder.cc @@ -3820,23 +3820,23 @@ TEST_CASE("WorldBuilder Parameters") CHECK(rapidjson::Pointer("/test/oneOf/0/oneOf/1/default value").Get(prm_temp.declarations)->GetDouble() == Approx(102.)); CHECK(rapidjson::Pointer("/test/oneOf/1/default value").Get(prm_temp.declarations)->GetDouble() == Approx(103.)); } - prm.declare_entry("one value at points one value string",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101.))), + prm.declare_entry("one value at points one value string",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101., 2.))), "Documentation"); - prm.declare_entry("array value at points one value string",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101.))), + prm.declare_entry("array value at points one value string",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101., 2.))), "Documentation"); - prm.declare_entry("value at points set ap val string",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101.))), + prm.declare_entry("value at points set ap val string",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101., 2.))), "Documentation"); - prm.declare_entry("value at points set ap p1 string",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101.))), + prm.declare_entry("value at points set ap p1 string",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101., 2.))), "Documentation"); - prm.declare_entry("value at points set ap p2 string",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101.))), + prm.declare_entry("value at points set ap p2 string",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101., 2.))), "Documentation"); - prm.declare_entry("one value at points one value",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101.))), + prm.declare_entry("one value at points one value",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101., 2.))), "Documentation"); - prm.declare_entry("array value at points one value",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101.))), + prm.declare_entry("array value at points one value",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101., 2.))), "Documentation"); - prm.declare_entry("value at points",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101.))), + prm.declare_entry("value at points",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101., 2.))), "Documentation"); - prm.declare_entry("value at points default ap",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101.))), + prm.declare_entry("value at points default ap",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101., 2.))), "Documentation"); prm.leave_subsection(); const std::vector > additional_points = {Point<2>(-10,-10,cartesian),Point<2>(-10,10,cartesian), From 97ea335a56347f143b5d5711c91f89ade98e2b05 Mon Sep 17 00:00:00 2001 From: danieldouglas92 Date: Fri, 16 Feb 2024 11:02:58 -0600 Subject: [PATCH 2/2] Add changelog --- CHANGELOG.md | 1 + .../subducting_plate_models/temperature/mass_conserving.cc | 2 -- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 49e6f28e3..75ec1e5cc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,7 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm - Added an option to use the plate model as the reference model for the mass conserving temperature of the slab. \[Haoyuan Li; 2024-02-02; [#471](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/471)\] - Added a cookbook for making a transform fault and using this model in ASPECT. \[Juliane Dannberg; 2024-02-14; [#563](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/563)\] - Added a system which allows users to tag features. The tag index can then be written out throught the gwb-dat program. \[Menno Fraters and Timo Heister; 2024-02-15; [[#598](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/598)]\] +- Added variable spreading for mid oceanic ridges. \[Daniel Douglas; 2024-02-16; [#617](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/617)\] ### Changed - Unified the directories `cookbooks/` and `doc/sphinx/user_manual/cookbooks`. All information about cookbooks including the documentation is now bundled in the top-level `cookbooks/` directory. \[Rene Gassmoeller; 2024-02-14; [#558](github.com/GeodynamicWorldBuilder/WorldBuilder/pull/558)\] 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 5667df28a..3d80cb585 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 @@ -254,8 +254,6 @@ namespace WorldBuilder { const double distance_from_plane = distance_from_planes.distance_from_plane; - const int plate_model_summation_number = 100; // for the plate model - if (distance_from_plane <= max_depth && distance_from_plane >= min_depth) {