Fastscape c++ coupling [WIP] #5323

Minerallo · 2023-07-14T18:55:17Z

Pull Request Checklist. Please read and check each box with an X. Delete any part not applicable. Ask on the forum if you need help with any step.

The existing coupling between ASPECT and FastScape, developed by Derek Neuharth and Anne Glerum, utilizes the Fortran version of FastScape. However, this Fortran version is no longer being actively developed. I want to couple ASPECT with the actively developed C++ version of FastScape (Fastscapelib), which will incorporate an unstructured grid and be parallelized. This will make ASPECT - FastScape working with large scale spherical models.
I have integrated the Fastscape C++ library to the contrib folder and added the paths to the Cmakelist. I re-use some of the fastscape fortran interface from PR #4396 to develop the initial FastScapecc (.cc, .h) which call the Stream Power law erosion function of fastscapelib.

I added a test case repository in tests repository with a output_screen but it is not working yet. It runs for few timesteps and crashes sending malloc errors.

Before your first pull request:

I have read the guidelines in our CONTRIBUTING.md document.

For all pull requests:

I have followed the instructions for indenting my code.

For new features/models or changes of existing features:

I have tested my new feature locally to ensure it is correct.
I have created a testcase for the new feature/benchmark in the tests/ directory.
I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

gassmoeller

Here are some initial comments, mostly on the header file and variable names.

gassmoeller · 2023-07-14T22:56:12Z

include/aspect/mesh_deformation/fastscapecc.h

@@ -0,0 +1,311 @@
+/*
+  Copyright (C) 2022 by the authors of the ASPECT code.


No, not actually corrected :-)

gassmoeller · 2023-07-14T22:57:39Z

include/aspect/mesh_deformation/fastscapecc.h

+        // /**
+        //  * Destructor for FastScape.
+        //  */
+        // ~FastScape() override;


remove commented lines

To be specific, if you don't need the destructor at all, just remove the whole block here.

gassmoeller · 2023-07-14T22:58:23Z

include/aspect/mesh_deformation/fastscapecc.h

+        void parse_parameters (ParameterHandler &prm);
+
+      private:
+


Unify the number of empty lines, we usually use 1 empty line if you want to separate lines logically. Also check the rest of the header file.

okay thanks corrected

gassmoeller · 2023-07-14T22:59:08Z

include/aspect/mesh_deformation/fastscapecc.h

+        /**
+         * We want to check that every point within the FastScape mesh has been given a value.
+         * However, this is done only on processor 0, so this allows us to send an error later on
+         * if necessary.
+         */
+        // mutable bool fastscape_mesh_filled = true;


Remove commented lines.

gassmoeller · 2023-07-14T22:59:55Z

include/aspect/mesh_deformation/fastscapecc.h

+        /**
+         * How many levels FastScape should be refined above the maximum ASPECT surface resolution.
+         */
+        unsigned int additional_refinement;


Just a suggestion

Suggested change

unsigned int additional_refinement;

unsigned int additional_refinement_levels;

gassmoeller · 2023-07-14T23:10:13Z

include/aspect/mesh_deformation/fastscapecc.h

+         */
+        unsigned int left;
+
+


be consistent with the number of empty lines between parameters

gassmoeller · 2023-07-14T23:11:17Z

include/aspect/mesh_deformation/fastscapecc.h

+         * for the value to be transferred. This is only necessary if use_v is set to 0
+         * and the free surface is used to advect the surface with a normal projection.
+         */
+        double precision;


This could be more self explanatory:

Suggested change

double precision;

double node_tolerance;

gassmoeller · 2023-07-14T23:12:12Z

source/mesh_deformation/fastscapecc.cc

+      // Sub intervals are 3 less than points, if including the ghost nodes. Otherwise 1 less.
+      table_intervals[0] = nx ;
+      table_intervals[dim-1] = 1;
+      table_intervals[1] = ny;


The comment seems to contradict the code

ah yes its corrected, I did not add ghost nodes yet and will see with Benoit how this could be directly done from fastscapelib.

