Merge pull request #36 from flowy-code/intermediate_dems

Intermediate dems

Author: MSallermann

flowy/include/asc_file.hpp

@@ -12,7 +12,7 @@ class AscFile : public TopographyFile
 {
 public:
     AscFile() = default;
-    AscFile( const Topography & topography, OutputQuantitiy output )
+    AscFile( const Topography & topography, OutputQuantity output )
             : TopographyFile::TopographyFile( topography, output )
     {
     }

flowy/include/config.hpp

@@ -39,6 +39,9 @@ class InputParams
     bool print_remaining_time = false;
     bool save_final_dem       = false;
 
+    // If this std::optional has a value `n`, write the lava thickness (not the dem!) to a file every `n` lobes
+    std::optional<int> write_thickness_every_n_lobes{};
+
     // The tolerance in the volume ratio when finding the threshold thickness for masking
     double masking_tolerance = 1e-5;
     // The maximum number of bisection search iterations when finding the threshold thickness for masking

flowy/include/lobe.hpp

@@ -130,27 +130,29 @@ class Lobe
         const double a2 = a * a;
         const double b2 = b * b;
 
+        // Compute coefficients of quadratic equation
         const double alpha = 1.0 / ( a2 ) * ( diff[0] * diff[0] ) + 1.0 / b2 * ( diff[1] * diff[1] );
         const double beta  = 2.0 * ( x1_prime[0] * diff[0] / a2 + x1_prime[1] * diff[1] / b2 );
         const double gamma = x1_prime[0] * x1_prime[0] / a2 + x1_prime[1] * x1_prime[1] / b2 - 1.0;
 
-        // The solution to this quadratic equation is
+        // The solution to the quadratic equation is
         // t = (-beta +- sqrt(beta^2 - 4*alpha*gamma)) / (2*alpha)
-
-        // Therefore, if beta^2 - 4*alpha*gamma < 0, the line se\beta\frgment misses the ellipse
+        // Therefore, if beta^2 - 4*alpha*gamma < 0, the line segment misses the ellipse
         const double radicand = beta * beta - 4 * alpha * gamma;
         if( radicand < 0 )
             return std::nullopt;
 
         const double sqrt_r = std::sqrt( radicand );
 
-        // Else, we compute the two points of intersection and check if they fall into the interval [0, 1]
+        // Else (if the line segment does not miss the ellipse) we compute the two points of intersection t1 and t2 and check if
         const double t1 = ( -beta - sqrt_r ) / ( 2.0 * alpha );
         const double t2 = ( -beta + sqrt_r ) / ( 2.0 * alpha );
 
-        // This condition determines if any of the points of the line segment lie within the lobe
-        // This is the condition for the intersection of [0,1] and [t1, t2] to not be empty
+        // The following condition determines if any of the points of the line segment lie within the lobe
+        // It can be thought of as the the condition for the intersection of [0,1] and [t1,t2] to not be empty
         const bool condition = !( ( t1 < 0 && t2 < 0 ) || ( t1 > 1 && t2 > 1 ) );
+
+        // If the conditions if fulfilled, we compute the intersection points
         if( condition )
         {
             // We clamp t1 and t2 to the interval [0,1] in case the line segment intersects the ellipse only once
@@ -165,20 +167,22 @@ class Lobe
             const std::array<Vector2, 2> res = { v1 + center, v2 + center };
             return res;
         }
-
-        return std::nullopt;
+        else // else we return nullopt
+        {
+            return std::nullopt;
+        }
     }
 
     // Gives a point on the perimeter of the ellipse
     // The angle is relative to the semi major axis angle
-    inline Vector2 point_at_angle( const double phi ) const
+    inline Vector2 point_at_angle( double phi ) const
     {
         return point_at_angle( std::sin( phi ), std::cos( phi ) );
     }
 
     // Gives a point on the perimeter of the ellipse
     // The angle is relative to the semi major axis angle
-    inline Vector2 point_at_angle( const double sin_phi, const double cos_phi ) const
+    inline Vector2 point_at_angle( double sin_phi, double cos_phi ) const
     {
         const double a      = semi_axes[0]; // major axis
         const double b      = semi_axes[1]; // minor axis
@@ -187,20 +191,8 @@ class Lobe
         return coord + center;
     }
 
