Initial Landlab coupling

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
+  <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifndef _aspect_mesh_deformation_landlab_h
+#define _aspect_mesh_deformation_landlab_h
+
+#include <aspect/mesh_deformation/interface.h>
+#include <aspect/simulator_access.h>
+#include <aspect/geometry_model/initial_topography_model/interface.h>
+
+#include <Python.h>
+
+namespace aspect
+{
+  namespace MeshDeformation
+  {
+    /**
+     * A plugin that computes the deformation of surface
+     * vertices by coupling with LandLab.
+     *
+     * @ingroup MeshDeformation
+     */
+    template <int dim>
+    class LandLab : public Interface<dim>, public SimulatorAccess<dim>
+    {
+      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<dim> &mesh_deformation_dof_handler,
+                                                 AffineConstraints<double> &mesh_velocity_constraints,
+                                                 const std::set<types::boundary_id> &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<dim> &free_surface_dof_handler,
+                               const IndexSet &mesh_locally_owned,
+                               const IndexSet &mesh_locally_relevant,
+                               LinearAlgebra::Vector &output,
+                               const std::set<types::boundary_id> &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<dim> &mesh_deformation_dof_handler,
+                                        const std::set<types::boundary_id> &boundary_ids) const;
+
+
+        /**
+        */
+        std::vector<Point<dim>> 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<types::boundary_id> additional_tangential_mesh_boundary_indicators;
+
+        PyObject *pModule;
+    };
+  }
+}
+
+
+#endif

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

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