gassmoeller · 2023-07-14T23:12:50Z

source/mesh_deformation/fastscapecc.cc

+      const GeometryModel::Box<dim> *geometry
+        = dynamic_cast<const GeometryModel::Box<dim>*> (&this->get_geometry_model());


As far as I see you require a box geometry? This would be the place to add an Assert to ensure that only box geometries are used.

Yes it is only box geometry for now I am adding back the assert. I will start to have a look how to integrate the use of the unstructured mesh (called trimesh) from fastscape for spherical models this week.

gassmoeller · 2023-07-14T23:13:02Z

source/mesh_deformation/fastscapecc.cc

+      if (this->get_timestep_number() == 0)
+        return;
+
+      // std::vector<std::vector<double>> temporary_variables;


remove commented code

benbovy

@Minerallo I'm very glad to see your effort of coupling fastscapelib (C++) with aspect! Thanks for working in this.

I left a few comments. I hope it will be helpful.

benbovy · 2023-07-15T10:16:02Z

CMakeLists.txt

+INCLUDE_DIRECTORIES("${CMAKE_SOURCE_DIR}/contrib/fastscape/include/")
+INCLUDE_DIRECTORIES("${CMAKE_SOURCE_DIR}/contrib/xtensor/include/")
+INCLUDE_DIRECTORIES("${CMAKE_SOURCE_DIR}/contrib/xtensor/include/xtensor/")
+INCLUDE_DIRECTORIES("${CMAKE_SOURCE_DIR}/contrib/xtl/include/")


I'm not sure how you usually deal with external dependencies in aspect, but I think that using find_package(fastscapelib REQUIRED) then target_link_libraries(${ASPECT_TARGET} INTERFACE fastscapelib) would be a much better alternative. Fastscapelib, xtensor and xtl all provide CMake configuration files when installed so normally finding and linking against fastscapelib should also take care of transitive dependencies (xtensor and xtl).

Or you could use FetchContent or ExternalProject if you want to download those dependencies as part of the aspect's cmake configuration/build.

If for some reason you still want to vendor fastscapelib, xtensor and xtl sources in aspect, please include their license.

Hi @benbovy, thanks for your comments, it's true that it will be a better way of doing it, I'll give it a try, it also makes more sense for the license problem. This was a quick way to get started and will have to be modified.

benbovy · 2023-07-15T10:27:49Z

source/mesh_deformation/fastscapecc.cc

+                            // FastScape indexes from 1 to n, starting at X and Y = 0, and increases
+                            // across the X row. At the end of the row, it jumps back to X = 0
+                            // and up to the next X row in increasing Y direction. We track
+                            // this to correctly place the variables later on.
+                            // Nx*ys effectively tells us what row we are in
+                            // and then indx tells us what position in that row.


I guess this has been copied from #4396?

Beware that there are some major differences between fastscapelib-fortran and fastscapelib (C++). E.g., Fastscapelib-fortran uses a column-major layout (Fortran-style arrays) for representing 2D raster grid nodes on 1D arrays, while fastscapelib (C++) uses a row-major layout (C-style arrays).

yes it was thanks for pointing it out

benbovy · 2023-07-15T10:30:15Z

source/mesh_deformation/fastscapecc.cc

+      for (unsigned int i=0; i<array_size; ++i)
+        {
+          kf[i] = kff;
+          kd[i] = kdd;
+        }


If you only support spatially uniform values for these two parameters in aspect / fastscape coupling, you don't need to create kf and kd vector/array containers. Fastscapelib-fortran explicitly requires arrays but fastscapelib (C++) works with both scalar and array values.

Right I didn't add it in this quick preliminary version of the interface, but these values may change spatially for example when using erosion depending of a prescribed wind direction (example in #4396).

benbovy · 2023-07-15T10:38:42Z

source/mesh_deformation/fastscapecc.cc

