diff --git a/.github/workflows/pre_commit.yml b/.github/workflows/pre_commit.yml index 2895110..c8284d3 100644 --- a/.github/workflows/pre_commit.yml +++ b/.github/workflows/pre_commit.yml @@ -9,7 +9,9 @@ jobs: steps: - uses: actions/checkout@v3 - name: Install Conda environment - uses: mamba-org/provision-with-micromamba@main + uses: mamba-org/setup-micromamba@main + with: + environment-file: environment.yml - name: Run precommit shell: bash -l {0} run: | diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index f533170..2df43f8 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -18,7 +18,7 @@ repos: name: cpplint entry: cpplint language: system - args: ["--exclude=thirdparty/", "--filter=-whitespace/comments,-runtime/references,-whitespace/indent,-whitespace/parens,-whitespace/braces,-whitespace/line_length,-whitespace/newline,-build/include_order,-readability/todo,-build/namespaces"] + args: ["--exclude=thirdparty/", "--filter=-whitespace/comments,-runtime/references,-whitespace/indent,-whitespace/parens,-whitespace/braces,-whitespace/line_length,-whitespace/newline,-build/include_order,-readability/todo,-build/namespaces,-build/c++11"] - id: clang-format name: clang-format entry: clang-format diff --git a/flowy/include/models/mr_lava_loba.hpp b/flowy/include/models/mr_lava_loba.hpp index ee01f2e..d21d022 100644 --- a/flowy/include/models/mr_lava_loba.hpp +++ b/flowy/include/models/mr_lava_loba.hpp @@ -5,7 +5,6 @@ #include "flowy/include/config.hpp" #include "flowy/include/definitions.hpp" #include "flowy/include/lobe.hpp" -#include "flowy/include/simulation.hpp" #include "flowy/include/vent_flags.hpp" #include "pdf_cpplib/include/probability_dist.hpp" #include @@ -63,214 +62,204 @@ class CommonLobeDimensions double exp_lobe_exponent = 1; }; -class MrLavaLoba +namespace MrLavaLoba { -public: - CommonLobeDimensions lobe_dimensions; - Config::InputParams input; - - MrLavaLoba( const Config::InputParams & input, std::mt19937 & gen ) - : lobe_dimensions( input ), input( input ), gen( gen ) +// Calculate n_lobes +inline int compute_n_lobes( int idx_flow, const Config::InputParams & input, std::mt19937 & gen ) +{ + int n_lobes{}; + // Number of lobes in the flow is a random number between the min and max values + if( input.a_beta == 0 && input.b_beta == 0 ) { + std::uniform_int_distribution<> dist_num_lobes( input.min_n_lobes, input.max_n_lobes ); + n_lobes = dist_num_lobes( gen ); } - - explicit MrLavaLoba( Simulation & simulation ) - : lobe_dimensions( simulation.input ), input( simulation.input ), gen( simulation.gen ) + // Deterministic number of lobes, such that a beta probability density distribution is used (not a beta + // distribution). However this means that n_lobes could potentially be greater than max_n_lobes + else { + const double x_beta = ( 1.0 * idx_flow ) / ( input.n_flows - 1.0 ); + const double random_number = Math::beta_pdf( x_beta, input.a_beta, input.b_beta ); + n_lobes = static_cast( + std::round( input.min_n_lobes + 0.5 * ( input.max_n_lobes - input.min_n_lobes ) * random_number ) ); } + return n_lobes; +} - // Calculate n_lobes - int compute_n_lobes( int idx_flow ) const - { - int n_lobes{}; - // Number of lobes in the flow is a random number between the min and max values - if( input.a_beta == 0 && input.b_beta == 0 ) - { - std::uniform_int_distribution<> dist_num_lobes( input.min_n_lobes, input.max_n_lobes ); - n_lobes = dist_num_lobes( gen ); - } - // Deterministic number of lobes, such that a beta probability density distribution is used (not a beta - // distribution). However this means that n_lobes could potentially be greater than min_n_lobes - else - { - double x_beta = ( 1.0 * idx_flow ) / ( input.n_flows - 1.0 ); - double random_number = Math::beta_pdf( x_beta, input.a_beta, input.b_beta ); - n_lobes = int( - std::round( input.min_n_lobes + 0.5 * ( input.max_n_lobes - input.min_n_lobes ) * random_number ) ); - } - return n_lobes; - } +// Calculate the current thickness of a lobe +inline double compute_current_lobe_thickness( + int idx_lobe, int n_lobes, const Config::InputParams & input, const CommonLobeDimensions & lobe_dimensions ) +{ + const double delta_lobe_thickness + = 2.0 * ( lobe_dimensions.avg_lobe_thickness - lobe_dimensions.thickness_min ) / ( n_lobes - 1.0 ); - // Calculate the current thickness of a lobe - double compute_current_lobe_thickness( int idx_lobe, int n_lobes ) const - { - double delta_lobe_thickness - = 2.0 * ( lobe_dimensions.avg_lobe_thickness - lobe_dimensions.thickness_min ) / ( n_lobes - 1.0 ); + // Calculated for each flow with n_lobes number of lobes + FLOWY_CHECK( delta_lobe_thickness ); - // Calculated for each flow with n_lobes number of lobes - FLOWY_CHECK( delta_lobe_thickness ); + return ( 1.0 - input.thickening_parameter ) * ( lobe_dimensions.thickness_min + idx_lobe * delta_lobe_thickness ); +} - return ( 1.0 - input.thickening_parameter ) - * ( lobe_dimensions.thickness_min + idx_lobe * delta_lobe_thickness ); - } +// calculates the initial lobe position +inline void +compute_initial_lobe_position( int idx_flow, Lobe & lobe, const Config::InputParams & input, std::mt19937 & gen ) +{ + std::unique_ptr f{}; - // calculates the initial lobe position - void compute_initial_lobe_position( int idx_flow, Lobe & lobe ) const + // Initial lobes are on the vent and flow starts from the first vent, second vent and so on + if( input.vent_flag == 0 ) { - std::unique_ptr f{}; - - // Initial lobes are on the vent and flow starts from the first vent, second vent and so on - if( input.vent_flag == 0 ) - { - f = std::make_unique( idx_flow, input.n_flows, input.vent_coordinates, gen ); - } - else if( input.vent_flag == 1 ) - { - f = std::make_unique( - input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); - } - else if( input.vent_flag == 2 ) - { - f = std::make_unique( - input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); - } - else if( input.vent_flag == 3 ) - { - f = std::make_unique( - input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); - } - else if( input.vent_flag == 4 ) - { - f = std::make_unique( - input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); - } - else if( input.vent_flag == 5 ) - { - f = std::make_unique( - input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); - } - else if( input.vent_flag == 6 ) - { - f = std::make_unique( - input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); - } - else if( input.vent_flag == 7 ) - { - f = std::make_unique( - input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); - } - else if( input.vent_flag == 8 ) - { - f = std::make_unique( - input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); - } - - lobe.center = f->get_position(); + f = std::make_unique( idx_flow, input.n_flows, input.vent_coordinates, gen ); } - - // perturbes the initial azimuthal angle of the lobe, which is - // computed from the terrain slope - void compute_lobe_axes( Lobe & lobe, double slope ) const + else if( input.vent_flag == 1 ) { - // Factor for the lobe eccentricity - double aspect_ratio = std::min( input.max_aspect_ratio, 1.0 + input.aspect_ratio_coeff * slope ); - - // Compute the semi-axes of the lobe - double semi_major_axis = std::sqrt( lobe_dimensions.lobe_area / Math::pi ) * std::sqrt( aspect_ratio ); - double semi_minor_axis = std::sqrt( lobe_dimensions.lobe_area / Math::pi ) / std::sqrt( aspect_ratio ); - // Set the semi-axes - lobe.semi_axes = { semi_major_axis, semi_minor_axis }; + f = std::make_unique( + input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); } - - void compute_descendent_lobe_position( Lobe & lobe, const Lobe & parent, const Vector2 & final_budding_point ) const + else if( input.vent_flag == 2 ) { - Vector2 direction_to_new_lobe - = ( final_budding_point - parent.center ) / xt::linalg::norm( final_budding_point - parent.center ); - Vector2 new_lobe_center = final_budding_point + input.dist_fact * direction_to_new_lobe * lobe.semi_axes[0]; - lobe.center = new_lobe_center; + f = std::make_unique( + input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); } - - void perturb_lobe_angle( Lobe & lobe, double slope ) const + else if( input.vent_flag == 3 ) { - const double slope_deg = 180.0 / Math::pi * std::atan( slope ); - - if( input.max_slope_prob < 1 ) - { - if( slope_deg > 0.0 && input.max_slope_prob > 0 ) - { - // Since we use radians instead of degrees, max_slope_prob has to be rescaled accordingly - const double sigma = ( 1.0 - input.max_slope_prob ) / input.max_slope_prob * ( 90.0 - slope_deg ) - / slope_deg * Math::pi / 180.0; - - ProbabilityDist::truncated_normal_distribution dist_truncated( 0, sigma, -Math::pi, Math::pi ); - const double angle_perturbation = dist_truncated( gen ); - lobe.set_azimuthal_angle( lobe.get_azimuthal_angle() + angle_perturbation ); - } - else - { - std::uniform_real_distribution dist_uniform( -Math::pi / 2, Math::pi / 2 ); - const double angle_perturbation = dist_uniform( gen ); - lobe.set_azimuthal_angle( lobe.get_azimuthal_angle() + angle_perturbation ); - } - } + f = std::make_unique( + input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); } - - // Select which lobe amongst the existing lobes will be the parent for the new descendent lobe - int select_parent_lobe( int idx_descendant, std::vector & lobes ) const + else if( input.vent_flag == 4 ) + { + f = std::make_unique( + input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); + } + else if( input.vent_flag == 5 ) + { + f = std::make_unique( + input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); + } + else if( input.vent_flag == 6 ) + { + f = std::make_unique( + input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); + } + else if( input.vent_flag == 7 ) { - Lobe & lobe_descendent = lobes[idx_descendant]; + f = std::make_unique( + input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); + } + else if( input.vent_flag == 8 ) + { + f = std::make_unique( + input.vent_coordinates, input.fissure_probabilities, input.fissure_end_coordinates, gen ); + } - int idx_parent{}; + lobe.center = f->get_position(); +} - // Generate from the last lobe - if( input.lobe_exponent <= 0 ) - { - idx_parent = idx_descendant - 1; - } - else if( input.lobe_exponent >= 1 ) // Draw from a uniform random distribution if exponent is 1 +// perturbes the initial azimuthal angle of the lobe, which is +// computed from the terrain slope +inline void compute_lobe_axes( + Lobe & lobe, double slope, const Config::InputParams & input, const CommonLobeDimensions & lobe_dimensions ) +{ + // Factor for the lobe eccentricity + const double aspect_ratio = std::min( input.max_aspect_ratio, 1.0 + input.aspect_ratio_coeff * slope ); + + // Compute the semi-axes of the lobe + const double semi_major_axis = std::sqrt( lobe_dimensions.lobe_area / Math::pi ) * std::sqrt( aspect_ratio ); + const double semi_minor_axis = std::sqrt( lobe_dimensions.lobe_area / Math::pi ) / std::sqrt( aspect_ratio ); + // Set the semi-axes + lobe.semi_axes = { semi_major_axis, semi_minor_axis }; +} + +inline void compute_descendent_lobe_position( + Lobe & lobe, const Lobe & parent, const Vector2 & final_budding_point, const Config::InputParams & input ) +{ + const Vector2 direction_to_new_lobe + = ( final_budding_point - parent.center ) / xt::linalg::norm( final_budding_point - parent.center ); + const Vector2 new_lobe_center = final_budding_point + input.dist_fact * direction_to_new_lobe * lobe.semi_axes[0]; + lobe.center = new_lobe_center; +} + +inline void perturb_lobe_angle( Lobe & lobe, double slope, const Config::InputParams & input, std::mt19937 & gen ) +{ + const double slope_deg = 180.0 / Math::pi * std::atan( slope ); + + if( input.max_slope_prob < 1 ) + { + if( slope_deg > 0.0 && input.max_slope_prob > 0 ) { - std::uniform_int_distribution dist_int( 0, idx_descendant - 1 ); - idx_parent = dist_int( gen ); + // Since we use radians instead of degrees, max_slope_prob has to be rescaled accordingly + const double sigma = ( 1.0 - input.max_slope_prob ) / input.max_slope_prob * ( 90.0 - slope_deg ) + / slope_deg * Math::pi / 180.0; + + ProbabilityDist::truncated_normal_distribution dist_truncated( 0, sigma, -Math::pi, Math::pi ); + const double angle_perturbation = dist_truncated( gen ); + lobe.set_azimuthal_angle( lobe.get_azimuthal_angle() + angle_perturbation ); } else { - std::uniform_real_distribution dist( 0, 1 ); - const double idx0 = dist( gen ); - const auto idx1 = std::pow( idx0, input.lobe_exponent ); - idx_parent = idx_descendant * idx1; + std::uniform_real_distribution dist_uniform( -Math::pi / 2, Math::pi / 2 ); + const double angle_perturbation = dist_uniform( gen ); + lobe.set_azimuthal_angle( lobe.get_azimuthal_angle() + angle_perturbation ); } + } +} + +// Select which lobe amongst the existing lobes will be the parent for the new descendent lobe +inline int select_parent_lobe( + int idx_descendant, std::vector & lobes, const Config::InputParams & input, + const CommonLobeDimensions & lobe_dimensions, std::mt19937 & gen ) +{ + Lobe & lobe_descendent = lobes[idx_descendant]; - // Update the lobe information - lobe_descendent.idx_parent = idx_parent; - lobe_descendent.dist_n_lobes = lobes[idx_parent].dist_n_lobes + 1; - lobe_descendent.parent_weight *= lobe_dimensions.exp_lobe_exponent; + int idx_parent{}; - return idx_parent; + // Generate from the last lobe + if( input.lobe_exponent <= 0 ) + { + idx_parent = idx_descendant - 1; } - - void add_inertial_contribution( Lobe & lobe, const Lobe & parent, double slope ) const + else if( input.lobe_exponent >= 1 ) // Draw from a uniform random distribution if exponent is 1 + { + std::uniform_int_distribution dist_int( 0, idx_descendant - 1 ); + idx_parent = dist_int( gen ); + } + else { - double cos_angle_parent = parent.get_cos_azimuthal_angle(); - double sin_angle_parent = parent.get_sin_azimuthal_angle(); - double cos_angle_lobe = lobe.get_cos_azimuthal_angle(); - double sin_angle_lobe = lobe.get_sin_azimuthal_angle(); + std::uniform_real_distribution dist( 0, 1 ); + const double idx0 = dist( gen ); + const auto idx1 = std::pow( idx0, input.lobe_exponent ); + idx_parent = idx_descendant * idx1; + } - double alpha_inertial = 0.0; + // Update the lobe information + lobe_descendent.idx_parent = idx_parent; + lobe_descendent.dist_n_lobes = lobes[idx_parent].dist_n_lobes + 1; + lobe_descendent.parent_weight *= lobe_dimensions.exp_lobe_exponent; - const double eta = input.inertial_exponent; - if( eta > 0 ) - { - alpha_inertial = std::pow( ( 1.0 - std::pow( 2.0 * std::atan( slope ) / Math::pi, eta ) ), ( 1.0 / eta ) ); - } + return idx_parent; +} + +inline void +add_inertial_contribution( Lobe & lobe, const Lobe & parent, double slope, const Config::InputParams & input ) +{ + const double cos_angle_parent = parent.get_cos_azimuthal_angle(); + const double sin_angle_parent = parent.get_sin_azimuthal_angle(); + const double cos_angle_lobe = lobe.get_cos_azimuthal_angle(); + const double sin_angle_lobe = lobe.get_sin_azimuthal_angle(); - const double x_avg = ( 1.0 - alpha_inertial ) * cos_angle_lobe + alpha_inertial * cos_angle_parent; - const double y_avg = ( 1.0 - alpha_inertial ) * sin_angle_lobe + alpha_inertial * sin_angle_parent; + double alpha_inertial = 0.0; - lobe.set_azimuthal_angle( std::atan2( y_avg, x_avg ) ); + const double eta = input.inertial_exponent; + if( eta > 0 ) + { + alpha_inertial = std::pow( ( 1.0 - std::pow( 2.0 * std::atan( slope ) / Math::pi, eta ) ), ( 1.0 / eta ) ); } -private: - std::mt19937 & gen; -}; + const double x_avg = ( 1.0 - alpha_inertial ) * cos_angle_lobe + alpha_inertial * cos_angle_parent; + const double y_avg = ( 1.0 - alpha_inertial ) * sin_angle_lobe + alpha_inertial * sin_angle_parent; + + lobe.set_azimuthal_angle( std::atan2( y_avg, x_avg ) ); +} + +} // namespace MrLavaLoba } // namespace Flowy diff --git a/flowy/include/simulation.hpp b/flowy/include/simulation.hpp index 0d72955..ee8d7fc 100644 --- a/flowy/include/simulation.hpp +++ b/flowy/include/simulation.hpp @@ -5,21 +5,51 @@ #include "flowy/include/config.hpp" #include "flowy/include/definitions.hpp" #include "flowy/include/lobe.hpp" +#include "flowy/include/models/mr_lava_loba.hpp" #include "flowy/include/topography.hpp" #include "flowy/include/topography_file.hpp" +#include #include #include +#include #include #include namespace Flowy { -class MrLavaLoba; + +/** + * @brief Track the current state of the Simulation run + * + */ +struct SimulationState +{ + std::chrono::time_point t_run_start{}; + std::vector lobes{}; + + int n_lobes_processed = 0; + int n_lobes = 0; + int idx_flow = 0; + int idx_lobe = 0; + + bool beginning_of_new_flow = true; +}; + +enum class FlowStatus +{ + Finished, + Ongoing +}; + +enum class RunStatus +{ + Finished, + Ongoing +}; class Simulation { public: - friend class MrLavaLoba; Simulation( const Config::InputParams & input, std::optional rng_seed ); Config::InputParams input; @@ -39,20 +69,42 @@ class Simulation void write_avg_thickness_file(); - // Check if the dem has to be written (because of the input.write_dem_every_n_lobes_setting) and, if yes, writes the topography - void write_thickness_if_necessary(int n_lobes_processed); + // Check if the dem has to be written (because of the input.write_dem_every_n_lobes_setting) and, if yes, writes the + // topography + void write_thickness_if_necessary( int n_lobes_processed ); - // Computes the topography_thickness field by substacting the initial topography and dividing by (1.0 - filling_parameter) + // Computes the topography_thickness field by subtracting the initial topography and dividing by (1.0 - filling_parameter) void compute_topography_thickness(); std::unique_ptr get_file_handle( const Topography & topography, OutputQuantity output_quantity ) const; + void save_post_run_output(); + void run(); + // Perform `n_steps` steps of the simulation (per step, a single lobe is added to the topography) + RunStatus steps( int n_steps ); + + void reset_simulation_state() + { + simulation_state = std::nullopt; + } + + std::optional get_simulation_state() const + { + return simulation_state; + } + private: int rng_seed; std::mt19937 gen{}; + std::optional simulation_state = std::nullopt; + + void post_flow_hook( int idx_flow, std::vector & lobes ); + void process_initial_lobe( int idx_flow, Lobe & lobe_cur, const CommonLobeDimensions & common_lobe_dimensions ); + FlowStatus process_descendent_lobe( + int idx_lobe, std::vector & lobes, const CommonLobeDimensions & common_lobe_dimensions ); }; } // namespace Flowy diff --git a/meson.build b/meson.build index 3ce9ca0..8bef871 100644 --- a/meson.build +++ b/meson.build @@ -39,8 +39,9 @@ _sources = [ _deps = [] pdflib_dep = dependency('pdf_cpplib', fallback : ['pdf_cpplib', 'pdflib_dep']) +with_netcdf = get_option('with_netcdf') ## Netcdf -if get_option('with_netcdf') +if with_netcdf # On windows the conda-forge netcdf pkg-config file tells us to link against a mysterious debug.lib # Therefore, we use cmake as the preferred method, which does not make problems netcdf_dep = dependency('netcdf', language : 'cpp', method: 'cmake', required: false) diff --git a/src/simulation.cpp b/src/simulation.cpp index 7f53e8c..ddc9b65 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -2,6 +2,7 @@ // Copyright 2023--present Flowy developers #include "flowy/include/simulation.hpp" #include "flowy/include/asc_file.hpp" +#include "flowy/include/config.hpp" #include "flowy/include/definitions.hpp" #include "flowy/include/lobe.hpp" #include "flowy/include/math.hpp" @@ -15,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -345,7 +347,7 @@ Simulation::get_file_handle( const Topography & topography, OutputQuantity outpu void Simulation::compute_topography_thickness() { - // Compute the thickness by substracting the initial topography and correcting for the thickening parametr + // Compute the thickness by subtracting the initial topography and correcting for the thickening parameter topography_thickness = topography; topography_thickness.no_data_value = DEFAULT_NO_DATA_VALUE_THICKNESS; topography_thickness.height_data -= topography_initial.height_data; @@ -370,189 +372,268 @@ void Simulation::write_thickness_if_necessary( int n_lobes_processed ) } } -void Simulation::run() +void Simulation::save_post_run_output() { - // Initialize MrLavaLoba method - auto mr_lava_loba = MrLavaLoba( input, gen ); + // Save initial topography to asc file + auto file_initial = get_file_handle( topography_initial, OutputQuantity::Height ); + file_initial->save( input.output_folder / fmt::format( "{}_DEM", input.run_name ) ); - int n_lobes_processed = 0; + // Save final topography to asc file + if( input.save_final_dem ) + { + auto file_final = get_file_handle( topography, OutputQuantity::Height ); + file_final->save( input.output_folder / fmt::format( "{}_DEM_final", input.run_name ) ); + } - // Make a copy of the initial topography - auto t_run_start = std::chrono::high_resolution_clock::now(); + // Save full thickness to asc file + compute_topography_thickness(); + auto file_thick = get_file_handle( topography_thickness, OutputQuantity::Height ); + file_thick->save( input.output_folder / fmt::format( "{}_thickness_full", input.run_name ) ); - for( int idx_flow = 0; idx_flow < input.n_flows; idx_flow++ ) + // Save the full hazard map + if( input.save_hazard_data ) { - // Determine n_lobes - int n_lobes = mr_lava_loba.compute_n_lobes( idx_flow ); + auto file_hazard = get_file_handle( topography, OutputQuantity::Hazard ); + file_hazard->save( input.output_folder / fmt::format( "{}_hazard_full", input.run_name ) ); + } - lobes = std::vector{}; - lobes.reserve( n_lobes ); + write_avg_thickness_file(); +} - // set the intersection cache - topography.reset_intersection_cache( n_lobes ); +void Simulation::process_initial_lobe( + int idx_flow, Lobe & lobe_cur, const CommonLobeDimensions & common_lobe_dimensions ) +{ + MrLavaLoba::compute_initial_lobe_position( idx_flow, lobe_cur, input, gen ); - // Build initial lobes which do not propagate descendents - for( int idx_lobe = 0; idx_lobe < input.n_init; idx_lobe++ ) - { - lobes.emplace_back(); - Lobe & lobe_cur = lobes.back(); + const auto [height_lobe_center, slope] = topography.height_and_slope( lobe_cur.center ); - mr_lava_loba.compute_initial_lobe_position( idx_flow, lobe_cur ); + if( height_lobe_center == topography.no_data_value ) + { + throw std::runtime_error( + "The initial lobe center has been placed on a no_data value point in the topography." ); + } - // Compute the thickness of the lobe - lobe_cur.thickness = mr_lava_loba.compute_current_lobe_thickness( idx_lobe, n_lobes ); + // Perturb the angle (and set it) + lobe_cur.set_azimuthal_angle( std::atan2( slope[1], slope[0] ) ); // Sets the angle prior to perturbation + const double slope_norm = xt::linalg::norm( slope, 2 ); // Similar to np.linalg.norm + MrLavaLoba::perturb_lobe_angle( lobe_cur, slope_norm, input, gen ); - auto [height_lobe_center, slope] = topography.height_and_slope( lobe_cur.center ); + // compute lobe axes + MrLavaLoba::compute_lobe_axes( lobe_cur, slope_norm, input, common_lobe_dimensions ); +} - if( height_lobe_center == topography.no_data_value ) - { - throw std::runtime_error( - "The initial lobe center has been placed on a no_data value point in the topography." ); - } +[[nodiscard]] FlowStatus Simulation::process_descendent_lobe( + int idx_lobe, std::vector & lobes, const CommonLobeDimensions & common_lobe_dimensions ) +{ + Lobe & lobe_cur = lobes.back(); - // Perturb the angle (and set it) - lobe_cur.set_azimuthal_angle( std::atan2( slope[1], slope[0] ) ); // Sets the angle prior to perturbation - const double slope_norm = xt::linalg::norm( slope, 2 ); // Similar to np.linalg.norm - mr_lava_loba.perturb_lobe_angle( lobe_cur, slope_norm ); + // Select which of the previously created lobes is the parent lobe + // from which the new descendent lobe will bud + const auto idx_parent = MrLavaLoba::select_parent_lobe( idx_lobe, lobes, input, common_lobe_dimensions, gen ); + const Lobe & lobe_parent = lobes[idx_parent]; - // compute lobe axes - mr_lava_loba.compute_lobe_axes( lobe_cur, slope_norm ); + // stopping condition (parent lobe close the domain boundary or at a not defined z value) + if( stop_condition( lobe_parent.center, lobe_parent.semi_axes[0] ) ) + { + lobes.pop_back(); + return FlowStatus::Finished; + } - // Add rasterized lobe - topography.add_lobe( lobe_cur, input.volume_correction, idx_lobe ); - n_lobes_processed++; - write_thickness_if_necessary( n_lobes_processed ); - } + // Find the preliminary budding point on the perimeter of the parent lobe (npoints is the number of raster + // points on the ellipse) + const Flowy::Vector2 budding_point = topography.find_preliminary_budding_point( lobe_parent, input.npoints ); - // Loop over the rest of the lobes (skipping the initial ones). - // Each lobe is a descendant of a parent lobe - for( int idx_lobe = input.n_init; idx_lobe < n_lobes; idx_lobe++ ) - { - lobes.emplace_back(); - Lobe & lobe_cur = lobes.back(); + const auto [height_parent, slope_parent] = topography.height_and_slope( lobe_parent.center ); - // Select which of the previously created lobes is the parent lobe - // from which the new descendent lobe will bud - auto idx_parent = mr_lava_loba.select_parent_lobe( idx_lobe, lobes ); - Lobe & lobe_parent = lobes[idx_parent]; + const Vector2 diff = ( budding_point - lobe_parent.center ); - // stopping condition (parent lobe close the domain boundary or at a not defined z value) - if( stop_condition( lobe_parent.center, lobe_parent.semi_axes[0] ) ) - { - lobes.pop_back(); - break; - } + // Perturb the angle and set it (not on the parent anymore) + lobe_cur.set_azimuthal_angle( std::atan2( diff[1], diff[0] ) ); // Sets the angle prior to perturbation + const double slope_parent_norm = topography.slope_between_points( lobe_parent.center, budding_point ); + MrLavaLoba::perturb_lobe_angle( lobe_cur, slope_parent_norm, input, gen ); - // Find the preliminary budding point on the perimeter of the parent lobe (npoints is the number of raster - // points on the ellipse) - Flowy::Vector2 budding_point = topography.find_preliminary_budding_point( lobe_parent, input.npoints ); + // Add the inertial contribution + MrLavaLoba::add_inertial_contribution( lobe_cur, lobe_parent, slope_parent_norm, input ); - auto [height_lobe_center, slope_parent] = topography.height_and_slope( lobe_parent.center ); - auto [height_bp, slope_bp] = topography.height_and_slope( budding_point ); + // Compute the final budding point + // It is defined by the point on the perimeter of the parent lobe closest to the center of the new lobe + const auto angle_diff = lobe_parent.get_azimuthal_angle() - lobe_cur.get_azimuthal_angle(); + const Vector2 final_budding_point = lobe_parent.point_at_angle( -angle_diff ); - Vector2 diff = ( budding_point - lobe_parent.center ); + // final_budding_point = budding_point; + if( stop_condition( final_budding_point, lobe_parent.semi_axes[0] ) ) + { + lobes.pop_back(); + return FlowStatus::Finished; + } + // Get the slope at the final budding point + const double slope_budding_point = topography.slope_between_points( lobe_parent.center, final_budding_point ); - // Perturb the angle and set it (not on the parent anymore) - lobe_cur.set_azimuthal_angle( std::atan2( diff[1], diff[0] ) ); // Sets the angle prior to perturbation - const double slope_parent_norm = topography.slope_between_points( lobe_parent.center, budding_point ); - mr_lava_loba.perturb_lobe_angle( lobe_cur, slope_parent_norm ); + // compute the new lobe axes + MrLavaLoba::compute_lobe_axes( lobe_cur, slope_budding_point, input, common_lobe_dimensions ); - // Add the inertial contribution - mr_lava_loba.add_inertial_contribution( lobe_cur, lobe_parent, slope_parent_norm ); + // Get new lobe center + MrLavaLoba::compute_descendent_lobe_position( lobe_cur, lobe_parent, final_budding_point, input ); - // Compute the final budding point - // It is defined by the point on the perimeter of the parent lobe closest to the center of the new lobe - auto angle_diff = lobe_parent.get_azimuthal_angle() - lobe_cur.get_azimuthal_angle(); - Vector2 final_budding_point = lobe_parent.point_at_angle( -angle_diff ); + if( stop_condition( lobe_cur.center, lobe_cur.semi_axes[0] ) ) + { + lobes.pop_back(); + return FlowStatus::Finished; + } - // final_budding_point = budding_point; - if( stop_condition( final_budding_point, lobe_parent.semi_axes[0] ) ) - { - lobes.pop_back(); - break; - } - // Get the slope at the final budding point - double slope_budding_point = topography.slope_between_points( lobe_parent.center, final_budding_point ); + return FlowStatus::Ongoing; +} - // compute the new lobe axes - mr_lava_loba.compute_lobe_axes( lobe_cur, slope_budding_point ); +void Simulation::post_flow_hook( int idx_flow, std::vector & lobes ) +{ + if( input.save_hazard_data ) + { + compute_cumulative_descendents( lobes ); + topography.compute_hazard_flow( lobes ); + } - // Get new lobe center - mr_lava_loba.compute_descendent_lobe_position( lobe_cur, lobe_parent, final_budding_point ); + if( input.write_lobes_csv ) + { + write_lobe_data_to_file( lobes, input.output_folder / fmt::format( "lobes_{}.csv", idx_flow ) ); + } - if( stop_condition( lobe_cur.center, lobe_cur.semi_axes[0] ) ) - { - lobes.pop_back(); - break; - } + if( input.print_remaining_time ) + { + const std::chrono::time_point t_cur + = std::chrono::high_resolution_clock::now(); - // Compute the thickness of the lobe - lobe_cur.thickness = mr_lava_loba.compute_current_lobe_thickness( idx_lobe, n_lobes ); + const auto remaining_time = std::chrono::duration_cast( + ( input.n_flows - idx_flow - 1 ) * ( t_cur - simulation_state->t_run_start ) / ( idx_flow + 1 ) ); + fmt::print( " remaining_time = {:%Hh %Mm %Ss}\n", remaining_time ); + } +} + +RunStatus Simulation::steps( int n_steps ) +{ + RunStatus run_status = RunStatus::Ongoing; + + // If the simulation state has no value, we initialize a new simulation + if( !simulation_state.has_value() ) + { + simulation_state = SimulationState(); + simulation_state->t_run_start = std::chrono::high_resolution_clock::now(); + } + const auto common_lobe_dimensions = CommonLobeDimensions( input ); + + // Define some aliases for convenience + auto & lobes = simulation_state->lobes; + auto & idx_flow = simulation_state->idx_flow; + auto & idx_lobe = simulation_state->idx_lobe; + auto & n_lobes_processed = simulation_state->n_lobes_processed; + auto & n_lobes = simulation_state->n_lobes; + + int n_lobes_processed_initial = n_lobes_processed; + + // Now we have to figure out the value of idx_flow and if we are past n_init or not + while( n_lobes_processed < n_lobes_processed_initial + n_steps ) + { + bool is_last_step_of_flow = idx_lobe == n_lobes - 1; + const bool is_an_initial_lobe = idx_lobe < input.n_init; + const bool is_last_step = ( idx_flow == input.n_flows - 1 ) && is_last_step_of_flow; + + if( simulation_state->beginning_of_new_flow ) + { + simulation_state->n_lobes = MrLavaLoba::compute_n_lobes( simulation_state->idx_flow, input, gen ); + simulation_state->lobes = std::vector{}; + simulation_state->lobes.reserve( simulation_state->n_lobes ); + // set the intersection cache + topography.reset_intersection_cache( simulation_state->n_lobes ); + + // set idx_lobe to zero, because we are at the beginning of a flow + idx_lobe = 0; + + // switch back the flag + simulation_state->beginning_of_new_flow = false; + } - // Add rasterized lobe + if( idx_lobe >= n_lobes || idx_flow >= input.n_flows ) + { + run_status = RunStatus::Finished; + break; + } + + // Add a new lobe + lobes.emplace_back(); + Lobe & lobe_cur = lobes.back(); + // Compute the thickness of the lobe + lobe_cur.thickness + = MrLavaLoba::compute_current_lobe_thickness( idx_lobe, n_lobes, input, common_lobe_dimensions ); + + if( is_an_initial_lobe ) + { + process_initial_lobe( idx_flow, lobe_cur, common_lobe_dimensions ); + + // Rasterized lobe and add it to the topography topography.add_lobe( lobe_cur, input.volume_correction, idx_lobe ); n_lobes_processed++; write_thickness_if_necessary( n_lobes_processed ); } - - if( input.save_hazard_data ) + else { - compute_cumulative_descendents( lobes ); - topography.compute_hazard_flow( lobes ); + const auto flow_status = process_descendent_lobe( idx_lobe, lobes, common_lobe_dimensions ); + if( flow_status == FlowStatus::Ongoing ) + { + // Rasterized lobe and add it to the topography + topography.add_lobe( lobe_cur, input.volume_correction, idx_lobe ); + n_lobes_processed++; + write_thickness_if_necessary( n_lobes_processed ); + } + else + { + is_last_step_of_flow = true; + } } - if( input.write_lobes_csv ) + idx_lobe++; + + if( is_last_step_of_flow ) { - write_lobe_data_to_file( lobes, input.output_folder / fmt::format( "lobes_{}.csv", idx_flow ) ); + post_flow_hook( idx_flow, lobes ); + simulation_state->beginning_of_new_flow = true; + idx_flow++; + idx_lobe = 0; } - if( input.print_remaining_time ) + if( is_last_step ) { - auto t_cur = std::chrono::high_resolution_clock::now(); - auto remaining_time = std::chrono::duration_cast( - ( input.n_flows - idx_flow - 1 ) * ( t_cur - t_run_start ) / ( idx_flow + 1 ) ); - fmt::print( " remaining_time = {:%Hh %Mm %Ss}\n", remaining_time ); + run_status = RunStatus::Finished; } } - auto t_cur = std::chrono::high_resolution_clock::now(); - auto total_time = std::chrono::duration_cast( ( t_cur - t_run_start ) ); - fmt::print( "total_time = {:%Hh %Mm %Ss}\n", total_time ); - - fmt::print( "Total number of processed lobes = {}\n", n_lobes_processed ); - - if( total_time.count() > 0 ) + if( run_status == RunStatus::Finished ) { - auto lobes_per_ms = n_lobes_processed / total_time.count(); - fmt::print( "n_lobes/ms = {}\n", lobes_per_ms ); - } + const auto t_cur = std::chrono::high_resolution_clock::now(); + const auto total_time + = std::chrono::duration_cast( ( t_cur - simulation_state->t_run_start ) ); + fmt::print( "total_time = {:%Hh %Mm %Ss}\n", total_time ); - fmt::print( "Used RNG seed: {}\n", rng_seed ); + fmt::print( "Total number of processed lobes = {}\n", n_lobes_processed ); - // Save initial topography to asc file - auto file_initial = get_file_handle( topography_initial, OutputQuantity::Height ); - file_initial->save( input.output_folder / fmt::format( "{}_DEM", input.run_name ) ); + if( total_time.count() > 0 ) + { + const auto lobes_per_ms = n_lobes_processed / total_time.count(); + fmt::print( "n_lobes/ms = {}\n", lobes_per_ms ); + } - // Save final topography to asc file - if( input.save_final_dem ) - { - auto file_final = get_file_handle( topography, OutputQuantity::Height ); - file_final->save( input.output_folder / fmt::format( "{}_DEM_final", input.run_name ) ); + fmt::print( "Used RNG seed: {}\n", rng_seed ); + save_post_run_output(); } - // Save full thickness to asc file - compute_topography_thickness(); - auto file_thick = get_file_handle( topography_thickness, OutputQuantity::Height ); - file_thick->save( input.output_folder / fmt::format( "{}_thickness_full", input.run_name ) ); + return run_status; +} - // Save the full hazard map - if( input.save_hazard_data ) +void Simulation::run() +{ + // This number is completely arbitrary, it should just be big enough to mitigate the overhead of creating + // CommonLobeDimensions + constexpr int NSTEPS = 1000; + while( steps( NSTEPS ) != RunStatus::Finished ) { - auto file_hazard = get_file_handle( topography, OutputQuantity::Hazard ); - file_hazard->save( input.output_folder / fmt::format( "{}_hazard_full", input.run_name ) ); - } - - write_avg_thickness_file(); + }; } - } // namespace Flowy diff --git a/src/topography.cpp b/src/topography.cpp index 62df190..cd6ca0b 100644 --- a/src/topography.cpp +++ b/src/topography.cpp @@ -503,8 +503,9 @@ void Topography::add_lobe( const Lobe & lobe, bool volume_correction, std::optio Vector2 Topography::find_preliminary_budding_point( const Lobe & lobe, size_t npoints ) { - bool compute_cache = cos_phi_lobe_perimeter == std::nullopt || sin_phi_lobe_perimeter == std::nullopt - || npoints != cos_phi_lobe_perimeter->size() || npoints != sin_phi_lobe_perimeter->size(); + const bool compute_cache = cos_phi_lobe_perimeter == std::nullopt || sin_phi_lobe_perimeter == std::nullopt + || npoints != cos_phi_lobe_perimeter->size() + || npoints != sin_phi_lobe_perimeter->size(); if( compute_cache ) { @@ -516,7 +517,7 @@ Vector2 Topography::find_preliminary_budding_point( const Lobe & lobe, size_t np // First, we rasterize the perimeter of the ellipse const auto sin = std::span( sin_phi_lobe_perimeter->begin(), sin_phi_lobe_perimeter->end() ); const auto cos = std::span( cos_phi_lobe_perimeter->begin(), cos_phi_lobe_perimeter->end() ); - std::vector perimeter = lobe.rasterize_perimeter( sin, cos ); + const std::vector perimeter = lobe.rasterize_perimeter( sin, cos ); // Then, we find the point of minimal elevation amongst the rasterized points on the perimeter auto min_elevation_point_it = std::min_element( diff --git a/test/test_simulation.cpp b/test/test_simulation.cpp index aee0266..c46a24f 100644 --- a/test/test_simulation.cpp +++ b/test/test_simulation.cpp @@ -31,9 +31,8 @@ TEST_CASE( "perturb_angle", "[perturb_angle]" ) input.prescribed_avg_lobe_thickness = 1.0; INFO( fmt::format( "{}", fmt::streamed( input.source ) ) ); - auto simulation = Simulation( input, 0 ); - auto mr_lava_loba = MrLavaLoba( simulation ); - + auto simulation = Simulation( input, 0 ); + std::mt19937 gen{}; Lobe my_lobe; int n_samples = 100000; @@ -53,7 +52,7 @@ TEST_CASE( "perturb_angle", "[perturb_angle]" ) { my_lobe.set_azimuthal_angle( std::atan2( slope[1], slope[0] ) ); // Sets the angle prior to perturbation const double slope_norm = xt::linalg::norm( slope, 2 ); // Similar to np.linalg.norm - mr_lava_loba.perturb_lobe_angle( my_lobe, slope_norm ); + MrLavaLoba::perturb_lobe_angle( my_lobe, slope_norm, input, gen ); angle_samples[i] = my_lobe.get_azimuthal_angle(); REQUIRE( ( ( my_lobe.get_azimuthal_angle() > -Math::pi ) && ( my_lobe.get_azimuthal_angle() < Math::pi ) ) ); @@ -91,9 +90,10 @@ TEST_CASE( "budding_point", "[budding_point]" ) Vector2 final_budding_point = lobe_parent.point_at_angle( 0 ); - auto simulation = Simulation( input_params, std::nullopt ); - auto mr_lava_loba = MrLavaLoba( simulation ); - mr_lava_loba.compute_descendent_lobe_position( lobe_cur, lobe_parent, final_budding_point ); + auto simulation = Simulation( input_params, std::nullopt ); + std::mt19937 gen{}; + + MrLavaLoba::compute_descendent_lobe_position( lobe_cur, lobe_parent, final_budding_point, input_params ); INFO( fmt::format( "budding_point = {}\n", fmt::streamed( final_budding_point ) ) ); INFO( fmt::format( "lobe_cur.center = {}\n", fmt::streamed( lobe_cur.center ) ) );