Merge pull request #6151 from bangerth/heat-flux

Document the heat flux plugin system.

doc/sphinx/user/extending/plugin-types/boundary-heat-flux.md

@@ -0,0 +1,34 @@
+
+# Boundary heat flux
+
+In contrast to prescribing the actual temperature at a boundary, it is
+also possible to prescribe the *heat flux* across a boundary. A number
+of heat flux models are already implemented, but if you want to
+implement a new boundary heat flux model, you need to overload the
+`aspect::BoundaryHeatFlux::Interface` class and use the
+`ASPECT_REGISTER_BOUNDARY_HEAT_FLUX_MODEL` macro to register your new class.
+The implementation of the new class should be in namespace
+`aspect::BoundaryTemperature`.
+
+The `Interface` class that you need to overload provides the
+prescribed heat flux at a given set of points via the `heat_flux()`
+function. The function receives information about material parameters
+at these points as well. For historical reasons, the function is asked
+to provide the heat flux as a vector, even though the place where the
+heat flux is used only uses the component of this vector that is
+*normal* to the boundary (which it computes by taking the dot product
+between the returned vector and the normal vector). Because there are
+situations where all you can do is compute the normal heat flux as a
+scalar, the `heat_flux()` function also receives the normal vector as
+an input argument. As a consequence, one way for the function to
+compute the required heat flux vector is to compute the scalar heat
+flux and multiply it by the normal vector.
+
+Finally, the function also receives the boundary indicator of the particular
+piece of boundary on which the point is located, as a hint in
+determining where this point may be located; this may, for example, be used to
+determine if a point is on the inner or outer boundary of a spherical shell.
+Models may also need to query the geometry model (via the `SimulatorAccess` class)
+to understand what a specific boundary indicator is supposed to mean – e.g.,
+if the given boundary indicator corresponds to an inner or outer surface of a shell
+geometry.

doc/sphinx/user/extending/plugin-types/index.md

@@ -9,6 +9,7 @@ gravity-models.md
 initial-conditions.md
 velocity-bc.md
 temp-bc.md
+boundary-heat-flux.md
 postprocessors.md
 vis-postprocessors.md
 mesh-refinement.md

include/aspect/boundary_heat_flux/interface.h

@@ -59,7 +59,18 @@ namespace aspect
          * @param material_model_inputs The material property inputs.
          * @param material_model_outputs The material property outputs.
          * @param normal_vectors The normal vector at each quadrature point.
+         *
          * @return A vector of heat flux vectors at the evaluation points.
+         *   For historical reasons, the function is asked
+         *   to provide the heat flux as a vector, even though the place where the
+         *   heat flux is used only uses the component of this vector that is
+         *   to the boundary (which it computes by taking the dot product *normal*
+         *   between the returned vector and the normal vector). Because there are
+         *   situations where all you can do is compute the normal heat flux as a
+         *   scalar, the `heat_flux()` function also receives the normal vector as
+         *   an input argument. As a consequence, one way for the function to
+         *   compute the required heat flux vector is to compute the scalar heat
+         *   flux and multiply it by the normal vector.
          */
         virtual
         std::vector<Tensor<1,dim>>

source/simulator/assemblers/advection.cc

@@ -461,13 +461,12 @@ namespace aspect
           // We are in the case of a Neumann temperature boundary.
           // Impose the Neumann value weakly using a RHS term.
 
-          std::vector<Tensor<1,dim>> heat_flux(n_face_q_points);
-          heat_flux = this->get_boundary_heat_flux().heat_flux(
-                        face->boundary_id(),
-                        scratch.face_material_model_inputs,
-                        scratch.face_material_model_outputs,
-                        scratch.face_finite_element_values->get_normal_vectors()
-                      );
+          const std::vector<Tensor<1,dim>> heat_flux
+            = this->get_boundary_heat_flux().heat_flux(
+                face->boundary_id(),
+                scratch.face_material_model_inputs,
+                scratch.face_material_model_outputs,
+                scratch.face_finite_element_values->get_normal_vectors());
 
           for (unsigned int q=0; q<n_face_q_points; ++q)
             {
@@ -479,7 +478,7 @@ namespace aspect
                 {
                   if (fe.system_to_component_index(i).first == solution_component)
                     {
-                      scratch.face_phi_field[i_advection]      = (*scratch.face_finite_element_values)[solution_field].value (i, q);
+                      scratch.face_phi_field[i_advection] = (*scratch.face_finite_element_values)[solution_field].value (i, q);
                       ++i_advection;
                     }
                   ++i;