+      fastscapelib::raster_boundary_status bs(fastscapelib::node_status::fixed_value);
+
+      auto grid = fastscapelib::raster_grid::from_length({ static_cast<unsigned long>(nx), static_cast<unsigned long>(ny) }, { x_extent, y_extent}, bs);
+
+      // flow graph with single direction flow routing
+      fastscapelib::flow_graph<fastscapelib::raster_grid> flow_graph(
+        grid, { fastscapelib::single_flow_router(), fastscapelib::mst_sink_resolver() });
+      std::cout << "here it works 2"<<std::endl;
+
+      // Setup eroders
+      auto spl_eroder = fastscapelib::make_spl_eroder(flow_graph, 2e-4, 0.4, 1, 1e-5);
+      auto diffusion_eroder = fastscapelib::make_diffusion_adi_eroder(grid, 0.01);


You might want to move the creation of the fastscapelib grid (as well as the flow_graph and eroder objects) to FastScapecc<dim>::initialize. This would be more efficient, especially if FastScapecc<dim>::compute_velocity_constraints_on_boundary is called multiple times and the size of the fastscapelib grid is fixed (as it seems to be the case looking at the values currently set in FastScapecc<dim>::initialize).

Those fastscapelib C++ classes don't provide default constructors, so if you want to have instances of those as members of FastScapecc you will probably need to wrap them in std::unique_ptr.

benbovy · 2023-07-15T10:42:12Z

source/mesh_deformation/fastscapecc.cc

+      //xt::xarray<double> uplift_rate(vec_shape.shape(), 0);
+      //  xt::xarray<double> uplift_rate = xt::adapt(uplift_rate_in_m_year,vec_shape.shape());
+
+      xt::xarray<double> uplift_rate(vec_shape.shape(), 0);


Although I haven't tested it yet, setting both the uplift rate and the (initial) elevation values to zero (vec_shape above) may work but may also make the simulation crash (undefined flow paths compute from the topographic surface).

(I assume you did this for quick testing here but your intent is to eventually set them from other aspect state variables?).

benbovy · 2023-07-15T10:55:35Z

source/mesh_deformation/fastscapecc.cc

+          // apply channel erosion then hillslope diffusion
+          auto spl_erosion = spl_eroder.erode(uplifted_elevation, drainage_area, fastscape_timestep_in_years);
+          std::cout << "here it works 8e "<<std::endl;
+          auto diff_erosion = diffusion_eroder.erode(uplifted_elevation - spl_erosion, fastscape_timestep_in_years);


uplifted_elevation - spl_erosion returns an xtensor (lazy) expression rather than an xtensor container. Maybe this is the reason why the simulation crashes at this point (as I see in the "output_screen" logs)?

I'm actually surprised that it compiles and that it runs fine for a few iterations / time steps (xtensor's lazy expression system is powerful but it might be also tricky and confusing sometimes).

Could you try saving those elevation values in a temporary container before passing it to erode? I.e.,

xt::xarray<double> temp_elevation = uplifted_elevation - spl_erosion; auto diff_erosion = diffusion_eroder.erode(temp_elevation, fastscape_timestep_in_years);

I should probably update the example in fastscapelib documentation too. Creating those temporary containers is not very nice here, better would be to update elevation inplace in each fastscapelib eroder erode() (I plan to add eventually add some overloads to support that).

benbovy · 2023-07-15T11:04:40Z

source/mesh_deformation/fastscapecc.cc

+      std::vector<double> V(array_size);
+
+      // Initialize the variables that will be sent to FastScape.
+      std::vector<double> h(array_size, std::numeric_limits<double>::max());


Is there a reason why you set initial default values to std::numeric_limits<double>::max() for the elevation container?

It doesn't seem like the safest way to initialize it regarding model stability (in case any element is left to its initial value later), but I might be missing something.

Actually, I'm not sure that you even need to create h and h_old if the reason of those two variables was to interface it with fastscapelib-fortran.

Fastscapelib-fortran has its own array of elevation values (among other internal state variables defined in a context), which you normally cannot access directly using its API (you need a copy). With Fastscapelib C++ it is different: you can define your own (xtensor) array container and pass it to the C++ API.

