diff --git a/CMakeLists.txt b/CMakeLists.txt index 571808902f1..86d34c5e295 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -231,6 +231,29 @@ if(ASPECT_WITH_FASTSCAPE) endif() +# Fastscapelib (C++) +set(ASPECT_WITH_FASTSCAPELIB OFF CACHE BOOL "Whether the user wants to compile ASPECT with the landscape evolution C++ code Fastscape(lib), or not.") +message(STATUS "Using ASPECT_WITH_FASTSCAPELIB = '${ASPECT_WITH_FASTSCAPE}'") +if(ASPECT_WITH_FASTSCAPELIB) + find_package(fastscapelib CONFIG) + if(${fastscapelib_FOUND}) + set(FS_WITH_HEALPIX OFF CACHE BOOL "Turn off Fastscapelib support of Healpix grid" FORCE) + message(STATUS "Using ASPECT_WITH_FASTSCAPELIB = '${ASPECT_WITH_FASTSCAPELIB}'") + message(STATUS " fastscapelib_INCLUDE_DIR: ${fastscapelib_INCLUDE_DIRS}") + message(STATUS " fastscapelib_VERSION: ${fastscapelib_VERSION}") + # fastscapelib dependencies + find_package(xtensor REQUIRED) + message(STATUS " xtensor_INCLUDE_DIR: ${xtensor_INCLUDE_DIRS}") + message(STATUS " xtensor_VERSION: ${xtensor_VERSION}") + else() + message(STATUS "Fastscapelib not found. Disabling ASPECT_WITH_FASTSCAPELIB. You can specify a hint to your installation directory with fastscapelib_DIR.") + set(ASPECT_WITH_FASTSCAPELIB OFF CACHE BOOL "" FORCE) + endif() +else() + message(STATUS "Using ASPECT_WITH_FASTSCAPELIB = 'OFF'") +endif() + + # NetCDF (c including parallel) set(ASPECT_WITH_NETCDF ON CACHE BOOL "Check if the user wants to compile ASPECT with the NetCDF libraries.") @@ -358,6 +381,7 @@ if(ASPECT_WITH_WORLD_BUILDER) endif() message(STATUS "") +#include "datatypes.h" # # Other stuff about external tools and how we interact with the system: @@ -882,6 +906,15 @@ if (FASTSCAPE) endforeach() endif() +# Fastscapelib +if (${fastscapelib_FOUND}) + message(STATUS "Linking ASPECT against Fastscapelib (header-only)") + foreach(_T ${TARGET_EXECUTABLES}) + target_include_directories(${_T} PRIVATE ${fastscapelib_INCLUDE_DIRS}) + target_include_directories(${_T} PRIVATE ${xtensor_INCLUDE_DIRS}) + endforeach() +endif() + # NetCDF if(${NETCDF_FOUND}) message(STATUS "Linking ASPECT against NetCDF") diff --git a/WITH_FASTSCAPELIB.md b/WITH_FASTSCAPELIB.md new file mode 100644 index 00000000000..5d15216f0bf --- /dev/null +++ b/WITH_FASTSCAPELIB.md @@ -0,0 +1,117 @@ +# Using Aspect with Fastscapelib (C++) + +## Install Fastscapelib using conda + +You can use conda to install fastscapelib and all of its dependencies (we assume +that you are familiar with managing conda environments): + +``` +conda install fastscapelib -c conda-forge +``` + +## Install Fastscapelib from source + +You can follow the steps below to install Fastscapelib and its dependencies from +source. + +### Set the installation prefix path for Fastscapelib and its dependencies + +Choose a common location where you want to download the source files of +Fastscapelib and its dependencies: + +``` +$ export FS_SOURCE_PATH=/path/to/fastscapelib/source +``` + +Choose a common location where you want to install Fastscapelib and its +dependencies: + +``` +export CMAKE_INSTALL_PREFIX=/path/to/fastscapelib/install/dir +``` + +### Install xtl and xtensor + +Download the `xtl` git repository and checkout the last stable version: + +``` +$ cd $FS_SOURCE_PATH +$ git clone https://github.com/xtensor-stack/xtl +$ cd xtl +$ git checkout 0.7.7 +``` + +Build and install `xtl`: + +``` +$ cmake -S. -Bbuild +$ cmake --install build +``` + +Download the `xtensor` git repository and checkout the last stable version: + +``` +$ cd $FS_SOURCE_PATH +$ git clone https://github.com/xtensor-stack/xtensor +$ cd xtensor +$ git checkout 0.25.0 +``` + +Build and install `xtensor`: + +``` +$ cmake -S. -Bbuild +$ cmake --install build +``` + +### Install Fastscapelib + +Download the `fastscapelib` git repository and checkout the last stable version: + +``` +$ cd $FS_SOURCE_PATH +$ git clone https://github.com/fastscape-lem/fastscapelib +$ cd fastscapelib +$ git checkout v0.2.2 +``` + +IMPORTANT: edit the `cmake/fastscapelibConfig.cmake` file with the editor of +your choice and add the following line just below `@PACKAGE_INIT@` (this will be +required every time you checkout another version or branch and re-install +Fastscapelib, but this won't be required for Fastscapelib version >0.2.2): + +``` +include(CMakeFindDependencyMacro) +``` + +Build and install `fastscapelib`: + +``` +$ cmake -S. -Bbuild +$ cmake --install build +``` + +## Compile Aspect with Fastscapelib enabled + +If you have installed Fastscapelib using conda, don't forget to activate the +environment where you have installed Fastscapelib. Once activated, configuring +`aspect` with CMake should then find Fastscapelib automatically. + +If you have installed Fastscapelib from source like described above, you can +help CMake find the Fastscapelib installation path using: + +``` +export CMAKE_PREFIX_PATH=/path/to/fastscapelib/install/dir +``` + +where the given path corresponds to the path that you set during installation +via the `CMAKE_INSTALL_PREFIX` environment variable (see instructions above). + +Configuring `aspect` should then find Fastscapelib. Check the output of the +following command run from aspect's source root directory (note: you might need +to set additional arguments or run extra commands in order to find the other +dependencies of `aspect`, we skipped it here): + +``` +cmake -S. -Bbuild -DASPECT_WITH_FASTSCAPELIB=ON +``` diff --git a/contrib/fastscape/fastscape_rewrite_VTK_with_time.py b/contrib/fastscape/fastscape_rewrite_VTK_with_time.py index ffde75c5bb6..f39d656de85 100644 --- a/contrib/fastscape/fastscape_rewrite_VTK_with_time.py +++ b/contrib/fastscape/fastscape_rewrite_VTK_with_time.py @@ -33,7 +33,7 @@ def get_times_pvd(filename): times[i] = float(temp_str) - return times + return times #%% Get file paths (absolute, not relative paths!) @@ -54,73 +54,73 @@ def absoluteFilePaths(directory): # trace generated using paraview version 5.8.0 # -# To ensure correct image size when batch processing, please search +# To ensure correct image size when batch processing, please search # for and uncomment the line `# renderView*.ViewSize = [*,*]` for File in FileNames: #### disable automatic camera reset on 'Show' paraview.simple._DisableFirstRenderCameraReset() - + # create a new 'Legacy VTK Reader' topography0000 = LegacyVTKReader(FileNames=File) - + # get animation scene animationScene1 = GetAnimationScene() - + # get the time-keeper timeKeeper1 = GetTimeKeeper() - + # update animation scene based on data timesteps animationScene1.UpdateAnimationUsingDataTimeSteps() - + # get active view renderView1 = GetActiveViewOrCreate('RenderView') # uncomment following to set a specific view size # renderView1.ViewSize = [882, 554] - + # get layout layout1 = GetLayout() - + # show data in view topography0000Display = Show(topography0000, renderView1, 'StructuredGridRepresentation') - + # trace defaults for the display properties. topography0000Display.Representation = 'Surface' - + # reset view to fit data renderView1.ResetCamera() - + # get the material library materialLibrary1 = GetMaterialLibrary() - + # show color bar/color legend topography0000Display.SetScalarBarVisibility(renderView1, True) - + # update the view to ensure updated data information renderView1.Update() - + # get color transfer function/color map for 'H' hLUT = GetColorTransferFunction('H') - + # get opacity transfer function/opacity map for 'H' hPWF = GetOpacityTransferFunction('H') - + # save data SaveData(File[:-4] + '.vts', proxy=topography0000, ChooseArraysToWrite=1, PointDataArrays=['HHHHH','basement','catchment','drainage_area','erosion_rate','topography','total_erosion']) - + #### saving camera placements for all active views - + # current camera placement for renderView1 renderView1.CameraPosition = [225625.0, 25625.0, 878497.0779980461] renderView1.CameraFocalPoint = [225625.0, 25625.0, 1135.7277145385742] renderView1.CameraParallelScale = 227077.82689023562 - + #### uncomment the following to render all views # RenderAllViews() # alternatively, if you want to write images, you can use SaveScreenshot(...). - + Delete(topography0000) del topography0000 @@ -140,17 +140,9 @@ def absoluteFilePaths(directory): for n, File in enumerate(FileNames): dataset_attributes = {"timestep": str(times[n]), "group": "", "group": "", "file": os.path.basename(File)[:-4]+'.vts'} - + ET.SubElement(doc, "DataSet", attrib=dataset_attributes).text="\n" tree = ET.ElementTree(root) -tree.write('./VTK/topography.pvd', xml_declaration=True, encoding='utf-8', method="xml") - - - - - - - - +tree.write('./VTK/topography.pvd', xml_declaration=True, encoding='utf-8', method="xml") diff --git a/include/aspect/mesh_deformation/fastscapecc.h b/include/aspect/mesh_deformation/fastscapecc.h new file mode 100644 index 00000000000..6cbd3d6fbee --- /dev/null +++ b/include/aspect/mesh_deformation/fastscapecc.h @@ -0,0 +1,267 @@ +/* + Copyright (C) 2023 by the authors of the ASPECT code. + This file is part of ASPECT. + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + +#ifndef _aspect_mesh_deformation_fastscapecc_h +#define _aspect_mesh_deformation_fastscapecc_h + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +namespace aspect +{ + using namespace dealii; + + namespace MeshDeformation + { + + /** + * A plugin that utilizes the landscape evolution code Fastscapelib (C++) + * to deform the ASPECT boundary through erosion. + * + */ + template + class FastScapecc : public Interface, public SimulatorAccess + { + public: + FastScapecc(); + + /** + * Initialize variables for FastScape. + */ + void initialize () override; + + /** + * A function that creates constraints for the velocity of certain mesh + * vertices (e.g. the surface vertices) for a specific boundary. + * The calling class will respect + * these constraints when computing the new vertex positions. + */ + void + compute_velocity_constraints_on_boundary(const DoFHandler &mesh_deformation_dof_handler, + AffineConstraints &mesh_velocity_constraints, + const std::set &boundary_id) const override; + + + /** + * Returns whether or not the plugin requires surface stabilization + */ + bool needs_surface_stabilization () const override; + + /** + * Declare parameters for the FastScape plugin. + */ + static + void declare_parameters (ParameterHandler &prm); + + /** + * Parse parameters for the FastScape plugin. + */ + void parse_parameters (ParameterHandler &prm) override; + + private: + /** + * Surface mesh and solution. + */ + using SurfaceMeshType = Triangulation; + + SurfaceMeshType surface_mesh; + DoFHandler surface_mesh_dof_handler; + LinearAlgebra::Vector surface_solution; + mutable AffineConstraints surface_constraints; // Constraints for hanging nodes + dealii::LinearAlgebra::distributed::Vector boundary_solution; + + template + void init_surface_mesh(M &geom_model); + + template >::value>> + void init_surface_mesh(M &geom_model); + + template >::value>> + void init_surface_mesh(M &geom_model); + + /** + * Pointers to Fastscapelib objects + */ + using GridAdapterType = typename fastscapelib::dealii_grid; + using FlowGraphType = typename fastscapelib::flow_graph; + std::unique_ptr grid; + std::unique_ptr flow_graph; + std::unique_ptr> spl_eroder; + + void project_surface_solution(const std::set &boundary_ids); + + /** + * Project the Stokes velocity solution onto the + * free surface. Called by make_constraints() + */ + // void project_velocity_onto_boundary (const DoFHandler &free_surface_dof_handler, + // const IndexSet &mesh_locally_owned, + // const IndexSet &mesh_locally_relevant, + // LinearAlgebra::Vector &output) const; + + + /** + * Fill velocity data table to be interpolated back onto the ASPECT mesh. + */ + Table fill_data_table(std::vector &values, + TableIndices &size_idx, + const int &array_size) const; + + /** + * Suggestion for the number of FastScape steps to run for every ASPECT timestep, + * where the FastScape timestep is determined by ASPECT_timestep_length divided by + * this parameter. + */ + unsigned int fastscape_steps_per_aspect_step; + + /** + * Maximum timestep allowed for FastScape, if the suggested timestep exceeds this + * limit it is repeatedly divided by 2 until the final timestep is smaller than this parameter. + */ + double maximum_fastscape_timestep; + + /** + * Maximum timestep allowed for FastScape, if the suggested timestep exceeds this + * limit it is repeatedly divided by 2 until the final timestep is smaller than this parameter. + */ + + /** + * Check whether FastScape needs to be restarted. + */ + mutable bool restart; + + /** + * ASPECT end time to check if we should destroy FastScape. + */ + double end_time; + + /** + * Total number of Fastscapelib grid nodes. + */ + unsigned int n_grid_nodes; + + /** + * How many levels FastScape should be refined above the maximum ASPECT surface resolution. + */ + unsigned int additional_refinement_levels; + + /** + * Maximum expected refinement level at ASPECT's surface. + * This and resolution_difference are required to properly transfer node data from + * ASPECT to FastScape. + */ + int maximum_surface_refinement_level; + + /** + * Difference in refinement levels expected at the ASPECT surface, + * where this would be set to 2 if 3 refinement leves are set at the surface. + * This and surface_resolution are required to properly transfer node data from + * ASPECT to FastScape. + * + * TODO: Should this be kept this way, or make it so the input is the expected levels + * of refinement at the surface, and we can subtract one within the code? Also, + * it would be good to find a way to check these are correct, because they are a + * common source of errors. + */ + int surface_refinement_difference; + + /** + * Seed number for initial topography noise in FastScape. + */ + int fs_seed; + + /** + * Table for interpolating FastScape surface velocities back to ASPECT. + */ + std::array< unsigned int, dim > table_intervals; + + /** + * Magnitude (m) of the initial noise applied to FastScape. + * Applied as either a + or - value to the topography + * such that the total difference can be up to 2*noise_h. + */ + double noise_h; + + // FastScape erosional parameters // + /** + * Drainage area exponent for the stream power law (m variable in FastScape surface equation). + */ + double m; + + /** + * Slope exponent for the steam power law (n variable in FastScape surface equation). + */ + double n; + + /** + * Slope exponent for multi-direction flow, where 0 is uniform, and 10 is steepest descent. (-1 varies with slope) + * (p variable in FastScape surface equation). + */ + double p; + + /** + * Bedrock river incision rate for the stream power law + * (meters^(1-2m)/yr, kf variable in FastScape surface equation). + */ + double kff; + + /** + * Sediment river incision rate for the stream power law (meters^(1-2m)/yr). + * When set to -1 this is identical to the bedrock value. + * (kf variable in FastScape surface equation applied to sediment). + */ + double kfsed; + + /** + * Bedrock transport coefficient for hillslope diffusion (m^2/yr, kd in FastScape surface equation). + */ + double kdd; + + /** + * Bedrock transport coefficient for hillslope diffusion (m^2/yr). When set to -1 this is + * identical to the bedrock value. + * (kd in FastScape surface equation applied to sediment). + */ + double kdsed; + + /** + * Precision value for how close a ASPECT node must be to the FastScape node + * 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 node_tolerance; + + /** + * FastScape X extent (ASPECT X extent + 2*dx for ghost nodes). + */ + double precision; + }; + } +} + +#endif diff --git a/include/aspect/mesh_deformation/fastscapecc_adapter.h b/include/aspect/mesh_deformation/fastscapecc_adapter.h new file mode 100644 index 00000000000..0e2023d9bbf --- /dev/null +++ b/include/aspect/mesh_deformation/fastscapecc_adapter.h @@ -0,0 +1,234 @@ +#ifndef FASTSCAPELIB_DEALII_GRID_H_ +#define FASTSCAPELIB_DEALII_GRID_H_ + +#include +#include + +#include +#include + +#include "fastscapelib/grid/base.hpp" +#include "fastscapelib/utils/xtensor_utils.hpp" + + +namespace fastscapelib +{ + template + class dealii_grid; + + /** + * Dealii grid specialized types + */ + template + struct grid_inner_types> + { + static constexpr bool is_structured = false; + static constexpr bool is_uniform = false; + + using grid_data_type = double; + + using xt_selector = S; + static constexpr std::size_t xt_ndims = 1; + + static constexpr uint8_t n_neighbors_max = 8; + using neighbors_cache_type = neighbors_no_cache; + }; + + /** + * @brief Grid adapter class for using a DEAL.II surface or boundary mesh with + * Fastscapelib + * + * This adapter should work with various kinds of DEAL.II meshes (e.g., box, + * sphere, etc.). The only requirement is that the mesh must be uniform (i.e., + * all cells have an equal area) at a given refinement level. + * + * Note: a ``dealii_grid`` node cooresponds to a DEAL.II mesh cell. + * + * @tparam T The DEAL.II mesh type. + * @tparam S The container selector for data array members. + */ + template + class dealii_grid : public grid> + { + public: + using self_type = dealii_grid; + using base_type = grid; + using inner_types = grid_inner_types; + + using grid_data_type = typename base_type::grid_data_type; + + using xt_selector = typename base_type::xt_selector; + using container_type = xt_tensor_t; + using neighbors_type = typename base_type::neighbors_type; + using neighbors_indices_type = typename base_type::neighbors_indices_type; + using neighbors_distances_type = typename base_type::neighbors_distances_type; + + using nodes_status_type = typename base_type::nodes_status_type; + using nodes_status_array_type = xt_tensor_t; + + using size_type = typename base_type::size_type; + using shape_type = typename base_type::shape_type; + + dealii_grid(T &triangulation, double cell_area); + + void set_nodes_status(const nodes_status_array_type &nodes_status); + + protected: + T &m_triangulation; + + shape_type m_shape; + size_type m_size; + double m_node_area; + + nodes_status_type m_nodes_status; + + using neighbors_distances_impl_type = typename base_type::neighbors_distances_impl_type; + using neighbors_indices_impl_type = typename base_type::neighbors_indices_impl_type; + + std::vector m_neighbors_count; + std::vector m_neighbors_indices; + std::vector m_neighbors_distances; + + void compute_connectivity(); + + inline container_type nodes_areas_impl() const; + inline grid_data_type nodes_areas_impl(const size_type &idx) const noexcept; + + inline size_type neighbors_count_impl(const size_type &idx) const; + + void neighbors_indices_impl(neighbors_indices_impl_type &neighbors, + const size_type &idx) const; + + inline const neighbors_distances_impl_type &neighbors_distances_impl( + const size_type &idx) const; + + static constexpr std::size_t dimension_impl() noexcept; + + friend class grid; + }; + + + /** + * Creates a new adapter from an existing DEAL.II surface mesh. + * + * @param triangulation The DEAL.II triangulation object. + * @param cell_area the uniform area of the mesh cells. + */ + template + dealii_grid::dealii_grid(T &triangulation, double cell_area) + : base_type(0) + , m_triangulation(triangulation) + , m_node_area(cell_area) + { + + m_size = triangulation.n_global_active_cells(); + m_shape = { static_cast(m_size) }; + + // pre-compute explicit grid connectivity + compute_connectivity(); + + // Set all grid nodes as "active" + // - Fastscapelib doesn't support yet labelling grid nodes as ghost nodes + // - we assume that all the erosion processes that we apply here depend on flow routing + // (i.e., they rely on a `fastscapelib::flow_graph` object built on top of this grid, + // such as `fastscapelib::spl_eroder`) + // - It is better (more explicit) to set base level nodes via the flow graph + // object (i.e., `fastscapelib::flow_graph::set_base_levels`) + // - Masks (e.g., oceans) can be defined via the flow graph too + // (i.e., `fastscapelib::flow_graph::set_mask`) + // - TODO: periodic boundary conditions (if any) should be handled in `compute_connectivity` + // (there is no need to label nodes as periodic boundaries since the grid connectivity + // is computed explicitly here) + m_nodes_status = node_status::core; + } + + template + void dealii_grid::set_nodes_status(const nodes_status_array_type &nodes_status) + { + if (!xt::same_shape(nodes_status.shape(), m_shape)) + { + throw std::invalid_argument( + "invalid shape for nodes_status array (expects shape [N] where N is the total number of nodes)"); + } + m_nodes_status = nodes_status; + } + + template + void dealii_grid::compute_connectivity() + { + m_neighbors_count.resize(m_size); + m_neighbors_indices.resize(m_size); + m_neighbors_distances.resize(m_size); + unsigned int counter = 0; + std::unordered_map global_to_local_index; + for (const auto &cell : m_triangulation->get_dof_handler().active_cell_iterators()) + { + m_neighbors_count[counter] = cell->n_faces(); + auto center = cell.center(); + for (auto face_i = 0; face_i < m_neighbors_count[counter]; ++face_i) + { + unsigned int global_neighbor_index = 0; + for (const auto &neighbor_cell : cell.neighbor(face_i)) + { + global_neighbor_index = neighbor_cell.global_active_cell_index(); + m_neighbors_distances[counter][face_i] = center.distance(neighbor_cell.center()); + // TODO: I assume we only have one neighbor per face, since we start with a uniform grid, we will only have one, but cwe could generatlize this + } + // TODO: add assert to check if neighbor_index is correclty set + global_to_local_index.emplace(global_neighbor_index, counter); + } + + counter++; + } + counter = 0; + for (const auto &cell : this->get_dof_handler().active_cell_iterators()) + { + for (auto k = 0; k <= m_neighbors_count[counter]; k++) + { + auto global_nb_index = m_neighbors_indices[counter][k]; + m_neighbors_indices[counter][k] = global_to_local_index.at(global_nb_index); + } + counter++; + } + } + + template + inline auto dealii_grid::nodes_areas_impl() const -> container_type + { + return xt::broadcast(m_node_area, m_shape); + } + + template + inline auto dealii_grid::nodes_areas_impl(const size_type & /*idx*/) const noexcept + -> grid_data_type + { + return m_node_area; + } + + template + inline auto dealii_grid::neighbors_count_impl(const size_type &idx) const -> size_type + { + return m_neighbors_count[idx]; + } + + template + void dealii_grid::neighbors_indices_impl(neighbors_indices_impl_type &neighbors, + const size_type &idx) const + { + const auto &size = m_neighbors_count[idx]; + + for (size_type i = 0; i < size; i++) + { + neighbors[i] = m_neighbors_indices[idx][i]; + } + } + + template + auto dealii_grid::neighbors_distances_impl(const size_type &idx) const + -> const neighbors_distances_impl_type & + { + return m_neighbors_distances[idx]; + } +} + +#endif // FASTSCAPELIB_DEALII_GRID_H__ diff --git a/source/mesh_deformation/fastscapecc.cc b/source/mesh_deformation/fastscapecc.cc new file mode 100644 index 00000000000..fb94642dc3a --- /dev/null +++ b/source/mesh_deformation/fastscapecc.cc @@ -0,0 +1,537 @@ +/* + Copyright (C) 2011 - 2022 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . + */ + + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + + +namespace aspect +{ + namespace MeshDeformation + { + template + FastScapecc::FastScapecc() + : surface_mesh_dof_handler(surface_mesh) // Link DoFHandler to surface_mesh + {} + + template + void FastScapecc::initialize() + { + init_surface_mesh(this->get_geometry_model()); + n_grid_nodes = surface_mesh.n_active_cells(); + } + + template + template + void FastScapecc::init_surface_mesh(M & /*geom_model*/) + { + AssertThrow(false, ExcMessage("FastScapecc plugin only supports 3D Box or Spherical Shell geometries.")); + } + + template + template >::value>> + void FastScapecc::init_surface_mesh(M &geom_model) + { + + this->get_pcout() << "Box geometry detected. Initializing FastScape for Box geometry..." << std::endl; + + // Create the rectangle surface mesh + const auto origin = geom_model->get_origin(); + const auto extent = geom_model->get_extents(); + + const Point<2> p1(origin[0], origin[1]); + const Point<2> p2(origin[0] + extent[0], origin[1] + extent[1]); + + auto rep = geom_model->get_repetitions(); + std::vector repetitions(rep.begin(), rep.end()); + + GridGenerator::subdivided_hyper_rectangle(surface_mesh, repetitions, p1, p2); + } + + template + template >::value>> + void FastScapecc::init_surface_mesh(M &geom_model) + { + + this->get_pcout() << "Spherical Shell geometry detected. Initializing FastScape for Spherical Shell geometry..." << std::endl; + + // Create the spherical surface mesh + const Point<3> center(0, 0, 0); // Center at the origin + GridGenerator::hyper_sphere(surface_mesh, center, geom_model->outer_radius()); + surface_mesh.refine_global(3); + + DoFTools::make_hanging_node_constraints(surface_mesh_dof_handler, surface_constraints); + surface_constraints.close(); + } + + template + void FastScapecc::project_surface_solution(const std::set &boundary_ids) + { + TimerOutput::Scope timer_section(this->get_computing_timer(), "Project surface solution"); + + // Define the quadrature rule for projection + QGauss<2> quadrature(surface_mesh_dof_handler.get_fe().degree + 1); + + // Initialize the surface solution vector + IndexSet locally_relevant_dofs_surface; + DoFTools::extract_locally_relevant_dofs(surface_mesh_dof_handler, locally_relevant_dofs_surface); + surface_solution.reinit(surface_mesh_dof_handler.locally_owned_dofs(), + locally_relevant_dofs_surface, + this->get_mpi_communicator()); + + // Initialize the boundary solution vector + IndexSet locally_relevant_dofs_main; + DoFTools::extract_locally_relevant_dofs(this->get_dof_handler(), locally_relevant_dofs_main); + boundary_solution.reinit(this->get_dof_handler().locally_owned_dofs(), + locally_relevant_dofs_main, + this->get_mpi_communicator()); + + const types::boundary_id relevant_boundary = this->get_geometry_model().translate_symbolic_boundary_name_to_id("top"); + + for (const auto &cell : this->get_dof_handler().active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary() && + cell->face(face)->boundary_id() == relevant_boundary) + { + for (unsigned int vertex = 0; vertex < GeometryInfo::vertices_per_face; ++vertex) + { + const unsigned int dof_index = cell->face(face)->vertex_dof_index(vertex, 0); + boundary_solution[dof_index] = this->get_solution()[dof_index]; + } + } + } + } + } + + boundary_solution.compress(VectorOperation::insert); + + // dealii::Functions::FEFieldFunction<2, dealii::LinearAlgebra::distributed::Vector,3> + // fe_function(this->get_dof_handler(), boundary_solution, this->get_mapping()); + + // VectorTools::interpolate_boundary_values (surface_mesh_dof_handler, + // *boundary_ids.begin(), + // boundary_solution, + // surface_constraints); + + // VectorTools::interpolate_to_different_mesh( + // this->get_dof_handler(), // Source DoFHandler + // this->get_solution(), // Source solution vector + // surface_mesh_dof_handler, // Target DoFHandler + // surface_solution // Target solution vector + // ); + + // Project the boundary solution onto the 2D surface mesh + // VectorTools::project( + // this->template get_mapping<2,3>(), + // surface_mesh_dof_handler, + // surface_constraints, + // quadrature, + // fe_function, + // surface_solution); + } + + + + template + void FastScapecc::compute_velocity_constraints_on_boundary(const DoFHandler &mesh_deformation_dof_handler, + AffineConstraints &mesh_velocity_constraints, + const std::set &boundary_ids) const + { + if (this->get_timestep_number() == 0) + return; + + TimerOutput::Scope timer_section(this->get_computing_timer(), "FastScape plugin"); + + const_cast *>(this)->project_surface_solution(boundary_ids); + + // VectorTools::interpolate_to_different_mesh( + // this->get_dof_handler(), // Source DoFHandler + // this->get_solution(), // Source solution vector + // surface_mesh_dof_handler, // Target DoFHandler + // surface_solution // Target solution vector + // ); + + // Apply hanging node constraints to ensure continuity + // TODO: this yields a link error + // surface_constraints.distribute(surface_solution); + + // Temporary storage for computation + std::vector> temporary_variables(dim + 2, std::vector()); + + // Iterate over active cells of the surface mesh + for (const auto &cell : surface_mesh_dof_handler.active_cell_iterators()) + { + for (unsigned int vertex_index = 0; vertex_index < GeometryInfo<2>::vertices_per_cell; ++vertex_index) + { + const Point vertex = cell->vertex(vertex_index); // Access the vertex + + // Access the DoF index for the current vertex + unsigned int vertex_dof_index = cell->vertex_dof_index(vertex_index, 0); // Assume component 0 + + // Retrieve the solution value at the DoF index + double interpolated_value = surface_solution[vertex_dof_index]; + + // Compute radial velocity + double radial_velocity = (vertex(0) * interpolated_value + + vertex(1) * interpolated_value + + vertex(2) * interpolated_value) / + vertex.norm(); + + // Surface elevation + double elevation = vertex(2); + + auto spherical_model = dynamic_cast*>(&this->get_geometry_model()); + if (spherical_model != nullptr) + { + elevation -= spherical_model->outer_radius(); + } + + // Store results + temporary_variables[0].push_back(elevation); + temporary_variables[2].push_back(radial_velocity * year_in_seconds); + } + } + + // int array_size = healpix_grid.Npix(); + std::vector V(n_grid_nodes); + + if (Utilities::MPI::this_mpi_process(this->get_mpi_communicator()) == 0) + { + // Initialize the variables that will be sent to FastScape. + std::vector h(n_grid_nodes, std::numeric_limits::max()); + std::vector vz(n_grid_nodes); + std::vector h_old(n_grid_nodes); + + for (unsigned int i = 0; i < temporary_variables[1].size(); ++i) + { + int index = static_cast(temporary_variables[1][i]); + h[index] = temporary_variables[0][i]; + vz[index] = temporary_variables[2][i]; + } + + for (unsigned int p = 1; p < Utilities::MPI::n_mpi_processes(this->get_mpi_communicator()); ++p) + { + MPI_Status status; + MPI_Probe(p, 42, this->get_mpi_communicator(), &status); + int incoming_size = 0; + MPI_Get_count(&status, MPI_DOUBLE, &incoming_size); + + for (unsigned int i = 0; i < temporary_variables.size(); ++i) + { + temporary_variables[i].resize(incoming_size); + } + + for (unsigned int i = 0; i < temporary_variables.size(); ++i) + MPI_Recv(&temporary_variables[i][0], incoming_size, MPI_DOUBLE, p, 42, this->get_mpi_communicator(), &status); + + for (unsigned int i = 0; i < temporary_variables[1].size(); ++i) + { + int index = static_cast(temporary_variables[1][i]); + h[index] = temporary_variables[0][i]; + vz[index] = temporary_variables[2][i]; + } + } + + for (unsigned int i = 0; i < n_grid_nodes; ++i) + { + h_old[i] = h[i]; + } + + const double aspect_timestep_in_years = this->get_timestep() / year_in_seconds; + + unsigned int fastscape_iterations = fastscape_steps_per_aspect_step; + double fastscape_timestep_in_years = aspect_timestep_in_years / fastscape_iterations; + while (fastscape_timestep_in_years > maximum_fastscape_timestep) + { + fastscape_iterations *= 2; + fastscape_timestep_in_years *= 0.5; + } + + // xt::xarray node_status_array = xt::zeros({ array_size }); + // // TODO: replace Healpix grid by dealii grid adapter + // auto grid = fastscapelib::healpix_grid<>(nsides, node_status_array, 6.371e6); + // auto flow_graph = fastscapelib::flow_graph>( + // grid, + // { + // fastscapelib::single_flow_router() + // } + // ); + // auto spl_eroder = fastscapelib::make_spl_eroder(flow_graph, 2e-4, 0.4, 1, 1e-5); + + // xt::xarray uplifted_elevation = xt::zeros(grid.shape()); + // xt::xarray drainage_area = xt::zeros(grid.shape()); + // xt::xarray sediment_flux = xt::zeros(grid.shape()); + + // std::vector shape = { static_cast(array_size) }; + // auto uplift_rate = xt::adapt(vz, shape); + // auto elevation = xt::adapt(h, shape); + // auto elevation_old = xt::adapt(h_old, shape); + + // for (unsigned int fastscape_iteration = 0; fastscape_iteration < fastscape_iterations; ++fastscape_iteration) + // { + // uplifted_elevation = elevation + fastscape_timestep_in_years * uplift_rate; + // flow_graph.update_routes(uplifted_elevation); + // flow_graph.accumulate(drainage_area, 1.0); + // auto spl_erosion = spl_eroder.erode(uplifted_elevation, drainage_area, fastscape_timestep_in_years); + // sediment_flux = flow_graph.accumulate(spl_erosion); + // elevation = uplifted_elevation - spl_erosion; + // } + + // std::vector elevation_std(elevation.begin(), elevation.end()); + // std::vector elevation_old_std(elevation_old.begin(), elevation.end()); + + // for (unsigned int i = 0; i < array_size; ++i) + // { + // V[i] = (elevation_std[i] - elevation_old_std[i]) / (this->get_timestep() / year_in_seconds); + // } + // MPI_Bcast(&V[0], array_size, MPI_DOUBLE, 0, this->get_mpi_communicator()); + } + else + { + for (unsigned int i = 0; i < temporary_variables.size(); ++i) + MPI_Ssend(&temporary_variables[i][0], temporary_variables[1].size(), MPI_DOUBLE, 0, 42, this->get_mpi_communicator()); + + MPI_Bcast(&V[0], n_grid_nodes, MPI_DOUBLE, 0, this->get_mpi_communicator()); + } + + // auto healpix_velocity_function = [&](const Point &p) -> double + // { + // int index = healpix_grid.vec2pix({p(0), p(1), p(2)}); + // return V[index]; + // }; + + // VectorFunctionFromScalarFunctionObject vector_function_object( + // healpix_velocity_function, + // dim - 1, + // dim); + + // VectorTools::interpolate_boundary_values (mesh_deformation_dof_handler, + // *boundary_ids.begin(), + // vector_function_object, + // mesh_velocity_constraints); + } + + + template + Table + FastScapecc::fill_data_table(std::vector &values, + TableIndices &size_idx, + const int &array_size) const + { + // Create data table based off of the given size. + Table data_table; + data_table.TableBase::reinit(size_idx); + TableIndices idx; + + // Loop through the data table and fill it with the velocities from FastScape. + + // Indexes through z, y, and then x. + for (unsigned int k=0; k + bool + FastScapecc:: + needs_surface_stabilization () const + { + return true; + } + + + template + void FastScapecc::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection ("Mesh deformation"); + { + prm.enter_subsection ("Fastscapecc"); + { + prm.declare_entry("Number of steps", "10", + Patterns::Integer(), + "Number of steps per ASPECT timestep"); + prm.declare_entry("Maximum timestep", "10e3", + Patterns::Double(0), + "Maximum timestep for FastScape. Units: $\\{yrs}$"); + prm.declare_entry("Additional fastscape refinement levels", "0", + Patterns::Integer(), + "How many levels above ASPECT FastScape should be refined."); + prm.declare_entry("Fastscape seed", "1000", + Patterns::Integer(), + "Seed used for adding an initial noise to FastScape topography based on the initial noise magnitude."); + prm.declare_entry("Maximum surface refinement level", "1", + Patterns::Integer(), + "This should be set to the highest ASPECT refinement level expected at the surface."); + prm.declare_entry("Surface refinement difference", "0", + Patterns::Integer(), + "The difference between the lowest and highest refinement level at the surface. E.g., if three resolution " + "levels are expected, this would be set to 2."); + prm.declare_entry ("Use velocities", "true", + Patterns::Bool (), + "Flag to use FastScape advection and uplift."); + prm.declare_entry("Precision", "0.001", + Patterns::Double(), + "Precision value for how close a ASPECT node must be to the FastScape node for the value to be transferred."); + prm.declare_entry("Initial noise magnitude", "5", + Patterns::Double(), + "Maximum topography change from the initial noise. Units: $\\{m}$"); + + prm.enter_subsection ("Erosional parameters"); + { + prm.declare_entry("Drainage area exponent", "0.4", + Patterns::Double(), + "Exponent for drainage area."); + prm.declare_entry("Slope exponent", "1", + Patterns::Double(), + "The slope exponent for SPL (n). Generally m/n should equal approximately 0.4"); + prm.declare_entry("Bedrock river incision rate", "1e-5", + Patterns::Double(), + "River incision rate for bedrock in the Stream Power Law. Units: $\\{m^(1-2*drainage_area_exponent)/yr}$"); + prm.declare_entry("Sediment river incision rate", "-1", + Patterns::Double(), + "River incision rate for sediment in the Stream Power Law. -1 sets this to the bedrock river incision rate. Units: $\\{m^(1-2*drainage_area_exponent)/yr}$ "); + prm.declare_entry("Bedrock diffusivity", "1e-2", + Patterns::Double(), + "Transport coefficient (diffusivity) for bedrock. Units: $\\{m^2/yr}$ "); + prm.declare_entry("Sediment diffusivity", "-1", + Patterns::Double(), + "Transport coefficient (diffusivity) for sediment. -1 sets this to the bedrock diffusivity. Units: $\\{m^2/yr}$"); + prm.declare_entry("Elevation factor", "1", + Patterns::Double(), + "Amount to multiply kf and kd by past given orographic elevation control."); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + prm.leave_subsection (); + } + + + template + void FastScapecc::parse_parameters(ParameterHandler &prm) + { + end_time = prm.get_double ("End time"); + if (prm.get_bool ("Use years in output instead of seconds") == true) + end_time *= year_in_seconds; + + prm.enter_subsection ("Mesh deformation"); + { + prm.enter_subsection("Fastscapecc"); + { + fastscape_steps_per_aspect_step = prm.get_integer("Number of steps"); + maximum_fastscape_timestep = prm.get_double("Maximum timestep"); + additional_refinement_levels = prm.get_integer("Additional fastscape refinement levels"); + fs_seed = prm.get_integer("Fastscape seed"); + maximum_surface_refinement_level = prm.get_integer("Maximum surface refinement level"); + surface_refinement_difference = prm.get_integer("Surface refinement difference"); + precision = prm.get_double("Precision"); + noise_h = prm.get_double("Initial noise magnitude"); + + if (!this->convert_output_to_years()) + { + maximum_fastscape_timestep /= year_in_seconds; + } + + prm.enter_subsection("Erosional parameters"); + { + m = prm.get_double("Drainage area exponent"); + n = prm.get_double("Slope exponent"); + kfsed = prm.get_double("Sediment river incision rate"); + kff = prm.get_double("Bedrock river incision rate"); + kdsed = prm.get_double("Sediment diffusivity"); + kdd = prm.get_double("Bedrock diffusivity"); + + // if (!this->convert_output_to_years()) + // { + // kff *= year_in_seconds; + // kdd *= year_in_seconds; + // kfsed *= year_in_seconds; + // kdsed *= year_in_seconds; + // } + + } + prm.leave_subsection(); + + } + prm.leave_subsection(); + } + prm.leave_subsection (); + } + } +} + + +// explicit instantiation of the functions we implement in this file +namespace aspect +{ + namespace MeshDeformation + { + ASPECT_REGISTER_MESH_DEFORMATION_MODEL(FastScapecc, + "fastscapecc", + "A plugin, which prescribes the surface mesh to " + "deform according to an analytically prescribed " + "function. Note that the function prescribes a " + "deformation velocity, i.e. the return value of " + "this plugin is later multiplied by the time step length " + "to compute the displacement increment in this time step. " + "The format of the " + "functions follows the syntax understood by the " + "muparser library, see Section~\\ref{sec:muparser-format}.") + } +} diff --git a/tests/fastscape_eroding_box_cc.prm b/tests/fastscape_eroding_box_cc.prm new file mode 100644 index 00000000000..c42661ecf99 --- /dev/null +++ b/tests/fastscape_eroding_box_cc.prm @@ -0,0 +1,195 @@ +# At the top, we define the number of space dimensions we would like to +# work in: +set Dimension = 3 + +set Use years in output instead of seconds = true +set End time = 1e6 +set Maximum time step = 1000 +set Output directory = eroding_box + +set Pressure normalization = no +set Surface pressure = 0 + +# A box which has one half raised by 1 km. +subsection Geometry model + set Model name = box + + subsection Box + set X extent = 100e3 + set Y extent = 100e3 + set Z extent = 25e3 + set X repetitions = 5 + set Y repetitions = 5 + end +end + +# subsection Initial topography model +# set Model name = function +# subsection Function +# set Function expression = if(y>20e3 && y<80e3 && x>20e3 && x<80e3, 2500, 0) +# end +# end +# end + +# We do not consider temperature in this setup +subsection Initial temperature model + set Model name = function + + subsection Function + set Function expression = 0 + end +end + +subsection Boundary temperature model + set Fixed temperature boundary indicators = bottom, top + set List of model names = initial temperature +end + +# X Set all boundaries except the top to freeslip. +subsection Boundary velocity model + set Tangential velocity boundary indicators = bottom, front, back, right, left +# set Prescribed velocity boundary indicators = left x:function +# subsection Function +# set Variable names = x,y,t +# set Function constants = cm=0.01,year=1 +# set Function expression = if(x<50e3, 1*cm/year,0) ;0 +# end +end + + + + +# Here we setup the top boundary with fastscape. +subsection Mesh deformation + set Mesh deformation boundary indicators = top : fastscapecc + + subsection Fastscapecc + # As the highest resolution at the surface is 4 in the mesh refinement function, we set this the same. + set Maximum surface refinement level = 0 + + # Because we only have the same level across the entire surface we set this to zero. + # If we had multiple resolutions, e.g. 3 different cell sizes at the surface, + # this would be set to 2 (number of cell sizes at surface - 1). + set Surface refinement difference = 0 + + # We set FastScape as 1 level more resolved than ASPECT's initial global surface resolution. + set Additional fastscape refinement = 0 + + # Number of FastScape timesteps we want per ASPECT timestep. + set Number of steps = 5 + + # Set the maximum allowable FastScape timestep. If the (ASPECT timestep)/(Number of steps) + # is greater than this value, then the number of steps is automatically doubled. + # This is continued until the FastScape timestep is less than this value. + set Maximum timestep = 10000 + + # Seed number for initial topography noise + set Fastscape seed = 0 + + + # Fix all boundaries in FastScape. + subsection Boundary conditions + set Front = 1 + set Right = 1 + set Back = 1 + set Left = 1 + end + + # We use both the Steam Power Law (SPL) and hillslope diffusion. The bedrock deposition + # variable allows deposition of sediment, however, because + # the sediment variables are not set, they are equal to + # bedrock values. + subsection Erosional parameters + set Drainage area exponent = 0.4 + set Bedrock diffusivity = 1e-2 + set Bedrock river incision rate = 1e-4 + set Slope exponent = 1 + end + end + + set Additional tangential mesh velocity boundary indicators = left, right, front, back + +end + +subsection Gravity model + set Model name = vertical + + subsection Vertical + set Magnitude = 0 + end +end + +# Dimensionless +subsection Material model + set Model name = simple + + subsection Simple model + set Reference density = 1 + set Reference specific heat = 1 + set Reference temperature = 0 + set Thermal conductivity = 1 + set Thermal expansion coefficient = 1 + set Viscosity = 1 + end +end + + +# We also have to specify that we want to use the Boussinesq +# approximation (assuming the density in the temperature +# equation to be constant, and incompressibility). +subsection Formulation + set Formulation = Boussinesq approximation +end + + +# We set a maximum surface refinement level of 4 using the max/minimum refinement level. +# This will make it so that the ASPECT surface refinement is the same as what we set +# for our maximum surface refinement level for FastScape. +# subsection Mesh refinement +# set Initial global refinement = 4 +# set Initial adaptive refinement = 1 +# set Time steps between mesh refinement = 0 +# set Strategy = minimum refinement function, maximum refinement function + +# subsection Maximum refinement function +# set Coordinate system = cartesian +# set Variable names = x,y,z +# set Function expression = if(z>=21000, 4, 2) +# end + +# subsection Minimum refinement function +# set Coordinate system = cartesian +# set Variable names = x,y,z +# set Function expression = if(z>=21000, 4, 2) +# end +# end + +subsection Mesh refinement + set Initial global refinement = 0 + set Initial adaptive refinement = 0 + set Time steps between mesh refinement = 0 + # set Strategy = minimum refinement function, maximum refinement function + + # subsection Maximum refinement function + # set Coordinate system = cartesian + # set Variable names = x,y,z + # set Function expression = if(z>=21000, 4, 2) + # end + + # subsection Minimum refinement function + # set Coordinate system = cartesian + # set Variable names = x,y,z + # set Function expression = if(z>=21000, 4, 2) + # end +end + +subsection Postprocess + set List of postprocessors = visualization + + subsection Visualization + set Time between graphical output = 0 + set List of output variables = strain rate + set Interpolate output = false + set Write higher order output = false + end +end diff --git a/tests/fastscape_eroding_box_cc/output_screen b/tests/fastscape_eroding_box_cc/output_screen new file mode 100644 index 00000000000..22d10c18c6e --- /dev/null +++ b/tests/fastscape_eroding_box_cc/output_screen @@ -0,0 +1,277 @@ +----------------------------------------------------------------------------- +-- This is ASPECT, the Advanced Solver for Problems in Earth's ConvecTion. +-- . version 2.6.0-pre (coupling_fs, 19ad382da) +-- . using deal.II 9.4.0 +-- . with 32 bit indices and vectorization level 0 (64 bits) +-- . using Trilinos 13.2.0 +-- . using p4est 2.3.2 +-- . running in DEBUG mode +-- . running with 1 MPI process +----------------------------------------------------------------------------- + + +----------------------------------------------------------------------------- +The output directory provided in the input file appears not to exist. +ASPECT will create it for you. +----------------------------------------------------------------------------- + + +----------------------------------------------------------------------------- +-- For information on how to cite ASPECT, see: +-- https://aspect.geodynamics.org/citing.html?ver=2.6.0-pre&sha=19ad382da&src=code +----------------------------------------------------------------------------- +Number of active cells: 25 (on 1 levels) +Number of degrees of freedom: 1,524 (1,089+72+363) + +Number of mesh deformation degrees of freedom: 216 + Solving mesh displacement system... 0 iterations. +*** Timestep 0: t=0 years, dt=0 years + Solving mesh displacement system... 0 iterations. + Skipping temperature solve because RHS is zero. + Rebuilding Stokes preconditioner... + Solving Stokes system... 0+0 iterations. + + Postprocessing: + Writing graphical output: eroding_box/solution/solution-00000 + +*** Timestep 1: t=1000 years, dt=1000 years + Initializing FastScape... 1 levels, cell size: 20000 m. +here it works 2 +here it works 3 +here it works 4 +here it works 5 +here it works 6 +here it works 7 +here it works 8 +Fastscape ite 0 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +here it works 9 +here it works 9b +here it works 8 +Fastscape ite 1 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +here it works 9 +here it works 9b +here it works 8 +Fastscape ite 2 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +here it works 9 +here it works 9b +here it works 8 +Fastscape ite 3 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +here it works 9 +here it works 9b +here it works 8 +Fastscape ite 4 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +here it works 9 +here it works 9b +here it works 10 +here it works 11 +here it works 12 + Solving mesh displacement system... 0 iterations. + Skipping temperature solve because RHS is zero. + Rebuilding Stokes preconditioner... + Solving Stokes system... 0+0 iterations. + + Postprocessing: + Writing graphical output: eroding_box/solution/solution-00001 + +*** Timestep 2: t=2000 years, dt=1000 years + Initializing FastScape... 1 levels, cell size: 20000 m. +here it works 2 +here it works 3 +here it works 4 +here it works 5 +here it works 6 +here it works 7 +here it works 8 +Fastscape ite 0 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +here it works 9 +here it works 9b +here it works 8 +Fastscape ite 1 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +here it works 9 +here it works 9b +here it works 8 +Fastscape ite 2 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +here it works 9 +here it works 9b +here it works 8 +Fastscape ite 3 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +here it works 9 +here it works 9b +here it works 8 +Fastscape ite 4 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +here it works 9 +here it works 9b +here it works 10 +here it works 11 +here it works 12 + Solving mesh displacement system... 0 iterations. + Skipping temperature solve because RHS is zero. + Rebuilding Stokes preconditioner... + Solving Stokes system... 0+0 iterations. + + Postprocessing: + Writing graphical output: eroding_box/solution/solution-00002 + +*** Timestep 3: t=3000 years, dt=1000 years + Initializing FastScape... 1 levels, cell size: 20000 m. +here it works 2 +here it works 3 +here it works 4 +here it works 5 +here it works 6 +here it works 7 +here it works 8 +Fastscape ite 0 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +here it works 9 +here it works 9b +here it works 8 +Fastscape ite 1 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +here it works 9 +here it works 9b +here it works 8 +Fastscape ite 2 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +here it works 9 +here it works 9b +here it works 8 +Fastscape ite 3 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +here it works 9 +here it works 9b +here it works 8 +Fastscape ite 4 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +here it works 9 +here it works 9b +here it works 10 +here it works 11 +here it works 12 + Solving mesh displacement system... 0 iterations. + Skipping temperature solve because RHS is zero. + Rebuilding Stokes preconditioner... + Solving Stokes system... 0+0 iterations. + + Postprocessing: + Writing graphical output: eroding_box/solution/solution-00003 + +*** Timestep 4: t=4000 years, dt=1000 years + Initializing FastScape... 1 levels, cell size: 20000 m. +here it works 2 +here it works 3 +here it works 4 +here it works 5 +here it works 6 +here it works 7 +here it works 8 +Fastscape ite 0 +here it works 8b +here it works 8c +here it works 8d +here it works 8e +[MacBook-Pro:08191] *** Process received signal *** +[MacBook-Pro:08191] Signal: Abort trap: 6 (6) +[MacBook-Pro:08191] Signal code: (0) +[MacBook-Pro:08191] [ 0] 0 libsystem_platform.dylib 0x000000019ea52a84 _sigtramp + 56 +[MacBook-Pro:08191] [ 1] 0 libsystem_pthread.dylib 0x000000019ea23c28 pthread_kill + +288 +[MacBook-Pro:08191] [ 2] 0 libsystem_c.dylib 0x000000019e931ae8 abort + 180 +[MacBook-Pro:08191] [ 3] 0 libsystem_malloc.dylib 0x000000019e852e28 malloc_vreport ++ 908 +[MacBook-Pro:08191] [ 4] 0 libsystem_malloc.dylib 0x000000019e8695d4 +malloc_zone_error + 104 +[MacBook-Pro:08191] [ 5] 0 libsystem_malloc.dylib 0x000000019e861148 +nanov2_guard_corruption_detected + 44 +[MacBook-Pro:08191] [ 6] 0 libsystem_malloc.dylib 0x000000019e860344 _nanov2_free + +0 +[MacBook-Pro:08191] [ 7] 0 libc++abi.dylib 0x000000019e9de924 _Znwm + 32 +[MacBook-Pro:08191] [ 8] 0 aspect 0x0000000105070738 +_ZN3xtl3mpl9static_ifIZN2xt18assign_xexpressionINS2_17xtensor_containerINS2_7uvectorIdNSt3__19allocatorIdEEEELm2ELNS2_11layout_typeE1ENS2_22xtensor_expression_tagEEENS2_13xstrided_viewIRSC_NS6_5arrayImLm2EEELSA_0ENS2_6detail20inner_storage_getterISE_EEEEEEvRNS2_11xexpressionIT_EERKNSL_IT0_EEEUlSM_E_ZNS3_ISC_SK_EEvSO_SS_EUlSM_E0_EEDcNS6_17integral_constantIbLb0EEERKSM_RKSP_ ++ 188 +[MacBook-Pro:08191] [ 9] 0 aspect 0x000000010506e1b8 +_ZN12fastscapelib20diffusion_adi_eroderINS_14raster_grid_xtINS_11xt_selectorELNS_14raster_connectE1ENS_15neighbors_cacheILh8EEEEES2_E13solve_adi_rowIN2xt13xstrided_viewIRNS9_17xtensor_containerINS9_7uvectorIdNSt3__19allocatorIdEEEELm2ELNS9_11layout_typeE1ENS9_22xtensor_expression_tagEEENSD_5arrayImLm2EEELSH_0ENS9_6detail20inner_storage_getterISK_EEEENSA_IRNSB_ISG_Lm3ELSH_1ESI_EENSL_ImLm3EEELSH_0ENSO_ISS_EEEESV_EEDaOT_OT0_OT1_mmd ++ 144 +[MacBook-Pro:08191] [10] 0 aspect 0x000000010501a780 +_ZN12fastscapelib20diffusion_adi_eroderINS_14raster_grid_xtINS_11xt_selectorELNS_14raster_connectE1ENS_15neighbors_cacheILh8EEEEES2_E5erodeERKN2xt16xarray_containerINS8_7uvectorIdNSt3__19allocatorIdEEEELNS8_11layout_typeE1ENS8_7svectorImLm4ENSC_ImEELb1EEENS8_22xtensor_expression_tagEEEd ++ 400 +[MacBook-Pro:08191] [11] 0 aspect 0x000000010501f064 +_ZNK6aspect15MeshDeformation11FastScapeccILi3EE40compute_velocity_constraints_on_boundaryERKN6dealii10DoFHandlerILi3ELi3EEERNS3_17AffineConstraintsIdEERKNSt3__13setIjNSB_4lessIjEENSB_9allocatorIjEEEE ++ 6844 +[MacBook-Pro:08191] [12] 0 aspect 0x0000000105039b38 +_ZN6aspect15MeshDeformation22MeshDeformationHandlerILi3EE16make_constraintsEv + 1036 +[MacBook-Pro:08191] [13] 0 aspect 0x0000000105039618 +_ZN6aspect15MeshDeformation22MeshDeformationHandlerILi3EE7executeEv + 152 +[MacBook-Pro:08191] [14] 0 aspect 0x00000001048a16f8 +_ZN6aspect9SimulatorILi3EE14solve_timestepEv + 44 +[MacBook-Pro:08191] [15] 0 aspect 0x00000001048a0084 +_ZN6aspect9SimulatorILi3EE3runEv + 892 +[MacBook-Pro:08191] [16] 0 aspect 0x00000001051b4f44 +_Z13run_simulatorILi3EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb + +868 +[MacBook-Pro:08191] [17] 0 aspect 0x00000001051b41a8 main + 2580 +[MacBook-Pro:08191] [18] 0 dyld 0x000000019e6cbf28 start + 2236 +[MacBook-Pro:08191] *** End of error message *** +-------------------------------------------------------------------------- +Primary job terminated normally, but 1 process returned +a non-zero exit code. Per user-direction, the job has been aborted. +-------------------------------------------------------------------------- +-------------------------------------------------------------------------- +mpirun noticed that process rank 0 with PID 0 on node MacBook-Pro exited on signal 6 (Abort trap: +6). +--------------------------------------------------------------------------