Skip to content

Commit

Permalink
Move calculation viscous dissipation into fill_elastic_outputs
Browse files Browse the repository at this point in the history
  • Loading branch information
@anne-glerum
anne-glerum committed May 14, 2024
1 parent cae39e6 commit e63b401
Show file tree
Hide file tree
Showing 4 changed files with 23 additions and 73 deletions.
13 changes: 2 additions & 11 deletions include/aspect/material_model/rheology/elasticity.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,8 @@ namespace aspect
* Given the stress of the previous time step in the material model inputs @p in,
* the elastic shear moduli @p average_elastic_shear_moduli a each point,
* and the (viscous) viscosities given in the material model outputs object @p out,
* fill a material model outputs objects with the elastic force terms.
* fill a material model outputs objects with the elastic force terms, viscoelastic
* strain rate and viscous dissipation.
*/
void
fill_elastic_outputs (const MaterialModel::MaterialModelInputs<dim> &in,
Expand Down Expand Up @@ -117,16 +118,6 @@ namespace aspect
const std::vector<double> &average_elastic_shear_moduli,
MaterialModel::MaterialModelOutputs<dim> &out) const;

/**
* A function that fills the elastic output in the
* MaterialModelOutputs object that is handed over, if it exists.
* Does nothing otherwise.
*/
void fill_viscous_dissipation(const unsigned int point_index,
const MaterialModel::MaterialModelInputs<dim> &in,
const double average_elastic_shear_modulus,
MaterialModel::MaterialModelOutputs<dim> &out) const;