So maybe all you need is an xt::xarray<double> elevation and copy its elements to/from other std::vector<double> variables defined in aspect?

bangerth

It turns out that I have 3 pending comments still. Let me flush them out.

bangerth · 2023-07-15T04:38:05Z

include/aspect/mesh_deformation/fastscapecc.h

+  {
+
+    /**
+     * A plugin that utilizes the landscape evolution code FastScape


Suggested change

* A plugin that utilizes the landscape evolution code FastScape

* A plugin that utilizes the landscape evolution code FastScape++

bangerth · 2023-07-15T04:38:55Z

include/aspect/mesh_deformation/fastscapecc.h

+     *
+     */
+    template<int dim>
+    class FastScapecc : public Interface<dim>, public SimulatorAccess<dim>


How do the FastScape++ people write their software when they can't use anything other than letters? FastScapecc? Or FastScapepp? or FastScapeXX?

bangerth · 2023-07-15T04:41:34Z

include/aspect/mesh_deformation/fastscapecc.h

+        // /**
+        //  * Destructor for FastScape.
+        //  */
+        // ~FastScape() override;


To be specific, if you don't need the destructor at all, just remove the whole block here.

bangerth

@Minerallo There are a number of comments that suggest that you fixed something, but it doesn't show up in github. Did you forget to push these changes? Would you mind doing so?

bangerth · 2023-10-08T00:09:27Z

include/aspect/mesh_deformation/fastscapecc.h

@@ -0,0 +1,311 @@
+/*
+  Copyright (C) 2022 by the authors of the ASPECT code.


No, not actually corrected :-)

bangerth

Well, my time's up. I'll see that I come back to this sometime soon, but perhaps these comments are already useful in some way!

bangerth · 2023-10-08T00:11:32Z

source/mesh_deformation/fastscapecc.cc

@@ -0,0 +1,698 @@
+#include <iostream>


This file is missing its customary header. Just take it from one of the other files, with the copyright year adjusted.

bangerth · 2023-10-08T00:12:40Z

source/mesh_deformation/fastscapecc.cc

+#include <xtensor/xarray.hpp>
+#include <xtensor/xmath.hpp>
+#include <xtensor/xview.hpp>
+#include <xtensor/xadapt.hpp>
+
+#include <fastscapelib/flow/flow_graph.hpp>
+#include <fastscapelib/flow/sink_resolver.hpp>
+#include <fastscapelib/flow/flow_router.hpp>
+#include <fastscapelib/grid/raster_grid.hpp>
+#include <fastscapelib/eroders/diffusion_adi.hpp>
+#include <fastscapelib/eroders/spl.hpp>


We're going to have to figure out a way to disable compilation of this file when fastscape wasn't found during configuration. The way this is typically done is by #includeing global.h and then having a #ifdef ASPECT_WITH_FASTSCAPE_CC or some such.

bangerth · 2023-10-08T00:13:11Z

source/mesh_deformation/fastscapecc.cc

+#include <aspect/geometry_model/box.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <aspect/postprocess/visualization.h>
+#include <ctime>


Move this C++ header include down to where the other C++ headers are (specifically <algorithm>).

bangerth · 2023-10-08T00:13:52Z

source/mesh_deformation/fastscapecc.cc

+    FastScapecc<dim>::initialize ()
+    {
+        AssertThrow(Plugins::plugin_type_matches<const GeometryModel::Box<dim>>(this->get_geometry_model()),
+                         ExcMessage("FastScape can only be run with a box geometry model."));
+
+             const GeometryModel::Box<dim> *geometry
+               = dynamic_cast<const GeometryModel::Box<dim>*> (&this->get_geometry_model());


You will have to run ./contrib/utilities/indent to get a consistent indentation of the entire file. You'll like the end result, it will be much easier to read!

bangerth · 2023-10-08T00:15:13Z

source/mesh_deformation/fastscapecc.cc

+             // Get the deformation type names called for each boundary.
+             std::map<types::boundary_id, std::vector<std::string>> mesh_deformation_boundary_indicators_map
+               = this->get_mesh_deformation_handler().get_active_mesh_deformation_names();
+
+             // Loop over each mesh deformation boundary, and make sure FastScape is only called on the surface.
+             for (std::set<types::boundary_id>::const_iterator p = mesh_deformation_boundary_ids.begin();
+                  p != mesh_deformation_boundary_ids.end(); ++p)
+               {
+                 const std::vector<std::string> names = mesh_deformation_boundary_indicators_map[*p];
+                 for (unsigned int i = 0; i < names.size(); ++i )
+                   {
+                     if (names[i] == "fastscape")
+                       AssertThrow((*p == top_boundary),
+                                   ExcMessage("FastScape can only be called on the surface boundary."));


Suggested change

// Get the deformation type names called for each boundary.

std::map<types::boundary_id, std::vector<std::string>> mesh_deformation_boundary_indicators_map

= this->get_mesh_deformation_handler().get_active_mesh_deformation_names();

// Loop over each mesh deformation boundary, and make sure FastScape is only called on the surface.

for (std::set<types::boundary_id>::const_iterator p = mesh_deformation_boundary_ids.begin();

p != mesh_deformation_boundary_ids.end(); ++p)

{

const std::vector<std::string> names = mesh_deformation_boundary_indicators_map[*p];

for (unsigned int i = 0; i < names.size(); ++i )

{

if (names[i] == "fastscape")

AssertThrow((*p == top_boundary),

ExcMessage("FastScape can only be called on the surface boundary."));

// Get the deformation type names called for each boundary.

const std::map<types::boundary_id, std::vector<std::string>> mesh_deformation_boundary_indicators_map

= this->get_mesh_deformation_handler().get_active_mesh_deformation_names();

// Loop over each mesh deformation boundary, and make sure FastScape is only called on the surface.

for (const auto p : mesh_deformation_boundary_ids)

{

const std::vector<std::string> names = mesh_deformation_boundary_indicators_map[p];

for (unsigned int i = 0; i < names.size(); ++i )

{

if (names[i] == "fastscape")

AssertThrow((p == top_boundary),

ExcMessage("FastScape can only be called on the surface boundary."));

bangerth · 2023-10-08T00:16:11Z

source/mesh_deformation/fastscapecc.cc

+
+
+
+
+
+
+
+
+//
+//      dx = 20000;
+//      dy = dx ;
+//      x_extent = 100e3;
+//      y_extent = 100e3;
+//      nx = 5 ;
+//      ny = 5 ;
+//      array_size = nx*ny;
+//
+//      // Sub intervals are 3 less than points, if including the ghost nodes. Otherwise 1 less.
+//      table_intervals[0] = nx ;
+//      table_intervals[dim-1] = 1;
+//      table_intervals[1] = ny;
+
+//      const GeometryModel::Box<dim> *geometry
+//        = dynamic_cast<const GeometryModel::Box<dim>*> (&this->get_geometry_model());
+//
+//      // The first entry represents the minimum coordinates of the model domain, the second the model extent.
+//      for (unsigned int i=0; i<dim; ++i)
+//        {
+//          grid_extent[i].first = geometry->get_origin()[i];
+//          grid_extent[i].second = geometry->get_extents()[i];
+//        }


If you have code that is no longer necessary, just remove it -- or document why you want to keep it.

bangerth · 2023-10-08T00:16:45Z

source/mesh_deformation/fastscapecc.cc

+         dx = (grid_extent[0].second)/( repetitions[0]);
+
+         // FastScape X extent, which is generally ASPECT's extent unless the ghost nodes are used,
+         // in which case 2 cells are added on either side.
+         x_extent = (grid_extent[0].second) ;


Suggested change

dx = (grid_extent[0].second)/( repetitions[0]);

// FastScape X extent, which is generally ASPECT's extent unless the ghost nodes are used,

// in which case 2 cells are added on either side.

x_extent = (grid_extent[0].second) ;

dx = (grid_extent[0].second)/repetitions[0];

// FastScape X extent, which is generally ASPECT's extent unless the ghost nodes are used,

// in which case 2 cells are added on either side.

x_extent = grid_extent[0].second;

bangerth · 2023-10-08T00:20:20Z

source/mesh_deformation/fastscapecc.cc

+      const types::boundary_id relevant_boundary = this->get_geometry_model().translate_symbolic_boundary_name_to_id ("top");
+      std::vector<std::vector<double>> temporary_variables(dim+2, std::vector<double>());
+
+      // Get a quadrature rule that exists only on the corners, and increase the refinement if specified.


Suggested change

// Get a quadrature rule that exists only on the corners, and increase the refinement if specified.

// Define a quadrature rule that has points only on the vertices of

// a face, then subdivide it a number of times. This way, we can evaluate

// the solution at the surface nodes and points inbetween.

bangerth · 2023-10-08T00:21:49Z

source/mesh_deformation/fastscapecc.cc

+      std::vector<std::vector<double>> temporary_variables(dim+2, std::vector<double>());
+
+      // Get a quadrature rule that exists only on the corners, and increase the refinement if specified.
+      const QIterated<dim-1> face_corners (QTrapez<1>(),


Can you give this variable a better name? Perhaps uniform_face_subdivision?

bangerth · 2023-10-08T00:24:57Z

source/mesh_deformation/fastscapecc.cc

+      for (const auto &cell : this->get_dof_handler().active_cell_iterators())
+        if (cell->is_locally_owned() && cell->at_boundary())
+          for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell; ++face_no)
+            if (cell->face(face_no)->at_boundary())
+              {
+                if ( cell->face(face_no)->boundary_id() != relevant_boundary)
+                  continue;
+
+                std::vector<Tensor<1,dim>> vel(face_corners.size());


This is perhaps a more general comment:

The last line here is not very efficient. You are allocating memory on every face you encounter, but the array you allocate always has the same size. Two lines below you then fill it with some value. At the end of the code block, the memory is then released again.

A better scheme would have declared the variable just before the outer loop, because then the memory is allocated only once. Within the loop, you are then just filling the already allocated array, but the expensive memory allocation and deallocation only happens once in that case.

Minerallo added 5 commits July 8, 2023 11:52

Add fastscapelib xtensor xtl to contrib and Cmakelist

8ce3de4

add interface fastcapecc.cc .h and prm in cookbook

19ad382

replaced hard variables back

14836e0

add test case fastscape_eroding_box_cc.prm

82af2fe

add test repository and output_screen

dc15e7b

gassmoeller reviewed Jul 15, 2023

View reviewed changes

benbovy reviewed Jul 15, 2023

View reviewed changes

correction of V in m/s

9745811

bangerth reviewed Oct 8, 2023

View reviewed changes

		@@ -0,0 +1,311 @@
		/*
		Copyright (C) 2022 by the authors of the ASPECT code.

	unsigned int additional_refinement;
	unsigned int additional_refinement_levels;

		const GeometryModel::Box<dim> *geometry
		= dynamic_cast<const GeometryModel::Box<dim>*> (&this->get_geometry_model());

	* A plugin that utilizes the landscape evolution code FastScape
	* A plugin that utilizes the landscape evolution code FastScape++

-      // Get a quadrature rule that exists only on the corners, and increase the refinement if specified.
+      // Define a quadrature rule that has points only on the vertices of
+      // a face, then subdivide it a number of times. This way, we can evaluate
+      // the solution at the surface nodes and points inbetween.

Fastscape c++ coupling [WIP] #5323

Are you sure you want to change the base?

Fastscape c++ coupling [WIP] #5323

Conversation

Minerallo commented Jul 14, 2023

Before your first pull request:

For all pull requests:

For new features/models or changes of existing features:

gassmoeller left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

benbovy left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Minerallo Jul 24, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bangerth left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bangerth left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bangerth left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Minerallo Jul 24, 2023 •

edited