From 5ed666d0ed81c13dde7decccc5f5d5c193412537 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Thu, 23 Oct 2025 16:53:03 +0200 Subject: [PATCH 1/6] Move anisotropic viscosity assemblers into main code --- .../assemblers/stokes_anisotropic_viscosity.h | 172 ++++++ .../stokes_anisotropic_viscosity.cc | 399 ++++++++++++++ tests/anisotropic_viscosity.cc | 490 +----------------- 3 files changed, 572 insertions(+), 489 deletions(-) create mode 100644 include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h create mode 100644 source/simulator/assemblers/stokes_anisotropic_viscosity.cc 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..e3cd494baf7 --- /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/simulator/assemblers/stokes_anisotropic_viscosity.cc b/source/simulator/assemblers/stokes_anisotropic_viscosity.cc new file mode 100644 index 00000000000..3949e371549 --- /dev/null +++ b/source/simulator/assemblers/stokes_anisotropic_viscosity.cc @@ -0,0 +1,399 @@ +/* + 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 + +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..fa73dd4d2f0 100644 --- a/tests/anisotropic_viscosity.cc +++ b/tests/anisotropic_viscosity.cc @@ -25,6 +25,7 @@ #include #include #include +#include #include #include @@ -32,478 +33,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); - - 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 @@ -784,23 +313,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, From c5f2883e54fe2d2a29bb91b29629aed93eb35d4f Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Thu, 23 Oct 2025 17:20:08 +0200 Subject: [PATCH 2/6] Add test for anisotropic viscosity cookbook --- tests/anisotropic_viscosity_cookbook.cc | 21 +++++++ tests/anisotropic_viscosity_cookbook.prm | 9 +++ .../screen-output | 62 +++++++++++++++++++ .../anisotropic_viscosity_cookbook/statistics | 31 ++++++++++ 4 files changed, 123 insertions(+) create mode 100644 tests/anisotropic_viscosity_cookbook.cc create mode 100644 tests/anisotropic_viscosity_cookbook.prm create mode 100644 tests/anisotropic_viscosity_cookbook/screen-output create mode 100644 tests/anisotropic_viscosity_cookbook/statistics 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 From 80347274c446e2a2aa27a1adecf08327c074cb8c Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Thu, 23 Oct 2025 17:20:22 +0200 Subject: [PATCH 3/6] fix anisotropic viscosity parameter file --- cookbooks/anisotropic_viscosity/AV_Rayleigh_Taylor.prm | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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; From 6b74606db83ddc7d845a05c2f8c4252b665d0c67 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Thu, 23 Oct 2025 17:27:20 +0200 Subject: [PATCH 4/6] Remove duplicated assembler from cookbook --- .../anisotropic_viscosity/av_material.cc | 323 +----------------- .../assemblers/stokes_anisotropic_viscosity.h | 38 +-- .../stokes_anisotropic_viscosity.cc | 2 + 3 files changed, 22 insertions(+), 341 deletions(-) diff --git a/cookbooks/anisotropic_viscosity/av_material.cc b/cookbooks/anisotropic_viscosity/av_material.cc index 7bc25da38e3..13e7a4d2dae 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 @@ -60,317 +60,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 @@ -792,16 +481,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, diff --git a/include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h b/include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h index e3cd494baf7..a5173ca85d0 100644 --- a/include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h +++ b/include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h @@ -27,7 +27,7 @@ namespace aspect { - namespace MaterialModel + namespace MaterialModel { /** * Additional output fields for anisotropic viscosities to be added to @@ -40,7 +40,7 @@ namespace aspect public: AnisotropicViscosity(const unsigned int n_points); - std::vector + std::vector get_nth_output(const unsigned int idx) const override; /** @@ -59,19 +59,19 @@ namespace aspect 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 + 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; + } } - } @@ -99,9 +99,9 @@ namespace aspect namespace Assemblers { - /** - * A class containing the functions to assemble the Stokes preconditioner. - */ + /** + * A class containing the functions to assemble the Stokes preconditioner. + */ template class StokesPreconditionerAnisotropicViscosity : public Assemblers::Interface, public SimulatorAccess @@ -114,7 +114,7 @@ namespace aspect /** * Create AnisotropicViscosities. */ - void + void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const override; }; @@ -148,7 +148,7 @@ namespace aspect /** * Create AdditionalMaterialOutputsStokesRHS if we need to do so. */ - void + void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &outputs) const override; }; diff --git a/source/simulator/assemblers/stokes_anisotropic_viscosity.cc b/source/simulator/assemblers/stokes_anisotropic_viscosity.cc index 3949e371549..5b9cffe66eb 100644 --- a/source/simulator/assemblers/stokes_anisotropic_viscosity.cc +++ b/source/simulator/assemblers/stokes_anisotropic_viscosity.cc @@ -20,6 +20,8 @@ #include +#include + #include namespace aspect From e544695de644a5dca9b0c59694f74f4dcb5b0875 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Wed, 29 Oct 2025 13:57:50 +0100 Subject: [PATCH 5/6] Add shear heating plugin --- .../anisotropic_viscosity/av_material.cc | 113 ---------------- .../shear_heating_anisotropic_viscosity.h | 61 +++++++++ .../shear_heating_anisotropic_viscosity.cc | 123 ++++++++++++++++++ tests/anisotropic_viscosity.cc | 119 +---------------- 4 files changed, 190 insertions(+), 226 deletions(-) create mode 100644 include/aspect/heating_model/shear_heating_anisotropic_viscosity.h create mode 100644 source/heating_model/shear_heating_anisotropic_viscosity.cc diff --git a/cookbooks/anisotropic_viscosity/av_material.cc b/cookbooks/anisotropic_viscosity/av_material.cc index 13e7a4d2dae..1770bdd3fba 100644 --- a/cookbooks/anisotropic_viscosity/av_material.cc +++ b/cookbooks/anisotropic_viscosity/av_material.cc @@ -48,8 +48,6 @@ #include #include -#include -#include #include #include #include @@ -60,103 +58,6 @@ namespace aspect { - 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 @@ -481,20 +382,6 @@ namespace aspect // explicit instantiations namespace aspect { - 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/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/tests/anisotropic_viscosity.cc b/tests/anisotropic_viscosity.cc index fa73dd4d2f0..45916434ea4 100644 --- a/tests/anisotropic_viscosity.cc +++ b/tests/anisotropic_viscosity.cc @@ -33,120 +33,25 @@ namespace aspect { - namespace HeatingModel - { - template - 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 &, @@ -313,18 +218,6 @@ namespace aspect // explicit instantiations namespace aspect { - 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, From 51c647a20f2984c2a54daf489c7dc8c207f1c38e Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Tue, 11 Nov 2025 15:31:57 +0100 Subject: [PATCH 6/6] Improve anisotropic viscosity code --- .../anisotropic_viscosity/av_material.cc | 2 +- .../shear_heating_anisotropic_viscosity.h | 4 +- .../anisotropic_viscosity.h | 61 ++++++++++ include/aspect/simulator/assemblers/stokes.h | 2 +- .../assemblers/stokes_anisotropic_viscosity.h | 82 ++------------ .../shear_heating_anisotropic_viscosity.cc | 14 +-- .../anisotropic_viscosity.cc | 80 ++++++++++++++ .../stokes_anisotropic_viscosity.cc | 104 +++++++----------- tests/anisotropic_viscosity.cc | 1 + 9 files changed, 198 insertions(+), 152 deletions(-) create mode 100644 include/aspect/material_model/additional_outputs/anisotropic_viscosity.h create mode 100644 source/material_model/additional_outputs/anisotropic_viscosity.cc diff --git a/cookbooks/anisotropic_viscosity/av_material.cc b/cookbooks/anisotropic_viscosity/av_material.cc index 1770bdd3fba..4a548ae8bcd 100644 --- a/cookbooks/anisotropic_viscosity/av_material.cc +++ b/cookbooks/anisotropic_viscosity/av_material.cc @@ -20,7 +20,7 @@ #include #include -#include +#include #include #include #include diff --git a/include/aspect/heating_model/shear_heating_anisotropic_viscosity.h b/include/aspect/heating_model/shear_heating_anisotropic_viscosity.h index 5614ff9334e..ad86514757e 100644 --- a/include/aspect/heating_model/shear_heating_anisotropic_viscosity.h +++ b/include/aspect/heating_model/shear_heating_anisotropic_viscosity.h @@ -31,8 +31,8 @@ namespace aspect { /** * 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 + * 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 diff --git a/include/aspect/material_model/additional_outputs/anisotropic_viscosity.h b/include/aspect/material_model/additional_outputs/anisotropic_viscosity.h new file mode 100644 index 00000000000..a48c30f5467 --- /dev/null +++ b/include/aspect/material_model/additional_outputs/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 doc/COPYING. If not see + . +*/ + +#ifndef _aspect_material_model_additional_outputs_anisotropic_viscosity_h +#define _aspect_material_model_additional_outputs_anisotropic_viscosity_h + + +#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; + }; + } +} + +#endif diff --git a/include/aspect/simulator/assemblers/stokes.h b/include/aspect/simulator/assemblers/stokes.h index 9ec2a5549fa..a45558e78e5 100644 --- a/include/aspect/simulator/assemblers/stokes.h +++ b/include/aspect/simulator/assemblers/stokes.h @@ -78,7 +78,7 @@ namespace aspect }; /** - * This class assembles the term that arises in the viscosity term of Stokes matrix for + * This class assembles the term that arises in the viscosity term of the Stokes matrix for * compressible models, because the divergence of the velocity is not longer zero. */ template diff --git a/include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h b/include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h index a5173ca85d0..a766150f83b 100644 --- a/include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h +++ b/include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h @@ -27,80 +27,11 @@ 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. + * A class containing the functions to assemble the Stokes preconditioner for the + * case of anisotropic viscosities. */ template class StokesPreconditionerAnisotropicViscosity : public Assemblers::Interface, @@ -120,7 +51,7 @@ namespace aspect /** * A class containing the functions to assemble the compressible adjustment - * to the Stokes preconditioner. + * to the Stokes preconditioner for the case of anisotropic viscosities. */ template class StokesCompressiblePreconditionerAnisotropicViscosity : public Assemblers::Interface, @@ -134,7 +65,7 @@ namespace aspect /** * This class assembles the terms for the matrix and right-hand-side of the incompressible - * Stokes equation for the current cell. + * Stokes equation for the case of anisotropic viscosities. */ template class StokesIncompressibleTermsAnisotropicViscosity : public Assemblers::Interface, @@ -153,8 +84,9 @@ namespace aspect }; /** - * 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. + * This class assembles the term that arises in the viscosity term of the Stokes matrix for + * compressible models for the case of anisotropic viscosities, because the divergence + * of the velocity is not longer zero. */ template class StokesCompressibleStrainRateViscosityTermAnisotropicViscosity : public Assemblers::Interface, diff --git a/source/heating_model/shear_heating_anisotropic_viscosity.cc b/source/heating_model/shear_heating_anisotropic_viscosity.cc index fed7445582e..da10e37bac8 100644 --- a/source/heating_model/shear_heating_anisotropic_viscosity.cc +++ b/source/heating_model/shear_heating_anisotropic_viscosity.cc @@ -21,8 +21,8 @@ #include #include +#include #include -#include #include #include @@ -49,15 +49,15 @@ namespace aspect const std::shared_ptr> anisotropic_viscosity = material_model_outputs.template get_additional_output_object>(); + Assert(anisotropic_viscosity != nullptr, + ExcMessage("This heating plugin should only be used with material models that provide " + "an anisotropic viscosity tensor, but none was provided.")); + 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> &directed_strain_rate = anisotropic_viscosity->stress_strain_directors[q] + * material_model_inputs.strain_rate[q]; const SymmetricTensor<2,dim> stress = 2 * material_model_outputs.viscosities[q] * diff --git a/source/material_model/additional_outputs/anisotropic_viscosity.cc b/source/material_model/additional_outputs/anisotropic_viscosity.cc new file mode 100644 index 00000000000..3f84b86dbc6 --- /dev/null +++ b/source/material_model/additional_outputs/anisotropic_viscosity.cc @@ -0,0 +1,80 @@ +/* + 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 + . +*/ + +#include + +namespace aspect +{ + namespace MaterialModel + { + namespace + { + template + std::vector make_anisotropic_viscosity_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_anisotropic_viscosity_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; + } + } +} + +// explicit instantiations +namespace aspect +{ + namespace MaterialModel + { +#define INSTANTIATE(dim) \ + template class AnisotropicViscosity; + + ASPECT_INSTANTIATE(INSTANTIATE) + +#undef INSTANTIATE + } +} diff --git a/source/simulator/assemblers/stokes_anisotropic_viscosity.cc b/source/simulator/assemblers/stokes_anisotropic_viscosity.cc index 5b9cffe66eb..0ec00b751c7 100644 --- a/source/simulator/assemblers/stokes_anisotropic_viscosity.cc +++ b/source/simulator/assemblers/stokes_anisotropic_viscosity.cc @@ -21,6 +21,7 @@ #include #include +#include #include @@ -37,9 +38,13 @@ namespace aspect internal::Assembly::Scratch::StokesPreconditioner &scratch = dynamic_cast&> (scratch_base); internal::Assembly::CopyData::StokesPreconditioner &data = dynamic_cast&> (data_base); - std::shared_ptr> anisotropic_viscosity = + const std::shared_ptr> anisotropic_viscosity = scratch.material_model_outputs.template get_additional_output_object>(); + Assert(anisotropic_viscosity != nullptr, + ExcMessage("This assembler should only be used with material models that provide " + "an anisotropic viscosity tensor, but none was provided.")); + const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); @@ -78,30 +83,16 @@ namespace aspect 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 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) += (( - 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])) + 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] @@ -140,6 +131,10 @@ namespace aspect const std::shared_ptr> anisotropic_viscosity = scratch.material_model_outputs.template get_additional_output_object>(); + Assert(anisotropic_viscosity != nullptr, + ExcMessage("This assembler should only be used with material models that provide " + "an anisotropic viscosity tensor, but none was provided.")); + const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); @@ -174,26 +169,15 @@ namespace aspect } 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 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) += (- (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]) - )) + data.local_matrix(i, j) += (- eta_two_thirds * + (scratch.div_phi_u[i] * trace(stress_strain_director * scratch.grads_phi_u[j]))) * JxW; } } @@ -212,6 +196,10 @@ namespace aspect const std::shared_ptr> anisotropic_viscosity = scratch.material_model_outputs.template get_additional_output_object>(); + Assert(anisotropic_viscosity != nullptr, + ExcMessage("This assembler should only be used with material models that provide " + "an anisotropic viscosity tensor, but none was provided.")); + const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); @@ -247,13 +235,7 @@ namespace aspect : numbers::signaling_nan()); - 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 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)); @@ -274,18 +256,15 @@ namespace aspect if (scratch.rebuild_stokes_matrix) for (unsigned int j=0; j> anisotropic_viscosity = scratch.material_model_outputs.template get_additional_output_object>(); + Assert(anisotropic_viscosity != nullptr, + ExcMessage("This assembler should only be used with material models that provide " + "an anisotropic viscosity tensor, but none was provided.")); + const Introspection &introspection = this->introspection(); const FiniteElement &fe = this->get_fe(); const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); @@ -357,25 +340,14 @@ namespace aspect // Viscosity scalar 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 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 #include +#include #include #include #include