Merge pull request #5752 from KerrMadeleine/grav_model_coords

Added coordinate systems to Gravity model function plug-in

doc/modules/changes/20240602_KerrMadeleine

@@ -0,0 +1,3 @@
+Added: The Function expression plug-in for the Gravity model can now evaluate functions in multiple coordinate systems. The plug-in was amended to read in entries to the "Coordinate system" parameter in the Function subsection and Gravity Model section of the parameter file. There are three added tests for the Gravity model function plug-in to ensure the cartesian, spherical, and depth coordinate systems evaluate functions as expected in a 2D cylindrical annulus.
+<br>
+(Madeleine Kerr, 2024/06/02)

include/aspect/gravity_model/function.h

@@ -24,6 +24,7 @@
 
 #include <aspect/gravity_model/interface.h>
 #include <aspect/simulator_access.h>
+#include <aspect/utilities.h>
 
 #include <deal.II/base/parsed_function.h>
 
@@ -85,6 +86,12 @@ namespace aspect
          * A function object representing the gravity.
          */
         Functions::ParsedFunction<dim> function;
+
+        /**
+         * The coordinate representation to evaluate the function. Possible
+         * choices are depth, cartesian and spherical.
+         */
+        Utilities::Coordinates::CoordinateSystem coordinate_system;
     };
   }
 }

source/gravity_model/function.cc

@@ -20,6 +20,7 @@
 
 
 #include <aspect/gravity_model/function.h>
+#include <aspect/geometry_model/interface.h>
 
 #include <iostream>
 
@@ -39,8 +40,11 @@ namespace aspect
     gravity_vector (const Point<dim> &position) const
     {
       Tensor<1,dim> gravity;
+      const Utilities::NaturalCoordinate<dim> point =
+        this->get_geometry_model().cartesian_to_other_coordinates(position, coordinate_system);
+      const Point<dim> point_in_coordinate_system = Utilities::convert_array_to_point<dim>(point.get_coordinates());
       for (unsigned int d=0; d<dim; ++d)
-        gravity[d] = function.value(position,d);
+        gravity[d] = function.value(point_in_coordinate_system,d);
       return gravity;
     }
 
@@ -66,6 +70,17 @@ namespace aspect
       {
         prm.enter_subsection("Function");
         {
+          prm.declare_entry ("Coordinate system", "cartesian",
+                             Patterns::Selection ("cartesian|spherical|depth"),
+                             "A selection that determines the assumed coordinate "
+                             "system for the function variables. Allowed values "
+                             "are `cartesian', `spherical', and `depth'. `spherical' coordinates "
+                             "are interpreted as r,phi or r,phi,theta in 2d/3d "
+                             "respectively with theta being the polar angle. `depth' "
+                             "will create a function, in which only the first "
+                             "parameter is non-zero, which is interpreted to "
+                             "be the depth of the point.");
+
           Functions::ParsedFunction<dim>::declare_parameters (prm, dim);
         }
         prm.leave_subsection();
@@ -81,6 +96,9 @@ namespace aspect
       prm.enter_subsection("Gravity model");
       {
         prm.enter_subsection("Function");
+        {
+          coordinate_system = Utilities::Coordinates::string_to_coordinate_system(prm.get("Coordinate system"));
+        }
         try
           {
             function.parse_parameters (prm);

tests/gravity_model_cartesian.prm

@@ -0,0 +1,98 @@
+# This is a test based on benchmarks/onset-of-convection/convection-box-base.prm.
+# It makes sure that the gravity plugin 'function' works correctly
+# for a Cartesian coordinate system.
+
+set Dimension                              = 2
+set Use years in output instead of seconds = true
+set Output directory                       = output
+set Pressure normalization                 = surface
+set Surface pressure                       = 0
+set Use conduction timestep                = true
+
+subsection Termination criteria
+  set Termination criteria                 = end step
+  set End step                             = 0
+end
+
+subsection Formulation
+  set Formulation = Boussinesq approximation
+end
+
+subsection Geometry model
+  set Model name = spherical shell
+
+  subsection Spherical shell
+    set Outer radius = 6371000
+    set Inner radius = 3481000
+  end
+end
+
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Coordinate system = spherical
+    set Variable names = r,phi
+    set Function expression = 2500 + 10*sin(10*phi)
+  end
+end
+
+subsection Boundary temperature model
+  set Fixed temperature boundary indicators = bottom, top
+  set List of model names = spherical constant
+
+  subsection Spherical constant
+    set Inner temperature = 3500
+    set Outer temperature = 0
+  end
+end
+
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = bottom, top
+end
+
+#subsection Gravity model
+#  set Model name = radial constant
+#
+#  subsection Radial constant
+#    set Magnitude = 10.0
+#  end
+#end
+
+subsection Gravity model
+  set Model name = function
+  subsection Function
+    set Coordinate system = cartesian
+    set Variable names = x,y
+    set Function expression = x;y
+  end
+end
+
+subsection Material model
+  set Model name = simple
+
+  subsection Simple model
+    set Reference density             = 4000
+    set Reference specific heat       = 1250
+    set Reference temperature         = 0
+    set Thermal conductivity          = 4.0
+    set Thermal expansion coefficient = 3e-5
+    set Viscosity                     = 1e23
+  end
+end
+
+subsection Mesh refinement
+  set Initial global refinement                = 0
+  set Initial adaptive refinement              = 0
+  set Time steps between mesh refinement       = 0
+end
+
+subsection Postprocess
+  set List of postprocessors = velocity statistics, visualization
+
+  subsection Visualization
+    set Time steps between graphical output = 5
+    set List of output variables = gravity
+    set Output format = gnuplot
+  end
+end

tests/gravity_model_cartesian/screen-output

@@ -0,0 +1,17 @@
+
+Number of active cells: 12 (on 1 levels)
+Number of degrees of freedom: 240 (144+24+72)
+
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system... 13+0 iterations.
+
+   Postprocessing:
+     RMS, max velocity:        98.6 m/year, 201 m/year
+     Writing graphical output: output-gravity_model_cartesian/solution/solution-00000
+
+Termination requested by criterion: end step
+
+
+

tests/gravity_model_cartesian/statistics

@@ -0,0 +1,14 @@
+# 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: Iterations for temperature solver
+# 8: Iterations for Stokes solver
+# 9: Velocity iterations in Stokes preconditioner
+# 10: Schur complement iterations in Stokes preconditioner
+# 11: RMS velocity (m/year)
+# 12: Max. velocity (m/year)
+# 13: Visualization file name
+0 0.000000000000e+00 0.000000000000e+00 12 168 72 0 12 14 39 9.86096651e+01 2.00683097e+02 output-gravity_model_cartesian/solution/solution-00000