Skip to content

Commit

Permalink
Merge 1148dd9 into 7b227fc
Browse files Browse the repository at this point in the history
  • Loading branch information
@danieldouglas92
danieldouglas92 committed Feb 23, 2024
2 parents 7b227fc + 1148dd9 commit 93f5610
Show file tree
Hide file tree
Showing 8 changed files with 298 additions and 78 deletions.
5 changes: 3 additions & 2 deletions include/world_builder/features/subducting_plate_models/temperature/mass_conserving.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,8 +110,9 @@ namespace WorldBuilder
double min_depth;
double max_depth;
double density;
std::pair<std::vector<double>,std::vector<double>> plate_velocities;
std::vector<std::vector<double>> plate_velocities_at_each_ridge_point;
std::vector<std::vector<double>> subducting_velocities;
std::pair<std::vector<double>,std::vector<double>> ridge_spreading_velocities;
std::vector<std::vector<double>> ridge_spreading_velocities_at_each_ridge_point;
double mantle_coupling_depth;
double forearc_cooling_factor;
double thermal_conductivity;
Expand Down
4 changes: 3 additions & 1 deletion include/world_builder/parameters.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,8 @@ namespace WorldBuilder
template<class T>
std::vector<T> get_vector(const std::string &name);

std::vector<std::vector<double>> get_vector_or_double(const std::string &name);

/**
* A specialized version of get which can return a value at points type.
* \param name The name of the entry to retrieved
Expand All @@ -131,7 +133,7 @@ namespace WorldBuilder
std::pair<std::vector<double>,std::vector<double>> get_value_at_array(const std::string &name);

/**
* A specialized verions of get which can return vectors/arrays.
* A specialized version of get which can return vectors/arrays.
* This version is designed for the plugin system.
* \param name The name of the entry to retrieved
*/
Expand Down
17 changes: 15 additions & 2 deletions include/world_builder/utilities.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -460,11 +460,24 @@ namespace WorldBuilder
* @param position_in_natural_coordinates_at_min_depth the current position in natural_coordinates
* @return The content of the file.
*/
std::pair<double, double>
std::vector<double>
calculate_ridge_distance_and_spreading(std::vector<std::vector<Point<2>>> mid_oceanic_ridges,
std::vector<std::vector<double>> mid_oceanic_spreading_velocities,
const std::unique_ptr<WorldBuilder::CoordinateSystems::Interface> &coordinate_system,
const Objects::NaturalCoordinate &position_in_natural_coordinates_at_min_depth);
const Objects::NaturalCoordinate &position_in_natural_coordinates_at_min_depth,
std::vector<std::vector<double>> subducting_plate_velocities = {{0.0}},
std::vector<double> ridge_migration_times = {0.0});

// todo_effective
/**
* Calculate the effective plate ages of a point on the slab surface, and also calculates
* the effective trench ages at the start of subduction.
* @param ridge_parameters The distance and spreading velocity relative to a mid ocean ridge
* @param distance_along_plane The distance along the slab surface plane
* @return The effective plate age and the trench age
*/
std::vector<double>
calculate_effective_trench_and_plate_ages(std::vector<double> ridge_parameters, double distance_along_plane);

} // namespace Utilities
} // namespace WorldBuilder
Expand Down
24 changes: 14 additions & 10 deletions source/world_builder/features/oceanic_plate_models/temperature/half_space_model.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.01),Types::Array(Types::ValueAtPoints(0.01, std::numeric_limits<uint64_t>::max()))),
prm.declare_entry("spreading velocity", Types::OneOf(Types::Double(0.05),Types::Array(Types::ValueAtPoints(0.05, std::numeric_limits<size_t>::max()))),
"The spreading velocity of the plate in meter per year. "
"This is the velocity with which one side moves away from the ridge.");

Expand Down Expand Up @@ -123,9 +123,13 @@ namespace WorldBuilder
for (unsigned int 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]);
{
spreading_rates_for_ridge.push_back(spreading_velocities.second[0]);
}
else
spreading_rates_for_ridge.push_back(spreading_velocities.second[ridge_point_index]);
{
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);
Expand Down Expand Up @@ -160,15 +164,15 @@ namespace WorldBuilder
this->world->specific_heat) * depth);
}

