[cig-commits] commit 1997 by heister to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Mon Oct 28 12:59:27 PDT 2013
Revision 1997
add visualization plugin to output entropy viscosity
U trunk/aspect/doc/modules/changes.h
A trunk/aspect/include/aspect/postprocess/visualization/artificial_viscosity.h
U trunk/aspect/include/aspect/simulator.h
U trunk/aspect/include/aspect/simulator_access.h
A trunk/aspect/source/postprocess/visualization/artificial_viscosity.cc
U trunk/aspect/source/simulator/assembly.cc
U trunk/aspect/source/simulator/simulator_access.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=1997&peg=1997
Diff:
Modified: trunk/aspect/doc/modules/changes.h
===================================================================
--- trunk/aspect/doc/modules/changes.h 2013-10-25 21:17:39 UTC (rev 1996)
+++ trunk/aspect/doc/modules/changes.h 2013-10-28 19:58:43 UTC (rev 1997)
@@ -8,6 +8,11 @@
</p>
<ol>
+ <li>New: add a visual postprocessor that outputs the artificial
+ viscosity parameter for the temperature equation on each cell.
+ <br>
+ (Timo Heister 2013/10/28)
+
<li>Fixed: moved particle generation to a class, changed particle
integration and generation to be factory patterned classes. There
should be no effect on the user but this will allow for easier
Added: trunk/aspect/include/aspect/postprocess/visualization/artificial_viscosity.h
===================================================================
--- trunk/aspect/include/aspect/postprocess/visualization/artificial_viscosity.h (rev 0)
+++ trunk/aspect/include/aspect/postprocess/visualization/artificial_viscosity.h 2013-10-28 19:58:43 UTC (rev 1997)
@@ -0,0 +1,70 @@
+/*
+ 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: error_indicator.h 1491 2012-12-11 06:05:16Z bangerth $ */
+
+
+#ifndef __aspect__postprocess_visualization_artificial_viscosity_h
+#define __aspect__postprocess_visualization_artificial_viscosity_h
+
+#include <aspect/postprocess/visualization.h>
+#include <aspect/simulator_access.h>
+
+
+namespace aspect
+{
+ namespace Postprocess
+ {
+ namespace VisualizationPostprocessors
+ {
+ /**
+ * A class derived that implements a function that provides
+ * the artificial viscosity for the temperature equation
+ * on each cell for graphical output.
+ */
+ template <int dim>
+ class ArtificialViscosity : public CellDataVectorCreator<dim>,
+ public SimulatorAccess<dim>
+ {
+ public:
+ /**
+ * The function classes have to implement that want to output
+ * cellwise data.
+ * @return A pair of values with the following meaning:
+ * - The first element provides the name by which this
+ * data should be written to the output file.
+ * - The second element is a pointer to a vector with
+ * one element per active cell on the current processor.
+ * Elements corresponding to active cells that are either
+ * artificial or ghost cells (in deal.II language, see the
+ * deal.II glossary) will be ignored but must nevertheless
+ * exist in the returned vector. While implementations of this
+ * function must create this vector, ownership is taken over
+ * by the caller of this function and the caller will take
+ * care of destroying the vector pointed to.
+ */
+ virtual
+ std::pair<std::string, Vector<float> *>
+ execute () const;
+ };
+ }
+ }
+}
+
+#endif
Modified: trunk/aspect/include/aspect/simulator.h
===================================================================
--- trunk/aspect/include/aspect/simulator.h 2013-10-25 21:17:39 UTC (rev 1996)
+++ trunk/aspect/include/aspect/simulator.h 2013-10-28 19:58:43 UTC (rev 1997)
@@ -736,6 +736,11 @@
void make_pressure_rhs_compatible(LinearAlgebra::BlockVector &vector);
/**
+ * Fills a vector with the artificial viscosity for the temperature on each local cell
+ */
+ void get_artificial_viscosity (Vector<float> &viscosity_per_cell) const;
+
+ /**
* Internal routine to compute the depth average of a certain quantitiy.
* The functor @p fctr should implement:
* 1. bool need_material_properties()
Modified: trunk/aspect/include/aspect/simulator_access.h
===================================================================
--- trunk/aspect/include/aspect/simulator_access.h 2013-10-25 21:17:39 UTC (rev 1996)
+++ trunk/aspect/include/aspect/simulator_access.h 2013-10-28 19:58:43 UTC (rev 1997)
@@ -202,6 +202,14 @@
*/
void
get_refinement_criteria(Vector<float> &estimated_error_per_cell) const;
+
+ /**
+ * Returns the entropy viscosity on each locally owned cell as it is used
+ * to stabilize the temperature equation.
+ */
+ void
+ get_artificial_viscosity(Vector<float> &viscosity_per_cell) const;
+
/** @} */
Added: trunk/aspect/source/postprocess/visualization/artificial_viscosity.cc
===================================================================
--- trunk/aspect/source/postprocess/visualization/artificial_viscosity.cc (rev 0)
+++ trunk/aspect/source/postprocess/visualization/artificial_viscosity.cc 2013-10-28 19:58:43 UTC (rev 1997)
@@ -0,0 +1,62 @@
+/*
+ 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: error_indicator.cc 1491 2012-12-11 06:05:16Z bangerth $ */
+
+
+#include <aspect/postprocess/visualization/artificial_viscosity.h>
+#include <aspect/simulator_access.h>
+
+
+namespace aspect
+{
+ namespace Postprocess
+ {
+ namespace VisualizationPostprocessors
+ {
+ template <int dim>
+ std::pair<std::string, Vector<float> *>
+ ArtificialViscosity<dim>::execute() const
+ {
+ std::pair<std::string, Vector<float> *>
+ return_value ("artificial_viscosity",
+ new Vector<float>(this->get_triangulation().n_active_cells()));
+ this->get_artificial_viscosity(*return_value.second);
+ return return_value;
+ }
+ }
+ }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+ namespace Postprocess
+ {
+ namespace VisualizationPostprocessors
+ {
+ ASPECT_REGISTER_VISUALIZATION_POSTPROCESSOR(ArtificialViscosity,
+ "artificial viscosity",
+ "A visualization output object that generates output "
+ "showing the value of the artificial viscosity on each "
+ "cell.")
+ }
+ }
+}
Modified: trunk/aspect/source/simulator/assembly.cc
===================================================================
--- trunk/aspect/source/simulator/assembly.cc 2013-10-25 21:17:39 UTC (rev 1996)
+++ trunk/aspect/source/simulator/assembly.cc 2013-10-28 19:58:43 UTC (rev 1997)
@@ -681,7 +681,173 @@
}
}
+ template <int dim>
+ void
+ Simulator<dim>::
+ get_artificial_viscosity (Vector<float> &viscosity_per_cell) const
+ {
+ viscosity_per_cell = 0.0;
+ //Assert length
+ Assert(viscosity_per_cell.size()==triangulation.n_active_cells(), ExcInternalError());
+
+ TemperatureOrComposition torc(TemperatureOrComposition::temperature_field);
+ const std::pair<double,double>
+ global_field_range = get_extrapolated_temperature_or_composition_range (torc);
+ double global_entropy_variation = get_entropy_variation ((global_field_range.first +
+ global_field_range.second) / 2, torc);
+ double global_max_velocity = get_maximal_velocity(old_solution);
+
+
+ internal::Assembly::Scratch::
+ AdvectionSystem<dim> scratch (finite_element,
+ finite_element.base_element(torc.is_temperature()
+ ?
+ introspection.block_indices.temperature
+ :
+ introspection.block_indices.compositional_fields[0]),
+ mapping,
+ QGauss<dim>((torc.is_temperature()
+ ?
+ parameters.temperature_degree
+ :
+ parameters.composition_degree)
+ +
+ (parameters.stokes_velocity_degree+1)/2),
+ parameters.n_compositional_fields);
+
+ typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active();
+ for (unsigned int cellidx=0;cellidx<triangulation.n_active_cells();++cellidx, ++cell)
+ {
+ if (!cell->is_locally_owned())
+ {
+ viscosity_per_cell[cellidx]=-1;
+ continue;
+ }
+
+ const bool use_bdf2_scheme = (timestep_number > 1);
+
+ const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points;
+
+ // also have the number of dofs that correspond just to the element for
+ // the system we are currently trying to assemble
+ const unsigned int advection_dofs_per_cell = scratch.phi_field.size();
+
+ Assert (advection_dofs_per_cell < scratch.finite_element_values.get_fe().dofs_per_cell, ExcInternalError());
+ Assert (scratch.grad_phi_field.size() == advection_dofs_per_cell, ExcInternalError());
+ Assert (scratch.phi_field.size() == advection_dofs_per_cell, ExcInternalError());
+
+ const unsigned int solution_component
+ = (torc.is_temperature()
+ ?
+ introspection.component_indices.temperature
+ :
+ introspection.component_indices.compositional_fields[torc.compositional_variable]
+ );
+
+ const FEValuesExtractors::Scalar solution_field
+ = (torc.is_temperature()
+ ?
+ introspection.extractors.temperature
+ :
+ introspection.extractors.compositional_fields[torc.compositional_variable]
+ );
+
+ scratch.finite_element_values.reinit (cell);
+
+ // get all dof indices on the current cell, then extract those
+ // that correspond to the solution_field we are interested in
+ cell->get_dof_indices (scratch.local_dof_indices);
+
+ if (torc.is_temperature())
+ {
+ scratch.finite_element_values[introspection.extractors.temperature].get_function_values (old_solution,
+ scratch.old_temperature_values);
+ scratch.finite_element_values[introspection.extractors.temperature].get_function_values (old_old_solution,
+ scratch.old_old_temperature_values);
+
+ scratch.finite_element_values[introspection.extractors.velocities].get_function_symmetric_gradients (old_solution,
+ scratch.old_strain_rates);
+ scratch.finite_element_values[introspection.extractors.velocities].get_function_symmetric_gradients (old_old_solution,
+ scratch.old_old_strain_rates);
+
+ scratch.finite_element_values[introspection.extractors.pressure].get_function_values (old_solution,
+ scratch.old_pressure);
+ scratch.finite_element_values[introspection.extractors.pressure].get_function_values (old_old_solution,
+ scratch.old_old_pressure);
+
+ for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
+ {
+ scratch.finite_element_values[introspection.extractors.compositional_fields[c]].get_function_values(old_solution,
+ scratch.old_composition_values[c]);
+ scratch.finite_element_values[introspection.extractors.compositional_fields[c]].get_function_values(old_old_solution,
+ scratch.old_old_composition_values[c]);
+ }
+ }
+ else
+ {
+ scratch.finite_element_values[introspection.extractors.compositional_fields[torc.compositional_variable]].get_function_values(old_solution,
+ scratch.old_composition_values[torc.compositional_variable]);
+ scratch.finite_element_values[introspection.extractors.compositional_fields[torc.compositional_variable]].get_function_values(old_old_solution,
+ scratch.old_old_composition_values[torc.compositional_variable]);
+ }
+
+ scratch.finite_element_values[introspection.extractors.velocities].get_function_values (old_solution,
+ scratch.old_velocity_values);
+ scratch.finite_element_values[introspection.extractors.velocities].get_function_values (old_old_solution,
+ scratch.old_old_velocity_values);
+ scratch.finite_element_values[introspection.extractors.velocities].get_function_values(current_linearization_point,
+ scratch.current_velocity_values);
+
+
+ scratch.old_field_values = ((torc.is_temperature()) ? &scratch.old_temperature_values : &scratch.old_composition_values[torc.compositional_variable]);
+ scratch.old_old_field_values = ((torc.is_temperature()) ? &scratch.old_old_temperature_values : &scratch.old_old_composition_values[torc.compositional_variable]);
+
+ scratch.finite_element_values[solution_field].get_function_gradients (old_solution,
+ scratch.old_field_grads);
+ scratch.finite_element_values[solution_field].get_function_gradients (old_old_solution,
+ scratch.old_old_field_grads);
+
+ scratch.finite_element_values[solution_field].get_function_laplacians (old_solution,
+ scratch.old_field_laplacians);
+ scratch.finite_element_values[solution_field].get_function_laplacians (old_old_solution,
+ scratch.old_old_field_laplacians);
+
+ if (torc.is_temperature())
+ {
+ compute_material_model_input_values (current_linearization_point,
+ scratch.finite_element_values,
+ true,
+ scratch.material_model_inputs);
+ material_model->evaluate(scratch.material_model_inputs,scratch.material_model_outputs);
+
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ scratch.explicit_material_model_inputs.temperature[q] = (scratch.old_temperature_values[q] + scratch.old_old_temperature_values[q]) / 2;
+ scratch.explicit_material_model_inputs.position[q] = scratch.finite_element_values.quadrature_point(q);
+ scratch.explicit_material_model_inputs.pressure[q] = (scratch.old_pressure[q] + scratch.old_old_pressure[q]) / 2;
+ for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
+ scratch.explicit_material_model_inputs.composition[q][c] = (scratch.old_composition_values[c][q] + scratch.old_old_composition_values[c][q]) / 2;
+ scratch.explicit_material_model_inputs.strain_rate[q] = (scratch.old_strain_rates[q] + scratch.old_old_strain_rates[q]) / 2;
+ }
+ material_model->evaluate(scratch.explicit_material_model_inputs,scratch.explicit_material_model_outputs);
+ }
+
+
+ viscosity_per_cell[cellidx] = compute_viscosity(scratch,
+ global_max_velocity,
+ global_field_range.second - global_field_range.first,
+ 0.5 * (global_field_range.second + global_field_range.first),
+ global_entropy_variation,
+ cell->diameter(),
+ torc);
+ }
+
+
+ }
+
+
+
template <int dim>
void
Simulator<dim>::
@@ -1523,6 +1689,7 @@
template void Simulator<dim>::copy_local_to_global_stokes_system ( \
const internal::Assembly::CopyData::StokesSystem<dim> &data); \
template void Simulator<dim>::assemble_stokes_system (); \
+ template void Simulator<dim>::get_artificial_viscosity (Vector<float> &viscosity_per_cell) const; \
template void Simulator<dim>::build_advection_preconditioner (const TemperatureOrComposition &, \
std_cxx1x::shared_ptr<aspect::LinearAlgebra::PreconditionILU> &preconditioner); \
template void Simulator<dim>::local_assemble_advection_system ( \
Modified: trunk/aspect/source/simulator/simulator_access.cc
===================================================================
--- trunk/aspect/source/simulator/simulator_access.cc 2013-10-25 21:17:39 UTC (rev 1996)
+++ trunk/aspect/source/simulator/simulator_access.cc 2013-10-28 19:58:43 UTC (rev 1997)
@@ -158,8 +158,13 @@
simulator->mesh_refinement_manager.execute (estimated_error_per_cell);
}
+ template <int dim>
+ void
+ SimulatorAccess<dim>::get_artificial_viscosity (Vector<float> &viscosity_per_cell) const
+ {
+ simulator->get_artificial_viscosity(viscosity_per_cell);
+ }
-
template <int dim>
const LinearAlgebra::BlockVector &
SimulatorAccess<dim>::get_solution () const
More information about the CIG-COMMITS
mailing list