Stepper func

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
+
+
+

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")

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")

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")