Merge pull request #5902 from danieldouglas92/fluid_velocity_postproc…

…essor

Add fluid velocity postprocessor

doc/modules/changes/20240613_danieldouglas92

@@ -0,0 +1,7 @@
+Added: A postprocessor 'darcy velocity' that visualizes the fluid
+velocity in models with a compositional field named 'porosity'. 
+This postprocessor is independent of the method the porosity is actually 
+advected with, i.e. it could be static, advected with the Darcy velocity, 
+or advected according to the two-phase flow equations.
+<br>
+(Daniel Douglas, 2024/06/13)

include/aspect/postprocess/visualization/darcy_velocity.h

@@ -0,0 +1,57 @@
+/*
+  Copyright (C) 2011 - 2024 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/>.
+*/
+
+
+#ifndef _aspect_postprocess_visualization_darcy_velocity_h
+#define _aspect_postprocess_visualization_darcy_velocity_h
+
+#include <aspect/postprocess/visualization.h>
+#include <aspect/simulator_access.h>
+
+#include <deal.II/numerics/data_postprocessor.h>
+
+
+namespace aspect
+{
+  namespace Postprocess
+  {
+    namespace VisualizationPostprocessors
+    {
+      template <int dim>
+      class DarcyVelocity
+        : public DataPostprocessorVector<dim>,
+          public SimulatorAccess<dim>,
+          public Interface<dim>
+      {
+        public:
+          DarcyVelocity ();
+
+          void
+          evaluate_vector_field(const DataPostprocessorInputs::Vector<dim> &input_data,
+                                std::vector<Vector<double>> &computed_quantities) const override;
+
+          std::string
+          get_physical_units () const override;
+      };
+    }
+  }
+}
+
+#endif

source/postprocess/visualization/darcy_velocity.cc

