Skip to content

Commit 5c249f4

Browse files
bangerthbangerth
authored
Merge pull request #6728 from gassmoeller/unify_anisotropic_viscosity_code
Unify anisotropic viscosity Stokes assemblers / Fix anisotropic viscosity cookbook
2 parents f50254f + 51c647a commit 5c249f4

File tree

14 files changed

+937
-1040
lines changed

14 files changed

+937
-1040
lines changed

‎cookbooks/anisotropic_viscosity/AV_Rayleigh_Taylor.prm‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -62,7 +62,7 @@ subsection Initial composition model
6262

6363
subsection Function
6464
set Variable names = x,z
65-
set Function constants = pi=3.1415926;
65+
set Function constants = pi=3.1415926
6666
set Function expression = 0.5*(1+tanh((z-0.85-0.02*cos(2*pi*x/2))/0.02)); \
6767
(0.5*(1+tanh((z-0.85-0.02*cos(2*pi*x/2))/0.02))) > 0.8 ? sin(45*pi/180) : 0.0; \
6868
(0.5*(1+tanh((z-0.85-0.02*cos(2*pi*x/2))/0.02))) > 0.8 ? cos(45*pi/180) : 0.0;

‎cookbooks/anisotropic_viscosity/av_material.cc‎

Lines changed: 2 additions & 436 deletions
Large diffs are not rendered by default.

‎include/aspect/heating_model/shear_heating_anisotropic_viscosity.h‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
/*
2+
Copyright (C) 2025 by the authors of the ASPECT code.
3+
4+
This file is part of ASPECT.
5+
6+
ASPECT is free software; you can redistribute it and/or modify
7+
it under the terms of the GNU General Public License as published by
8+
the Free Software Foundation; either version 2, or (at your option)
9+
any later version.
10+
11+
ASPECT is distributed in the hope that it will be useful,
12+
but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
GNU General Public License for more details.
15+
16+
You should have received a copy of the GNU General Public License
17+
along with ASPECT; see the file LICENSE. If not see
18+
<http://www.gnu.org/licenses/>.
19+
*/
20+
21+
#ifndef _aspect_heating_model_shear_heating_anisotropic_viscosity_h
22+
#define _aspect_heating_model_shear_heating_anisotropic_viscosity_h
23+
24+
#include <aspect/heating_model/interface.h>
25+
26+
#include <aspect/simulator_access.h>
27+
28+
namespace aspect
29+
{
30+
namespace HeatingModel
31+
{
32+
/**
33+
* A class that implements a standard model for shear heating extended for an
34+
* anisotropic viscosity tensor. If the material model provides a stress-strain
35+
* director tensor, then the strain-rate is multiplied with this
36+
* tensor to compute the stress that is used when computing the shear heating.
37+
*
38+
* @ingroup HeatingModels
39+
*/
40+
template <int dim>
41+
class ShearHeatingAnisotropicViscosity : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
42+
{
43+
public:
44+
/**
45+
* Compute the heating model outputs for this class.
46+
*/
47+
void
48+
evaluate (const MaterialModel::MaterialModelInputs<dim> &material_model_inputs,
49+
const MaterialModel::MaterialModelOutputs<dim> &material_model_outputs,
50+
HeatingModel::HeatingModelOutputs &heating_model_outputs) const override;
51+
52+
/**
53+
* Allow the heating model to attach additional material model outputs.
54+
*/
55+
void
56+
create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &material_model_outputs) const override;
57+
};
58+
}
59+
}
60+
61+
#endif

‎include/aspect/material_model/additional_outputs/anisotropic_viscosity.h‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
/*
2+
Copyright (C) 2025 - by the authors of the ASPECT code.
3+
4+
This file is part of ASPECT.
5+
6+
ASPECT is free software; you can redistribute it and/or modify
7+
it under the terms of the GNU General Public License as published by
8+
the Free Software Foundation; either version 2, or (at your option)
9+
any later version.
10+
11+
ASPECT is distributed in the hope that it will be useful,
12+
but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
GNU General Public License for more details.
15+
16+
You should have received a copy of the GNU General Public License
17+
along with ASPECT; see the file doc/COPYING. If not see
18+
<http://www.gnu.org/licenses/>.
19+
*/
20+
21+
#ifndef _aspect_material_model_additional_outputs_anisotropic_viscosity_h
22+
#define _aspect_material_model_additional_outputs_anisotropic_viscosity_h
23+
24+
25+
#include <aspect/material_model/interface.h>
26+
27+
namespace aspect
28+
{
29+
namespace MaterialModel
30+
{
31+
/**
32+
* Additional output fields for anisotropic viscosities to be added to
33+
* the MaterialModel::MaterialModelOutputs structure and filled in the
34+
* MaterialModel::Interface::evaluate() function.
35+
*/
36+
template <int dim>
37+
class AnisotropicViscosity : public NamedAdditionalMaterialOutputs<dim>
38+
{
39+
public:
40+
AnisotropicViscosity(const unsigned int n_points);
41+
42+
std::vector<double>
43+
get_nth_output(const unsigned int idx) const override;
44+
45+
/**
46+
* Stress-strain "director" tensors at the given positions. This
47+
* variable is used to implement anisotropic viscosity.
48+
*
49+
* @note The strain rate term in equation (1) of the manual will be
50+
* multiplied by this tensor *and* the viscosity scalar ($\eta$), as
51+
* described in the manual section titled "Constitutive laws". This
52+
* variable is assigned the rank-four identity tensor by default.
53+
* This leaves the isotropic constitutive law unchanged if the material
54+
* model does not explicitly assign a value.
55+
*/
56+
std::vector<SymmetricTensor<4,dim>> stress_strain_directors;
57+
};
58+
}
59+
}
60+
61+
#endif

