diff --git a/include/aspect/mesh_deformation/landlab.h b/include/aspect/mesh_deformation/landlab.h
new file mode 100644
index 00000000000..51cdeaec92f
--- /dev/null
+++ b/include/aspect/mesh_deformation/landlab.h
@@ -0,0 +1,144 @@
+/*
+ Copyright (C) 2020 - 2024 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_landlab_h
+#define _aspect_mesh_deformation_landlab_h
+
+#include
+#include
+#include
+
+#include
+
+namespace aspect
+{
+ namespace MeshDeformation
+ {
+ /**
+ * A plugin that computes the deformation of surface
+ * vertices by coupling with LandLab.
+ *
+ * @ingroup MeshDeformation
+ */
+ template
+ class LandLab : public Interface, public SimulatorAccess
+ {
+ public:
+ LandLab();
+
+ /**
+ * Initialize function, which sets the start time and
+ * start timestep of diffusion.
+ */
+ void initialize() override;
+
+ /**
+ * The update function sets the current time and determines
+ * whether diffusion should be applied in this timestep.
+ */
+ void update() 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 diffusion of the surface.
+ */
+ static
+ void declare_parameters (ParameterHandler &prm);
+
+ /**
+ * Parse parameters for the diffusion of the surface.
+ */
+ void parse_parameters (ParameterHandler &prm) override;
+
+ private:
+ /**
+ * Compute the surface velocity from a difference
+ * in surface height given by the solution of
+ * the hillslope diffusion problem.
+ */
+ void diffuse_boundary (const DoFHandler &free_surface_dof_handler,
+ const IndexSet &mesh_locally_owned,
+ const IndexSet &mesh_locally_relevant,
+ LinearAlgebra::Vector &output,
+ const std::set &boundary_id) const;
+
+ /**
+ * Check that the size of the next time step is not larger than the conduction
+ * timestep based on the mesh size and the diffusivity on each cell along the
+ * surface. The computed time step has to satisfy the CFL
+ * number chosen in the input parameter file on each cell of the mesh.
+ */
+ void check_diffusion_time_step (const DoFHandler &mesh_deformation_dof_handler,
+ const std::set &boundary_ids) const;
+
+
+ /**
+ */
+ std::vector> surface_point_locations;
+
+ /**
+ * The hillslope transport coefficient or diffusivity [m2/s]
+ * used in the hillslope diffusion of the deformed
+ * surface.
+ */
+ double diffusivity;
+
+ /**
+ * Maximum number of steps between the application of diffusion.
+ */
+ unsigned int timesteps_between_diffusion;
+
+ /**
+ * Whether or not in the current timestep, diffusion
+ * should be applied.
+ */
+ bool apply_diffusion;
+
+ /**
+ * Boundaries along which the mesh is allowed to move tangentially
+ * despite of the Stokes velocity boundary conditions.
+ */
+ std::set additional_tangential_mesh_boundary_indicators;
+
+ PyObject *pModule;
+ };
+ }
+}
+
+
+#endif
diff --git a/landlab/landlab.py b/landlab/landlab.py
new file mode 100644
index 00000000000..cfdc9b50efe
--- /dev/null
+++ b/landlab/landlab.py
@@ -0,0 +1,121 @@
+from mpi4py import MPI
+import numpy as np
+
+print("loading landlab.py...")
+
+current_time = 0
+
+comm = None
+
+def initialize(comm_handle):
+ # Convert the handle back to an MPI communicator
+ global comm
+ comm = MPI.Comm.f2py(comm_handle)
+
+ rank = comm.Get_rank()
+ size = comm.Get_size()
+
+ print(f"Python: Hello from Rank {rank} of {size}")
+
+ data = 1
+ globalsum = comm.allreduce(data, op=MPI.SUM)
+ if comm.rank == 0:
+ print(f"\tPython: testing communication; sum {globalsum}")
+
+def finalize():
+ pass
+
+
+# Run the Landlab simulation from the current time to end_time and return
+# the new topographic elevation (in m) at each local node.
+# dict_variable_name_to_value_in_nodes is a dictionary mapping variables
+# (x velocity, y velocity, temperature, etc.) to an array of values in each
+# node.
+def update_until(end_time, dict_variable_name_to_value_in_nodes):
+ dt = end_time-current_time
+ x_velocity = dict_variable_name_to_value_in_nodes["x velocity"]
+ y_velocity = dict_variable_name_to_value_in_nodes["y velocity"]
+ z_velocity = dict_variable_name_to_value_in_nodes["z velocity"]
+
+ if dt>0:
+ n_substeps = 10
+ sub_dt = dt / n_substeps
+ for _ in range(n_substeps):
+
+ uplift(z_velocity * sub_dt)
+ advect(x_velocity * sub_dt, y_velocity * sub_dt)
+ elevation_before = "topographic__elevation"
+ run all components
+ deposition_erosion += "topographic__elevation" - elevation_before
+ pass
+
+ current_time = end_time
+
+ return deposition_erosion
+
+
+mg = None
+
+# Notify the code to create the parallel grid. Done
+# once at the beginning of the simulation and later,
+# if the ASPECT volume mesh changes. This can be
+# ignored if no adaptivity is supported on the
+# Landlab side.
+def set_mesh_information(dict_grid_information):
+ if not mg:
+ if rank==0:
+ global_grid = HexGrid()
+ partition_grid()
+ broadcast_pieces()
+ else
+ mg = receive_my_piece()
+
+# Return the x coordinates of the locally owned nodes on this
+# MPI rank. grid_id is always 0.
+def get_grid_x(grid_id):
+ return mg.node_x
+
+# Return the y coordinates of the locally owned nodes on this
+# MPI rank. grid_id is always 0.
+def get_grid_y(grid_id):
+ return mg.node_y
+
+# Return the initial topography at the start of the simulation
+# in each node.
+def get_initial_topography(grid_id):
+ return topographic__elevation
+
+
+
+# old:
+
+def get_grid_nodes(dict_grid_information):
+ numpoints = 5
+ xcoords = np.array([np.random.uniform(0,1) for i in range(numpoints)])
+ ycoords = np.array([np.random.uniform(0,1) for i in range(numpoints)])
+ return xcoords, ycoords
+
+def do_timestep():
+ print("hello")
+
+def start_update(t, vec):
+ print(f"update at time {t}, {vec}")
+
+def write_output():
+ pass
+
+def update_single(x,y,old):
+ print(f"{x} {y}")
+ n = old + np.random.uniform(-1,1)*0.1
+ return n
+
+if __name__ == "__main__":
+ comm = MPI.COMM_WORLD
+ initialize(MPI.comm.py2f(comm))
+
+ set_mesh_information({})
+
+ dt = 0.1
+ for n in range(10):
+ update_until(n*dt)
+ write_output()
diff --git a/landlab/test1.prm b/landlab/test1.prm
new file mode 100644
index 00000000000..702d3ca32af
--- /dev/null
+++ b/landlab/test1.prm
@@ -0,0 +1,115 @@
+set Dimension = 3
+set Use years in output instead of seconds = false
+set End time = 1e6
+set Maximum time step = 1e2
+set Nonlinear solver scheme = no Advection, no Stokes
+set Pressure normalization = surface
+set Surface pressure = 0
+
+# 1x1 box with an initial hill topography
+subsection Geometry model
+ set Model name = box
+
+ subsection Box
+ set X extent = 1
+ set Y extent = 1
+ set Z extent = 1
+ end
+
+ subsection Initial topography model
+ set Model name = function
+
+ subsection Function
+ set Function constants = A=0.3
+ set Function expression = \
+ A * sin(x*pi) * sin(y*pi)
+ end
+ end
+end
+
+# Temperature effects are ignored
+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
+
+# Free slip on all boundaries
+subsection Boundary velocity model
+ set Tangential velocity boundary indicators = left, right, bottom, top, front, back
+end
+
+# The mesh will be deformed according to the displacement
+# of the surface nodes due to diffusion of the topography.
+# The mesh is allowed to move vertical along the left and
+# right boundary.
+subsection Mesh deformation
+ set Mesh deformation boundary indicators = top: landlab
+ set Additional tangential mesh velocity boundary indicators = left, right, front, back
+
+ subsection Diffusion
+ # The diffusivity
+ set Hillslope transport coefficient = 0.25
+ end
+end
+
+# Vertical gravity
+subsection Gravity model
+ set Model name = vertical
+
+ subsection Vertical
+ set Magnitude = 1
+ end
+end
+
+# One material with unity properties
+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 here use a globally refined mesh without
+# adaptive mesh refinement.
+subsection Mesh refinement
+ set Initial global refinement = 3
+ set Initial adaptive refinement = 0
+ set Time steps between mesh refinement = 0
+end
+
+subsection Postprocess
+ set List of postprocessors = visualization
+
+ subsection Visualization
+ set Time between graphical output = 0
+ set Output mesh velocity = true
+ set Interpolate output = false
+ end
+end
+
+subsection Solver parameters
+ subsection Stokes solver parameters
+ set Linear solver tolerance = 1e-7
+ end
+end
diff --git a/source/mesh_deformation/landlab.cc b/source/mesh_deformation/landlab.cc
new file mode 100644
index 00000000000..2b63452cbf0
--- /dev/null
+++ b/source/mesh_deformation/landlab.cc
@@ -0,0 +1,725 @@
+/*
+ Copyright (C) 2020 - 2024 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
+
+
+#define PY_SSIZE_T_CLEAN
+#include
+#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
+#include
+
+PyObject *vector_to_numpy(const std::vector &vec)
+{
+ npy_intp size = vec.size();
+ double *data = const_cast(vec.data());
+ return PyArray_SimpleNewFromData(1, &size, NPY_DOUBLE, data);
+}
+
+
+
+namespace aspect
+{
+
+ namespace MeshDeformation
+ {
+ template
+ LandLab::LandLab()
+ :
+ diffusivity(0.),
+ timesteps_between_diffusion (1),
+ apply_diffusion(false)
+ {}
+
+
+
+ template
+ void
+ LandLab::initialize ()
+ {
+ AssertThrow(Plugins::plugin_type_matches>(this->get_geometry_model()) ||
+ Plugins::plugin_type_matches>(this->get_geometry_model()),
+ ExcMessage("The surface diffusion mesh deformation plugin only works for Box geometries."));
+
+ Py_Initialize();
+ _import_array(); // required for Numpy interop
+
+ PyRun_SimpleString("import sys; sys.path.append(\"" ASPECT_SOURCE_DIR "/landlab\")");
+
+ // Import Python module
+ pModule = PyImport_ImportModule("landlab");
+ AssertThrow(pModule, ExcMessage("Failed to load Python module"));
+
+ // Get function from module
+ PyObject *pFunc = PyObject_GetAttrString(pModule, "initialize");
+ if (!pFunc || !PyCallable_Check(pFunc))
+ {
+ std::cerr << "Failed to load function" << std::endl;
+ Py_DECREF(pModule);
+ Py_Finalize();
+ return;
+ }
+
+ // Call Python function with communicator handle
+ PyObject *pArgs = PyTuple_Pack(1, PyLong_FromLong(MPI_Comm_c2f(this->get_mpi_communicator())));
+ //PyObject* pValue =
+ PyObject_CallObject(pFunc, pArgs);
+
+ // Clean up
+ Py_DECREF(pArgs);
+ Py_DECREF(pFunc);
+// Py_DECREF(pModule);
+ }
+
+
+
+ template
+ void
+ LandLab::update ()
+ {
+ // Set the current timestep number
+ const unsigned int current_timestep_number = this->get_timestep_number();
+
+ // Determine whether we need to apply diffusion based
+ // on the timestep interval between applications.
+ // Because mesh deformation is computed at the beginning
+ // of each timestep, we do not apply diffusion in the first timestep.
+ if (current_timestep_number != 0)
+ {
+ if (current_timestep_number % timesteps_between_diffusion == 0)
+ apply_diffusion = true;
+ else
+ apply_diffusion = false;
+ }
+
+
+ {
+ PyObject *pFunc = PyObject_GetAttrString(pModule, "define_mesh");
+ PyObject *result = PyObject_CallObject(pFunc, nullptr);
+ PyArrayObject *xcoord_obj = reinterpret_cast(PyTuple_GetItem(result, 0));
+ PyArrayObject *ycoord_obj = reinterpret_cast(PyTuple_GetItem(result, 1));
+
+ double *xcoord = static_cast(PyArray_DATA(xcoord_obj));
+ double *ycoord = static_cast(PyArray_DATA(ycoord_obj));
+
+ npy_intp size1 = PyArray_SIZE(xcoord_obj);
+ npy_intp size2 = PyArray_SIZE(ycoord_obj);
+ AssertThrow(size1 == size2, ExcMessage("return values do not match."));
+
+ surface_point_locations.resize(size1);
+ for (unsigned i=0; i nat_position = {{xcoord[i], ycoord[i]}};
+ nat_position[dim-1] = 1.0; // depth=0 is z=1
+ surface_point_locations[i] = this->get_geometry_model().natural_to_cartesian_coordinates(nat_position);
+
+ std::cout << xcoord[i] << " - " << ycoord[i] <
+ void LandLab::diffuse_boundary(const DoFHandler &mesh_deformation_dof_handler,
+ const IndexSet &mesh_locally_owned,
+ const IndexSet &mesh_locally_relevant,
+ LinearAlgebra::Vector &output,
+ const std::set &boundary_ids) const
+ {
+
+
+ {
+ LinearAlgebra::Vector displacements = this->get_mesh_deformation_handler().get_mesh_displacements();
+
+ Utilities::MPI::RemotePointEvaluation rpe;
+ MappingQ mapping(1);
+ rpe.reinit(surface_point_locations, this->get_triangulation(), mapping);
+
+ const unsigned int myrank = Utilities::MPI::this_mpi_process(this->get_mpi_communicator());
+
+ std::vector values1 = VectorTools::point_values<1>(rpe,
+ mesh_deformation_dof_handler,
+ displacements,
+ dealii::VectorTools::EvaluationFlags::avg, /*component*/ 2);
+ std::vector values2 = VectorTools::point_values<1>(rpe,
+ mesh_deformation_dof_handler,
+ this->get_mesh_deformation_handler().get_initial_topography(),
+ dealii::VectorTools::EvaluationFlags::avg, /*component*/ 2);
+
+ for (unsigned int i=0; i &values,
+ const typename Utilities::MPI::RemotePointEvaluation::CellData &cell_data)
+ {
+// TODO
+ };
+
+// rpe.template process_and_evaluate(values, eval_func, /*sort_data*/ true);
+
+
+ }
+
+
+ {
+ //def start_update(t):
+ PyObject *pFunc = PyObject_GetAttrString(pModule, "start_update");
+ if (!pFunc || !PyCallable_Check(pFunc))
+ {
+ std::cerr << "Failed to load function" << std::endl;
+ Py_DECREF(pModule);
+ Py_Finalize();
+ return;
+ }
+ std::vector vec = {1.0, 2.0, 3.0};
+ PyObject *arr = vector_to_numpy(vec);
+
+ // Call Python function with communicator handle
+ PyObject *pArgs = PyTuple_Pack(2, PyFloat_FromDouble(this->get_time()), arr);
+ //PyObject* pValue =
+ PyObject_CallObject(pFunc, pArgs);
+
+
+
+
+ Py_DECREF(arr);
+
+ // Clean up
+ Py_DECREF(pArgs);
+ Py_DECREF(pFunc);
+
+ }
+
+
+ // Check that the current timestep does not exceed the diffusion timestep
+ check_diffusion_time_step(mesh_deformation_dof_handler, boundary_ids);
+
+ // Set up constraints
+#if DEAL_II_VERSION_GTE(9,6,0)
+ AffineConstraints matrix_constraints(mesh_locally_relevant, mesh_locally_relevant);
+#else
+ AffineConstraints matrix_constraints(mesh_locally_relevant);
+#endif
+
+ DoFTools::make_hanging_node_constraints(mesh_deformation_dof_handler, matrix_constraints);
+
+ std::set periodic_boundary_indicators;
+ using periodic_boundary_pairs = std::set, unsigned int>>;
+ const periodic_boundary_pairs pbp = this->get_geometry_model().get_periodic_boundary_pairs();
+ for (const auto &p : pbp)
+ {
+ DoFTools::make_periodicity_constraints(mesh_deformation_dof_handler,
+ p.first.first, p.first.second, p.second, matrix_constraints);
+ periodic_boundary_indicators.insert(p.first.first);
+ periodic_boundary_indicators.insert(p.first.second);
+ }
+
+ // A boundary can be:
+ // 1) A mesh deformation boundary with its deformation computed by this plugin.
+ // 2) A zero mesh displacement boundary:
+ // a) boundaries with zero Stokes velocity that are not an Additional tangential mesh boundary or mesh deformation boundary,
+ // b) boundaries with prescribed Stokes velocity that are not an Additional tangential mesh boundary or mesh deformation boundary.
+ // 3) A tangential mesh displacement boundary (no normal flux):
+ // a) tangential Stokes boundaries that are not a mesh deformation boundary,
+ // b) additional tangential mesh boundaries other than those in a) and c),
+ // c) periodic boundaries that are not a mesh deformation boundary.
+ // We obtain the zero mesh displacement boundaries by subtracting the mesh deformation ones (1+3) from all used boundary indicators. In this
+ // case, when no Stokes velocity boundary conditions are set (e.g. when using the single Advection,
+ // no Stokes solver scheme, we do end up with mesh displacement boundary conditions.
+
+ // All boundary indicators used by the current geometry.
+ const std::set all_boundaries = this->get_geometry_model().get_used_boundary_indicators();
+
+ // Get the tangential Stokes velocity boundary indicators.
+ const std::set tangential_velocity_boundary_indicators =
+ this->get_boundary_velocity_manager().get_tangential_boundary_velocity_indicators();
+
+ // The list of all tangential mesh deformation boundary indicators.
+ std::set tmp_tangential_boundaries, tmp2_tangential_boundaries, all_tangential_boundaries;
+ std::set_union(tangential_velocity_boundary_indicators.begin(),tangential_velocity_boundary_indicators.end(),
+ additional_tangential_mesh_boundary_indicators.begin(),additional_tangential_mesh_boundary_indicators.end(),
+ std::inserter(tmp_tangential_boundaries, tmp_tangential_boundaries.end()));
+ std::set_union(tmp_tangential_boundaries.begin(),tmp_tangential_boundaries.end(),
+ periodic_boundary_indicators.begin(),periodic_boundary_indicators.end(),
+ std::inserter(tmp2_tangential_boundaries, tmp2_tangential_boundaries.end()));
+ // Tangential Stokes velocity boundaries can also be mesh deformation boundaries,
+ // so correct for that.
+ std::set_difference(tmp2_tangential_boundaries.begin(),tmp2_tangential_boundaries.end(),
+ boundary_ids.begin(),boundary_ids.end(),
+ std::inserter(all_tangential_boundaries, all_tangential_boundaries.end()));
+
+ // The list of mesh deformation boundary indicators of boundaries that are allowed to move either
+ // tangentially or in all directions.
+ std::set all_deformable_boundaries;
+ std::set_union(all_tangential_boundaries.begin(),all_tangential_boundaries.end(),
+ boundary_ids.begin(),boundary_ids.end(),
+ std::inserter(all_deformable_boundaries, all_deformable_boundaries.end()));
+
+ // The list of mesh deformation boundary indicators of boundaries that are not allowed to move.
+ std::set all_fixed_boundaries;
+ std::set_difference(all_boundaries.begin(),all_boundaries.end(),
+ all_deformable_boundaries.begin(),all_deformable_boundaries.end(),
+ std::inserter(all_fixed_boundaries, all_fixed_boundaries.end()));
+
+ // Make the no flux boundary constraints
+ for (const types::boundary_id &boundary_id : all_fixed_boundaries)
+ {
+ VectorTools::interpolate_boundary_values (this->get_mapping(),
+ mesh_deformation_dof_handler,
+ boundary_id,
+ dealii::Functions::ZeroFunction(dim),
+ matrix_constraints);
+ }
+
+ // Make the no normal flux boundary constraints
+ VectorTools::compute_no_normal_flux_constraints (mesh_deformation_dof_handler,
+ /* first_vector_component= */
+ 0,
+ all_tangential_boundaries,
+ matrix_constraints, this->get_mapping());
+
+ matrix_constraints.close();
+
+ // Set up the system to solve
+ LinearAlgebra::SparseMatrix matrix;
+
+ // Sparsity of the matrix
+ TrilinosWrappers::SparsityPattern sp (mesh_locally_owned,
+ mesh_locally_owned,
+ mesh_locally_relevant,
+ this->get_mpi_communicator());
+ DoFTools::make_sparsity_pattern (mesh_deformation_dof_handler, sp, matrix_constraints, false,
+ Utilities::MPI::this_mpi_process(this->get_mpi_communicator()));
+
+ sp.compress();
+ matrix.reinit (sp);
+
+ LinearAlgebra::Vector system_rhs, solution;
+ system_rhs.reinit(mesh_locally_owned, this->get_mpi_communicator());
+ solution.reinit(mesh_locally_owned, this->get_mpi_communicator());
+
+ // Initialize Gauss-Legendre quadrature for degree+1 quadrature points of the surface faces
+ const QGauss face_quadrature(mesh_deformation_dof_handler.get_fe().degree+1);
+ // Update shape function values and gradients, the quadrature points and the Jacobian x quadrature weights.
+ const UpdateFlags update_flags = UpdateFlags(update_values | update_gradients | update_quadrature_points | update_normal_vectors | update_JxW_values);
+ // We want to extract the displacement at the free surface faces of the mesh deformation element.
+ FEFaceValues fs_fe_face_values (this->get_mapping(), mesh_deformation_dof_handler.get_fe(), face_quadrature, update_flags);
+ // and to solve on the whole mesh deformation mesh
+ // The number of quadrature points on a mesh deformation surface face
+ const unsigned int n_fs_face_q_points = fs_fe_face_values.n_quadrature_points;
+
+ // What we need to build our system on the mesh deformation element
+
+ // The nr of shape functions per mesh deformation element
+ const unsigned int dofs_per_cell = mesh_deformation_dof_handler.get_fe().dofs_per_cell;
+
+ // Map of local to global cell dof indices
+ std::vector cell_dof_indices (dofs_per_cell);
+
+ // The local rhs vector
+ Vector cell_vector (dofs_per_cell);
+ // The local matrix
+ FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell);
+
+ // Vector for getting the local dim displacement values
+ std::vector> displacement_values(n_fs_face_q_points);
+
+ // Vector for getting the local dim initial topography values
+ std::vector> initial_topography_values(n_fs_face_q_points);
+
+ // The global displacements on the MeshDeformation FE
+ LinearAlgebra::Vector displacements = this->get_mesh_deformation_handler().get_mesh_displacements();
+
+ // The global initial topography on the MeshDeformation FE
+ // TODO Once the initial mesh deformation is ready, this
+ // can be removed.
+ LinearAlgebra::Vector initial_topography = this->get_mesh_deformation_handler().get_initial_topography();
+
+ // Do nothing at time zero
+ if (this->get_timestep_number() < 1)
+ return;
+
+ // An extractor for the dim-valued displacement vectors
+ // Later on we will compute the gravity-parallel displacement
+ FEValuesExtractors::Vector extract_vertical_displacements(0);
+
+ // An extractor for the dim-valued initial topography vectors
+ // Later on we will compute the gravity-parallel displacement
+ FEValuesExtractors::Vector extract_initial_topography(0);
+
+ // Cell iterator over the MeshDeformation FE
+ typename DoFHandler::active_cell_iterator
+ fscell = mesh_deformation_dof_handler.begin_active(),
+ fsendc = mesh_deformation_dof_handler.end();
+
+ // Iterate over all cells to find those at the mesh deformation boundary
+ for (; fscell!=fsendc; ++fscell)
+ if (fscell->at_boundary() && fscell->is_locally_owned())
+ for (const unsigned int face_no : fscell->face_indices())
+ if (fscell->face(face_no)->at_boundary())
+ {
+ // Boundary indicator of current cell face
+ const types::boundary_id boundary_indicator
+ = fscell->face(face_no)->boundary_id();
+
+ // Only apply diffusion to the requested boundaries
+ if (boundary_ids.find(boundary_indicator) == boundary_ids.end())
+ continue;
+
+ // Recompute values, gradients, etc on the faces
+ fs_fe_face_values.reinit (fscell, face_no);
+
+ // Get the global numbers of the local DoFs of the mesh deformation cell
+ fscell->get_dof_indices (cell_dof_indices);
+
+ // Extract the displacement values
+ fs_fe_face_values[extract_vertical_displacements].get_function_values (displacements, displacement_values);
+
+ // Extract the initial topography values
+ fs_fe_face_values[extract_initial_topography].get_function_values (initial_topography, initial_topography_values);
+
+ // Reset local rhs and matrix
+ cell_vector = 0;
+ cell_matrix = 0;
+
+ // Loop over the quadrature points of the current face
+ for (unsigned int q=0; q direction = -(this->get_gravity_model().gravity_vector(fs_fe_face_values.quadrature_point(q)));
+ // Normalize direction vector
+ if (direction.norm() > 0.0)
+ direction *= 1./direction.norm();
+ // TODO this is only correct for box geometries
+ else
+ {
+ AssertThrow(Plugins::plugin_type_matches>(this->get_geometry_model()) ||
+ Plugins::plugin_type_matches>(this->get_geometry_model()),
+ ExcMessage("The surface diffusion mesh deformation plugin only works for Box geometries."));
+ direction[dim-1] = 1.;
+ }
+
+ // Compute the total displacement in the gravity direction,
+ // i.e. the initial topography + any additional mesh displacement.
+ const double displacement = direction * (displacement_values[q] + initial_topography_values[q]);
+
+ // To project onto the tangent space of the surface,
+ // we define the projection P:= I- n x n,
+ // with I the unit tensor and n the unit normal to the surface.
+ // The surface gradient then is P times the usual gradient of the shape functions.
+ const Tensor<2, dim, double> projection = unit_symmetric_tensor() -
+ outer_product(fs_fe_face_values.normal_vector(q), fs_fe_face_values.normal_vector(q));
+
+ const double JxW = fs_fe_face_values.JxW(q);
+
+ // The shape values for the i-loop
+ std::vector phi(dofs_per_cell);
+
+ // The projected gradients of the shape values for the i-loop
+ std::vector> projected_grad_phi(dofs_per_cell);
+
+ // Loop over the shape functions
+ for (unsigned int i=0; i>(this->get_geometry_model()) ||
+ Plugins::plugin_type_matches>(this->get_geometry_model()),
+ ExcMessage("The surface diffusion mesh deformation plugin only works for Box geometries."));
+ if (mesh_deformation_dof_handler.get_fe().system_to_component_index(i).first == dim-1)
+ {
+ // Precompute shape values and projected shape value gradients
+ phi[i] = fs_fe_face_values.shape_value (i, q);
+ projected_grad_phi[i] = projection * fs_fe_face_values.shape_grad(i, q);
+
+ // Assemble the RHS
+ // RHS = M*H_old
+ cell_vector(i) += phi[i] * displacement * JxW;
+
+ for (unsigned int j=0; jget_timestep() * diffusivity *
+ projected_grad_phi[i] *
+ (projection * fs_fe_face_values.shape_grad(j, q))
+ )
+ * JxW;
+ }
+ }
+ }
+ }
+ }
+
+ matrix_constraints.distribute_local_to_global (cell_matrix, cell_vector,
+ cell_dof_indices, matrix, system_rhs, false);
+ }
+
+ system_rhs.compress (VectorOperation::add);
+ matrix.compress(VectorOperation::add);
+
+ // Jacobi seems to be fine here. Other preconditioners (ILU, IC) run into trouble
+ // because the matrix is mostly empty, since we don't touch internal vertices.
+ LinearAlgebra::PreconditionJacobi preconditioner_mass;
+ LinearAlgebra::PreconditionJacobi::AdditionalData preconditioner_control(1,1e-16,1);
+ preconditioner_mass.initialize(matrix, preconditioner_control);
+
+ this->get_pcout() << " Solving mesh surface diffusion" << std::endl;
+ SolverControl solver_control(5*system_rhs.size(), this->get_parameters().linear_stokes_solver_tolerance*system_rhs.l2_norm());
+ SolverCG cg(solver_control);
+ cg.solve (matrix, solution, system_rhs, preconditioner_mass);
+
+ // Distribute constraints on mass matrix
+ matrix_constraints.distribute (solution);
+
+ // The solution contains the new displacements, but we need to return a velocity.
+ // Therefore, we compute v=d_displacement/d_t.
+ // d_displacement are the new mesh node locations
+ // minus the old locations, which are initial_topography + displacements.
+ LinearAlgebra::Vector velocity(mesh_locally_owned, mesh_locally_relevant, this->get_mpi_communicator());
+ velocity = solution;
+ velocity -= initial_topography;
+ velocity -= displacements;
+
+ // The velocity
+ if (this->get_timestep() > 0.)
+ velocity /= this->get_timestep();
+ else
+ AssertThrow(false, ExcZero());
+
+ output = velocity;
+ }
+
+
+
+ template
+ void LandLab::check_diffusion_time_step (const DoFHandler &mesh_deformation_dof_handler,
+ const std::set &boundary_ids) const
+ {
+ double min_local_conduction_timestep = std::numeric_limits::max();
+
+ for (const auto &fscell : mesh_deformation_dof_handler.active_cell_iterators())
+ if (fscell->at_boundary() && fscell->is_locally_owned())
+ for (const unsigned int face_no : fscell->face_indices())
+ if (fscell->face(face_no)->at_boundary())
+ {
+ // Get the boundary indicator of current cell face
+ const types::boundary_id boundary_indicator
+ = fscell->face(face_no)->boundary_id();
+
+ // Only consider the requested boundaries
+ if (boundary_ids.find(boundary_indicator) == boundary_ids.end())
+ continue;
+
+ // Calculate the corresponding conduction timestep
+ min_local_conduction_timestep = std::min(min_local_conduction_timestep,
+ this->get_parameters().CFL_number*Utilities::fixed_power<2>(fscell->face(face_no)->minimum_vertex_distance())
+ / diffusivity);
+ }
+
+ // Get the global minimum timestep
+ const double min_conduction_timestep = Utilities::MPI::min (min_local_conduction_timestep, this->get_mpi_communicator());
+
+ double conduction_timestep = min_conduction_timestep;
+ if (this->convert_output_to_years())
+ conduction_timestep /= year_in_seconds;
+
+ AssertThrow (this->get_timestep() <= min_conduction_timestep,
+ ExcMessage("The numerical timestep is too large for diffusion of the surface. Although the "
+ "diffusion scheme is stable, note that the error increases linearly with the timestep. "
+ "The diffusion timestep is: " + std::to_string(conduction_timestep) + ", while "
+ "the advection timestep is: " + std::to_string (this->get_timestep ())));
+ }
+
+
+
+ template
+ void
+ LandLab::compute_velocity_constraints_on_boundary(const DoFHandler &mesh_deformation_dof_handler,
+ AffineConstraints &mesh_velocity_constraints,
+ const std::set &boundary_id) const
+ {
+ if (!apply_diffusion)
+ return;
+
+ LinearAlgebra::Vector boundary_velocity;
+
+ const IndexSet &mesh_locally_owned = mesh_deformation_dof_handler.locally_owned_dofs();
+ const IndexSet mesh_locally_relevant = DoFTools::extract_locally_relevant_dofs (mesh_deformation_dof_handler);
+ boundary_velocity.reinit(mesh_locally_owned, mesh_locally_relevant,
+ this->get_mpi_communicator());
+
+ // Determine the mesh velocity at the surface based on diffusion of
+ // the topography
+ diffuse_boundary(mesh_deformation_dof_handler, mesh_locally_owned,
+ mesh_locally_relevant, boundary_velocity, boundary_id);
+
+ // now insert the relevant part of the solution into the mesh constraints
+ const IndexSet constrained_dofs =
+ DoFTools::extract_boundary_dofs(mesh_deformation_dof_handler,
+ ComponentMask(dim, true),
+ boundary_id);
+
+ for (unsigned int i = 0; i < constrained_dofs.n_elements(); ++i)
+ {
+ types::global_dof_index index = constrained_dofs.nth_index_in_set(i);
+ if (mesh_velocity_constraints.can_store_line(index))
+ if (mesh_velocity_constraints.is_constrained(index)==false)
+ {
+ mesh_velocity_constraints.add_line(index);
+ mesh_velocity_constraints.set_inhomogeneity(index, boundary_velocity[index]);
+ }
+ }
+ }
+
+
+
+ template
+ bool
+ LandLab::
+ needs_surface_stabilization () const
+ {
+ return false;
+ }
+
+
+
+ template
+ void LandLab::declare_parameters(ParameterHandler &prm)
+ {
+ prm.enter_subsection("Mesh deformation");
+ {
+ prm.enter_subsection("LandLab");
+ {
+ prm.declare_entry("Hillslope transport coefficient", "1e-6",
+ Patterns::Double(0),
+ "The hillslope transport coefficient $\\kappa$ used to "
+ "diffuse the free surface, either as a "
+ "stabilization step or to mimic erosional "
+ "and depositional processes. Units: $\\si{m^2/s}$. ");
+ prm.declare_entry("Time steps between diffusion", "1",
+ Patterns::Integer(0,std::numeric_limits::max()),
+ "The number of time steps between each application of "
+ "diffusion.");
+ }
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+ }
+
+
+
+ template
+ void LandLab::parse_parameters(ParameterHandler &prm)
+ {
+ prm.enter_subsection ("Mesh deformation");
+ {
+ // Create the list of tangential mesh movement boundary indicators.
+ try
+ {
+ const std::vector x_additional_tangential_mesh_boundary_indicators
+ = this->get_geometry_model().translate_symbolic_boundary_names_to_ids(Utilities::split_string_list
+ (prm.get ("Additional tangential mesh velocity boundary indicators")));
+
+ additional_tangential_mesh_boundary_indicators.insert(x_additional_tangential_mesh_boundary_indicators.begin(),
+ x_additional_tangential_mesh_boundary_indicators.end());
+ }
+ catch (const std::string &error)
+ {
+ AssertThrow (false, ExcMessage ("While parsing the entry , there was an error. Specifically, "
+ "the conversion function complained as follows:\n\n"
+ + error));
+ }
+
+ prm.enter_subsection ("LandLab");
+ {
+ diffusivity = prm.get_double("Hillslope transport coefficient");
+ timesteps_between_diffusion = prm.get_integer("Time steps between diffusion");
+ }
+ 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(LandLab,
+ "landlab",
+ "A plugin that computes the deformation of surface "
+ "vertices according to the solution of the hillslope diffusion problem. "
+ "Specifically, at the end of each timestep, or after a specific number "
+ "of timesteps, this plugin solves the following equation: "
+ "\\begin{align*}"
+ " \\frac{\\partial h}{\\partial t} = \\kappa \\left( \\frac{\\partial^{2} h}{\\partial x^{2}} "
+ "+ \\frac{\\partial^{2} h}{\\partial y^{2}} \\right), "
+ "\\end{align*} "
+ "where $\\kappa$ is the hillslope diffusion coefficient (diffusivity), and $h(x,y)$ the "
+ "height of a point along the top boundary with respect to the surface of the unperturbed domain. "
+ "\n\n"
+ "Using this definition, the plugin then solves for one time step, i.e., "
+ "using as initial condition $h(t_{n-1})$ the current surface elevation, "
+ "and computing $h(t_n)$ from it by solving the equation above over "
+ "the time interval $t_n-t_{n-1}$. From this, one can then compute "
+ "a surface velocity $v = \\frac{h(t_n)-h(t_{n-1})}{t_n-t_{n-1}}$. "
+ "\n\n"
+ "This surface velocity is used to deform the surface and as a boundary condition "
+ "for solving the Laplace equation to determine the mesh velocity in the "
+ "domain interior. "
+ "Diffusion can be applied every timestep, mimicking surface processes of erosion "
+ "and deposition, or at a user-defined timestep interval to purely smooth the surface "
+ "topography to avoid too great a distortion of mesh elements when a free "
+ "surface is also used.")
+ }
+}