Update point_is_in_domain function

source/geometry_model/box.cc

@@ -21,6 +21,7 @@
 
 #include <aspect/geometry_model/box.h>
 #include <aspect/geometry_model/initial_topography_model/zero_topography.h>
+#include <aspect/mesh_deformation/interface.h>
 #include <aspect/simulator_signals.h>
 
 #include <deal.II/grid/grid_generator.h>
@@ -316,20 +317,77 @@ namespace aspect
     bool
     Box<dim>::point_is_in_domain(const Point<dim> &point) const
     {
-      AssertThrow(!this->get_parameters().mesh_deformation_enabled ||
-                  this->simulator_is_past_initialization() == false,
-                  ExcMessage("After displacement of the free surface, this function can no longer be used to determine whether a point lies in the domain or not."));
+      // If mesh deformation occurs over time, and we're past the
+      // first timestep (during which no mesh deformation is applied),
+      // we can no longer use the original extents of the model domain
+      // to check whether the point lies in the domain.
+      // Initial mesh deformation is applied during the first
+      // timestep, but this function might be called before the
+      // simulator is initialized and thus before we can ask whether initial
+      // or time-dependent mesh deformation is applied.
+      // Therefore, when any mesh deformation is enabled,
+      // we loop over all the current grid cells to see if the point
+      // lies within the domain.
+      // Do note that this function can be called before mesh deformation
+      // is applied in the first timestep, and therefore there is no guarantee
+      // that the point will lie in the domain after initial mesh deformation.
+      // For example, boundary traction plugins are initialized before
+      // mesh deformation plugins and the initial_lithostatic_pressure
+      // plugin calls this function while the simulator is being initialized.
+      // However, if the simulator has not been initialized,
+      // no mesh deformation has been applied and we can use the geometry model's
+      // extents to check whether a point lies in the domain.
+      if (this->get_parameters().mesh_deformation_enabled &&
+          this->simulator_is_past_initialization())
+        {
+          bool cell_found = false;
+          std::pair<const typename parallel::distributed::Triangulation<dim>::active_cell_iterator,
+              Point<dim>> it =
+                GridTools::find_active_cell_around_point<> (this->get_mapping(), this->get_triangulation(), point);
+
+          // If the cell the point is in is on this processor, we have found the right cell.
+          if (it.first.state() == IteratorState::valid && it.first->is_locally_owned())
+            cell_found = true;
+
+          // Compute how many processes found the cell.
+          const int n_procs_cell_found = Utilities::MPI::sum(cell_found ? 1 : 0, this->get_mpi_communicator());
+
+          // If at least one process found the cell, the point is in the domain.
+          if (n_procs_cell_found > 0)
+            return true;
+          else
+            return false;
+        }
+      // Without mesh deformation enabled, it is much cheaper to check whether the point lies in the domain.
+      else
+        {
 
-      AssertThrow(Plugins::plugin_type_matches<const InitialTopographyModel::ZeroTopography<dim>>(this->get_initial_topography_model()),
-                  ExcMessage("After adding topography, this function can no longer be used "
-                             "to determine whether a point lies in the domain or not."));
+          // The maximal extents of the unperturbed box domain.
+          Point<dim> max_point = extents+box_origin;
 
-      for (unsigned int d = 0; d < dim; ++d)
-        if (point[d] > extents[d]+box_origin[d]+std::numeric_limits<double>::epsilon()*extents[d] ||
-            point[d] < box_origin[d]-std::numeric_limits<double>::epsilon()*extents[d])
-          return false;
+          // If mesh deformation is not enabled, but initial topography
+          // was/will be applied to the mesh, include this topography in the
+          // extent of the domain.
+          if (!Plugins::plugin_type_matches<const InitialTopographyModel::ZeroTopography<dim>>(this->get_initial_topography_model()))
+            {
+              // Get the surface x (,y) point
+              Point<dim-1> surface_point;
+              for (unsigned int d=0; d<dim-1; ++d)
+                surface_point[d] = point[d];
+
+              // Get the surface topography at this point
+              const double topo = topo_model->value(surface_point);
+              max_point[dim-1] += topo;
+            }
+
+          // Check whether point lies within the min/max coordinates of the domain including initial topography.
+          for (unsigned int d = 0; d < dim; ++d)
+            if (point[d] > max_point[d]+std::numeric_limits<double>::epsilon()*extents[d] ||
+                point[d] < box_origin[d]-std::numeric_limits<double>::epsilon()*extents[d])
+              return false;
 
-      return true;
+          return true;
+        }
     }
 
     template <int dim>

