Skip to content

Commit

Permalink
Use ElasticOutputs in shear_heating
Browse files Browse the repository at this point in the history
  • Loading branch information
@anne-glerum
anne-glerum committed Feb 18, 2024
1 parent 718794b commit 83bcac5
Show file tree
Hide file tree
Showing 2 changed files with 4 additions and 66 deletions.
66 changes: 2 additions & 64 deletions source/heating_model/shear_heating.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,75 +80,13 @@ namespace aspect

// If elasticity is used, the stress should include the elastic stresses
// and only the visco-plastic (non-elastic) strain rate should contribute.
// Directly ask the material model for the dissipation.
if (this->get_parameters().enable_elasticity == true)
{
// Visco-elastic stresses are stored on the fields
SymmetricTensor<2, dim> stress_0, stress_old;
stress_0[0][0] = material_model_inputs.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")];
stress_0[1][1] = material_model_inputs.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")];
stress_0[0][1] = material_model_inputs.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")];

stress_old[0][0] = material_model_inputs.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx_old")];
stress_old[1][1] = material_model_inputs.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy_old")];
stress_old[0][1] = material_model_inputs.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy_old")];

if (dim == 3)
{
stress_0[2][2] = material_model_inputs.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")];
stress_0[0][2] = material_model_inputs.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")];
stress_0[1][2] = material_model_inputs.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")];

stress_old[2][2] = material_model_inputs.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz_old")];
stress_old[0][2] = material_model_inputs.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz_old")];
stress_old[1][2] = material_model_inputs.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz_old")];
}

const MaterialModel::ElasticAdditionalOutputs<dim> *elastic_out = material_model_outputs.template get_additional_output<MaterialModel::ElasticAdditionalOutputs<dim>>();
AssertThrow(elastic_out != nullptr, ExcMessage("ElasticAdditionalOutputs are requested, but have not been computed."));

const double shear_modulus = elastic_out->elastic_shear_moduli[q];

// Retrieve the timestep ratio dtc/dte and the elastic timestep and viscosity,
// only two material models support elasticity.
double timestep_ratio = 0.;
double elastic_timestep = 0.;
double elastic_viscosity = 0.;
if (Plugins::plugin_type_matches<MaterialModel::ViscoPlastic<dim>>(this->get_material_model()))
{
const MaterialModel::ViscoPlastic<dim> &vp = Plugins::get_plugin_as_type<const MaterialModel::ViscoPlastic<dim>>(this->get_material_model());
elastic_viscosity = vp.get_elastic_viscosity(shear_modulus);
elastic_timestep = vp.get_elastic_timestep();
timestep_ratio = vp.get_timestep_ratio();
}
else if (Plugins::plugin_type_matches<MaterialModel::Viscoelastic<dim>>(this->get_material_model()))
{
const MaterialModel::Viscoelastic<dim> &ve = Plugins::get_plugin_as_type<const MaterialModel::Viscoelastic<dim>>(this->get_material_model());
elastic_viscosity = ve.get_elastic_viscosity(shear_modulus);
elastic_timestep = ve.get_elastic_timestep();
timestep_ratio = ve.get_timestep_ratio();
}
else
AssertThrow(false, ExcMessage("Shear heating cannot be used with elasticity for material models other than ViscoPlastic and Viscoelastic."));

// Apply the stress update to get the total stress of timestep t.
// Both the viscous and the elastic viscosity have been scaled with the timestep ratio.
stress = 2. * material_model_outputs.viscosities[q] * deviatoric_strain_rate + material_model_outputs.viscosities[q] / elastic_viscosity * stress_0 +
(1. - timestep_ratio) * (1. - material_model_outputs.viscosities[q] / elastic_viscosity) * stress_old;

// Obtain the computational timestep by multiplying the ratio between the computational
// and elastic timestep $\frac{\Delta t_c}{\Delta t_{el}}$ with the elastic timestep.
const double dtc = timestep_ratio * elastic_timestep;

SymmetricTensor<2, dim> visco_plastic_strain_rate = material_model_inputs.strain_rate[q] - ((stress - stress_0) / (2. * dtc * shear_modulus));

if (this->get_material_model().is_compressible())
visco_plastic_strain_rate = visco_plastic_strain_rate -
1. / 3. * trace(visco_plastic_strain_rate) * unit_symmetric_tensor<dim>();

// TODO only use elastic_out
std::cout << "Old and new viscous dissipation: " << stress*visco_plastic_strain_rate << " " << elastic_out->viscous_dissipation[q] << std::endl;

heating_model_outputs.heating_source_terms[q] = stress * visco_plastic_strain_rate;
heating_model_outputs.heating_source_terms[q] = elastic_out->viscous_dissipation[q];
}

// If shear heating work fractions are provided, reduce the
Expand Down
4 changes: 2 additions & 2 deletions source/material_model/rheology/elasticity.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -584,7 +584,7 @@ namespace aspect

if (elastic_out != nullptr)
{
const SymmetricTensor<2,dim> deviatoric_strain_rate = in.strain_rate[q];
const SymmetricTensor<2, dim> deviatoric_strain_rate = in.strain_rate[q];

SymmetricTensor<2, dim> stress_0;
// The current stress is stored on the first n_independent_components fields.
Expand All @@ -594,7 +594,7 @@ namespace aspect
// The old stress is stored on the second set of n_independent_components fields.
SymmetricTensor<2, dim> stress_old;
for (unsigned int j = 0; j < SymmetricTensor<2, dim>::n_independent_components; ++j)
stress_old[SymmetricTensor<2, dim>::unrolled_to_component_indices(j)] = in.composition[q][SymmetricTensor<2, dim>::n_independent_components+j];
stress_old[SymmetricTensor<2, dim>::unrolled_to_component_indices(j)] = in.composition[q][SymmetricTensor<2, dim>::n_independent_components + j];

// Retrieve the timestep ratio dtc/dte and the elastic viscosity,
// only two material models support elasticity.
Expand Down

0 comments on commit 83bcac5

Please sign in to comment.