From 5a487aee80764e417f9dc2b8bcf4936b2135d113 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Wed, 26 Mar 2025 15:24:12 +0000 Subject: [PATCH 01/15] CHR: fix typo in comment --- flowy/include/simulation.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flowy/include/simulation.hpp b/flowy/include/simulation.hpp index 0d72955..af68fc6 100644 --- a/flowy/include/simulation.hpp +++ b/flowy/include/simulation.hpp @@ -42,7 +42,7 @@ class Simulation // 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 From 108952d6cec768a391c80b12cdf37773bc211142 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Wed, 26 Mar 2025 15:42:24 +0000 Subject: [PATCH 02/15] ENH: Turned the MrLavaLoba class into a namespace Co-authored-by: Amrita Goswami --- flowy/include/models/mr_lava_loba.hpp | 336 +++++++++++++------------- flowy/include/simulation.hpp | 2 - src/simulation.cpp | 34 ++- test/test_simulation.cpp | 14 +- 4 files changed, 186 insertions(+), 200 deletions(-) diff --git a/flowy/include/models/mr_lava_loba.hpp b/flowy/include/models/mr_lava_loba.hpp index ee01f2e..f752118 100644 --- a/flowy/include/models/mr_lava_loba.hpp +++ b/flowy/include/models/mr_lava_loba.hpp @@ -63,214 +63,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 min_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 + = int( 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, 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, 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, + 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 af68fc6..82c4be2 100644 --- a/flowy/include/simulation.hpp +++ b/flowy/include/simulation.hpp @@ -14,12 +14,10 @@ namespace Flowy { -class MrLavaLoba; class Simulation { public: - friend class MrLavaLoba; Simulation( const Config::InputParams & input, std::optional rng_seed ); Config::InputParams input; diff --git a/src/simulation.cpp b/src/simulation.cpp index 7f53e8c..60d1a58 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -231,8 +231,7 @@ void Simulation::write_avg_thickness_file() // This lambda performs bisection search to find the threshold thickness at which a // relative volume proportion of `thresh` is contained within cells with greater thickness than the threshold thickness - auto bisection_search = [&]( double thresh, double tol, int max_iter ) - { + auto bisection_search = [&]( double thresh, double tol, int max_iter ) { int idx_lo = 0; int idx_hi = n_cells - 1; @@ -372,18 +371,15 @@ void Simulation::write_thickness_if_necessary( int n_lobes_processed ) void Simulation::run() { - // Initialize MrLavaLoba method - auto mr_lava_loba = MrLavaLoba( input, gen ); - int n_lobes_processed = 0; - // Make a copy of the initial topography - auto t_run_start = std::chrono::high_resolution_clock::now(); + auto t_run_start = std::chrono::high_resolution_clock::now(); + auto common_lobe_dimensions = CommonLobeDimensions( input ); for( int idx_flow = 0; idx_flow < input.n_flows; idx_flow++ ) { // Determine n_lobes - int n_lobes = mr_lava_loba.compute_n_lobes( idx_flow ); + int n_lobes = MrLavaLoba::compute_n_lobes( idx_flow, input, gen ); lobes = std::vector{}; lobes.reserve( n_lobes ); @@ -397,10 +393,11 @@ void Simulation::run() lobes.emplace_back(); Lobe & lobe_cur = lobes.back(); - mr_lava_loba.compute_initial_lobe_position( idx_flow, lobe_cur ); + MrLavaLoba::compute_initial_lobe_position( idx_flow, lobe_cur, input, gen ); // Compute the thickness of the lobe - lobe_cur.thickness = mr_lava_loba.compute_current_lobe_thickness( idx_lobe, n_lobes ); + lobe_cur.thickness + = MrLavaLoba::compute_current_lobe_thickness( idx_lobe, n_lobes, input, common_lobe_dimensions ); auto [height_lobe_center, slope] = topography.height_and_slope( lobe_cur.center ); @@ -413,10 +410,10 @@ void Simulation::run() // 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 ); + MrLavaLoba::perturb_lobe_angle( lobe_cur, slope_norm, input, gen ); // compute lobe axes - mr_lava_loba.compute_lobe_axes( lobe_cur, slope_norm ); + MrLavaLoba::compute_lobe_axes( lobe_cur, slope_norm, input, common_lobe_dimensions ); // Add rasterized lobe topography.add_lobe( lobe_cur, input.volume_correction, idx_lobe ); @@ -433,7 +430,7 @@ void Simulation::run() // 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 ); + auto idx_parent = MrLavaLoba::select_parent_lobe( idx_lobe, lobes, input, common_lobe_dimensions, gen ); Lobe & lobe_parent = lobes[idx_parent]; // stopping condition (parent lobe close the domain boundary or at a not defined z value) @@ -455,10 +452,10 @@ void Simulation::run() // 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 ); + MrLavaLoba::perturb_lobe_angle( lobe_cur, slope_parent_norm, input, gen ); // Add the inertial contribution - mr_lava_loba.add_inertial_contribution( lobe_cur, lobe_parent, slope_parent_norm ); + MrLavaLoba::add_inertial_contribution( lobe_cur, lobe_parent, slope_parent_norm, 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 @@ -475,10 +472,10 @@ void Simulation::run() double slope_budding_point = topography.slope_between_points( lobe_parent.center, final_budding_point ); // compute the new lobe axes - mr_lava_loba.compute_lobe_axes( lobe_cur, slope_budding_point ); + MrLavaLoba::compute_lobe_axes( lobe_cur, slope_budding_point, input, common_lobe_dimensions ); // Get new lobe center - mr_lava_loba.compute_descendent_lobe_position( lobe_cur, lobe_parent, final_budding_point ); + MrLavaLoba::compute_descendent_lobe_position( lobe_cur, lobe_parent, final_budding_point, input ); if( stop_condition( lobe_cur.center, lobe_cur.semi_axes[0] ) ) { @@ -487,7 +484,8 @@ void Simulation::run() } // Compute the thickness of the lobe - lobe_cur.thickness = mr_lava_loba.compute_current_lobe_thickness( idx_lobe, n_lobes ); + lobe_cur.thickness + = MrLavaLoba::compute_current_lobe_thickness( idx_lobe, n_lobes, input, common_lobe_dimensions ); // Add rasterized lobe topography.add_lobe( lobe_cur, input.volume_correction, idx_lobe ); 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 ) ) ); From fa505601b97d5a55a39c6f5938a11bbd4f9fc761 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Wed, 26 Mar 2025 16:01:12 +0000 Subject: [PATCH 03/15] ENH: `save_post_run_output` function a bit more concise Co-authored-by: Amrita Goswami --- flowy/include/models/mr_lava_loba.hpp | 1 - flowy/include/simulation.hpp | 7 +++- src/simulation.cpp | 55 +++++++++++++++------------ 3 files changed, 35 insertions(+), 28 deletions(-) diff --git a/flowy/include/models/mr_lava_loba.hpp b/flowy/include/models/mr_lava_loba.hpp index f752118..0b50150 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 diff --git a/flowy/include/simulation.hpp b/flowy/include/simulation.hpp index 82c4be2..7cd4f3e 100644 --- a/flowy/include/simulation.hpp +++ b/flowy/include/simulation.hpp @@ -37,8 +37,9 @@ 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 subtracting the initial topography and dividing by (1.0 - filling_parameter) void compute_topography_thickness(); @@ -46,6 +47,8 @@ class Simulation std::unique_ptr get_file_handle( const Topography & topography, OutputQuantity output_quantity ) const; + void save_post_run_output(); + void run(); private: diff --git a/src/simulation.cpp b/src/simulation.cpp index 60d1a58..ce47db0 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -344,7 +344,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 parametr topography_thickness = topography; topography_thickness.no_data_value = DEFAULT_NO_DATA_VALUE_THICKNESS; topography_thickness.height_data -= topography_initial.height_data; @@ -369,6 +369,34 @@ void Simulation::write_thickness_if_necessary( int n_lobes_processed ) } } +void Simulation::save_post_run_output() +{ + // 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 ) ); + + // 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 ) ); + } + + // 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 ) ); + + // Save the full hazard map + if( input.save_hazard_data ) + { + 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(); +} + void Simulation::run() { int n_lobes_processed = 0; @@ -527,30 +555,7 @@ void Simulation::run() fmt::print( "Used RNG seed: {}\n", rng_seed ); - // 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 ) ); - - // 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 ) ); - } - - // 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 ) ); - - // Save the full hazard map - if( input.save_hazard_data ) - { - 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(); + save_post_run_output(); } } // namespace Flowy From 54f45377de151bbd0d41d4345d7cb3a9901e8349 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Wed, 26 Mar 2025 16:39:11 +0000 Subject: [PATCH 04/15] Simulation: Started to write `steps` function (WIP) Co-authored-by: Amrita Goswami --- flowy/include/simulation.hpp | 29 ++++++ src/simulation.cpp | 190 +++++++++++++++++++++++++++++++++++ 2 files changed, 219 insertions(+) diff --git a/flowy/include/simulation.hpp b/flowy/include/simulation.hpp index 7cd4f3e..00befbf 100644 --- a/flowy/include/simulation.hpp +++ b/flowy/include/simulation.hpp @@ -7,14 +7,39 @@ #include "flowy/include/lobe.hpp" #include "flowy/include/topography.hpp" #include "flowy/include/topography_file.hpp" +#include #include #include +#include #include #include namespace Flowy { +/** + * @brief Track the current state of the Simulation run + * + */ +struct SimulationState +{ + int n_lobes_processed = 0; + std::chrono::time_point t_run_start{}; + int n_lobes{}; + std::vector lobes{}; + + int step = 0; + int idx_flow = 0; + int idx_lobe = 0; + int n_lobes_current_flow = 0; +}; + +enum class RunStatus +{ + Finished, + Ongoing +}; + class Simulation { public: @@ -51,9 +76,13 @@ class Simulation void run(); + // Perform `n_steps` steps of the simulation (per step, a single lobe is added to the topography) + RunStatus steps( int n_steps ); + private: int rng_seed; std::mt19937 gen{}; + std::optional simulation_state = std::nullopt; }; } // namespace Flowy diff --git a/src/simulation.cpp b/src/simulation.cpp index ce47db0..b682294 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -397,6 +397,196 @@ void Simulation::save_post_run_output() write_avg_thickness_file(); } +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(); + } + 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; + + // Now we have to figure out the value of idx_flow and if we are past n_init or not + const int step_initial = simulation_state->step; + + for( int step = step_initial; step < step_initial + n_steps; step++ ) + { + bool is_first_step_of_new_flow = {}; // TODO: figure this out + bool is_last_step_of_flow = {}; + bool is_an_initial_lobe = {}; + bool is_last_step = {}; + + if( is_first_step_of_new_flow ) + { + // simulation_state->idx_flow++; + simulation_state->n_lobes_current_flow + = MrLavaLoba::compute_n_lobes( simulation_state->idx_flow, input, gen ); + simulation_state->lobes = std::vector{}; + simulation_state->lobes.reserve( simulation_state->n_lobes_current_flow ); + // set the intersection cache + topography.reset_intersection_cache( simulation_state->n_lobes_current_flow ); + + // set idx_lobe to zero, because we are at the beginning of a flow + idx_lobe = 0; + } + + if( is_an_initial_lobe ) + { + lobes.emplace_back(); + Lobe & lobe_cur = lobes.back(); + + MrLavaLoba::compute_initial_lobe_position( idx_flow, lobe_cur, input, gen ); + + // Compute the thickness of the lobe + lobe_cur.thickness + = MrLavaLoba::compute_current_lobe_thickness( idx_lobe, n_lobes, input, common_lobe_dimensions ); + + auto [height_lobe_center, slope] = topography.height_and_slope( lobe_cur.center ); + + 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." ); + } + + // 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 ); + + // compute lobe axes + MrLavaLoba::compute_lobe_axes( lobe_cur, slope_norm, input, common_lobe_dimensions ); + + // Add rasterized lobe + topography.add_lobe( lobe_cur, input.volume_correction, idx_lobe ); + n_lobes_processed++; + write_thickness_if_necessary( n_lobes_processed ); + } + else + { + lobes.emplace_back(); + Lobe & lobe_cur = lobes.back(); + + // Select which of the previously created lobes is the parent lobe + // from which the new descendent lobe will bud + auto idx_parent = MrLavaLoba::select_parent_lobe( idx_lobe, lobes, input, common_lobe_dimensions, gen ); + Lobe & lobe_parent = lobes[idx_parent]; + + // 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(); + run_status = RunStatus::Finished; + break; + } + + // 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 ); + + auto [height_lobe_center, slope_parent] = topography.height_and_slope( lobe_parent.center ); + auto [height_bp, slope_bp] = topography.height_and_slope( budding_point ); + + Vector2 diff = ( budding_point - lobe_parent.center ); + + // 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 ); + + // Add the inertial contribution + MrLavaLoba::add_inertial_contribution( lobe_cur, lobe_parent, slope_parent_norm, 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 ); + + // final_budding_point = budding_point; + if( stop_condition( final_budding_point, lobe_parent.semi_axes[0] ) ) + { + lobes.pop_back(); + run_status = RunStatus::Finished; + break; + } + // Get the slope at the final budding point + double slope_budding_point = topography.slope_between_points( lobe_parent.center, final_budding_point ); + + // compute the new lobe axes + MrLavaLoba::compute_lobe_axes( lobe_cur, slope_budding_point, input, common_lobe_dimensions ); + + // Get new lobe center + MrLavaLoba::compute_descendent_lobe_position( lobe_cur, lobe_parent, final_budding_point, input ); + + if( stop_condition( lobe_cur.center, lobe_cur.semi_axes[0] ) ) + { + lobes.pop_back(); + run_status = RunStatus::Finished; + break; + } + + // Compute the thickness of the lobe + lobe_cur.thickness + = MrLavaLoba::compute_current_lobe_thickness( idx_lobe, n_lobes, input, common_lobe_dimensions ); + + // Add rasterized lobe + topography.add_lobe( lobe_cur, input.volume_correction, idx_lobe ); + n_lobes_processed++; + write_thickness_if_necessary( n_lobes_processed ); + } + + if( is_last_step_of_flow ) + { + // Increment idx_flow + idx_flow++; + + if( input.save_hazard_data ) + { + compute_cumulative_descendents( lobes ); + topography.compute_hazard_flow( lobes ); + } + + if( input.write_lobes_csv ) + { + write_lobe_data_to_file( lobes, input.output_folder / fmt::format( "lobes_{}.csv", idx_flow ) ); + } + + if( input.print_remaining_time ) + { + auto t_cur = std::chrono::high_resolution_clock::now(); + 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 ); + } + } + + if( is_last_step ) + { + run_status = RunStatus::Finished; + } + + n_lobes_processed++; + idx_lobe++; + } + + if( run_status == RunStatus::Finished ) + { + save_post_run_output(); + } + + return run_status; +} + void Simulation::run() { int n_lobes_processed = 0; From 96112dae7bac594539910626e931c8555ee31b2a Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Wed, 26 Mar 2025 20:15:00 +0000 Subject: [PATCH 05/15] ENH: pre-commit: filter out build/c++11 --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From da10e5bf100f77b4236ef9578c614bef70ff5f21 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Wed, 26 Mar 2025 20:15:50 +0000 Subject: [PATCH 06/15] mr_lava_loba: int(...) -> static_cast(...) cpplint was complaining --- flowy/include/models/mr_lava_loba.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/flowy/include/models/mr_lava_loba.hpp b/flowy/include/models/mr_lava_loba.hpp index 0b50150..8ffd423 100644 --- a/flowy/include/models/mr_lava_loba.hpp +++ b/flowy/include/models/mr_lava_loba.hpp @@ -80,8 +80,8 @@ inline int compute_n_lobes( int idx_flow, const Config::InputParams & input, std { 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 - = int( std::round( input.min_n_lobes + 0.5 * ( input.max_n_lobes - input.min_n_lobes ) * random_number ) ); + 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; } From 9328b4f77c9121f47fec8a7b52ce1198b2710475 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Wed, 26 Mar 2025 21:04:25 +0000 Subject: [PATCH 07/15] ENH: simulation: wrap up implementation of `steps` ENH: simulation: better const correctness :) Co-authored-by: Amrita Goswami --- flowy/include/models/mr_lava_loba.hpp | 6 +-- flowy/include/simulation.hpp | 6 ++- src/simulation.cpp | 68 ++++++++++++++------------- 3 files changed, 43 insertions(+), 37 deletions(-) diff --git a/flowy/include/models/mr_lava_loba.hpp b/flowy/include/models/mr_lava_loba.hpp index 8ffd423..979f908 100644 --- a/flowy/include/models/mr_lava_loba.hpp +++ b/flowy/include/models/mr_lava_loba.hpp @@ -88,7 +88,7 @@ inline int compute_n_lobes( int idx_flow, const Config::InputParams & input, std // Calculate the current thickness of a lobe inline double compute_current_lobe_thickness( - int idx_lobe, int n_lobes, const Config::InputParams & input, CommonLobeDimensions & lobe_dimensions ) + 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 ); @@ -157,7 +157,7 @@ compute_initial_lobe_position( int idx_flow, Lobe & lobe, const Config::InputPar // 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, CommonLobeDimensions & lobe_dimensions ) + 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 ); @@ -206,7 +206,7 @@ inline void perturb_lobe_angle( Lobe & lobe, double slope, const Config::InputPa // 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, - CommonLobeDimensions & lobe_dimensions, std::mt19937 & gen ) + const CommonLobeDimensions & lobe_dimensions, std::mt19937 & gen ) { Lobe & lobe_descendent = lobes[idx_descendant]; diff --git a/flowy/include/simulation.hpp b/flowy/include/simulation.hpp index 00befbf..5b6fcc2 100644 --- a/flowy/include/simulation.hpp +++ b/flowy/include/simulation.hpp @@ -23,15 +23,17 @@ namespace Flowy */ struct SimulationState { - int n_lobes_processed = 0; std::chrono::time_point t_run_start{}; - int n_lobes{}; std::vector lobes{}; + int n_lobes_processed = 0; + int n_lobes = 0; int step = 0; int idx_flow = 0; int idx_lobe = 0; int n_lobes_current_flow = 0; + + bool beginning_of_new_flow = true; }; enum class RunStatus diff --git a/src/simulation.cpp b/src/simulation.cpp index b682294..ea0de38 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -407,7 +407,7 @@ RunStatus Simulation::steps( int n_steps ) simulation_state = SimulationState(); simulation_state->t_run_start = std::chrono::high_resolution_clock::now(); } - auto common_lobe_dimensions = CommonLobeDimensions( input ); + const auto common_lobe_dimensions = CommonLobeDimensions( input ); // Define some aliases for convenience auto & lobes = simulation_state->lobes; @@ -421,10 +421,10 @@ RunStatus Simulation::steps( int n_steps ) for( int step = step_initial; step < step_initial + n_steps; step++ ) { - bool is_first_step_of_new_flow = {}; // TODO: figure this out - bool is_last_step_of_flow = {}; - bool is_an_initial_lobe = {}; - bool is_last_step = {}; + const bool is_last_step_of_flow = idx_lobe == n_lobes - 1; + const bool is_first_step_of_new_flow = simulation_state->beginning_of_new_flow; + const bool is_an_initial_lobe = idx_lobe < input.n_init; + const bool is_last_step = ( idx_flow == input.n_flows ) && is_last_step_of_flow; if( is_first_step_of_new_flow ) { @@ -438,6 +438,9 @@ RunStatus Simulation::steps( int n_steps ) // 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; } if( is_an_initial_lobe ) @@ -451,7 +454,7 @@ RunStatus Simulation::steps( int n_steps ) lobe_cur.thickness = MrLavaLoba::compute_current_lobe_thickness( idx_lobe, n_lobes, input, common_lobe_dimensions ); - auto [height_lobe_center, slope] = topography.height_and_slope( lobe_cur.center ); + const auto [height_lobe_center, slope] = topography.height_and_slope( lobe_cur.center ); if( height_lobe_center == topography.no_data_value ) { @@ -479,8 +482,8 @@ RunStatus Simulation::steps( int n_steps ) // Select which of the previously created lobes is the parent lobe // from which the new descendent lobe will bud - auto idx_parent = MrLavaLoba::select_parent_lobe( idx_lobe, lobes, input, common_lobe_dimensions, gen ); - Lobe & lobe_parent = lobes[idx_parent]; + const auto idx_parent = MrLavaLoba::select_parent_lobe( idx_lobe, lobes, input, common_lobe_dimensions, gen ); + const Lobe & lobe_parent = lobes[idx_parent]; // 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] ) ) @@ -492,12 +495,12 @@ RunStatus Simulation::steps( int n_steps ) // 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 ); + const Flowy::Vector2 budding_point = topography.find_preliminary_budding_point( lobe_parent, input.npoints ); - auto [height_lobe_center, slope_parent] = topography.height_and_slope( lobe_parent.center ); - auto [height_bp, slope_bp] = topography.height_and_slope( budding_point ); + const auto [height_lobe_center, slope_parent] = topography.height_and_slope( lobe_parent.center ); + const auto [height_bp, slope_bp] = topography.height_and_slope( budding_point ); - Vector2 diff = ( budding_point - lobe_parent.center ); + const Vector2 diff = ( budding_point - lobe_parent.center ); // 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 @@ -509,8 +512,8 @@ RunStatus Simulation::steps( int n_steps ) // 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 ); + 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 ); // final_budding_point = budding_point; if( stop_condition( final_budding_point, lobe_parent.semi_axes[0] ) ) @@ -520,7 +523,7 @@ RunStatus Simulation::steps( int n_steps ) break; } // Get the slope at the final budding point - double slope_budding_point = topography.slope_between_points( lobe_parent.center, final_budding_point ); + const double slope_budding_point = topography.slope_between_points( lobe_parent.center, final_budding_point ); // compute the new lobe axes MrLavaLoba::compute_lobe_axes( lobe_cur, slope_budding_point, input, common_lobe_dimensions ); @@ -549,6 +552,7 @@ RunStatus Simulation::steps( int n_steps ) { // Increment idx_flow idx_flow++; + simulation_state->beginning_of_new_flow = true; if( input.save_hazard_data ) { @@ -563,8 +567,8 @@ RunStatus Simulation::steps( int n_steps ) if( input.print_remaining_time ) { - auto t_cur = std::chrono::high_resolution_clock::now(); - auto remaining_time = std::chrono::duration_cast( + const auto t_cur = std::chrono::high_resolution_clock::now(); + 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 ); } @@ -592,7 +596,7 @@ void Simulation::run() int n_lobes_processed = 0; auto t_run_start = std::chrono::high_resolution_clock::now(); - auto common_lobe_dimensions = CommonLobeDimensions( input ); + const auto common_lobe_dimensions = CommonLobeDimensions( input ); for( int idx_flow = 0; idx_flow < input.n_flows; idx_flow++ ) { @@ -617,7 +621,7 @@ void Simulation::run() lobe_cur.thickness = MrLavaLoba::compute_current_lobe_thickness( idx_lobe, n_lobes, input, common_lobe_dimensions ); - auto [height_lobe_center, slope] = topography.height_and_slope( lobe_cur.center ); + const auto [height_lobe_center, slope] = topography.height_and_slope( lobe_cur.center ); if( height_lobe_center == topography.no_data_value ) { @@ -648,8 +652,8 @@ void Simulation::run() // Select which of the previously created lobes is the parent lobe // from which the new descendent lobe will bud - auto idx_parent = MrLavaLoba::select_parent_lobe( idx_lobe, lobes, input, common_lobe_dimensions, gen ); - Lobe & lobe_parent = lobes[idx_parent]; + const auto idx_parent = MrLavaLoba::select_parent_lobe( idx_lobe, lobes, input, common_lobe_dimensions, gen ); + const Lobe & lobe_parent = lobes[idx_parent]; // 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] ) ) @@ -662,10 +666,10 @@ void Simulation::run() // points on the ellipse) Flowy::Vector2 budding_point = topography.find_preliminary_budding_point( lobe_parent, input.npoints ); - auto [height_lobe_center, slope_parent] = topography.height_and_slope( lobe_parent.center ); - auto [height_bp, slope_bp] = topography.height_and_slope( budding_point ); + const auto [height_lobe_center, slope_parent] = topography.height_and_slope( lobe_parent.center ); + const auto [height_bp, slope_bp] = topography.height_and_slope( budding_point ); - Vector2 diff = ( budding_point - lobe_parent.center ); + const Vector2 diff = ( budding_point - lobe_parent.center ); // 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 @@ -677,8 +681,8 @@ void Simulation::run() // 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 ); + 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 ); // final_budding_point = budding_point; if( stop_condition( final_budding_point, lobe_parent.semi_axes[0] ) ) @@ -687,7 +691,7 @@ void Simulation::run() break; } // Get the slope at the final budding point - double slope_budding_point = topography.slope_between_points( lobe_parent.center, final_budding_point ); + const double slope_budding_point = topography.slope_between_points( lobe_parent.center, final_budding_point ); // compute the new lobe axes MrLavaLoba::compute_lobe_axes( lobe_cur, slope_budding_point, input, common_lobe_dimensions ); @@ -724,22 +728,22 @@ void Simulation::run() if( input.print_remaining_time ) { - auto t_cur = std::chrono::high_resolution_clock::now(); - auto remaining_time = std::chrono::duration_cast( + const auto t_cur = std::chrono::high_resolution_clock::now(); + const 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 ); } } - auto t_cur = std::chrono::high_resolution_clock::now(); - auto total_time = std::chrono::duration_cast( ( t_cur - t_run_start ) ); + const auto t_cur = std::chrono::high_resolution_clock::now(); + const 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 ) { - auto lobes_per_ms = n_lobes_processed / total_time.count(); + const auto lobes_per_ms = n_lobes_processed / total_time.count(); fmt::print( "n_lobes/ms = {}\n", lobes_per_ms ); } From 2e99a359eecc2c2db2b4bead32e0d2b46045cd70 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Wed, 26 Mar 2025 22:15:15 +0000 Subject: [PATCH 08/15] WIP: Simulation::steps function Co-authored-by: Amrita Goswami --- flowy/include/models/mr_lava_loba.hpp | 2 +- flowy/include/simulation.hpp | 10 ++-- src/simulation.cpp | 73 ++++++++++++++++----------- 3 files changed, 49 insertions(+), 36 deletions(-) diff --git a/flowy/include/models/mr_lava_loba.hpp b/flowy/include/models/mr_lava_loba.hpp index 979f908..d21d022 100644 --- a/flowy/include/models/mr_lava_loba.hpp +++ b/flowy/include/models/mr_lava_loba.hpp @@ -75,7 +75,7 @@ inline int compute_n_lobes( int idx_flow, const Config::InputParams & input, std 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 + // 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 ); diff --git a/flowy/include/simulation.hpp b/flowy/include/simulation.hpp index 5b6fcc2..7ae7412 100644 --- a/flowy/include/simulation.hpp +++ b/flowy/include/simulation.hpp @@ -26,12 +26,10 @@ struct SimulationState std::chrono::time_point t_run_start{}; std::vector lobes{}; - int n_lobes_processed = 0; - int n_lobes = 0; - int step = 0; - int idx_flow = 0; - int idx_lobe = 0; - int n_lobes_current_flow = 0; + int n_lobes_processed = 0; + int n_lobes = 0; + int idx_flow = 0; + int idx_lobe = 0; bool beginning_of_new_flow = true; }; diff --git a/src/simulation.cpp b/src/simulation.cpp index ea0de38..19f717c 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -231,7 +231,8 @@ void Simulation::write_avg_thickness_file() // This lambda performs bisection search to find the threshold thickness at which a // relative volume proportion of `thresh` is contained within cells with greater thickness than the threshold thickness - auto bisection_search = [&]( double thresh, double tol, int max_iter ) { + auto bisection_search = [&]( double thresh, double tol, int max_iter ) + { int idx_lo = 0; int idx_hi = n_cells - 1; @@ -417,32 +418,43 @@ RunStatus Simulation::steps( int n_steps ) auto & n_lobes = simulation_state->n_lobes; // Now we have to figure out the value of idx_flow and if we are past n_init or not - const int step_initial = simulation_state->step; - - for( int step = step_initial; step < step_initial + n_steps; step++ ) + for( int step = 0; step < n_steps; step++ ) { - const bool is_last_step_of_flow = idx_lobe == n_lobes - 1; - const bool is_first_step_of_new_flow = simulation_state->beginning_of_new_flow; - const bool is_an_initial_lobe = idx_lobe < input.n_init; - const bool is_last_step = ( idx_flow == input.n_flows ) && is_last_step_of_flow; + std::cout << "step " << step << "\n"; + std::cout << "idx_flow " << idx_flow << "\n"; + std::cout << "idx_lobe " << idx_lobe << "\n"; + std::cout << "n_lobes " << n_lobes << "\n"; + std::cout << "n_lobes_processed " << n_lobes_processed << "\n"; + + const 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( is_first_step_of_new_flow ) + std::cout << "is_last_step_of_flow " << is_last_step_of_flow << "\n"; + std::cout << "is_an_initial_lobe " << is_an_initial_lobe << "\n"; + std::cout << "is_last_step " << is_last_step << "\n ===== \n\n"; + + if( simulation_state->beginning_of_new_flow ) { - // simulation_state->idx_flow++; - simulation_state->n_lobes_current_flow - = MrLavaLoba::compute_n_lobes( simulation_state->idx_flow, input, gen ); - simulation_state->lobes = std::vector{}; - simulation_state->lobes.reserve( simulation_state->n_lobes_current_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_current_flow ); + 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 + // switch back the flag simulation_state->beginning_of_new_flow = false; } + if( idx_lobe >= n_lobes || idx_flow >= input.n_flows ) + { + run_status = RunStatus::Finished; + break; + } + if( is_an_initial_lobe ) { lobes.emplace_back(); @@ -472,7 +484,6 @@ RunStatus Simulation::steps( int n_steps ) // Add rasterized lobe topography.add_lobe( lobe_cur, input.volume_correction, idx_lobe ); - n_lobes_processed++; write_thickness_if_necessary( n_lobes_processed ); } else @@ -482,7 +493,8 @@ RunStatus Simulation::steps( int n_steps ) // 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 auto idx_parent + = MrLavaLoba::select_parent_lobe( idx_lobe, lobes, input, common_lobe_dimensions, gen ); const Lobe & lobe_parent = lobes[idx_parent]; // stopping condition (parent lobe close the domain boundary or at a not defined z value) @@ -495,7 +507,8 @@ RunStatus Simulation::steps( int n_steps ) // 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 ); + const Flowy::Vector2 budding_point + = topography.find_preliminary_budding_point( lobe_parent, input.npoints ); const auto [height_lobe_center, slope_parent] = topography.height_and_slope( lobe_parent.center ); const auto [height_bp, slope_bp] = topography.height_and_slope( budding_point ); @@ -523,7 +536,8 @@ RunStatus Simulation::steps( int n_steps ) break; } // Get the slope at the final budding point - const double slope_budding_point = topography.slope_between_points( lobe_parent.center, final_budding_point ); + const double slope_budding_point + = topography.slope_between_points( lobe_parent.center, final_budding_point ); // compute the new lobe axes MrLavaLoba::compute_lobe_axes( lobe_cur, slope_budding_point, input, common_lobe_dimensions ); @@ -544,14 +558,14 @@ RunStatus Simulation::steps( int n_steps ) // Add rasterized lobe topography.add_lobe( lobe_cur, input.volume_correction, idx_lobe ); - n_lobes_processed++; write_thickness_if_necessary( n_lobes_processed ); } + n_lobes_processed++; + idx_lobe++; + if( is_last_step_of_flow ) { - // Increment idx_flow - idx_flow++; simulation_state->beginning_of_new_flow = true; if( input.save_hazard_data ) @@ -572,15 +586,14 @@ RunStatus Simulation::steps( int n_steps ) ( 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 ); } + idx_flow++; + idx_lobe = 0; } if( is_last_step ) { run_status = RunStatus::Finished; } - - n_lobes_processed++; - idx_lobe++; } if( run_status == RunStatus::Finished ) @@ -595,7 +608,7 @@ void Simulation::run() { int n_lobes_processed = 0; - auto t_run_start = std::chrono::high_resolution_clock::now(); + auto t_run_start = std::chrono::high_resolution_clock::now(); const auto common_lobe_dimensions = CommonLobeDimensions( input ); for( int idx_flow = 0; idx_flow < input.n_flows; idx_flow++ ) @@ -652,7 +665,8 @@ void Simulation::run() // 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 auto idx_parent + = MrLavaLoba::select_parent_lobe( idx_lobe, lobes, input, common_lobe_dimensions, gen ); const Lobe & lobe_parent = lobes[idx_parent]; // stopping condition (parent lobe close the domain boundary or at a not defined z value) @@ -691,7 +705,8 @@ void Simulation::run() break; } // Get the slope at the final budding point - const double slope_budding_point = topography.slope_between_points( lobe_parent.center, final_budding_point ); + const double slope_budding_point + = topography.slope_between_points( lobe_parent.center, final_budding_point ); // compute the new lobe axes MrLavaLoba::compute_lobe_axes( lobe_cur, slope_budding_point, input, common_lobe_dimensions ); From 469d0124ef0570cd6afe7b8c85be4dd4bbf72ad4 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Wed, 26 Mar 2025 23:40:27 +0000 Subject: [PATCH 09/15] ENH: simulation: removed some duplication Co-authored-by: Amrita Goswami --- src/simulation.cpp | 39 ++++++++++++--------------------------- 1 file changed, 12 insertions(+), 27 deletions(-) diff --git a/src/simulation.cpp b/src/simulation.cpp index 19f717c..e6f76c1 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -231,8 +231,7 @@ void Simulation::write_avg_thickness_file() // This lambda performs bisection search to find the threshold thickness at which a // relative volume proportion of `thresh` is contained within cells with greater thickness than the threshold thickness - auto bisection_search = [&]( double thresh, double tol, int max_iter ) - { + auto bisection_search = [&]( double thresh, double tol, int max_iter ) { int idx_lo = 0; int idx_hi = n_cells - 1; @@ -455,17 +454,17 @@ RunStatus Simulation::steps( int n_steps ) 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 ) { - lobes.emplace_back(); - Lobe & lobe_cur = lobes.back(); - MrLavaLoba::compute_initial_lobe_position( idx_flow, lobe_cur, input, gen ); - // Compute the thickness of the lobe - lobe_cur.thickness - = MrLavaLoba::compute_current_lobe_thickness( idx_lobe, n_lobes, input, common_lobe_dimensions ); - const auto [height_lobe_center, slope] = topography.height_and_slope( lobe_cur.center ); if( height_lobe_center == topography.no_data_value ) @@ -481,16 +480,9 @@ RunStatus Simulation::steps( int n_steps ) // compute lobe axes MrLavaLoba::compute_lobe_axes( lobe_cur, slope_norm, input, common_lobe_dimensions ); - - // Add rasterized lobe - topography.add_lobe( lobe_cur, input.volume_correction, idx_lobe ); - write_thickness_if_necessary( n_lobes_processed ); } else { - lobes.emplace_back(); - Lobe & lobe_cur = lobes.back(); - // Select which of the previously created lobes is the parent lobe // from which the new descendent lobe will bud const auto idx_parent @@ -510,8 +502,7 @@ RunStatus Simulation::steps( int n_steps ) const Flowy::Vector2 budding_point = topography.find_preliminary_budding_point( lobe_parent, input.npoints ); - const auto [height_lobe_center, slope_parent] = topography.height_and_slope( lobe_parent.center ); - const auto [height_bp, slope_bp] = topography.height_and_slope( budding_point ); + const auto [height_parent, slope_parent] = topography.height_and_slope( lobe_parent.center ); const Vector2 diff = ( budding_point - lobe_parent.center ); @@ -551,16 +542,11 @@ RunStatus Simulation::steps( int n_steps ) run_status = RunStatus::Finished; break; } - - // Compute the thickness of the lobe - lobe_cur.thickness - = MrLavaLoba::compute_current_lobe_thickness( idx_lobe, n_lobes, input, common_lobe_dimensions ); - - // Add rasterized lobe - topography.add_lobe( lobe_cur, input.volume_correction, idx_lobe ); - write_thickness_if_necessary( n_lobes_processed ); } + // Add rasterized lobe + topography.add_lobe( lobe_cur, input.volume_correction, idx_lobe ); + write_thickness_if_necessary( n_lobes_processed ); n_lobes_processed++; idx_lobe++; @@ -681,7 +667,6 @@ void Simulation::run() Flowy::Vector2 budding_point = topography.find_preliminary_budding_point( lobe_parent, input.npoints ); const auto [height_lobe_center, slope_parent] = topography.height_and_slope( lobe_parent.center ); - const auto [height_bp, slope_bp] = topography.height_and_slope( budding_point ); const Vector2 diff = ( budding_point - lobe_parent.center ); From d52cb63058a78065629e6a02f2f336f51e27d885 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 27 Mar 2025 00:35:48 +0000 Subject: [PATCH 10/15] ENH: refactored simulation::steps function moved some functional units into new class methods Co-authored-by: Amrita Goswami --- flowy/include/simulation.hpp | 17 +++ src/simulation.cpp | 235 +++++++++++++++++++---------------- 2 files changed, 144 insertions(+), 108 deletions(-) diff --git a/flowy/include/simulation.hpp b/flowy/include/simulation.hpp index 7ae7412..6821a5b 100644 --- a/flowy/include/simulation.hpp +++ b/flowy/include/simulation.hpp @@ -5,6 +5,7 @@ #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 @@ -34,6 +35,12 @@ struct SimulationState bool beginning_of_new_flow = true; }; +enum class FlowStatus +{ + Finished, + Ongoing +}; + enum class RunStatus { Finished, @@ -79,10 +86,20 @@ class Simulation // Perform `n_steps` steps of the simulation (per step, a single lobe is added to the topography) RunStatus steps( int n_steps ); + std::optional get_simulation_state() + { + 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/src/simulation.cpp b/src/simulation.cpp index e6f76c1..20d5ab2 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" @@ -397,6 +398,112 @@ void Simulation::save_post_run_output() write_avg_thickness_file(); } +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 ); + + const auto [height_lobe_center, slope] = topography.height_and_slope( lobe_cur.center ); + + 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." ); + } + + // 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 ); + + // compute lobe axes + MrLavaLoba::compute_lobe_axes( lobe_cur, slope_norm, input, common_lobe_dimensions ); +} + +[[nodiscard]] FlowStatus Simulation::process_descendent_lobe( + int idx_lobe, std::vector & lobes, const CommonLobeDimensions & common_lobe_dimensions ) +{ + Lobe & lobe_cur = lobes.back(); + + // 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]; + + // 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; + } + + // 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 ); + + const auto [height_parent, slope_parent] = topography.height_and_slope( lobe_parent.center ); + + const Vector2 diff = ( budding_point - lobe_parent.center ); + + // 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 ); + + // Add the inertial contribution + MrLavaLoba::add_inertial_contribution( lobe_cur, lobe_parent, slope_parent_norm, 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 + 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 ); + + // 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 ); + + // compute the new lobe axes + MrLavaLoba::compute_lobe_axes( lobe_cur, slope_budding_point, input, common_lobe_dimensions ); + + // Get new lobe center + MrLavaLoba::compute_descendent_lobe_position( lobe_cur, lobe_parent, final_budding_point, input ); + + if( stop_condition( lobe_cur.center, lobe_cur.semi_axes[0] ) ) + { + lobes.pop_back(); + return FlowStatus::Finished; + } + + return FlowStatus::Ongoing; +} + +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 ); + } + + if( input.write_lobes_csv ) + { + write_lobe_data_to_file( lobes, input.output_folder / fmt::format( "lobes_{}.csv", idx_flow ) ); + } + + if( input.print_remaining_time ) + { + const auto t_cur = std::chrono::high_resolution_clock::now(); + 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; @@ -416,22 +523,14 @@ RunStatus Simulation::steps( int n_steps ) 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 - for( int step = 0; step < n_steps; step++ ) + while( n_lobes_processed < n_lobes_processed_initial + n_steps ) { - std::cout << "step " << step << "\n"; - std::cout << "idx_flow " << idx_flow << "\n"; - std::cout << "idx_lobe " << idx_lobe << "\n"; - std::cout << "n_lobes " << n_lobes << "\n"; - std::cout << "n_lobes_processed " << n_lobes_processed << "\n"; - - const 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; - - std::cout << "is_last_step_of_flow " << is_last_step_of_flow << "\n"; - std::cout << "is_an_initial_lobe " << is_an_initial_lobe << "\n"; - std::cout << "is_last_step " << is_last_step << "\n ===== \n\n"; + 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 ) { @@ -463,115 +562,35 @@ RunStatus Simulation::steps( int n_steps ) if( is_an_initial_lobe ) { - MrLavaLoba::compute_initial_lobe_position( idx_flow, lobe_cur, input, gen ); + process_initial_lobe( idx_flow, lobe_cur, common_lobe_dimensions ); - const auto [height_lobe_center, slope] = topography.height_and_slope( lobe_cur.center ); - - 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." ); - } - - // 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 ); - - // compute lobe axes - MrLavaLoba::compute_lobe_axes( lobe_cur, slope_norm, input, 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 ); } else { - // 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]; - - // 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(); - run_status = RunStatus::Finished; - break; - } - - // 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 ); - - const auto [height_parent, slope_parent] = topography.height_and_slope( lobe_parent.center ); - - const Vector2 diff = ( budding_point - lobe_parent.center ); - - // 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 ); - - // Add the inertial contribution - MrLavaLoba::add_inertial_contribution( lobe_cur, lobe_parent, slope_parent_norm, 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 - 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 ); - - // final_budding_point = budding_point; - if( stop_condition( final_budding_point, lobe_parent.semi_axes[0] ) ) + const auto flow_status = process_descendent_lobe( idx_lobe, lobes, common_lobe_dimensions ); + if( flow_status == FlowStatus::Ongoing ) { - lobes.pop_back(); - run_status = RunStatus::Finished; - break; + // 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 ); } - // Get the slope at the final budding point - const double slope_budding_point - = topography.slope_between_points( lobe_parent.center, final_budding_point ); - - // compute the new lobe axes - MrLavaLoba::compute_lobe_axes( lobe_cur, slope_budding_point, input, common_lobe_dimensions ); - - // Get new lobe center - MrLavaLoba::compute_descendent_lobe_position( lobe_cur, lobe_parent, final_budding_point, input ); - - if( stop_condition( lobe_cur.center, lobe_cur.semi_axes[0] ) ) + else { - lobes.pop_back(); - run_status = RunStatus::Finished; - break; + is_last_step_of_flow = true; } } - // Add rasterized lobe - topography.add_lobe( lobe_cur, input.volume_correction, idx_lobe ); - write_thickness_if_necessary( n_lobes_processed ); - n_lobes_processed++; idx_lobe++; if( is_last_step_of_flow ) { + post_flow_hook( idx_flow, lobes ); simulation_state->beginning_of_new_flow = true; - - if( input.save_hazard_data ) - { - compute_cumulative_descendents( lobes ); - topography.compute_hazard_flow( lobes ); - } - - if( input.write_lobes_csv ) - { - write_lobe_data_to_file( lobes, input.output_folder / fmt::format( "lobes_{}.csv", idx_flow ) ); - } - - if( input.print_remaining_time ) - { - const auto t_cur = std::chrono::high_resolution_clock::now(); - 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 ); - } idx_flow++; idx_lobe = 0; } From 8407b8a5c44bd286836d0f7a234f240cede6aa80 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 27 Mar 2025 00:46:26 +0000 Subject: [PATCH 11/15] simulation: new function reset_simulation_state simulation: implement `run` in terms of `steps` --- flowy/include/simulation.hpp | 7 +- src/simulation.cpp | 176 ++++------------------------------- 2 files changed, 25 insertions(+), 158 deletions(-) diff --git a/flowy/include/simulation.hpp b/flowy/include/simulation.hpp index 6821a5b..2c823ae 100644 --- a/flowy/include/simulation.hpp +++ b/flowy/include/simulation.hpp @@ -86,7 +86,12 @@ class Simulation // Perform `n_steps` steps of the simulation (per step, a single lobe is added to the topography) RunStatus steps( int n_steps ); - std::optional get_simulation_state() + void reset_simulation_state() + { + simulation_state = std::nullopt; + } + + std::optional get_simulation_state() const { return simulation_state; } diff --git a/src/simulation.cpp b/src/simulation.cpp index 20d5ab2..2b6496e 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -603,172 +603,34 @@ RunStatus Simulation::steps( int n_steps ) if( run_status == RunStatus::Finished ) { - save_post_run_output(); - } - - return run_status; -} -void Simulation::run() -{ - int n_lobes_processed = 0; + 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 ); - auto t_run_start = std::chrono::high_resolution_clock::now(); - const auto common_lobe_dimensions = CommonLobeDimensions( input ); - - for( int idx_flow = 0; idx_flow < input.n_flows; idx_flow++ ) - { - // Determine n_lobes - int n_lobes = MrLavaLoba::compute_n_lobes( idx_flow, input, gen ); + fmt::print( "Total number of processed lobes = {}\n", n_lobes_processed ); - lobes = std::vector{}; - lobes.reserve( n_lobes ); - - // set the intersection cache - topography.reset_intersection_cache( n_lobes ); - - // Build initial lobes which do not propagate descendents - for( int idx_lobe = 0; idx_lobe < input.n_init; idx_lobe++ ) + if( total_time.count() > 0 ) { - lobes.emplace_back(); - Lobe & lobe_cur = lobes.back(); - - MrLavaLoba::compute_initial_lobe_position( idx_flow, lobe_cur, input, gen ); - - // Compute the thickness of the lobe - lobe_cur.thickness - = MrLavaLoba::compute_current_lobe_thickness( idx_lobe, n_lobes, input, common_lobe_dimensions ); - - const auto [height_lobe_center, slope] = topography.height_and_slope( lobe_cur.center ); - - 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." ); - } - - // 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 ); - - // compute lobe axes - MrLavaLoba::compute_lobe_axes( lobe_cur, slope_norm, input, common_lobe_dimensions ); - - // Add rasterized lobe - topography.add_lobe( lobe_cur, input.volume_correction, idx_lobe ); - n_lobes_processed++; - write_thickness_if_necessary( n_lobes_processed ); + const auto lobes_per_ms = n_lobes_processed / total_time.count(); + fmt::print( "n_lobes/ms = {}\n", lobes_per_ms ); } - // 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(); - - // 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]; - - // 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; - } - - // 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 ); - - const auto [height_lobe_center, slope_parent] = topography.height_and_slope( lobe_parent.center ); - - const Vector2 diff = ( budding_point - lobe_parent.center ); - - // 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 ); - - // Add the inertial contribution - MrLavaLoba::add_inertial_contribution( lobe_cur, lobe_parent, slope_parent_norm, 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 - 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 ); - - // 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 - const double slope_budding_point - = topography.slope_between_points( lobe_parent.center, final_budding_point ); - - // compute the new lobe axes - MrLavaLoba::compute_lobe_axes( lobe_cur, slope_budding_point, input, common_lobe_dimensions ); - - // Get new lobe center - MrLavaLoba::compute_descendent_lobe_position( lobe_cur, lobe_parent, final_budding_point, input ); - - if( stop_condition( lobe_cur.center, lobe_cur.semi_axes[0] ) ) - { - lobes.pop_back(); - break; - } - - // Compute the thickness of the lobe - lobe_cur.thickness - = MrLavaLoba::compute_current_lobe_thickness( idx_lobe, n_lobes, input, common_lobe_dimensions ); - - // Add rasterized lobe - 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 ) - { - compute_cumulative_descendents( lobes ); - topography.compute_hazard_flow( lobes ); - } - - if( input.write_lobes_csv ) - { - write_lobe_data_to_file( lobes, input.output_folder / fmt::format( "lobes_{}.csv", idx_flow ) ); - } - - if( input.print_remaining_time ) - { - const auto t_cur = std::chrono::high_resolution_clock::now(); - const 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 ); - } + fmt::print( "Used RNG seed: {}\n", rng_seed ); + save_post_run_output(); } - const auto t_cur = std::chrono::high_resolution_clock::now(); - const 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 ); + return run_status; +} - if( total_time.count() > 0 ) +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 ) { - const auto lobes_per_ms = n_lobes_processed / total_time.count(); - fmt::print( "n_lobes/ms = {}\n", lobes_per_ms ); - } - - fmt::print( "Used RNG seed: {}\n", rng_seed ); - - save_post_run_output(); + }; } - } // namespace Flowy From ba654c77773da309f63b874537c296c847732133 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 27 Mar 2025 21:12:25 +0000 Subject: [PATCH 12/15] BLD: save `with_netcdf` option as a variable In pyflowy we can then fetch it from the flowy subproject with get_variable. Co-authored-by: Amrita Goswami --- meson.build | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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) From 5c35ddab0d5be68b2fa8afcd28c58f9d19301d8e Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 27 Mar 2025 22:08:40 +0000 Subject: [PATCH 13/15] FIX: no auto types for chrono timepoints Got some errors with clang... Co-authored-by: Amrita Goswami --- flowy/include/simulation.hpp | 2 +- src/simulation.cpp | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/flowy/include/simulation.hpp b/flowy/include/simulation.hpp index 2c823ae..ee8d7fc 100644 --- a/flowy/include/simulation.hpp +++ b/flowy/include/simulation.hpp @@ -24,7 +24,7 @@ namespace Flowy */ struct SimulationState { - std::chrono::time_point t_run_start{}; + std::chrono::time_point t_run_start{}; std::vector lobes{}; int n_lobes_processed = 0; diff --git a/src/simulation.cpp b/src/simulation.cpp index 2b6496e..04e4844 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -345,7 +346,7 @@ Simulation::get_file_handle( const Topography & topography, OutputQuantity outpu void Simulation::compute_topography_thickness() { - // Compute the thickness by subtracting 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; @@ -497,7 +498,9 @@ void Simulation::post_flow_hook( int idx_flow, std::vector & lobes ) if( input.print_remaining_time ) { - const auto t_cur = std::chrono::high_resolution_clock::now(); + const std::chrono::time_point t_cur + = std::chrono::high_resolution_clock::now(); + 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 ); From 3a0db72a541d3daabebd3e30f8b81500d2e7f549 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 27 Mar 2025 22:18:11 +0000 Subject: [PATCH 14/15] FIX: minor changes to make precommit happy Co-authored-by: Amrita Goswami --- src/simulation.cpp | 4 ++-- src/topography.cpp | 7 ++++--- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/simulation.cpp b/src/simulation.cpp index 04e4844..ddc9b65 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -233,7 +233,8 @@ void Simulation::write_avg_thickness_file() // This lambda performs bisection search to find the threshold thickness at which a // relative volume proportion of `thresh` is contained within cells with greater thickness than the threshold thickness - auto bisection_search = [&]( double thresh, double tol, int max_iter ) { + auto bisection_search = [&]( double thresh, double tol, int max_iter ) + { int idx_lo = 0; int idx_hi = n_cells - 1; @@ -606,7 +607,6 @@ RunStatus Simulation::steps( int n_steps ) if( run_status == RunStatus::Finished ) { - 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 ) ); 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( From afdb90aa743e0d3f994242f2ff478760b255fd23 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 27 Mar 2025 22:26:44 +0000 Subject: [PATCH 15/15] FIX: pre_commit update to setup_micromamba action provision-with-micromamba has been deprecated --- .github/workflows/pre_commit.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) 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: |