Merge pull request #5662 from gassmoeller/add_coordinate_system_to_he…

…ating_function

Add coordinate system to heating function

include/aspect/heating_model/function.h

@@ -83,6 +83,12 @@ namespace aspect
          * A function object representing the components of the velocity.
          */
         Functions::ParsedFunction<dim> heating_model_function;
+
+        /**
+         * The coordinate representation to evaluate the function. Possible
+         * choices are depth, cartesian and spherical.
+         */
+        Utilities::Coordinates::CoordinateSystem coordinate_system;
     };
   }
 }

source/heating_model/function.cc

@@ -20,7 +20,8 @@
 
 
 #include <aspect/heating_model/function.h>
-#include <aspect/global.h>
+
+#include <aspect/geometry_model/interface.h>
 
 namespace aspect
 {
@@ -42,9 +43,13 @@ namespace aspect
     {
       for (unsigned int q=0; q<heating_model_outputs.heating_source_terms.size(); ++q)
         {
-          // return a constant value
+          // convert the position into the selected coordinate system
           const Point<dim> position = material_model_inputs.position[q];
-          heating_model_outputs.heating_source_terms[q] = heating_model_function.value(position)
+          const Utilities::NaturalCoordinate<dim> point =
+            this->get_geometry_model().cartesian_to_other_coordinates(position, coordinate_system);
+
+          // then compute the heating function value at this position
+          heating_model_outputs.heating_source_terms[q] = heating_model_function.value(Utilities::convert_array_to_point<dim>(point.get_coordinates()))
                                                           * material_model_outputs.densities[q];
           heating_model_outputs.lhs_latent_heat_terms[q] = 0.0;
         }
@@ -73,6 +78,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, 1);
         }
         prm.leave_subsection();
@@ -88,21 +104,23 @@ namespace aspect
       prm.enter_subsection("Heating model");
       {
         prm.enter_subsection("Function");
-        try
-          {
-            heating_model_function.parse_parameters (prm);
-          }
-        catch (...)
-          {
-            std::cerr << "ERROR: FunctionParser failed to parse\n"
-                      << "\t'Heating model.Function'\n"
-                      << "with expression\n"
-                      << "\t'" << prm.get("Function expression") << "'"
-                      << "More information about the cause of the parse error \n"
-                      << "is shown below.\n";
-            throw;
-          }
         {
+          try
+            {
+              heating_model_function.parse_parameters (prm);
+            }
+          catch (...)
+            {
+              std::cerr << "ERROR: FunctionParser failed to parse\n"
+                        << "\t'Heating model.Function'\n"
+                        << "with expression\n"
+                        << "\t'" << prm.get("Function expression") << "'"
+                        << "More information about the cause of the parse error \n"
+                        << "is shown below.\n";
+              throw;
+            }
+
+          coordinate_system = Utilities::Coordinates::string_to_coordinate_system(prm.get("Coordinate system"));
         }
         prm.leave_subsection();
       }

tests/heating_function.prm

@@ -0,0 +1,67 @@
+# A test for the function heating plugin.
+# The heating function is given in cartesian coordinates.
+
+set Dimension                              = 2
+set End time                               = 0.0
+set Use years in output instead of seconds = false
+
+subsection Geometry model
+  set Model name = spherical shell
+
+  subsection Spherical shell
+    set Inner radius = 0.55
+    set Outer radius  = 1.0
+    set Opening angle = 360
+  end
+end
+
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = bottom, top
+end
+
+subsection Material model
+  set Model name = simple
+
+  subsection Simple model
+    set Reference density             = 1
+    set Viscosity                     = 1
+    set Thermal expansion coefficient = 0
+  end
+end
+
+subsection Gravity model
+  set Model name = radial constant
+end
+
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 0
+  end
+end
+
+subsection Heating model
+  set List of model names = function
+
+  subsection Function
+    set Function expression = x
+  end
+end
+
+subsection Mesh refinement
+  set Initial adaptive refinement        = 0
+  set Strategy                           = thermal energy density
+  set Initial global refinement          = 2
+  set Time steps between mesh refinement = 0
+end
+
+subsection Postprocess
+  set List of postprocessors = heating statistics, visualization
+
+  subsection Visualization
+    set List of output variables = heating
+    set Output format = gnuplot
+    set Time between graphical output = 0.0
+  end
+end

tests/heating_function/screen-output

@@ -0,0 +1,17 @@
+
+Number of active cells: 192 (on 3 levels)
+Number of degrees of freedom: 2,832 (1,728+240+864)
+
+*** Timestep 0:  t=0 seconds, dt=0 seconds
+   Skipping temperature solve because RHS is zero.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system... 34+0 iterations.
+
+   Postprocessing:
+     Heating rate (average/total):  3.083e-17 W/kg, 6.755e-17 W
+     Writing graphical output:      output-heating_function/solution/solution-00000
+
+Termination requested by criterion: end time
+
+
+