@@ -0,0 +1,120 @@
+/*
+  Copyright (C) 2011 - 2024 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/visualization/darcy_velocity.h>
+#include <aspect/utilities.h>
+#include <aspect/melt.h>
+#include <aspect/simulator.h>
+
+namespace aspect
+{
+  namespace Postprocess
+  {
+    namespace VisualizationPostprocessors
+    {
+      template <int dim>
+      DarcyVelocity<dim>::
+      DarcyVelocity ()
+        :
+        DataPostprocessorVector<dim> ("darcy_velocity",
+                                      update_values | update_quadrature_points | update_gradients),
+        Interface<dim>()
+      {}
+
+
+
+      template <int dim>
+      std::string
+      DarcyVelocity<dim>::
+      get_physical_units () const
+      {
+        if (this->convert_output_to_years())
+          return "m/year";
+        else
+          return "m/s";
+      }
+
+
+
+      template <int dim>
+      void
+      DarcyVelocity<dim>::
+      evaluate_vector_field(const DataPostprocessorInputs::Vector<dim> &input_data,
+                            std::vector<Vector<double>> &computed_quantities) const
+      {
+        AssertThrow(this->introspection().compositional_name_exists("porosity"),
+                    ExcMessage("The 'darcy velocity' postprocessor requires a compositional "
+                               "field called porosity."));
+
+        const unsigned int porosity_idx = this->introspection().find_composition_type(CompositionalFieldDescription::porosity);
+
+        const double velocity_scaling_factor =
+          this->convert_output_to_years() ? year_in_seconds : 1.0;
+
+        MaterialModel::MaterialModelInputs<dim> in(input_data, this->introspection());
+        MaterialModel::MaterialModelOutputs<dim> out(in.n_evaluation_points(), this->n_compositional_fields());
+        MeltHandler<dim>::create_material_model_outputs(out);
+        this->get_material_model().evaluate(in, out);
+        MaterialModel::MeltOutputs<dim> *fluid_out = out.template get_additional_output<MaterialModel::MeltOutputs<dim>>();
+        AssertThrow(std::isfinite(fluid_out->fluid_viscosities[0]),
+                    ExcMessage("To compute the Darcy velocity the material model needs to provide the melt material model "
+                               "outputs. At least the fluid viscosity was not computed, or is not a number."));
+
+        for (unsigned int q=0; q<in.n_evaluation_points(); ++q)
+          {
+            const double porosity = std::max(in.composition[q][porosity_idx], 1e-10);
+            const Tensor<1,dim> gravity = this->get_gravity_model().gravity_vector(in.position[q]);
+            const double solid_density = out.densities[q];
+            const double fluid_viscosity = fluid_out->fluid_viscosities[q];
+            const double fluid_density = fluid_out->fluid_densities[q];
+            const double permeability = fluid_out->permeabilities[q];
+            const Tensor<1,dim> solid_velocity = in.velocity[q];
+            const Tensor<1,dim> darcy_velocity = solid_velocity -
+                                                 permeability / fluid_viscosity / porosity * gravity *
+                                                 (solid_density - fluid_density) * velocity_scaling_factor;
+
+            for (unsigned int k=0; k<dim; ++k)
+              computed_quantities[q](k) = darcy_velocity[k];
+          }
+      }
+    }
+  }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+  namespace Postprocess
+  {
+    namespace VisualizationPostprocessors
+    {
+      ASPECT_REGISTER_VISUALIZATION_POSTPROCESSOR(DarcyVelocity,
+                                                  "darcy velocity",
+                                                  "A visualization output object that outputs the Darcy velocity "
+                                                  "vector. This postprocessor requires a compositional field named "
+                                                  "'porosity'."
+                                                  "\n\n"
+                                                  "Physical units: $\\frac{\\text{m}}{\\text{s}}$ or "
+                                                  "$\\frac{\\text{m}}{\\text{year}}$, depending on settings in the input file.")
+    }
+  }
+}

tests/darcy_convection_step.prm

@@ -27,7 +27,8 @@
 # pemeability is 1e-6, and gravity is set to -10. This results in a fluid
 # velocity of 24.25 m/yr. The mesh is 2.5 km x 2.5 km, and with CFL = 1,
 # this means that the time step should be 2500 m / 24.25 m/yr = 103.11
-
+# This also serves as a test for the 'darcy velocity' postprocessor, by
+# confirming that the darcy velocity is 24.25 m/yr.
 ############### Global parameters
 
 set Dimension                              = 2
@@ -154,6 +155,11 @@ subsection Melt settings
   set Include melt transport             = false
 end
 
+# Output the darcy velocity
 subsection Postprocess
-  set List of postprocessors =
+  set List of postprocessors = visualization
+  subsection Visualization
+    set List of output variables = darcy velocity
+    set Output format            = gnuplot
+  end
 end

tests/darcy_convection_step/log.txt

@@ -0,0 +1,44 @@
+
+Number of active cells: 4 (on 2 levels)
+Number of degrees of freedom: 134 (50+9+25+25+25)
+
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+   Solving porosity system ... 0 iterations.
+   Skipping bound_fluid composition solve because RHS is zero.
+   Solving Stokes system... 7+0 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 6.2054e-17, 1.34414e-16, 0, 5.09696e-16
+      Relative nonlinear residual (total system) after nonlinear iteration 1: 5.09696e-16
+
+
+   Postprocessing:
+     Writing graphical output: output-darcy_convection_step/solution/solution-00000
+
+*** Timestep 1:  t=103.11 years, dt=103.11 years
+   Solving composition reactions... in 1 substep(s).
+   Solving temperature system... 0 iterations.
+   Solving porosity system ... 0 iterations.
+   Skipping bound_fluid composition solve because RHS is zero.
+   Solving Stokes system... 6+0 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 3.51395e-16, 1.5726e-16, 0, 2.2185e-16
+      Relative nonlinear residual (total system) after nonlinear iteration 1: 3.51395e-16
+
+
+   Postprocessing:
+
+*** Timestep 2:  t=200 years, dt=96.8895 years
+   Solving composition reactions... in 1 substep(s).
+   Solving temperature system... 0 iterations.
+   Solving porosity system ... 0 iterations.
+   Skipping bound_fluid composition solve because RHS is zero.
+   Solving Stokes system... 7+0 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.14854e-16, 1.32969e-16, 0, 3.48944e-16
+      Relative nonlinear residual (total system) after nonlinear iteration 1: 3.48944e-16
+
+
+   Postprocessing:
+
+Termination requested by criterion: end time
+
+
+

tests/darcy_convection_step/screen-output

@@ -12,6 +12,7 @@ Number of degrees of freedom: 134 (50+9+25+25+25)
 
 
    Postprocessing:
+     Writing graphical output: output-darcy_convection_step/solution/solution-00000
 
 *** Timestep 1:  t=103.11 years, dt=103.11 years
    Solving composition reactions... in 1 substep(s).

tests/darcy_convection_step/statistics

@@ -12,6 +12,7 @@
 # 12: Iterations for Stokes solver
 # 13: Velocity iterations in Stokes preconditioner
 # 14: Schur complement iterations in Stokes preconditioner
-0 0.000000000000e+00 0.000000000000e+00 4 59 25 50 1 0 0 0 6 8 8 
-1 1.031104829590e+02 1.031104829590e+02 4 59 25 50 1 0 0 0 5 7 7 
-2 2.000000000000e+02 9.688951704104e+01 4 59 25 50 1 0 0 0 6 8 8 
+# 15: Visualization file name
+0 0.000000000000e+00 0.000000000000e+00 4 59 25 50 1 0 0 0 6 8 8 output-darcy_convection_step/solution/solution-00000 
+1 1.031104829590e+02 1.031104829590e+02 4 59 25 50 1 0 0 0 5 7 7                                                   "" 
+2 2.000000000000e+02 9.688951704104e+01 4 59 25 50 1 0 0 0 6 8 8                                                   ""