-    inline std::vector<Vector2> rasterize_perimeter( int n_raster_points ) const
-    {
-        VectorX phi_list = xt::linspace<double>( 0.0, 2.0 * Math::pi, n_raster_points, false );
-        auto res         = std::vector<Vector2>( n_raster_points );
-
-        for( int idx_phi = 0; idx_phi < n_raster_points; idx_phi++ )
-        {
-            res[idx_phi] = point_at_angle( phi_list[idx_phi] );
-        }
-
-        return res;
-    }
-
-    inline std::vector<Vector2> rasterize_perimeter( std::span<double> sin_phi, std::span<double> cos_phi ) const
+    inline std::vector<Vector2>
+    rasterize_perimeter( const std::span<double> sin_phi, const std::span<double> cos_phi ) const
     {
         const int n_raster_points = sin_phi.size();
         auto res                  = std::vector<Vector2>( n_raster_points );

flowy/include/models/mr_lava_loba.hpp

@@ -80,7 +80,7 @@ class MrLavaLoba
     }
 
     // Calculate n_lobes
-    int compute_n_lobes( int idx_flow )
+    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
@@ -115,7 +115,7 @@ class MrLavaLoba
     }
 
     // calculates the initial lobe position
-    void compute_initial_lobe_position( int idx_flow, Lobe & lobe )
+    void compute_initial_lobe_position( int idx_flow, Lobe & lobe ) const
     {
         std::unique_ptr<VentFlag> f{};
 
@@ -182,15 +182,15 @@ class MrLavaLoba
         lobe.semi_axes = { semi_major_axis, semi_minor_axis };
     }
 
-    void compute_descendent_lobe_position( Lobe & lobe, const Lobe & parent, Vector2 final_budding_point )
+    void compute_descendent_lobe_position( Lobe & lobe, const Lobe & parent, const Vector2 & final_budding_point ) const
     {
         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;
     }
 
-    void perturb_lobe_angle( Lobe & lobe, double slope )
+    void perturb_lobe_angle( Lobe & lobe, double slope ) const
     {
         const double slope_deg = 180.0 / Math::pi * std::atan( slope );
 
@@ -216,7 +216,7 @@ class MrLavaLoba
     }
 
     // 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<Lobe> & lobes )
+    int select_parent_lobe( int idx_descendant, std::vector<Lobe> & lobes ) const
     {
         Lobe & lobe_descendent = lobes[idx_descendant];
 

flowy/include/netcdf_file.hpp

@@ -24,7 +24,7 @@ class NetCDFFile : public TopographyFile
 
     StorageDataType data_type = StorageDataType::Short;
     NetCDFFile()              = default;
-    NetCDFFile( const Topography & topography, OutputQuantitiy output )
+    NetCDFFile( const Topography & topography, OutputQuantity output )
             : TopographyFile::TopographyFile( topography, output )
     {
     }

flowy/include/simulation.hpp

@@ -39,8 +39,14 @@ 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);
+
+    // Computes the topography_thickness field by substacting the initial topography and dividing by (1.0 - filling_parameter)
+    void compute_topography_thickness();
+
     std::unique_ptr<TopographyFile>
-    get_file_handle( const Topography & topography, OutputQuantitiy output_quantity ) const;
+    get_file_handle( const Topography & topography, OutputQuantity output_quantity ) const;
 
     void run();
 

flowy/include/topography_file.hpp

@@ -19,7 +19,7 @@ struct TopographyCrop
     double y_max;
 };
 
-enum class OutputQuantitiy
+enum class OutputQuantity
 {
     Hazard,
     Height
@@ -51,14 +51,14 @@ class TopographyFile
     TopographyFile() = default;
 
     // Constructor that takes topography
-    TopographyFile( const Topography & topography, OutputQuantitiy output )
+    TopographyFile( const Topography & topography, OutputQuantity output )
             : x_data( topography.x_data ), y_data( topography.y_data ), no_data_value( topography.no_data_value )
     {
-        if( output == OutputQuantitiy::Height )
+        if( output == OutputQuantity::Height )
         {
             data = topography.height_data;
         }
-        else if( output == OutputQuantitiy::Hazard )
+        else if( output == OutputQuantity::Hazard )
         {
             data = topography.hazard;
         }

src/config_parser.cpp

@@ -52,6 +52,7 @@ InputParams parse_config( const std::filesystem::path & path )
     set_if_specified( params.write_lobes_csv, tbl["write_lobes_csv"] );
     set_if_specified( params.print_remaining_time, tbl["print_remaining_time"] );
     set_if_specified( params.save_final_dem, tbl["save_final_dem"] );
+    params.write_thickness_every_n_lobes = tbl["write_thickness_every_n_lobes"].value<int>();
 
     std::optional<std::string> output_folder_string{};
     output_folder_string = tbl["output_folder"].value<std::string>();
@@ -270,6 +271,11 @@ void validate_settings( const InputParams & options )
         "Allowed values of the vent_flag are 0 to 8, inclusive." );
 
     // Output settings validation
+    if( options.write_thickness_every_n_lobes.has_value() )
+    {
+        check( name_and_var( options.write_thickness_every_n_lobes.value() ), g_zero );
+    }
+
     check(
         name_and_var( options.output_settings.compression_level ), []( auto x ) { return x >= 0 && x <= 9; },
         "The compression level can only be between 0 and 9, inclusive." );

src/simulation.cpp

@@ -296,13 +296,13 @@ void Simulation::write_avg_thickness_file()
         xt::filter( topography_masked.hazard, topography_thickness.height_data < threshold_thickness )      = 0.0;
 
         // Write the masked thickness and the masked hazard maps
-        auto file_thick = get_file_handle( topography_masked, OutputQuantitiy::Height );
+        auto file_thick = get_file_handle( topography_masked, OutputQuantity::Height );
         file_thick->save(
             input.output_folder / fmt::format( "{}_thickness_masked_{:.2f}", input.run_name, threshold ) );
 
         if( input.save_hazard_data )
         {
-            auto file_hazard = get_file_handle( topography_masked, OutputQuantitiy::Hazard );
+            auto file_hazard = get_file_handle( topography_masked, OutputQuantity::Hazard );
             file_hazard->save(
                 input.output_folder / fmt::format( "{}_hazard_masked_{:.2f}", input.run_name, threshold ) );
         }
@@ -311,7 +311,7 @@ void Simulation::write_avg_thickness_file()
 }
 
 std::unique_ptr<TopographyFile>
