diff --git a/flowy/include/asc_file.hpp b/flowy/include/asc_file.hpp index 79dbd63..50aead7 100644 --- a/flowy/include/asc_file.hpp +++ b/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 ) { } diff --git a/flowy/include/config.hpp b/flowy/include/config.hpp index 3e37c3b..e72f59e 100644 --- a/flowy/include/config.hpp +++ b/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 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 diff --git a/flowy/include/lobe.hpp b/flowy/include/lobe.hpp index 5b17dc6..50f4277 100644 --- a/flowy/include/lobe.hpp +++ b/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 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 rasterize_perimeter( int n_raster_points ) const - { - VectorX phi_list = xt::linspace( 0.0, 2.0 * Math::pi, n_raster_points, false ); - auto res = std::vector( 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 rasterize_perimeter( std::span sin_phi, std::span cos_phi ) const + inline std::vector + rasterize_perimeter( const std::span sin_phi, const std::span cos_phi ) const { const int n_raster_points = sin_phi.size(); auto res = std::vector( n_raster_points ); diff --git a/flowy/include/models/mr_lava_loba.hpp b/flowy/include/models/mr_lava_loba.hpp index cbb00a2..ee01f2e 100644 --- a/flowy/include/models/mr_lava_loba.hpp +++ b/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 f{}; @@ -182,7 +182,7 @@ 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 ); @@ -190,7 +190,7 @@ class MrLavaLoba 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 & lobes ) + int select_parent_lobe( int idx_descendant, std::vector & lobes ) const { Lobe & lobe_descendent = lobes[idx_descendant]; diff --git a/flowy/include/netcdf_file.hpp b/flowy/include/netcdf_file.hpp index 6cd6894..48c0f4f 100644 --- a/flowy/include/netcdf_file.hpp +++ b/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 ) { } diff --git a/flowy/include/simulation.hpp b/flowy/include/simulation.hpp index bc26925..0d72955 100644 --- a/flowy/include/simulation.hpp +++ b/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 - get_file_handle( const Topography & topography, OutputQuantitiy output_quantity ) const; + get_file_handle( const Topography & topography, OutputQuantity output_quantity ) const; void run(); diff --git a/flowy/include/topography_file.hpp b/flowy/include/topography_file.hpp index 182b87d..1a06bb2 100644 --- a/flowy/include/topography_file.hpp +++ b/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; } diff --git a/src/config_parser.cpp b/src/config_parser.cpp index 5c9eb83..d54ed51 100644 --- a/src/config_parser.cpp +++ b/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(); std::optional output_folder_string{}; output_folder_string = tbl["output_folder"].value(); @@ -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." ); diff --git a/src/simulation.cpp b/src/simulation.cpp index e36e36b..7f53e8c 100644 --- a/src/simulation.cpp +++ b/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 -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 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 ) ); } diff --git a/src/topography.cpp b/src/topography.cpp index 11d3652..62df190 100644 --- a/src/topography.cpp +++ b/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( 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 ); - // std::vector 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 ) diff --git a/test/test_file_io.cpp b/test/test_file_io.cpp index 56e74fc..eceed6f 100644 --- a/test/test_file_io.cpp +++ b/test/test_file_io.cpp @@ -58,41 +58,41 @@ TEST_CASE( "file_io", "[netcdf]" ) auto write_asc_height = [&]( const Flowy::Topography & topography, const fs::path & path ) { - auto file = Flowy::AscFile( topography, Flowy::OutputQuantitiy::Height ); + auto file = Flowy::AscFile( topography, Flowy::OutputQuantity::Height ); file.save( path ); }; auto write_asc_hazard = [&]( const Flowy::Topography & topography, const fs::path & path ) { - auto file = Flowy::AscFile( topography, Flowy::OutputQuantitiy::Hazard ); + auto file = Flowy::AscFile( topography, Flowy::OutputQuantity::Hazard ); file.save( path ); }; auto read_asc = [&]( const fs::path & path ) { return Flowy::AscFile( path ); }; auto write_netcdf_height = [&]( const Flowy::Topography & topography, const fs::path & path ) { - auto file = Flowy::NetCDFFile( topography, Flowy::OutputQuantitiy::Height ); + auto file = Flowy::NetCDFFile( topography, Flowy::OutputQuantity::Height ); file.data_type = Flowy::StorageDataType::Double; file.save( path ); }; auto write_netcdf_height_float = [&]( const Flowy::Topography & topography, const fs::path & path ) { - auto file = Flowy::NetCDFFile( topography, Flowy::OutputQuantitiy::Height ); + auto file = Flowy::NetCDFFile( topography, Flowy::OutputQuantity::Height ); file.data_type = Flowy::StorageDataType::Float; file.save( path ); }; auto write_netcdf_height_short = [&]( const Flowy::Topography & topography, const fs::path & path ) { - auto file = Flowy::NetCDFFile( topography, Flowy::OutputQuantitiy::Height ); + auto file = Flowy::NetCDFFile( topography, Flowy::OutputQuantity::Height ); file.data_type = Flowy::StorageDataType::Short; file.save( path ); }; auto write_netcdf_hazard = [&]( const Flowy::Topography & topography, const fs::path & path ) { - auto file = Flowy::NetCDFFile( topography, Flowy::OutputQuantitiy::Hazard ); + auto file = Flowy::NetCDFFile( topography, Flowy::OutputQuantity::Hazard ); file.save( path ); }; diff --git a/test/test_topography.cpp b/test/test_topography.cpp index aa66de1..4a1d36d 100644 --- a/test/test_topography.cpp +++ b/test/test_topography.cpp @@ -140,13 +140,6 @@ TEST_CASE( "find_preliminary_budding_point", "[budding_point]" ) my_lobe.semi_axes = { 0.8, 0.8 }; my_lobe.set_azimuthal_angle( 0.0 ); - auto perimeter = my_lobe.rasterize_perimeter( 32 ); - for( auto & p : perimeter ) - { - INFO( fmt::format( "p = {}\n", fmt::streamed( p ) ) ); - INFO( fmt::format( "height = {}\n\n", topography.height_and_slope( p ).first ) ); - } - Flowy::Vector2 budding_point = topography.find_preliminary_budding_point( my_lobe, 32 ); INFO( fmt::format( "budding_point = {}", budding_point ) );