diff --git a/include/world_builder/features/fault_models/composition/smooth.h b/include/world_builder/features/fault_models/composition/smooth.h new file mode 100644 index 000000000..19fce097a --- /dev/null +++ b/include/world_builder/features/fault_models/composition/smooth.h @@ -0,0 +1,96 @@ +/* + Copyright (C) 2018 - 2020 by the authors of the World Builder code. + + This file is part of the World Builder. + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published + by the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program. If not, see . +*/ + +#ifndef _world_builder_features_fault_composition_smooth_h +#define _world_builder_features_fault_composition_smooth_h + +#include +#include +#include + + +namespace WorldBuilder +{ + namespace Features + { + namespace FaultModels + { + namespace Composition + { + /** + * This class represents a subducting plate and can implement submodules + * for temperature and composition. These submodules determine what + * the returned temperature or composition of the temperature and composition + * functions of this class will be. + */ + class Smooth: public Interface + { + public: + /** + * constructor + */ + Smooth(WorldBuilder::World *world); + + /** + * Destructor + */ + ~Smooth(); + + /** + * declare and read in the world builder file into the parameters class + */ + static + void declare_entries(Parameters &prm, const std::string &parent_name = ""); + + /** + * declare and read in the world builder file into the parameters class + */ + void parse_entries(Parameters &prm) override final; + + + /** + * Returns a composition based on the given position, depth in the model, + * and current compostion. + */ + double get_composition(const Point<3> &position, + const double depth, + const unsigned int composition_number, + double composition, + const double feature_min_depth, + const double feature_max_depth, + const WorldBuilder::Utilities::PointDistanceFromCurvedPlanes &distance_from_planes, + const Utilities::AdditionalParameters &additional_paramters) const override final; + + private: + // linear temperature submodule parameters + double min_distance; + double side_distance; + double center_composition; + // currenly not using the side composition, but maybe usefu if you want another composition towards the end + double side_composition; + unsigned int compositions; + Utilities::Operations operation; + + }; + } + } + } +} + +#endif diff --git a/include/world_builder/features/utilities.h b/include/world_builder/features/utilities.h index 227932614..3cd21cb65 100644 --- a/include/world_builder/features/utilities.h +++ b/include/world_builder/features/utilities.h @@ -33,7 +33,7 @@ namespace WorldBuilder { enum class Operations { - REPLACE,ADD,SUBTRACT + REPLACE,ADD,SUBTRACT,MAX }; /** @@ -64,6 +64,9 @@ namespace WorldBuilder case Utilities::Operations::SUBTRACT: return old_value - new_value; + case Utilities::Operations::MAX: + return std::max(old_value, new_value); + default: WBAssert(false,"Operation not found."); } diff --git a/include/world_builder/utilities.h b/include/world_builder/utilities.h index db3c610bc..72289548e 100644 --- a/include/world_builder/utilities.h +++ b/include/world_builder/utilities.h @@ -412,6 +412,16 @@ namespace WorldBuilder */ std::array,3> euler_angles_to_rotation_matrix(double phi1, double theta, double phi2); + + /** + * Read a file and distribute the content over all MPI processes. + * If WB_WITH_MPI is not defined, this function will just read the file. + * + * @param filename The name of the file to read. + * @return The content of the file. + */ + std::string + read_and_distribute_file_content(const std::string &filename); } // namespace Utilities } // namespace WorldBuilder diff --git a/source/features/fault.cc b/source/features/fault.cc index 7222b34b2..2d4bd12cc 100644 --- a/source/features/fault.cc +++ b/source/features/fault.cc @@ -612,13 +612,34 @@ namespace WorldBuilder section_fraction * (total_fault_length[next_section] - total_fault_length[current_section]); + const double rounded_corner = sqrt(abs(0.25*thickness_local*thickness_local - distance_from_plane*distance_from_plane)); - // Because both sides return positve values, we have to - // devide the thickness_local by two + // We need 3d points in order to conert them from natural to cartesian coordinates + const bool bool_cartesian = (this->world->parameters.coordinate_system->natural_coordinate_system() == cartesian); + + const Point<3> P1 (bool_cartesian ? coordinates[current_section][0] : starting_radius, + bool_cartesian ? coordinates[current_section][1] : coordinates[current_section][0], + bool_cartesian ? starting_radius : coordinates[current_section][1], + this->world->parameters.coordinate_system->natural_coordinate_system()); + + const Point<3> P2 (bool_cartesian ? coordinates[next_section][0] : starting_radius, + bool_cartesian ? coordinates[next_section][1] : coordinates[current_section][0], + bool_cartesian ? starting_radius : coordinates[next_section][1], + this->world->parameters.coordinate_system->natural_coordinate_system()); + + Point<3> P1_cartesian(this->world->parameters.coordinate_system->natural_to_cartesian_coordinates(P1.get_array()),cartesian); + Point<3> P2_cartesian(this->world->parameters.coordinate_system->natural_to_cartesian_coordinates(P2.get_array()),cartesian); + + const double fault_length_at_surface = (P2_cartesian - P1_cartesian).norm(); + + // Because both sides return positive values, we have to + // divide the thickness_local by two if (std::fabs(distance_from_plane) > 0 && std::fabs(distance_from_plane) <= thickness_local * 0.5 && distance_along_plane > 0 && - distance_along_plane <= max_fault_length) + distance_along_plane <= max_fault_length && + section_fraction * fault_length_at_surface > - rounded_corner && + (section_fraction - 1) * fault_length_at_surface < rounded_corner) { // Inside the fault! const Features::Utilities::AdditionalParameters additional_parameters = {max_fault_length,thickness_local}; diff --git a/source/features/fault_models/composition/smooth.cc b/source/features/fault_models/composition/smooth.cc new file mode 100644 index 000000000..2199b083b --- /dev/null +++ b/source/features/fault_models/composition/smooth.cc @@ -0,0 +1,124 @@ +/* + Copyright (C) 2018 - 2020 by the authors of the World Builder code. + + This file is part of the World Builder. + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published + by the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with this program. If not, see . +*/ + +#include +#include +#include +#include + +#include +#include +#include +#include + + +namespace WorldBuilder +{ + using namespace Utilities; + + namespace Features + { + namespace FaultModels + { + namespace Composition + { + Smooth::Smooth(WorldBuilder::World *world_) + : + min_distance(NaN::DSNAN), + side_distance(NaN::DSNAN), + center_composition(NaN::DSNAN), + side_composition(NaN::DSNAN), + operation(Utilities::Operations::MAX) + { + this->world = world_; + this->name = "smooth"; + } + + Smooth::~Smooth() + { } + + void + Smooth::declare_entries(Parameters &prm, const std::string &) + { + // Add compositions to the required parameters. + prm.declare_entry("", Types::Object({"compositions"}), "Compositional model object"); + + prm.declare_entry("min distance fault center", Types::Double(0), + "The distance in meters from which the composition of this feature is present."); + + prm.declare_entry("side distance fault center", Types::Double(std::numeric_limits::max()), + "The distance over which the composition is reduced from 1 to 0."); + + prm.declare_entry("center composition", Types::Double(1.), + "The composition at the center of the fault."); + + prm.declare_entry("side composition", Types::Double(0), + "The composition at the sides of this feature."); + + prm.declare_entry("operation", Types::String("max", std::vector {"max"}), + "Whether the value should replace any value previously defined at this location (replace) or " + "add the value to the previously define value (add, not implemented). Replacing implies that all values not " + "explicitly defined are set to zero."); + } + + void + Smooth::parse_entries(Parameters &prm) + { + min_distance = prm.get("min distance fault center"); + side_distance = prm.get("side distance fault center"); + WBAssert(side_distance >= min_distance, "distance at the side needs to be larger or equal than the min distance."); + operation = Utilities::string_operations_to_enum(prm.get("operation")); + center_composition = prm.get("center composition"); + side_composition = prm.get("side composition"); + compositions = prm.get("compositions"); + } + + + double + Smooth::get_composition( const Point<3> &, + const double , + const unsigned int composition_number, + double composition_, + const double , + const double , + const WorldBuilder::Utilities::PointDistanceFromCurvedPlanes &distance_from_planes, + const Utilities::AdditionalParameters &additional_paramters) const + { + double composition = 0; + + if (compositions == composition_number) + { + + // Hyperbolic tangent goes from 0 to 1 over approximately x=(0, 2) without any arguements. The function is written + // so that the composition returned 1 to 0 over the side_distance on either sides. + composition = (center_composition - std::tanh(10 * (distance_from_planes.distance_from_plane - side_distance/2)/side_distance ) )/2 ; + + return Utilities::apply_operation(operation,composition_,composition); + } + + return composition; + + } + + WB_REGISTER_FEATURE_FAULT_COMPOSITION_MODEL (Smooth, smooth) + } + } + } +} + diff --git a/source/features/utilities.cc b/source/features/utilities.cc index 7f52efe0e..d76f52be5 100644 --- a/source/features/utilities.cc +++ b/source/features/utilities.cc @@ -31,6 +31,7 @@ namespace WorldBuilder { if (operation == "add") return Operations::ADD; if (operation == "subtract") return Operations::SUBTRACT; + if (operation == "max") return Operations::MAX; WBAssert(operation == "replace", "Could not find operation: " << operation << "."); return Operations::REPLACE; diff --git a/source/parameters.cc b/source/parameters.cc index 6525cbe93..a5649cb71 100644 --- a/source/parameters.cc +++ b/source/parameters.cc @@ -80,16 +80,9 @@ namespace WorldBuilder } path_level =0; - // Now read in the world builder file into a file stream and - // put it into a the rapidjason document - std::ifstream json_input_stream(filename.c_str()); - - // Get world builder file and check whether it exists - WBAssertThrow(json_input_stream.good(), - "Could not find the world builder file at the specified location: " + filename); - - WBAssert(json_input_stream, "Could not read the world builder file."); - + // Now read in the world builder file into a stringstream and + // put it into a the rapidjson document + std::stringstream json_input_stream(WorldBuilder::Utilities::read_and_distribute_file_content(filename)); rapidjson::IStreamWrapper isw(json_input_stream); // relaxing sytax by allowing comments () for now, maybe also allow trailing commas and (kParseTrailingCommasFlag) and nan's, inf etc (kParseNanAndInfFlag)? @@ -120,7 +113,6 @@ namespace WorldBuilder )); WBAssertThrow(parameters.IsObject(), "World builder file is is not an object."); - json_input_stream.close(); SchemaDocument schema(declarations); diff --git a/source/utilities.cc b/source/utilities.cc index b5c274e7e..19801123f 100644 --- a/source/utilities.cc +++ b/source/utilities.cc @@ -19,6 +19,15 @@ #include #include +#include +#include +#include +#include + +#ifdef WB_WITH_MPI +#define OMPI_SKIP_MPICXX 1 +#include +#endif #include "world_builder/nan.h" #include "world_builder/utilities.h" @@ -505,11 +514,16 @@ namespace WorldBuilder "angle_at_begin_segment_with_surface are implemented"); double min_distance_check_point_surface_2d_line = INFINITY; + double min_distance_outside_section = INFINITY; size_t i_section_min_distance = 0; Point<2> closest_point_on_line_2d(NaN::DSNAN,NaN::DSNAN,natural_coordinate_system); + size_t i_section_min_distance_outside_section = 0; Point<2> closest_point_on_line_2d_temp(0,0,natural_coordinate_system); + Point<2> closest_point_outside_section(0,0,natural_coordinate_system); double fraction_CPL_P1P2_strict = INFINITY; // or NAN? + double fraction_CPL_P1P2_outside_section = INFINITY; double fraction_CPL_P1P2 = INFINITY; + bool is_within_a_section = false; bool continue_computation = false; if (interpolation_type != InterpolationType::ContinuousMonotoneSpline) @@ -552,11 +566,37 @@ namespace WorldBuilder i_section_min_distance = i_section; closest_point_on_line_2d = closest_point_on_line_2d_temp; fraction_CPL_P1P2_strict = fraction_CPL_P1P2_strict_temp; + is_within_a_section = true; + } + else if (((fraction_CPL_P1P2_strict_temp >= 1. && fraction_CPL_P1P2_strict_temp <= 2.) || + (fraction_CPL_P1P2_strict_temp >= -1. && fraction_CPL_P1P2_strict_temp < 0.) ) + && fabs(min_distance_check_point_surface_2d_line_temp) < fabs(min_distance_outside_section)) + { + min_distance_outside_section = min_distance_check_point_surface_2d_line_temp; + i_section_min_distance_outside_section = i_section; + closest_point_outside_section = closest_point_on_line_2d_temp; + fraction_CPL_P1P2_outside_section = fraction_CPL_P1P2_strict_temp; } } + + if (!is_within_a_section) + { + min_distance_check_point_surface_2d_line = min_distance_outside_section; + i_section_min_distance = i_section_min_distance_outside_section; + closest_point_on_line_2d = closest_point_outside_section; + fraction_CPL_P1P2_strict = fraction_CPL_P1P2_outside_section; + } + + if (!is_within_a_section) + { + min_distance_check_point_surface_2d_line = min_distance_outside_section; + i_section_min_distance = i_section_min_distance_outside_section; + closest_point_on_line_2d = closest_point_outside_section; + fraction_CPL_P1P2_strict = fraction_CPL_P1P2_outside_section; + } // If the point on the line does not lay between point P1 and P2 // then ignore it. Otherwise continue. - continue_computation = (fabs(fraction_CPL_P1P2_strict) < INFINITY && fraction_CPL_P1P2_strict >= 0. && fraction_CPL_P1P2_strict <= 1.); + continue_computation = (fabs(fraction_CPL_P1P2_strict) < INFINITY && fraction_CPL_P1P2_strict >= -1. && fraction_CPL_P1P2_strict <= 2.); fraction_CPL_P1P2 = global_x_list[i_section_min_distance] - static_cast(global_x_list[i_section_min_distance]) + (global_x_list[i_section_min_distance+1]-global_x_list[i_section_min_distance]) * fraction_CPL_P1P2_strict; @@ -615,12 +655,23 @@ namespace WorldBuilder } } double solution = min_estimate_solution; - - continue_computation = (solution > 0 && floor(solution) <= global_x_list[point_list.size()-2] && floor(solution) >= 0); - - closest_point_on_line_2d = Point<2>(x_spline(solution),y_spline(solution),natural_coordinate_system); - i_section_min_distance = static_cast(floor(solution)); - fraction_CPL_P1P2 = solution-floor(solution); + bool is_outside_section = (floor(solution) > global_x_list[point_list.size()-2]); + + continue_computation = ((solution > 0 && floor(solution) <= global_x_list[point_list.size()-2] && floor(solution) >= 0) || is_outside_section); + + if (!is_outside_section) + { + closest_point_on_line_2d = Point<2>(x_spline(solution),y_spline(solution),natural_coordinate_system); + i_section_min_distance = static_cast(floor(solution)); + fraction_CPL_P1P2 = solution-floor(solution); + } + + else if (is_outside_section) + { + closest_point_on_line_2d = Point<2>(x_spline(solution),y_spline(solution),natural_coordinate_system); + i_section_min_distance = closest_point_on_line_2d.cheap_relative_distance(check_point_surface_2d); + fraction_CPL_P1P2 = solution-1; + } } // We now need 3d points from this point on, so make them. // The order of a Cartesian coordinate is x,y,z and the order of @@ -1376,6 +1427,88 @@ namespace WorldBuilder return rot_matrix; } + std::string + read_and_distribute_file_content(const std::string &filename) + { + std::string data_string; + + const unsigned int invalid_unsigned_int = static_cast(-1); + + const MPI_Comm comm = MPI_COMM_WORLD; + int my_rank = invalid_unsigned_int; + MPI_Comm_rank(comm, &my_rank); + if (my_rank == 0) + { + std::size_t filesize; + std::ifstream filestream; + filestream.open(filename.c_str()); + + // Need to convert to unsigned int, because MPI_Bcast does not support + // size_t or const unsigned int + unsigned int invalid_file_size = invalid_unsigned_int; + + if (!filestream) + { + // broadcast failure state, then throw + const int ierr = MPI_Bcast(&invalid_file_size, 1, MPI_UNSIGNED, 0, comm); + WBAssertThrow (ierr, + std::string("Could not open file <") + filename + ">."); + } + + // Read data from disk + std::stringstream datastream; + + try + { + datastream << filestream.rdbuf(); + } + catch (const std::ios::failure &) + { + // broadcast failure state, then throw + const int ierr = MPI_Bcast(&invalid_file_size, 1, MPI_UNSIGNED, 0, comm); + WBAssertThrow(ierr == 0, "MPI_Bcast failed."); + WBAssertThrow (false, + std::string("Could not read file content from <") + filename + ">."); + } + + data_string = datastream.str(); + filesize = data_string.size(); + + // Distribute data_size and data across processes + int ierr = MPI_Bcast(&filesize, 1, MPI_UNSIGNED, 0, comm); + WBAssertThrow(ierr == 0, "MPI_Bcast failed."); + + // Receive and store data + ierr = MPI_Bcast(&data_string[0], + filesize, + MPI_CHAR, + 0, + comm); + WBAssertThrow(ierr == 0, "MPI_Bcast failed."); + } + else + { + // Prepare for receiving data + unsigned int filesize = 0; + int ierr = MPI_Bcast(&filesize, 1, MPI_UNSIGNED, 0, comm); + WBAssertThrow(ierr == 0, "MPI_Bcast failed."); + WBAssertThrow(filesize != invalid_unsigned_int, + std::string("Could not open file <") + filename + ">."); + + data_string.resize(filesize); + + // Receive and store data + ierr = MPI_Bcast(&data_string[0], + filesize, + MPI_CHAR, + 0, + comm); + WBAssertThrow(ierr == 0, "MPI_Bcast failed."); + } + + return data_string; + } + template std::array convert_point_to_array<2>(const Point<2> &point_); template std::array convert_point_to_array<3>(const Point<3> &point_); } // namespace Utilities