ENH: simulation: wrap up implementation of steps

ENH: simulation: better const correctness :)

Co-authored-by: Amrita Goswami <amrita16thaug646@gmail.com>

Author: MSallermann

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<Lobe> & lobes, const Config::InputParams & input,
-    CommonLobeDimensions & lobe_dimensions, std::mt19937 & gen )
+    const CommonLobeDimensions & lobe_dimensions, std::mt19937 & gen )
 {
     Lobe & lobe_descendent = lobes[idx_descendant];
 

flowy/include/simulation.hpp

@@ -23,15 +23,17 @@ namespace Flowy
  */
 struct SimulationState
 {
-    int n_lobes_processed = 0;
     std::chrono::time_point<std::chrono::system_clock> t_run_start{};
-    int n_lobes{};
     std::vector<Lobe> 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

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<std::chrono::milliseconds>(
+                const auto t_cur          = std::chrono::high_resolution_clock::now();
+                const auto remaining_time = std::chrono::duration_cast<std::chrono::milliseconds>(
                     ( 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<std::chrono::milliseconds>(
+            const auto t_cur          = std::chrono::high_resolution_clock::now();
+            const auto remaining_time = std::chrono::duration_cast<std::chrono::milliseconds>(
                 ( 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<std::chrono::milliseconds>( ( t_cur - t_run_start ) );
+    const auto t_cur      = std::chrono::high_resolution_clock::now();
+    const auto total_time = std::chrono::duration_cast<std::chrono::milliseconds>( ( 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 );
     }