Unify anisotropic viscosity Stokes assemblers / Fix anisotropic viscosity cookbook #6728
Conversation
gassmoeller
force-pushed
the
unify_anisotropic_viscosity_code
branch
from
7d29297 to
dbf80e1
Compare
|
/rebuild |
gassmoeller
force-pushed
the
unify_anisotropic_viscosity_code
branch
from
dbf80e1 to
e544695
Compare
|
/rebuild |
| subsection Function | ||
| set Variable names = x,z | ||
| set Function constants = pi=3.1415926; | ||
| set Function constants = pi=3.1415926 |
There was a problem hiding this comment.
Did that yield an error? If not, our parsing is not robust enough...
bangerth
left a comment
bangerth
left a comment
There was a problem hiding this comment.
I still have to look at the .cc file for the assembly, but here are a few comments already.
| * anisotropic viscosity tensor. If the material model provides a stress- | ||
| * strain director tensor, then the strain-rate is multiplied with this |
There was a problem hiding this comment.
| * 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 |
| const std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity = | ||
| material_model_outputs.template get_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>(); | ||
|
|
||
| for (unsigned int q=0; q<heating_model_outputs.heating_source_terms.size(); ++q) | ||
| { | ||
| // If there is an anisotropic viscosity, use it to compute the correct stress | ||
| const SymmetricTensor<2,dim> &directed_strain_rate = ((anisotropic_viscosity != nullptr) | ||
| ? | ||
| anisotropic_viscosity->stress_strain_directors[q] | ||
| * material_model_inputs.strain_rate[q] | ||
| : | ||
| material_model_inputs.strain_rate[q]); |
There was a problem hiding this comment.
Would it make sense to just always expect that a director is given when using this heating model?
| 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 <int dim> | ||
| class AnisotropicViscosity : public NamedAdditionalMaterialOutputs<dim> | ||
| { | ||
| public: | ||
| AnisotropicViscosity(const unsigned int n_points); |
There was a problem hiding this comment.
The assemblers directory is a bit of a funny place to find the declaration of a material model (or named additional output, as it may be). Have you thought moving it into the respective directory?
| namespace | ||
| { | ||
| template <int dim> | ||
| std::vector<std::string> make_AnisotropicViscosity_additional_outputs_names() | ||
| { | ||
| std::vector<std::string> 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; | ||
| } | ||
| } |
There was a problem hiding this comment.
Modules don't allow for anonymous namespaces in header files, and it's generally a bad idea anyway. I think that you're using this function only immediately below in one place, so it should be easy to move this code into its only caller.
| template <int dim> | ||
| AnisotropicViscosity<dim>::AnisotropicViscosity (const unsigned int n_points) | ||
| : | ||
| NamedAdditionalMaterialOutputs<dim>(make_AnisotropicViscosity_additional_outputs_names<dim>()), |
There was a problem hiding this comment.
Just convert the function above into a lambda function that you can declare and call in-place here.
| template <int dim> | ||
| std::vector<double> | ||
| AnisotropicViscosity<dim>::get_nth_output(const unsigned int idx) const | ||
| { | ||
| std::vector<double> 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; | ||
| } |
There was a problem hiding this comment.
This function could probably just as well live in a .cc file.
| namespace Assemblers | ||
| { | ||
| /** | ||
| * A class containing the functions to assemble the Stokes preconditioner. |
There was a problem hiding this comment.
| * 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. |
| * A class containing the functions to assemble the compressible adjustment | ||
| * to the Stokes preconditioner. |
There was a problem hiding this comment.
| * A class containing the functions to assemble the compressible adjustment | |
| * to the Stokes preconditioner. | |
| * A class containing the functions to assemble the compressible adjustment | |
| * to the Stokes preconditioner for the case of anisotropic viscosities. |
| * This class assembles the terms for the matrix and right-hand-side of the incompressible | ||
| * Stokes equation for the current cell. |
There was a problem hiding this comment.
same here, and perhaps also below
| }; | ||
|
|
||
| /** | ||
| * This class assembles the term that arises in the viscosity term of Stokes matrix for |
There was a problem hiding this comment.
| * 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 |
This is a copy-paste mistake -- the grammatically wrong code is already in line include/aspect/simulator/assemblers/stokes.h:81.
| std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity = | ||
| scratch.material_model_outputs.template get_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>(); |
There was a problem hiding this comment.
Could this variable be const?
| const double eta = scratch.material_model_outputs.viscosities[q]; | ||
| const double one_over_eta = 1. / eta; | ||
|
|
||
| const bool use_tensor = (anisotropic_viscosity != nullptr); |
There was a problem hiding this comment.
This variable should move out of the loop. But more generally, shouldn't this assemble only be called for the anisotropic case? We could assert that at the top, and then simplify a lot of the code below. (I recognize the history of the code as what we used to have, where we catered to both the anisotropic and isotropic case, but perhaps we don't need to carry this provenance around with us for forever.)
| template <int dim> | ||
| void | ||
| StokesCompressiblePreconditionerAnisotropicViscosity<dim>:: | ||
| execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch_base, | ||
| internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const |
There was a problem hiding this comment.
same comments for this function as for the one above.
| template <int dim> | ||
| void | ||
| StokesIncompressibleTermsAnisotropicViscosity<dim>:: | ||
| execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch_base, | ||
| internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const |
| template <int dim> | ||
| void | ||
| StokesCompressibleStrainRateViscosityTermAnisotropicViscosity<dim>:: | ||
| execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch_base, | ||
| internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const |
2 participants
In #2139 we had moved the anisotropic viscosity assemblers out of the main code and into the tests, since there was no realistic application for them at the time. #3066 then introduced a cookbook that duplicated the existing code from the test in a plugin. #6615 now needs the same code again, so I think it is time to move the code back into the main program and unify the code. The assemblers are separated from the normal Stokes assemblers, because they change things in the innermost loop of the Stokes assembler and are significantly more expensive than the isotropic viscosity ones (that was the motivation in #2139 to split them from the normal assemblers).
To make sure nothing was changed between the cookbook and the test assemblers, I first created a test for the cookbook (meanwhile fixing a bug in the cookbook prm), and then unified the assembler code. The results of the test remained unchanged.
The anisotropic viscosity assemblers are still not accessible from a simple parameter file change, they need a plugin to replace the isotropic assemblers with the anisotropic ones. But at least we now dont need to duplicate the code anymore.
I have now also unified the shear heating plugins into a new heating plugin.
@Wang-yijun this will simplify your work in #6615