source/geometry_model/chunk.cc

@@ -697,21 +697,47 @@ namespace aspect
     bool
     Chunk<dim>::point_is_in_domain(const Point<dim> &point) const
     {
-      AssertThrow(!this->get_parameters().mesh_deformation_enabled ||
-                  this->simulator_is_past_initialization() == false,
-                  ExcMessage("After displacement of the free surface, this function can no longer be used to determine whether a point lies in the domain or not."));
-
-      AssertThrow(Plugins::plugin_type_matches<const InitialTopographyModel::ZeroTopography<dim>>(this->get_initial_topography_model()),
-                  ExcMessage("After adding topography, this function can no longer be used to determine whether a point lies in the domain or not."));
-
-      const Point<dim> spherical_point = manifold->pull_back(point);
+      // If mesh deformation has possibly been applied,
+      // e.g. initial mesh deformation after the simulator
+      // has been initialized or time-dependent mesh deformation
+      // after timestep 0, we cannot use information from the
+      // geometry model to determine whether the point lies in
+      // the model domain. Instead, we loop over all cells
+      // to check whether the point lies in one of them.
+      if (this->get_parameters().mesh_deformation_enabled &&
+          this->simulator_is_past_initialization())
+        {
+          bool cell_found = false;
+          std::pair<const typename parallel::distributed::Triangulation<dim>::active_cell_iterator,
+              Point<dim>>
+              it =
+                GridTools::find_active_cell_around_point<>(this->get_mapping(), this->get_triangulation(), point);
+
+          // If the cell the point is in is on this processor, we have found the right cell.
+          if (it.first.state() == IteratorState::valid && it.first->is_locally_owned())
+            cell_found = true;
+
+          // Compute how many processes found the cell.
+          const int n_procs_cell_found = Utilities::MPI::sum(cell_found ? 1 : 0, this->get_mpi_communicator());
+
+          // If at least one process found the cell, the point is in the domain.
+          if (n_procs_cell_found > 0)
+            return true;
+          else
+            return false;
+        }
+      // Without mesh deformation enabled, it is much cheaper to check whether the point lies in the domain.
+      else
+        {
+          const Point<dim> spherical_point = manifold->pull_back(point);
 
-      for (unsigned int d = 0; d < dim; ++d)
-        if (spherical_point[d] > point2[d]+std::numeric_limits<double>::epsilon()*std::abs(point2[d]) ||
-            spherical_point[d] < point1[d]-std::numeric_limits<double>::epsilon()*std::abs(point2[d]))
-          return false;
+          for (unsigned int d = 0; d < dim; ++d)
+            if (spherical_point[d] > point2[d]+std::numeric_limits<double>::epsilon()*std::abs(point2[d]) ||
+                spherical_point[d] < point1[d]-std::numeric_limits<double>::epsilon()*std::abs(point2[d]))
+              return false;
 
-      return true;
+          return true;
+        }
     }
 
 

tests/airy_isostasy_initial_topo.cc

@@ -0,0 +1,77 @@
+/*
+  Copyright (C) 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
+  <http://www.gnu.org/licenses/>.
+*/
+
+#include <aspect/postprocess/interface.h>
+#include <aspect/geometry_model/interface.h>
+#include <aspect/simulator_access.h>
+
+
+namespace aspect
+{
+  namespace Postprocess
+  {
+    template <int dim>
+    class PointIsInDomain : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
+    {
+      public:
+        virtual
+        std::pair<std::string,std::string>
+        execute (TableHandler &);
+    };
+  }
+}
+
+
+namespace aspect
+{
+  namespace Postprocess
+  {
+    template <int dim>
+    std::pair<std::string,std::string>
+    PointIsInDomain<dim>::execute(TableHandler &)
+    {
+      // The airy isostasy box domain is 1x1 with max 0.05 topography.
+      const Point<dim> point_in_domain(0.1, 0.1);
+      const Point<dim> point_not_in_domain(0.1, 2.0);
+      const Point<dim> point_not_in_domain_2(2.0, 0.01);
+
+      std::ostringstream screen_text;
+      screen_text << std::boolalpha << this->get_geometry_model().point_is_in_domain(point_in_domain) ;
+      screen_text << std::boolalpha << this->get_geometry_model().point_is_in_domain(point_not_in_domain);
+      screen_text << std::boolalpha << this->get_geometry_model().point_is_in_domain(point_not_in_domain_2);
+
+      return std::pair<std::string, std::string> ("Points lie in domain: ",
+                                                  screen_text.str());
+    }
+  }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+  namespace Postprocess
+  {
+    ASPECT_REGISTER_POSTPROCESSOR(PointIsInDomain,
+                                  "point is in domain",
+                                  "A postprocessor that tests whether certain points "
+                                  "lie within the domain.")
+  }
+}

