[cig-commits] r1278 - in branches/active_compositions: include/aspect include/aspect/material_model source/material_model source/postprocess/visualization source/simulator
dannberg at dealii.org
dannberg at dealii.org
Mon Oct 15 15:40:08 PDT 2012
Author: dannberg
Date: 2012-10-15 16:40:07 -0600 (Mon, 15 Oct 2012)
New Revision: 1278
Modified:
branches/active_compositions/include/aspect/material_model/duretz_et_al.h
branches/active_compositions/include/aspect/material_model/interface.h
branches/active_compositions/include/aspect/material_model/simple.h
branches/active_compositions/include/aspect/material_model/steinberger.h
branches/active_compositions/include/aspect/material_model/table.h
branches/active_compositions/include/aspect/material_model/tan_gurnis.h
branches/active_compositions/include/aspect/simulator.h
branches/active_compositions/source/material_model/duretz_et_al.cc
branches/active_compositions/source/material_model/interface.cc
branches/active_compositions/source/material_model/simple.cc
branches/active_compositions/source/material_model/steinberger.cc
branches/active_compositions/source/material_model/table.cc
branches/active_compositions/source/material_model/tan_gurnis.cc
branches/active_compositions/source/postprocess/visualization/friction_heating.cc
branches/active_compositions/source/postprocess/visualization/viscosity.cc
branches/active_compositions/source/postprocess/visualization/viscosity_ratio.cc
branches/active_compositions/source/simulator/assembly.cc
branches/active_compositions/source/simulator/helper_functions.cc
branches/active_compositions/source/simulator/parameters.cc
branches/active_compositions/source/simulator/simulator_access.cc
Log:
Introduce composition dependence of viscosity and viscosity ratio for all material models. Change structure of compositional fields vectors in the material model (input). Write some general routines to extract a vector with the values of the compositional field, which is needed for the material model.
Modified: branches/active_compositions/include/aspect/material_model/duretz_et_al.h
===================================================================
--- branches/active_compositions/include/aspect/material_model/duretz_et_al.h 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/include/aspect/material_model/duretz_et_al.h 2012-10-15 22:40:07 UTC (rev 1278)
@@ -76,6 +76,7 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
@@ -244,6 +245,7 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
@@ -365,6 +367,7 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
Modified: branches/active_compositions/include/aspect/material_model/interface.h
===================================================================
--- branches/active_compositions/include/aspect/material_model/interface.h 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/include/aspect/material_model/interface.h 2012-10-15 22:40:07 UTC (rev 1278)
@@ -86,7 +86,7 @@
*/
/**
* Return the viscosity $\eta$ of the model as a function of temperature,
- * pressure, strain rate, and position.
+ * pressure, composition, strain rate, and position.
*
* @note The strain rate given as the third argument of this function
* is computed as $\varepsilon(\mathbf u)=\frac 12 (\nabla \mathbf u +
@@ -98,6 +98,7 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const = 0;
@@ -107,6 +108,7 @@
*/
virtual double viscosity_ratio (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strainrate,
const Point<dim> &position) const;
Modified: branches/active_compositions/include/aspect/material_model/simple.h
===================================================================
--- branches/active_compositions/include/aspect/material_model/simple.h 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/include/aspect/material_model/simple.h 2012-10-15 22:40:07 UTC (rev 1278)
@@ -53,6 +53,7 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
Modified: branches/active_compositions/include/aspect/material_model/steinberger.h
===================================================================
--- branches/active_compositions/include/aspect/material_model/steinberger.h 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/include/aspect/material_model/steinberger.h 2012-10-15 22:40:07 UTC (rev 1278)
@@ -56,6 +56,7 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
Modified: branches/active_compositions/include/aspect/material_model/table.h
===================================================================
--- branches/active_compositions/include/aspect/material_model/table.h 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/include/aspect/material_model/table.h 2012-10-15 22:40:07 UTC (rev 1278)
@@ -48,11 +48,13 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
virtual double viscosity_ratio (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
Modified: branches/active_compositions/include/aspect/material_model/tan_gurnis.h
===================================================================
--- branches/active_compositions/include/aspect/material_model/tan_gurnis.h 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/include/aspect/material_model/tan_gurnis.h 2012-10-15 22:40:07 UTC (rev 1278)
@@ -51,6 +51,7 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
Modified: branches/active_compositions/include/aspect/simulator.h
===================================================================
--- branches/active_compositions/include/aspect/simulator.h 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/include/aspect/simulator.h 2012-10-15 22:40:07 UTC (rev 1278)
@@ -362,6 +362,19 @@
const BoundaryTemperature::Interface<dim> &
get_boundary_temperature () const;
+ /**
+ * Copy the values of the compositional fields at the quadrature point
+ * q given as input parameter to the output vector
+ * composition_values_at_q_point.
+ *
+ * This function is implemented in
+ * <code>source/simulator/helper_functions.cc</code>.
+ */
+ void
+ get_composition_values_at_q_point (const std::vector<std::vector<double>> &composition_values,
+ const unsigned int q,
+ std::vector<double> &composition_values_at_q_point) const;
+
/** @} */
private:
@@ -536,7 +549,6 @@
*/
unsigned int n_compositional_fields;
std::vector<unsigned int> normalized_fields;
- std::vector<double> chemical_diffusivities;
/**
* @}
*/
@@ -1248,8 +1260,42 @@
const unsigned int composition_index,
double &max_residual,
double &max_velocity) const;
+ /**
+ * Extract the values of temperature, pressure, composition and
+ * optional strain rate for the current linearization point.
+ * These values are stored as input arguments for the material
+ * model. The compositional fields are extracted with the
+ * individual compositional fields as outer vectors and the values
+ * at each quadrature point as inner vectors, but the material
+ * model needs it the other way round. Hence, this vector of vectors
+ * is transposed.
+ *
+ * @param compute_strainrate If the strain rate should be computed
+ * or not
+ *
+ * This function is implemented in
+ * <code>source/simulator/assembly.cc</code>.
+ */
+ void
+ compute_material_model_input_values (const TrilinosWrappers::MPI::BlockVector ¤t_linearization_point,
+ const FEValues<dim,dim> &input_finite_element_values,
+ const bool compute_strainrate,
+ typename MaterialModel::Interface<dim>::MaterialModelInputs &material_model_inputs) const;
/**
+ * Copy the values of the compositional fields at the quadrature point
+ * q given as input parameter to the output vector
+ * composition_values_at_q_point.
+ *
+ * This function is implemented in
+ * <code>source/simulator/helper_functions.cc</code>.
+ */
+ void
+ extract_composition_values_at_q_point (const std::vector<std::vector<double>> &composition_values,
+ const unsigned int q,
+ std::vector<double> &composition_values_at_q_point) const;
+
+ /**
* Generate and output some statistics like timing information and memory consumption.
* Whether this function does anything or not is controlled through the
* variable aspect::output_parallel_statistics.
Modified: branches/active_compositions/source/material_model/duretz_et_al.cc
===================================================================
--- branches/active_compositions/source/material_model/duretz_et_al.cc 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/source/material_model/duretz_et_al.cc 2012-10-15 22:40:07 UTC (rev 1278)
@@ -38,6 +38,7 @@
SolCx<dim>::
viscosity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const SymmetricTensor<2,dim> &,
const Point<dim> &p) const
{
@@ -257,6 +258,7 @@
SolKz<dim>::
viscosity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const SymmetricTensor<2,dim> &,
const Point<dim> &p) const
{
@@ -420,6 +422,7 @@
Inclusion<dim>::
viscosity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const SymmetricTensor<2,dim> &,
const Point<dim> &p) const
{
Modified: branches/active_compositions/source/material_model/interface.cc
===================================================================
--- branches/active_compositions/source/material_model/interface.cc 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/source/material_model/interface.cc 2012-10-15 22:40:07 UTC (rev 1278)
@@ -21,6 +21,7 @@
#include <aspect/global.h>
+#include <aspect/simulator.h>
#include <aspect/material_model/interface.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/std_cxx1x/tuple.h>
@@ -171,6 +172,7 @@
Interface<dim>::
viscosity_ratio (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const
{
@@ -247,9 +249,9 @@
position.resize(n_points);
temperature.resize(n_points);
pressure.resize(n_points);
- composition.resize(n_comp);
- for(unsigned int i=0;i<n_comp;++i)
- composition[i].resize(n_points);
+ composition.resize(n_points);
+ for(unsigned int i=0;i<n_points;++i)
+ composition[i].resize(n_comp);
strain_rate.resize(n_points);
}
@@ -273,7 +275,7 @@
{
for (unsigned int i=0; i < in.temperature.size(); ++i)
{
- out.viscosities[i] = viscosity (in.temperature[i], in.pressure[i], in.strain_rate[i], in.position[i]);
+ out.viscosities[i] = viscosity (in.temperature[i], in.pressure[i], in.composition[i], in.strain_rate[i], in.position[i]);
out.densities[i] = density (in.temperature[i], in.pressure[i], in.position[i]);
out.thermal_expansion_coefficients[i] = thermal_expansion_coefficient (in.temperature[i], in.pressure[i], in.position[i]);
out.seismic_Vp[i] = seismic_Vp (in.temperature[i], in.pressure[i], in.position[i]);
Modified: branches/active_compositions/source/material_model/simple.cc
===================================================================
--- branches/active_compositions/source/material_model/simple.cc 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/source/material_model/simple.cc 2012-10-15 22:40:07 UTC (rev 1278)
@@ -35,10 +35,12 @@
Simple<dim>::
viscosity (const double,
const double,
+ const std::vector<double> &composition, /*composition*/
const SymmetricTensor<2,dim> &,
const Point<dim> &) const
{
- return eta;
+// return eta;
+ return (4.0*composition[0]+1) * eta;
}
@@ -141,11 +143,10 @@
{
for (unsigned int i=0; i < in.temperature.size(); ++i)
{
- out.viscosities[i] = viscosity (in.temperature[i], in.pressure[i], in.strain_rate[i], in.position[i]);
-// out.viscosities[i] = (29.0*in.composition[0][i]+1) * eta;
+ out.viscosities[i] = viscosity (in.temperature[i], in.pressure[i], in.composition[i], in.strain_rate[i], in.position[i]);
out.densities[i] = (this->n_compositional_fields()>0)
?
- (100*in.composition[0][i]) + reference_rho *(1 - thermal_alpha * (in.temperature[i] - reference_T))
+ (100*in.composition[i][0]) + reference_rho *(1 - thermal_alpha * (in.temperature[i] - reference_T))
:
density (in.temperature[i], in.pressure[i], in.position[i]);
out.thermal_expansion_coefficients[i] = thermal_expansion_coefficient (in.temperature[i], in.pressure[i], in.position[i]);
Modified: branches/active_compositions/source/material_model/steinberger.cc
===================================================================
--- branches/active_compositions/source/material_model/steinberger.cc 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/source/material_model/steinberger.cc 2012-10-15 22:40:07 UTC (rev 1278)
@@ -391,6 +391,7 @@
Steinberger<dim>::
viscosity (const double temperature,
const double,
+ const std::vector<double> &, /*composition*/
const SymmetricTensor<2,dim> &,
const Point<dim> &position) const
{
Modified: branches/active_compositions/source/material_model/table.cc
===================================================================
--- branches/active_compositions/source/material_model/table.cc 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/source/material_model/table.cc 2012-10-15 22:40:07 UTC (rev 1278)
@@ -347,6 +347,7 @@
Table<dim>::
viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &, /*composition*/
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const
{
@@ -413,6 +414,7 @@
Table<dim>::
viscosity_ratio (const double temperature,
const double pressure,
+ const std::vector<double> &, /*composition*/
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const
{
Modified: branches/active_compositions/source/material_model/tan_gurnis.cc
===================================================================
--- branches/active_compositions/source/material_model/tan_gurnis.cc 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/source/material_model/tan_gurnis.cc 2012-10-15 22:40:07 UTC (rev 1278)
@@ -50,6 +50,7 @@
TanGurnis<dim>::
viscosity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const SymmetricTensor<2,dim> &,
const Point<dim> &pos) const
{
Modified: branches/active_compositions/source/postprocess/visualization/friction_heating.cc
===================================================================
--- branches/active_compositions/source/postprocess/visualization/friction_heating.cc 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/source/postprocess/visualization/friction_heating.cc 2012-10-15 22:40:07 UTC (rev 1278)
@@ -82,6 +82,7 @@
computed_quantities[q](0) = 2 * this->get_material_model().viscosity(temperature,
pressure,
+ composition,
strain_rate,
evaluation_points[q]) *
compressible_strain_rate * compressible_strain_rate;
Modified: branches/active_compositions/source/postprocess/visualization/viscosity.cc
===================================================================
--- branches/active_compositions/source/postprocess/visualization/viscosity.cc 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/source/postprocess/visualization/viscosity.cc 2012-10-15 22:40:07 UTC (rev 1278)
@@ -82,6 +82,7 @@
computed_quantities[q](0) = this->get_material_model().viscosity(temperature,
pressure,
+ composition,
strain_rate,
evaluation_points[q]);
}
Modified: branches/active_compositions/source/postprocess/visualization/viscosity_ratio.cc
===================================================================
--- branches/active_compositions/source/postprocess/visualization/viscosity_ratio.cc 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/source/postprocess/visualization/viscosity_ratio.cc 2012-10-15 22:40:07 UTC (rev 1278)
@@ -82,6 +82,7 @@
computed_quantities[q](0) = this->get_material_model().viscosity_ratio(temperature,
pressure,
+ composition,
strain_rate,
evaluation_points[q]);
}
Modified: branches/active_compositions/source/simulator/assembly.cc
===================================================================
--- branches/active_compositions/source/simulator/assembly.cc 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/source/simulator/assembly.cc 2012-10-15 22:40:07 UTC (rev 1278)
@@ -677,47 +677,47 @@
const unsigned int n_q_points = old_temperature.size();
- typename MaterialModel::Interface<dim>::MaterialModelInputs material_model_inputs (n_q_points, parameters.n_compositional_fields);
- typename MaterialModel::Interface<dim>::MaterialModelOutputs material_model_outputs (n_q_points);
+ typename MaterialModel::Interface<dim>::MaterialModelInputs inputs (n_q_points, parameters.n_compositional_fields);
+ typename MaterialModel::Interface<dim>::MaterialModelOutputs outputs (n_q_points);
for (unsigned int q=0; q<n_q_points; ++q) {
- material_model_inputs.temperature[q] = (old_temperature[q] + old_old_temperature[q]) / 2;
- material_model_inputs.position[q] = evaluation_points[q];
- material_model_inputs.pressure[q] = (old_pressure[q] + old_old_pressure[q]) / 2;
+ inputs.temperature[q] = (old_temperature[q] + old_old_temperature[q]) / 2;
+ inputs.position[q] = evaluation_points[q];
+ inputs.pressure[q] = (old_pressure[q] + old_old_pressure[q]) / 2;
for (unsigned int i=0; i<parameters.n_compositional_fields; ++i)
- material_model_inputs.composition[i][q] = (old_composition[i][q] + old_old_composition[i][q]) / 2;
- material_model_inputs.strain_rate[q] = (old_strain_rates[q] + old_old_strain_rates[q]) / 2;
+ inputs.composition[q][i] = (old_composition[i][q] + old_old_composition[i][q]) / 2;
+ inputs.strain_rate[q] = (old_strain_rates[q] + old_old_strain_rates[q]) / 2;
}
- material_model->compute_parameters(material_model_inputs,material_model_outputs);
+ material_model->compute_parameters(inputs,outputs);
for (unsigned int q=0; q < n_q_points; ++q)
{
const Tensor<1,dim> u = (old_velocity_values[q] +
old_old_velocity_values[q]) / 2;
- const SymmetricTensor<2,dim> strain_rate = material_model_inputs.strain_rate[q];
+ const SymmetricTensor<2,dim> strain_rate = inputs.strain_rate[q];
const double dT_dt = (old_temperature[q] - old_old_temperature[q])
/ old_time_step;
const double u_grad_T = u * (old_temperature_grads[q] +
old_old_temperature_grads[q]) / 2;
- const double alpha = material_model_outputs.thermal_expansion_coefficients[q];
- const double density = material_model_outputs.densities[q];
- const double thermal_conductivity = material_model_outputs.thermal_conductivities[q];
- const double c_P = material_model_outputs.specific_heat[q];
+ const double alpha = outputs.thermal_expansion_coefficients[q];
+ const double density = outputs.densities[q];
+ const double thermal_conductivity = outputs.thermal_conductivities[q];
+ const double c_P = outputs.specific_heat[q];
const double k_Delta_T = thermal_conductivity
* (old_temperature_laplacians[q] +
old_old_temperature_laplacians[q]) / 2;
// verify correctness of the heating term
- const double viscosity = material_model_outputs.viscosities[q];
- const bool is_compressible = material_model_outputs.is_compressible;
+ const double viscosity = outputs.viscosities[q];
+ const bool is_compressible = outputs.is_compressible;
const double compressibility
= (is_compressible
?
- material_model_outputs.compressibilities[q]
+ outputs.compressibilities[q]
:
std::numeric_limits<double>::quiet_NaN() );
const Tensor<1,dim> gravity = gravity_model->gravity_vector (evaluation_points[q] );
@@ -758,14 +758,14 @@
// alpha = - 1/rho drho/dT
(parameters.include_adiabatic_heating
?
- alpha * density * (u*gravity) * material_model_inputs.temperature[q]
+ alpha * density * (u*gravity) * inputs.temperature[q]
:
0)
);
double residual
= std::abs(density * c_P * (dT_dt + u_grad_T) - k_Delta_T - gamma);
if (parameters.stabilization_alpha == 2)
- residual *= std::abs(material_model_inputs.temperature[q] - average_temperature);
+ residual *= std::abs(inputs.temperature[q] - average_temperature);
max_residual = std::max (residual, max_residual);
max_velocity = std::max (std::sqrt (u*u), max_velocity);
@@ -802,7 +802,9 @@
const double u_grad_C = u * (old_composition_grads[q] +
old_old_composition_grads[q]) / 2;
- const double kappa = parameters.chemical_diffusivities[composition_index];
+ // the chemical diffusivity kappa is always smaller than the numerical diffusivity
+ // therefore we set it to zero
+ const double kappa = 0;
const double kappa_Delta_C = kappa
* (old_composition_laplacians[q] +
@@ -919,7 +921,54 @@
}
+ template <int dim>
+ void
+ Simulator<dim>::
+ compute_material_model_input_values (const TrilinosWrappers::MPI::BlockVector ¤t_linearization_point,
+ const FEValues<dim> &input_finite_element_values,
+ const bool compute_strainrate,
+ typename MaterialModel::Interface<dim>::MaterialModelInputs &material_model_inputs) const
+ {
+ const FEValuesExtractors::Vector input_velocities (0);
+ const FEValuesExtractors::Scalar input_pressure (dim);
+ const FEValuesExtractors::Scalar input_temperature (dim+1);
+ std::vector<FEValuesExtractors::Scalar> input_composition;
+
+ unsigned int n_q_points = material_model_inputs.temperature.size();
+
+ for (unsigned int q=0;q<parameters.n_compositional_fields;++q)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+q);
+ input_composition.push_back(temp);
+ }
+
+ input_finite_element_values[input_temperature].get_function_values (current_linearization_point,
+ material_model_inputs.temperature);
+ input_finite_element_values[input_pressure].get_function_values(current_linearization_point,
+ material_model_inputs.pressure);
+ if(compute_strainrate)
+ input_finite_element_values[input_velocities].get_function_symmetric_gradients(current_linearization_point,
+ material_model_inputs.strain_rate);
+
+ // 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>> composition_values (parameters.n_compositional_fields,
+ std::vector<double> (n_q_points));
+
+ for(unsigned int i=0;i<parameters.n_compositional_fields;++i)
+ input_finite_element_values[input_composition[i]].get_function_values(current_linearization_point,
+ composition_values[i]);
+
+ // then we copy these values to exchange the inner and outer vector, because for the material
+ // model we need a vector with values of all the compositional fields for every quadrature point
+ for(unsigned int q=0;q<n_q_points;++q)
+ for(unsigned int i=0;i<parameters.n_compositional_fields;++i)
+ material_model_inputs.composition[q][i] = composition_values[i][q];
+
+ }
+
+
template <int dim>
void
Simulator<dim>::
@@ -932,43 +981,22 @@
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
- const FEValuesExtractors::Scalar temperature (dim+1);
- std::vector<FEValuesExtractors::Scalar> compositional_fields;
- for (unsigned int q=0;q<parameters.n_compositional_fields;++q)
- {
- const FEValuesExtractors::Scalar temp(dim+2+q);
- compositional_fields.push_back(temp);
- }
-
scratch.finite_element_values.reinit (cell);
- scratch.finite_element_values[temperature].get_function_values (current_linearization_point,
- scratch.temperature_values);
- scratch.finite_element_values[pressure].get_function_values(current_linearization_point,
- scratch.pressure_values);
- scratch.finite_element_values[velocities].get_function_symmetric_gradients(current_linearization_point,
- scratch.strain_rates);
- for(unsigned int q=0;q<parameters.n_compositional_fields;++q)
- scratch.finite_element_values[compositional_fields[q]].get_function_values(current_linearization_point,
- scratch.composition_values[q]);
-
data.local_matrix = 0;
- scratch.material_model_inputs.temperature = scratch.temperature_values;
- for (unsigned int q=0; q<n_q_points; ++q)
- scratch.material_model_inputs.position[q] = scratch.finite_element_values.quadrature_point(q);
- scratch.material_model_inputs.pressure = scratch.pressure_values;
- for (unsigned int i=0; i<parameters.n_compositional_fields; ++i)
- scratch.material_model_inputs.composition[i] = scratch.composition_values[i];
- scratch.material_model_inputs.strain_rate = scratch.strain_rates;
+ compute_material_model_input_values (current_linearization_point,
+ scratch.finite_element_values,
+ true,
+ scratch.material_model_inputs);
material_model->compute_parameters(scratch.material_model_inputs,scratch.material_model_outputs);
for (unsigned int q=0; q<n_q_points; ++q)
{
- const double current_temperature = scratch.temperature_values[q];
- const double current_pressure = scratch.pressure_values[q];
+ const double current_temperature = scratch.material_model_inputs.temperature[q];
+ const double current_pressure = scratch.material_model_inputs.pressure[q];
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
@@ -1110,57 +1138,31 @@
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
- const FEValuesExtractors::Scalar temperature (dim+1);
- std::vector<FEValuesExtractors::Scalar> compositional_fields;
- for (unsigned int q=0;q<parameters.n_compositional_fields;++q)
- {
- const FEValuesExtractors::Scalar temp(dim+2+q);
- compositional_fields.push_back(temp);
- }
-
scratch.finite_element_values.reinit (cell);
- //scratch.finite_element_values[temperature].get_function_values (old_solution,
- // scratch.old_temperature_values);
- // Assuming we already have the temperature for the current time step:
- scratch.finite_element_values[temperature].get_function_values (current_linearization_point,
- scratch.temperature_values);
- scratch.finite_element_values[pressure].get_function_values(current_linearization_point,
- scratch.pressure_values);
- scratch.finite_element_values[velocities].get_function_values(current_linearization_point,
- scratch.velocity_values);
- for(unsigned int q=0;q<parameters.n_compositional_fields;++q)
- scratch.finite_element_values[compositional_fields[q]].get_function_values(current_linearization_point,
- scratch.composition_values[q]);
- // we only need the strain rates for the viscosity,
- // which we only need when rebuilding the matrix
if (rebuild_stokes_matrix)
- scratch.finite_element_values[velocities]
- .get_function_symmetric_gradients(current_linearization_point,
- scratch.strain_rates);
-
- if (rebuild_stokes_matrix)
data.local_matrix = 0;
data.local_rhs = 0;
if (material_model->is_compressible())
data.local_pressure_shape_function_integrals = 0;
- scratch.material_model_inputs.temperature = scratch.temperature_values;
- for (unsigned int q=0; q<n_q_points; ++q)
- scratch.material_model_inputs.position[q] = scratch.finite_element_values.quadrature_point(q);
- scratch.material_model_inputs.pressure = scratch.pressure_values;
- for (unsigned int i=0; i<parameters.n_compositional_fields; ++i)
- scratch.material_model_inputs.composition[i] = scratch.composition_values[i];
- scratch.material_model_inputs.strain_rate = scratch.strain_rates;
+ // we only need the strain rates for the viscosity,
+ // which we only need when rebuilding the matrix
+ compute_material_model_input_values (current_linearization_point,
+ scratch.finite_element_values,
+ rebuild_stokes_matrix,
+ scratch.material_model_inputs);
material_model->compute_parameters(scratch.material_model_inputs,scratch.material_model_outputs);
+ scratch.finite_element_values[velocities].get_function_values(current_linearization_point,
+ scratch.velocity_values);
for (unsigned int q=0; q<n_q_points; ++q)
{
- const double current_temperature = scratch.temperature_values[q];
- const double current_pressure = scratch.pressure_values[q];
+ const double current_temperature = scratch.material_model_inputs.temperature[q];
+ const double current_pressure = scratch.material_model_inputs.pressure[q];
const Tensor<1,dim> current_u = scratch.velocity_values[q];
for (unsigned int k=0; k<dofs_per_cell; ++k)
@@ -1387,18 +1389,10 @@
scratch.finite_element_values[pressure].get_function_values (old_old_solution,
scratch.old_old_pressure);
- scratch.finite_element_values[temperature].get_function_values(current_linearization_point,
- scratch.current_temperature_values);
scratch.finite_element_values[velocities].get_function_values(current_linearization_point,
scratch.current_velocity_values);
- scratch.finite_element_values[velocities].get_function_symmetric_gradients(current_linearization_point,
- scratch.current_strain_rates);
- scratch.finite_element_values[pressure].get_function_values(current_linearization_point,
- scratch.current_pressure_values);
for(unsigned int q=0;q<parameters.n_compositional_fields;++q) {
- scratch.finite_element_values[compositional_fields[q]].get_function_values(current_linearization_point,
- scratch.current_composition_values[q]);
scratch.finite_element_values[compositional_fields[q]].get_function_values(old_solution,
scratch.old_composition_values[q]);
scratch.finite_element_values[compositional_fields[q]].get_function_values(old_old_solution,
@@ -1432,13 +1426,10 @@
cell->diameter(),
0); //index for temperature field
- scratch.material_model_inputs.temperature = scratch.current_temperature_values;
- for (unsigned int q=0; q<n_q_points; ++q)
- scratch.material_model_inputs.position[q] = scratch.finite_element_values.quadrature_point(q);
- scratch.material_model_inputs.pressure = scratch.current_pressure_values;
- for (unsigned int i=0; i<parameters.n_compositional_fields; ++i)
- scratch.material_model_inputs.composition[i] = scratch.current_composition_values[i];
- scratch.material_model_inputs.strain_rate = scratch.current_strain_rates;
+ compute_material_model_input_values (current_linearization_point,
+ scratch.finite_element_values,
+ true,
+ scratch.material_model_inputs);
material_model->compute_parameters(scratch.material_model_inputs,scratch.material_model_outputs);
@@ -1461,10 +1452,10 @@
:
scratch.old_temperature_values[q]);
- const double current_T = scratch.current_temperature_values[q];
- const SymmetricTensor<2,dim> current_strain_rate = scratch.current_strain_rates[q];
+ const double current_T = scratch.material_model_inputs.temperature[q];
+ const SymmetricTensor<2,dim> current_strain_rate = scratch.material_model_inputs.strain_rate[q];
const Tensor<1,dim> current_u = scratch.current_velocity_values[q];
- const double current_p = scratch.current_pressure_values[q];
+ const double current_p = scratch.material_model_inputs.pressure[q];
const double alpha = scratch.material_model_outputs.thermal_expansion_coefficients[q];
const double density = scratch.material_model_outputs.densities[q];
@@ -1729,7 +1720,9 @@
const Tensor<1,dim> current_u = scratch.current_velocity_values[q];
- const double kappa = parameters.chemical_diffusivities[composition_index];
+ // the chemical diffusivity kappa is always smaller than the numerical diffusivity
+ // therefore we set it to zero
+ const double kappa = 0.0;
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
Modified: branches/active_compositions/source/simulator/helper_functions.cc
===================================================================
--- branches/active_compositions/source/simulator/helper_functions.cc 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/source/simulator/helper_functions.cc 2012-10-15 22:40:07 UTC (rev 1278)
@@ -254,7 +254,16 @@
}
+ template <int dim>
+ void Simulator<dim>::extract_composition_values_at_q_point (const std::vector<std::vector<double>> &composition_values,
+ const unsigned int q,
+ std::vector<double> &composition_values_at_q_point) const
+ {
+ for(unsigned int k=0; k < composition_values_at_q_point.size(); ++k)
+ composition_values_at_q_point[k] = composition_values[k][q];
+ }
+
template <int dim>
std::pair<double,double>
Simulator<dim>::get_extrapolated_temperature_or_composition_range (const unsigned int index) const
@@ -616,10 +625,19 @@
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
const FEValuesExtractors::Scalar temperature (dim+1);
+ std::vector<FEValuesExtractors::Scalar> compositional_fields;
+ for (unsigned int q=0;q<parameters.n_compositional_fields;++q)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+q);
+ compositional_fields.push_back(temp);
+ }
+
std::vector<SymmetricTensor<2,dim> > strain_rates(n_q_points);
std::vector<double> pressure_values(n_q_points);
std::vector<double> temperature_values(n_q_points);
+ std::vector<std::vector<double>> composition_values (parameters.n_compositional_fields,std::vector<double> (n_q_points));
+ std::vector<double> composition_values_at_q_point (parameters.n_compositional_fields);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
@@ -637,6 +655,9 @@
temperature_values);
fe_values[velocities].get_function_symmetric_gradients (this->solution,
strain_rates);
+ for(unsigned int i=0;i<parameters.n_compositional_fields;++i)
+ fe_values[compositional_fields[i]].get_function_values(this->solution,
+ composition_values[i]);
for (unsigned int q = 0; q < n_q_points; ++q)
{
@@ -644,8 +665,13 @@
const unsigned int idx = static_cast<unsigned int>((depth*num_slices)/max_depth);
Assert(idx<num_slices, ExcInternalError());
+ extract_composition_values_at_q_point (composition_values,
+ q,
+ composition_values_at_q_point);
+
const double viscosity = material_model->viscosity(temperature_values[q],
pressure_values[q],
+ composition_values_at_q_point,
strain_rates[q],
fe_values.quadrature_point(q));
++counts[idx];
@@ -1021,6 +1047,9 @@
template void Simulator<dim>::normalize_pressure(LinearAlgebra::BlockVector &vector); \
template double Simulator<dim>::get_maximal_velocity (const LinearAlgebra::BlockVector &solution) const; \
template std::pair<double,double> Simulator<dim>::get_extrapolated_temperature_or_composition_range (const unsigned int index) const; \
+ template void Simulator<dim>::extract_composition_values_at_q_point (const std::vector<std::vector<double>> &composition_values, \
+ const unsigned int q, \
+ std::vector<double> &composition_values_at_q_point) const; \
template double Simulator<dim>::compute_time_step () const; \
template void Simulator<dim>::make_pressure_rhs_compatible(LinearAlgebra::BlockVector &vector); \
template void Simulator<dim>::compute_depth_average_field(const unsigned int index, std::vector<double> &values) const; \
Modified: branches/active_compositions/source/simulator/parameters.cc
===================================================================
--- branches/active_compositions/source/simulator/parameters.cc 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/source/simulator/parameters.cc 2012-10-15 22:40:07 UTC (rev 1278)
@@ -421,11 +421,6 @@
"at every point and the global maximum is determined."
"Second, the compositional fields to be normalized are "
"divided by this maximum.");
- prm.declare_entry ("List of chemical diffusivities", "",
- Patterns::List (Patterns::Double(0)),
- "A list of doubles with as many entries as there are"
- "compositional fields (see Number of fields). These are"
- "the chemical diffusivities of each compositional field. ");
}
prm.leave_subsection ();
}
@@ -599,16 +594,9 @@
(Utilities::split_string_list(prm.get ("List of normalized fields")));
normalized_fields = std::vector<unsigned int> (n_normalized_fields.begin(),
n_normalized_fields.end());
- const std::vector<double> n_chemical_diffusivities = Utilities::string_to_double
- (Utilities::split_string_list(prm.get ("List of chemical diffusivities")));
- chemical_diffusivities = std::vector<double> (n_chemical_diffusivities.begin(),
- n_chemical_diffusivities.end());
Assert (normalized_fields.size() <= n_compositional_fields,
ExcMessage("Invalid input parameter file: Too many entries in List of normalized fields"));
- Assert (chemical_diffusivities.size() == n_compositional_fields,
- ExcMessage("Invalid input parameter file: number of chemical diffusivities must be equal "
- "to the number of compositional fields!"));
}
prm.leave_subsection ();
}
Modified: branches/active_compositions/source/simulator/simulator_access.cc
===================================================================
--- branches/active_compositions/source/simulator/simulator_access.cc 2012-10-14 07:51:45 UTC (rev 1277)
+++ branches/active_compositions/source/simulator/simulator_access.cc 2012-10-15 22:40:07 UTC (rev 1278)
@@ -287,6 +287,18 @@
{
return *simulator->adiabatic_conditions.get();
}
+
+ template <int dim>
+ void
+ SimulatorAccess<dim>::get_composition_values_at_q_point (const std::vector<std::vector<double>> &composition_values,
+ const unsigned int q,
+ std::vector<double> &composition_values_at_q_point) const
+ {
+ simulator->extract_composition_values_at_q_point (composition_values,
+ q,
+ composition_values_at_q_point);
+ }
+
}
More information about the CIG-COMMITS
mailing list