std::pair<double, double> 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);
std::vector<double> 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);



const double thermal_diffusivity = this->world->thermal_diffusivity;
const double age = ridge_parameters.second / ridge_parameters.first;
const double age = ridge_parameters[1] / ridge_parameters[0];

double temperature = bottom_temperature_local;

Expand All @@ -178,12 +182,12 @@ namespace WorldBuilder
<< ". Relevant variables: bottom_temperature_local = " << bottom_temperature_local
<< ", top_temperature = " << top_temperature
<< ", max_depth = " << max_depth
<< ", spreading_velocity = " << ridge_parameters.first
<< ", spreading_velocity = " << ridge_parameters[0]
<< ", 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 = " << ridge_parameters.first
<< ", spreading_velocity = " << ridge_parameters[0]
<< ", thermal_diffusivity = " << thermal_diffusivity
<< ", age = " << age << '.');

Expand Down
30 changes: 17 additions & 13 deletions source/world_builder/features/oceanic_plate_models/temperature/plate_model.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.01),Types::Array(Types::ValueAtPoints(0.01, std::numeric_limits<uint64_t>::max()))),
prm.declare_entry("spreading velocity", Types::OneOf(Types::Double(0.05),Types::Array(Types::ValueAtPoints(0.05, std::numeric_limits<size_t>::max()))),
"The spreading velocity of the plate in meter per year. "
"This is the velocity with which one side moves away from the ridge.");

Expand Down Expand Up @@ -118,9 +118,13 @@ namespace WorldBuilder
for (unsigned int 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]);
{
spreading_rates_for_ridge.push_back(spreading_velocities.second[0]);
}
else
spreading_rates_for_ridge.push_back(spreading_velocities.second[ridge_point_index]);
{
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);
Expand Down Expand Up @@ -158,13 +162,13 @@ namespace WorldBuilder

const int summation_number = 100;

std::pair<double, double> 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);
std::vector<double> 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);

const double thermal_diffusivity = this->world->thermal_diffusivity;
const double age = ridge_parameters.second / ridge_parameters.first;
const double age = ridge_parameters[1] / ridge_parameters[0];
double temperature = top_temperature + (bottom_temperature_local - top_temperature) * (depth / max_depth);

// This formula addresses the horizontal heat transfer by having the spreading velocity and distance to the ridge in it.
Expand All @@ -173,23 +177,23 @@ 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((((ridge_parameters.first * max_depth)/(2 * thermal_diffusivity)) -
std::sqrt(((ridge_parameters.first*ridge_parameters.first*max_depth*max_depth) /
std::exp((((ridge_parameters[0] * max_depth)/(2 * thermal_diffusivity)) -
std::sqrt(((ridge_parameters[0]*ridge_parameters[0]*max_depth*max_depth) /
(4*thermal_diffusivity*thermal_diffusivity)) + double(i) * double(i) * Consts::PI * Consts::PI)) *
((ridge_parameters.first * age) / max_depth)));
((ridge_parameters[0] * age) / max_depth)));

}

WBAssert(!std::isnan(temperature), "Temperature inside plate model is not a number: " << temperature
<< ". Relevant variables: bottom_temperature_local = " << bottom_temperature_local
<< ", top_temperature = " << top_temperature
<< ", max_depth = " << max_depth
<< ", spreading_velocity = " << ridge_parameters.first
<< ", spreading_velocity = " << ridge_parameters[0]
<< ", thermal_diffusivity = " << thermal_diffusivity
<< ", age = " << age << '.');
WBAssert(std::isfinite(temperature), "Temperature inside plate model is not a finite: " << temperature << ". Relevant variables: bottom_temperature_local = " << bottom_temperature_local
<< ", top_temperature = " << top_temperature
<< ", spreading_velocity = " << ridge_parameters.first
<< ", spreading_velocity = " << ridge_parameters[0]
<< ", thermal_diffusivity = " << thermal_diffusivity
<< ", age = " << age << '.');

Expand Down
Loading

0 comments on commit 93f5610

Please sign in to comment.