From 15379cd6c35c63bd3ab97573229c780cc319c2b6 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Wed, 26 Mar 2025 20:37:53 +0000 Subject: [PATCH 1/9] bindings: added Simulation::steps bindings: mr_lava_loba in new submodule --- python_bindings/bindings.cpp | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/python_bindings/bindings.cpp b/python_bindings/bindings.cpp index 52beb6d..7e6604a 100644 --- a/python_bindings/bindings.cpp +++ b/python_bindings/bindings.cpp @@ -219,16 +219,8 @@ PYBIND11_MODULE( flowycpp, m ) .def_readwrite( "topography", &Flowy::Simulation::topography ) .def_readwrite( "lobes", &Flowy::Simulation::lobes ) .def( "stop_condition", &Flowy::Simulation::stop_condition ) - .def( "run", &Flowy::Simulation::run ); - - py::class_( m, "MrLavaLoba" ) - .def_readwrite( "lobe_dimensions", &Flowy::MrLavaLoba::lobe_dimensions ) - .def( "compute_initial_lobe_position", &Flowy::MrLavaLoba::compute_initial_lobe_position ) - .def( "compute_lobe_axes", &Flowy::MrLavaLoba::compute_lobe_axes ) - .def( "compute_descendent_lobe_position", &Flowy::MrLavaLoba::compute_descendent_lobe_position ) - .def( "perturb_lobe_angle", &Flowy::MrLavaLoba::perturb_lobe_angle ) - .def( "select_parent_lobe", &Flowy::MrLavaLoba::select_parent_lobe ) - .def( "add_inertial_contribution", &Flowy::MrLavaLoba::add_inertial_contribution ); + .def( "run", &Flowy::Simulation::run ) + .def( "steps", &Flowy::Simulation::steps ); m.def( "parse_config", &Flowy::Config::parse_config, "A function to parse input settings from a TOML file.", @@ -236,4 +228,13 @@ PYBIND11_MODULE( flowycpp, m ) m.def( "validate_settings", &Flowy::Config::validate_settings, "A function to validate the Flowy config settings"_a ); + + auto mr_lava_loba_m = m.def_submodule( "mr_lava_loba" ); + mr_lava_loba_m.def( "compute_initial_lobe_position", &Flowy::MrLavaLoba::compute_initial_lobe_position ); + mr_lava_loba_m.def( "compute_lobe_axes", &Flowy::MrLavaLoba::compute_lobe_axes ); + mr_lava_loba_m.def( "compute_descendent_lobe_position", &Flowy::MrLavaLoba::compute_descendent_lobe_position ); + mr_lava_loba_m.def( "perturb_lobe_angle", &Flowy::MrLavaLoba::perturb_lobe_angle ); + mr_lava_loba_m.def( "select_parent_lobe", &Flowy::MrLavaLoba::select_parent_lobe ); + mr_lava_loba_m.def( "add_inertial_contribution", &Flowy::MrLavaLoba::add_inertial_contribution ); + } From b7956f5bf53f64fb0d2760ab5718681ade15f392 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 27 Mar 2025 01:13:48 +0000 Subject: [PATCH 2/9] ENH: added bindings for - SimulationState - get_simulation_state - reset_simulation_state ENH: unit test for simualiton with steps ENH: use pytest fixtures Co-authored-by: Amrita Goswami --- python_bindings/bindings.cpp | 13 +++++++++++- subprojects/flowy.wrap | 2 +- tests/conftest.py | 12 ++++++++++- tests/test_simulation.py | 39 +++++++++++++++++++++++++++++------- 4 files changed, 56 insertions(+), 10 deletions(-) diff --git a/python_bindings/bindings.cpp b/python_bindings/bindings.cpp index 7e6604a..e09a79e 100644 --- a/python_bindings/bindings.cpp +++ b/python_bindings/bindings.cpp @@ -213,11 +213,23 @@ PYBIND11_MODULE( flowycpp, m ) .def_readwrite( "max_semiaxis", &Flowy::CommonLobeDimensions::max_semiaxis ) .def_readwrite( "thickness_min", &Flowy::CommonLobeDimensions::thickness_min ); + py::enum_( m, "RunStatus" ) + .value( "Finished", Flowy::RunStatus::Finished ) + .value( "Ongoing", Flowy::RunStatus::Ongoing ); + + py::class_( m, "SimulationState" ) + .def_readonly( "n_lobes_processed", &Flowy::SimulationState::n_lobes_processed ) + .def_readonly( "n_lobes", &Flowy::SimulationState::n_lobes ) + .def_readonly( "idx_flow", &Flowy::SimulationState::idx_flow ) + .def_readonly( "idx_lobe", &Flowy::SimulationState::idx_lobe ); + py::class_( m, "Simulation" ) .def( py::init>() ) .def_readwrite( "input", &Flowy::Simulation::input ) .def_readwrite( "topography", &Flowy::Simulation::topography ) .def_readwrite( "lobes", &Flowy::Simulation::lobes ) + .def( "reset_simulation_state", &Flowy::Simulation::reset_simulation_state ) + .def( "get_simulation_state", &Flowy::Simulation::get_simulation_state ) .def( "stop_condition", &Flowy::Simulation::stop_condition ) .def( "run", &Flowy::Simulation::run ) .def( "steps", &Flowy::Simulation::steps ); @@ -236,5 +248,4 @@ PYBIND11_MODULE( flowycpp, m ) mr_lava_loba_m.def( "perturb_lobe_angle", &Flowy::MrLavaLoba::perturb_lobe_angle ); mr_lava_loba_m.def( "select_parent_lobe", &Flowy::MrLavaLoba::select_parent_lobe ); mr_lava_loba_m.def( "add_inertial_contribution", &Flowy::MrLavaLoba::add_inertial_contribution ); - } diff --git a/subprojects/flowy.wrap b/subprojects/flowy.wrap index 9e1d7f0..c3e5a0a 100644 --- a/subprojects/flowy.wrap +++ b/subprojects/flowy.wrap @@ -1,5 +1,5 @@ [wrap-git] directory=flowy url=https://github.com/flowy-code/flowy.git -revision=8520829e141edda73547eeaf063003e269c7c44e +revision=8407b8a5c44bd286836d0f7a234f240cede6aa80 depth=1 \ No newline at end of file diff --git a/tests/conftest.py b/tests/conftest.py index 3ae039c..5a2a2ff 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -10,4 +10,14 @@ def hawaii_example() -> tuple[Path, Path]: asc_filepath = Path(f"{FILEPATH}/hawaii/test20m.asc") util.download_file_from_github( toml_filepath, relative_file_path="examples/KILAUEA2014-2015/input.toml" ) util.download_file_from_github( asc_filepath, relative_file_path="examples/KILAUEA2014-2015/test20m.asc" ) - return toml_filepath, asc_filepath \ No newline at end of file + return toml_filepath, asc_filepath + +@pytest.fixture +def output_folder() -> Path: + output_folder = Path(f"{FILEPATH}/output_test") + # Remove any previous files because otherwise you just get more extra files and the old ones persist + if output_folder.exists(): + [f.unlink() for f in output_folder.glob("*")] + output_folder.rmdir() + + return output_folder \ No newline at end of file diff --git a/tests/test_simulation.py b/tests/test_simulation.py index db8bbaf..2c93be3 100644 --- a/tests/test_simulation.py +++ b/tests/test_simulation.py @@ -6,7 +6,7 @@ FILEPATH = Path(__file__).parent -def test_simulation(hawaii_example): +def test_simulation(hawaii_example, output_folder): toml_filepath, asc_filepath = hawaii_example @@ -14,12 +14,6 @@ def test_simulation(hawaii_example): input_params.write_thickness_every_n_lobes = 5000 pfy.flowycpp.validate_settings(input_params) - output_folder = Path(f"{FILEPATH}/output_test") - # Remove any previous files because otherwise you just get more extra files and the old ones persist - if output_folder.exists(): - [f.unlink() for f in output_folder.glob("*")] - output_folder.rmdir() - input_params.output_folder = output_folder input_params.source = asc_filepath @@ -39,3 +33,34 @@ def test_simulation(hawaii_example): assert Path(output_folder/f"{input_params.run_name}_hazard_masked_{input_params.masking_threshold[-1]}.nc").exists() # Make sure that this summary file was written assert Path(output_folder / f"{input_params.run_name}_avg_thick.txt").exists() + + +def test_simulation_with_steps(hawaii_example, output_folder): + toml_filepath, asc_filepath = hawaii_example + + input_params = pfy.flowycpp.parse_config( toml_filepath ) + input_params.write_thickness_every_n_lobes = 100 + pfy.flowycpp.validate_settings(input_params) + + input_params.output_folder = output_folder + input_params.source = asc_filepath + + # Change some of the output settings + input_params.output_settings.use_netcdf = True + input_params.output_settings.crop_to_content = True + input_params.output_settings.data_type = pfy.flowycpp.StorageDataType.Short + + simulation = pfy.flowycpp.Simulation(input_params, 0) + + simulation.input.n_flows = 10 + simulation.input.min_n_lobes = 50 + simulation.input.max_n_lobes = 50 + simulation.input.max_slope_prob = 0.0 + simulation.input.lobe_exponent = 0.0 + simulation.input.inertial_exponent = 0.0 + simulation.input.thickening_parameter = 0.0 + + while simulation.steps(100) == pfy.flowycpp.RunStatus.Ongoing: + simulation.input.max_slope_prob += 0.1 + + assert Path(output_folder/f"{input_params.run_name}_thickness_after_{input_params.write_thickness_every_n_lobes}_lobes.nc").exists() \ No newline at end of file From 0c74cead55fe316859c181d12310613e90076165 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 27 Mar 2025 21:55:47 +0000 Subject: [PATCH 3/9] ENH: bind mt19937 and DEFAULT_NO_DATA_VALUE_* the mt19937 from the cpp std lib is needed in the restructured mr_lava_loba functions and the DEFAULT_NO_DATA_VALUE_ are used in the constructor of the Topography class. Co-authored-by: Amrita Goswami --- python_bindings/bindings.cpp | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/python_bindings/bindings.cpp b/python_bindings/bindings.cpp index e09a79e..6f5c1ca 100644 --- a/python_bindings/bindings.cpp +++ b/python_bindings/bindings.cpp @@ -9,6 +9,7 @@ #include "flowy/include/topography_file.hpp" #include "pybind11/detail/common.h" #include "pybind11/pytypes.h" +#include #ifdef WITH_NETCDF #include "flowy/include/netcdf_file.hpp" @@ -70,7 +71,7 @@ PYBIND11_MODULE( flowycpp, m ) #ifdef WITH_NETCDF py::class_( m, "NetCDFFile" ) .def( py::init<>() ) - .def( py::init() ) + .def( py::init() ) .def( py::init(), "file_path"_a ) .def( py::init() ) .def( "save", &Flowy::NetCDFFile::save ) @@ -112,6 +113,9 @@ PYBIND11_MODULE( flowycpp, m ) .def_readwrite( "cells_intersecting", &Flowy::LobeCells::cells_intersecting ) .def_readwrite( "cells_enclosed", &Flowy::LobeCells::cells_enclosed ); + m.attr( "DEFAULT_NO_DATA_VALUE_HEIGHT" ) = py::float_( Flowy::DEFAULT_NO_DATA_VALUE_HEIGHT ); + m.attr( "DEFAULT_NO_DATA_VALUE_THICKNESS" ) = py::float_( Flowy::DEFAULT_NO_DATA_VALUE_THICKNESS ); + py::class_( m, "Topography" ) .def( py::init<>() ) .def( py::init() ) @@ -241,6 +245,11 @@ PYBIND11_MODULE( flowycpp, m ) m.def( "validate_settings", &Flowy::Config::validate_settings, "A function to validate the Flowy config settings"_a ); + py::class_( m, "mt19937" ) + .def( py::init<>() ) + .def( py::init() ) + .def( "__call__", []( std::mt19937 & p ) { return p(); } ); + auto mr_lava_loba_m = m.def_submodule( "mr_lava_loba" ); mr_lava_loba_m.def( "compute_initial_lobe_position", &Flowy::MrLavaLoba::compute_initial_lobe_position ); mr_lava_loba_m.def( "compute_lobe_axes", &Flowy::MrLavaLoba::compute_lobe_axes ); From bad213a563760e44d52513db9ae3ed241f5968e6 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 27 Mar 2025 21:57:27 +0000 Subject: [PATCH 4/9] BLD: read the `with_netcdf` value from flowy Since the compilation of the netcdf file classes is controlled with via the preprocessor (`#ifdef WITH_NETCDF`), we also have to define this preprocessor variable in the bindings code. This is needed to decide if we have to generate the bindings for the netcdf files. Co-authored-by: Amrita Goswami --- meson.build | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/meson.build b/meson.build index 9a7e456..f683b5e 100644 --- a/meson.build +++ b/meson.build @@ -13,6 +13,10 @@ flowy_proj = subproject('flowy', default_options: ['default_library=static', flowylib = flowy_proj.get_variable('flowylib_static_dep') deps += [flowylib] +if flowy_proj.get_variable('with_netcdf') + args += ['-DWITH_NETCDF'] +endif + py = import('python').find_installation('python', modules: ['numpy']) deps += py.dependency() From 565a7607f219deb77ff46d5d5f2e960f0f4b3f6d Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 27 Mar 2025 21:59:59 +0000 Subject: [PATCH 5/9] ENH: update examples So they work with the new bindings. Co-authored-by: Amrita Goswami --- examples/add_lobe_example.py | 25 ++++-- examples/compute_cell_intersection.py | 16 ++-- examples/create_toy_topographies.py | 60 +++++-------- examples/descendent_lobe.py | 104 +++++++++++++++------- examples/line_segment_intersects.py | 4 +- examples/raster.py | 12 ++- examples/simulation.py | 122 ++++++++++++++++++-------- examples/topography_interpolation.py | 7 +- 8 files changed, 222 insertions(+), 128 deletions(-) diff --git a/examples/add_lobe_example.py b/examples/add_lobe_example.py index 938dfbc..1af93fe 100644 --- a/examples/add_lobe_example.py +++ b/examples/add_lobe_example.py @@ -10,17 +10,21 @@ extent = lobe.extent_xy() -perimeter = np.array(lobe.rasterize_perimeter(30)) +perimeter = np.array( + [lobe.point_at_angle(phi) for phi in np.linspace(0, 2 * np.pi, 30, endpoint=True)] +) -x_data = np.linspace(0, 40, 40) -y_data = np.linspace(0, 20, 20) +x_data = np.linspace(0, 40, 40, endpoint=False) +y_data = np.linspace(0, 20, 20, endpoint=False) height_data = np.zeros(shape=(len(x_data), len(y_data))) height_data = np.array( [[i + j for j in range(len(y_data))] for i in range(len(x_data))] ) -topography = pfy.flowycpp.Topography(height_data, x_data, y_data) +topography = pfy.flowycpp.Topography( + height_data, x_data, y_data, pfy.flowycpp.DEFAULT_NO_DATA_VALUE_HEIGHT +) for p in perimeter: print(topography.height_and_slope(p)) @@ -33,15 +37,20 @@ new_lobe.thickness = 20 new_lobe.center = [20, 10] new_lobe.semi_axes = [8, 2] -new_lobe_perimeter = np.array(new_lobe.rasterize_perimeter(30)) +new_lobe_perimeter = np.array( + [ + new_lobe.point_at_angle(phi) + for phi in np.linspace(0, 2 * np.pi, 30, endpoint=True) + ] +) new_lobe_extent = new_lobe.extent_xy() -topography.add_lobe(lobe, None) -topography.add_lobe(new_lobe, None) +topography.add_lobe(lobe, False, None) +topography.add_lobe(new_lobe, False, None) cell = topography.cell_size() -plt.pcolormesh(x_data+0.5*cell, y_data+0.5*cell, height_data.T) +plt.pcolormesh(x_data + 0.5 * cell, y_data + 0.5 * cell, height_data.T) plt.axvline(lobe.center[0] + extent[0], color="black") plt.axvline(lobe.center[0] - extent[0], color="black") diff --git a/examples/compute_cell_intersection.py b/examples/compute_cell_intersection.py index 6dc7377..266c6bf 100644 --- a/examples/compute_cell_intersection.py +++ b/examples/compute_cell_intersection.py @@ -8,6 +8,7 @@ lobe.set_azimuthal_angle(np.pi / 4) lobe.center = [20, 10] + def add_lobe(lobe, topography, height_data): # height_data = np.zeros_like(topography.height_data) lobecells = topography.get_cells_intersecting_lobe(lobe, None) @@ -19,19 +20,24 @@ def add_lobe(lobe, topography, height_data): height_data[c[0], c[1]] += 30 return height_data + extent = lobe.extent_xy() -perimeter = np.array(lobe.rasterize_perimeter(30)) +perimeter = np.array( + [lobe.point_at_angle(phi) for phi in np.linspace(0, 2 * np.pi, 30, endpoint=True)] +) -x_data = np.linspace(0, 40, 40) -y_data = np.linspace(0, 20, 20) +x_data = np.linspace(0, 40, 40, endpoint=False) +y_data = np.linspace(0, 20, 20, endpoint=False) height_data = np.zeros(shape=(len(x_data), len(y_data))) height_data = np.array( [[i + j for j in range(len(y_data))] for i in range(len(x_data))] ) -topography = pfy.flowycpp.Topography(height_data, x_data, y_data) +topography = pfy.flowycpp.Topography( + height_data, x_data, y_data, pfy.flowycpp.DEFAULT_NO_DATA_VALUE_HEIGHT +) bbox = topography.bounding_box(lobe.center, extent[0], extent[1]) @@ -41,7 +47,7 @@ def add_lobe(lobe, topography, height_data): # Plot cell = topography.cell_size() # plt.pcolormesh(x_data+0.5*cell, y_data+0.5*cell, height_data.T) -plt.pcolormesh(x_data+0.5*cell, y_data+0.5*cell, new_heights.T) +plt.pcolormesh(x_data + 0.5 * cell, y_data + 0.5 * cell, new_heights.T) plt.axvline(x_data[bbox.idx_x_lower], color="black") plt.axvline(x_data[bbox.idx_x_higher], color="black") plt.axhline(y_data[bbox.idx_y_lower], color="black") diff --git a/examples/create_toy_topographies.py b/examples/create_toy_topographies.py index 7dc9518..e81d2c9 100644 --- a/examples/create_toy_topographies.py +++ b/examples/create_toy_topographies.py @@ -1,60 +1,40 @@ import pyflowy as pfy import numpy as np -import matplotlib.pyplot as plt -# First we create a topograhy with constant sloe -x_data = np.linspace(0, 10000, 1000) -y_data = np.linspace(0, 10000, 1000) +def create_topography( + height_cb, filename, xdim=10000, ydim=10000, res=10, lower_left_corner=[0, 0] +): -def height(x, y): - return 5e-2 * x - + # First we create a topograhy with constant slope + x_data = np.arange(lower_left_corner[0], xdim + lower_left_corner[0], res) + y_data = np.arange(lower_left_corner[1], ydim + lower_left_corner[1], res) -height_data = np.array([[height(x, y) for y in y_data] for x in x_data]) + height_data = np.array([[height_cb(x, y) for y in y_data] for x in x_data]) -asc_file = pfy.flowycpp.AscFile() -asc_file.x_data = x_data -asc_file.y_data = y_data -asc_file.height_data = height_data -asc_file.cell_size = x_data[1] - x_data[0] -asc_file.save("constant_slope.asc") + file = pfy.flowycpp.NetCDFFile() + file.x_data = x_data + file.y_data = y_data + file.data = height_data + file.save(filename) -# Then we create a saddle topography -x_data = np.linspace(0, 10000, 1000) -y_data = np.linspace(0, 10000, 1000) +def height_const_slope(x, y): + return 5e-2 * x -def height(x, y, center): - return -1e-4 * (x - center[0]) * (y - center[1]) +create_topography(height_const_slope, filename="constant_slope") -center = np.array([np.mean(x_data), np.mean(y_data)]) -height_data = np.array([[height(x, y, center) for y in y_data] for x in x_data]) +def height_parabola(x, y): + return -1e-4 * (x - 5000) * (y - 5000) -asc_file = pfy.flowycpp.AscFile() -asc_file.x_data = x_data -asc_file.y_data = y_data -asc_file.height_data = height_data -asc_file.cell_size = x_data[1] - x_data[0] -asc_file.save("saddle.asc") -# Then we create a flat topography -x_data = np.linspace(0, 50000, 5000) -y_data = np.linspace(0, 50000, 5000) +create_topography(height_parabola, filename="parabola") -def height(x, y, center): +def height_flat(x, y): return 0.0 -center = np.array([np.mean(x_data), np.mean(y_data)]) -height_data = np.array([[height(x, y, center) for y in y_data] for x in x_data]) - -asc_file = pfy.flowycpp.AscFile() -asc_file.x_data = x_data -asc_file.y_data = y_data -asc_file.height_data = height_data -asc_file.cell_size = x_data[1] - x_data[0] -asc_file.save("flat.asc") +create_topography(height_flat, filename="flat") diff --git a/examples/descendent_lobe.py b/examples/descendent_lobe.py index 573251d..fe74b3d 100644 --- a/examples/descendent_lobe.py +++ b/examples/descendent_lobe.py @@ -3,24 +3,31 @@ import matplotlib.pyplot as plt import math -def compute_descendent_lobe_position( descendent_lobe, parent_lobe, final_budding_point ): - direction_to_new_lobe = ( final_budding_point - parent_lobe.center ) / np.linalg.norm( final_budding_point - parent_lobe.center ) - new_lobe_center = final_budding_point + direction_to_new_lobe * descendent_lobe.semi_axes[0] + +def compute_descendent_lobe_position(descendent_lobe, parent_lobe, final_budding_point): + direction_to_new_lobe = (final_budding_point - parent_lobe.center) / np.linalg.norm( + final_budding_point - parent_lobe.center + ) + new_lobe_center = ( + final_budding_point + direction_to_new_lobe * descendent_lobe.semi_axes[0] + ) return new_lobe_center -def compute_lobe_axes( lobe, slope, input, lobe_dimensions ): - slope_norm = np.linalg.norm( slope, 2 ) + +def compute_lobe_axes(lobe, slope, input, lobe_dimensions): + slope_norm = np.linalg.norm(slope, 2) # Factor for the lobe eccentricity - aspect_ratio = min( input.max_aspect_ratio, 1.0 + input.aspect_ratio_coeff * slope_norm ) - # breakpoint() - # aspect_ratio = 2.0 + aspect_ratio = min( + input.max_aspect_ratio, 1.0 + input.aspect_ratio_coeff * slope_norm + ) # Compute the semi-axes of the lobe - semi_major_axis = np.sqrt( lobe_dimensions.lobe_area / np.pi ) * np.sqrt( aspect_ratio ) - semi_minor_axis = np.sqrt( lobe_dimensions.lobe_area / np.pi ) / np.sqrt( aspect_ratio ) + semi_major_axis = np.sqrt(lobe_dimensions.lobe_area / np.pi) * np.sqrt(aspect_ratio) + semi_minor_axis = np.sqrt(lobe_dimensions.lobe_area / np.pi) / np.sqrt(aspect_ratio) return semi_major_axis, semi_minor_axis + # Input parameters (we just construct this here) input = pfy.flowycpp.InputParams() input.source = "file.asc" @@ -40,13 +47,15 @@ def compute_lobe_axes( lobe, slope, input, lobe_dimensions ): # Create a Simulation object simulation = pfy.flowycpp.Simulation(input, None) # Topography data -x_data = np.linspace(0, 40, 40) -y_data = np.linspace(0, 20, 20) +x_data = np.linspace(0, 40, 40, endpoint=False) +y_data = np.linspace(0, 20, 20, endpoint=False) height_data = np.zeros(shape=(len(x_data), len(y_data))) height_data = np.array( [[i + j for j in range(len(y_data))] for i in range(len(x_data))] ) -topography = pfy.flowycpp.Topography(height_data, x_data, y_data) +topography = pfy.flowycpp.Topography( + height_data, x_data, y_data, pfy.flowycpp.DEFAULT_NO_DATA_VALUE_HEIGHT +) simulation.topography = topography # Parent lobe @@ -62,42 +71,75 @@ def compute_lobe_axes( lobe, slope, input, lobe_dimensions ): height_parent_center, slope_parent_center = topography.height_and_slope( parent_lobe.center ) -simulation.perturb_lobe_angle(parent_lobe, slope_parent_center) + +gen = pfy.flowycpp.mt19937() + +pfy.flowycpp.mr_lava_loba.perturb_lobe_angle( + parent_lobe, np.linalg.norm(slope_parent_center), simulation.input, gen +) + # Perimeter of the parent lobe, rasterized so you can see it -perimeter = np.array(parent_lobe.rasterize_perimeter(30)) -print("Azimuthal angle of the parent in pi", parent_lobe.get_azimuthal_angle()/np.pi) +perimeter = np.array( + [ + parent_lobe.point_at_angle(phi) + for phi in np.linspace(0, 2 * np.pi, 30, endpoint=True) + ] +) +print("Azimuthal angle of the parent in pi", parent_lobe.get_azimuthal_angle() / np.pi) #### Now we are done changing the parent lobe budding_point = topography.find_preliminary_budding_point(parent_lobe, 30) -simulation.perturb_lobe_angle(descendent_lobe, slope_parent_center) # This should have set the azimuthal angle +pfy.flowycpp.mr_lava_loba.perturb_lobe_angle( + descendent_lobe, np.linalg.norm(slope_parent_center), simulation.input, gen +) # This should have set the azimuthal angle -simulation.add_inertial_contribution( descendent_lobe, parent_lobe, slope_parent_center ) +pfy.flowycpp.mr_lava_loba.add_inertial_contribution( + descendent_lobe, parent_lobe, np.linalg.norm(slope_parent_center), simulation.input +) -angle_diff = parent_lobe.get_azimuthal_angle() - descendent_lobe.get_azimuthal_angle() -final_budding_point = parent_lobe.point_at_angle( np.pi - angle_diff ) -height_budding_point, slope_budding_point = topography.height_and_slope( final_budding_point ) -print("Azimuthal angle of the descendent in pi ", descendent_lobe.get_azimuthal_angle()/np.pi) +angle_diff = parent_lobe.get_azimuthal_angle() - descendent_lobe.get_azimuthal_angle() +final_budding_point = parent_lobe.point_at_angle(np.pi - angle_diff) +height_budding_point, slope_budding_point = topography.height_and_slope( + final_budding_point +) +print( + "Azimuthal angle of the descendent in pi ", + descendent_lobe.get_azimuthal_angle() / np.pi, +) + +common_lobe_dimensions = pfy.flowycpp.CommonLobeDimensions(simulation.input) # Compute the lobe axes and center of the new descendent lobe -simulation.compute_lobe_axes( descendent_lobe, slope_budding_point ) -# major_axis, minor_axis = compute_lobe_axes( descendent_lobe, slope_budding_point, input, lobe_dimensions ) +pfy.flowycpp.mr_lava_loba.compute_lobe_axes( + descendent_lobe, + np.linalg.norm(slope_budding_point), + simulation.input, + common_lobe_dimensions, +) + # descendent_lobe.semi_axes = [major_axis, minor_axis] print("Axes = ", descendent_lobe.semi_axes) -# Compute the final descendent center -new_lobe_center = compute_descendent_lobe_position( descendent_lobe, parent_lobe, final_budding_point ) -descendent_lobe.center = new_lobe_center +# Compute the final descendent center +new_lobe_center = compute_descendent_lobe_position( + descendent_lobe, parent_lobe, final_budding_point +) +descendent_lobe.center = new_lobe_center print("new lobe center ", new_lobe_center) -perimeter_descendent = np.array(descendent_lobe.rasterize_perimeter(30)) +perimeter_descendent = np.array( + [descendent_lobe.point_at_angle(phi) for phi in np.linspace(0, 2 * np.pi, 30)] +) # Representation of the topography plt.pcolormesh(x_data, y_data, topography.height_data.T) -plt.plot(perimeter[:, 0], perimeter[:, 1], color="black") # parent lobe -plt.plot(perimeter_descendent[:, 0], perimeter_descendent[:, 1], color="red") # descendent lobe +plt.plot(perimeter[:, 0], perimeter[:, 1], color="black") # parent lobe +plt.plot( + perimeter_descendent[:, 0], perimeter_descendent[:, 1], color="red" +) # descendent lobe # Budding points plt.plot(budding_point[0], budding_point[1], marker="o", color="black", ms=12) plt.plot(final_budding_point[0], final_budding_point[1], marker=".", color="red", ms=12) plt.plot(new_lobe_center[0], new_lobe_center[1], marker=".", color="green", ms=12) plt.gca().set_box_aspect(1) -plt.show() plt.savefig("descendent.png", dpi=300) +plt.show() diff --git a/examples/line_segment_intersects.py b/examples/line_segment_intersects.py index e42e678..c75ffe4 100644 --- a/examples/line_segment_intersects.py +++ b/examples/line_segment_intersects.py @@ -30,8 +30,10 @@ def plot_lines_intersect(x1, x2): # extent = lobe.extent_xy() -perimeter = np.array(lobe.rasterize_perimeter(30)) +perimeter = np.array( + [lobe.point_at_angle(phi) for phi in np.linspace(0, 2 * np.pi, 30, endpoint=True)] +) # Plot plt.plot(perimeter[:, 0], perimeter[:, 1], color="black") diff --git a/examples/raster.py b/examples/raster.py index 9b1a81c..0f519fb 100644 --- a/examples/raster.py +++ b/examples/raster.py @@ -19,15 +19,19 @@ def add_lobe(lobe, topography, new_height_data): extent = lobe.extent_xy() -perimeter = np.array(lobe.rasterize_perimeter(30)) +perimeter = np.array( + [lobe.point_at_angle(phi) for phi in np.linspace(0, 2 * np.pi, 30, endpoint=True)] +) -x_data = np.linspace(0, 40, 40) -y_data = np.linspace(0, 20, 20) +x_data = np.linspace(0, 40, 40, endpoint=False) +y_data = np.linspace(0, 20, 20, endpoint=False) height_data = np.zeros(shape=(len(x_data), len(y_data))) height_data = np.array([[0 for j in range(len(y_data))] for i in range(len(x_data))]) -topography = pfy.flowycpp.Topography(height_data, x_data, y_data) +topography = pfy.flowycpp.Topography( + height_data, x_data, y_data, pfy.flowycpp.DEFAULT_NO_DATA_VALUE_HEIGHT +) bbox = topography.bounding_box(lobe.center, extent[0], extent[1]) diff --git a/examples/simulation.py b/examples/simulation.py index d35d177..4c5a1b3 100644 --- a/examples/simulation.py +++ b/examples/simulation.py @@ -4,24 +4,31 @@ from matplotlib.colors import Normalize from matplotlib.cm import ScalarMappable + def select_parent_lobe(idx_descendent): - return idx_descendent-1 + return idx_descendent - 1 + -def compute_descendent_lobe_position( descendent_lobe, parent_lobe, final_budding_point ): - direction_to_new_lobe = ( final_budding_point - parent_lobe.center ) / np.linalg.norm( final_budding_point - parent_lobe.center ) - new_lobe_center = final_budding_point + direction_to_new_lobe * descendent_lobe.semi_axes[0] +def compute_descendent_lobe_position(descendent_lobe, parent_lobe, final_budding_point): + direction_to_new_lobe = (final_budding_point - parent_lobe.center) / np.linalg.norm( + final_budding_point - parent_lobe.center + ) + new_lobe_center = ( + final_budding_point + direction_to_new_lobe * descendent_lobe.semi_axes[0] + ) return new_lobe_center + # Input parameters (we just construct this here) input = pfy.flowycpp.InputParams() input.source = "file.asc" input.total_volume = 20 input.prescribed_lobe_area = 20 -input.vent_coordinates = np.array([[20,10]]) -input.npoints = 30 # ellipse rasterization +input.vent_coordinates = np.array([[20, 10]]) +input.npoints = 30 # ellipse rasterization input.prescribed_avg_lobe_thickness = 1 -input.n_init = 1 # One initial lobe -# Required for calculating the thickness +input.n_init = 1 # One initial lobe +# Required for calculating the thickness input.n_flows = 1 input.min_n_lobes = 1 input.max_n_lobes = 1 @@ -37,7 +44,7 @@ def compute_descendent_lobe_position( descendent_lobe, parent_lobe, final_buddin # max slope direction; # inertial_exponent > 0 => the max probability direction for the new lobe takes # into account also the direction of the parent lobe and -# the inertia increaes with increasing exponent +# the inertia increase with increasing exponent input.inertial_exponent = 0.0 # Aspect ratio info required to create the descendent lobe axes input.max_aspect_ratio = 2.5 @@ -45,78 +52,119 @@ def compute_descendent_lobe_position( descendent_lobe, parent_lobe, final_buddin # Create a Simulation object simulation = pfy.flowycpp.Simulation(input, None) +gen = pfy.flowycpp.mt19937() + # Topography data -x_data = np.linspace(0, 40, 40) -y_data = np.linspace(0, 20, 20) +x_data = np.linspace(0, 40, 40, endpoint=False) +y_data = np.linspace(0, 20, 20, endpoint=False) height_data = np.zeros(shape=(len(x_data), len(y_data))) height_data = np.array( [[i + j for j in range(len(y_data))] for i in range(len(x_data))] ) -topography = pfy.flowycpp.Topography(height_data, x_data, y_data) +topography = pfy.flowycpp.Topography( + height_data, x_data, y_data, pfy.flowycpp.DEFAULT_NO_DATA_VALUE_HEIGHT +) simulation.topography = topography -# Reproduce the simulation.run loop +# Reproduce the simulation.run loop n_lobes = 3 -lobes = [] # List, equivalent to simulation.lobes +lobes = [] # List, equivalent to simulation.lobes budding_point_list = [] -thickness = 1 # All lobes have the same thickness +thickness = 1 # All lobes have the same thickness + +common_lobe_dimensions = pfy.flowycpp.CommonLobeDimensions(simulation.input) for idx_lobe in range(input.n_init): lobe_cur = pfy.flowycpp.Lobe() - - simulation.compute_initial_lobe_position( 0, lobe_cur ) + + pfy.flowycpp.mr_lava_loba.compute_initial_lobe_position( + 0, lobe_cur, simulation.input, gen + ) lobe_cur_thickness = thickness - height_lobe_center, slope_lobe_center = topography.height_and_slope( lobe_cur.center ) - simulation.perturb_lobe_angle( lobe_cur, slope_lobe_center ) - simulation.compute_lobe_axes( lobe_cur, slope_lobe_center ) + height_lobe_center, slope_lobe_center = topography.height_and_slope(lobe_cur.center) + + pfy.flowycpp.mr_lava_loba.perturb_lobe_angle( + lobe_cur, np.linalg.norm(slope_lobe_center), simulation.input, gen + ) + pfy.flowycpp.mr_lava_loba.compute_lobe_axes( + lobe_cur, + np.linalg.norm(slope_lobe_center), + simulation.input, + common_lobe_dimensions, + ) + lobes.append(lobe_cur) # Skip initial lobes and go over the rest for idx_lobe in range(input.n_init, n_lobes): lobe_cur = pfy.flowycpp.Lobe() - idx_parent = select_parent_lobe( idx_lobe ) + idx_parent = select_parent_lobe(idx_lobe) lobe_parent = lobes[idx_parent] - preliminary_budding_point = topography.find_preliminary_budding_point( lobe_parent, input.npoints ) - height_lobe_center, slope_parent = topography.height_and_slope( lobe_parent.center ) + preliminary_budding_point = topography.find_preliminary_budding_point( + lobe_parent, input.npoints + ) + height_lobe_center, slope_parent = topography.height_and_slope(lobe_parent.center) # Perturb and add inertial contribution (also set the azimuthal angle of the descendent lobe in these functions) - simulation.perturb_lobe_angle( lobe_cur, slope_parent ) - simulation.add_inertial_contribution( lobe_cur, lobe_parent, slope_parent ) + pfy.flowycpp.mr_lava_loba.perturb_lobe_angle( + lobe_cur, np.linalg.norm(slope_parent), simulation.input, gen + ) + pfy.flowycpp.mr_lava_loba.add_inertial_contribution( + lobe_cur, lobe_parent, np.linalg.norm(slope_parent), simulation.input + ) # Final budding point - angle_diff = lobe_parent.get_azimuthal_angle() - lobe_cur.get_azimuthal_angle() - final_budding_point = lobe_parent.point_at_angle( - angle_diff ) + angle_diff = lobe_parent.get_azimuthal_angle() - lobe_cur.get_azimuthal_angle() + final_budding_point = lobe_parent.point_at_angle(-angle_diff) budding_point_list.append(final_budding_point) - height_budding_point, slope_budding_point = topography.height_and_slope( final_budding_point ) + height_budding_point, slope_budding_point = topography.height_and_slope( + final_budding_point + ) # compute the new lobe axes and get the lobe center and thickness - simulation.compute_lobe_axes( lobe_cur, slope_budding_point ) + pfy.flowycpp.mr_lava_loba.compute_lobe_axes( + lobe_cur, + np.linalg.norm(slope_budding_point), + simulation.input, + common_lobe_dimensions, + ) + # simulation.compute_descendent_lobe_position( lobe_cur, lobe_parent, final_budding_point ) # does not work!! - new_lobe_center = compute_descendent_lobe_position( lobe_cur, lobe_parent, final_budding_point ) - lobe_cur.center = new_lobe_center - lobe_cur.thickness = thickness # Set to constant for now + new_lobe_center = compute_descendent_lobe_position( + lobe_cur, lobe_parent, final_budding_point + ) + lobe_cur.center = new_lobe_center + lobe_cur.thickness = thickness # Set to constant for now lobes.append(lobe_cur) # Representation of the topography plt.pcolormesh(x_data, y_data, topography.height_data.T) -# Plot the lobes +# Plot the lobes colormap = plt.cm.bwr -norm = Normalize(vmin=0, vmax=n_lobes - 1) # Normalize index to range [0, 1] for colormap +norm = Normalize( + vmin=0, vmax=n_lobes - 1 +) # Normalize index to range [0, 1] for colormap scalar_map = ScalarMappable(norm=norm, cmap=colormap) for i in range(n_lobes): color = scalar_map.to_rgba(i) - perimeter = np.array(lobes[i].rasterize_perimeter(input.npoints)) - plt.plot(perimeter[:, 0], perimeter[:, 1], color=color, label=f'Lobe {i}') + + perimeter = np.array( + [ + lobes[i].point_at_angle(phi) + for phi in np.linspace(0, 2 * np.pi, 30, endpoint=True) + ] + ) + plt.plot(perimeter[:, 0], perimeter[:, 1], color=color, label=f"Lobe {i}") plt.plot(lobes[i].center[0], lobes[i].center[1], marker="o", color="black", ms=12) for budding_point in budding_point_list: plt.plot(budding_point[0], budding_point[1], marker=".", color="red", ms=12) plt.gca().set_box_aspect(1) -plt.show() plt.savefig("simulation_test.png", dpi=300) +plt.show() diff --git a/examples/topography_interpolation.py b/examples/topography_interpolation.py index 8a9760b..fd436af 100644 --- a/examples/topography_interpolation.py +++ b/examples/topography_interpolation.py @@ -3,12 +3,15 @@ import matplotlib.pyplot as plt asc_file = pfy.flowycpp.AscFile("topo.asc") -height_data = asc_file.height_data +height_data = asc_file.data x_data = asc_file.x_data y_data = asc_file.y_data topography = pfy.flowycpp.Topography( - asc_file.height_data, asc_file.x_data, asc_file.y_data + asc_file.data, + asc_file.x_data, + asc_file.y_data, + pfy.flowycpp.DEFAULT_NO_DATA_VALUE_HEIGHT, ) x_data_interp = np.linspace(np.min(x_data), np.max(x_data), 500) From ba47a6635add4d5b782f9feedf2d05bf78ecefbb Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 27 Mar 2025 22:00:48 +0000 Subject: [PATCH 6/9] DOC: Update README with more information Co-authored-by: Amrita Goswami --- README.md | 77 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 77 insertions(+) diff --git a/README.md b/README.md index 05a71e6..1c369cd 100644 --- a/README.md +++ b/README.md @@ -23,3 +23,80 @@ git restore subprojects meson setup build --wipe pip install -e . --no-build-isolation ``` + +## How to Run a Simulation + +First, we will briefly describe the building blocks of the simulation, before describing how to run `PyFlowy` from `Python`. +Every simulation consists of flows, which are, in turn made up of lobes. Every lobe changes the topography, and therefore, subsequent flows are affected by previously deposited flows and lobes. + +This beggars the question: why not just have one flow with many lobes? What's the difference between 1 flow with 50 lobes and 10 flows with 5 lobes each? + +- 1 flow with 50 lobes: There is only one lobe with the largest thickness (the thickness of each lobe varies within each flow). +- 10 flows with 5 lobes: Since there are 10 chains of lobes, there are 10 lobes with the largest thickness + +Furthermore (perhaps more crucially; we do not know), parent lobes are only selected within the same flow. Practically, this means that the simulation consisting of a single large flow, will, on average, spread out further from the vent. + +This is enough to start with, to understand how to run `Flowy` and `PyFlowy`. The test `./tests/test_simulation.py` has a working example. + +```python +# A TOML file is parsed which creates an object with several input settings +# An example of the TOML file can be found at this link: +# https://raw.githubusercontent.com/flowy-code/flowy/main/examples/KILAUEA2014-2015/input.toml +input_params = pfy.flowycpp.parse_config( toml_filepath ) +# you can access and write to these settings +input_params.write_thickness_every_n_lobes = 5000 # Writes out the thickness every time 5000 lobes have been processsed. Note: the initial DEM is ignored when writing out intermediate thicknesses. You can still place the intermediate thicknesses on top of the initial DEM + +# These user-defined input settings should be validated (for instance, when the thickening_parameter is set to 1.0 then Flowy will crash. See Issue # 17 for details) +pfy.flowycpp.validate_settings(input_params) +input_params.output_folder = output_folder +# Asc and NetCDF file formats are supported by Flowy +input_params.source = asc_filepath # this example uses an .asc file as input. The topography is constructed from the initial DEM (which is the source here) + +# Outputs can be created in the .asc file format, or the NetCDF file format (generally better), if you have compiled with NetCDF. The NetCDF file format can be compressed, and the user can control the data format in which the output will be saved (only in the case of NetCDF files). Also, some HDF5 viewers can be used to look at NetCDF files! +input_params.output_settings.use_netcdf = True # this uses the NetCDF file format + +# The crop to content functionality lets you "crop" parts of the output which have no lava. This can save a lot of space because DEMs are usually very big and the lava doesn't cover all of it +input_params.output_settings.crop_to_content = True + +# You can save the NetCDF files with data types Short, Float or Double +input_params.output_settings.data_type = pfy.flowycpp.StorageDataType.Short + +# The Simulation class. The first argument is the InputParams class. The second argument is an optional seed. Just set it to None if you don't want to specify it +simulation = pfy.flowycpp.Simulation(input_params, None) +# The run function actually runs the entire simulation (for all flows and lobes) +simulation.run() +``` + +The example above creates and runs a simulation. However, it is also possible to change some input settings in the middle of a simulation and then restart it. Instead of using the `run` function of the simulation class, you would use the `steps` function. The `steps` function takes the number of steps as an argument and returns an enum (`Ongoing` or `Finished`) to indicate whether the simulation was completed or not. Note that the number of steps means the total number of lobes deposited (over all flows). The maximum number of steps that you can perform, therefore, is the product of the number of flows and the number of lobes in each flow. Note that although there are variables for `min_n_lobes` and `max_n_lobes`, in practice, we have found that this is always set to the same value (keeping the number of lobes in each flow constant). + +The following code snippet illustrates this functionality (see also `./tests/test_simulation.py`): + +```python +# Setup for the simulation +input_params = pfy.flowycpp.parse_config( toml_filepath ) +input_params.output_folder = output_folder +input_params.source = asc_filepath +# Create the Simulation object (with a seed of 0) +simulation = pfy.flowycpp.Simulation(input_params, 0) +# Some simulation parameters +simulation.input.n_flows = 10 # Number of flows +simulation.input.min_n_lobes = 50 # in practice the number of lobes is kept constant by keeping min_n_lobes equal to max_n_lobes +simulation.input.max_n_lobes = 50 +# Various options that control the spreading of the lava +simulation.input.max_slope_prob = 0.0 +simulation.input.lobe_exponent = 0.0 +simulation.input.inertial_exponent = 0.0 +simulation.input.thickening_parameter = 0.0 + +# This loop keeps running the simulation in chunks of 100 steps, modifying the max_slop_prob each time +# until the simulation is completed +while simulation.steps(100) == pfy.flowycpp.RunStatus.Ongoing: + simulation.input.max_slope_prob += 0.1 +``` + +## Citation + +https://doi.org/10.48550/arXiv.2405.20144 + + + From 829caa68151bdda2279f3ee62a139de46e849ef8 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 27 Mar 2025 22:01:12 +0000 Subject: [PATCH 7/9] BLD: pin flowy.wrap to new revision --- subprojects/flowy.wrap | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/subprojects/flowy.wrap b/subprojects/flowy.wrap index c3e5a0a..e89db97 100644 --- a/subprojects/flowy.wrap +++ b/subprojects/flowy.wrap @@ -1,5 +1,5 @@ [wrap-git] directory=flowy url=https://github.com/flowy-code/flowy.git -revision=8407b8a5c44bd286836d0f7a234f240cede6aa80 +revision=ba654c77773da309f63b874537c296c847732133 depth=1 \ No newline at end of file From 27051eb8d8b86806ec3d1d402ebfea9867e479d5 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 27 Mar 2025 22:40:28 +0000 Subject: [PATCH 8/9] BLD: new revision for flowy.wrap Co-authored-by: Amrita Goswami --- subprojects/flowy.wrap | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/subprojects/flowy.wrap b/subprojects/flowy.wrap index e89db97..fc03d1f 100644 --- a/subprojects/flowy.wrap +++ b/subprojects/flowy.wrap @@ -1,5 +1,5 @@ [wrap-git] directory=flowy url=https://github.com/flowy-code/flowy.git -revision=ba654c77773da309f63b874537c296c847732133 +revision=6beb3ba141e8ab413e27fc5c57ce4ea1b9b48efd depth=1 \ No newline at end of file From e4c18ae8e72e54485cba7df23d71c446ecb957d1 Mon Sep 17 00:00:00 2001 From: Moritz Sallermann Date: Thu, 27 Mar 2025 22:42:16 +0000 Subject: [PATCH 9/9] Bump version to 1.0.0 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 71c5d9e..02c9f5a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "pyflowy" -version = "0.0.1" +version = "1.0.0" description = "Python bindings for flowy" authors = [ {name = "Amrita Goswami", email = "amrita16thaug646@gmail.com"},