diff --git a/cookbooks/anisotropic_viscosity/AV_Rayleigh_Taylor.prm b/cookbooks/anisotropic_viscosity/AV_Rayleigh_Taylor.prm index 57ca6cc1027..57a1a3d5b28 100644 --- a/cookbooks/anisotropic_viscosity/AV_Rayleigh_Taylor.prm +++ b/cookbooks/anisotropic_viscosity/AV_Rayleigh_Taylor.prm @@ -62,7 +62,7 @@ subsection Initial composition model subsection Function set Variable names = x,z - set Function constants = pi=3.1415926; + set Function constants = pi=3.1415926 set Function expression = 0.5*(1+tanh((z-0.85-0.02*cos(2*pi*x/2))/0.02)); \ (0.5*(1+tanh((z-0.85-0.02*cos(2*pi*x/2))/0.02))) > 0.8 ? sin(45*pi/180) : 0.0; \ (0.5*(1+tanh((z-0.85-0.02*cos(2*pi*x/2))/0.02))) > 0.8 ? cos(45*pi/180) : 0.0; diff --git a/cookbooks/anisotropic_viscosity/av_material.cc b/cookbooks/anisotropic_viscosity/av_material.cc index 7bc25da38e3..1770bdd3fba 100644 --- a/cookbooks/anisotropic_viscosity/av_material.cc +++ b/cookbooks/anisotropic_viscosity/av_material.cc @@ -21,7 +21,7 @@ #include #include #include -#include +#include #include #include #include @@ -48,8 +48,6 @@ #include #include -#include -#include #include #include #include @@ -60,414 +58,6 @@ namespace aspect { - namespace MaterialModel - { - /** - * Additional output fields for anisotropic viscosities to be added to - * the MaterialModel::MaterialModelOutputs structure and filled in the - * MaterialModel::Interface::evaluate() function. - */ - template - class AnisotropicViscosity : public NamedAdditionalMaterialOutputs - { - public: - AnisotropicViscosity(const unsigned int n_points); - - std::vector get_nth_output(const unsigned int idx) const override; - - /** - * Stress-strain "director" tensors at the given positions. This - * variable is used to implement anisotropic viscosity. - * - * @note The strain rate term in equation (1) of the manual will be - * multiplied by this tensor *and* the viscosity scalar ($\eta$), as - * described in the manual section titled "Constitutive laws". This - * variable is assigned the rank-four identity tensor by default. - * This leaves the isotropic constitutive law unchanged if the material - * model does not explicitly assign a value. - */ - std::vector> stress_strain_directors; - }; - - namespace - { - - - - template - std::vector make_AnisotropicViscosity_additional_outputs_names() - { - std::vector names; - - for (unsigned int i = 0; i < Tensor<4,dim>::n_independent_components ; ++i) - { - TableIndices<4> indices(Tensor<4,dim>::unrolled_to_component_indices(i)); - names.push_back("anisotropic_viscosity"+std::to_string(indices[0])+std::to_string(indices[1])+std::to_string(indices[2])+std::to_string(indices[3])); - } - return names; - } - } - - - - template - AnisotropicViscosity::AnisotropicViscosity (const unsigned int n_points) - : - NamedAdditionalMaterialOutputs(make_AnisotropicViscosity_additional_outputs_names()), - stress_strain_directors(n_points, dealii::identity_tensor ()) - {} - - - - template - std::vector - AnisotropicViscosity::get_nth_output(const unsigned int idx) const - { - std::vector output(stress_strain_directors.size()); - for (unsigned int i = 0; i < stress_strain_directors.size() ; ++i) - { - output[i]= stress_strain_directors[i][Tensor<4,dim>::unrolled_to_component_indices(idx)]; - } - return output; - } - } -} - -namespace aspect -{ - namespace Assemblers - { - /** - * A class containing the functions to assemble the Stokes preconditioner. - */ - template - class StokesPreconditionerAnisotropicViscosity : public Assemblers::Interface, - public SimulatorAccess - { - public: - void - execute(internal::Assembly::Scratch::ScratchBase &scratch, - internal::Assembly::CopyData::CopyDataBase &data) const override; - - /** - * Create AnisotropicViscosities. - */ - void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const override; - }; - - /** - * This class assembles the terms for the matrix and right-hand-side of the incompressible - * Stokes equation for the current cell. - */ - template - class StokesIncompressibleTermsAnisotropicViscosity : public Assemblers::Interface, - public SimulatorAccess - { - public: - void - execute(internal::Assembly::Scratch::ScratchBase &scratch, - internal::Assembly::CopyData::CopyDataBase &data) const override; - - /** - * Create AdditionalMaterialOutputsStokesRHS if we need to do so. - */ - void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const override; - }; - - - - template - void - StokesPreconditionerAnisotropicViscosity:: - execute (internal::Assembly::Scratch::ScratchBase &scratch_base, - internal::Assembly::CopyData::CopyDataBase &data_base) const - { - internal::Assembly::Scratch::StokesPreconditioner &scratch = dynamic_cast&> (scratch_base); - internal::Assembly::CopyData::StokesPreconditioner &data = dynamic_cast&> (data_base); - - const std::shared_ptr> anisotropic_viscosity - = scratch.material_model_outputs.template get_additional_output_object>(); - - const Introspection &introspection = this->introspection(); - const FiniteElement &fe = this->get_fe(); - const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); - const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; - const double pressure_scaling = this->get_pressure_scaling(); - - // First loop over all dofs and find those that are in the Stokes system - // save the component (pressure and dim velocities) each belongs to. - for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/) - { - if (introspection.is_stokes_component(fe.system_to_component_index(i).first)) - { - scratch.dof_component_indices[i_stokes] = fe.system_to_component_index(i).first; - ++i_stokes; - } - ++i; - } - - // Loop over all quadrature points and assemble their contributions to - // the preconditioner matrix - for (unsigned int q = 0; q < n_q_points; ++q) - { - for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/) - { - if (introspection.is_stokes_component(fe.system_to_component_index(i).first)) - { - scratch.grads_phi_u[i_stokes] = - scratch.finite_element_values[introspection.extractors - .velocities].symmetric_gradient(i, q); - scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection - .extractors.pressure].value(i, q); - ++i_stokes; - } - ++i; - } - - const double eta = scratch.material_model_outputs.viscosities[q]; - const double one_over_eta = 1. / eta; - const SymmetricTensor<4, dim> &stress_strain_director = anisotropic_viscosity->stress_strain_directors[q]; - const double JxW = scratch.finite_element_values.JxW(q); - - for (unsigned int i = 0; i < stokes_dofs_per_cell; ++i) - for (unsigned int j = 0; j < stokes_dofs_per_cell; ++j) - if (scratch.dof_component_indices[i] == - scratch.dof_component_indices[j]) - data.local_matrix(i, j) += (2.0 * eta * (scratch.grads_phi_u[i] - * stress_strain_director - * scratch.grads_phi_u[j]) - + one_over_eta * pressure_scaling - * pressure_scaling - * (scratch.phi_p[i] - * scratch.phi_p[j])) - * JxW; - } - } - - - - template - void - StokesPreconditionerAnisotropicViscosity:: - create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const - { - const unsigned int n_points = outputs.viscosities.size(); - - if (outputs.template has_additional_output_object>() == false) - { - outputs.additional_outputs.push_back( - std::make_unique> (n_points)); - } - } - - - - template - void - StokesIncompressibleTermsAnisotropicViscosity:: - execute (internal::Assembly::Scratch::ScratchBase &scratch_base, - internal::Assembly::CopyData::CopyDataBase &data_base) const - { - internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast&> (scratch_base); - internal::Assembly::CopyData::StokesSystem &data = dynamic_cast&> (data_base); - - const std::shared_ptr> anisotropic_viscosity - = scratch.material_model_outputs.template get_additional_output_object>(); - - const Introspection &introspection = this->introspection(); - const FiniteElement &fe = this->get_fe(); - const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); - const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; - const double pressure_scaling = this->get_pressure_scaling(); - - const std::shared_ptr> force - = scratch.material_model_outputs.template get_additional_output_object>(); - - for (unsigned int q=0; q()); - - const SymmetricTensor<4, dim> &stress_strain_director = anisotropic_viscosity->stress_strain_directors[q]; - - const Tensor<1,dim> - gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q)); - - const double density = scratch.material_model_outputs.densities[q]; - const double JxW = scratch.finite_element_values.JxW(q); - - for (unsigned int i=0; irhs_u[q] * scratch.phi_u[i] - + pressure_scaling * force->rhs_p[q] * scratch.phi_p[i]) - * JxW; - - if (scratch.rebuild_stokes_matrix) - for (unsigned int j=0; j - void - StokesIncompressibleTermsAnisotropicViscosity:: - create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const - { - const unsigned int n_points = outputs.viscosities.size(); - - if (outputs.template has_additional_output_object>() == false) - { - outputs.additional_outputs.push_back( - std::make_unique> (n_points)); - } - - if (this->get_parameters().enable_additional_stokes_rhs - && outputs.template has_additional_output_object>() == false) - { - outputs.additional_outputs.push_back( - std::make_unique> (n_points)); - } - Assert(!this->get_parameters().enable_additional_stokes_rhs - || - outputs.template get_additional_output_object>()->rhs_u.size() - == n_points, ExcInternalError()); - } - } - - namespace HeatingModel - { - template - class ShearHeatingAnisotropicViscosity : public Interface, public ::aspect::SimulatorAccess - { - public: - /** - * Compute the heating model outputs for this class. - */ - void - evaluate (const MaterialModel::MaterialModelInputs &material_model_inputs, - const MaterialModel::MaterialModelOutputs &material_model_outputs, - HeatingModel::HeatingModelOutputs &heating_model_outputs) const override; - - /** - * Allow the heating model to attach additional material model outputs. - */ - void - create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &material_model_outputs) const override; - }; - - - - template - void - ShearHeatingAnisotropicViscosity:: - evaluate (const MaterialModel::MaterialModelInputs &material_model_inputs, - const MaterialModel::MaterialModelOutputs &material_model_outputs, - HeatingModel::HeatingModelOutputs &heating_model_outputs) const - { - Assert(heating_model_outputs.heating_source_terms.size() == material_model_inputs.position.size(), - ExcMessage ("Heating outputs need to have the same number of entries as the material model inputs.")); - - // Some material models provide dislocation viscosities and boundary area work fractions - // as additional material outputs. If they are attached, use them. - const std::shared_ptr> shear_heating_out - = material_model_outputs.template get_additional_output_object>(); - - const std::shared_ptr> anisotropic_viscosity - = material_model_outputs.template get_additional_output_object>(); - - for (unsigned int q=0; q &directed_strain_rate = ((anisotropic_viscosity != nullptr) - ? - anisotropic_viscosity->stress_strain_directors[q] - * material_model_inputs.strain_rate[q] - : - material_model_inputs.strain_rate[q]); - - const SymmetricTensor<2,dim> stress = - 2 * material_model_outputs.viscosities[q] * - (this->get_material_model().is_compressible() - ? - directed_strain_rate - 1./3. * trace(directed_strain_rate) * unit_symmetric_tensor() - : - directed_strain_rate); - - const SymmetricTensor<2,dim> deviatoric_strain_rate = - (this->get_material_model().is_compressible() - ? - material_model_inputs.strain_rate[q] - - 1./3. * trace(material_model_inputs.strain_rate[q]) * unit_symmetric_tensor() - : - material_model_inputs.strain_rate[q]); - - heating_model_outputs.heating_source_terms[q] = stress * deviatoric_strain_rate; - - // If shear heating work fractions are provided, reduce the - // overall heating by this amount (which is assumed to be converted into other forms of energy) - if (shear_heating_out != nullptr) - heating_model_outputs.heating_source_terms[q] *= shear_heating_out->shear_heating_work_fractions[q]; - - heating_model_outputs.lhs_latent_heat_terms[q] = 0.0; - } - } - - - - template - void - ShearHeatingAnisotropicViscosity:: - create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &material_model_outputs) const - { - const unsigned int n_points = material_model_outputs.viscosities.size(); - - if (material_model_outputs.template has_additional_output_object>() == false) - { - material_model_outputs.additional_outputs.push_back( - std::make_unique> (n_points)); - } - - this->get_material_model().create_additional_named_outputs(material_model_outputs); - } - } - namespace MaterialModel { // The AV material model calculates an anisotropic viscosity tensor from director vectors and the normal and shear @@ -792,30 +382,6 @@ namespace aspect // explicit instantiations namespace aspect { - namespace Assemblers - { -#define INSTANTIATE(dim) \ - template class StokesPreconditioner; \ - template class StokesIncompressibleTerms; \ - template class StokesBoundaryTraction; - - ASPECT_INSTANTIATE(INSTANTIATE) - } - - namespace HeatingModel - { - ASPECT_REGISTER_HEATING_MODEL(ShearHeatingAnisotropicViscosity, - "anisotropic shear heating", - "Implementation of a standard model for shear heating. " - "Adds the term: " - "$ 2 \\eta \\left( \\varepsilon - \\frac{1}{3} \\text{tr} " - "\\varepsilon \\mathbf 1 \\right) : \\left( \\varepsilon - \\frac{1}{3} " - "\\text{tr} \\varepsilon \\mathbf 1 \\right)$ to the " - "right-hand side of the temperature equation.") - } - - - namespace MaterialModel { ASPECT_REGISTER_MATERIAL_MODEL(AV, diff --git a/include/aspect/heating_model/shear_heating_anisotropic_viscosity.h b/include/aspect/heating_model/shear_heating_anisotropic_viscosity.h new file mode 100644 index 00000000000..5614ff9334e --- /dev/null +++ b/include/aspect/heating_model/shear_heating_anisotropic_viscosity.h @@ -0,0 +1,61 @@ +/* + Copyright (C) 2025 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . + */ + +#ifndef _aspect_heating_model_shear_heating_anisotropic_viscosity_h +#define _aspect_heating_model_shear_heating_anisotropic_viscosity_h + +#include + +#include + +namespace aspect +{ + namespace HeatingModel + { + /** + * A class that implements a standard model for shear heating extended for an + * anisotropic viscosity tensor. If the material model provides a stress- + * strain director tensor, then the strain-rate is multiplied with this + * tensor to compute the stress that is used when computing the shear heating. + * + * @ingroup HeatingModels + */ + template + class ShearHeatingAnisotropicViscosity : public Interface, public ::aspect::SimulatorAccess + { + public: + /** + * Compute the heating model outputs for this class. + */ + void + evaluate (const MaterialModel::MaterialModelInputs &material_model_inputs, + const MaterialModel::MaterialModelOutputs &material_model_outputs, + HeatingModel::HeatingModelOutputs &heating_model_outputs) const override; + + /** + * Allow the heating model to attach additional material model outputs. + */ + void + create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &material_model_outputs) const override; + }; + } +} + +#endif diff --git a/include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h b/include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h new file mode 100644 index 00000000000..a5173ca85d0 --- /dev/null +++ b/include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h @@ -0,0 +1,172 @@ +/* + Copyright (C) 2025 - by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file doc/COPYING. If not see + . +*/ + +#ifndef _aspect_simulator_assemblers_stokes_anisotropic_viscosity_h +#define _aspect_simulator_assemblers_stokes_anisotropic_viscosity_h + + +#include +#include + +namespace aspect +{ + namespace MaterialModel + { + /** + * Additional output fields for anisotropic viscosities to be added to + * the MaterialModel::MaterialModelOutputs structure and filled in the + * MaterialModel::Interface::evaluate() function. + */ + template + class AnisotropicViscosity : public NamedAdditionalMaterialOutputs + { + public: + AnisotropicViscosity(const unsigned int n_points); + + std::vector + get_nth_output(const unsigned int idx) const override; + + /** + * Stress-strain "director" tensors at the given positions. This + * variable is used to implement anisotropic viscosity. + * + * @note The strain rate term in equation (1) of the manual will be + * multiplied by this tensor *and* the viscosity scalar ($\eta$), as + * described in the manual section titled "Constitutive laws". This + * variable is assigned the rank-four identity tensor by default. + * This leaves the isotropic constitutive law unchanged if the material + * model does not explicitly assign a value. + */ + std::vector> stress_strain_directors; + }; + + namespace + { + template + std::vector make_AnisotropicViscosity_additional_outputs_names() + { + std::vector names; + + for (unsigned int i = 0; i < Tensor<4,dim>::n_independent_components ; ++i) + { + TableIndices<4> indices(Tensor<4,dim>::unrolled_to_component_indices(i)); + names.push_back("anisotropic_viscosity"+std::to_string(indices[0])+std::to_string(indices[1])+std::to_string(indices[2])+std::to_string(indices[3])); + } + return names; + } + } + + + + template + AnisotropicViscosity::AnisotropicViscosity (const unsigned int n_points) + : + NamedAdditionalMaterialOutputs(make_AnisotropicViscosity_additional_outputs_names()), + stress_strain_directors(n_points, dealii::identity_tensor ()) + {} + + + + template + std::vector + AnisotropicViscosity::get_nth_output(const unsigned int idx) const + { + std::vector output(stress_strain_directors.size()); + for (unsigned int i = 0; i < stress_strain_directors.size() ; ++i) + { + output[i]= stress_strain_directors[i][Tensor<4,dim>::unrolled_to_component_indices(idx)]; + } + return output; + } + } + + namespace Assemblers + { + /** + * A class containing the functions to assemble the Stokes preconditioner. + */ + template + class StokesPreconditionerAnisotropicViscosity : public Assemblers::Interface, + public SimulatorAccess + { + public: + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const override; + + /** + * Create AnisotropicViscosities. + */ + void + create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const override; + }; + + /** + * A class containing the functions to assemble the compressible adjustment + * to the Stokes preconditioner. + */ + template + class StokesCompressiblePreconditionerAnisotropicViscosity : public Assemblers::Interface, + public SimulatorAccess + { + public: + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const override; + }; + + /** + * This class assembles the terms for the matrix and right-hand-side of the incompressible + * Stokes equation for the current cell. + */ + template + class StokesIncompressibleTermsAnisotropicViscosity : public Assemblers::Interface, + public SimulatorAccess + { + public: + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const override; + + /** + * Create AdditionalMaterialOutputsStokesRHS if we need to do so. + */ + void + create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const override; + }; + + /** + * This class assembles the term that arises in the viscosity term of Stokes matrix for + * compressible models, because the divergence of the velocity is not longer zero. + */ + template + class StokesCompressibleStrainRateViscosityTermAnisotropicViscosity : public Assemblers::Interface, + public SimulatorAccess + { + public: + void + execute(internal::Assembly::Scratch::ScratchBase &scratch, + internal::Assembly::CopyData::CopyDataBase &data) const override; + }; + } +} + + +#endif diff --git a/source/heating_model/shear_heating_anisotropic_viscosity.cc b/source/heating_model/shear_heating_anisotropic_viscosity.cc new file mode 100644 index 00000000000..fed7445582e --- /dev/null +++ b/source/heating_model/shear_heating_anisotropic_viscosity.cc @@ -0,0 +1,123 @@ +/* + Copyright (C) 2015 - 2023 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . + */ + +#include + +#include +#include +#include + +#include +#include + +namespace aspect +{ + namespace HeatingModel + { + template + void + ShearHeatingAnisotropicViscosity:: + evaluate (const MaterialModel::MaterialModelInputs &material_model_inputs, + const MaterialModel::MaterialModelOutputs &material_model_outputs, + HeatingModel::HeatingModelOutputs &heating_model_outputs) const + { + Assert(heating_model_outputs.heating_source_terms.size() == material_model_inputs.position.size(), + ExcMessage ("Heating outputs need to have the same number of entries as the material model inputs.")); + + // Some material models provide dislocation viscosities and boundary area work fractions + // as additional material outputs. If they are attached, use them. + const std::shared_ptr> shear_heating_out = + material_model_outputs.template get_additional_output_object>(); + + const std::shared_ptr> anisotropic_viscosity = + material_model_outputs.template get_additional_output_object>(); + + for (unsigned int q=0; q &directed_strain_rate = ((anisotropic_viscosity != nullptr) + ? + anisotropic_viscosity->stress_strain_directors[q] + * material_model_inputs.strain_rate[q] + : + material_model_inputs.strain_rate[q]); + + const SymmetricTensor<2,dim> stress = + 2 * material_model_outputs.viscosities[q] * + (this->get_material_model().is_compressible() + ? + directed_strain_rate - 1./3. * trace(directed_strain_rate) * unit_symmetric_tensor() + : + directed_strain_rate); + + const SymmetricTensor<2,dim> deviatoric_strain_rate = + (this->get_material_model().is_compressible() + ? + material_model_inputs.strain_rate[q] + - 1./3. * trace(material_model_inputs.strain_rate[q]) * unit_symmetric_tensor() + : + material_model_inputs.strain_rate[q]); + + heating_model_outputs.heating_source_terms[q] = stress * deviatoric_strain_rate; + + // If shear heating work fractions are provided, reduce the + // overall heating by this amount (which is assumed to be converted into other forms of energy) + if (shear_heating_out != nullptr) + heating_model_outputs.heating_source_terms[q] *= shear_heating_out->shear_heating_work_fractions[q]; + + heating_model_outputs.lhs_latent_heat_terms[q] = 0.0; + } + } + + + + template + void + ShearHeatingAnisotropicViscosity:: + create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &material_model_outputs) const + { + const unsigned int n_points = material_model_outputs.viscosities.size(); + + if (material_model_outputs.template has_additional_output_object>() == false) + { + material_model_outputs.additional_outputs.push_back( + std::make_unique> (n_points)); + } + + this->get_material_model().create_additional_named_outputs(material_model_outputs); + } + } +} + + + +// explicit instantiations +namespace aspect +{ + namespace HeatingModel + { + ASPECT_REGISTER_HEATING_MODEL(ShearHeatingAnisotropicViscosity, + "anisotropic shear heating", + "Implementation of a standard model for shear heating extended for an " + "anisotropic viscosity tensor. If the material model provides a stress-" + "strain director tensor, then the strain-rate is multiplied with this " + "tensor to compute the stress that is used when computing the shear heating.") + } +} diff --git a/source/simulator/assemblers/stokes_anisotropic_viscosity.cc b/source/simulator/assemblers/stokes_anisotropic_viscosity.cc new file mode 100644 index 00000000000..5b9cffe66eb --- /dev/null +++ b/source/simulator/assemblers/stokes_anisotropic_viscosity.cc @@ -0,0 +1,401 @@ +/* + Copyright (C) 2025 - by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + +#include + +#include + +#include + +namespace aspect +{ + namespace Assemblers + { + template + void + StokesPreconditionerAnisotropicViscosity:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const + { + internal::Assembly::Scratch::StokesPreconditioner &scratch = dynamic_cast&> (scratch_base); + internal::Assembly::CopyData::StokesPreconditioner &data = dynamic_cast&> (data_base); + + std::shared_ptr> anisotropic_viscosity = + scratch.material_model_outputs.template get_additional_output_object>(); + + const Introspection &introspection = this->introspection(); + const FiniteElement &fe = this->get_fe(); + const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); + const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; + const double pressure_scaling = this->get_pressure_scaling(); + + // First loop over all dofs and find those that are in the Stokes system + // save the component (pressure and dim velocities) each belongs to. + for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/) + { + if (introspection.is_stokes_component(fe.system_to_component_index(i).first)) + { + scratch.dof_component_indices[i_stokes] = fe.system_to_component_index(i).first; + ++i_stokes; + } + ++i; + } + + // Loop over all quadrature points and assemble their contributions to + // the preconditioner matrix + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/) + { + if (introspection.is_stokes_component(fe.system_to_component_index(i).first)) + { + scratch.grads_phi_u[i_stokes] = + scratch.finite_element_values[introspection.extractors + .velocities].symmetric_gradient(i, q); + scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection + .extractors.pressure].value(i, q); + ++i_stokes; + } + ++i; + } + + const double eta = scratch.material_model_outputs.viscosities[q]; + const double one_over_eta = 1. / eta; + + const bool use_tensor = (anisotropic_viscosity != nullptr); + + const SymmetricTensor<4, dim> &stress_strain_director = (use_tensor) + ? + anisotropic_viscosity->stress_strain_directors[q] + : + dealii::identity_tensor(); + + + + const double JxW = scratch.finite_element_values.JxW(q); + + for (unsigned int i = 0; i < stokes_dofs_per_cell; ++i) + for (unsigned int j = 0; j < stokes_dofs_per_cell; ++j) + if (scratch.dof_component_indices[i] == + scratch.dof_component_indices[j]) + data.local_matrix(i, j) += (( + use_tensor ? + 2.0 * eta * (scratch.grads_phi_u[i] + * stress_strain_director + * scratch.grads_phi_u[j]) : + 2.0 * eta * (scratch.grads_phi_u[i] + * scratch.grads_phi_u[j])) + + one_over_eta * pressure_scaling + * pressure_scaling + * (scratch.phi_p[i] + * scratch.phi_p[j])) + * JxW; + } + } + + + + template + void + StokesPreconditionerAnisotropicViscosity:: + create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const + { + const unsigned int n_points = outputs.viscosities.size(); + + if (outputs.template has_additional_output_object>() == false) + { + outputs.additional_outputs.push_back( + std::make_unique> (n_points)); + } + } + + + + template + void + StokesCompressiblePreconditionerAnisotropicViscosity:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const + { + internal::Assembly::Scratch::StokesPreconditioner &scratch = dynamic_cast&> (scratch_base); + internal::Assembly::CopyData::StokesPreconditioner &data = dynamic_cast&> (data_base); + + const std::shared_ptr> anisotropic_viscosity = + scratch.material_model_outputs.template get_additional_output_object>(); + + const Introspection &introspection = this->introspection(); + const FiniteElement &fe = this->get_fe(); + const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); + const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; + + // First loop over all dofs and find those that are in the Stokes system + // save the component (pressure and dim velocities) each belongs to. + for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/) + { + if (introspection.is_stokes_component(fe.system_to_component_index(i).first)) + { + scratch.dof_component_indices[i_stokes] = fe.system_to_component_index(i).first; + ++i_stokes; + } + ++i; + } + + // Loop over all quadrature points and assemble their contributions to + // the preconditioner matrix + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/) + { + if (introspection.is_stokes_component(fe.system_to_component_index(i).first)) + { + scratch.grads_phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].symmetric_gradient(i,q); + scratch.div_phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].divergence (i, q); + + ++i_stokes; + } + ++i; + } + + const double eta_two_thirds = scratch.material_model_outputs.viscosities[q] * 2.0 / 3.0; + + const bool use_tensor = (anisotropic_viscosity != nullptr); + + const SymmetricTensor<4, dim> &stress_strain_director = (use_tensor) + ? + anisotropic_viscosity->stress_strain_directors[q] + : + dealii::identity_tensor(); + + const double JxW = scratch.finite_element_values.JxW(q); + + for (unsigned int i = 0; i < stokes_dofs_per_cell; ++i) + for (unsigned int j = 0; j < stokes_dofs_per_cell; ++j) + if (scratch.dof_component_indices[i] == + scratch.dof_component_indices[j]) + data.local_matrix(i, j) += (- (use_tensor ? + eta_two_thirds * (scratch.div_phi_u[i] * trace(stress_strain_director * scratch.grads_phi_u[j])) + : + eta_two_thirds * (scratch.div_phi_u[i] * scratch.div_phi_u[j]) + )) + * JxW; + } + } + + + + template + void + StokesIncompressibleTermsAnisotropicViscosity:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const + { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast&> (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast&> (data_base); + + const std::shared_ptr> anisotropic_viscosity = + scratch.material_model_outputs.template get_additional_output_object>(); + + const Introspection &introspection = this->introspection(); + const FiniteElement &fe = this->get_fe(); + const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); + const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; + const double pressure_scaling = this->get_pressure_scaling(); + + const std::shared_ptr> force + = scratch.material_model_outputs.template get_additional_output_object>(); + + for (unsigned int q=0; q()); + + const bool use_tensor = (anisotropic_viscosity != nullptr); + + const SymmetricTensor<4, dim> &stress_strain_director = (use_tensor) + ? + anisotropic_viscosity->stress_strain_directors[q] + : + dealii::identity_tensor(); + + const Tensor<1,dim> + gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q)); + + const double density = scratch.material_model_outputs.densities[q]; + const double JxW = scratch.finite_element_values.JxW(q); + + for (unsigned int i=0; irhs_u[q] * scratch.phi_u[i] + + pressure_scaling * force->rhs_p[q] * scratch.phi_p[i]) + * JxW; + + if (scratch.rebuild_stokes_matrix) + for (unsigned int j=0; j + void + StokesIncompressibleTermsAnisotropicViscosity:: + create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const + { + const unsigned int n_points = outputs.viscosities.size(); + + if (outputs.template has_additional_output_object>() == false) + { + outputs.additional_outputs.push_back( + std::make_unique> (n_points)); + } + + if (this->get_parameters().enable_additional_stokes_rhs + && outputs.template has_additional_output_object>() == false) + { + outputs.additional_outputs.push_back( + std::make_unique> (n_points)); + } + Assert(!this->get_parameters().enable_additional_stokes_rhs + || + outputs.template get_additional_output_object>()->rhs_u.size() + == n_points, ExcInternalError()); + } + + + + template + void + StokesCompressibleStrainRateViscosityTermAnisotropicViscosity:: + execute (internal::Assembly::Scratch::ScratchBase &scratch_base, + internal::Assembly::CopyData::CopyDataBase &data_base) const + { + internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast&> (scratch_base); + internal::Assembly::CopyData::StokesSystem &data = dynamic_cast&> (data_base); + + if (!scratch.rebuild_stokes_matrix) + return; + + const std::shared_ptr> anisotropic_viscosity = + scratch.material_model_outputs.template get_additional_output_object>(); + + const Introspection &introspection = this->introspection(); + const FiniteElement &fe = this->get_fe(); + const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); + const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; + + for (unsigned int q=0; q &stress_strain_director = (use_tensor) + ? + anisotropic_viscosity->stress_strain_directors[q] + : + dealii::identity_tensor(); + + const double JxW = scratch.finite_element_values.JxW(q); + + for (unsigned int i=0; i; \ + template class StokesCompressiblePreconditionerAnisotropicViscosity; \ + template class StokesIncompressibleTermsAnisotropicViscosity; \ + template class StokesCompressibleStrainRateViscosityTermAnisotropicViscosity; + + ASPECT_INSTANTIATE(INSTANTIATE) + +#undef INSTANTIATE + } +} diff --git a/tests/anisotropic_viscosity.cc b/tests/anisotropic_viscosity.cc index a097b2985e5..45916434ea4 100644 --- a/tests/anisotropic_viscosity.cc +++ b/tests/anisotropic_viscosity.cc @@ -25,6 +25,7 @@ #include #include #include +#include #include #include @@ -32,592 +33,25 @@ namespace aspect { - namespace MaterialModel - { - /** - * Additional output fields for anisotropic viscosities to be added to - * the MaterialModel::MaterialModelOutputs structure and filled in the - * MaterialModel::Interface::evaluate() function. - */ - template - class AnisotropicViscosity : public NamedAdditionalMaterialOutputs - { - public: - AnisotropicViscosity(const unsigned int n_points); - - virtual std::vector get_nth_output(const unsigned int idx) const; - - /** - * Stress-strain "director" tensors at the given positions. This - * variable is used to implement anisotropic viscosity. - * - * @note The strain rate term in equation (1) of the manual will be - * multiplied by this tensor *and* the viscosity scalar ($\eta$), as - * described in the manual section titled "Constitutive laws". This - * variable is assigned the rank-four identity tensor by default. - * This leaves the isotropic constitutive law unchanged if the material - * model does not explicitly assign a value. - */ - std::vector> stress_strain_directors; - }; - - template - AnisotropicViscosity::AnisotropicViscosity (const unsigned int n_points) - : - NamedAdditionalMaterialOutputs(std::vector(1,"anisotropic_viscosity")), - stress_strain_directors(n_points, dealii::identity_tensor ()) - {} - - - - template - std::vector - AnisotropicViscosity::get_nth_output(const unsigned int /*idx*/) const - { - return std::vector(); - } - } - - namespace Assemblers - { - /** - * A class containing the functions to assemble the Stokes preconditioner. - */ - template - class StokesPreconditionerAnisotropicViscosity : public Assemblers::Interface, - public SimulatorAccess - { - public: - virtual - void - execute(internal::Assembly::Scratch::ScratchBase &scratch, - internal::Assembly::CopyData::CopyDataBase &data) const; - - /** - * Create AnisotropicViscosities. - */ - virtual void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const; - }; - - /** - * A class containing the functions to assemble the compressible adjustment - * to the Stokes preconditioner. - */ - template - class StokesCompressiblePreconditionerAnisotropicViscosity : public Assemblers::Interface, - public SimulatorAccess - { - public: - virtual - void - execute(internal::Assembly::Scratch::ScratchBase &scratch, - internal::Assembly::CopyData::CopyDataBase &data) const; - }; - - /** - * This class assembles the terms for the matrix and right-hand-side of the incompressible - * Stokes equation for the current cell. - */ - template - class StokesIncompressibleTermsAnisotropicViscosity : public Assemblers::Interface, - public SimulatorAccess - { - public: - virtual - void - execute(internal::Assembly::Scratch::ScratchBase &scratch, - internal::Assembly::CopyData::CopyDataBase &data) const; - - /** - * Create AdditionalMaterialOutputsStokesRHS if we need to do so. - */ - virtual void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const; - }; - - /** - * This class assembles the term that arises in the viscosity term of Stokes matrix for - * compressible models, because the divergence of the velocity is not longer zero. - */ - template - class StokesCompressibleStrainRateViscosityTermAnisotropicViscosity : public Assemblers::Interface, - public SimulatorAccess - { - public: - virtual - void - execute(internal::Assembly::Scratch::ScratchBase &scratch, - internal::Assembly::CopyData::CopyDataBase &data) const; - }; - - template - void - StokesPreconditionerAnisotropicViscosity:: - execute (internal::Assembly::Scratch::ScratchBase &scratch_base, - internal::Assembly::CopyData::CopyDataBase &data_base) const - { - internal::Assembly::Scratch::StokesPreconditioner &scratch = dynamic_cast&> (scratch_base); - internal::Assembly::CopyData::StokesPreconditioner &data = dynamic_cast&> (data_base); - - std::shared_ptr> anisotropic_viscosity = - scratch.material_model_outputs.template get_additional_output_object>(); - - const Introspection &introspection = this->introspection(); - const FiniteElement &fe = this->get_fe(); - const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); - const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; - const double pressure_scaling = this->get_pressure_scaling(); - - // First loop over all dofs and find those that are in the Stokes system - // save the component (pressure and dim velocities) each belongs to. - for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/) - { - if (introspection.is_stokes_component(fe.system_to_component_index(i).first)) - { - scratch.dof_component_indices[i_stokes] = fe.system_to_component_index(i).first; - ++i_stokes; - } - ++i; - } - - // Loop over all quadrature points and assemble their contributions to - // the preconditioner matrix - for (unsigned int q = 0; q < n_q_points; ++q) - { - for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/) - { - if (introspection.is_stokes_component(fe.system_to_component_index(i).first)) - { - scratch.grads_phi_u[i_stokes] = - scratch.finite_element_values[introspection.extractors - .velocities].symmetric_gradient(i, q); - scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection - .extractors.pressure].value(i, q); - ++i_stokes; - } - ++i; - } - - const double eta = scratch.material_model_outputs.viscosities[q]; - const double one_over_eta = 1. / eta; - - const bool use_tensor = (anisotropic_viscosity != nullptr); - - const SymmetricTensor<4, dim> &stress_strain_director = (use_tensor) - ? - anisotropic_viscosity->stress_strain_directors[q] - : - dealii::identity_tensor(); - - - - const double JxW = scratch.finite_element_values.JxW(q); - - for (unsigned int i = 0; i < stokes_dofs_per_cell; ++i) - for (unsigned int j = 0; j < stokes_dofs_per_cell; ++j) - if (scratch.dof_component_indices[i] == - scratch.dof_component_indices[j]) - data.local_matrix(i, j) += (( - use_tensor ? - 2.0 * eta * (scratch.grads_phi_u[i] - * stress_strain_director - * scratch.grads_phi_u[j]) : - 2.0 * eta * (scratch.grads_phi_u[i] - * scratch.grads_phi_u[j])) - + one_over_eta * pressure_scaling - * pressure_scaling - * (scratch.phi_p[i] - * scratch.phi_p[j])) - * JxW; - } - } - - - - template - void - StokesPreconditionerAnisotropicViscosity:: - create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const - { - const unsigned int n_points = outputs.viscosities.size(); - - if (outputs.template has_additional_output_object>() == false) - { - outputs.additional_outputs.push_back( - std::make_unique> (n_points)); - } - } - - - - template - void - StokesCompressiblePreconditionerAnisotropicViscosity:: - execute (internal::Assembly::Scratch::ScratchBase &scratch_base, - internal::Assembly::CopyData::CopyDataBase &data_base) const - { - internal::Assembly::Scratch::StokesPreconditioner &scratch = dynamic_cast&> (scratch_base); - internal::Assembly::CopyData::StokesPreconditioner &data = dynamic_cast&> (data_base); - - const std::shared_ptr> anisotropic_viscosity = - scratch.material_model_outputs.template get_additional_output_object>(); - - const Introspection &introspection = this->introspection(); - const FiniteElement &fe = this->get_fe(); - const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); - const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; - - // First loop over all dofs and find those that are in the Stokes system - // save the component (pressure and dim velocities) each belongs to. - for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/) - { - if (introspection.is_stokes_component(fe.system_to_component_index(i).first)) - { - scratch.dof_component_indices[i_stokes] = fe.system_to_component_index(i).first; - ++i_stokes; - } - ++i; - } - - // Loop over all quadrature points and assemble their contributions to - // the preconditioner matrix - for (unsigned int q = 0; q < n_q_points; ++q) - { - for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/) - { - if (introspection.is_stokes_component(fe.system_to_component_index(i).first)) - { - scratch.grads_phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].symmetric_gradient(i,q); - scratch.div_phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].divergence (i, q); - - ++i_stokes; - } - ++i; - } - - const double eta_two_thirds = scratch.material_model_outputs.viscosities[q] * 2.0 / 3.0; - - const bool use_tensor = (anisotropic_viscosity != nullptr); - - const SymmetricTensor<4, dim> &stress_strain_director = (use_tensor) - ? - anisotropic_viscosity->stress_strain_directors[q] - : - dealii::identity_tensor(); - - const double JxW = scratch.finite_element_values.JxW(q); - - for (unsigned int i = 0; i < stokes_dofs_per_cell; ++i) - for (unsigned int j = 0; j < stokes_dofs_per_cell; ++j) - if (scratch.dof_component_indices[i] == - scratch.dof_component_indices[j]) - data.local_matrix(i, j) += (- (use_tensor ? - eta_two_thirds * (scratch.div_phi_u[i] * trace(stress_strain_director * scratch.grads_phi_u[j])) - : - eta_two_thirds * (scratch.div_phi_u[i] * scratch.div_phi_u[j]) - )) - * JxW; - } - } - - - - template - void - StokesIncompressibleTermsAnisotropicViscosity:: - execute (internal::Assembly::Scratch::ScratchBase &scratch_base, - internal::Assembly::CopyData::CopyDataBase &data_base) const - { - internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast&> (scratch_base); - internal::Assembly::CopyData::StokesSystem &data = dynamic_cast&> (data_base); - - const std::shared_ptr> anisotropic_viscosity = - scratch.material_model_outputs.template get_additional_output_object>(); - - const Introspection &introspection = this->introspection(); - const FiniteElement &fe = this->get_fe(); - const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); - const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; - const double pressure_scaling = this->get_pressure_scaling(); - - const std::shared_ptr> force - = scratch.material_model_outputs.template get_additional_output_object>(); - - for (unsigned int q=0; q()); - - const bool use_tensor = (anisotropic_viscosity != nullptr); - - const SymmetricTensor<4, dim> &stress_strain_director = (use_tensor) - ? - anisotropic_viscosity->stress_strain_directors[q] - : - dealii::identity_tensor(); - - const Tensor<1,dim> - gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q)); - - const double density = scratch.material_model_outputs.densities[q]; - const double JxW = scratch.finite_element_values.JxW(q); - - for (unsigned int i=0; irhs_u[q] * scratch.phi_u[i] - + pressure_scaling * force->rhs_p[q] * scratch.phi_p[i]) - * JxW; - - if (scratch.rebuild_stokes_matrix) - for (unsigned int j=0; j - void - StokesIncompressibleTermsAnisotropicViscosity:: - create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const - { - const unsigned int n_points = outputs.viscosities.size(); - - if (outputs.template has_additional_output_object>() == false) - { - outputs.additional_outputs.push_back( - std::make_unique> (n_points)); - } - - if (this->get_parameters().enable_additional_stokes_rhs - && outputs.template has_additional_output_object>() == false) - { - outputs.additional_outputs.push_back( - std::make_unique> (n_points)); - } - Assert(!this->get_parameters().enable_additional_stokes_rhs - || - outputs.template get_additional_output_object>()->rhs_u.size() - == n_points, ExcInternalError()); - } - - - - template - void - StokesCompressibleStrainRateViscosityTermAnisotropicViscosity:: - execute (internal::Assembly::Scratch::ScratchBase &scratch_base, - internal::Assembly::CopyData::CopyDataBase &data_base) const - { - internal::Assembly::Scratch::StokesSystem &scratch = dynamic_cast&> (scratch_base); - internal::Assembly::CopyData::StokesSystem &data = dynamic_cast&> (data_base); - - if (!scratch.rebuild_stokes_matrix) - return; - - const std::shared_ptr> anisotropic_viscosity = - scratch.material_model_outputs.template get_additional_output_object>(); - - const Introspection &introspection = this->introspection(); - const FiniteElement &fe = this->get_fe(); - const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); - const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; - - for (unsigned int q=0; q &stress_strain_director = (use_tensor) - ? - anisotropic_viscosity->stress_strain_directors[q] - : - dealii::identity_tensor(); - - const double JxW = scratch.finite_element_values.JxW(q); - - for (unsigned int i=0; i - class ShearHeatingAnisotropicViscosity : public Interface, public ::aspect::SimulatorAccess - { - public: - /** - * Compute the heating model outputs for this class. - */ - virtual - void - evaluate (const MaterialModel::MaterialModelInputs &material_model_inputs, - const MaterialModel::MaterialModelOutputs &material_model_outputs, - HeatingModel::HeatingModelOutputs &heating_model_outputs) const; - - /** - * Allow the heating model to attach additional material model outputs. - */ - virtual - void - create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &material_model_outputs) const; - }; - - template - void - ShearHeatingAnisotropicViscosity:: - evaluate (const MaterialModel::MaterialModelInputs &material_model_inputs, - const MaterialModel::MaterialModelOutputs &material_model_outputs, - HeatingModel::HeatingModelOutputs &heating_model_outputs) const - { - Assert(heating_model_outputs.heating_source_terms.size() == material_model_inputs.position.size(), - ExcMessage ("Heating outputs need to have the same number of entries as the material model inputs.")); - - // Some material models provide dislocation viscosities and boundary area work fractions - // as additional material outputs. If they are attached, use them. - const std::shared_ptr> shear_heating_out = - material_model_outputs.template get_additional_output_object>(); - - const std::shared_ptr> anisotropic_viscosity = - material_model_outputs.template get_additional_output_object>(); - - for (unsigned int q=0; q &directed_strain_rate = ((anisotropic_viscosity != nullptr) - ? - anisotropic_viscosity->stress_strain_directors[q] - * material_model_inputs.strain_rate[q] - : - material_model_inputs.strain_rate[q]); - - const SymmetricTensor<2,dim> stress = - 2 * material_model_outputs.viscosities[q] * - (this->get_material_model().is_compressible() - ? - directed_strain_rate - 1./3. * trace(directed_strain_rate) * unit_symmetric_tensor() - : - directed_strain_rate); - - const SymmetricTensor<2,dim> deviatoric_strain_rate = - (this->get_material_model().is_compressible() - ? - material_model_inputs.strain_rate[q] - - 1./3. * trace(material_model_inputs.strain_rate[q]) * unit_symmetric_tensor() - : - material_model_inputs.strain_rate[q]); - - heating_model_outputs.heating_source_terms[q] = stress * deviatoric_strain_rate; - - // If shear heating work fractions are provided, reduce the - // overall heating by this amount (which is assumed to be converted into other forms of energy) - if (shear_heating_out != nullptr) - heating_model_outputs.heating_source_terms[q] *= shear_heating_out->shear_heating_work_fractions[q]; - - heating_model_outputs.lhs_latent_heat_terms[q] = 0.0; - } - } - - template - void - ShearHeatingAnisotropicViscosity:: - create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &material_model_outputs) const - { - const unsigned int n_points = material_model_outputs.viscosities.size(); - - if (material_model_outputs.template has_additional_output_object>() == false) - { - material_model_outputs.additional_outputs.push_back( - std::make_unique> (n_points)); - } - - this->get_material_model().create_additional_named_outputs(material_model_outputs); - } - } - namespace MaterialModel { template class Anisotropic : public MaterialModel::Simple { public: - virtual void initialize(); + void initialize() override; - virtual void evaluate(const MaterialModel::MaterialModelInputs &in, - MaterialModel::MaterialModelOutputs &out) const; + void evaluate(const MaterialModel::MaterialModelInputs &in, + MaterialModel::MaterialModelOutputs &out) const override; static void declare_parameters(ParameterHandler &prm); - virtual void parse_parameters(ParameterHandler &prm); + void parse_parameters(ParameterHandler &prm) override; /** * Return true if the compressibility() function returns something that * is not zero. */ - virtual bool - is_compressible () const; + bool + is_compressible () const override; private: void set_assemblers(const SimulatorAccess &, @@ -784,35 +218,6 @@ namespace aspect // explicit instantiations namespace aspect { - namespace Assemblers - { -#define INSTANTIATE(dim) \ - template class StokesPreconditioner; \ - template class StokesCompressiblePreconditioner; \ - template class StokesIncompressibleTerms; \ - template class StokesCompressibleStrainRateViscosityTerm; \ - template class StokesReferenceDensityCompressibilityTerm; \ - template class StokesImplicitReferenceDensityCompressibilityTerm; \ - template class StokesIsentropicCompressionTerm; \ - template class StokesHydrostaticCompressionTerm; \ - template class StokesPressureRHSCompatibilityModification; \ - template class StokesBoundaryTraction; - - ASPECT_INSTANTIATE(INSTANTIATE) - } - - namespace HeatingModel - { - ASPECT_REGISTER_HEATING_MODEL(ShearHeatingAnisotropicViscosity, - "anisotropic shear heating", - "Implementation of a standard model for shear heating. " - "Adds the term: " - "$ 2 \\eta \\left( \\varepsilon - \\frac{1}{3} \\text{tr} " - "\\varepsilon \\mathbf 1 \\right) : \\left( \\varepsilon - \\frac{1}{3} " - "\\text{tr} \\varepsilon \\mathbf 1 \\right)$ to the " - "right-hand side of the temperature equation.") - } - namespace MaterialModel { ASPECT_REGISTER_MATERIAL_MODEL(Anisotropic, diff --git a/tests/anisotropic_viscosity_cookbook.cc b/tests/anisotropic_viscosity_cookbook.cc new file mode 100644 index 00000000000..891cfb86af2 --- /dev/null +++ b/tests/anisotropic_viscosity_cookbook.cc @@ -0,0 +1,21 @@ +/* + Copyright (C) 2022 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + +#include "../cookbooks/anisotropic_viscosity/av_material.cc" diff --git a/tests/anisotropic_viscosity_cookbook.prm b/tests/anisotropic_viscosity_cookbook.prm new file mode 100644 index 00000000000..da8766c281d --- /dev/null +++ b/tests/anisotropic_viscosity_cookbook.prm @@ -0,0 +1,9 @@ +# A test for the anisotropic viscosity cookbook + +include $ASPECT_SOURCE_DIR/cookbooks/anisotropic_viscosity/AV_Rayleigh_Taylor.prm + +set End time = 1 + +subsection Mesh refinement + set Initial global refinement = 4 +end diff --git a/tests/anisotropic_viscosity_cookbook/screen-output b/tests/anisotropic_viscosity_cookbook/screen-output new file mode 100644 index 00000000000..97a600963cc --- /dev/null +++ b/tests/anisotropic_viscosity_cookbook/screen-output @@ -0,0 +1,62 @@ + +Loading shared library <./libanisotropic_viscosity_cookbook.debug.so> + +Number of active cells: 256 (on 5 levels) +Number of degrees of freedom: 6,823 (2,178+289+1,089+1,089+1,089+1,089) + +*** Timestep 0: t=0 seconds, dt=0 seconds + Skipping temperature solve because RHS is zero. + Solving gamma system ... 0 iterations. + Solving ni system ... 11 iterations. + Solving nj system ... 11 iterations. + Solving Stokes system (GMG)... 12+0 iterations. + +Number of active cells: 124 (on 6 levels) +Number of degrees of freedom: 3,618 (1,154+156+577+577+577+577) + +*** Timestep 0: t=0 seconds, dt=0 seconds + Skipping temperature solve because RHS is zero. + Solving gamma system ... 0 iterations. + Solving ni system ... 10 iterations. + Solving nj system ... 10 iterations. + Solving Stokes system (GMG)... 15+0 iterations. + +Number of active cells: 112 (on 6 levels) +Number of degrees of freedom: 3,416 (1,090+146+545+545+545+545) + +*** Timestep 0: t=0 seconds, dt=0 seconds + Skipping temperature solve because RHS is zero. + Solving gamma system ... 0 iterations. + Solving ni system ... 10 iterations. + Solving nj system ... 10 iterations. + Solving Stokes system (GMG)... 16+0 iterations. + + Postprocessing: + Writing graphical output: output-anisotropic_viscosity_cookbook/solution/solution-00000 + RMS, max velocity: 0.00115 m/s, 0.00323 m/s + Compositions min/max/mass: 0/1/0.2988 // -0.09361/0.7853/0.192 // -0.09361/0.7853/0.192 + Average density / Average viscosity / Total mass: 0.1329 kg/m^3, 1 Pa s, 0.2658 kg + +Number of active cells: 130 (on 6 levels) +Number of degrees of freedom: 3,966 (1,266+168+633+633+633+633) + +*** Timestep 1: t=1 seconds, dt=1 seconds + Skipping temperature solve because RHS is zero. + Solving gamma system ... 9 iterations. + Solving ni system ... 10 iterations. + Solving nj system ... 9 iterations. + Solving Stokes system (GMG)... 16+0 iterations. + + Postprocessing: + Writing graphical output: output-anisotropic_viscosity_cookbook/solution/solution-00001 + RMS, max velocity: 0.00119 m/s, 0.00324 m/s + Compositions min/max/mass: -0.0005129/1.014/0.2987 // -0.08678/0.787/0.1925 // -0.08677/0.7925/0.1922 + Average density / Average viscosity / Total mass: 0.1329 kg/m^3, 1 Pa s, 0.2657 kg + +Number of active cells: 145 (on 6 levels) +Number of degrees of freedom: 4,435 (1,416+187+708+708+708+708) + +Termination requested by criterion: end time + + + diff --git a/tests/anisotropic_viscosity_cookbook/statistics b/tests/anisotropic_viscosity_cookbook/statistics new file mode 100644 index 00000000000..f18f1daf3e3 --- /dev/null +++ b/tests/anisotropic_viscosity_cookbook/statistics @@ -0,0 +1,31 @@ +# 1: Time step number +# 2: Time (seconds) +# 3: Time step size (seconds) +# 4: Number of mesh cells +# 5: Number of Stokes degrees of freedom +# 6: Number of temperature degrees of freedom +# 7: Number of degrees of freedom for all compositions +# 8: Iterations for temperature solver +# 9: Iterations for composition solver 1 +# 10: Iterations for composition solver 2 +# 11: Iterations for composition solver 3 +# 12: Iterations for Stokes solver +# 13: Velocity iterations in Stokes preconditioner +# 14: Schur complement iterations in Stokes preconditioner +# 15: Visualization file name +# 16: RMS velocity (m/s) +# 17: Max. velocity (m/s) +# 18: Minimal value for composition gamma +# 19: Maximal value for composition gamma +# 20: Global mass for composition gamma +# 21: Minimal value for composition ni +# 22: Maximal value for composition ni +# 23: Global mass for composition ni +# 24: Minimal value for composition nj +# 25: Maximal value for composition nj +# 26: Global mass for composition nj +# 27: Average density (kg/m^3) +# 28: Average viscosity (Pa s) +# 29: Total mass (kg) +0 0.000000000000e+00 0.000000000000e+00 112 1236 545 1635 0 0 10 10 16 17 17 output-anisotropic_viscosity_cookbook/solution/solution-00000 1.14927553e-03 3.23431850e-03 0.00000000e+00 9.99999959e-01 2.98754046e-01 -9.36092543e-02 7.85323275e-01 1.92049844e-01 -9.36092543e-02 7.85323275e-01 1.92049844e-01 1.32914979e-01 1.00000000e+00 2.65829958e-01 +1 1.000000000000e+00 1.000000000000e+00 130 1434 633 1899 0 9 10 9 16 17 17 output-anisotropic_viscosity_cookbook/solution/solution-00001 1.18953522e-03 3.24176496e-03 -5.12893159e-04 1.01443242e+00 2.98668555e-01 -8.67764149e-02 7.86971956e-01 1.92470152e-01 -8.67669395e-02 7.92497164e-01 1.92244641e-01 1.32854697e-01 1.00000000e+00 2.65709394e-01