[cig-commits] [commit] master: Initial version of spherical_velocity_statistics postprocessor (1d0cb79)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue May 20 11:26:35 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/19fcaae1dae23206c4e6f1f4e9965965f0128967...6204d1b7337a0741858189f8e4ee9cf2f447ebd6
>---------------------------------------------------------------
commit 1d0cb79359042ee835f22f018e71a0f8b12b2dec
Author: anne-glerum <a.c.glerum at uu.nl>
Date: Mon May 19 23:11:58 2014 +0200
Initial version of spherical_velocity_statistics postprocessor
>---------------------------------------------------------------
1d0cb79359042ee835f22f018e71a0f8b12b2dec
...tatistics.h => spherical_velocity_statistics.h} | 12 +-
.../postprocess/spherical_velocity_statistics.cc | 182 +++++++++++++++++++++
2 files changed, 188 insertions(+), 6 deletions(-)
diff --git a/include/aspect/postprocess/velocity_statistics.h b/include/aspect/postprocess/spherical_velocity_statistics.h
similarity index 77%
copy from include/aspect/postprocess/velocity_statistics.h
copy to include/aspect/postprocess/spherical_velocity_statistics.h
index 1792b47..30d4213 100644
--- a/include/aspect/postprocess/velocity_statistics.h
+++ b/include/aspect/postprocess/spherical_velocity_statistics.h
@@ -17,11 +17,11 @@
along with ASPECT; see the file doc/COPYING. If not see
<http://www.gnu.org/licenses/>.
*/
-/* $Id$ */
+/* $Id: spherical_velocity_statistics.h 1433 2012-12-08 08:24:55Z bangerth $ */
-#ifndef __aspect__postprocess_velocity_statistics_h
-#define __aspect__postprocess_velocity_statistics_h
+#ifndef __aspect__postprocess_spherical_velocity_statistics_h
+#define __aspect__postprocess_spherical_velocity_statistics_h
#include <aspect/postprocess/interface.h>
#include <aspect/simulator_access.h>
@@ -32,17 +32,17 @@ namespace aspect
{
/**
- * A postprocessor that computes some statistics about the velocity.
+ * A postprocessor that computes some statistics about the velocity of spherical/cylindrical models.
*
* @ingroup Postprocessing
*/
template <int dim>
- class VelocityStatistics : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
+ class SphericalVelocityStatistics : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
{
public:
/**
* Evaluate the solution for some velocity statistics.
- */
+ **/
virtual
std::pair<std::string,std::string>
execute (TableHandler &statistics);
diff --git a/source/postprocess/spherical_velocity_statistics.cc b/source/postprocess/spherical_velocity_statistics.cc
new file mode 100644
index 0000000..5a61bef
--- /dev/null
+++ b/source/postprocess/spherical_velocity_statistics.cc
@@ -0,0 +1,182 @@
+/*
+ Copyright (C) 2011, 2012 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/>.
+*/
+/* $Id: spherical_velocity_statistics.cc 1969 2013-10-16 19:32:59Z heien $ */
+
+
+#include <aspect/postprocess/spherical_velocity_statistics.h>
+#include <aspect/geometry_model/spherical_shell.h>
+#include <aspect/simulator_access.h>
+#include <aspect/global.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>
+ SphericalVelocityStatistics<dim>::execute (TableHandler &statistics)
+ {
+ Assert (dynamic_cast<const GeometryModel::SphericalShell<dim>*> (&this->get_geometry_model()) != 0,
+ ExcMessage ("This postprocessor can only be used if the geometry "
+ "is a spherical shell."));
+
+ const QGauss<dim> quadrature_formula (this->get_fe()
+ .base_element(this->introspection().base_elements.velocities).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<Tensor<1,dim> > velocity_values(n_q_points);
+ std::vector<Point<dim> > position_point(n_q_points);
+
+ double local_rad_velocity_square_integral = 0;
+ double local_tan_velocity_square_integral = 0;
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = this->get_dof_handler().begin_active(),
+ endc = this->get_dof_handler().end();
+ for (; cell!=endc; ++cell)
+ if (cell->is_locally_owned())
+ {
+ fe_values.reinit (cell);
+ fe_values[this->introspection().extractors.velocities].get_function_values (this->get_solution(),
+ velocity_values);
+ position_point = fe_values.get_quadrature_points();
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ {
+ // create unit vector in radial direction
+ Tensor<1,dim> position;
+// for (unsigned int i = 0; i<dim; ++i)
+// position[i] = position_point[q][i];
+ position = position_point[q];
+ position /= position.norm();
+ // compute the radial velocity by multiplying with the unit vector in the radial direction twice
+ Tensor<1,dim> vel = (velocity_values[q] * position) * position;
+ local_rad_velocity_square_integral += (vel * vel) *
+ fe_values.JxW(q);
+ // compute the tangential velocity by subtracting the radial velocity from the velocity
+ vel /= -1;
+ vel += velocity_values[q];
+ local_tan_velocity_square_integral += (vel * vel) *
+ fe_values.JxW(q);
+ }
+ }
+
+ const double global_rad_velocity_square_integral
+ = Utilities::MPI::sum (local_rad_velocity_square_integral, this->get_mpi_communicator());
+ const double global_tan_velocity_square_integral
+ = Utilities::MPI::sum (local_tan_velocity_square_integral, this->get_mpi_communicator());
+
+ const double rad_vrms = std::sqrt(global_rad_velocity_square_integral) /
+ std::sqrt(this->get_volume());
+ const double tan_vrms = std::sqrt(global_tan_velocity_square_integral) /
+ std::sqrt(this->get_volume());
+ const double vrms = std::sqrt(rad_vrms*rad_vrms + tan_vrms*tan_vrms);
+
+ if (this->convert_output_to_years() == true)
+ {
+ statistics.add_value ("Radial RMS velocity (m/yr)",
+ rad_vrms * year_in_seconds);
+ statistics.add_value ("Tangential RMS velocity (m/yr)",
+ tan_vrms * year_in_seconds);
+ statistics.add_value ("Total RMS velocity (m/yr)",
+ vrms * year_in_seconds);
+
+ // 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[] = { "Radial RMS velocity (m/yr)",
+ "Tangential RMS velocity (m/yr)",
+ "Total RMS velocity (m/yr)"
+ };
+ for (unsigned int i=0; i<sizeof(columns)/sizeof(columns[0]); ++i)
+ {
+ statistics.set_precision (columns[i], 8);
+ statistics.set_scientific (columns[i], true);
+ }
+ }
+ }
+ else
+ {
+ statistics.add_value ("Radial RMS velocity (m/s)", rad_vrms);
+ statistics.add_value ("Tangential RMS velocity (m/s)", tan_vrms);
+ statistics.add_value ("Total RMS velocity (m/s)", vrms);
+
+ // 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[] = { "Radial RMS velocity (m/s)",
+ "Tangential RMS velocity (m/s)",
+ "Total RMS velocity (m/s)"
+ };
+ 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(3);
+ if (this->convert_output_to_years() == true)
+ output << rad_vrms *year_in_seconds
+ << " m/yr, "
+ << tan_vrms *year_in_seconds
+ << " m/yr, "
+ << vrms *year_in_seconds
+ << " m/yr";
+ else
+ output << rad_vrms
+ << " m/s, "
+ << tan_vrms
+ << " m/s, "
+ << vrms
+ << " m/s";
+
+ return std::pair<std::string, std::string> ("Radial RMS, tangential RMS, total RMS velocity:",
+ output.str());
+ }
+ }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+ namespace Postprocess
+ {
+ ASPECT_REGISTER_POSTPROCESSOR(SphericalVelocityStatistics,
+ "spherical velocity statistics",
+ "A postprocessor that computes radial, tangential and total RMS "
+ "velocity.")
+ }
+}
More information about the CIG-COMMITS
mailing list