[cig-commits] commit 2340 by bangerth to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Thu Mar 20 20:33:38 PDT 2014
Revision 2340
Add a postprocessor that computes dynamic topography. This is inspired by Jacky Austermann's needs.
U trunk/aspect/doc/modules/changes.h
A trunk/aspect/include/aspect/postprocess/dynamic_topography.h
A trunk/aspect/source/postprocess/dynamic_topography.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=2340&peg=2340
Diff:
Modified: trunk/aspect/doc/modules/changes.h
===================================================================
--- trunk/aspect/doc/modules/changes.h 2014-03-19 15:53:16 UTC (rev 2339)
+++ trunk/aspect/doc/modules/changes.h 2014-03-21 03:33:34 UTC (rev 2340)
@@ -8,6 +8,11 @@
</p>
<ol>
+ <li>New: There is now a postprocess that computes a measure of
+ dynamic topography, based on the vertical component of the stress.
+ <br>
+ (Wolfgang Bangerth, Jacqueline Austermann 2014/03/03)
+
<li>New: Aspect now installs a file <code>AspectConfig.cmake</code>
into the same directory as the executable that can be used by
plugins to set up compiler flags, include paths, etc, to compile
Copied: trunk/aspect/include/aspect/postprocess/dynamic_topography.h (from rev 2331, trunk/aspect/include/aspect/postprocess/heat_flux_statistics.h)
===================================================================
--- trunk/aspect/include/aspect/postprocess/dynamic_topography.h (rev 0)
+++ trunk/aspect/include/aspect/postprocess/dynamic_topography.h 2014-03-21 03:33:34 UTC (rev 2340)
@@ -0,0 +1,55 @@
+/*
+ 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$ */
+
+
+#ifndef __aspect__postprocess_surface_topography_h
+#define __aspect__postprocess_surface_topography_h
+
+#include <aspect/postprocess/interface.h>
+#include <aspect/simulator_access.h>
+
+
+namespace aspect
+{
+ namespace Postprocess
+ {
+
+ /**
+ * A postprocessor that computes dynamic topography at the surface.
+ *
+ * @ingroup Postprocessing
+ */
+ template <int dim>
+ class DynamicTopography : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
+ {
+ public:
+ /**
+ * Evaluate the solution for the dynamic topography.
+ **/
+ virtual
+ std::pair<std::string,std::string>
+ execute (TableHandler &statistics);
+ };
+ }
+}
+
+
+#endif
Copied: trunk/aspect/source/postprocess/dynamic_topography.cc (from rev 2331, trunk/aspect/source/postprocess/heat_flux_statistics.cc)
===================================================================
--- trunk/aspect/source/postprocess/dynamic_topography.cc (rev 0)
+++ trunk/aspect/source/postprocess/dynamic_topography.cc 2014-03-21 03:33:34 UTC (rev 2340)
@@ -0,0 +1,153 @@
+/*
+ Copyright (C) 2011, 2012, 2013, 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/>.
+*/
+/* $Id$ */
+
+
+#include <aspect/postprocess/dynamic_topography.h>
+#include <aspect/simulator_access.h>
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/fe/fe_values.h>
+
+
+namespace aspect
+{
+ namespace Postprocess
+ {
+ template <int dim>
+ std::pair<std::string,std::string>
+ DynamicTopography<dim>::execute (TableHandler &statistics)
+ {
+ const QMidpoint<dim> quadrature_formula;
+
+ FEValues<dim> fe_values (this->get_mapping(),
+ this->get_fe(),
+ quadrature_formula,
+ update_values |
+ update_gradients |
+ update_q_points);
+
+ typename MaterialModel::Interface<dim>::MaterialModelInputs in(fe_values.n_quadrature_points, this->n_compositional_fields());
+ typename MaterialModel::Interface<dim>::MaterialModelOutputs out(fe_values.n_quadrature_points, this->n_compositional_fields());
+
+ std::vector<std::vector<double> > composition_values (this->n_compositional_fields(),std::vector<double> (quadrature_formula.size()));
+
+ std::string filename = this->get_output_directory() +
+ "surface_topography." +
+ Utilities::int_to_string
+ (this->get_triangulation().locally_owned_subdomain(), 4);
+ std::ofstream output (filename.c_str());
+
+ // loop over all of the surface cells and if one less than h/2 away from
+ // the top surface, evaluate the stress at its center
+ 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())
+ if (cell->at_boundary() &&
+ (this->get_geometry_model().depth (cell->center()) < cell->diameter()/2))
+ {
+ fe_values.reinit (cell);
+
+ // get the various components of the solution, then
+ // evaluate the material properties there
+ fe_values[this->introspection().extractors.temperature]
+ .get_function_values (this->get_solution(), in.temperature);
+ fe_values[this->introspection().extractors.pressure]
+ .get_function_values (this->get_solution(), in.pressure);
+ fe_values[this->introspection().extractors.velocities]
+ .get_function_symmetric_gradients (this->get_solution(), in.strain_rate);
+
+ in.position = fe_values.get_quadrature_points();
+
+ for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
+ fe_values[this->introspection().extractors.compositional_fields[c]]
+ .get_function_values(this->get_solution(),
+ composition_values[c]);
+ for (unsigned int i=0; i<fe_values.n_quadrature_points; ++i)
+ {
+ for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
+ in.composition[i][c] = composition_values[c][i];
+ }
+
+ this->get_material_model().evaluate(in, out);
+
+
+ // for each of the quadrature points, evaluate the
+ // stress and compute the component in direction of the
+ // gravity vector
+ for (unsigned int q=0; q<quadrature_formula.size(); ++q)
+ {
+ const Point<dim> location = fe_values.quadrature_point(q);
+ const double viscosity = out.viscosities[q];
+ const double density = out.densities[q];
+
+ const SymmetricTensor<2,dim> stress = viscosity * in.strain_rate[q];
+
+ const Tensor<1,dim> gravity = this->get_gravity_model().gravity_vector(location);
+ const Tensor<1,dim> gravity_direction = gravity/gravity.norm();
+
+ const double sigma_rr = gravity_direction * (stress * gravity_direction);
+ const double dynamic_topography = sigma_rr / gravity.norm() / density;
+
+ output << location
+ << ' '
+ << dynamic_topography
+ << std::endl;
+ }
+ }
+
+ return std::pair<std::string,std::string>("Writing surface topography...",
+ filename);
+ }
+ }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+ namespace Postprocess
+ {
+ ASPECT_REGISTER_POSTPROCESSOR(DynamicTopography,
+ "dynamic topography",
+ "A postprocessor that computes a measure of dynamic topograph "
+ "based on the stress at the surface."
+ "
"
+ "The exact approach works as follows: At the centers of all cells "
+ "that sit along the top surface, we evaluate the stress and "
+ "evaluate the component of it in the direction in which "
+ "gravity acts. In other words, we compute "
+ "$\sigma_{rr}={\hat g}^T(\eta \varepsilon(\mathbf u))\hat g$ "
+ "where $\hat g = \mathbf g/\|\mathbf g\|$ is the direction of "
+ "the gravity vector $\mathbf g$. From this, the dynamic "
+ "topography is computed using the formula "
+ "$h=\frac{\sigma_{rr}}{\|\mathbf g\| \rho}$ where $\rho$ "
+ "is the density at the cell center."
+ "
"
+ "(As a side note, the postprocessor chooses the cell center "
+ "instead of the center of the cell face at the surface, where we "
+ "really are interested in the quantity, since "
+ "this often gives better accuracy. The results should in essence "
+ "be the same, though.)")
+ }
+}
More information about the CIG-COMMITS
mailing list