[cig-commits] commit 1882 by dannberg to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Wed Sep 11 08:02:28 PDT 2013
Revision 1882
add the mesh refinement criterion viscosity
A branches/j-dannberg/include/aspect/mesh_refinement/viscosity.h
A branches/j-dannberg/source/mesh_refinement/viscosity.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=1882&peg=1882
Diff:
Added: branches/j-dannberg/include/aspect/mesh_refinement/viscosity.h
===================================================================
--- branches/j-dannberg/include/aspect/mesh_refinement/viscosity.h (rev 0)
+++ branches/j-dannberg/include/aspect/mesh_refinement/viscosity.h 2013-09-11 15:00:37 UTC (rev 1882)
@@ -0,0 +1,63 @@
+/*
+ Copyright (C) 2013 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: density.h 1433 2012-12-08 08:24:55Z bangerth $ */
+
+
+#ifndef __aspect__mesh_refinement_viscosity_h
+#define __aspect__mesh_refinement_viscosity_h
+
+#include <aspect/mesh_refinement/interface.h>
+#include <aspect/simulator_access.h>
+
+namespace aspect
+{
+ namespace MeshRefinement
+ {
+
+ /**
+ * A class that implements a mesh refinement criterion based on
+ * a field representing the viscosity $\eta$ as a
+ * spatially variable function and that we compute from the solution
+ * vector (assuming the viscosity actually depends on the solution).
+ *
+ * @ingroup MeshRefinement
+ */
+ template <int dim>
+ class Viscosity : public Interface<dim>,
+ public SimulatorAccess<dim>
+ {
+ public:
+ /**
+ * Execute this mesh refinement criterion.
+ *
+ * @param[out] error_indicators A vector that for every active
+ * cell of the current mesh
+ * (which may be a partition of a distributed mesh) provides an error
+ * indicator. This vector will already have the correct size when the
+ * function is called.
+ */
+ virtual
+ void
+ execute (Vector<float> &error_indicators) const;
+ };
+ }
+}
+
+#endif
Added: branches/j-dannberg/source/mesh_refinement/viscosity.cc
===================================================================
--- branches/j-dannberg/source/mesh_refinement/viscosity.cc (rev 0)
+++ branches/j-dannberg/source/mesh_refinement/viscosity.cc 2013-09-11 15:00:37 UTC (rev 1882)
@@ -0,0 +1,168 @@
+/*
+ Copyright (C) 2013 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: density.cc 1624 2013-04-28 21:15:28Z heister $ */
+
+
+#include <aspect/mesh_refinement/viscosity.h>
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/derivative_approximation.h>
+
+namespace aspect
+{
+ namespace MeshRefinement
+ {
+ template <int dim>
+ void
+ Viscosity<dim>::execute(Vector<float> &indicators) const
+ {
+ indicators = 0;
+
+//TODO: if the viscosity doesn't actually depend on the solution
+ // then we can get away with simply interpolating it spatially
+
+
+ // create a vector in which we set the temperature block to
+ // be a finite element interpolation of the viscosity.
+ // we do so by setting up a quadrature formula with the
+ // temperature unit support points, then looping over these
+ // points, compute the output quantity at them, and writing
+ // the result into the output vector in the same order
+ // (because quadrature points and temperature dofs are,
+ // by design of the quadrature formula, numbered in the
+ // same way)
+ LinearAlgebra::BlockVector vec_distributed (this->introspection().index_sets.system_partitioning,
+ this->get_mpi_communicator());
+
+ const Quadrature<dim> quadrature(this->get_fe().base_element(2).get_unit_support_points());
+ std::vector<unsigned int> local_dof_indices (this->get_fe().dofs_per_cell);
+ FEValues<dim> fe_values (this->get_mapping(),
+ this->get_fe(),
+ quadrature,
+ update_quadrature_points | update_values);
+
+ // the values of the compositional fields are stored as blockvectors for each field
+ // we have to extract them in this structure
+ std::vector<std::vector<double> > prelim_composition_values (this->n_compositional_fields(),
+ std::vector<double> (quadrature.size()));
+
+ typename MaterialModel::Interface<dim>::MaterialModelInputs in(quadrature.size(),
+ this->n_compositional_fields());
+ typename MaterialModel::Interface<dim>::MaterialModelOutputs out(quadrature.size(),
+ this->n_compositional_fields());
+
+ 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.pressure].get_function_values (this->get_solution(),
+ in.pressure);
+ fe_values[this->introspection().extractors.temperature].get_function_values (this->get_solution(),
+ in.temperature);
+ 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(),
+ prelim_composition_values[c]);
+
+ in.position = fe_values.get_quadrature_points();
+ in.strain_rate.resize(0);// we are not reading the viscosity
+ for (unsigned int i=0;i<quadrature.size();++i)
+ {
+ for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
+ in.composition[i][c] = prelim_composition_values[c][i];
+ }
+ this->get_material_model().evaluate(in, out);
+
+ cell->get_dof_indices (local_dof_indices);
+
+ // for each temperature dof, write into the output
+ // vector the viscosity. note that quadrature points and
+ // dofs are enumerated in the same order
+ for (unsigned int i=0; i<this->get_fe().base_element(2).dofs_per_cell; ++i)
+ {
+ const unsigned int system_local_dof
+ = this->get_fe().component_to_system_index(/*temperature component=*/dim+1,
+ /*dof index within component=*/i);
+
+ vec_distributed(local_dof_indices[system_local_dof])
+ = std::log(out.viscosities[i]);
+ }
+ }
+
+ // now create a vector with the requisite ghost elements
+ // and use it for estimating the gradients
+ LinearAlgebra::BlockVector vec (this->introspection().index_sets.system_relevant_partitioning,
+ this->get_mpi_communicator());
+ vec = vec_distributed;
+
+ DerivativeApproximation::approximate_gradient (this->get_mapping(),
+ this->get_dof_handler(),
+ vec,
+ indicators,
+//TODO: replace by the appropriate component mask
+ dim+1);
+
+ // Scale gradient in each cell with the correct power of h. Otherwise,
+ // error indicators do not reduce when refined if there is a viscosity
+ // jump. We need at least order 1 for the error not to grow when
+ // refining, so anything >1 should work. (note that the gradient
+ // itself scales like 1/h, so multiplying it with any factor h^s, s>1
+ // will yield convergence of the error indicators to zero as h->0)
+ const double power = 1.0 + dim/2.0;
+ {
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = this->get_dof_handler().begin_active(),
+ endc = this->get_dof_handler().end();
+ unsigned int i=0;
+ for (; cell!=endc; ++cell, ++i)
+ if (cell->is_locally_owned())
+ indicators(i) *= std::pow(cell->diameter(), power);
+ }
+ }
+ }
+}
+
+// explicit instantiations
+namespace aspect
+{
+ namespace MeshRefinement
+ {
+ ASPECT_REGISTER_MESH_REFINEMENT_CRITERION(Viscosity,
+ "viscosity",
+ "A mesh refinement criterion that computes "
+ "refinement indicators from a field that describes "
+ "the spatial variability of the viscosity, $\eta$. "
+ "Because this quantity may not be a continuous function ($\eta$ "
+ "may be a discontinuous function along discontinuities in the "
+ "medium, for example due to phase changes), we approximate the "
+ "gradient of this quantity to refine the mesh. The error indicator "
+ "defined here takes the magnitude of the approximate gradient "
+ "and scales it by $h_K^{1+d/2}$ where $h_K$ is the diameter of each cell "
+ "and $d$ is the dimension. "
+ "This scaling ensures that the error indicators converge to zero as "
+ "$h_K\rightarrow 0$ even if the energy density is discontinuous, since "
+ "the gradient of a discontinuous function grows like $1/h_K$.")
+ }
+}
More information about the CIG-COMMITS
mailing list