‎include/aspect/simulator/assemblers/stokes.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -78,7 +78,7 @@ namespace aspect
7878
};
7979

8080
/**
81-
* This class assembles the term that arises in the viscosity term of Stokes matrix for
81+
* This class assembles the term that arises in the viscosity term of the Stokes matrix for
8282
* compressible models, because the divergence of the velocity is not longer zero.
8383
*/
8484
template <int dim>

‎include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h‎

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
/*
2+
Copyright (C) 2025 - by the authors of the ASPECT code.
3+
4+
This file is part of ASPECT.
5+
6+
ASPECT is free software; you can redistribute it and/or modify
7+
it under the terms of the GNU General Public License as published by
8+
the Free Software Foundation; either version 2, or (at your option)
9+
any later version.
10+
11+
ASPECT is distributed in the hope that it will be useful,
12+
but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
GNU General Public License for more details.
15+
16+
You should have received a copy of the GNU General Public License
17+
along with ASPECT; see the file doc/COPYING. If not see
18+
<http://www.gnu.org/licenses/>.
19+
*/
20+
21+
#ifndef _aspect_simulator_assemblers_stokes_anisotropic_viscosity_h
22+
#define _aspect_simulator_assemblers_stokes_anisotropic_viscosity_h
23+
24+
25+
#include <aspect/simulator/assemblers/interface.h>
26+
#include <aspect/simulator_access.h>
27+
28+
namespace aspect
29+
{
30+
namespace Assemblers
31+
{
32+
/**
33+
* A class containing the functions to assemble the Stokes preconditioner for the
34+
* case of anisotropic viscosities.
35+
*/
36+
template <int dim>
37+
class StokesPreconditionerAnisotropicViscosity : public Assemblers::Interface<dim>,
38+
public SimulatorAccess<dim>
39+
{
40+
public:
41+
void
42+
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch,
43+
internal::Assembly::CopyData::CopyDataBase<dim> &data) const override;
44+
45+
/**
46+
* Create AnisotropicViscosities.
47+
*/
48+
void
49+
create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const override;
50+
};
51+
52+
/**
53+
* A class containing the functions to assemble the compressible adjustment
54+
* to the Stokes preconditioner for the case of anisotropic viscosities.
55+
*/
56+
template <int dim>
57+
class StokesCompressiblePreconditionerAnisotropicViscosity : public Assemblers::Interface<dim>,
58+
public SimulatorAccess<dim>
59+
{
60+
public:
61+
void
62+
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch,
63+
internal::Assembly::CopyData::CopyDataBase<dim> &data) const override;
64+
};
65+
66+
/**
67+
* This class assembles the terms for the matrix and right-hand-side of the incompressible
68+
* Stokes equation for the case of anisotropic viscosities.
69+
*/
70+
template <int dim>
71+
class StokesIncompressibleTermsAnisotropicViscosity : public Assemblers::Interface<dim>,
72+
public SimulatorAccess<dim>
73+
{
74+
public:
75+
void
76+
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch,
77+
internal::Assembly::CopyData::CopyDataBase<dim> &data) const override;
78+
79+
/**
80+
* Create AdditionalMaterialOutputsStokesRHS if we need to do so.
81+
*/
82+
void
83+
create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const override;
84+
};
85+
86+
/**
87+
* This class assembles the term that arises in the viscosity term of the Stokes matrix for
88+
* compressible models for the case of anisotropic viscosities, because the divergence
89+
* of the velocity is not longer zero.
90+
*/
91+
template <int dim>
92+
class StokesCompressibleStrainRateViscosityTermAnisotropicViscosity : public Assemblers::Interface<dim>,
93+
public SimulatorAccess<dim>
94+
{
95+
public:
96+
void
97+
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch,
98+
internal::Assembly::CopyData::CopyDataBase<dim> &data) const override;
99+
};
100+
}
101+
}
102+
103+
104+
#endif

‎source/heating_model/shear_heating_anisotropic_viscosity.cc‎