tests/airy_isostasy_initial_topo.prm

@@ -30,3 +30,8 @@ subsection Boundary traction model
     set Function expression = 0;10
   end
 end
+
+# Also test the point_is_in_domain function of the geometry model.
+subsection Postprocess
+  set List of postprocessors = velocity statistics, pressure statistics, topography, point is in domain
+end

tests/airy_isostasy_initial_topo_lithostatic_pressure.prm

@@ -0,0 +1,40 @@
+# We setup three columns, of which the middle column has a higher density
+# and viscosity. By choosing a free surface in combination with a
+# prescribed lithostatic pressure,
+# we see the middle column sink until isostasy is reached (if end time is increased)
+# Compared to airy_isostasy.prm, we start off with a nonzero initial topography.
+# This test in particular tests the point_is_in_domain function of the geometry
+# model. It should take into account the initial topography and find that the
+# Representative point lies within the domain. At x = 0.2, the initial topography
+# is namely 0.02948.
+
+include $ASPECT_SOURCE_DIR/tests/airy_isostasy.prm
+
+set Dimension = 2
+# The point_is_in_domain function is only called in
+# timestep 0, so we do not need to run the test further.
+set End time                               = 0
+
+subsection Geometry model
+  subsection Initial topography model
+    set Model name = function
+
+    subsection Function
+      set Function constants = L=1, R=0.05, C=0.5
+      set Function expression = R * cos(pi*(x-C)/(L))
+      set Maximum topography value = 0.05
+    end
+  end
+end
+
+# The representative point determines where the pressure
+# profile is calculated. We pick a point in the left
+# column which has some initial topography prescribed.
+subsection Boundary traction model
+  set Prescribed traction boundary indicators = 2 y: initial lithostatic pressure
+
+  subsection Initial lithostatic pressure
+    set Number of integration points = 260
+    set Representative point         = 0.2,1.02
+  end
+end

tests/airy_isostasy_initial_topo_lithostatic_pressure/screen-output

@@ -0,0 +1,46 @@
+
+Number of active cells: 1,024 (on 6 levels)
+Number of degrees of freedom: 17,989 (8,450+1,089+4,225+4,225)
+
+Number of mesh deformation degrees of freedom: 2,178
+   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.
+   Solving C_1 system ... 0 iterations.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system... 200+8 iterations.
+
+Number of active cells: 1,204 (on 7 levels)
+Number of degrees of freedom: 21,360 (10,034+1,292+5,017+5,017)
+
+Number of mesh deformation degrees of freedom: 2,584
+   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.
+   Solving C_1 system ... 0 iterations.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system... 167+0 iterations.
+
+Number of active cells: 1,948 (on 8 levels)
+Number of degrees of freedom: 34,621 (16,270+2,081+8,135+8,135)
+
+Number of mesh deformation degrees of freedom: 4,162
+   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.
+   Solving C_1 system ... 0 iterations.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system... 200+6 iterations.
+
+   Postprocessing:
+     RMS, max velocity:    6.63e+05 m/year, 8.88e+05 m/year
+     Pressure min/avg/max: -0.6507 Pa, 5.005 Pa, 10 Pa
+     Topography min/max:   0 m, 0 m
+
+Termination requested by criterion: end time
+
+
+

tests/airy_isostasy_initial_topo_lithostatic_pressure/statistics

@@ -0,0 +1,20 @@
+# 1: Time step number
+# 2: Time (years)
+# 3: Time step size (years)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Number of degrees of freedom for all compositions
+# 8: Iterations for temperature solver
+# 9: Iterations for composition solver 1
+# 10: Iterations for Stokes solver
+# 11: Velocity iterations in Stokes preconditioner
+# 12: Schur complement iterations in Stokes preconditioner
+# 13: RMS velocity (m/year)
+# 14: Max. velocity (m/year)
+# 15: Minimal pressure (Pa)
+# 16: Average pressure (Pa)
+# 17: Maximal pressure (Pa)
+# 18: Minimum topography (m)
+# 19: Maximum topography (m)
+0 0.000000000000e+00 0.000000000000e+00 1948 18351 8135 8135 0 0 206 444 840 6.62545995e+05 8.87572504e+05 -6.50724076e-01 5.00501899e+00 1.00000000e+01 0.00000000e+00 0.00000000e+00