[cig-commits] [commit] master: Add a postprocessor that computes statistics about the pressure. Add a simple Poisseuille flow test. (c1fc7e2)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Fri May 23 07:27:30 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/19b0c92b6eac4e5192d682ac769218eb3fa1db05...1d916725840b7e74c3611154de2dcabdfa3dfdce
>---------------------------------------------------------------
commit c1fc7e2ec50dcca1d024e26b9e85e2c4ab40147b
Author: Wolfgang Bangerth <bangerth at math.tamu.edu>
Date: Fri May 23 09:13:06 2014 -0500
Add a postprocessor that computes statistics about the pressure. Add a simple Poisseuille flow test.
>---------------------------------------------------------------
c1fc7e2ec50dcca1d024e26b9e85e2c4ab40147b
...eat_flux_statistics.h => pressure_statistics.h} | 13 +-
source/postprocess/pressure_statistics.cc | 158 +++++++++++++++++++++
.../{global_refine_amr.prm => poisseuille-2d.prm} | 40 +++---
tests/poisseuille-2d/screen-output | 0
4 files changed, 187 insertions(+), 24 deletions(-)
diff --git a/include/aspect/postprocess/heat_flux_statistics.h b/include/aspect/postprocess/pressure_statistics.h
similarity index 74%
copy from include/aspect/postprocess/heat_flux_statistics.h
copy to include/aspect/postprocess/pressure_statistics.h
index a9c7b98..3c7608a 100644
--- a/include/aspect/postprocess/heat_flux_statistics.h
+++ b/include/aspect/postprocess/pressure_statistics.h
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2011, 2012 by the authors of the ASPECT code.
+ Copyright (C) 2011, 2012, 2014 by the authors of the ASPECT code.
This file is part of ASPECT.
@@ -19,29 +19,28 @@
*/
-#ifndef __aspect__postprocess_heat_flux_statistics_h
-#define __aspect__postprocess_heat_flux_statistics_h
+#ifndef __aspect__postprocess_pressure_statistics_h
+#define __aspect__postprocess_pressure_statistics_h
#include <aspect/postprocess/interface.h>
#include <aspect/simulator_access.h>
-
namespace aspect
{
namespace Postprocess
{
/**
- * A postprocessor that computes some statistics about the heat_flux.
+ * A postprocessor that computes some statistics about the pressure.
*
* @ingroup Postprocessing
*/
template <int dim>
- class HeatFluxStatistics : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
+ class PressureStatistics : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
{
public:
/**
- * Evaluate the solution for some heat_flux statistics.
+ * Evaluate the solution for some pressure statistics.
*/
virtual
std::pair<std::string,std::string>
diff --git a/source/postprocess/pressure_statistics.cc b/source/postprocess/pressure_statistics.cc
new file mode 100644
index 0000000..025b818
--- /dev/null
+++ b/source/postprocess/pressure_statistics.cc
@@ -0,0 +1,158 @@
+/*
+ Copyright (C) 2011 - 2014 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 doc/COPYING. If not see
+ <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <aspect/postprocess/pressure_statistics.h>
+#include <aspect/simulator_access.h>
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <boost/archive/text_oarchive.hpp>
+#include <boost/archive/text_iarchive.hpp>
+
+
+namespace aspect
+{
+ namespace Postprocess
+ {
+ template <int dim>
+ std::pair<std::string,std::string>
+ PressureStatistics<dim>::execute (TableHandler &statistics)
+ {
+ // create a quadrature formula based on the pressure element alone.
+ const QGauss<dim> quadrature_formula (this->get_fe().base_element(this->introspection().base_elements.pressure).degree+1);
+ const unsigned int n_q_points = quadrature_formula.size();
+
+ FEValues<dim> fe_values (this->get_mapping(),
+ this->get_fe(),
+ quadrature_formula,
+ update_values |
+ update_quadrature_points |
+ update_JxW_values);
+
+ std::vector<double> pressure_values(n_q_points);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = this->get_dof_handler().begin_active(),
+ endc = this->get_dof_handler().end();
+
+ double local_pressure_integral = 0;
+
+ // compute the integral quantities by quadrature
+ for (; cell!=endc; ++cell)
+ if (cell->is_locally_owned())
+ {
+ fe_values.reinit (cell);
+ fe_values[this->introspection().extractors.pressure].get_function_values (this->get_solution(),
+ pressure_values);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ local_pressure_integral += pressure_values[q]*fe_values.JxW(q);
+ }
+ }
+
+ // compute min/max by simply
+ // looping over the elements of the
+ // solution vector. the reason is
+ // that minimum and maximum are
+ // usually attained at the
+ // boundary, and so taking their
+ // values at Gauss quadrature
+ // points gives an inaccurate
+ // picture of their true values
+ double local_min_pressure = std::numeric_limits<double>::max();
+ double local_max_pressure = -std::numeric_limits<double>::max();
+ const unsigned int pressure_block = this->introspection().block_indices.pressure;
+ IndexSet range = this->get_solution().block(pressure_block).locally_owned_elements();
+ for (unsigned int i=0; i<range.n_elements(); ++i)
+ {
+ const unsigned int idx = range.nth_index_in_set(i);
+ const double val = this->get_solution().block(pressure_block)(idx);
+
+ local_min_pressure = std::min<double> (local_min_pressure, val);
+ local_max_pressure = std::max<double> (local_max_pressure, val);
+ }
+
+ const double global_pressure_integral
+ = Utilities::MPI::sum (local_pressure_integral, this->get_mpi_communicator());
+ double global_min_pressure = 0;
+ double global_max_pressure = 0;
+
+ // now do the reductions that are
+ // min/max operations. do them in
+ // one communication by multiplying
+ // one value by -1
+ {
+ double local_values[2] = { -local_min_pressure, local_max_pressure };
+ double global_values[2];
+
+ Utilities::MPI::max (local_values, this->get_mpi_communicator(), global_values);
+
+ global_min_pressure = -global_values[0];
+ global_max_pressure = global_values[1];
+ }
+
+ double global_mean_pressure = global_pressure_integral / this->get_volume();
+ statistics.add_value ("Minimal pressure (Pa)",
+ global_min_pressure);
+ statistics.add_value ("Average pressure (Pa)",
+ global_mean_pressure);
+ statistics.add_value ("Maximal pressure (Pa)",
+ global_max_pressure);
+
+ // also make sure that the other columns filled by the this object
+ // all show up with sufficient accuracy and in scientific notation
+ {
+ const char *columns[] = { "Minimal pressure (Pa)",
+ "Average pressure (Pa)",
+ "Maximal pressure (Pa)"
+ };
+ for (unsigned int i=0; i<sizeof(columns)/sizeof(columns[0]); ++i)
+ {
+ statistics.set_precision (columns[i], 8);
+ statistics.set_scientific (columns[i], true);
+ }
+ }
+
+ std::ostringstream output;
+ output.precision(4);
+ output << global_min_pressure << " Pa, "
+ << global_pressure_integral / this->get_volume() << " Pa, "
+ << global_max_pressure << " Pa";
+
+ return std::pair<std::string, std::string> ("Pressure min/avg/max:",
+ output.str());
+ }
+ }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+ namespace Postprocess
+ {
+ ASPECT_REGISTER_POSTPROCESSOR(PressureStatistics,
+ "pressure statistics",
+ "A postprocessor that computes some statistics about "
+ "the pressure field.")
+ }
+}
diff --git a/tests/global_refine_amr.prm b/tests/poisseuille-2d.prm
similarity index 57%
copy from tests/global_refine_amr.prm
copy to tests/poisseuille-2d.prm
index a49a772..269a582 100644
--- a/tests/global_refine_amr.prm
+++ b/tests/poisseuille-2d.prm
@@ -16,6 +16,9 @@ end
subsection Gravity model
set Model name = vertical
+ subsection Vertical
+ set Magnitude = 0
+ end
end
@@ -23,15 +26,14 @@ subsection Geometry model
set Model name = box
subsection Box
- set X extent = 1.2 # default: 1
+ set X extent = 2
set Y extent = 1
- set Z extent = 1
end
end
subsection Initial conditions
- set Model name = perturbed box
+ set Model name = function
end
@@ -39,33 +41,37 @@ subsection Material model
set Model name = simple
subsection Simple model
- set Reference density = 1 # default: 3300
- set Reference specific heat = 1250
- set Reference temperature = 1 # default: 293
- set Thermal conductivity = 1e-6 # default: 4.7
- set Thermal expansion coefficient = 2e-5
- set Viscosity = 1 # default: 5e24
+ set Reference density = 1
+ set Reference specific heat = 1
+ set Reference temperature = 1
+ set Thermal conductivity = 0
+ set Thermal expansion coefficient = 0
+ set Viscosity = 1
end
end
subsection Mesh refinement
set Initial adaptive refinement = 0
- set Initial global refinement = 6
- set Strategy = AMRLeft
-
-
+ set Initial global refinement = 2
end
subsection Model settings
set Fixed temperature boundary indicators = 0, 1
- set Prescribed velocity boundary indicators =
- set Tangential velocity boundary indicators = 1
- set Zero velocity boundary indicators = 0, 2, 3
+ set Prescribed velocity boundary indicators = 0: function, 1: function
+ set Tangential velocity boundary indicators =
+ set Zero velocity boundary indicators = 2, 3
+end
+
+subsection Boundary velocity model
+ subsection Function
+ set Variable names = x,z
+ set Function expression = z*(1-z);0
+ end
end
subsection Postprocess
- set List of postprocessors = visualization, velocity statistics, temperature statistics, heat flux statistics
+ set List of postprocessors = visualization, velocity statistics, pressure statistics
end
diff --git a/tests/poisseuille-2d/screen-output b/tests/poisseuille-2d/screen-output
new file mode 100644
index 0000000..e69de29
More information about the CIG-COMMITS
mailing list