diff --git a/CHANGELOG.md b/CHANGELOG.md index d0087d950..5fdeaa792 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,7 +11,7 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm ## [Unreleased] ### Added -- Implemented a method for modifying the mass conserving temperature model based on the movement of a spreading center through time. \[Daniel Douglas; Haoyuan Li; 2024-02-29; [#654](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/654)\] +- Implemented the framework that will allow the mass conserving temperature model to account for the the movement of a spreading center through time. \[Daniel Douglas; Haoyuan Li; 2024-02-29; [#654](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/654)\] - Added an option to apply a spline in the mass conserving temperature of the slab. \[Haoyuan Li; 2024-02-27; [#659](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/659)\] - Added a cookbook tutorial for a simple subduction model in 2D Chunk geometry. \[Magali Billen; 2024-02-16; [#535](https://github.com/GeodynamicWorldBuilder/WorldBuilder/issues/535)\] - Added a cookbook tutorial for a simple subduction model in 2D Cartesian geometry. \[Magali Billen; 2024-02-14; [#535](https://github.com/GeodynamicWorldBuilder/WorldBuilder/issues/535)\] @@ -109,7 +109,7 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm - Changed the underlying function which determines when you are inside a slab or a fault. This may slightly change the result, but the new result is more realistic I think. \[Menno Fraters; 2021-05-01; [#229](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/229)\] - The Vizualizer now uses compressed output by default. This decreases the file size and increases performance. \[Menno Fraters; 2021-05-22; [#239](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/239)\] - The Vizualizer buffers output before it writes it to a file. This increases performance. \[Menno Fraters; 2021-05-22; [#239](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/239)\] -- In the fault temperature model linear, the entiry top temperature is changed to center temperature and the entry bottom temperature is changed to side temperature, since this represents the actual sides more accurately. \[Menno Fraters; 2021-07-09; [#260](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/260)\] +- In the fault temperature model linear, the entry top temperature is changed to center temperature and the entry bottom temperature is changed to side temperature, since this represents the actual sides more accurately. \[Menno Fraters; 2021-07-09; [#260](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/260)\] - A large overhoal of the distance_point_from_curved_planes function improving the accuracy in spherical coordinates. This may slightly change the results, but it should be an improvement both in accuracy and performance. Also makes available a coordinate system independent distrance function and some fast trigonometric functions. \[Menno Fraters; 2021-06-11; [#255](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/255)\] ### Fixed diff --git a/include/world_builder/utilities.h b/include/world_builder/utilities.h index fb112aa04..62f31b2c0 100644 --- a/include/world_builder/utilities.h +++ b/include/world_builder/utilities.h @@ -455,9 +455,17 @@ namespace WorldBuilder * TODO: make the spreading velocity spatially/temporally variable * * @param mid_oceanic_ridges The coordinates of the mid oceanic ridges - * @param spreading_velocity The spreading rate of the mid oceanic ridges + * @param mid_oceanic_spreading_velocities The spreading rate of the mid oceanic ridges at each ridge coordinate * @param coordinate_system The coordinate system * @param position_in_natural_coordinates_at_min_depth the current position in natural_coordinates + * @param subducting_plate_velocities the subducting plate velocities, currently this is only an optional parameter + * that can be used in the mass conserving subducting plate temperature model. This parameter allows the user to + * track the effect of ridge migration on the slab thermal structure + * @param ridge_migration_times the times that the corresponding section of the ridge has been moving, in years. This + * is used in combination with subducting_plate_velocities, and mid_oceanic_spreading_velocities to compute the distance + * that the spreading center has migrated. This vector is obtained from the input parameter "spreading velocity" in the + * mass conserving model when "spreading velocity" has the form: [ [t1,[[v11, v12, ...]], [t2,[[v21, v22, ...]], ... ]. + * where tn is the time that ridge section n has been moving. * @return The content of the file. */ std::vector 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 5fb70f615..ad064978d 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 @@ -258,6 +258,24 @@ namespace WorldBuilder apply_spline = prm.get("apply spline"); spline_n_points = prm.get("number of points in spline"); + + if (subducting_velocities[0].size() > 1) + { + for (unsigned int ridge_index = 0; ridge_index < mid_oceanic_ridges.size(); ridge_index++) + { + WBAssertThrow(subducting_velocities.size() == mid_oceanic_ridges.size() && + subducting_velocities[ridge_index].size() == mid_oceanic_ridges[ridge_index].size(), + "subducting velocity must have the same dimension as the ridge coordinates parameter."); + for (unsigned int point_index = 0; point_index < mid_oceanic_ridges[ridge_index].size(); point_index++) + { + WBAssertThrow(subducting_velocities[ridge_index][point_index] == + ridge_spreading_velocities_at_each_ridge_point[ridge_index][point_index], + "Currently, subducting velocity must equal the spreading velocity to satisfy conservation of mass. " + "This will be changed in the future to allow for the spreading center to move depending on the relative " + "difference between the spreading velocity and the subducting velocity, thereby conserving mass."); + } + } + } } double diff --git a/source/world_builder/utilities.cc b/source/world_builder/utilities.cc index 3288e50c9..01d4a6092 100644 --- a/source/world_builder/utilities.cc +++ b/source/world_builder/utilities.cc @@ -1358,12 +1358,12 @@ namespace WorldBuilder // 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."); + WBAssert(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"); + WBAssert(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]; @@ -1452,6 +1452,8 @@ namespace WorldBuilder double spreading_velocity_at_ridge_pt = spreading_velocity_at_ridge_pt1; double subducting_velocity_at_trench_pt = subducting_velocity_at_trench_pt1; + // This is required in spherical coordinates to ensure that the distance + // returned is the shortest distance around the sphere. if (compare_distance2 < compare_distance1) { compare_distance = compare_distance2; @@ -1476,7 +1478,7 @@ namespace WorldBuilder return result; } - // todo_effective + // TODO: implement method for modifying the age of the slab based on ridge/trench migration. std::vector calculate_effective_trench_and_plate_ages(std::vector ridge_parameters, double distance_along_plane) { @@ -1484,33 +1486,15 @@ namespace WorldBuilder const double seconds_in_year = 60.0 * 60.0 * 24.0 * 365.25; // sec/y const double spreading_velocity = ridge_parameters[0] * seconds_in_year; // m/yr double subducting_velocity = ridge_parameters[2] * seconds_in_year; // m/yr - const double ridge_velocity = spreading_velocity - subducting_velocity; // m/yr - const int sign_v = (ridge_velocity > 0) ? 1 : ((ridge_velocity < 0) ? -1 : 0); - const double time_of_ridge_migration = ridge_parameters[3]; //yr - const double ridge_drift_distance = abs(ridge_velocity) * time_of_ridge_migration; // m - - double effective_age_shift = 0; - double trench_age_shift = 0; if (subducting_velocity <= 0) subducting_velocity = spreading_velocity; - else - { - if (distance_along_plane < ridge_drift_distance) - { - effective_age_shift = (1 - (distance_along_plane / ridge_drift_distance)) * sign_v * \ - distance_along_plane / subducting_velocity; // yr - trench_age_shift = (1 - (distance_along_plane / ridge_drift_distance)) * sign_v * \ - distance_along_plane / spreading_velocity; // yr - } - } - - const double age_at_trench = ridge_parameters[1] / spreading_velocity + trench_age_shift * 0; // m/(m/y) = yr + const double age_at_trench = ridge_parameters[1] / spreading_velocity; // m/(m/y) = yr const double plate_age_sec = age_at_trench * seconds_in_year; // y --> seconds // Plate age increases with distance along the slab in the mantle - double effective_plate_age = plate_age_sec + (distance_along_plane / subducting_velocity + effective_age_shift) * seconds_in_year; // m/(m/y) = y(seconds_in_year) + double effective_plate_age = plate_age_sec + (distance_along_plane / subducting_velocity) * seconds_in_year; // m/(m/y) = y(seconds_in_year) WBAssertThrow(effective_plate_age >= 0, "The age of the subducting plate is less than or equal to 0. " "Effective plate age: " << effective_plate_age); std::vector result; @@ -1519,7 +1503,6 @@ namespace WorldBuilder return result; } - } // namespace Utilities } // namespace WorldBuilder diff --git a/tests/data/type_data.json b/tests/data/type_data.json index 925754d35..41415018e 100644 --- a/tests/data/type_data.json +++ b/tests/data/type_data.json @@ -13,6 +13,8 @@ "value at points set ap val string": [[100,[[1,2],[3,4]]],[105],["test3",[[5,6]]],[300,[[-10,10]]]], "value at points set ap p1 string": [[100,[[1,2],[3,4]]],[105],[200,[[5,6]]],[300,[["test4",10]]]], "value at points set ap p2 string": [[100,[[1,2],[3,4]]],[105],[200,[[5,6]]],[300,[[-10,"test5"]]]], +"value at array full": [[10,[[1, 2, 3]]], [20,[[4, 5]]]], +"vector for vector or double": [[10, 20, 30, 40], [100, 200]], "unsigned int array" :[25,26,27], "bool array" :[true,false,true,false,true,false], "bool array nob" :[true,false,true,false,"True","false"], diff --git a/tests/gwb-dat/cartesian_no_ridge_drift.dat b/tests/gwb-dat/cartesian_no_ridge_drift.dat new file mode 100644 index 000000000..1804f3775 --- /dev/null +++ b/tests/gwb-dat/cartesian_no_ridge_drift.dat @@ -0,0 +1,20 @@ +# This is a comment in the data +# file. +# Now define parameters: +# dim = 2 +# compositions = 0 +# x y d T C0 +500e3 380e3 20e3 +525e3 380e3 20e3 +550e3 380e3 20e3 +575e3 380e3 20e3 +600e3 380e3 20e3 +625e3 380e3 20e3 +650e3 380e3 20e3 +675e3 380e3 20e3 +700e3 380e3 20e3 +600e3 360e3 40e3 +650e3 340e3 60e3 +700e3 320e3 80e3 +750e3 280e3 120e3 +800e3 240e3 160e3 \ No newline at end of file diff --git a/tests/gwb-dat/cartesian_no_ridge_drift.wb b/tests/gwb-dat/cartesian_no_ridge_drift.wb new file mode 100644 index 000000000..8bc20540c --- /dev/null +++ b/tests/gwb-dat/cartesian_no_ridge_drift.wb @@ -0,0 +1,87 @@ +{ + "version":"0.6", + "coordinate system":{"model":"cartesian"}, + "gravity model":{"model":"uniform", "magnitude":10}, + "cross section":[[0,0],[10000e3,0]],"surface temperature":273, "force surface temperature":true, + "potential mantle temperature":1673, "thermal expansion coefficient":3.1e-5, + "specific heat":1000, "thermal diffusivity":1.0e-6, + "features": + [ + { + "model": "oceanic plate", + "name": "Subducting Plate plate", + "max depth": 100e3, + "min depth": 0.0, + "coordinates": [ [0, 1000], [0, -1000], [500e3, -1000], [500e3, 1000] ], + "temperature models": [ + { + "model": "half space model", + "min depth": 0.0, + "max depth": 100e3, + "spreading velocity": 0.05, + "ridge coordinates": [[[-100e3, 1000.0], [-100e3, -1000.0]]] + } + ], + "composition models": + [ + {"model": "uniform", "min depth": 0.0, "max depth": 100e3, "compositions": [0]} + ] + }, + + { + "model": "oceanic plate", + "name": "Overriding Plate", + "max depth": 100e3, + "min depth": 0.0, + "coordinates": [ [1000e3, 1000], [1000e3, -1000], [500e3, -1000], [500e3, 1000] ], + "temperature models": + [ + { + "model": "half space model", "min depth": 0.0, "max depth": 100e3, "spreading velocity": 0.03, + "ridge coordinates": [[[-100e3, 1000.0], [-100e3, -1000.0]]] + } + ] + }, + + { + "model": "subducting plate", "name": "initial slab", "coordinates": [[500e3, -1000.0], [500e3, 1000.0]], + "dip point": [400e6, 0], + "segments": [ + {"length": 418880.0, + "thickness": [100e3], + "top truncation": [-100000.0], + "angle": [0,60], + "composition models": [ + { + "model": "uniform", + "compositions": [0], + "max distance slab top": 100e3 + } + ] + }, + { + "length": 100000.0, + "thickness": [100e3], + "top truncation": [-100000.0], + "angle": [60, 60] + } + ], + "temperature models": [ + { + "model": "mass conserving", + "density": 3300, + "thermal conductivity": 3.3, + "adiabatic heating": true, + "plate velocity": 0.05, + "subducting velocity": [[0.05, 0.05]], + "ridge coordinates": [[[-100e3, 1000.0], [-100e3, -1000.0]]], + "coupling depth": 50000.0, + "taper distance": 100000.0, + "min distance slab top": -100000.0, + "max distance slab top": 100e3, + "reference model name": "half space model" + } + ] + } + ] +} diff --git a/tests/gwb-dat/cartesian_no_ridge_drift/screen-output.log b/tests/gwb-dat/cartesian_no_ridge_drift/screen-output.log new file mode 100644 index 000000000..d4dc961a0 --- /dev/null +++ b/tests/gwb-dat/cartesian_no_ridge_drift/screen-output.log @@ -0,0 +1,15 @@ +# x z d T tag +500e3 380e3 20e3 1005.87 1 +525e3 380e3 20e3 969.259 1 +550e3 380e3 20e3 882.454 1 +575e3 380e3 20e3 744.612 1 +600e3 380e3 20e3 558.663 1 +625e3 380e3 20e3 635.761 1 +650e3 380e3 20e3 809.857 1 +675e3 380e3 20e3 814.935 1 +700e3 380e3 20e3 804.603 1 +600e3 360e3 40e3 1138.89 1 +650e3 340e3 60e3 1169.63 1 +700e3 320e3 80e3 1035.12 1 +750e3 280e3 120e3 1091.76 1 +800e3 240e3 160e3 926.557 1 diff --git a/tests/gwb-dat/spherical_ridge_distance_check.dat b/tests/gwb-dat/spherical_ridge_distance_check.dat new file mode 100644 index 000000000..d11d51a5f --- /dev/null +++ b/tests/gwb-dat/spherical_ridge_distance_check.dat @@ -0,0 +1,9 @@ +# 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 20 50e3 +6321e3 22.5 22 50e3 diff --git a/tests/gwb-dat/spherical_ridge_distance_check.wb b/tests/gwb-dat/spherical_ridge_distance_check.wb new file mode 100644 index 000000000..b2a75f137 --- /dev/null +++ b/tests/gwb-dat/spherical_ridge_distance_check.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" :[[-380, 15],[-380, 30],[380, 30],[380, 15]], + "temperature models":[ + {"model":"half space model", "min depth":0, "max depth":100e3, "spreading velocity":0.05, + "top temperature":273, + "ridge coordinates": [[[-370,10], [-371, 16], [-373, 22.5]], [[371, 22.5], [375, 31]]]}]} + ] +} diff --git a/tests/gwb-dat/spherical_ridge_distance_check/screen-output.log b/tests/gwb-dat/spherical_ridge_distance_check/screen-output.log new file mode 100644 index 000000000..c01acb7ca --- /dev/null +++ b/tests/gwb-dat/spherical_ridge_distance_check/screen-output.log @@ -0,0 +1,3 @@ +# x y z d g T tag +6321e3 22.5 20 50e3 1045.9 0 +6321e3 22.5 22 50e3 1044.78 0 diff --git a/tests/unit_tests/approval_tests/approved/unit_test_world_builder.WorldBuilder_Parameters.txt b/tests/unit_tests/approval_tests/approved/unit_test_world_builder.WorldBuilder_Parameters.txt index 94bb65f65..dbae6f707 100644 --- a/tests/unit_tests/approval_tests/approved/unit_test_world_builder.WorldBuilder_Parameters.txt +++ b/tests/unit_tests/approval_tests/approved/unit_test_world_builder.WorldBuilder_Parameters.txt @@ -108,4 +108,33 @@ TITLE [104] = 7 [105] = 8 [106] = 9 +[107] = 2 +[108] = 10 +[109] = 20 +[110] = 5 +[111] = 1 +[112] = 2 +[113] = 3 +[114] = 4 +[115] = 5 +[116] = 1 +[117] = 0 +[118] = 1 +[119] = 200 +[120] = 1 +[121] = 0 +[122] = 1 +[123] = 101 +[124] = 2 +[125] = 4 +[126] = 10 +[127] = 20 +[128] = 30 +[129] = 40 +[130] = 2 +[131] = 100 +[132] = 200 +[133] = 1 +[134] = 1 +[135] = 200 diff --git a/tests/unit_tests/unit_test_world_builder.cc b/tests/unit_tests/unit_test_world_builder.cc index 5d879736d..9db97b789 100644 --- a/tests/unit_tests/unit_test_world_builder.cc +++ b/tests/unit_tests/unit_test_world_builder.cc @@ -3849,6 +3849,12 @@ TEST_CASE("WorldBuilder Parameters") "Documentation"); prm.declare_entry("value at points default ap",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101., 2.))), "Documentation"); + prm.declare_entry("value at array full",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101., std::numeric_limits::max()))), + "Documentation"); + prm.declare_entry("value at array",Types::OneOf(Types::Double(101.),Types::Array(Types::ValueAtPoints(101., std::numeric_limits::max()))), + "Documentation"); + prm.declare_entry("vector for vector or double", Types::OneOf(Types::Double(-1), Types::Array(Types::Array(Types::Double(-1), 1), 1)), + "Documentation"); prm.leave_subsection(); const std::vector > additional_points = {Point<2>(-10,-10,cartesian),Point<2>(-10,10,cartesian), Point<2>(10,10,cartesian),Point<2>(10,-10,cartesian) @@ -4092,7 +4098,47 @@ TEST_CASE("WorldBuilder Parameters") Contains("Could not convert values of /vector of vectors of points<2> nan/1 into doubles")); - + std::pair,std::vector> value_at_array = prm.get_value_at_array("value at array full"); + approval_tests.emplace_back((double)value_at_array.first.size()); + approval_tests.emplace_back(value_at_array.first[0]); + approval_tests.emplace_back(value_at_array.first[1]); + + approval_tests.emplace_back((double)value_at_array.second.size()); + approval_tests.emplace_back(value_at_array.second[0]); + approval_tests.emplace_back(value_at_array.second[1]); + approval_tests.emplace_back(value_at_array.second[2]); + approval_tests.emplace_back(value_at_array.second[3]); + approval_tests.emplace_back(value_at_array.second[4]); + + + std::pair,std::vector> double_value_at_array = prm.get_value_at_array("one value at points one value"); + approval_tests.emplace_back((double)double_value_at_array.first.size()); + approval_tests.emplace_back(double_value_at_array.first[0]); + approval_tests.emplace_back((double)double_value_at_array.second.size()); + approval_tests.emplace_back(double_value_at_array.second[0]); + + std::pair,std::vector> default_value_at_array = prm.get_value_at_array("value at array"); + approval_tests.emplace_back((double)default_value_at_array.first.size()); + approval_tests.emplace_back(default_value_at_array.first[0]); + approval_tests.emplace_back((double)default_value_at_array.second.size()); + approval_tests.emplace_back(default_value_at_array.second[0]); + + std::vector> vector_for_vector_or_double = prm.get_vector_or_double("vector for vector or double"); + approval_tests.emplace_back((double)vector_for_vector_or_double.size()); + approval_tests.emplace_back((double)vector_for_vector_or_double[0].size()); + approval_tests.emplace_back((double)vector_for_vector_or_double[0][0]); + approval_tests.emplace_back((double)vector_for_vector_or_double[0][1]); + approval_tests.emplace_back((double)vector_for_vector_or_double[0][2]); + approval_tests.emplace_back((double)vector_for_vector_or_double[0][3]); + + approval_tests.emplace_back((double)vector_for_vector_or_double[1].size()); + approval_tests.emplace_back((double)vector_for_vector_or_double[1][0]); + approval_tests.emplace_back((double)vector_for_vector_or_double[1][1]); + + std::vector> double_for_vector_or_double = prm.get_vector_or_double("one value at points one value"); + approval_tests.emplace_back((double)double_for_vector_or_double.size()); + approval_tests.emplace_back((double)double_for_vector_or_double[0].size()); + approval_tests.emplace_back((double)double_for_vector_or_double[0][0]); /*CHECK_THROWS_WITH(prm.get_vector("non existent string vector"), Contains("internal error: could not retrieve the default value at"));