Merge pull request #6374 from alarshi/modify_depths_chunk

Modify depths in presence of topography in chunk.

Author: gassmoeller

include/aspect/geometry_model/chunk.h

@@ -119,6 +119,11 @@ namespace aspect
           std::unique_ptr<Manifold<dim,dim>>
           clone() const override;
 
+          /**
+           * Return the topography at a given point
+           */
+          double topography_for_point(const Point<dim> &x_y_z) const;
+
         private:
           /**
            * A pointer to the topography model.

source/geometry_model/chunk.cc

@@ -433,6 +433,34 @@ namespace aspect
 
         return r_phi_theta;
       }
+
+
+
+      template <int dim>
+      double
+      ChunkGeometry<dim>::
+      topography_for_point(const Point<dim> &x_y_z) const
+      {
+        if (dynamic_cast<const InitialTopographyModel::ZeroTopography<dim>*>(topo) != nullptr)
+          return 0;
+        else
+          {
+            // The natural coordinate system in chunk geometry is r/lon/lat.
+            // So, after converting the point into the spherical system we need to
+            // change the theta (co-latitudes) into latitudes.
+            const std::array<double, dim> r_phi_theta = Utilities::Coordinates::cartesian_to_spherical_coordinates(x_y_z);
+
+            // Grab lon,lat coordinates
+            Point<dim-1> surface_point;
+            // Convert latitude to colatitude
+            if (dim == 3)
+              surface_point[1] = 0.5*numbers::PI - surface_point[1];
+
+            for (unsigned int d=0; d<dim-1; ++d)
+              surface_point[d] = r_phi_theta[d+1];
+            return topo->value(surface_point);
+          }
+      }
     }
 
 
@@ -555,8 +583,12 @@ namespace aspect
     Chunk<dim>::depth(const Point<dim> &position) const
     {
       // depth is defined wrt the reference surface point2[0]
-      // negative depth is not allowed
-      return std::max (0., std::min (point2[0]-position.norm(), maximal_depth()));
+      // plus initial topography. Negative depth is not allowed.
+      if (this->simulator_is_past_initialization() &&
+          !Plugins::plugin_type_matches<const InitialTopographyModel::ZeroTopography<dim>>(this->get_initial_topography_model()))
+        return std::min(std::max(point2[0]+ manifold->topography_for_point(position) - position.norm(), 0.), maximal_depth());
+      else
+        return std::min(std::max(point2[0]-position.norm(), 0.), maximal_depth());
     }
 
 
@@ -582,10 +614,9 @@ namespace aspect
       // Choose a point at the mean longitude (and latitude)
       Point<dim> p = 0.5*(point2+point1);
       // at a depth beneath the top surface
-      p[0] = point2[0]-depth;
+      p[0] = point2[0]+manifold->topography_for_point(p) - depth;
 
-      // Now convert to Cartesian coordinates. This ignores the surface topography,
-      // but that is as documented.
+      // Now convert to Cartesian coordinates.
       return manifold->push_forward_sphere(p);
     }
 
@@ -658,12 +689,12 @@ namespace aspect
     double
     Chunk<dim>::maximal_depth() const
     {
-      // The depth is defined as relative to a reference surface (without
+      // The depth is defined as relative to a reference surface (with
       // topography) and since we don't apply topography on the CMB,
       // the maximal depth really is the formula below, unless one applies a
       // topography that is always strictly below zero (i.e., where the
       // actual surface lies strictly below the reference surface).
-      return point2[0]-point1[0];
+      return point2[0] + this->get_initial_topography_model().max_topography() - point1[0];
     }
 
 

tests/chunk_initial_topography_ascii_data/screen-output

@@ -7,11 +7,11 @@ Number of degrees of freedom: 3,556 (2,178+289+1,089)
 
 *** Timestep 0:  t=0 years, dt=0 years
    Solving temperature system... 0 iterations.
-   Solving Stokes system (GMG)... 29+0 iterations.
+   Solving Stokes system (GMG)... 28+0 iterations.
 
    Postprocessing:
 
-     Model domain depth (m):                        3e+06
+     Model domain depth (m):                        3.09e+06
      Temperature contrast across model domain (K):  3500
      Reference depth (m):                           0
      Reference temperature (K):                     1613
@@ -23,7 +23,7 @@ Number of degrees of freedom: 3,556 (2,178+289+1,089)
      Reference thermal conductivity (W/(m*K)):      4.7
      Reference viscosity (Pa*s):                    1e+21
      Reference thermal diffusivity (m^2/s):         1.22101e-06
-     Rayleigh number:                               7.4479e+07
+     Rayleigh number:                               8.13852e+07
 
      RMS, max velocity:       0.0119 m/year, 0.0933 m/year
      Temperature min/avg/max: 1613 K, 1613 K, 1613 K

tests/chunk_initial_topography_ascii_data/statistics

@@ -16,4 +16,4 @@
 # 16: Average nondimensional temperature (K)
 # 17: Minimum topography (m)
 # 18: Maximum topography (m)
-0 0.000000000000e+00 0.000000000000e+00 256 2467 1089 0 29 30 30 1.18604659e-02 9.33132004e-02 1.61300000e+03 1.61300000e+03 1.61300000e+03 4.40406633e-10 -4.75663606e+04 9.00000000e+04 
+0 0.000000000000e+00 0.000000000000e+00 256 2467 1089 0 28 29 29 1.18604723e-02 9.33132938e-02 1.61300000e+03 1.61300000e+03 1.61300000e+03 4.40406633e-10 -4.75663606e+04 9.00000000e+04 

tests/chunk_initial_topography_ascii_data_3d/screen-output

@@ -11,7 +11,7 @@ Number of degrees of freedom: 3,041 (2,187+125+729)
 
    Postprocessing:
 
-     Model domain depth (m):                        3e+06
+     Model domain depth (m):                        3.009e+06
      Temperature contrast across model domain (K):  3500
      Reference depth (m):                           0
      Reference temperature (K):                     1613
@@ -23,7 +23,7 @@ Number of degrees of freedom: 3,041 (2,187+125+729)
      Reference thermal conductivity (W/(m*K)):      4.7
      Reference viscosity (Pa*s):                    1e+21
      Reference thermal diffusivity (m^2/s):         1.22101e-06
-     Rayleigh number:                               7.4479e+07
+     Rayleigh number:                               7.51513e+07
 
      RMS, max velocity:       0.0442 m/year, 0.113 m/year
      Temperature min/avg/max: 1613 K, 1613 K, 1613 K

tests/chunk_initial_topography_ascii_data_3d/statistics

@@ -16,4 +16,4 @@
 # 16: Average nondimensional temperature (K)
 # 17: Minimum topography (m)
 # 18: Maximum topography (m)
-0 0.000000000000e+00 0.000000000000e+00 64 2312 729 0 23 24 24 4.42347040e-02 1.13166975e-01 1.61300000e+03 1.61300006e+03 1.61300000e+03 1.71059328e-08 -1.98131486e+03 8.99532925e+03 
+0 0.000000000000e+00 0.000000000000e+00 64 2312 729 0 23 24 24 4.42347089e-02 1.13166977e-01 1.61300000e+03 1.61300006e+03 1.61300000e+03 1.71059328e-08 -1.98131486e+03 8.99532925e+03 

tests/chunk_initial_topography_ascii_data_3d_colatitude/screen-output

@@ -7,11 +7,11 @@ Number of degrees of freedom: 3,041 (2,187+125+729)
 
 *** Timestep 0:  t=0 years, dt=0 years
    Solving temperature system... 0 iterations.
-   Solving Stokes system (GMG)... 21+0 iterations.
+   Solving Stokes system (GMG)... 20+0 iterations.
 
    Postprocessing:
 
-     Model domain depth (m):                        3e+06
+     Model domain depth (m):                        3.009e+06
      Temperature contrast across model domain (K):  3500
      Reference depth (m):                           0
      Reference temperature (K):                     1613
@@ -23,7 +23,7 @@ Number of degrees of freedom: 3,041 (2,187+125+729)
      Reference thermal conductivity (W/(m*K)):      4.7
      Reference viscosity (Pa*s):                    1e+21
      Reference thermal diffusivity (m^2/s):         1.22101e-06
-     Rayleigh number:                               7.4479e+07
+     Rayleigh number:                               7.51513e+07
 
      RMS, max velocity:       0.0124 m/year, 0.0606 m/year
      Temperature min/avg/max: 1613 K, 1613 K, 1613 K

tests/chunk_initial_topography_ascii_data_3d_colatitude/statistics

@@ -16,4 +16,4 @@
 # 16: Average nondimensional temperature (K)
 # 17: Minimum topography (m)
 # 18: Maximum topography (m)
-0 0.000000000000e+00 0.000000000000e+00 64 2312 729 0 21 22 22 1.24375134e-02 6.06479500e-02 1.61300000e+03 1.61300000e+03 1.61300000e+03 -1.75439514e-10 -1.95928291e+03 8.99918928e+03 
+0 0.000000000000e+00 0.000000000000e+00 64 2312 729 0 20 21 21 1.24375096e-02 6.06479458e-02 1.61300000e+03 1.61300000e+03 1.61300000e+03 -1.75439514e-10 -1.95928291e+03 8.99918928e+03 

tests/chunk_lithostatic_pressure_2d_initial_topography/screen-output

@@ -11,8 +11,8 @@ Number of degrees of freedom: 109 (50+9+25+25)
    Solving Stokes system (GMG)... 13+0 iterations.
 
    Postprocessing:
-     RMS, max velocity:        0.261 m/year, 0.401 m/year
-     Pressure min/avg/max:     -2.139e+09 Pa, 3.927e+10 Pa, 9.249e+10 Pa
+     RMS, max velocity:        0.27 m/year, 0.427 m/year
+     Pressure min/avg/max:     -2.152e+09 Pa, 3.928e+10 Pa, 9.249e+10 Pa
      Writing graphical output: output-chunk_lithostatic_pressure_2d_initial_topography/solution/solution-00000
 
 Termination requested by criterion: end time