diff --git a/cookbooks/CPO_induced_anisotropic_viscosity/doc/CPO_induced_anisotropic_viscosity.md b/cookbooks/CPO_induced_anisotropic_viscosity/doc/CPO_induced_anisotropic_viscosity.md new file mode 100644 index 00000000000..61cf781e0f2 --- /dev/null +++ b/cookbooks/CPO_induced_anisotropic_viscosity/doc/CPO_induced_anisotropic_viscosity.md @@ -0,0 +1,141 @@ +```{tags} +category:cookbook +feature:3d +feature:cartesian +feature:modular-equations +feature:compositional-fields +feature:particles + +(sec:cookbooks:CPO_induced_anisotropic_viscosity)= +# CPO induced anisotropic viscosity + +*This section was contributed by Yijun Wang and Ágnes Király.* + +This cookbook explains how to use the CPO-induced anisotropic viscosity material model to set up a shear box model. + +## Introduction + +Individual crystals of the mineral olivine reorganize their orientations into the crystal-preferred orientations (CPO) under deformation. The viscous properties of olivine crystals are direction-dependent (anisotropic), which suggests that the effective viscosity for olivine rocks/aggregates is different when deformations occur in different directions relative to the CPO. This cookbook model computes an anisotropic viscosity based on the CPO evolution predicted by D-Rex ({cite}`fraters_billen_2021_cpo`; {cite}`kaminski2004`) and includes this information in the subsequent modeling process. + +Our constitutive equation for the relationship between the strain rate and stress using the anisotropic viscosity tensor is adapted from {cite:t}`signorelli:etal:2021`: + +```{math} +:label: eqn:anisotropic_general_stress +\dot{\varepsilon}_{ij} = \gamma J(\sigma_{ij})^{(n-1)} A_{ij} \sigma_{ij} +``` + +where $\gamma$ is the part of fluidity (the inverse of viscosity) which is temperature- and grain-size dependent: + +```{math} +:label: eqn:fluidity +\gamma=\gamma_0 exp \left(\frac{-Q}{RT} \right) /d^m +``` + +$\gamma_0=1.1\times 10^{5}$ is the isotropic fluidity, $Q=530$ $kJ/mol$ is the activation energy, $R=8.314 m^3 \cdot Pa \cdot K^{−1} \cdot$ $mol^{−1}$ is the gas constant, $d=0.001$ $m$ is the grain size, and $m=0.73$ is the grain size exponent. These values for olivine are taken from rock experiments performed by {cite:t}`hansen:etal:2016` and {cite:t}`HK04`. $J(\sigma_{ij})$ is the equivalent yield stress, where $\sigma_{ij}$ is the deviatoric (anisotropic) stress computed using the tensorial and scalar component of the anisotropic viscosity: + +```{math} +:label: eqn:equivalent_yield_stress +J(\sigma_{ij})=(F(\sigma_{11} - \sigma_{22})^2+G(\sigma_{22} - \sigma_{33})^2+H(\sigma_{33} - \sigma_{11})^2+2L\sigma_{23}^2+2M\sigma_{13}^2+2N\sigma_{12}^2)^{1/2} +``` + +and $A_{ij}$ is the anisotropic tensor of fluidity in Kelvin notation: + +```{math} +:label: eqn:anisotropic_fluidity +A_{ij}=\frac{2}{3} \left[ +\begin{matrix} +F+H & -F & -H & 0 & 0 & 0 \\ +-F & G+F & -G & 0 & 0 & 0 \\ +-H & -G & H+G & 0 & 0 & 0 \\ +0 & 0 & 0 & L & 0 & 0 \\ +0 & 0 & 0 & 0 & M & 0 \\ +0 & 0 & 0 & 0 & 0 & N +\end{matrix} \right] +``` + +$J(\sigma_{ij})$ and $A_{ij}$ are computed using Hill coefficients $H, J, K, L, M,$ and $N$ {cite}`hill:1948`, which are obtained from regression analysis of a texture database constructed with olivine textures from laboratory experiments, shear box models, and subduction models (Kiraly et al., in rev.). + +In this material model plugin, strain rate, density, temperature, and other parameters are taken as input to compute the anisotropic viscosity, which is passed into the Stokes system to compute the stress. As a result, we adapt {math:numref}`eqn:anisotropic_general_stress` to be: + +```{math} +:label: eqn:anisotropic_stress +\sigma_{ij} = \frac{1}{\gamma J(\sigma_{ij})^{(n-1)}} * A_{ij}^{-1} * \dot\varepsilon_{ij} +``` + +Since the Hill coefficients are defined in the microscopic CPO reference frame, and parameters computed in ASPECT are in the macroscopic model reference frame, several reference frame conversions are needed. First, we need to rotate $\sigma_{ij}$ in {math:numref}`eqn:equivalent_yield_stress` from the model reference frame to the CPO reference frame so that $J(\sigma_{ij})$ is in the CPO reference frame. This is achieved by constructing a matrix from the eigenvectors corresponding with the largest eigenvalues of the covariance matrix for the a-, b-, and c-axis of olivine textures and then we assign the nearest orthogonal matrix to be the rotation matrix R: + +```{math} +:label: eqn:rotation_matrix +R = \left[ +\begin{matrix} +\verb|max_eigenvector|\_{a1} & \verb|max_eigenvector|\_{b1} & \verb|max_eigenvector|\_{c1} \\ +\verb|max_eigenvector|\_{a2} & \verb|max_eigenvector|\_{b2} & \verb|max_eigenvector|\_{c2} \\ +\verb|max_eigenvector|\_{a3} & \verb|max_eigenvector|\_{b3} & \verb|max_eigenvector|\_{c3} \\ +\end{matrix} \right] +``` + +We compute the rotation matrix R on the particles and further convert it to Euler angles for computation and memory efficiency. These properties need to be interpolated from particles to fields to be used in the material model. As a result, the anisotropic viscosity material model requires at least one particle in each cell so that all cells can have the texture parameters (Euler angles and eigenvalues) for constructing the rotation matrix R and compute the Hill coefficients. In the material model, the interpolated Euler angles are converted to the rotation matrix again. We use the same notation R to describe the rotation matrix used in the material model in the following paragraphs. + +The inverse of the A tensor then needs to be rotated to the model reference frame. Since $A_{ij}^{-1}$ is the Kelvin notation of the rank-4 tensor, we apply the Kelvin notation representation of the R rotation matrix, $R_K$, on $A_{ij}^{-1}$: + +```{math} +:label: eqn:rotation_matrix_kelvin +R_K = \left[ +\begin{matrix} +R_{11}^2 & R_{12}^2 & R_{13}^2 & \sqrt2* R_{12}* R_{13} & \sqrt2* R_{11}* R_{13} & \sqrt2* R_{11}* R_{12} \\ +R_{21}^2 & R_{22}^2 & R_{23}^2 & \sqrt2* R_{22}* R_{23} & \sqrt2* R_{21}* R_{23} & \sqrt2* R_{21}* R_{22} \\ +R_{31}^2 & R_{32}^2 & R_{33}^2 & \sqrt2* R_{32}* R_{33} & \sqrt2* R_{31}* R_{33} & \sqrt2* R_{31}* R_{32} \\ +\sqrt2* R_{21}* R_{31} & \sqrt2* R_{23}* R_{32} & \sqrt2* R_{23}* R_{33} & R_{22}* R_{33}+R_{23}* R_{32} & R_{21}* R_{33}+R_{23}* R_{31} & R_{21}* R_{32}+R_{22}* R_{31} \\ +\sqrt2* R_{11}* R_{31} & \sqrt2* R_{12}* R_{32} & \sqrt2* R_{13}* R_{33} & R_{12}* R_{33}+R_{13}* R_{32} & R_{11}* R_{33}+R_{13}* R_{31} & R_{11}* R_{32}+R_{12}* R_{31} \\ +\sqrt2* R_{11}* R_{21} & \sqrt2* R_{12}* R_{22} & \sqrt2* R_{13}* R_{23} & R_{12}* R_{23}+R_{13}* R_{22} & R_{11}* R_{23}+R_{13}* R_{32} & R_{11}* R_{22}+R_{12}* R_{21} +\end{matrix} \right] +``` + +The final equation involving all reference frame conversions is: + +```{math} +:label: eqn:anisotropic_stress_final +\sigma_{ij} = \frac{1}{\gamma J(R'*\sigma_{ij}* R)^{(n-1)}} * (R_K * A_{ij}^{-1} * R_K') * \dot\varepsilon_i +``` + +$R'$ and $R_K'$ is the transpose of matrix $R$ and $R_K$ respectively. We save $\frac{1}{\gamma J(R'*\sigma_{ij}* R)^{(n-1)}}$ as the material model viscosity output, which we call the scalar viscosity. The tensorial part of anisotropic viscosity, $R_K * A_{ij}^{-1} * R_K'$ is called the stress-strain director and is stored in the additional outputs. Due to the dependence of the scalar viscosity on stress and the stress is determined by the scalar viscosity from the previous time step, the computation of the scalar viscosity is not stable and its value naturally oscillates. This fluctuation ultimately causes a numerical instability associated with the prediction of the anisotropic viscosity. Therefore, we damp the scalar viscosity using a non-linear Newton iteration, and in each iteration, only half of the change in the scalar viscosity is applied until the result is stable. + + + +## Compiling requirement + +Since the determinant of $A_{ij}$ is 0, $A_{ij}$ is a singular, non-invertible matrix. We find the MoorePenrose pseudoinverse of the matrix $A_{ij}$ as an approximate of the inverse of $A_{ij}$, using the SCALAPACK package provided in deal.II. Thus it is required to link ASPECT with a deal.II version with the scalapack package included in order to run this cookbook. When using candi you can enable the scalapack package by including `scalapack` in the list of installed packages or uncommenting the line in `candi.cfg` that corresponds to the scalapack installation. Alternatively, you can install ScaLAPACK yourself and enable the configuration variable `DEAL_II_WITH_SCALAPACK` during the cmake configuration of deal.II. + +## Model setup + +The usage of the AV material model is demonstrated with a 3d simple shear box model, where its dimension is $1 \times 1 \times 1 $ (non-dimensionalized). The shear strain rate is set to +$0.5$. The origin is the center of the box, and one Olivine particle with 1000 grains sits at the origin to track CPO developments for computation of anisotropic viscosity parameters. + +Since the AV material model computes viscosity based on the evolving CPO stored on particles, several setup requirements must be met: + +- **Particles per cell**: Each computational cell must contain at least one particle, to allow interpolation of the CPO particle property. This is achieved by setting (in the Particles subsection): + +```{literalinclude} min_particles_per_cell.part.prm +``` + +- **CPO particle property**: The CPO particle property must be stored for use by the AV model. This requires enabling the particle and crystal preferred orientation postprocessors and the relevant subsuctions for them, including the CPO Bingham Average plugin, which calculates the Hill coefficients: + +```{literalinclude} cpo_particle_property.part.prm +``` + +Note: These settings are similar to those used for simulations involving CPO alone. However, for the AV model, it is essential to set `Use rotation matrix = false` in the CPO Bingham Average subsection, so that the CPO is represented using Euler angles, as required. + +- **Compositional fields**: The eigenvalues and Euler angles of the CPO tensor are stored in compositional fields. This requires the following input file section: + +```{literalinclude} compositional_field.part.prm +``` + +In the `CPO induced Anisotropic Viscosity` material model subsection, all parameters have reasonable default values and do not need to be manually specified unless customization is needed. + +This shear box model uses an additional postprocessor, anisotropic stress, which is also implemented in this cookbook. It outputs a 3-by-3 matrix that can be visualized as a tensor, similar to the standard stress postprocessor. With the anisotropic viscosity material model, applying simple shear produces deformation in multiple directions. As a result, the anisotropic stress tensor appears as elongated and slightly tilted glyphs (indicating the principal stress directions), in contrast to the isotropic stress tensor (see figure below). + +```{figure-md} fig:anisotropic_stress_shearbox + + +Expected output of the shear box model using anisotropic viscosity material model, showing the anisotropic stress and stress postprocessor as tensor glyphs (blue disks) in Paraview. The arrows indicate the direction and magnitude of velocity. +``` \ No newline at end of file diff --git a/cookbooks/CPO_induced_anisotropic_viscosity/doc/anisotropic_stress.png b/cookbooks/CPO_induced_anisotropic_viscosity/doc/anisotropic_stress.png new file mode 100644 index 00000000000..3d0f2f8a1ed Binary files /dev/null and b/cookbooks/CPO_induced_anisotropic_viscosity/doc/anisotropic_stress.png differ diff --git a/cookbooks/CPO_induced_anisotropic_viscosity/doc/compositional_field.part.prm b/cookbooks/CPO_induced_anisotropic_viscosity/doc/compositional_field.part.prm new file mode 100644 index 00000000000..4ddadbf4bb0 --- /dev/null +++ b/cookbooks/CPO_induced_anisotropic_viscosity/doc/compositional_field.part.prm @@ -0,0 +1,32 @@ +subsection Initial composition model + set List of model names = function + subsection Function + set Function expression = 0;\ + 0; 0; 0; 0;\ + 0; 0; 0; 0;\ + 0; 0; 0; 0;\ + 0; 0 + end +end + +subsection Compositional fields + set Number of fields = 15 + set Names of fields = scalar_vis, \ + phi1, eigvalue_a1, eigvalue_a2, eigvalue_a3,\ + theta, eigvalue_b1, eigvalue_b2, eigvalue_b3,\ + phi2, eigvalue_c1, eigvalue_c2, eigvalue_c3,\ + D, water + set Compositional field methods = prescribed field, \ + particles, particles, particles, particles, \ + particles, particles, particles, particles, \ + particles, particles, particles, particles, \ + static, static + set Mapped particle properties = phi1:cpo mineral 0 phi1[0], eigvalue_a1:cpo mineral 0 eigenvalues a axis[0], eigvalue_a2:cpo mineral 0 eigenvalues a axis[1], eigvalue_a3:cpo mineral 0 eigenvalues a axis[2], \ + theta:cpo mineral 0 theta[0], eigvalue_b1:cpo mineral 0 eigenvalues b axis[0], eigvalue_b2:cpo mineral 0 eigenvalues b axis[1], eigvalue_b3:cpo mineral 0 eigenvalues b axis[2], \ + phi2:cpo mineral 0 phi2[0], eigvalue_c1:cpo mineral 0 eigenvalues c axis[0], eigvalue_c2:cpo mineral 0 eigenvalues c axis[1], eigvalue_c3:cpo mineral 0 eigenvalues c axis[2] + set Types of fields = generic, \ + generic, generic, generic, generic, \ + generic, generic, generic, generic, \ + generic, generic, generic, generic, \ + generic, generic +end diff --git a/cookbooks/CPO_induced_anisotropic_viscosity/doc/cpo_particle_property.part.prm b/cookbooks/CPO_induced_anisotropic_viscosity/doc/cpo_particle_property.part.prm new file mode 100644 index 00000000000..5b8ab6779a4 --- /dev/null +++ b/cookbooks/CPO_induced_anisotropic_viscosity/doc/cpo_particle_property.part.prm @@ -0,0 +1,60 @@ +subsection Postprocess + set List of postprocessors = velocity statistics, composition statistics, memory statistics, visualization, particles, crystal preferred orientation + + subsection Visualization + set Time between graphical output = 0.1 + set List of output variables = material properties, strain rate, named additional outputs, shear stress, stress + + subsection Material properties + set List of material properties = density, viscosity + end + end + subsection Particles + set Time between data output = 0.1 + set Data output format = gnuplot, vtu + set Exclude output properties = rotation_matrix, volume fraction + end + subsection Crystal Preferred Orientation + set Time between data output = 0.1 + set Write in background thread = true + set Compress cpo data files = false + set Write out raw cpo data = mineral 0: volume fraction, mineral 0: Euler angles #, mineral 1: volume fraction, mineral 1: Euler angles + set Write out draw volume weighted cpo data = mineral 0: volume fraction, mineral 0: Euler angles #, mineral 1: volume fraction, mineral 1: Euler angles + end +end + +subsection Particles + set List of particle properties = integrated strain, integrated strain invariant, velocity, pT path, strain rate, crystal preferred orientation, cpo bingham average, cpo elastic tensor #integrated strain, integrated strain invariant, velocity, pT path, strain rate, velocity gradient, stress, crystal preferred orientation, cpo bingham average, cpo bingham average euler angle, cpo elastic tensor + set Allow cells without particles = false + set Interpolation scheme = nearest neighbor + set Minimum particles per cell = 1 + set Maximum particles per cell = 5 + set Load balancing strategy = add particles + set Particle generator name = ascii file + subsection Generator + subsection Ascii file + set Data directory = ./ + set Data file name = particle_one.dat + end + end + subsection Crystal Preferred Orientation + subsection Initial grains + set Minerals = Olivine: A-fabric + set Volume fractions minerals = 1.0 + end + set Number of grains per particle = 1000 ## Annotation grain count + set Property advection method = Backward Euler ## Annotation Property advection method + set Property advection tolerance = 1e-15 + set CPO derivatives algorithm = D-Rex 2004 + subsection D-Rex 2004 + set Mobility = 10 ## Annotation LPO Mobility + set Stress exponents = 3.5 + set Exponents p = 1.5 + set Threshold GBS = 0.3 ## Annotation TGBS + end + end + subsection CPO Bingham Average + set Random number seed = 200 + set Use rotation matrix = false + end +end diff --git a/cookbooks/CPO_induced_anisotropic_viscosity/doc/min_particles_per_cell.part.prm b/cookbooks/CPO_induced_anisotropic_viscosity/doc/min_particles_per_cell.part.prm new file mode 100644 index 00000000000..9469d4c2aab --- /dev/null +++ b/cookbooks/CPO_induced_anisotropic_viscosity/doc/min_particles_per_cell.part.prm @@ -0,0 +1,2 @@ +set Minimum particles per cell = 1 +set Load balancing strategy = add particles diff --git a/cookbooks/CPO_induced_anisotropic_viscosity/particle_one.dat b/cookbooks/CPO_induced_anisotropic_viscosity/particle_one.dat new file mode 100644 index 00000000000..d7185572a95 --- /dev/null +++ b/cookbooks/CPO_induced_anisotropic_viscosity/particle_one.dat @@ -0,0 +1 @@ +0.5 0.5 0.5 diff --git a/cookbooks/CPO_induced_anisotropic_viscosity/plugin/CMakeLists.txt b/cookbooks/CPO_induced_anisotropic_viscosity/plugin/CMakeLists.txt new file mode 100644 index 00000000000..38bce63c4d5 --- /dev/null +++ b/cookbooks/CPO_induced_anisotropic_viscosity/plugin/CMakeLists.txt @@ -0,0 +1,39 @@ +# Copyright (C) 2011 - 2024 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 +# . + +CMAKE_MINIMUM_REQUIRED(VERSION 3.13.4) + +FIND_PACKAGE(Aspect 2.4.0 QUIET HINTS ${Aspect_DIR} ../ ../../ $ENV{ASPECT_DIR}) + +IF (NOT Aspect_FOUND) + MESSAGE(FATAL_ERROR "\n" + "Could not find a valid ASPECT build/installation directory. " + "Please specify the directory where you are building ASPECT by passing\n" + " -D Aspect_DIR=\n" + "to cmake or by setting the environment variable ASPECT_DIR in your shell " + "before calling cmake. See the section 'How to write a plugin' in the " + "manual for more information.") +ENDIF () + +DEAL_II_INITIALIZE_CACHED_VARIABLES() + +SET(TARGET "CPO_induced_anisotropic_viscosity") +PROJECT(${TARGET}) + +ADD_LIBRARY(${TARGET} SHARED cpo_induced_anisotropic_viscosity.h cpo_induced_anisotropic_viscosity.cc anisotropic_stress.cc anisotropic_stress.h) +ASPECT_SETUP_PLUGIN(${TARGET}) diff --git a/cookbooks/CPO_induced_anisotropic_viscosity/plugin/anisotropic_stress.cc b/cookbooks/CPO_induced_anisotropic_viscosity/plugin/anisotropic_stress.cc new file mode 100644 index 00000000000..2859ee76f8d --- /dev/null +++ b/cookbooks/CPO_induced_anisotropic_viscosity/plugin/anisotropic_stress.cc @@ -0,0 +1,127 @@ +/* + Copyright (C) 2011 - 2021 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 "anisotropic_stress.h" +#include "cpo_induced_anisotropic_viscosity.h" + + +namespace aspect +{ + namespace Postprocess + { + namespace VisualizationPostprocessors + { + template + AnisotropicStress:: + AnisotropicStress () + : + DataPostprocessorTensor ("anisotropic_stress", + update_values | update_gradients | update_quadrature_points), + Interface("Pa") + {} + + + + template + void + AnisotropicStress:: + evaluate_vector_field(const DataPostprocessorInputs::Vector &input_data, + std::vector> &computed_quantities) const + { + const unsigned int n_quadrature_points = input_data.solution_values.size(); + Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError()); + Assert ((computed_quantities[0].size() == Tensor<2,dim>::n_independent_components), + ExcInternalError()); + Assert (input_data.solution_values[0].size() == this->introspection().n_components, ExcInternalError()); + Assert (input_data.solution_gradients[0].size() == this->introspection().n_components, ExcInternalError()); + + MaterialModel::MaterialModelInputs in(input_data,this->introspection()); + MaterialModel::MaterialModelOutputs out(n_quadrature_points,this->n_compositional_fields()); + + this->get_material_model().create_additional_named_outputs(out); + this->get_material_model().evaluate(in, out); + + const std::shared_ptr> anisotropic_viscosity = out.template get_additional_output_object>(); + AssertThrow(anisotropic_viscosity != nullptr, + ExcMessage("Need anisotropic viscosity tensor from the anisotropic viscosity material model for computing the anisotropic stress.")); + + // ...and use it to compute the stresses + for (unsigned int q=0; q strain_rate = in.strain_rate[q]; + const SymmetricTensor<2,dim> deviatoric_strain_rate + = (this->get_material_model().is_compressible() + ? + strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor() + : + strain_rate); + + const double eta = out.viscosities[q]; + + SymmetricTensor<2,dim> aniso_stress; + aniso_stress= 2.*eta*deviatoric_strain_rate*anisotropic_viscosity->stress_strain_directors[q]; + + for (unsigned int d=0; d::component_to_unrolled_index(TableIndices<2>(d,e))] + = aniso_stress[d][e]; + } + + // average the values if requested + const auto &viz = this->get_postprocess_manager().template get_matching_active_plugin>(); + if (!viz.output_pointwise_stress_and_strain()) + average_quantities(computed_quantities); + } + + template + void + AnisotropicStress:: + create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &out) const + { + this->get_material_model().create_additional_named_outputs(out); + } + } + } +} + + +// explicit instantiations +namespace aspect +{ + namespace Postprocess + { + namespace VisualizationPostprocessors + { + ASPECT_REGISTER_VISUALIZATION_POSTPROCESSOR(AnisotropicStress, + "Anisotropic stress", + "A visualization output object that generates output " + "for the 6 (in 3d) components of the anisotropic stress " + "tensor. The anisotropic stress is defined as $2 \eta " + "(\varepsilon(\mathbf u) - \tfrac 13 \textrm{trace}\ " + "\varepsilon(\mathbf u) \mathbf 1) = 2\eta (\varepsilon(\mathbf u) - " + "\frac 13 (\nabla \cdot \mathbf u) \mathbf I)$ * stress_strain_directors, and differs from the " + "full stress by the absence of the pressure. The second term in the " + "difference is zero if the model is incompressible. " + "This particle property plugin should only be activated when using the " + "CPO induced anisotropic viscosity material model from the cookbook. ") + } + } +} diff --git a/cookbooks/CPO_induced_anisotropic_viscosity/plugin/anisotropic_stress.h b/cookbooks/CPO_induced_anisotropic_viscosity/plugin/anisotropic_stress.h new file mode 100644 index 00000000000..8cfac7d5849 --- /dev/null +++ b/cookbooks/CPO_induced_anisotropic_viscosity/plugin/anisotropic_stress.h @@ -0,0 +1,74 @@ +/* + Copyright (C) 2015 - 2021 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_postprocess_visualization_anisotropic_stress_h +#define _aspect_postprocess_visualization_anisotropic_stress_h + +#include +#include + +#include + + +namespace aspect +{ + namespace Postprocess + { + namespace VisualizationPostprocessors + { + /** + * A class derived from DataPostprocessor that takes an output vector + * and computes a variable that represents the 3 or 6 independent + * components (in 2d and 3d, respectively) of the anisotropic stress tensor at + * every point. The anisotropic stress is defined as $2 \eta + * (\varepsilon(\mathbf u) - \tfrac 13 \textrm{trace}\ + * \varepsilon(\mathbf u) \mathbf 1) = 2\eta (\varepsilon(\mathbf u) - + * \frac 13 (\nabla \cdot \mathbf u) \mathbf I)$ * stress_strain_directors, and differs from the + * full stress by the absence of the pressure. The second term in the + * difference is zero if the model is incompressible. + * + * The member functions are all implementations of those declared in the + * base class. See there for their meaning. + */ + template + class AnisotropicStress + : public DataPostprocessorTensor, + public SimulatorAccess, + public Interface + { + public: + /** + * Constructor. + */ + AnisotropicStress (); + + void + evaluate_vector_field(const DataPostprocessorInputs::Vector &input_data, + std::vector> &computed_quantities) const override; + + void + create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs &out) const; + }; + } + } +} + +#endif diff --git a/cookbooks/CPO_induced_anisotropic_viscosity/plugin/cpo_induced_anisotropic_viscosity.cc b/cookbooks/CPO_induced_anisotropic_viscosity/plugin/cpo_induced_anisotropic_viscosity.cc new file mode 100644 index 00000000000..651aa5cf28e --- /dev/null +++ b/cookbooks/CPO_induced_anisotropic_viscosity/plugin/cpo_induced_anisotropic_viscosity.cc @@ -0,0 +1,687 @@ + +/* + Copyright (C) 2015 - 2024 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 "cpo_induced_anisotropic_viscosity.h" +#include "../../cookbooks/anisotropic_viscosity/av_material.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +DEAL_II_DISABLE_EXTRA_DIAGNOSTICS +#include +DEAL_II_ENABLE_EXTRA_DIAGNOSTICS + + +namespace aspect +{ + namespace MaterialModel + { + using namespace dealii; + } + + 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); + } + } + +} + +namespace aspect +{ + +//Next session is a more evolved implementation of anisotropic viscosity in the material model based on Hansen et al 2016 and Kiraly et al 2020 + namespace MaterialModel + { + + template + void + CPO_AV_3D::set_assemblers(const SimulatorAccess &, + Assemblers::Manager &assemblers) const + { + for (unsigned int i=0; i>(*(assemblers.stokes_preconditioner[i]))) + assemblers.stokes_preconditioner[i] = std::make_unique> (); + } + + for (unsigned int i=0; i>(*(assemblers.stokes_system[i]))) + assemblers.stokes_system[i] = std::make_unique> (); + } + } + + + + template + void + CPO_AV_3D:: + initialize() + { + this->get_signals().set_assemblers.connect (std::bind(&CPO_AV_3D::set_assemblers, + std::cref(*this), + std::placeholders::_1, + std::placeholders::_2)); + AssertThrow((dim==3), + ExcMessage("Olivine has 3 independent slip systems, allowing for deformation in 3 independent directions, hence these models only work in 3D")); + + cpo_bingham_avg_a.push_back (this->introspection().compositional_index_for_name("phi1")); + cpo_bingham_avg_a.push_back (this->introspection().compositional_index_for_name("eigvalue_a1")); + cpo_bingham_avg_a.push_back (this->introspection().compositional_index_for_name("eigvalue_a2")); + cpo_bingham_avg_a.push_back (this->introspection().compositional_index_for_name("eigvalue_a3")); + + cpo_bingham_avg_b.push_back (this->introspection().compositional_index_for_name("theta")); + cpo_bingham_avg_b.push_back (this->introspection().compositional_index_for_name("eigvalue_b1")); + cpo_bingham_avg_b.push_back (this->introspection().compositional_index_for_name("eigvalue_b2")); + cpo_bingham_avg_b.push_back (this->introspection().compositional_index_for_name("eigvalue_b3")); + + cpo_bingham_avg_c.push_back (this->introspection().compositional_index_for_name("phi2")); + cpo_bingham_avg_c.push_back (this->introspection().compositional_index_for_name("eigvalue_c1")); + cpo_bingham_avg_c.push_back (this->introspection().compositional_index_for_name("eigvalue_c2")); + cpo_bingham_avg_c.push_back (this->introspection().compositional_index_for_name("eigvalue_c3")); + + + } + + + + template + Tensor<2,3> + AnisotropicViscosity::euler_angles_to_rotation_matrix(double phi1, double theta, double phi2) + { + Tensor<2,3> rot_matrix; + //R3*R2*R1 ZXZ rotation. Note it is not exactly the same as in utilities.cc + rot_matrix[0][0] = cos(phi2)*cos(phi1) - cos(theta)*sin(phi1)*sin(phi2); // + rot_matrix[0][1] = -cos(phi2)*sin(phi1) - cos(theta)*cos(phi1)*sin(phi2); //cos(phi2)*sin(phi1) + cos(theta)*cos(phi1)*sin(phi2); + rot_matrix[0][2] = sin(phi2)*sin(theta); + rot_matrix[1][0] = sin(phi2)*cos(phi1) + cos(theta)*sin(phi1)*cos(phi2); //-sin(phi2)*cos(phi1) - cos(theta)*sin(phi1)*cos(phi2); + rot_matrix[1][1] = -sin(phi2)*sin(phi1) + cos(theta)*cos(phi1)*cos(phi2); + rot_matrix[1][2] = -cos(phi2)*sin(theta); //cos(phi2)*sin(theta); + rot_matrix[2][0] = sin(theta)*sin(phi1); + rot_matrix[2][1] = sin(theta)*cos(phi1); //-sin(theta)*cos(phi1); + rot_matrix[2][2] = cos(theta); // + AssertThrow(rot_matrix[2][2] <= 1.0, ExcMessage("rot_matrix[2][2] > 1.0")); + return rot_matrix; + } + + + + template <> + void + CPO_AV_3D<2>::evaluate (const MaterialModel::MaterialModelInputs<2> &, + MaterialModel::MaterialModelOutputs<2> &) const + { + Assert (false, ExcNotImplemented()); + } + + + template <> + void + CPO_AV_3D<3>::evaluate (const MaterialModel::MaterialModelInputs<3> &in, + MaterialModel::MaterialModelOutputs<3> &out) const + { + const int dim=3; + const std::shared_ptr> anisotropic_viscosity = + out.template get_additional_output_object>(); + EquationOfStateOutputs eos_outputs (1); + + for (unsigned int q=0; q composition = in.composition[q]; + const std::vector volume_fractions = MaterialUtilities::compute_only_composition_fractions(composition, + this->introspection().chemical_composition_field_indices()); + + out.densities[q] = eos_outputs.densities[0]; + out.viscosities[q] = eta; + out.thermal_expansion_coefficients[q] = eos_outputs.thermal_expansion_coefficients[0]; + out.specific_heat[q] = eos_outputs.specific_heat_capacities[0]; + out.thermal_conductivities[q] = 1; + out.compressibilities[q] = eos_outputs.compressibilities[0]; + out.entropy_derivative_pressure[q] = eos_outputs.entropy_derivative_pressure[0]; + out.entropy_derivative_temperature[q] = eos_outputs.entropy_derivative_temperature[0]; + + // Calculate effective viscosity + const SymmetricTensor<2,dim> strain_rate = in.strain_rate[q]; + const SymmetricTensor<2,dim> deviatoric_strain_rate + = (this->get_material_model().is_compressible() + ? + strain_rate - 1./3. * trace(strain_rate) * unit_symmetric_tensor() + : + strain_rate); + + // The computation of the viscosity tensor is only necessary after the simulator has been initialized + // and when the condition allows dislocation creep + if ((this->simulator_is_past_initialization()) && (this->get_timestep_number() > 0) && (in.temperature[q]>1000) && (isfinite(determinant(deviatoric_strain_rate)))) + { + const unsigned int ind_vis = this->introspection().compositional_index_for_name("scalar_vis"); + + // Create constant value to use for AV + const double A_o = 1.1e5*exp(-530000/(8.314*in.temperature[q])); + const double n = 3.5; + const double Gamma = (A_o/(std::pow(grain_size,0.73)));// 3.5322e-15[1/(s*Pa^n)] if T=1600K and d=1000 microns + + // Get eigen values from compositional fields + const double eigvalue_a1 = composition[cpo_bingham_avg_a[1]]; + const double eigvalue_b1 = composition[cpo_bingham_avg_b[1]]; + const double eigvalue_c1 = composition[cpo_bingham_avg_c[1]]; + const double eigvalue_a2 = composition[cpo_bingham_avg_a[2]]; + const double eigvalue_b2 = composition[cpo_bingham_avg_b[2]]; + const double eigvalue_c2 = composition[cpo_bingham_avg_c[2]]; + const double eigvalue_a3 = composition[cpo_bingham_avg_a[3]]; + const double eigvalue_b3 = composition[cpo_bingham_avg_b[3]]; + const double eigvalue_c3 = composition[cpo_bingham_avg_c[3]]; + + // Calculate the rotation matrix from the euler angles + const double phi1 = composition[cpo_bingham_avg_a[0]]; + const double theta = composition[cpo_bingham_avg_b[0]]; + const double phi2 = composition[cpo_bingham_avg_c[0]]; + + Tensor<2,3> R = transpose(AnisotropicViscosity::euler_angles_to_rotation_matrix(phi1, theta, phi2)); + + // Compute Hill Parameters FGHLMN from the eigenvalues of a,b,c axis + double F, G, H, L, M, N; + // CPO2Hill v3 model (default coefficients are given for v5): + // F = std::pow(eigvalue_a1,2)*CnI_F[0] + eigvalue_a1*CnI_F[1] + eigvalue_a2*CnI_F[2] + (1/eigvalue_a3)*CnI_F[3] + std::pow(eigvalue_b1,2)*CnI_F[4] + eigvalue_b1*CnI_F[5] + eigvalue_b2*CnI_F[6] + (1/eigvalue_b3)*CnI_F[7] + std::pow(eigvalue_c1,2)*CnI_F[8] + eigvalue_c1*CnI_F[9] + eigvalue_c2*CnI_F[10] + (1/eigvalue_c3)*CnI_F[11] + CnI_F[12]; + // G = std::pow(eigvalue_a1,2)*CnI_G[0] + eigvalue_a1*CnI_G[1] + eigvalue_a2*CnI_G[2] + (1/eigvalue_a3)*CnI_G[3] + std::pow(eigvalue_b1,2)*CnI_G[4] + eigvalue_b1*CnI_G[5] + eigvalue_b2*CnI_G[6] + (1/eigvalue_b3)*CnI_G[7] + std::pow(eigvalue_c1,2)*CnI_G[8] + eigvalue_c1*CnI_G[9] + eigvalue_c2*CnI_G[10] + (1/eigvalue_c3)*CnI_G[11] + CnI_G[12]; + // H = std::pow(eigvalue_a1,2)*CnI_H[0] + eigvalue_a1*CnI_H[1] + eigvalue_a2*CnI_H[2] + (1/eigvalue_a3)*CnI_H[3] + std::pow(eigvalue_b1,2)*CnI_H[4] + eigvalue_b1*CnI_H[5] + eigvalue_b2*CnI_H[6] + (1/eigvalue_b3)*CnI_H[7] + std::pow(eigvalue_c1,2)*CnI_H[8] + eigvalue_c1*CnI_H[9] + eigvalue_c2*CnI_H[10] + (1/eigvalue_c3)*CnI_H[11] + CnI_H[12]; + // L = std::abs(std::pow(eigvalue_a1,2)*CnI_L[0] + eigvalue_a1*CnI_L[1] + eigvalue_a2*CnI_L[2] + (1/eigvalue_a3)*CnI_L[3] + std::pow(eigvalue_b1,2)*CnI_L[4] + eigvalue_b1*CnI_L[5] + eigvalue_b2*CnI_L[6] + (1/eigvalue_b3)*CnI_L[7] + std::pow(eigvalue_c1,2)*CnI_L[8] + eigvalue_c1*CnI_L[9] + eigvalue_c2*CnI_L[10] + (1/eigvalue_c3)*CnI_L[11] + CnI_L[12]); + // M = std::abs(std::pow(eigvalue_a1,2)*CnI_M[0] + eigvalue_a1*CnI_M[1] + eigvalue_a2*CnI_M[2] + (1/eigvalue_a3)*CnI_M[3] + std::pow(eigvalue_b1,2)*CnI_M[4] + eigvalue_b1*CnI_M[5] + eigvalue_b2*CnI_M[6] + (1/eigvalue_b3)*CnI_M[7] + std::pow(eigvalue_c1,2)*CnI_M[8] + eigvalue_c1*CnI_M[9] + eigvalue_c2*CnI_M[10] + (1/eigvalue_c3)*CnI_M[11] + CnI_M[12]); + // N = std::abs(std::pow(eigvalue_a1,2)*CnI_N[0] + eigvalue_a1*CnI_N[1] + eigvalue_a2*CnI_N[2] + (1/eigvalue_a3)*CnI_N[3] + std::pow(eigvalue_b1,2)*CnI_N[4] + eigvalue_b1*CnI_N[5] + eigvalue_b2*CnI_N[6] + (1/eigvalue_b3)*CnI_N[7] + std::pow(eigvalue_c1,2)*CnI_N[8] + eigvalue_c1*CnI_N[9] + eigvalue_c2*CnI_N[10] + (1/eigvalue_c3)*CnI_N[11] + CnI_N[12]); + + // CPO2Hill v5 model: + F = std::pow(eigvalue_a1,2)*CnI_F[0] + eigvalue_a2*CnI_F[1] + (1/eigvalue_a3)*CnI_F[2] + std::pow(eigvalue_b1,2)*CnI_F[3] + eigvalue_b2*CnI_F[4] + (1/eigvalue_b3)*CnI_F[5] + std::pow(eigvalue_c1,2)*CnI_F[6] + eigvalue_c2*CnI_F[7] + (1/eigvalue_c3)*CnI_F[8] + CnI_F[9]; + G = std::pow(eigvalue_a1,2)*CnI_G[0] + eigvalue_a2*CnI_G[1] + (1/eigvalue_a3)*CnI_G[2] + std::pow(eigvalue_b1,2)*CnI_G[3] + eigvalue_b2*CnI_G[4] + (1/eigvalue_b3)*CnI_G[5] + std::pow(eigvalue_c1,2)*CnI_G[6] + eigvalue_c2*CnI_G[7] + (1/eigvalue_c3)*CnI_G[8] + CnI_G[9]; + H = std::pow(eigvalue_a1,2)*CnI_H[0] + eigvalue_a2*CnI_H[1] + (1/eigvalue_a3)*CnI_H[2] + std::pow(eigvalue_b1,2)*CnI_H[3] + eigvalue_b2*CnI_H[4] + (1/eigvalue_b3)*CnI_H[5] + std::pow(eigvalue_c1,2)*CnI_H[6] + eigvalue_c2*CnI_H[7] + (1/eigvalue_c3)*CnI_H[8] + CnI_H[9]; + L = std::abs(std::pow(eigvalue_a1,2)*CnI_L[0] + eigvalue_a2*CnI_L[1] + (1/eigvalue_a3)*CnI_L[2] + std::pow(eigvalue_b1,2)*CnI_L[3] + eigvalue_b2*CnI_L[4] + (1/eigvalue_b3)*CnI_L[5] + std::pow(eigvalue_c1,2)*CnI_L[6] + eigvalue_c2*CnI_L[7] + (1/eigvalue_c3)*CnI_L[8] + CnI_L[9]); + M = std::abs(std::pow(eigvalue_a1,2)*CnI_M[0] + eigvalue_a2*CnI_M[1] + (1/eigvalue_a3)*CnI_M[2] + std::pow(eigvalue_b1,2)*CnI_M[3] + eigvalue_b2*CnI_M[4] + (1/eigvalue_b3)*CnI_M[5] + std::pow(eigvalue_c1,2)*CnI_M[6] + eigvalue_c2*CnI_M[7] + (1/eigvalue_c3)*CnI_M[8] + CnI_M[9]); + N = std::abs(std::pow(eigvalue_a1,2)*CnI_N[0] + eigvalue_a2*CnI_N[1] + (1/eigvalue_a3)*CnI_N[2] + std::pow(eigvalue_b1,2)*CnI_N[3] + eigvalue_b2*CnI_N[4] + (1/eigvalue_b3)*CnI_N[5] + std::pow(eigvalue_c1,2)*CnI_N[6] + eigvalue_c2*CnI_N[7] + (1/eigvalue_c3)*CnI_N[8] + CnI_N[9]); + + Tensor<2,6> R_CPO_K; + R_CPO_K[0][0] = std::pow(R[0][0],2); + R_CPO_K[0][1] = std::pow(R[0][1],2); + R_CPO_K[0][2] = std::pow(R[0][2],2); + R_CPO_K[0][3] = numbers::SQRT2*R[0][1]*R[0][2]; + R_CPO_K[0][4] = numbers::SQRT2*R[0][0]*R[0][2]; + R_CPO_K[0][5] = numbers::SQRT2*R[0][0]*R[0][1]; + + R_CPO_K[1][0] = std::pow(R[1][0],2); + R_CPO_K[1][1] = std::pow(R[1][1],2); + R_CPO_K[1][2] = std::pow(R[1][2],2); + R_CPO_K[1][3] = numbers::SQRT2*R[1][1]*R[1][2]; + R_CPO_K[1][4] = numbers::SQRT2*R[1][0]*R[1][2]; + R_CPO_K[1][5] = numbers::SQRT2*R[1][0]*R[1][1]; + + R_CPO_K[2][0] = std::pow(R[2][0],2); + R_CPO_K[2][1] = std::pow(R[2][1],2); + R_CPO_K[2][2] = std::pow(R[2][2],2); + R_CPO_K[2][3] = numbers::SQRT2*R[2][1]*R[2][2]; + R_CPO_K[2][4] = numbers::SQRT2*R[2][0]*R[2][2]; + R_CPO_K[2][5] = numbers::SQRT2*R[2][0]*R[2][1]; + + R_CPO_K[3][0] = numbers::SQRT2*R[1][0]*R[2][0]; + R_CPO_K[3][1] = numbers::SQRT2*R[1][1]*R[2][1]; + R_CPO_K[3][2] = numbers::SQRT2*R[1][2]*R[2][2]; + R_CPO_K[3][3] = R[1][1]*R[2][2]+R[1][2]*R[2][1]; + R_CPO_K[3][4] = R[1][0]*R[2][2]+R[1][2]*R[2][0]; + R_CPO_K[3][5] = R[1][0]*R[2][1]+R[1][1]*R[2][0]; + + R_CPO_K[4][0] = numbers::SQRT2*R[0][0]*R[2][0]; + R_CPO_K[4][1] = numbers::SQRT2*R[0][1]*R[2][1]; + R_CPO_K[4][2] = numbers::SQRT2*R[0][2]*R[2][2]; + R_CPO_K[4][3] = R[0][1]*R[2][2]+R[0][2]*R[2][1]; + R_CPO_K[4][4] = R[0][0]*R[2][2]+R[0][2]*R[2][0]; + R_CPO_K[4][5] = R[0][0]*R[2][1]+R[0][1]*R[2][0]; + + R_CPO_K[5][0] = numbers::SQRT2*R[0][0]*R[1][0]; + R_CPO_K[5][1] = numbers::SQRT2*R[0][1]*R[1][1]; + R_CPO_K[5][2] = numbers::SQRT2*R[0][2]*R[1][2]; + R_CPO_K[5][3] = R[0][1]*R[1][2]+R[0][2]*R[1][1]; + R_CPO_K[5][4] = R[0][0]*R[1][2]+R[0][2]*R[1][0]; + R_CPO_K[5][5] = R[0][0]*R[1][1]+R[0][1]*R[1][0]; + + SymmetricTensor<2,6> A; + A[0][0] = 2.0/3.0*(F+H); + A[0][1] = 2.0/3.0*(-F); + A[0][2] = 2.0/3.0*(-H); + A[1][1] = 2.0/3.0*(G+F); + A[1][2] = 2.0/3.0*(-G); + A[2][2] = 2.0/3.0*(H+G); + A[3][3] = 2.0/3.0*L; + A[4][4] = 2.0/3.0*M; + A[5][5] = 2.0/3.0*N; + + // Invert using ScaLAPACK in dealii + FullMatrix A_mat(6, 6); + for (unsigned int ai=0; ai<6; ++ai) + { + for (unsigned int aj=0; aj<6; ++aj) + { + A_mat(ai,aj) = A[ai][aj]; + } + } + + // A is the anisotropic tensor for the fluidity. We need its inverse, but it's not invertible due to singularity. + // Thus we use a pseudo inverse function imported from the scalapack package from deal.ii + const double ratio = 1e-8; + std::shared_ptr grid = std::make_shared(this->get_mpi_communicator(),6,6,4,4); + ScaLAPACKMatrix A_scalapack(6,6,grid,4,4); + A_scalapack = A_mat; + A_scalapack.pseudoinverse(ratio); + FullMatrix pinvA_mat(6,6); + A_scalapack.copy_to(pinvA_mat); + + SymmetricTensor<2,6> invA; + for (unsigned int ai=0; ai<6; ++ai) + { + for (unsigned int aj=0; aj<6; ++aj) + { + invA[ai][aj] = pinvA_mat(ai,aj); + } + } + + // Calculate the fluidity tensor in the CPO frame + Tensor<2,6> V = R_CPO_K * invA * transpose(R_CPO_K);//invA;// + + // Convert rank 2 viscosity tensor to rank 4 + FullMatrix V_mat(6,6); + for (unsigned int vi=0; vi<6; ++vi) + { + for (unsigned int vj=0; vj<6; ++vj) + { + V_mat[vi][vj] = V[vi][vj]; + } + } + SymmetricTensor<4,dim> V_r4; + dealii::Physics::Notation::Kelvin::to_tensor(V_mat, V_r4); + + if (anisotropic_viscosity != nullptr) + { + anisotropic_viscosity->stress_strain_directors[q] = V_r4; + } + + double scalar_viscosity = composition[ind_vis]; + + // In first time step using input viscosity can lead to convergence issue if the strainrate varies significantly within the model domain. + // Thus for the first timestep we calculate an initial viscosity based on the strain rate. Why not later: seems to cause unstable solution(?) + if (this->get_timestep_number() == 1) + { + const double edot_ii=std::max(std::sqrt(std::max(-second_invariant(deviator(strain_rate)), 0.)), + min_strain_rate); + scalar_viscosity= 1/Gamma * std::pow(edot_ii,((1. - n)/n)); + } + + double n_iterations = 1; + double max_iteration = 100; + double residual = scalar_viscosity; + double threshold = 0.0001*scalar_viscosity; + SymmetricTensor<2,dim> stress; + stress = 2*scalar_viscosity * V_r4 * deviatoric_strain_rate / 1e6; // Use stress in MPa + + while (std::abs(residual) > threshold && n_iterations < max_iteration) + { + stress = (1./2.) * (stress + 2*scalar_viscosity * V_r4 * deviatoric_strain_rate / 1e6); + + Tensor<2,dim> S_CPO=transpose(R)*stress*R; + + double Jhill = F*pow((S_CPO[0][0]-S_CPO[1][1]),2) + G*pow((S_CPO[1][1]-S_CPO[2][2]),2) + H*pow((S_CPO[2][2]-S_CPO[0][0]),2) + 2*L*pow(S_CPO[1][2],2) + 2*M*pow(S_CPO[0][2],2) + 2*N*pow(S_CPO[0][1],2); + if (Jhill < 0) + { + Jhill = std::abs(F)*pow((S_CPO[0][0]-S_CPO[1][1]),2) + std::abs(G)*pow((S_CPO[1][1]-S_CPO[2][2]),2) + std::abs(H)*pow((S_CPO[2][2]-S_CPO[0][0]),2) + 2*L*pow(S_CPO[1][2],2) + 2*M*pow(S_CPO[0][2],2) + 2*N*pow(S_CPO[0][1],2); + } + + AssertThrow(isfinite(Jhill), + ExcMessage("Jhill should be finite")); + AssertThrow(Jhill >= 0, + ExcMessage("Jhill should not be negative")); + + double scalar_viscosity_new = (1 / (Gamma * std::pow(Jhill,(n-1)/2))); // + residual = std::abs(scalar_viscosity_new - scalar_viscosity); + scalar_viscosity = scalar_viscosity_new; + threshold = 0.001*scalar_viscosity; + n_iterations += 1; + + } + // Overwrite the scalar viscosity with an effective viscosity + out.viscosities[q] = scalar_viscosity; + + AssertThrow(out.viscosities[q] > 0, + ExcMessage("Viscosity should be positive")); + AssertThrow(isfinite(out.viscosities[q]), + ExcMessage("Viscosity should be finite")); + + + } + else + { + if (anisotropic_viscosity != nullptr) + { + SymmetricTensor<2,6> V; + V[0][0] = 2.0/3.0; + V[0][1] = -1.0/3.0; + V[0][2] = -1.0/3.0; + V[1][1] = 2.0/3.0; + V[1][2] = -1.0/3.0; + V[2][2] = 2.0/3.0; + V[3][3] = 1; + V[4][4] = 1; + V[5][5] = 1; + + // Convert rank 2 viscosity tensor to rank 4 + FullMatrix V_mat(6,6); + for (unsigned int vi=0; vi<6; ++vi) + { + for (unsigned int vj=0; vj<6; ++vj) + { + V_mat[vi][vj] = V[vi][vj]; + } + } + SymmetricTensor<4,dim> V_r4; + dealii::Physics::Notation::Kelvin::to_tensor(V_mat, V_r4); + anisotropic_viscosity->stress_strain_directors[q] = V_r4; + + } + } + // Prescribe the stress strain directors and scalar viscosity to compositional field for access in the next time step + if (const std::shared_ptr> prescribed_field_out + = out.template get_additional_output_object>()) + { + const unsigned int ind_vis = this->introspection().compositional_index_for_name("scalar_vis"); + prescribed_field_out->prescribed_field_outputs[q][ind_vis] = out.viscosities[q]; + } + } + } + + + + template + bool + CPO_AV_3D::is_compressible () const + { + return false; + } + + + + template + void + CPO_AV_3D::declare_parameters (ParameterHandler &prm) + { + prm.enter_subsection("Material model"); + { + prm.enter_subsection("CPO-induced Anisotropic Viscosity"); + { + EquationOfState::LinearizedIncompressible::declare_parameters (prm); + + prm.declare_entry ("Coefficients and intercept for F", "1.0390459583037057, -0.767458622, 0.003066208, 0.19651133418307049, 0.413093763, 0.015463162, -0.935925291, -2.392877563, 0.051834768, 1.0799807050187482", + Patterns::List(Patterns::Double()), + "6 Coefficients and 1 intercept to compute the Hill Parameter F."); + prm.declare_entry ("Coefficients and intercept for G", "-2.836270315, -1.632453092, 0.000687606, 0.2671850239576621, -0.993392913, 0.002699241, 1.9689530759060374, 2.314442451425019, -0.018655905, 0.6887411607403755", + Patterns::List(Patterns::Double()), + "6 Coefficients and 1 intercept to compute the Hill Parameter G."); + prm.declare_entry ("Coefficients and intercept for H", "1.6687493021559732, 0.5797579293682223, 0.003241593, 0.701661336, 0.2513824481429968, 0.000229291, -2.003227619, -2.57032429, 0.071454541, 0.7490268673620638", + Patterns::List(Patterns::Double()), + "6 Coefficients and 1 intercept to compute the Hill Parameter H."); + prm.declare_entry ("Coefficients and intercept for L", "-0.325145943, 0.7284642859944138, 0.000404879, -0.665446098, 0.5152847961409479, 0.002722782, -1.026786493, -1.262574542, 0.009168498, 1.595422603", + Patterns::List(Patterns::Double()), + "6 Coefficients and 1 intercept to compute the Hill Parameter L."); + prm.declare_entry ("Coefficients and intercept for M", "1.6427437063774875, 0.8777500120437522, 0.004651732, 2.489417876177839, 0.8162729707609052, -0.010736521, -2.49420455, -0.511446494, -0.009362491, 0.893677343", + Patterns::List(Patterns::Double()), + "6 Coefficients and 1 intercept to compute the Hill Parameter M."); + prm.declare_entry ("Coefficients and intercept for N", "0.8122098589701904, 0.15663795996228266, 0.001500252, -1.648578168, 0.19362392490527092, -0.009650519, 1.6796559729985163, -0.103640482, 0.01971017, 1.2132200780065174", + Patterns::List(Patterns::Double()), + "6 Coefficients and 1 intercept to compute the Hill Parameter N."); + + prm.declare_entry ("Reference viscosity", "1e9", + Patterns::Double(), + "Magnitude of reference viscosity."); + prm.declare_entry ("Minimum strain rate", "1.4e-20", Patterns::Double(), + "Stabilizes strain dependent viscosity. Units: \\si{\\per\\second}"); + prm.declare_entry ("Grain size", "1000", + Patterns::Double(), + "Olivine anisotropic viscosity is dependent of grain size. Value is given in microns"); + prm.declare_entry ("Density differential for compositional field 1", "0.", + Patterns::Double(), + "If compositional fields are used, then one would frequently want " + "to make the density depend on these fields. In this simple material " + "model, we make the following assumptions: if no compositional fields " + "are used in the current simulation, then the density is simply the usual " + "one with its linear dependence on the temperature. If there are compositional " + "fields, then the density only depends on the first one in such a way that " + "the density has an additional term of the kind $+\\Delta \\rho \\; c_1(\\mathbf x)$. " + "This parameter describes the value of $\\Delta \\rho$. " + "Units: \\si{\\kilogram\\per\\meter\\cubed}/unit change in composition."); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + + } + + + + template + void + CPO_AV_3D::parse_parameters (ParameterHandler &prm) + { + prm.enter_subsection("Material model"); + { + prm.enter_subsection("CPO-induced Anisotropic Viscosity"); + { + + equation_of_state.parse_parameters (prm); + eta = prm.get_double("Reference viscosity"); + min_strain_rate = prm.get_double("Minimum strain rate"); + grain_size = prm.get_double("Grain size"); + CnI_F = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for F"))); + CnI_G = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for G"))); + CnI_H = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for H"))); + CnI_L = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for L"))); + CnI_M = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for M"))); + CnI_N = dealii::Utilities::string_to_double(dealii::Utilities::split_string_list(prm.get("Coefficients and intercept for N"))); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + + + + // Declare dependence + this->model_dependence.density = NonlinearDependence::compositional_fields; + } + + + + template + void + CPO_AV_3D::create_additional_named_outputs(MaterialModel::MaterialModelOutputs &out) const + { + if (out.template get_additional_output_object>() == nullptr) + { + const unsigned int n_points = out.n_evaluation_points(); + out.additional_outputs.push_back( + std::make_unique> (n_points)); + } + + if (out.template get_additional_output_object>() == NULL) + { + const unsigned int n_points = out.n_evaluation_points(); + out.additional_outputs.push_back( + std::make_unique> (n_points,this->n_compositional_fields())); + } + } + } +} + + + +// explicit instantiations +namespace aspect +{ + namespace Assemblers + { +#define INSTANTIATE(dim) \ + template class StokesPreconditionerAnisotropicViscosity; \ + template class StokesIncompressibleTermsAnisotropicViscosity; \ + + ASPECT_INSTANTIATE(INSTANTIATE) + } + + namespace HeatingModel + { + ASPECT_REGISTER_HEATING_MODEL(ShearHeatingAnisotropicViscosity, + "anisotropic shear heating for cpo_induced_anisotropic_viscosity", + "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(CPO_AV_3D, + "CPO-induced Anisotropic Viscosity", + "Olivine CPO related viscous anisotropy based on the Simple material model") + } +} diff --git a/cookbooks/CPO_induced_anisotropic_viscosity/plugin/cpo_induced_anisotropic_viscosity.h b/cookbooks/CPO_induced_anisotropic_viscosity/plugin/cpo_induced_anisotropic_viscosity.h new file mode 100644 index 00000000000..7a8649b1904 --- /dev/null +++ b/cookbooks/CPO_induced_anisotropic_viscosity/plugin/cpo_induced_anisotropic_viscosity.h @@ -0,0 +1,135 @@ +/* + Copyright (C) 2019 - 2020 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_material_model_CPO_AV_3D_h +#define _aspect_material_model_CPO_AV_3D_h + +#include +#include +#include +#include +#include + +namespace aspect +{ + namespace MaterialModel + { + using namespace dealii; + + /** + * A material model that implements the micromechanical behaviour of olivine grains to create anisotropic viscosity. + * Based on Hansen et al., 2016 (JGR) and Kiraly et al., 2020 (G3). + * The micromechanical model requires the euler angles of the olivine grains (now stored on 3 compositional field), + * the grainsize, temperature, and strain rate to calculate the stress that is needed to create the input strain rate. + * The material model otherwise is based on the Simple material model. + * @ingroup MaterialModels + */ + 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$ /i.e. effective viscosity), 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; + + static Tensor<2,3> euler_angles_to_rotation_matrix(double phi1, double theta, double phi2); + + + + }; + + + + template + class CPO_AV_3D : public MaterialModel::Simple + { + public: + void initialize() override; + + void evaluate (const MaterialModel::MaterialModelInputs &in, + MaterialModel::MaterialModelOutputs &out) const override; + + bool is_compressible () const override; + + static void declare_parameters (ParameterHandler &prm); + + void parse_parameters (ParameterHandler &prm) override; + + void create_additional_named_outputs(MaterialModel::MaterialModelOutputs &out) const override; + + private: + + /** + * Reference viscosity. + */ + double eta; + + /** + * Defining a minimum strain rate stabilizes the viscosity calculation, + * which involves a division by the strain rate. Units: 1/s. + */ + double min_strain_rate; + + /** + * These are arrays that store eigenvalues of the olivine textures in a, b, and c axis, and + * the olivine texture represented by Euler angles. For more details, please refer to the + * cpo_bingham_average particle property. To use the anisotropic viscosity plugin in this + * cookbook, the CPO Bingham Average particle property must be included and Use rotation matrix + * must be set to false. The resulting arrays are: + * cpo_bingham_avg_a = [phi, eigenvalue 1 for a-axis, eigenvalue 2 for a-axis, eigenvalue 3 for a-axis] + * cpo_bingham_avg_b = [phi, eigenvalue 1 for b-axis, eigenvalue 2 for b-axis, eigenvalue 3 for b-axis] + * cpo_bingham_avg_c = [phi, eigenvalue 1 for c-axis, eigenvalue 2 for c-axis, eigenvalue 3 for c-axis] + * They are used in computing rotation matrix with regards to the CPO reference frame, and + * the anisotropic Hill coefficients FGHLMN. + */ + std::vector cpo_bingham_avg_a, cpo_bingham_avg_b, cpo_bingham_avg_c; + + /** + * These are arrays that store coefficients used to compute the anisotropic Hill coefficients FGHLMN from + * a certain olivine texture represented with the eigenvalues of its a-, b-, and c-axis. Each array contains + * 9 coefficients and 1 constant. + */ + std::vector CnI_F, CnI_G, CnI_H, CnI_L, CnI_M, CnI_N; + + double grain_size; + + EquationOfState::LinearizedIncompressible equation_of_state; + + void set_assemblers(const SimulatorAccess &, + Assemblers::Manager &assemblers) const; + + }; + } +} + +#endif diff --git a/cookbooks/CPO_induced_anisotropic_viscosity/shearbox_cpo_av.prm b/cookbooks/CPO_induced_anisotropic_viscosity/shearbox_cpo_av.prm new file mode 100644 index 00000000000..1a7aa1d4f98 --- /dev/null +++ b/cookbooks/CPO_induced_anisotropic_viscosity/shearbox_cpo_av.prm @@ -0,0 +1,191 @@ +# This cookbook demonstrates how CPO induced anisotropic viscosity evolves during +# simple shearing. The model consists of 1 cell and one particle in the middle. +# The particle includes 500 grains with initially random orientations (isotropic +# texture), which then evolves into a CPO aligned with the shear direction. + +set Additional shared libraries = ./plugin/libCPO_induced_anisotropic_viscosity.debug.so +set CFL number = 0.1 +set Timing output frequency = 10000 +set Dimension = 3 +set Pressure normalization = surface +set Surface pressure = 0 +set Nonlinear solver scheme = single Advection, single Stokes +set End time = 0.1 +set Use years instead of seconds = false +set Output directory = output_Shearbox_CPO_AV + +subsection Initial composition model + set List of model names = function + + subsection Function + set Function expression = 0;\ + 0; 0; 0; 0;\ + 0; 0; 0; 0;\ + 0; 0; 0; 0;\ + 0 + end +end + +subsection Compositional fields + set Number of fields = 14 + set Names of fields = scalar_vis, \ + phi1, eigvalue_a1, eigvalue_a2, eigvalue_a3,\ + theta, eigvalue_b1, eigvalue_b2, eigvalue_b3,\ + phi2, eigvalue_c1, eigvalue_c2, eigvalue_c3,\ + water + set Compositional field methods = prescribed field, \ + particles, particles, particles, particles, \ + particles, particles, particles, particles, \ + particles, particles, particles, particles, \ + static + set Mapped particle properties = phi1:cpo mineral 0 phi1[0], eigvalue_a1:cpo mineral 0 eigenvalues a axis[0], eigvalue_a2:cpo mineral 0 eigenvalues a axis[1], eigvalue_a3:cpo mineral 0 eigenvalues a axis[2], \ + theta:cpo mineral 0 theta[0], eigvalue_b1:cpo mineral 0 eigenvalues b axis[0], eigvalue_b2:cpo mineral 0 eigenvalues b axis[1], eigvalue_b3:cpo mineral 0 eigenvalues b axis[2], \ + phi2:cpo mineral 0 phi2[0], eigvalue_c1:cpo mineral 0 eigenvalues c axis[0], eigvalue_c2:cpo mineral 0 eigenvalues c axis[1], eigvalue_c3:cpo mineral 0 eigenvalues c axis[2] + set Types of fields = generic, \ + generic, generic, generic, generic, \ + generic, generic, generic, generic, \ + generic, generic, generic, generic, \ + generic +end + +subsection Gravity model + set Model name = vertical + + subsection Vertical + set Magnitude = 0 + end +end + +subsection Geometry model + set Model name = box + + subsection Box + set X extent = 1 + set Y extent = 1 + set Z extent = 1 + set Box origin X coordinate = 0 + set Box origin Y coordinate = 0 + set Box origin Z coordinate = 0 + end +end + +subsection Initial temperature model + set Model name = function + + subsection Function + # Temperature function (anisotropic viscosity in this material is disabled at T<1000K) + set Function expression = 1600 + end +end + +subsection Boundary temperature model + set List of model names = box + set Fixed temperature boundary indicators = bottom, top, front, back, left, right + + subsection Box + set Bottom temperature = 1600 + set Top temperature = 1600 + set Front temperature = 1600 + set Back temperature = 1600 + set Left temperature = 1600 + set Right temperature = 1600 + end +end + +subsection Boundary velocity model + set Prescribed velocity boundary indicators = top:function, bottom:function, left:function, right:function, front:function, back:function + + subsection Function + set Variable names = x,y,z,t + set Function expression = (z-0.5);0;0 + end +end + + +subsection Material model + set Model name = CPO-induced Anisotropic Viscosity + + subsection CPO-induced Anisotropic Viscosity + set Reference viscosity = 1e9 + end +end + +subsection Mesh refinement + set Initial global refinement = 0 + set Initial adaptive refinement = 0 + # For this simple shear box test, there is no need for AMR + set Time steps between mesh refinement = 0 +end + +subsection Postprocess + set List of postprocessors = velocity statistics, composition statistics, memory statistics, visualization, particles, crystal preferred orientation + + subsection Visualization + set Time between graphical output = 0.1 + set List of output variables = material properties, strain rate, named additional outputs, shear stress, stress + + subsection Material properties + set List of material properties = density, viscosity + end +end + + subsection Particles + set Time between data output = 0.1 + set Data output format = gnuplot, vtu + set Exclude output properties = rotation_matrix, volume fraction + end + + subsection Crystal Preferred Orientation + set Time between data output = 0.1 + set Write in background thread = true + set Compress cpo data files = false + set Write out raw cpo data = mineral 0: volume fraction, mineral 0: Euler angles + set Write out draw volume weighted cpo data = mineral 0: volume fraction, mineral 0: Euler angles + end +end + +subsection Particles + set List of particle properties = integrated strain, integrated strain invariant, velocity, pT path, strain rate, crystal preferred orientation, cpo bingham average, cpo elastic tensor + set Allow cells without particles = false + set Interpolation scheme = nearest neighbor + set Particle generator name = ascii file + + subsection Generator + subsection Ascii file + set Data directory = ./ + set Data file name = particle_one.dat + end + end + + subsection Crystal Preferred Orientation + + subsection Initial grains + set Minerals = Olivine: A-fabric + set Volume fractions minerals = 1.0 + end + + set Number of grains per particle = 1000 + set Property advection method = Backward Euler + set Property advection tolerance = 1e-15 + set CPO derivatives algorithm = D-Rex 2004 + + subsection D-Rex 2004 + set Mobility = 10 + set Stress exponents = 3.5 + set Exponents p = 1.5 + set Threshold GBS = 0.3 + end + end + + subsection CPO Bingham Average + set Random number seed = 200 + set Use rotation matrix = false + end +end + +subsection Solver parameters + # set Temperature solver tolerance = 1e-10 + subsection Stokes solver parameters + set Stokes solver type = block AMG + end +end diff --git a/cookbooks/anisotropic_viscosity/CMakeLists.txt b/cookbooks/anisotropic_viscosity/CMakeLists.txt index 1c64d3d7194..83162e9149c 100644 --- a/cookbooks/anisotropic_viscosity/CMakeLists.txt +++ b/cookbooks/anisotropic_viscosity/CMakeLists.txt @@ -35,5 +35,5 @@ DEAL_II_INITIALIZE_CACHED_VARIABLES() set(TARGET "anisotropic_viscosity") project(${TARGET}) -add_library(${TARGET} SHARED av_material.cc) +add_library(${TARGET} SHARED av_material.cc av_material.h) ASPECT_SETUP_PLUGIN(${TARGET}) diff --git a/cookbooks/anisotropic_viscosity/av_material.cc b/cookbooks/anisotropic_viscosity/av_material.cc index 4a548ae8bcd..2a1698ede8b 100644 --- a/cookbooks/anisotropic_viscosity/av_material.cc +++ b/cookbooks/anisotropic_viscosity/av_material.cc @@ -18,6 +18,7 @@ . */ +#include "av_material.h" #include #include #include diff --git a/cookbooks/anisotropic_viscosity/av_material.h b/cookbooks/anisotropic_viscosity/av_material.h new file mode 100644 index 00000000000..a7d06c7b2ab --- /dev/null +++ b/cookbooks/anisotropic_viscosity/av_material.h @@ -0,0 +1,50 @@ +/* + Copyright (C) 2015 - 2024 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_material_model_AV_h +#define _aspect_material_model_AV_h + +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; + }; + } +} + + +#endif diff --git a/doc/modules/changes/20250827_Wang-yijun b/doc/modules/changes/20250827_Wang-yijun new file mode 100644 index 00000000000..c71de5ff97c --- /dev/null +++ b/doc/modules/changes/20250827_Wang-yijun @@ -0,0 +1,3 @@ +Added: A cookbook for CPO-induced anisotropic viscosity. +
+(Yijun Wang, 2025/08/27) diff --git a/tests/cpo_induced_anisotropic_viscosity.cc b/tests/cpo_induced_anisotropic_viscosity.cc new file mode 100644 index 00000000000..a901644740c --- /dev/null +++ b/tests/cpo_induced_anisotropic_viscosity.cc @@ -0,0 +1,24 @@ +/* + 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 "../cookbooks/CPO_induced_anisotropic_viscosity/plugin/cpo_induced_anisotropic_viscosity.cc" +#include "../cookbooks/CPO_induced_anisotropic_viscosity/plugin/anisotropic_stress.cc" +#include "../cookbooks/CPO_induced_anisotropic_viscosity/plugin/cpo_induced_anisotropic_viscosity.h" +#include "../cookbooks/CPO_induced_anisotropic_viscosity/plugin/anisotropic_stress.h" diff --git a/tests/cpo_induced_anisotropic_viscosity.prm b/tests/cpo_induced_anisotropic_viscosity.prm new file mode 100644 index 00000000000..c388f470681 --- /dev/null +++ b/tests/cpo_induced_anisotropic_viscosity.prm @@ -0,0 +1,16 @@ +include $ASPECT_SOURCE_DIR/cookbooks/CPO_induced_anisotropic_viscosity/shearbox_cpo_av.prm +set Additional shared libraries = ./libCPO_induced_anisotropic_viscosity.debug.so + +set End time = 0 + +subsection Particles + subsection Generator + subsection Ascii file + set Data directory = $ASPECT_SOURCE_DIR/cookbooks/CPO_induced_anisotropic_viscosity/ + set Data file name = particle_one.dat + end + end + subsection Crystal Preferred Orientation + set Number of grains per particle = 10 ## Annotation grain count + end +end diff --git a/tests/cpo_induced_anisotropic_viscosity/screen-output b/tests/cpo_induced_anisotropic_viscosity/screen-output new file mode 100644 index 00000000000..d3419b2bd3d --- /dev/null +++ b/tests/cpo_induced_anisotropic_viscosity/screen-output @@ -0,0 +1,35 @@ +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- + +Loading shared library <./libcpo_induced_anisotropic_viscosity.debug.so> + +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- +Number of active cells: 1 (on 1 levels) +Number of degrees of freedom: 494 (81+8+27+27+27+27+27+27+27+27+27+27+27+27+27+27+27) + +*** Timestep 0: t=0 seconds, dt=0 seconds + Solving temperature system... 0 iterations. + Advecting particles... done. + Copying properties into prescribed compositional field scalar_vis... done. + Rebuilding Stokes preconditioner... + Solving Stokes system (AMG)... 6+0 iterations. + + Postprocessing: + RMS, max velocity: 0.289 m/s, 0.5 m/s + Compositions min/max/mass: 1e+09/1e+09/1e+09 // 5.276/5.276/5.276 // 0.5203/0.5203/0.5203 // 0.3479/0.3479/0.3479 // 0.1318/0.1318/0.1318 // 1.736/1.736/1.736 // 0.5743/0.5743/0.5743 // 0.3808/0.3808/0.3808 // 0.04489/0.04489/0.04489 // 5.435/5.435/5.435 // 0.6079/0.6079/0.6079 // 0.3106/0.3106/0.3106 // 0.08157/0.08157/0.08157 // 0/0/0 + System matrix memory consumption: 0.36 MB + Writing graphical output: output-cpo_induced_anisotropic_viscosity/solution/solution-00000 + Writing particle output: output-cpo_induced_anisotropic_viscosity/particles/particles-00000 + Writing particle cpo output: output-cpo_induced_anisotropic_viscosity/particles_cpo/CPO-00000 + +Termination requested by criterion: end time + + ++------------------------------------------------+------------+------------+ ++------------------------------------+-----------+------------+------------+ ++------------------------------------+-----------+------------+------------+ + +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- diff --git a/tests/cpo_induced_anisotropic_viscosity/statistics b/tests/cpo_induced_anisotropic_viscosity/statistics new file mode 100644 index 00000000000..dc8a028f0e7 --- /dev/null +++ b/tests/cpo_induced_anisotropic_viscosity/statistics @@ -0,0 +1,69 @@ +# 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 14 +# 11: Iterations for Stokes solver +# 12: Velocity iterations in Stokes preconditioner +# 13: Schur complement iterations in Stokes preconditioner +# 14: RMS velocity (m/s) +# 15: Max. velocity (m/s) +# 16: Minimal value for composition scalar_vis +# 17: Maximal value for composition scalar_vis +# 18: Global mass for composition scalar_vis +# 19: Minimal value for composition phi1 +# 20: Maximal value for composition phi1 +# 21: Global mass for composition phi1 +# 22: Minimal value for composition eigvalue_a1 +# 23: Maximal value for composition eigvalue_a1 +# 24: Global mass for composition eigvalue_a1 +# 25: Minimal value for composition eigvalue_a2 +# 26: Maximal value for composition eigvalue_a2 +# 27: Global mass for composition eigvalue_a2 +# 28: Minimal value for composition eigvalue_a3 +# 29: Maximal value for composition eigvalue_a3 +# 30: Global mass for composition eigvalue_a3 +# 31: Minimal value for composition theta +# 32: Maximal value for composition theta +# 33: Global mass for composition theta +# 34: Minimal value for composition eigvalue_b1 +# 35: Maximal value for composition eigvalue_b1 +# 36: Global mass for composition eigvalue_b1 +# 37: Minimal value for composition eigvalue_b2 +# 38: Maximal value for composition eigvalue_b2 +# 39: Global mass for composition eigvalue_b2 +# 40: Minimal value for composition eigvalue_b3 +# 41: Maximal value for composition eigvalue_b3 +# 42: Global mass for composition eigvalue_b3 +# 43: Minimal value for composition phi2 +# 44: Maximal value for composition phi2 +# 45: Global mass for composition phi2 +# 46: Minimal value for composition eigvalue_c1 +# 47: Maximal value for composition eigvalue_c1 +# 48: Global mass for composition eigvalue_c1 +# 49: Minimal value for composition eigvalue_c2 +# 50: Maximal value for composition eigvalue_c2 +# 51: Global mass for composition eigvalue_c2 +# 52: Minimal value for composition eigvalue_c3 +# 53: Maximal value for composition eigvalue_c3 +# 54: Global mass for composition eigvalue_c3 +# 55: Minimal value for composition water +# 56: Maximal value for composition water +# 57: Global mass for composition water +# 58: System matrix memory consumption (MB) +# 59: Triangulation memory consumption (MB) +# 60: p4est memory consumption (MB) +# 61: DoFHandler memory consumption (MB) +# 62: AffineConstraints memory consumption (MB) +# 63: Solution vector memory consumption (MB) +# 64: Peak virtual memory usage (VmPeak) (MB) +# 65: Visualization file name +# 66: Number of advected particles +# 67: Particle file name +# 68: Particle CPO file name +0 0.000000000000e+00 0.000000000000e+00 1 89 27 378 0 0 0 6 7 7 2.88675135e-01 5.00000000e-01 1.00000000e+09 1.00000000e+09 1.00000000e+09 5.27630079e+00 5.27630079e+00 5.27630079e+00 5.20309402e-01 5.20309402e-01 5.20309402e-01 3.47923711e-01 3.47923711e-01 3.47923711e-01 1.31766887e-01 1.31766887e-01 1.31766887e-01 1.73622715e+00 1.73622715e+00 1.73622715e+00 5.74325744e-01 5.74325744e-01 5.74325744e-01 3.80780388e-01 3.80780388e-01 3.80780388e-01 4.48938682e-02 4.48938682e-02 4.48938682e-02 5.43504380e+00 5.43504380e+00 5.43504380e+00 6.07858069e-01 6.07858069e-01 6.07858069e-01 3.10569263e-01 3.10569263e-01 3.10569263e-01 8.15726682e-02 8.15726682e-02 8.15726682e-02 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.3605 0.0082 0.0050 0.0032 0.0040 0.0091 9507.7070 output-cpo_induced_anisotropic_viscosity/solution/solution-00000 1 output-cpo_induced_anisotropic_viscosity/particles/particles-00000 output-cpo_induced_anisotropic_viscosity/particles_cpo/CPO-00000