-Simulation::get_file_handle( const Topography & topography, OutputQuantitiy output_quantity ) const
+Simulation::get_file_handle( const Topography & topography, OutputQuantity output_quantity ) const
 {
     std::unique_ptr<TopographyFile> res{};
 
@@ -343,6 +343,33 @@ Simulation::get_file_handle( const Topography & topography, OutputQuantitiy outp
     return res;
 }
 
+void Simulation::compute_topography_thickness()
+{
+    // Compute the thickness by substracting 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;
+    topography_thickness.height_data /= ( 1.0 - input.thickening_parameter );
+}
+
+void Simulation::write_thickness_if_necessary( int n_lobes_processed )
+{
+    // If the optional does not have a value, we simply return without doing anything
+    if( !input.write_thickness_every_n_lobes.has_value() )
+    {
+        return;
+    }
+
+    // Else we check if we have to write a dem file according to our setting
+    if( n_lobes_processed % input.write_thickness_every_n_lobes.value() == 0 )
+    {
+        compute_topography_thickness();
+        auto file_thick = get_file_handle( topography_thickness, OutputQuantity::Height );
+        file_thick->save(
+            input.output_folder / fmt::format( "{}_thickness_after_{}_lobes", input.run_name, n_lobes_processed ) );
+    }
+}
+
 void Simulation::run()
 {
     // Initialize MrLavaLoba method
@@ -394,6 +421,7 @@ void Simulation::run()
             // Add rasterized lobe
             topography.add_lobe( lobe_cur, input.volume_correction, idx_lobe );
             n_lobes_processed++;
+            write_thickness_if_necessary( n_lobes_processed );
         }
 
         // Loop over the rest of the lobes (skipping the initial ones).
@@ -464,6 +492,7 @@ void Simulation::run()
             // 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 )
@@ -501,29 +530,25 @@ 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, OutputQuantitiy::Height );
+    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, OutputQuantitiy::Height );
+        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
-    topography_thickness               = topography;
-    topography_thickness.no_data_value = DEFAULT_NO_DATA_VALUE_THICKNESS;
-    topography_thickness.height_data -= topography_initial.height_data;
-    topography_thickness.height_data /= ( 1.0 - input.thickening_parameter );
-
-    auto file_thick = get_file_handle( topography_thickness, OutputQuantitiy::Height );
+    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, OutputQuantitiy::Hazard );
+        auto file_hazard = get_file_handle( topography, OutputQuantity::Hazard );
         file_hazard->save( input.output_folder / fmt::format( "{}_hazard_full", input.run_name ) );
     }
 

src/topography.cpp

@@ -516,11 +516,8 @@ 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<double>( sin_phi_lobe_perimeter->begin(), sin_phi_lobe_perimeter->end() );
     const auto cos = std::span<double>( cos_phi_lobe_perimeter->begin(), cos_phi_lobe_perimeter->end() );
-
     std::vector<Vector2> perimeter = lobe.rasterize_perimeter( sin, cos );
 
-    // std::vector<Vector2> perimeter = lobe.rasterize_perimeter( npoints );
-
     // Then, we find the point of minimal elevation amongst the rasterized points on the perimeter
     auto min_elevation_point_it = std::min_element(
         perimeter.begin(), perimeter.end(), [&]( const Vector2 & p1, const Vector2 & p2 )