[cig-commits] commit 2064 by buerg to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Wed Nov 27 08:43:55 PST 2013
Revision 2064
Add viscosity stabilisation for temperature equation.
U trunk/aspire/include/aspect/simulator.h
U trunk/aspire/source/material_model/interface.cc
U trunk/aspire/source/simulator/assembly.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=2064&peg=2064
Diff:
Modified: trunk/aspire/include/aspect/simulator.h
===================================================================
--- trunk/aspire/include/aspect/simulator.h 2013-11-26 23:01:30 UTC (rev 2063)
+++ trunk/aspire/include/aspect/simulator.h 2013-11-27 16:43:09 UTC (rev 2064)
@@ -864,7 +864,48 @@
* <code>source/simulator/helper_functions.cc</code>.
*/
void output_statistics();
+
/**
+ * Compute the heating term for the temperature system.
+ *
+ * This function is implemented in
+ * <code>source/simulator/assembly.cc</code>.
+ */
+ double compute_heating_term (const internal::Assembly::Scratch::TemperatureSystem<dim>& scratch,
+ const unsigned int q) const;
+
+ /**
+ * Compute the residual of the temperature equation to be used for
+ * the artificial diffusion coefficient value on a cell given the
+ * values and gradients of the solution passed as arguments.
+ *
+ * This function is implemented in
+ * <code>source/simulator/assembly.cc</code>.
+ */
+ void
+ compute_temperature_system_residual (const internal::Assembly::Scratch::TemperatureSystem<dim>& scratch,
+ const double average_temperature,
+ double& max_residual,
+ double& max_velocity,
+ double& max_density,
+ double& max_specific_heat) const;
+
+ /**
+ * Compute the residual of one advection equation to be used for the
+ * artificial diffusion coefficient value on a cell given the values
+ * and gradients of the solution passed as arguments.
+ *
+ * This function is implemented in
+ * <code>source/simulator/assembly.cc</code>.
+ */
+ double
+ compute_viscosity (internal::Assembly::Scratch::TemperatureSystem<dim>& scratch,
+ const double global_max_velocity,
+ const double global_field_variation,
+ const double average_field,
+ const double global_entropy_variation,
+ const double cell_diameter) const;
+ /**
* @}
*/
Modified: trunk/aspire/source/material_model/interface.cc
===================================================================
--- trunk/aspire/source/material_model/interface.cc 2013-11-26 23:01:30 UTC (rev 2063)
+++ trunk/aspire/source/material_model/interface.cc 2013-11-27 16:43:09 UTC (rev 2064)
@@ -317,12 +317,12 @@
template <int dim>
Interface<dim>::MaterialModelInputs::MaterialModelInputs(unsigned int n_points, unsigned int n_comp)
{
- position.resize(n_points);
- temperature.resize(n_points);
- composition.resize(n_points);
+ position.resize (n_points);
+ temperature.resize (n_points);
+ composition.resize (n_points);
for (unsigned int i=0; i<n_points; ++i)
- composition[i].resize(n_comp);
+ composition[i].resize (n_comp);
}
template <int dim>
Modified: trunk/aspire/source/simulator/assembly.cc
===================================================================
--- trunk/aspire/source/simulator/assembly.cc 2013-11-26 23:01:30 UTC (rev 2063)
+++ trunk/aspire/source/simulator/assembly.cc 2013-11-27 16:43:09 UTC (rev 2064)
@@ -189,24 +189,29 @@
const unsigned int n_compositional_fields);
TemperatureSystem (const TemperatureSystem &data);
- FEValues<dim> fe_values_velocity;
- FEValues<dim> fe_values_advection;
+ FEValues<dim> fe_values_velocity;
+ FEValues<dim> fe_values_advection;
- std::vector<double> phi_field;
- std::vector<Tensor<1, dim> > grad_phi_field;
+ std::vector<double> phi_field;
+ std::vector<Tensor<1, dim> > grad_phi_field;
- std::vector<double> current_temperature_values;
- std::vector<double> old_temperature_values;
- std::vector<double> old_old_temperature_values;
- std::vector<Tensor<1, dim> > old_temperature_grads;
- std::vector<Tensor<1, dim> > old_old_temperature_grads;
+ std::vector<double> current_temperature_values;
+ std::vector<double> old_temperature_values;
+ std::vector<double> old_old_temperature_values;
+ std::vector<double> old_temperature_laplacians;
+ std::vector<double> old_old_temperature_laplacians;
+ std::vector<Tensor<1, dim> > old_temperature_grads;
+ std::vector<Tensor<1, dim> > old_old_temperature_grads;
- std::vector<Tensor<1, dim> > current_velocity_values;
- std::vector<Tensor<1, dim> > old_velocity_values;
- std::vector<Tensor<1, dim> > old_old_velocity_values;
+ std::vector<Tensor<1, dim> > current_velocity_values;
+ std::vector<Tensor<1, dim> > old_velocity_values;
+ std::vector<Tensor<1, dim> > old_old_velocity_values;
+ std::vector<SymmetricTensor<2, dim> > current_strain_rates;
- std::vector<std::vector<double> > old_composition_values;
- std::vector<std::vector<double> > old_old_composition_values;
+ std::vector<std::vector<double> > old_composition_values;
+ std::vector<std::vector<double> > old_old_composition_values;
+ std::vector<std::vector<Tensor<1, dim> > > old_composition_grads;
+ std::vector<std::vector<Tensor<1, dim> > > old_old_composition_grads;
typename MaterialModel::Interface<dim>::MaterialModelInputs material_model_inputs;
typename MaterialModel::Interface<dim>::MaterialModelOutputs material_model_outputs;
@@ -256,8 +261,7 @@
:
fe_values_velocity (mapping,
fe_velocity, quadrature,
- update_values | update_gradients |
- update_quadrature_points),
+ update_values | update_gradients),
fe_values_advection (mapping,
fe_advection, quadrature,
update_gradients | update_JxW_values |
@@ -268,13 +272,18 @@
current_temperature_values(quadrature.size ()),
old_temperature_values (quadrature.size ()),
old_old_temperature_values (quadrature.size ()),
+ old_temperature_laplacians (quadrature.size ()),
+ old_old_temperature_laplacians (quadrature.size ()),
old_temperature_grads (quadrature.size ()),
old_old_temperature_grads (quadrature.size ()),
current_velocity_values (quadrature.size ()),
old_velocity_values (quadrature.size ()),
old_old_velocity_values (quadrature.size ()),
+ current_strain_rates (quadrature.size ()),
old_composition_values (n_compositional_fields, std::vector<double> (quadrature.size ())),
old_old_composition_values (n_compositional_fields, std::vector<double> (quadrature.size ())),
+ old_composition_grads (n_compositional_fields, std::vector<Tensor<1, dim> > (quadrature.size ())),
+ old_old_composition_grads (n_compositional_fields, std::vector<Tensor<1, dim> > (quadrature.size ())),
material_model_inputs (quadrature.size (), n_compositional_fields),
material_model_outputs (quadrature.size (), n_compositional_fields),
explicit_material_model_inputs (quadrature.size (), n_compositional_fields),
@@ -333,13 +342,18 @@
current_temperature_values (scratch.current_temperature_values),
old_temperature_values (scratch.old_temperature_values),
old_old_temperature_values (scratch.old_old_temperature_values),
+ old_temperature_laplacians (scratch.old_temperature_values),
+ old_old_temperature_laplacians (scratch.old_old_temperature_values),
old_temperature_grads (scratch.old_temperature_grads),
old_old_temperature_grads (scratch.old_old_temperature_grads),
current_velocity_values (scratch.current_velocity_values),
old_velocity_values (scratch.old_velocity_values),
old_old_velocity_values (scratch.old_old_velocity_values),
+ current_strain_rates (scratch.current_strain_rates),
old_composition_values (scratch.old_composition_values),
old_old_composition_values (scratch.old_old_composition_values),
+ old_composition_grads (scratch.old_composition_grads),
+ old_old_composition_grads (scratch.old_old_composition_grads),
material_model_inputs (scratch.material_model_inputs),
material_model_outputs (scratch.material_model_outputs),
explicit_material_model_inputs (scratch.explicit_material_model_inputs),
@@ -1136,7 +1150,127 @@
computing_timer.exit_section();
}
+
+ template<int dim>
+ double Simulator<dim>::
+ compute_heating_term (const internal::Assembly::Scratch::TemperatureSystem<dim>& scratch,
+ const unsigned int q) const
+ {
+ const SymmetricTensor<2, dim> current_strain_rate = scratch.current_strain_rates[q];
+ const Tensor<1, dim> current_u = scratch.current_velocity_values[q];
+ const double density = scratch.explicit_material_model_outputs.densities[q];
+ const double viscosity = scratch.explicit_material_model_outputs.viscosities[q];
+ double compressibility = scratch.explicit_material_model_outputs.density_changes[q][0]
+ * (current_u * (scratch.old_temperature_grads[q] + scratch.old_old_temperature_grads[q]));
+
+ for (unsigned int c = 0; c < parameters.n_compositional_fields; ++c)
+ compressibility += scratch.explicit_material_model_outputs.density_changes[q][c + 1]
+ * (current_u * (scratch.old_composition_grads[c][q] + scratch.old_old_composition_grads[c][q]));
+
+ const double gamma = // add the term 2*eta*(eps - 1/3*(tr eps)1):(eps - 1/3*(tr eps)1)
+ //
+ // we can multiply this out to obtain
+ // 2*eta*(eps:eps - 1/3*(tr eps)^2)
+ // and can then use that in the compressible case we have
+ // tr eps = div u
+ // = -1/rho u . grad rho
+ // and by the usual approximation we make,
+ // tr eps = -1/rho u . (drho/dT grad T
+ // + drho/dc_0 grad c_0
+ // + ...
+ // + drho/dc_n grad c_n)
+ (parameters.include_shear_heating ?
+ 2.0 * viscosity * (current_strain_rate * current_strain_rate
+ - (material_model->is_compressible () ?
+ std::pow (0.5 * compressibility / density, 2.0) :
+ 0.0)) :
+ 0.0);
+
+ return gamma;
+ }
+ template<int dim>
+ void Simulator<dim>::
+ compute_temperature_system_residual (const internal::Assembly::Scratch::TemperatureSystem<dim>& scratch,
+ const double average_temperature,
+ double& max_residual,
+ double& max_velocity,
+ double& max_density,
+ double& max_specific_heat) const
+ {
+ const unsigned int n_q_points = scratch.old_temperature_values.size ();
+
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ {
+ const Tensor<1, dim> u = 0.5 * (scratch.old_velocity_values[q] + scratch.old_old_velocity_values[q]);
+ const double dT_dt = (scratch.old_temperature_values[q] - scratch.old_old_temperature_values[q]) / old_time_step;
+ const double u_grad_T = 0.5 * (u * (scratch.old_temperature_grads[q] + scratch.old_old_temperature_grads[q]));
+ const double density = scratch.explicit_material_model_outputs.densities[q];
+ const double conductivity = scratch.explicit_material_model_outputs.thermal_conductivities[q];
+ const double c_P = scratch.explicit_material_model_outputs.specific_heat[q];
+ const double k_Delta_T = 0.5 * conductivity * (scratch.old_temperature_laplacians[q]
+ + scratch.old_old_temperature_laplacians[q]);
+ const double temperature = 0.5 * (scratch.old_temperature_values[q] + scratch.old_old_temperature_values[q]);
+ const double gamma = compute_heating_term (scratch, q);
+ double residual = std::abs (density * c_P * (dT_dt + u_grad_T) - k_Delta_T - gamma);
+
+ if (parameters.stabilization_alpha == 2)
+ residual *= std::abs (temperature - average_temperature);
+
+ max_residual = std::max (residual, max_residual);
+ max_velocity = std::max (std::sqrt (u * u), max_velocity);
+ max_density = std::max (density, max_density);
+ max_specific_heat = std::max (c_P, max_specific_heat);
+ }
+ }
+
+ template<int dim>
+ double Simulator<dim>::
+ compute_viscosity (internal::Assembly::Scratch::TemperatureSystem<dim>& scratch,
+ const double global_max_velocity,
+ const double global_field_variation,
+ const double average_field,
+ const double global_entropy_variation,
+ const double cell_diameter) const
+ {
+ if (std::min (std::min (std::abs (global_max_velocity), std::abs (global_field_variation)),
+ std::abs (global_entropy_variation)) < 1e-50)
+ return 5e-3 * cell_diameter;
+
+ double max_residual = 0.0;
+ double max_velocity = 0.0;
+ double max_density = 0.0;
+ double max_specific_heat = 0.0;
+
+ compute_temperature_system_residual (scratch, average_field, max_residual, max_velocity,
+ max_density, max_specific_heat);
+
+ const double max_viscosity = parameters.stabilization_beta * max_density * max_specific_heat
+ * max_velocity * cell_diameter;
+
+ if (timestep_number <= 1)
+ // we don't have sensible timesteps during the first two iterations
+ return max_viscosity;
+
+ else
+ {
+ Assert (old_time_step > 0, ExcInternalError ());
+
+ double entropy_viscosity;
+
+ if (parameters.stabilization_alpha == 2)
+ entropy_viscosity = parameters.stabilization_c_R * cell_diameter * cell_diameter
+ * max_residual / global_entropy_variation;
+
+ else
+ entropy_viscosity = parameters.stabilization_c_R * cell_diameter * global_Omega_diameter
+ * max_velocity * max_residual
+ / (global_max_velocity * global_field_variation);
+
+ return std::min (max_viscosity, entropy_viscosity);
+ }
+ }
+
template <int dim>
void Simulator<dim>::
local_assemble_temperature_system (const std::pair<double,double> global_field_range,
@@ -1164,6 +1298,23 @@
scratch.old_temperature_grads);
scratch.fe_values_advection.get_function_gradients (old_old_solution_temperature,
scratch.old_old_temperature_grads);
+ scratch.fe_values_advection.get_function_laplacians (old_solution_temperature,
+ scratch.old_temperature_laplacians);
+ scratch.fe_values_advection.get_function_laplacians (old_old_solution_temperature,
+ scratch.old_old_temperature_laplacians);
+
+ for (unsigned int c = 0; c < parameters.n_compositional_fields; ++c)
+ {
+ scratch.fe_values_advection.get_function_values (old_solution_composition.block (c),
+ scratch.old_composition_values[c]);
+ scratch.fe_values_advection.get_function_values (old_old_solution_composition.block (c),
+ scratch.old_old_composition_values[c]);
+ scratch.fe_values_advection.get_function_gradients (old_solution_composition.block (c),
+ scratch.old_composition_grads[c]);
+ scratch.fe_values_advection.get_function_gradients (old_old_solution_composition.block (c),
+ scratch.old_old_composition_grads[c]);
+ }
+
scratch.fe_values_velocity.reinit (std_cxx1x::get<0> (cell.iterators));
scratch.fe_values_velocity[introspection.extractors.velocities].get_function_values (current_linearization_point_velocity,
scratch.current_velocity_values);
@@ -1171,13 +1322,30 @@
scratch.old_velocity_values);
scratch.fe_values_velocity[introspection.extractors.velocities].get_function_values (old_old_solution_velocity,
scratch.old_old_velocity_values);
+ scratch.fe_values_velocity[introspection.extractors.velocities].get_function_symmetric_gradients (current_linearization_point_velocity,
+ scratch.current_strain_rates);
compute_material_model_input_values (current_linearization_point_temperature,
current_linearization_point_composition,
scratch.fe_values_advection,
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] = 0.5 * (scratch.old_temperature_values[q] + scratch.old_old_temperature_values[q]);
+ scratch.explicit_material_model_inputs.position[q] = scratch.fe_values_advection.quadrature_point (q);
+
+ for (unsigned int c = 0; c < parameters.n_compositional_fields; ++c)
+ scratch.explicit_material_model_inputs.composition[q][c] = 0.5 * (scratch.old_composition_values[c][q] + scratch.old_old_composition_values[c][q]);
+ }
+
+ material_model->evaluate (scratch.explicit_material_model_inputs, scratch.explicit_material_model_outputs);
+
const bool axisymmetric_terms = parameters.use_axisymmetric_discretization;
+ const double nu = compute_viscosity (scratch, global_max_velocity,
+ global_field_range.second - global_field_range.first,
+ 0.5 * (global_field_range.first + global_field_range.second),
+ global_entropy_variation, std_cxx1x::get<1> (cell.iterators)->diameter ());
for (unsigned int q = 0; q < n_q_points; ++q)
{
@@ -1188,7 +1356,7 @@
}
const double c_P = scratch.material_model_outputs.specific_heat[q];
- const double lambda = scratch.material_model_outputs.thermal_conductivities[q];
+ const double lambda = scratch.material_model_outputs.thermal_conductivities[q] + nu;
const double current_value = scratch.current_temperature_values[q];
const double rho = scratch.material_model_outputs.densities[q];
const Tensor<1, dim> current_u = scratch.current_velocity_values[q];
More information about the CIG-COMMITS
mailing list