Lines changed: 123 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,123 @@
1+
/*
2+
Copyright (C) 2015 - 2023 by the authors of the ASPECT code.
3+
4+
This file is part of ASPECT.
5+
6+
ASPECT is free software; you can redistribute it and/or modify
7+
it under the terms of the GNU General Public License as published by
8+
the Free Software Foundation; either version 2, or (at your option)
9+
any later version.
10+
11+
ASPECT is distributed in the hope that it will be useful,
12+
but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
GNU General Public License for more details.
15+
16+
You should have received a copy of the GNU General Public License
17+
along with ASPECT; see the file LICENSE. If not see
18+
<http://www.gnu.org/licenses/>.
19+
*/
20+
21+
#include <aspect/heating_model/shear_heating_anisotropic_viscosity.h>
22+
23+
#include <aspect/material_model/interface.h>
24+
#include <aspect/material_model/additional_outputs/anisotropic_viscosity.h>
25+
#include <aspect/heating_model/shear_heating.h>
26+
27+
#include <deal.II/base/symmetric_tensor.h>
28+
#include <deal.II/base/signaling_nan.h>
29+
30+
namespace aspect
31+
{
32+
namespace HeatingModel
33+
{
34+
template <int dim>
35+
void
36+
ShearHeatingAnisotropicViscosity<dim>::
37+
evaluate (const MaterialModel::MaterialModelInputs<dim> &material_model_inputs,
38+
const MaterialModel::MaterialModelOutputs<dim> &material_model_outputs,
39+
HeatingModel::HeatingModelOutputs &heating_model_outputs) const
40+
{
41+
Assert(heating_model_outputs.heating_source_terms.size() == material_model_inputs.position.size(),
42+
ExcMessage ("Heating outputs need to have the same number of entries as the material model inputs."));
43+
44+
// Some material models provide dislocation viscosities and boundary area work fractions
45+
// as additional material outputs. If they are attached, use them.
46+
const std::shared_ptr<const ShearHeatingOutputs<dim>> shear_heating_out =
47+
material_model_outputs.template get_additional_output_object<ShearHeatingOutputs<dim>>();
48+
49+
const std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity =
50+
material_model_outputs.template get_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>();
51+
52+
Assert(anisotropic_viscosity != nullptr,
53+
ExcMessage("This heating plugin should only be used with material models that provide "
54+
"an anisotropic viscosity tensor, but none was provided."));
55+
56+
for (unsigned int q=0; q<heating_model_outputs.heating_source_terms.size(); ++q)
57+
{
58+
// If there is an anisotropic viscosity, use it to compute the correct stress
59+
const SymmetricTensor<2,dim> &directed_strain_rate = anisotropic_viscosity->stress_strain_directors[q]
60+
* material_model_inputs.strain_rate[q];
61+
62+
const SymmetricTensor<2,dim> stress =
63+
2 * material_model_outputs.viscosities[q] *
64+
(this->get_material_model().is_compressible()
65+
?
66+
directed_strain_rate - 1./3. * trace(directed_strain_rate) * unit_symmetric_tensor<dim>()
67+
:
68+
directed_strain_rate);
69+
70+
const SymmetricTensor<2,dim> deviatoric_strain_rate =
71+
(this->get_material_model().is_compressible()
72+
?
73+
material_model_inputs.strain_rate[q]
74+
- 1./3. * trace(material_model_inputs.strain_rate[q]) * unit_symmetric_tensor<dim>()
75+
:
76+
material_model_inputs.strain_rate[q]);
77+
78+
heating_model_outputs.heating_source_terms[q] = stress * deviatoric_strain_rate;
79+
80+
// If shear heating work fractions are provided, reduce the
81+
// overall heating by this amount (which is assumed to be converted into other forms of energy)
82+
if (shear_heating_out != nullptr)
83+
heating_model_outputs.heating_source_terms[q] *= shear_heating_out->shear_heating_work_fractions[q];
84+
85+
heating_model_outputs.lhs_latent_heat_terms[q] = 0.0;
86+
}
87+
}
88+
89+
90+
91+
template <int dim>
92+
void
93+
ShearHeatingAnisotropicViscosity<dim>::
94+
create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &material_model_outputs) const
95+
{
96+
const unsigned int n_points = material_model_outputs.viscosities.size();
97+
98+
if (material_model_outputs.template has_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>() == false)
99+
{
100+
material_model_outputs.additional_outputs.push_back(
101+
std::make_unique<MaterialModel::AnisotropicViscosity<dim>> (n_points));
102+
}
103+
104+
this->get_material_model().create_additional_named_outputs(material_model_outputs);
105+
}
106+
}
107+
}
108+
109+
110+
111+
// explicit instantiations
112+
namespace aspect
113+
{
114+
namespace HeatingModel
115+
{
116+
ASPECT_REGISTER_HEATING_MODEL(ShearHeatingAnisotropicViscosity,
117+
"anisotropic shear heating",
118+
"Implementation of a standard model for shear heating extended for an "
119+
"anisotropic viscosity tensor. If the material model provides a stress-"
120+
"strain director tensor, then the strain-rate is multiplied with this "
121+
"tensor to compute the stress that is used when computing the shear heating.")
122+
}
123+
}

0 commit comments

Comments
 (0)