[cig-commits] [commit] master: Internal heating statistics computes mass averages. (f627ac0)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Sun Jun 22 04:22:27 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/8ad680900561ab32ad4d0420b88f2843ad3d74b5...48f5db061ecf9e78a0c7ef83891ac7743617afcd
>---------------------------------------------------------------
commit f627ac0812e6d6f2ed8409102917d78fa7b79485
Author: Rene Gassmoeller <R.Gassmoeller at mailbox.org>
Date: Fri Jun 20 15:58:05 2014 +0200
Internal heating statistics computes mass averages.
Internal heating statistics should compute mass averages and integrate
over mass instead of volume since the internal heating rate is also
given as specific heating rate (W/kg instead of W/m^3). In the non-
dimensional testcases this was not a problem so far, because the
density was commonly assumed to be 1 and constant.
>---------------------------------------------------------------
f627ac0812e6d6f2ed8409102917d78fa7b79485
source/postprocess/internal_heating_statistics.cc | 66 ++++++++++++++---------
1 file changed, 41 insertions(+), 25 deletions(-)
diff --git a/source/postprocess/internal_heating_statistics.cc b/source/postprocess/internal_heating_statistics.cc
index 908cd6b..5596b11 100644
--- a/source/postprocess/internal_heating_statistics.cc
+++ b/source/postprocess/internal_heating_statistics.cc
@@ -51,58 +51,74 @@ namespace aspect
update_quadrature_points |
update_JxW_values);
- std::vector<double> temperature_values (n_q_points);
- std::vector<double> pressure_values (n_q_points);
- std::vector<std::vector<double> > compositional_values (this->n_compositional_fields(),std::vector<double> (n_q_points));
- std::vector<double> composition_values_at_q_point (this->n_compositional_fields());
- std::vector<Point<dim> > quadrature_points(n_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());
+
+ in.strain_rate.resize(0); // we do not need the viscosity
+ std::vector<std::vector<double> > composition_values (this->n_compositional_fields(),std::vector<double> (quadrature_formula.size()));
typename DoFHandler<dim>::active_cell_iterator
cell = this->get_dof_handler().begin_active(),
endc = this->get_dof_handler().end();
double local_internal_heating_integrals = 0;
- double local_volume = 0;
+ double local_mass = 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.temperature].get_function_values (this->get_solution(),
- temperature_values);
- fe_values[this->introspection().extractors.pressure].get_function_values (this->get_solution(),
- pressure_values);
- for (unsigned c=0; c<this->n_compositional_fields(); c++)
- fe_values[this->introspection().extractors.compositional_fields[c]].get_function_values (this->get_solution(),
- compositional_values[c]);
- quadrature_points=fe_values.get_quadrature_points();
+ 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);
+ 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];
+ }
+
+ in.position = fe_values.get_quadrature_points();
+
+ this->get_material_model().evaluate(in, out);
for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned c=0; c<this->n_compositional_fields(); c++)
- composition_values_at_q_point[c]=compositional_values[c][q];
- local_internal_heating_integrals += heating_model.specific_heating_rate(temperature_values[q],
- pressure_values[q],
- composition_values_at_q_point,
- quadrature_points[q])*fe_values.JxW(q);
- local_volume+=fe_values.JxW(q);
+ in.composition[q][c] = composition_values[c][q];
+
+ local_internal_heating_integrals += heating_model.specific_heating_rate(in.temperature[q],
+ in.pressure[q],
+ in.composition[q],
+ in.position[q])
+ * out.densities[q] * fe_values.JxW(q);
+
+ local_mass += out.densities[q] * fe_values.JxW(q);
}
}
+
// compute the sum over all processors
std::vector<double> local_value;
std::vector<double> global_value;
local_value.push_back(local_internal_heating_integrals);
- local_value.push_back(local_volume);
+ local_value.push_back(local_mass);
+
Utilities::MPI::sum (local_value,
this->get_mpi_communicator(),
global_value);
- double global_internal_heating_integrals=global_value[0];
- double global_volume=global_value[1];
+ const double global_internal_heating_integrals = global_value[0];
+ const double global_mass = global_value[1];
// finally produce something for the statistics file
const std::string name1("Average internal heating rate (W/kg) ");
- statistics.add_value (name1, global_internal_heating_integrals/global_volume);
+ statistics.add_value (name1, global_internal_heating_integrals/global_mass);
// also make sure that the other columns filled by the this object
// all show up with sufficient accuracy and in scientific notation
statistics.set_precision (name1, 8);
@@ -120,7 +136,7 @@ namespace aspect
std::ostringstream output;
output.precision(4);
- output << global_internal_heating_integrals/global_volume << " W/kg, "
+ output << global_internal_heating_integrals/global_mass << " W/kg, "
<< global_internal_heating_integrals << " W";
return std::pair<std::string, std::string> ("Internal heating rate (average/total): ",
More information about the CIG-COMMITS
mailing list