/**
* Return the values of the elastic shear moduli for each composition used in the
* rheology model.
Expand Down
64 changes: 18 additions & 46 deletions source/material_model/rheology/elasticity.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -378,6 +378,8 @@ namespace aspect

for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
{
const SymmetricTensor<2, dim> deviatoric_strain_rate = in.strain_rate[i];

// Get stress from timestep $t$ rotated and advected into the current
// timestep $t+\Delta t_c$ from the compositional fields.
// This function is evaluated before the assembly of the Stokes equations
Expand Down Expand Up @@ -423,7 +425,22 @@ namespace aspect
if ((nonlinear_solver == Parameters<dim>::NonlinearSolver::iterated_Advection_and_Newton_Stokes) ||
(nonlinear_solver == Parameters<dim>::NonlinearSolver::single_Advection_iterated_Newton_Stokes))
elastic_out->viscoelastic_strain_rate[i] = calculate_viscoelastic_strain_rate(
in.strain_rate[i], stress_0_advected, stress_old, effective_creep_viscosity, average_elastic_shear_moduli[i]);
deviatoric_strain_rate, stress_0_advected, stress_old, effective_creep_viscosity, average_elastic_shear_moduli[i]);

// Apply the stress update to get the total stress of timestep t.
const SymmetricTensor<2, dim> stress = 2. * effective_creep_viscosity * deviatoric_strain_rate + viscosity_ratio * stress_0_advected +
(1. - timestep_ratio) * (1. - viscosity_ratio) * 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();

// Assume incompressibility. If compressible,
// visco_plastic_strain_rate = visco_plastic_strain_rate -
// 1. / 3. * trace(visco_plastic_strain_rate) * unit_symmetric_tensor<dim>();
const SymmetricTensor<2, dim> visco_plastic_strain_rate = deviatoric_strain_rate - ((stress - stress_0_advected) / (2. * dtc * average_elastic_shear_moduli[i]));

elastic_out->viscous_dissipation[i] = stress * visco_plastic_strain_rate;
}

}
Expand Down Expand Up @@ -650,52 +667,7 @@ namespace aspect
}
}

template <int dim>
void
Elasticity<dim>::
fill_viscous_dissipation(const unsigned int q,
const MaterialModel::MaterialModelInputs<dim> &in,
const double average_elastic_shear_modulus,
MaterialModel::MaterialModelOutputs<dim> &out) const
{
MaterialModel::ElasticOutputs<dim> *elastic_out = out.template get_additional_output<MaterialModel::ElasticOutputs<dim>>();

if (elastic_out != nullptr)
{
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.
for (unsigned int j = 0; j < SymmetricTensor<2, dim>::n_independent_components; ++j)
stress_0[SymmetricTensor<2, dim>::unrolled_to_component_indices(j)] = in.composition[q][stress_field_indices[j]];

// 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][stress_field_indices[SymmetricTensor<2, dim>::n_independent_components + j]];

// Retrieve the timestep ratio dtc/dte and the elastic viscosity,
// only two material models support elasticity.
const double timestep_ratio = calculate_timestep_ratio();
const double elastic_viscosity = calculate_elastic_viscosity(average_elastic_shear_modulus);

// 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.
const SymmetricTensor<2, dim> stress = 2. * out.viscosities[q] * deviatoric_strain_rate + out.viscosities[q] / elastic_viscosity * stress_0 +
(1. - timestep_ratio) * (1. - out.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();

// Assume incompressibility. If compressible,
// visco_plastic_strain_rate = visco_plastic_strain_rate -
// 1. / 3. * trace(visco_plastic_strain_rate) * unit_symmetric_tensor<dim>();
const SymmetricTensor<2, dim> visco_plastic_strain_rate = deviatoric_strain_rate - ((stress - stress_0) / (2. * dtc * average_elastic_shear_modulus));

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

template <int dim>
double
Expand Down
10 changes: 2 additions & 8 deletions source/material_model/visco_plastic.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -296,13 +296,6 @@ namespace aspect
elastic_additional_out->elastic_shear_moduli[i] = average_elastic_shear_moduli[i];

}

// Fill the material properties that are part of the elastic outputs other than the force term.
// At the moment, that means compute the viscous dissipation based on the viscoplastic strain rate.
if (ElasticOutputs<dim> *elastic_out = out.template get_additional_output<ElasticOutputs<dim>>())
{
rheology->elastic_rheology.fill_viscous_dissipation(i, in, average_elastic_shear_moduli[i], out);
}
}
}

Expand All @@ -311,7 +304,8 @@ namespace aspect

if (this->get_parameters().enable_elasticity)
{
// Fill the force outputs with the body force term for the RHS.
// Fill the elastic outputs with the body force term for the RHS, the viscoelastic strain rate
// and the viscous dissipation.
rheology->elastic_rheology.fill_elastic_outputs(in, average_elastic_shear_moduli, out);
// Fill the reaction terms that account for the rotation of the stresses.
rheology->elastic_rheology.fill_reaction_outputs(in, average_elastic_shear_moduli, out);
Expand Down
9 changes: 1 addition & 8 deletions source/material_model/viscoelastic.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,16 +83,9 @@ namespace aspect
// calculate_viscoelastic_viscosity function.
out.viscosities[i] = elastic_rheology.calculate_viscoelastic_viscosity(average_viscosity,
average_elastic_shear_moduli[i]);

// Fill the material properties that are part of the elastic additional outputs
if (ElasticAdditionalOutputs<dim> *elastic_out = out.template get_additional_output<ElasticAdditionalOutputs<dim>>())
{
elastic_out->elastic_shear_moduli[i] = average_elastic_shear_moduli[i];
elastic_rheology.fill_viscous_dissipation(i, in, average_elastic_shear_moduli[i], out);
}
}

// Fill the body force term.
// Fill the body force term, viscoelastic strain rate and viscous dissipation.
elastic_rheology.fill_elastic_outputs(in, average_elastic_shear_moduli, out);
// Fill the reaction terms to apply the rotation of the stresses into the current timestep.
elastic_rheology.fill_reaction_outputs(in, average_elastic_shear_moduli, out);
Expand Down

0 comments on commit e63b401

Please sign in to comment.