From: Jean-Paul Pelteret Date: Wed, 27 May 2020 19:47:04 +0000 (+0200) Subject: Add step-71 X-Git-Tag: v9.3.0-rc1~124^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0fa7eeea7cfb1149679df8cd5b3ab3fb3eb37c72;p=dealii.git Add step-71 --- diff --git a/doc/doxygen/headers/automatic_and_symbolic_differentiation.h b/doc/doxygen/headers/automatic_and_symbolic_differentiation.h index 06128e5c0c..2b003cafec 100644 --- a/doc/doxygen/headers/automatic_and_symbolic_differentiation.h +++ b/doc/doxygen/headers/automatic_and_symbolic_differentiation.h @@ -491,7 +491,7 @@ * - Classes designed to operate at the quadrature point level (or any general continuum point): * - Differentiation::AD::ScalarFunction: %Differentiation of a scalar-valued function. * One typical use would be the the development of constitutive laws directly from a strain - * energy function. + * energy function. An example of this exact use case is given in step-71. * - Differentiation::AD::VectorFunction: %Differentiation of a vector-valued function. * This could be used to linearize the kinematic variables of a constitutive law, or assist * in solving the evolution equations of local internal variables. @@ -591,6 +591,7 @@ * expressions using methods such as common subexpression elimination (CSE), as well as by generating * high performance code-paths to evaluate these expressions through the use of a custom-generated * `std::function` or by compiling the expression using the LLVM JIT compiler. + * The usage of the Differentiation::SD::BatchOptimizer class is exemplified in step-71. * * As a final note, it is important to recognize the remaining major deficiencies in deal.II's current * implementation of the interface to the supported symbolic library. @@ -598,6 +599,7 @@ * symbolic algebra to the traditional use case (i.e. scalar and tensor algebra, as might be useful to * define constitutive relations or complex functions for application as boundary conditions or * source terms). + * In fact, step-71 demonstrates how it can be used to implement challenging constitutive models. * In the future we will also implement classes to assist in performing assembly operations in * the same spirit as that which has been done in the Differentiation::AD namespace. * diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index c4a9a2fba0..4af204cd8e 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -877,6 +877,155 @@ year = {2008}, } +% ------------------------------------ +% Step 71 +% ------------------------------------ + + +@book{Pelteret2019a, + author = {Pelteret, J-P. V. and Steinmann, P.}, + publisher = {De Gruyter Mouton}, + title = {Magneto-active polymers: Fabrication, characterisation, modelling and simulation at the micro- and macro-scale}, + year = {2019}, + edition = {1}, + isbn = {9783110419511}, + comment = {In collaboration with: B. Brands, G. Chatzigeorgiou, D. Davydov, M. Hossain, A. Javili, D. Pivovarov, P. Saxena, F. Vogel, D.K. Vu, B. Walter and R. Zabihyan}, + doi = {10.1515/9783110418576}, + url = {https://www.degruyter.com/document/doi/10.1515/9783110418576/html}, +} + +@book{Logg2012a, + title = {Automated Solution of Differential Equations by the Finite Element Method}, + publisher = {Springer Berlin Heidelberg}, + year = {2012}, + editor = {A. Logg and K.-A. Mardal and G. Wells}, + isbn = {978-3-642-23098-1}, + doi = {10.1007/978-3-642-23099-8}, +} + +@book{Korelc2016a, + title = {Automation of Finite Element Methods}, + publisher = {Springer Nature}, + year = {2016}, + author = {J. Korelc and P. Wriggers}, + isbn = {978-3-319-39003-1}, + doi = {10.1007/978-3-319-39005-5}, +} + +@inbook{Truesdell1960a, + chapter = {2: The classical field theories}, + pages = {226--794}, + title = {Encyclopedia of Physics: Principles of Thermodynamics and Statics}, + publisher = {Springer-Verlag Berlin Heidelberg}, + year = {1960}, + author = {Truesdell, C. and Toupin, R.}, + editor = {Fl\"{u}gge, S.}, + volume = {1}, + isbn = {978-3-540-02547-4}, + doi = {10.1007/978-3-642-45943-6}, +} + +@Article{Coleman1963a, + author = {Coleman, B. D. and Noll, W.}, + title = {The thermodynamics of elastic materials with heat conduction and viscosity}, + journal = {Archive for Rational Mechanics and Analysis}, + year = {1963}, + volume = {13}, + number = {1}, + pages = {167--178}, + month = dec, + doi = {10.1007/BF01262690}, + publisher = {Springer}, +} + +@Article{Coleman1967a, + author = {Coleman, B. D. and Gurtin, M. E.}, + title = {Thermodynamics with internal state variables}, + year = {1967}, + volume = {47}, + number = {2}, + pages = {597--613}, + doi = {10.1063/1.1711937}, + journal = {The Journal of Chemical Physics}, + publisher = {AIP Publishing}, +} + +@InCollection{Pao1978a, + author = {Pao, Y. H.}, + title = {Electromagnetic Forces in Deformable Continua}, + booktitle = {Mechanics Today}, + publisher = {Elsevier}, + year = {1978}, + editor = {Nemat-Nasser, S.}, + volume = {4}, + series = {Pergamon Mechanics Today Series}, + chapter = {IV}, + pages = {209--305}, + address = {New York}, + isbn = {978-0-08-021792-5}, + doi = {10.1016/b978-0-08-021792-5.50012-4}, +} + +@Article{Holzapfel1996a, + author = {Holzapfel, G. A. and Simo, J. C.}, + title = {A new viscoelastic constitutive model for continuous media at finite thermomechanical changes}, + journal = {International Journal of Solids and Structures}, + year = {1996}, + volume = {33}, + number = {20--22}, + pages = {3019--3034}, + month = aug, + doi = {10.1016/0020-7683(95)00263-4}, + publisher = {Elsevier}, +} + +@Book{Holzapfel2007a, + title = {Nonlinear Solid Mechanics: A Continuum Approach for Engineering}, + publisher = {John Wiley \& Sons Ltd.}, + year = {2007}, + author = {Holzapfel, G. A.}, + address = {West Sussex, England}, + isbn = {0-471-82304-X} +} + +@Article{Linder2011a, + author = {Linder, C. and Tkachuk, M. and Miehe, C.}, + title = {A micromechanically motivated diffusion-based transient network model and its incorporation into finite rubber viscoelasticity}, + journal = {Journal of the Mechanics and Physics of Solids}, + year = {2011}, + volume = {59}, + number = {10}, + pages = {2134--2156}, + month = oct, + doi = {10.1016/j.jmps.2011.05.005}, + publisher = {Elsevier}, +} + +@Article{Pelteret2018a, + author = {Pelteret, J.-P. V. and Walter, B. and Steinmann, P.}, + title = {Application of metaheuristic algorithms to the identification of nonlinear magneto-viscoelastic constitutive parameters}, + journal = {Journal of Magnetism and Magnetic Materials}, + year = {2018}, + volume = {464}, + pages = {116--131}, + month = oct, + doi = {10.1016/j.jmmm.2018.02.094}, + publisher = {Elsevier {BV}}, +} + +@Article{Koprowski-Theiss2011a, + author = {Koprowski-Theiss, N. and Johlitz, M. and Diebels, S.}, + title = {Characterizing the time dependence of filled {EPDM}}, + journal = {Rubber Chemistry and Technology}, + year = {2011}, + volume = {84}, + number = {2}, + pages = {147--165}, + month = jun, + doi = {10.5254/1.3570527} +} + + % ------------------------------------ % Step 74 % ------------------------------------ @@ -1072,7 +1221,8 @@ eprint = {http://dx.doi.org/10.1137/0917003} title = {Code Complete}, publisher = {Microsoft Press}, year = 2004, - edition = {second}} + edition = {second} +} @article{gottlieb2001strong, title={Strong stability-preserving high-order time discretization methods}, diff --git a/doc/doxygen/tutorial/tutorial.h.in b/doc/doxygen/tutorial/tutorial.h.in index d83057834b..76715f8b34 100644 --- a/doc/doxygen/tutorial/tutorial.h.in +++ b/doc/doxygen/tutorial/tutorial.h.in @@ -598,6 +598,15 @@ * * * + * step-71 + * Constitutive modelling: a demonstration of how automatic and symbolic + * differentiation can be used to rapidly implement a complex coupled constitutive + * law. + *
Keywords: Automatic differentiation, Symbolic differentiation, + * Constitutive modelling + * + * + * * step-74 * The Symmetric interior penalty Galerkin (SIPG) method for Poisson's equation. *
Keywords: MeshWorker::mesh_loop(), FEInterfaceValues, ConvergenceTable @@ -720,7 +729,8 @@ * step-42, * step-43, * step-50, - * step-55 + * step-55, + * step-71 * * * @@ -798,7 +808,8 @@ * step-60, * step-62, * step-69, - * step-70 + * step-70, + * step-71 * * * @@ -877,9 +888,12 @@ * * * - * Computing Jacobians from residuals, automatic differentiation + * Computing Jacobians from residuals, automatic and symbolic differentiation + * + * + * step-33, + * step-71 * - * step-33 * * * diff --git a/examples/step-33/doc/builds-on b/examples/step-33/doc/builds-on index 20c3064395..396f9c7478 100644 --- a/examples/step-33/doc/builds-on +++ b/examples/step-33/doc/builds-on @@ -1 +1 @@ -step-12 +step-12 step-71 diff --git a/examples/step-71/CMakeLists.txt b/examples/step-71/CMakeLists.txt new file mode 100644 index 0000000000..7518476d40 --- /dev/null +++ b/examples/step-71/CMakeLists.txt @@ -0,0 +1,59 @@ +## +# CMake script +## + +# Set the name of the project and target: +SET(TARGET "step-71") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) + +FIND_PACKAGE(deal.II 9.3.0 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +# +# Are all dependencies fulfilled? +# +IF(NOT ((DEAL_II_WITH_TRILINOS AND DEAL_II_TRILINOS_WITH_SACADO) AND DEAL_II_WITH_SYMENGINE)) # keep in one line + MESSAGE(FATAL_ERROR " +Error! This tutorial requires a deal.II library that was configured with the following options: + DEAL_II_WITH_TRILINOS = ON + DEAL_II_TRILINOS_WITH_SACADO = ON +and + DEAL_II_WITH_SYMENGINE = ON +However, the deal.II library found at ${DEAL_II_PATH} was configured with these options + DEAL_II_WITH_TRILINOS = ${DEAL_II_WITH_TRILINOS} + DEAL_II_TRILINOS_WITH_SACADO = ${DEAL_II_TRILINOS_WITH_SACADO} + DEAL_II_WITH_SYMENGINE = ${DEAL_II_WITH_SYMENGINE} +which conflict with the requirements. +One or both of the aforementioned combinations of prerequisites are not met by your installation, but at least one is required for this tutorial step." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +SET(CLEAN_UP_FILES *.log *.gmv *.gnuplot *.gpl *.eps *.pov *.ucd *.d2 *.vtu *.pvtu) +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/examples/step-71/doc/builds-on b/examples/step-71/doc/builds-on new file mode 100644 index 0000000000..8b13789179 --- /dev/null +++ b/examples/step-71/doc/builds-on @@ -0,0 +1 @@ + diff --git a/examples/step-71/doc/intro.dox b/examples/step-71/doc/intro.dox new file mode 100644 index 0000000000..1ad1cfa8b9 --- /dev/null +++ b/examples/step-71/doc/intro.dox @@ -0,0 +1,678 @@ +
+ +This program was contributed by Jean-Paul Pelteret. + + + +

Introduction

+ +The aim of this tutorial is, quite simply, to introduce the fundamentals of both +[automatic](https://en.wikipedia.org/wiki/Automatic_differentiation) +and [symbolic differentiation](https://en.wikipedia.org/wiki/Computer_algebra) +(respectively abbreviated as AD +and SD): Ways in which one can, in source code, describe a function +$\mathbf f(\mathbf x)$ and automatically also obtain a representation of derivatives +$\nabla \mathbf f(\mathbf x)$ (the "Jacobian"), +$\nabla^2 \mathbf f(\mathbf x)$ (the "Hessian"), etc., without having +to write additional lines of code. Doing this is quite helpful in +solving nonlinear or optimization problems where one would like to +only describe the nonlinear equation or the objective function in the +code, without having to also provide their derivatives (which are +necessary for a Newton method for solving a nonlinear problem, or for +finding a minimizer). + +Since AD and SD tools are somewhat independent of finite elements and boundary value +problems, this tutorial is going to be different to the others that you may have +read beforehand. It will focus specifically on how these frameworks work and +the principles and thinking behind them, and will forgo looking at them in the +direct context of a finite element simulation. + +We will, in fact, look at two different sets of problems that have greatly +different levels of complexity, but when framed properly hold sufficient +similarity that the same AD and SD frameworks can be leveraged. With these +examples the aim is to build up an understanding of the steps that are required +to use the AD and SD tools, the differences between them, and hopefully identify +where they could be immediately be used in order to improve or simplify existing +code. + +It's plausible that you're wondering what AD and SD are, in the first place. Well, +that question is easy to answer but without context is not very insightful. So +we're not going to cover that in this introduction, but will rather defer this +until the first introductory example where we lay out the key points as this +example unfolds. To complement this, we should mention that the core theory for +both frameworks is extensively discussed in the @ref auto_symb_diff module, so +it bears little repeating here. + +Since we have to pick *some* sufficiently interesting topic to investigate +and identify where AD and SD can be used effectively, the main problem that's +implemented in the second half of the tutorial is one of modeling a coupled +constitutive law, specifically a magneto-active material (with hysteretic effects). +As a means of an introduction to that, later in the introduction some grounding +theory for that class of materials will be presented. +Naturally, this is not a field (or even a class of materials) that is of +interest to a wide audience. Therefore, the author wishes to express up front +that this theory and any subsequent derivations mustn't be considered the focus +of this tutorial. Instead, keep in mind the complexity of the problem that arises +from the relatively innocuous description of the constitutive law, and what we +might (in the context of a boundary value problem) need to derive from that. +We will perform some computations with these constitutive laws at the level of a +representative continuum point (so, remaining in the realm of continuum +mechanics), and will produce some benchmark results around which we can frame +a final discussion on the topic of computational performance. + +Once we have the foundation upon which we can build further concepts, we +will see how AD in particular can be exploited at a finite element (rather than +continuum) level: this is a topic that is covered in step-33. But before then, let's take a moment to +think about why we might want to consider using these sorts of tools, and what +benefits they can potentially offer you. + + +

A motivation: Why would I use these tools?

+ +The primary driver for using AD or SD is typically that there is some situation +that requires differentiation to be performed, and that doing so is sufficiently +challenging to make the prospect of using an external tool to perform that specific +task appealing. A broad categorization for the circumstances under which AD or +SD can be rendered most useful include (but are probably not limited to) the +following: +- Rapid prototyping: For a new class of problems where you're trying to + implement a solution quickly, and want to remove some of the intricate details + (in terms of both the mathematics as well as the organizational structure of + the code itself). You might be willing to justify any additional computational + cost, which would be offset by an increased agility in restructuring your code + or modifying the part of the problem that is introducing some complex nonlinearity + with minimal effort. +- Complex problems: It could very well be that some problems just happen to have + a nonlinearity that is incredibly challenging to linearize or formulate by hand. + Having this challenge taken care of for you by a tool that is, for the most part, + robust, reliable, and accurate may alleviate some of the pains in implementing + certain problems. Examples of this include step-15, where the + derivative of the nonlinear PDE we solve is not incredibly difficult + to derive, but sufficiently cumbersome that one has to pay attention + in doing so by hand, and where implementing the corresponding finite + element formulation of the Newton step takes more than just the few + lines that it generally takes to implement the bilinear form; + step-33 (where we actually use AD) is an even more extreme example. +- Verification: For materials and simulations that exhibit nonlinear response, + an accurate rather than only approximate material tangent (the term mechanical engineers use for + the derivative of a material law) can be the difference between convergent and + divergent behavior, especially at high external (or coupling) loads. + As the complexity of the problem increases, so do the opportunities to introduce + subtle (or, perhaps, not-so-subtle) errors that produce predictably negative + results. + Additionally, there is a lot to be gained by verifying that the implementation is + completely correct. For example, certain categories of problems are known to exhibit + instabilities, and therefore when you start to lose quadratic convergence in a + nonlinear solver (e.g., Newton's method) then this may not be a huge surprise to + the investigator. However, it is hard (if not impossible) to distinguish between + convergence behavior that is produced as you near an unstable solution and when + you simply have an error in the material or finite element linearization, and + start to drift off the optimal convergence path due to that. Having a + method of verifying the correctness of the implementation of a constitutive law + linearization, for example, is perhaps the only meaningful way that you can + use to catch such errors, assuming that you've got nobody else to scrutinize your code. + Thankfully, with some tactical programming it is quite straight-forward to structure + a code for reuse, such that you can use the same classes in production code and + directly verify them in, for instance, a unit-test framework. + +This tutorial program will have two parts: One where we just introduce +the basic ideas of automatic and symbolic differentiation support in +deal.II using a simple set of examples; and one where we apply this to +a realistic but much more complicated case. For that second half, the +next section will provide some background on magneto-mechanical +materials -- you can skip this section if all you want to learn +about is what AD and SD actually are, but you probably want to read +over this section if you are interested in how to apply AD and SD for +concrete situations. + + +

Theory for magneto-mechanical materials

+ +

Thermodynamic principles

+ +As a prelude to introducing the coupled magneto-mechanical material law that we'll use +to model a magneto-active polymer, we'll start with a very concise summary of +the salient thermodynamics to which these constitutive laws must subscribe. +The basis for the theory, as summarized here, is described in copious detail by +Truesdell and Toupin @cite Truesdell1960a and Coleman and Noll @cite Coleman1963a, +and follows the logic laid out by Holzapfel @cite Holzapfel2007a. + +Starting from the first law of thermodynamics, and following a few technical +assumptions, it can be shown the the balance between the kinetic plus internal +energy rates and the power supplied to the system from external +sources is given by the following relationship that equates the rate +of change of the energy in an (arbitrary) volume $V$ on the left, and +the sum of forces acting on that volume on the right: +@f[ + D_{t} \int\limits_{V} \left[ + \frac{1}{2} \rho_{0} \mathbf{v} \cdot \mathbf{v} + + U^{*}_{0} \right] dV += \int\limits_{V} \left[ + \rho_{0} \mathbf{v} \cdot \mathbf{a} + + \mathbf{P}^{\text{tot}} : \dot{\mathbf{F}} + + \boldsymbol{\mathbb{H}} \cdot \dot{\boldsymbol{\mathbb{B}}} + + \mathbb{E} \cdot \dot{\mathbb{D}} + - D_{t} M^{*}_{0} + - \nabla_{0} \cdot \mathbf{Q} + + R_{0} \right] dV . +@f] +Here $D_{t}$ represents the total time derivative, +$\rho_{0}$ is the material density as measured in the Lagrangian reference frame, +$\mathbf{v}$ is the material velocity and $\mathbf{a}$ its acceleration, +$U^{*}_{0}$ is the internal energy per unit reference volume, +$\mathbf{P}^{\text{tot}}$ is the total Piola stress tensor and $\dot{\mathbf{F}}$ is +the time rate of the deformation gradient tensor, +$\boldsymbol{\mathbb{H}}$ and $\boldsymbol{\mathbb{B}}$ are, respectively, the magnetic field vector and the +magnetic induction (or magnetic flux density) vector, +$\mathbb{E}$ and $\mathbb{D}$ are the electric field vector and electric +displacement vector, and +$\mathbf{Q}$ and $R_{0}$ represent the referential thermal flux vector and thermal +source. +The material differential operator +$\nabla_{0} (\bullet) \dealcoloneq \frac{d(\bullet)}{d\mathbf{X}}$ +where $\mathbf{X}$ is the material position vector. +With some rearrangement of terms, invoking the arbitrariness of the integration +volume $V$, the total internal energy density rate $\dot{E}_{0}$ can be identified as +@f[ + \dot{E}_{0} += \mathbf{P}^{\text{tot}} : \dot{\mathbf{F}} + + \boldsymbol{\mathbb{H}} \cdot \dot{\boldsymbol{\mathbb{B}}} + + \mathbb{E} \cdot \dot{\mathbb{D}} + - \nabla_{0} \cdot \mathbf{Q} + + R_{0} . +@f] +The total internal energy includes contributions that arise not only due to +mechanical deformation (the first term), and thermal fluxes and sources (the +fourth and fifth terms), but also due to the intrinsic energy stored in the +magnetic and electric fields themselves (the second and third terms, +respectively). + +The second law of thermodynamics, known also as the entropy inequality principle, +informs us that certain thermodynamic processes are irreversible. After accounting +for the total entropy and rate of entropy input, the Clausius-Duhem inequality +can be derived. In local form (and in the material configuration), this reads +@f[ + \theta \dot{\eta}_{0} + - R_{0} + + \nabla_{0} \cdot \mathbf{Q} + - \frac{1}{\theta} \nabla_{0} \theta \cdot \mathbf{Q} + \geq 0 . +@f] +The quantity $\theta$ is the absolute temperature, and +$\eta_{0}$ represents the entropy per unit reference volume. + +Using this to replace $R_{0} - \nabla_{0} \cdot \mathbf{Q}$ in the result +stemming from the first law of thermodynamics, we now have the relation +@f[ + \mathbf{P}^{\text{tot}} : \dot{\mathbf{F}} + + \boldsymbol{\mathbb{H}} \cdot \dot{\boldsymbol{\mathbb{B}}} + + \mathbb{E} \cdot \dot{\mathbb{D}} + + \theta \dot{\eta}_{0} + - \dot{E}_{0} + - \frac{1}{\theta} \nabla_{0} \theta \cdot \mathbf{Q} + \geq 0 . +@f] +On the basis of Fourier's law, which informs us that heat flows from regions +of high temperature to low temperature, the last term is always positive and +can be ignored. +This renders the local dissipation inequality +@f[ + \mathbf{P}^{\text{tot}} : \dot{\mathbf{F}} + + \boldsymbol{\mathbb{H}} \cdot \dot{\boldsymbol{\mathbb{B}}} + + \mathbb{E} \cdot \dot{\mathbb{D}} + - \left[ \dot{E}_{0} - \theta \dot{\eta}_{0} \right] + \geq 0 . +@f] +It is postulated @cite Holzapfel2007a that the Legendre transformation +@f[ + \psi^{*}_{0} += \psi^{*}_{0} \left( \mathbf{F}, \boldsymbol{\mathbb{B}}, \mathbb{D}, \theta \right) += E_{0} - \theta \eta_{0} , +@f] +from which we may define the free energy density function $\psi^{*}_{0}$ with the stated +parameterization, exists and is valid. +Taking the material rate of this equation and substituting it into the local +dissipation inequality results in the generic expression +@f[ + \mathcal{D}_{\text{int}} + = \mathbf{P}^{\text{tot}} : \dot{\mathbf{F}} + + \boldsymbol{\mathbb{H}} \cdot \dot{\boldsymbol{\mathbb{B}}} + + \mathbb{E} \cdot \dot{\mathbb{D}} + - \dot{\theta} \eta_{0} + - \dot{\psi}^{*}_{0} \left( \mathbf{F}, \boldsymbol{\mathbb{B}}, \mathbb{D}, \theta \right) + \geq 0 . +@f] +Under the assumption of isothermal conditions, and that the electric field does +not excite the material in a manner that is considered non-negligible, then this +dissipation inequality reduces to +@f[ + \mathcal{D}_{\text{int}} + = \mathbf{P}^{\text{tot}} : \dot{\mathbf{F}} + + \boldsymbol{\mathbb{H}} \cdot \dot{\boldsymbol{\mathbb{B}}} + - \dot{\psi}^{*}_{0} \left( \mathbf{F}, \boldsymbol{\mathbb{B}} \right) + \geq 0 . +@f] + +

Constitutive laws

+ +When considering materials that exhibit mechanically dissipative behavior, +it can be shown that this can be captured within the dissipation inequality +through the augmentation of the material free energy density function with additional +parameters that represent internal variables @cite Holzapfel1996a. Consequently, +we write it as +@f[ + \mathcal{D}_{\text{int}} + = \mathbf{P}^{\text{tot}} : \dot{\mathbf{F}} + + \boldsymbol{\mathbb{H}} \cdot \dot{\boldsymbol{\mathbb{B}}} + - \dot{\psi}^{*}_{0} \left( \mathbf{F}, \mathbf{F}_{v}^{i}, \boldsymbol{\mathbb{B}} \right) + \geq 0 . +@f] +where $\mathbf{F}_{v}^{i} = \mathbf{F}_{v}^{i} \left( t \right)$ represents the +internal variable (which acts like a measure of the deformation gradient) +associated with the `i`th mechanical dissipative (viscous) mechanism. +As can be inferred from its parameterization, each of these internal parameters +is considered to evolve in time. +Currently the free energy density function $\psi^{*}_{0}$ is parameterized in terms of +the magnetic induction $\boldsymbol{\mathbb{B}}$. This is the natural parameterization that +comes as a consequence of the considered balance laws. Should such a class of +materials to be incorporated within a finite-element model, it would be ascertained +that a certain formulation of the magnetic problem, known as the magnetic vector +potential formulation, would need to be adopted. This has its own set of challenges, +so where possible the more simple magnetic scalar potential formulation may be +preferred. In that case, the magnetic problem needs to be parameterized in terms +of the magnetic field $\boldsymbol{\mathbb{H}}$. To make this re-parameterization, we execute +a final Legendre transformation +@f[ + \tilde{\psi}_{0} \left( \mathbf{F}, \mathbf{F}_{v}^{i}, \boldsymbol{\mathbb{H}} \right) + = \psi^{*}_{0} \left( \mathbf{F}, \mathbf{F}_{v}^{i}, \boldsymbol{\mathbb{B}} \right) + - \boldsymbol{\mathbb{H}} \cdot \boldsymbol{\mathbb{B}} . +@f] +At the same time, we may take advantage of the principle of material frame +indifference in order to express the energy density function in terms of symmetric +deformation measures: +@f[ + \psi_{0} \left( \mathbf{C}, \mathbf{C}_{v}^{i}, \boldsymbol{\mathbb{H}} \right) + = \tilde{\psi}_{0} \left( \mathbf{F}, \mathbf{F}_{v}^{i}, \boldsymbol{\mathbb{H}} \right) . +@f] +The upshot of these two transformations (leaving out considerable explicit and +hidden details) renders the final expression for the reduced dissipation +inequality as +@f[ + \mathcal{D}_{\text{int}} + = \mathbf{S}^{\text{tot}} : \frac{1}{2} \dot{\mathbf{C}} + - \boldsymbol{\mathbb{B}} \cdot \dot{\boldsymbol{\mathbb{H}}} + - \dot{\psi}_{0} \left( \mathbf{C}, \mathbf{C}_{v}^{i}, \boldsymbol{\mathbb{H}} \right) + \geq 0 . +@f] +(Notice the sign change on the second term on the right hand side, and the +transfer of the time derivative to the magnetic induction vector.) +The stress quantity $\mathbf{S}^{\text{tot}}$ is known as the total Piola-Kirchhoff +stress tensor and its energy conjugate $\mathbf{C} = \mathbf{F}^{T} \cdot \mathbf{F}$ +is the right Cauchy-Green deformation tensor, and +$\mathbf{C}_{v}^{i} = \mathbf{C}_{v}^{i} \left( t \right)$ is the re-parameterized +internal variable associated with the `i`th mechanical dissipative (viscous) +mechanism. + +Expansion of the material rate of the energy density function, and rearrangement of the +various terms, results in the expression +@f[ + \mathcal{D}_{\text{int}} + = \left[ \mathbf{S}^{\text{tot}} - 2 \frac{\partial \psi_{0}}{\partial \mathbf{C}} \right] : \frac{1}{2} \dot{\mathbf{C}} + - \sum\limits_{i}\left[ 2 \frac{\partial \psi_{0}}{\partial \mathbf{C}_{v}^{i}} \right] : \frac{1}{2} \dot{\mathbf{C}}_{v}^{i} + + \left[ - \boldsymbol{\mathbb{B}} - \frac{\partial \psi_{0}}{\partial \boldsymbol{\mathbb{H}}} \right] \cdot \dot{\boldsymbol{\mathbb{H}}} + \geq 0 . +@f] +At this point, its worth noting the use of the +[partial derivatives](https://en.wikipedia.org/wiki/Partial_derivative) +$\partial \left( \bullet \right)$. This is an important detail that will be +fundamental to a certain design choice made within the tutorial. +As brief reminder of what this signifies, the partial derivative of a +multi-variate function returns the derivative of that function with respect +to one of those variables while holding the others constant: +@f[ + \frac{\partial f\left(x, y\right)}{\partial x} + = \frac{d f\left(x, y\right)}{d x} \vert_{y} . +@f] +More specific to what's encoded in the dissipation inequality (with the very general +free energy density function $\psi_{0}$ with its parameterization yet to be formalized), +if one of the input variables is a function of another, it is also held constant +and the chain rule does not propagate any further, while the computing total +derivative would imply judicious use of the chain rule. This can be better +understood by comparing the following two statements: +@f{align*} + \frac{\partial f\left(x, y\left(x\right)\right)}{\partial x} + &= \frac{d f\left(x, y\left(x\right)\right)}{d x} \vert_{y} \\ + \frac{d f\left(x, y\left(x\right)\right)}{d x} + &= \frac{d f\left(x, y\left(x\right)\right)}{d x} \vert_{y} + + \frac{d f\left(x, y\left(x\right)\right)}{d y} \vert_{x} \frac{d y\left(x\right)}{x} . +@f} + +Returning to the thermodynamics of the problem, we next exploit the arbitrariness +of the quantities $\dot{\mathbf{C}}$ and $\dot{\boldsymbol{\mathbb{H}}}$, +by application of the Coleman-Noll procedure @cite Coleman1963a, @cite Coleman1967a. +This leads to the identification of the kinetic conjugate quantities +@f[ + \mathbf{S}^{\text{tot}} + = \mathbf{S}^{\text{tot}} \left( \mathbf{C}, \mathbf{C}_{v}^{i}, \boldsymbol{\mathbb{H}} \right) + \dealcoloneq 2 \frac{\partial \psi_{0} \left( \mathbf{C}, \mathbf{C}_{v}^{i}, \boldsymbol{\mathbb{H}} \right)}{\partial \mathbf{C}} , \\ + \boldsymbol{\mathbb{B}} + = \boldsymbol{\mathbb{B}} \left( \mathbf{C}, \mathbf{C}_{v}^{i}, \boldsymbol{\mathbb{H}} \right) + \dealcoloneq - \frac{\partial \psi_{0} \left( \mathbf{C}, \mathbf{C}_{v}^{i}, \boldsymbol{\mathbb{H}} \right)}{\partial \boldsymbol{\mathbb{H}}} . +@f] +(Again, note the use of the partial derivatives to define the stress and magnetic +induction in this generalized setting.) +From what terms remain in the dissipative power (namely those related to the +mechanical dissipative mechanisms), if they are assumed to be independent of +one another then, for each mechanism `i`, +@f[ + \frac{\partial \psi_{0}}{\partial \mathbf{C}_{v}^{i}} : \dot{\mathbf{C}}_{v}^{i} + \leq 0 . +@f] +This constraint must be satisfies through the appropriate choice of free energy +function, as well as a carefully considered evolution law for the internal +variables. + +In the case that there are no dissipative mechanisms to be captured within the +constitutive model (e.g., if the material to be modelled is magneto-hyperelastic) +then the free energy density function +$\psi_{0} = \psi_{0} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right)$ reduces to a stored +energy density function, and the total stress and magnetic induction can be simplified +@f{align*}{ + \mathbf{S}^{\text{tot}} + = \mathbf{S}^{\text{tot}} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) + &\dealcoloneq 2 \frac{d \psi_{0} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right)}{d \mathbf{C}} , \\ + \boldsymbol{\mathbb{B}} + = \boldsymbol{\mathbb{B}} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) + &\dealcoloneq - \frac{d \psi_{0} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right)}{d \boldsymbol{\mathbb{H}}} , +@f} +where the operator $d$ denotes the total derivative operation. + +For completeness, the linearization of the stress tensor and magnetic induction +are captured within the fourth-order total referential elastic tangent tensor +$\mathcal{H}^{\text{tot}} $, the second-order magnetostatic tangent tensor $\mathbb{D}$ and the +third-order total referential magnetoelastic coupling tensor $\mathfrak{P}^{\text{tot}}$. +Irrespective of the parameterization of $\mathbf{S}^{\text{tot}}$ and $\boldsymbol{\mathbb{B}}$, +these quantities may be computed by +@f{align*}{ + \mathcal{H}^{\text{tot}} + &= 2 \frac{d \mathbf{S}^{\text{tot}}}{d \mathbf{C}} , \\ + \mathbb{D} + &= \frac{d \boldsymbol{\mathbb{B}}}{d \boldsymbol{\mathbb{H}}} , \\ + \mathfrak{P}^{\text{tot}} + &= - \frac{d \mathbf{S}^{\text{tot}}}{d \boldsymbol{\mathbb{H}}} , \\ + \left[ \mathfrak{P}^{\text{tot}} \right]^{T} + &= 2 \frac{d \boldsymbol{\mathbb{B}}}{d \mathbf{C}} . +@f} +For the case of rate-dependent materials, this expands to +@f{align*}{ + \mathcal{H}^{\text{tot}} \left( \mathbf{C}, \mathbf{C}_{v}^{i}, \boldsymbol{\mathbb{H}} \right) + &= 2 \frac{d^{2} \psi_{0} \left( \mathbf{C}, \mathbf{C}_{v}^{i}, \boldsymbol{\mathbb{H}} \right)}{d \mathbf{C} \otimes \partial \mathbf{C}} , \\ + \mathbb{D} \left( \mathbf{C}, \mathbf{C}_{v}^{i}, \boldsymbol{\mathbb{H}} \right) + &= -\frac{d^{2} \psi_{0} \left( \mathbf{C}, \mathbf{C}_{v}^{i}, \boldsymbol{\mathbb{H}} \right)}{d \boldsymbol{\mathbb{H}} \otimes \partial \boldsymbol{\mathbb{H}}} , \\ + \mathfrak{P}^{\text{tot}} \left( \mathbf{C}, \mathbf{C}_{v}^{i}, \boldsymbol{\mathbb{H}} \right) + &= - 2 \frac{d^{2} \psi_{0} \left( \mathbf{C}, \mathbf{C}_{v}^{i}, \boldsymbol{\mathbb{H}} \right)}{d \boldsymbol{\mathbb{H}} \otimes \partial \mathbf{C}} , \\ + \left[ \mathfrak{P}^{\text{tot}} \left( \mathbf{C}, \mathbf{C}_{v}^{i}, \boldsymbol{\mathbb{H}} \right) \right]^{T} + &= - 2 \frac{d^{2} \psi_{0} \left( \mathbf{C}, \mathbf{C}_{v}^{i}, \boldsymbol{\mathbb{H}} \right)}{d \mathbf{C} \otimes \partial \boldsymbol{\mathbb{H}}} , +@f} +while for rate-independent materials the linearizations are +@f{align*}{ + \mathcal{H}^{\text{tot}} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) + &= 2 \frac{d^{2} \psi_{0} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right)}{d \mathbf{C} \otimes d \mathbf{C}} , \\ + \mathbb{D} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) + &= -\frac{d^{2} \psi_{0} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right)}{d \boldsymbol{\mathbb{H}} \otimes d \boldsymbol{\mathbb{H}}} , \\ + \mathfrak{P}^{\text{tot}} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) + &= - 2 \frac{d^{2} \psi_{0} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right)}{d \boldsymbol{\mathbb{H}} \otimes d \mathbf{C}} , \\ + \left[ \mathfrak{P}^{\text{tot}} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) \right]^{T} + &= - 2 \frac{d^{2} \psi_{0} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right)}{d \mathbf{C} \otimes d \boldsymbol{\mathbb{H}}} . +@f} +The subtle difference between them is the application of a partial derivative during +the calculation of the first derivatives. We'll see later how this affects the choice +of AD versus SD for this specific application. For now, we'll simply introduce +the two specific materials that are implemented within this tutorial. + +
Magnetoelastic constitutive law
+ +The first material that we'll consider is one that is governed by a +magneto-hyperelastic constitutive law. This material responds to both +deformation as well as immersion in a magnetic field, but exhibits no +time- or history-dependent behavior (such as dissipation through viscous +damping or magnetic hysteresis, etc.). The *stored energy density +function* for such a material is only parameterized in terms of the +(current) field variables, but not their time derivatives or past values. + +We'll choose the energy density function, which captures both the energy +stored in the material due to deformation and magnetization, as well as +the energy stored in the magnetic field itself, to be +@f[ + \psi_{0} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) += \frac{1}{2} \mu_{e} f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} \right) + \left[ \text{tr}(\mathbf{C}) - d - 2 \ln (\text{det}(\mathbf{F})) + \right] ++ \lambda_{e} \ln^{2} \left(\text{det}(\mathbf{F}) \right) +- \frac{1}{2} \mu_{0} \mu_{r} \text{det}(\mathbf{F}) + \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + \boldsymbol{\mathbb{H}} \right] +@f] +with +@f[ + f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} \right) += 1 + \left[ \frac{\mu_{e}^{\infty}}{\mu_{e}} - 1 \right] + \tanh \left( 2 \frac{\boldsymbol{\mathbb{H}} \cdot + \boldsymbol{\mathbb{H}}} + {\left(h_{e}^{\text{sat}}\right)^{2}} \right) +@f] +and for which the variable $d = \text{tr}(\mathbf{I})$ ($\mathbf{I}$ +being the rank-2 identity tensor) represents the spatial dimension and +$\mathbf{F}$ is the deformation gradient tensor. To give some brief +background to the various components of $\psi_{0}$, the first two terms +bear a great resemblance to the stored energy density function for a +(hyperelastic) Neohookean material. The only difference between what's +used here and the Neohookean material is the scaling of the elastic shear +modulus by the magnetic field-sensitive saturation function $f_{\mu_{e}} +\left( \boldsymbol{\mathbb{H}} \right)$ (see @cite Pelteret2018a, equation +29). This function will, in effect, cause the material to stiffen in the +presence of a strong magnetic field. As it is governed by a sigmoid-type +function, the shear modulus will asymptotically converge on the specified +saturation shear modulus. It can also be shown that the last term in +$\psi_{0}$ is the stored energy density function for magnetic field (as +derived from first principles), scaled by the relative permeability +constant. This definition collectively implies that the material is +linearly magnetized, i.e., the magnetization vector and magnetic field +vector are aligned. (This is certainly not obvious with the magnetic energy +stated in its current form, but when the magnetic induction and magnetization +are derived from $\psi_{0}$ and all magnetic fields are expressed in the +current configuration then this correlation becomes clear.) +As for the specifics of what the magnetic induction, stress tensor, and the +various material tangents look like, we'll defer presenting these to the +tutorial body where the full, unassisted implementation of the constitutive +law is defined. + +
Magneto-viscoelastic constitutive law
+ +The second material that we'll formulate is one that for a +magneto-viscoelastic material with a single dissipative mechanism `i`. +The *free energy density function* that we'll be considering is defined as +@f{align*}{ + \psi_{0} \left( \mathbf{C}, \mathbf{C}_{v}, \boldsymbol{\mathbb{H}} + \right) +&= \psi_{0}^{ME} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) ++ \psi_{0}^{MVE} \left( \mathbf{C}, \mathbf{C}_{v}, +\boldsymbol{\mathbb{H}} \right) +\\ \psi_{0}^{ME} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) +&= \frac{1}{2} \mu_{e} f_{\mu_{e}^{ME}} \left( \boldsymbol{\mathbb{H}} +\right) + \left[ \text{tr}(\mathbf{C}) - d - 2 \ln (\text{det}(\mathbf{F})) + \right] ++ \lambda_{e} \ln^{2} \left(\text{det}(\mathbf{F}) \right) +- \frac{1}{2} \mu_{0} \mu_{r} \text{det}(\mathbf{F}) + \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + \boldsymbol{\mathbb{H}} \right] +\\ \psi_{0}^{MVE} \left( \mathbf{C}, \mathbf{C}_{v}, +\boldsymbol{\mathbb{H}} \right) +&= \frac{1}{2} \mu_{v} f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} +\right) + \left[ \mathbf{C}_{v} : \left[ + \left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + \mathbf{C} \right] - d - \ln\left( + \text{det}\left(\mathbf{C}_{v}\right) \right) \right] +@f} +with +@f[ + f_{\mu_{e}}^{ME} \left( \boldsymbol{\mathbb{H}} \right) += 1 + \left[ \frac{\mu_{e}^{\infty}}{\mu_{e}} - 1 \right] + \tanh \left( 2 \frac{\boldsymbol{\mathbb{H}} \cdot + \boldsymbol{\mathbb{H}}} + {\left(h_{e}^{\text{sat}}\right)^{2}} \right) +@f] +@f[ + f_{\mu_{v}}^{MVE} \left( \boldsymbol{\mathbb{H}} \right) += 1 + \left[ \frac{\mu_{v}^{\infty}}{\mu_{v}} - 1 \right] + \tanh \left( 2 \frac{\boldsymbol{\mathbb{H}} \cdot + \boldsymbol{\mathbb{H}}} + {\left(h_{v}^{\text{sat}}\right)^{2}} \right) +@f] +and the evolution law +@f[ + \dot{\mathbf{C}}_{v} \left( \mathbf{C} \right) += \frac{1}{\tau} \left[ + \left[\left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + \mathbf{C}\right]^{-1} + - \mathbf{C}_{v} \right] +@f] +for the internal viscous variable. +We've chosen the magnetoelastic part of the energy +$\psi_{0}^{ME} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right)$ +to match that of the first material model that we explored, so this part +needs no further explanation. As for the viscous part $\psi_{0}^{MVE}$, +this component of the free energy (in conjunction with the evolution law for +the viscous deformation tensor) is taken from @cite Linder2011a (with the +additional scaling by the viscous saturation function described in +@cite Pelteret2018a). It is derived in a thermodynamically consistent +framework that, at its core, models the movement of polymer chains on a +micro-scale level. + +To proceed beyond this point, we'll also need to consider the time +discretization of the evolution law. +Choosing the implicit first-order backwards difference scheme, then +@f[ + \dot{\mathbf{C}}_{v} +\approx \frac{\mathbf{C}_{v}^{(t)} - \mathbf{C}_{v}^{(t-1)}}{\Delta t} += \frac{1}{\tau} \left[ + \left[\left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + \mathbf{C}\right]^{-1} + - \mathbf{C}_{v}^{(t)} \right] +@f] +where the superscript $(t)$ denotes that the quantity is taken at the +current timestep, and $(t-1)$ denotes quantities taken at the previous +timestep (i.e., a history variable). The timestep size $\Delta t$ is the +difference between the current time and that of the previous timestep. +Rearranging the terms so that all internal variable quantities at the +current time are on the left hand side of the equation, we get +@f[ +\mathbf{C}_{v}^{(t)} += \frac{1}{1 + \frac{\Delta t}{\tau_{v}}} \left[ + \mathbf{C}_{v}^{(t-1)} + + \frac{\Delta t}{\tau_{v}} + \left[\left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + \mathbf{C} \right]^{-1} + \right] +@f] +that matches @cite Linder2011a equation 54. + +

Rheological experiment

+ +Now that we have shown all of these formulas for the thermodynamics and theory +governing magneto-mechanics and constitutive models, let us outline what the +program will do with all of this. +We wish to do something *meaningful* with the materials laws that we've formulated, +and so it makes sense to subject them to some mechanical and magnetic loading +conditions that are, in some way, representative of some conditions that might +be found either in an application or in a laboratory setting. One way to achieve +that aim would be to embed these constitutive laws in a finite element model to +simulate a device. In this instance, though, we'll keep things simple (we are +focussing on the automatic and symbolic differentiation concepts, after all) +and will find a concise way to faithfully replicate an industry-standard +rheological experiment using an analytical expression for the loading conditions. + +The rheological experiment that we'll reproduce, +which idealizes a laboratory experiment that was used to characterize +magneto-active polymers, is detailed in @cite Pelteret2018a +(as well as @cite Pelteret2019a, in which it is documented along with the +real-world experiments). The images below provide a visual description of +the problem set up. + + + + + + +
+ +

+ The basic functional geometry of the parallel-plate rotational + rheometer. The smooth rotor (blue) applies a torque to an + experimental sample (red) of radius $r$ and height $H$ while an + axially aligned magnetic field generated by a a + magneto-rheological device. Although the time-dependent + deformation profile of the may be varied, one common experiment + would be to subject the material to a harmonic torsional + deformation of constant amplitude and frequency $\omega$. +

+ 		+ +

+ Schematic of the kinematics of the problem, assuming no + preloading or compression of the sample. A point $\mathbf{P}$ + located at azimuth $\Theta$ is displaced to location $\mathbf{p}$ + at azimuth $\theta = \Theta + \alpha$. +

+
+ +Under the assumptions that an incompressible medium is being tested, +and that the deformation profile through the sample thickness is linear, +then the displacement at some measurement point $\mathbf{X}$ within +the sample, expressed in radial coordinates, is +@f{align*} + r(\mathbf{X}) + &= \frac{R(X_{1}, X_{2})}{\sqrt{\lambda_{3}}} , \\ + \theta(\mathbf{X}) + & = \Theta(X_{1}, X_{2}) + \underbrace{\tau(t) + \lambda_{3} X_{3}}_{\alpha(X_{3}, t)} , \\ + z(\mathbf{X}) + &= \lambda_{3} X_{3} +@f} +where +$R(X_{1}, X_{2})$ and $\Theta(X_{1}, X_{2})$ are the radius at +-- and angle of -- the sampling point, +$\lambda_{3}$ is the (constant) axial deformation, +$\tau(t) = \frac{A}{RH} \sin\left(\omega t\right)$ is the time-dependent +torsion angle per unit length that will be prescribed using a +sinusoidally repeating oscillation of fixed amplitude $A$. +The magnetic field is aligned axially, i.e., in the $X_{3}$ direction. + +This summarizes everything that we need to fully characterize the idealized +loading at any point within the rheological sample. We'll set up the problem +in such a way that we "pick" a representative point with this sample, and +subject it to a harmonic shear deformation at a constant axial deformation +(by default, a compressive load) and a constant, axially applied magnetic +field. We will record the stress and magnetic induction at this point, and +will output that data to file for post-processing. Although its not necessary +for this particular problem, we will also be computing the tangents as well. +Even though they are not directly used in this particular piece of work, these +second derivatives are needed to embed the constitutive law within a +finite element model (one possible extension to this work). We'll therefore +take the opportunity to checked our hand calculations for correctness using +the assisted differentiation frameworks. + +

Suggested literature

+ +In addition to the already mentioned @ref auto_symb_diff module, the following are a few +references that discuss in more detail +- magneto-mechanics, and some aspects of automated differentiation frameworks: @cite Pao1978a, @cite Pelteret2019a, and +- the automation of finite element frameworks using AD and/or SD: @cite Logg2012a, @cite Korelc2016a. + +
diff --git a/examples/step-71/doc/kind b/examples/step-71/doc/kind new file mode 100644 index 0000000000..c1d9154931 --- /dev/null +++ b/examples/step-71/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/examples/step-71/doc/results.dox b/examples/step-71/doc/results.dox new file mode 100644 index 0000000000..c3ab8b3c01 --- /dev/null +++ b/examples/step-71/doc/results.dox @@ -0,0 +1,333 @@ +

Results

+ +

Introductory example

+ +The first exploratory example produces the following output. It is verified that +all three implementations produce identical results. +@code +> ./step-71 +Simple example using automatic differentiation... +... all calculations are correct! +Simple example using symbolic differentiation. +... all calculations are correct! +@endcode + +

Constitutive modelling

+ +To help summarize the results from the virtual experiment itself, below are some +graphs showing the shear stress, plotted against the shear strain, at a select +location within the material sample. The plots show the stress-strain curves under +three different magnetic loads, and for the last cycle of the (mechanical) +loading profile, when the rate-dependent material reaches a repeatable +("steady-state") response. These types of graphs are often referred to as +[Lissajous plots](https://en.wikipedia.org/wiki/Lissajous_curve). The area +of the ellipse that the curve takes for viscoelastic materials provides some +measure of how much energy is dissipated by the material, and its ellipticity +indicates the phase shift of the viscous response with respect to the elastic +response. + + + + + + +
+ +

+ Lissajous plot for the magneto-elastic material. +

+ 		+ +

+ Lissajous plot for the magneto-viscoelastic material. +

+
+ +It is not surprising to see that the magneto-elastic material response has an unloading +curve that matches the loading curve -- the material is non-dissipative after all. +But here it's clearly noticeable how the gradient of the curve increases as the +applied magnetic field increases. The tangent at any point along this curve is +related to the instantaneous shear modulus and, due to the way that the energy +density function was defined, we expect that the shear modulus increases as the +magnetic field strength increases. +We observe much the same behavior for the magneto-viscoelastic material. The major +axis of the ellipse traced by the loading-unloading curve has a slope that increases +as a greater magnetic load is applied. At the same time, the more energy is +dissipated by the material. + +As for the code output, this is what is printed to the console for the part +pertaining to the rheological experiment conducted with the magnetoelastic +material: +@code +Coupled magnetoelastic constitutive law using automatic differentiation. +Timestep = 0 @ time = 0s. +Timestep = 125 @ time = 0.314159s. +Timestep = 250 @ time = 0.628318s. +Timestep = 375 @ time = 0.942477s. +... +Timestep = 12250 @ time = 30.7876s. +Timestep = 12375 @ time = 31.1018s. +Timestep = 12500 @ time = 31.4159s. +... all calculations are correct! +@endcode + +And this portion of the output pertains to the experiment performed with the +magneto-viscoelastic material: +@code +Coupled magneto-viscoelastic constitutive law using symbolic differentiation. +Using LLVM optimizer. +Timestep = 0 @ time = 0s. +Timestep = 125 @ time = 0.314159s. +Timestep = 250 @ time = 0.628318s. +Timestep = 375 @ time = 0.942477s. +... +Timestep = 12250 @ time = 30.7876s. +Timestep = 12375 @ time = 31.1018s. +Timestep = 12500 @ time = 31.4159s. +... all calculations are correct! +@endcode + +The timer output is also emitted to the console, so we can compare time taken +to perform the hand- and assisted- calculations and get some idea of the overhead +of using the AD and SD frameworks. +Here are the timings taken from the magnetoelastic experiment using +the AD framework, based on the Sacado component of the Trilinos library: +@code ++---------------------------------------------+------------+------------+ +| Total wallclock time elapsed since start | 3.2s | | +| | | | +| Section | no. calls | wall time | % of total | ++---------------------------------+-----------+------------+------------+ +| Assisted computation | 12501 | 3.02s | 95% | +| Hand calculated | 12501 | 0.0464s | 1.5% | ++---------------------------------+-----------+------------+------------+ +@endcode +With respect to the computations performed using automatic differentiation +(as a reminder, this is with two levels of differentiation using the Sacado +library in conjunction with dynamic forward auto-differentiable types), we +observe that the assisted computations takes about $65 \times$ longer to +compute the desired quantities. This does seem like quite a lot of overhead +but, as mentioned in the introduction, it's entirely subjective and +circumstance-dependent as to whether or not this is acceptable or not: +Do you value computer time more than human time for doing the +necessary hand-computations of derivatives, verify their correctness, +implement them, and verify the correctness of the implementation? If +you develop a research code that will only be run for a relatively +small number of experiments, you might value your own time more. If +you develop a production code that will be run over and over on +10,000-core clusters for hours, your considerations might be different. +In any case, the one nice feature +of the AD approach is the "drop in" capability when functions and classes are +templated on the scalar type. This means that minimal effort is required to +start working with it. + +In contrast, the timings for magneto-viscoelastic material as implemented using +just-in-time (JIT) compiled symbolic algebra indicate that, at some non-negligible cost during +initialization, the calculations themselves are a lot more efficiently executed: +@code ++---------------------------------------------+------------+------------+ +| Total wallclock time elapsed since start | 1.34s | | +| | | | +| Section | no. calls | wall time | % of total | ++---------------------------------+-----------+------------+------------+ +| Assisted computation | 12501 | 0.376s | 28% | +| Hand calculated | 12501 | 0.368s | 27% | +| Initialize symbolic CL | 1 | 0.466s | 35% | ++---------------------------------+-----------+------------+------------+ +@endcode +Since the initialization phase need, most likely, only be executed once per +thread, this initial expensive phase can be offset by the repeated use of a +single Differentiation::SD::BatchOptimizer instance. Even though the +magneto-viscoelastic constitutive law has more terms to calculate when compared +to its magnetoelastic counterpart, it still is a whole order of magnitude faster +to execute the computations of the kinetic variables and tangents. And when compared +to the hand computed variant that uses the caching scheme, the calculation time +is nearly equal. So although using the symbolic framework requires a paradigm +shift in terms of how one implements and manipulates the symbolic expressions, +it can offer good performance and flexibility that the AD frameworks lack. + +On the point of data caching, the added cost of value caching for the +magneto-viscoelastic material implementation is, in fact, about a $6\times$ +increase in the time spent in `update_internal_data()` when compared to the +implementation using intermediate values for the numerical experiments conducted +with this material. Here's a sample output of the timing comparison extracted for +the "hand calculated" variant when the caching data structure is removed: +@code ++---------------------------------------------+------------+------------+ +| Total wallclock time elapsed since start | 1.01s | | +| | | | +| Section | no. calls | wall time | % of total | ++---------------------------------+-----------+------------+------------+ +| Assisted computation | 12501 | 0.361s | 36% | +| Hand calculated | 12501 | 0.0562s | 5.6% | +| Initialize symbolic CL | 1 | 0.469s | 47% | ++---------------------------------+-----------+------------+------------+ +@endcode + +With some minor adjustment we can quite easily test the different optimization +schemes for the batch optimizer. So let's compare the computational expense +associated with the `LLVM` batch optimizer setting versus the alternatives. +Below are the timings reported for the `lambda` optimization method (retaining +the use of CSE): +@code ++---------------------------------------------+------------+------------+ +| Total wallclock time elapsed since start | 3.87s | | +| | | | +| Section | no. calls | wall time | % of total | ++---------------------------------+-----------+------------+------------+ +| Assisted computation | 12501 | 3.12s | 81% | +| Hand calculated | 12501 | 0.394s | 10% | +| Initialize symbolic CL | 1 | 0.209s | 5.4% | ++---------------------------------+-----------+------------+------------+ +@endcode +The primary observation here is that an order of magnitude greater time is spent +in the "Assisted computation" section when compared to the `LLVM` approach. + +Last of all we'll test how `dictionary` substitution, in conjunction with CSE, +performs. Dictionary substitution simply does all of the evaluation within the +native CAS framework itself, with no transformation of the underlying data +structures taking place. Only the use of CSE, which caches intermediate results, +will provide any "acceleration" in this instance. With that in mind, here are +the results from this selection: +@code ++---------------------------------------------+------------+------------+ +| Total wallclock time elapsed since start | 1.54e+03s | | +| | | | +| Section | no. calls | wall time | % of total | ++---------------------------------+-----------+------------+------------+ +| Assisted computation | 12501 | 1.54e+03s | 1e+02% | +| Hand calculated | 12501 | 0.563s | 0% | +| Initialize symbolic CL | 1 | 0.184s | 0% | ++---------------------------------+-----------+------------+------------+ +@endcode +Needless to say, compared to the other two methods, these results took quite +some time to produce... The `dictionary` substitution +method is perhaps only really viable for simple expressions or when the number +of calls is sufficiently small. + +

So, which framework should I use?

+ +Perhaps you've been convinced that these tools have some merit, and can be +of immediate help or use to you. The obvious question now is which one to +use. Focussing specifically at a continuum point level, where you would be +using these frameworks to compute derivatives of a constitutive law in +particular, we can say the following: +- Automatic differentiation probably provides the simplest entry point into + the world of assisted differentiation. +- Given a sufficiently generic implementation of a constitutive framework, + AD can often be used as a drop-in replacement for the intrinsic scalar types + and the helper classes can then be leveraged to compute first (and possibly + higher order) derivatives with minimal effort. +- As a qualification to the above point, being a "drop-in replacement" does not + mean that you must not be contentious of what the algorithms that these numbers + are being passed through are doing. It is possible to inadvertently perform + an operation that would, upon differentiating, return an incorrect result. + So this is definitely something that one should be aware of. + A concrete example: When computing the eigenvalues of a tensor, if the tensor + is diagonal then a short-cut to the result is simply to return the diagonal + entries directly (as extracted from the input tensor). This is completely + correct in terms of computing the eigenvalues themselves, but not going + through the algorithm that would otherwise compute the eigenvalues for a + non-diagonal tensor has had an unintended side-effect, namely that the + eigenvalues appear (to the AD framework) to be completely decoupled from + one another and their cross-sensitivities are not encoded in the returned + result. Upon differentiating, many entries of the derivative tensor will + be missing. To fix this issue, one has to ensure that the standard eigenvalue + solving algorithm is used so that the sensitivities of the returned eigenvalues + with respect to one another are encoded in the result. +- Computations involving AD number types may be expensive. The expense increases + (sometimes quite considerably) as the order of the differential operations + increases. This may be mitigated by computational complexity of surrounding + operations (such as a linear solve, for example), but is ultimately problem + specific. +- AD is restricted to the case where only total derivatives are required. If a + differential operation requires a partial derivative with respect to an + independent variable then it is not appropriate to use it. +- Each AD library has its own quirks (sad to say but, in the author's experience, + true), so it may take some trial and error to find the appropriate library and + choice of AD number to suit your purposes. The reason for these "quirks" + often boils down to the overall philosophy behind the library (data structures, + the use of template meta-programming, etc.) as well as the mathematical + implementation of the derivative computations (for example, manipulations of + results using logarithmic functions to change basis might restrict the domain + for the input values -- details all hidden from the user, of course). + Furthermore, one library might be able to compute the desired results quicker + than another, so some initial exploration might be beneficial in that regard. +- Symbolic differentiation (well, the use of a CAS in general) provides the most + flexible framework with which to perform assisted computations. +- The SD framework can do everything that the AD frameworks can, with the + additional benefit of having low-level control over when certain manipulations + and operations are performed. +- Acceleration of expression evaluation is possible, potentially leading to + near-native performance of the SD framework compared to some hand implementations + (this comparison being dependent on the overall program design, of course) at + the expense of the initial optimization call. +- Clever use of the Differentiation::SD::BatchOptimizer could minimize the + expense of the costly call that optimizes the dependent expressions. + The possibility to serialize the Differentiation::SD::BatchOptimizer + that often (but not always) this expensive call can be done once and then + reused in a later simulation. +- If two or more material laws differ by only their material parameters, for + instance, then a single batch optimizer can be shared between them as long + as those material parameters are considered to be symbolic. The implication + of this is that you can "differentiate once, evaluate in many contexts". +- The SD framework may partially be used as a "drop-in replacement" for scalar + types, but one (at the very least) has to add some more framework around it + to perform the value substitution step, converting symbolic types to their + numerical counterparts. +- It may not be possible to use SD numbers within some specialized algorithms. + For example, if an algorithm has an exit point or code branch based off of + some concrete, numerical value that the (symbolic) input argument should take, + then obviously this isn't going to work. One either has to reimplement the + algorithm specifically for SD number types (somewhat inconvenient, but + frequently possible as conditionals are supported by the + Differentiation::SD::Expression class), or one must use a creative means + around this specific issue (e.g., introduce a symbolic expression that + represents the result returned by this algorithm, perhaps declaring it + to be a + [symbolic function](https://dealii.org/developer/doxygen/deal.II/namespaceDifferentiation_1_1SD.html#a876041f6048705c7a8ad0855cdb1bd7a) + if that makes sense within the context in which it is to be used. This can + later be substituted by its numerical values, and if declared a symbolic + function then its deferred derivatives may also be incorporated into the + calculations as substituted results.). +- The biggest drawback to using SD is that using it requires a paradigm shift, + and that one has to frame most problems differently in order to take the + most advantage of it. (Careful consideration of how the data structures + are used and reused is also essential to get it to work effectively.) This may + mean that one needs to play around with it a bit and build up an understanding + of what the sequence of typical operations is and what specifically each step + does in terms of manipulating the underlying data. If one has the time and + inclination to do so, then the benefits of using this tool may be substantial. + +

Possibilities for extension

+ +There are a few logical ways in which this program could be extended: +- Perhaps the most obvious extension would be to implement and test other constitutive models. + This could still be within the realm of coupled magneto-mechanical problems, perhaps considering + alternatives to the "Neo-Hookean"-type elastic part of the energy functions, changing the + constitutive law for the dissipative energy (and its associated evolution law), or including + magnetic hysteretic effects or damage models for the composite polymer that these material + seek to model. +- Of course, the implemented models could be modified or completely replaced with models that are + focussed on other aspects of physics, such as electro-active polymers, biomechanical materials, + elastoplastic media, etc. +- Implement a different time-discretization scheme for the viscoelastic evolution law. +- Instead of deriving everything directly from an energy density function, use the + Differentiation::AD::VectorFunction to directly linearize the kinetic quantities. + This would mean that only a once-differentiable auto-differentiable number type + would be required, and would certainly improve the performance greatly. + Such an approach would also offer the opportunity for dissipative materials, + such as the magneto-viscoelastic one consider here, to be implemented in + conjunction with AD. This is because the linearization invokes the total + derivative of the dependent variables with respect to the field variables, which + is exactly what the AD frameworks can provide. +- Investigate using other auto-differentiable number types and frameworks (such as + ADOL-C). Since each AD library has its own implementation, the choice of which + to use could result in performance increases and, in the most unfortunate cases, + more stable computations. It can at least be said that for the AD libraries that + deal.II supports, the accuracy of results should be largely unaffected by this decision. +- Embed one of these constitutive laws within a finite element simulation. + +With less effort, one could think about re-writing nonlinear problem +solvers such as the one implemented in step-15 using AD or SD +approaches to compute the Newton matrix. diff --git a/examples/step-71/doc/tooltip b/examples/step-71/doc/tooltip new file mode 100644 index 0000000000..b66dc87e3f --- /dev/null +++ b/examples/step-71/doc/tooltip @@ -0,0 +1 @@ +Coupled constitutive modeling using automatic and symbolic differentiation diff --git a/examples/step-71/step-71.cc b/examples/step-71/step-71.cc new file mode 100644 index 0000000000..bc9e23043f --- /dev/null +++ b/examples/step-71/step-71.cc @@ -0,0 +1,4048 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2021 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE.md at + * the top level directory of deal.II. + * + * --------------------------------------------------------------------- + * + * Author: Jean-Paul Pelteret, 2021 + */ + + +// We start by including all the necessary deal.II header files and some C++ +// related ones. +// This first header will give us access to a data structure that will allow +// us to store arbitrary data within it. +#include + +// Next come some core classes, including one that provides an implementation +// for time-stepping. +#include +#include +#include +#include +#include +#include +#include + +// Then some headers that define some useful coordinate transformations and +// kinematic relationships that are often found in nonlinear elasticity. +#include +#include +#include + +// The following two headers provide all of the functionality that we need +// to perform automatic differentiation, and use the symbolic computer algebra +// system that deal.II can utilize. The headers of all automatic +// differentiation and symbolic differentiation wrapper classes, and any +// ancillary data structures that are required, are all collected inside these +// unifying headers. +#include +#include + +// Including this header allows us the capability to write output to a +// file stream. +#include + + +// As per usual, the entire tutorial program is defined within its own unique +// namespace. +namespace Step71 +{ + using namespace dealii; + + // @sect3{An introductory example: The fundamentals of automatic and symbolic differentiation} + + // Automatic and symbolic differentiation have some magical and mystical + // qualities. Although their use in a project can be beneficial for a + // multitude of reasons, the barrier to understanding how to use these + // frameworks or how they can be leveraged may exceed the patience of + // the developer that is trying to (reliably) integrate them into their work. + // + // Although it is the wish of the author to successfully illustrate how these + // tools can be integrated into workflows for finite element modelling, it + // might be best to first take a step back and start right from the basics. + // So to start off with, we'll first have a look at differentiating a "simple" + // mathematical function using both frameworks, so that the fundamental + // operations (both their sequence and function) can be firmly established and + // understood with minimal complication. In the second part of this tutorial + // we will put these fundamentals into practice and build on them further. + // + // Accompanying the description of the algorithmic steps to use the frameworks + // will be a simplified view as to what they *might* be doing in the + // background. This description will be very much one designed to aid + // understanding, and the reader is encouraged to view the @ref auto_symb_diff + // module documentation for a far more formal description into how these tools + // actually work. + // + // @sect4{An analytical function} + namespace SimpleExample + { + // In order to convince the reader that these tools are indeed useful in + // practice, let us choose a function for which it is not too difficult to + // compute the analytical derivatives by hand. It's just sufficiently + // complicated to make you think about whether or not you truly want to go + // through with this exercise, and might also make you question whether you + // are completely sure that your calculations and implementation for its + // derivatives are correct. The point, of course, is that differentiation of + // functions is in a sense relatively formulaic and should be something + // computers are good at -- if we could build on existing software that + // understands the rules, we wouldn't have to bother with doing it + // ourselves. + // + // We choose the two variable trigonometric function + // $f(x,y) = \cos\left(\frac{y}{x}\right)$ for this purpose. Notice that + // this function is templated on the number type. This is done because we + // can often (but not always) use special auto-differentiable and symbolic + // types as drop-in replacements for real or complex valued types, and these + // will then perform some elementary calculations, such as evaluate a + // function value along with its derivatives. We will exploit that property + // and make sure that we need only define our function once, and then it can + // be re-used in whichever context we wish to perform differential + // operations on it. + template + NumberType f(const NumberType &x, const NumberType &y) + { + return std::cos(y / x); + } + + // Rather than revealing this function's derivatives immediately, we'll + // forward declare functions that return them and defer their definition to + // later. As implied by the function names, they respectively return + // the derivatives $\frac{df(x,y)}{dx}$: + double df_dx(const double x, const double y); + + // $\frac{df(x,y)}{dy}$: + double df_dy(const double x, const double y); + + // $\frac{d^{2}f(x,y)}{dx^{2}}$: + double d2f_dx_dx(const double x, const double y); + + // $\frac{d^{2}f(x,y)}{dx dy}$: + double d2f_dx_dy(const double x, const double y); + + // $\frac{d^{2}f(x,y)}{dy dx}$: + double d2f_dy_dx(const double x, const double y); + + // and, lastly, $\frac{d^{2}f(x,y)}{dy^{2}}$: + double d2f_dy_dy(const double x, const double y); + + + // @sect4{Computing derivatives using automatic differentiation} + + // To begin, we'll use AD as the tool to automatically + // compute derivatives for us. We will evaluate the function with the + // arguments `x` and `y`, and expect the resulting value and all of the + // derivatives to match to within the given tolerance. + void + run_and_verify_ad(const double x, const double y, const double tol = 1e-12) + { + // Our function $f(x,y)$ is a scalar-valued function, with arguments that + // represent the typical input variables that one comes across in + // algebraic calculations or tensor calculus. For this reason, the + // Differentiation::AD::ScalarFunction class is the appropriate wrapper + // class to use to do the computations that we require. (As a point of + // comparison, if the function arguments represented finite element cell + // degrees-of-freedom, we'd want to treat them differently.) The spatial + // dimension of the problem is irrelevant since we have no vector- or + // tensor-valued arguments to accommodate, so the `dim` template argument + // is arbitrarily assigned a value of 1. The second template argument + // stipulates which AD framework will be used (deal.II has support for + // several external AD frameworks), and what the underlying number type + // provided by this framework is to be employed. This number type + // influences the maximum order of the differential operation, and the + // underlying algorithms that are used to compute them. Given its template + // nature, this choice is a compile-time decision because many (but not + // all) of the AD libraries exploit compile-time meta-programming to + // implement these special number types in an efficient manner. The third + // template parameter states what the result type is; in our case, we're + // working with `double`s. + constexpr unsigned int dim = 1; + constexpr Differentiation::AD::NumberTypes ADTypeCode = + Differentiation::AD::NumberTypes::sacado_dfad_dfad; + using ADHelper = + Differentiation::AD::ScalarFunction; + + // It is necessary that we pre-register with our @p ADHelper class how many + // arguments (what we will call "independent variables") the function + // $f(x,y)$ has. Those arguments are `x` and `y`, so obviously there + // are two of them. + constexpr unsigned int n_independent_variables = 2; + + // We now have sufficient information to create and initialize an + // instance of the helper class. We can also get the concrete + // number type that will be used in all subsequent calculations. + // This is useful, because we can write everything from here on by + // referencing this type, and if we ever want to change the framework + // used, or number type (e.g., if we need more differential operations) + // then we need only adjust the `ADTypeCode` template parameter. + ADHelper ad_helper(n_independent_variables); + using ADNumberType = typename ADHelper::ad_type; + + // The next step is to register the numerical values of the independent + // variables with the helper class. This is done because the function + // and its derivatives will be evaluated for exactly these arguments. + // Since we register them in the order `{x,y}`, the variable `x` will + // be assigned component number `0`, and `y` will be component `1` + // -- a detail that will be used in the next few lines. + ad_helper.register_independent_variables({x, y}); + + // We now ask for the helper class to give to us the independent variables + // with their auto-differentiable representation. These are termed + // "sensitive variables", because from this point on any operations that + // we do with the components `independent_variables_ad` are tracked and + // recorded by the AD framework, and will be considered + // when we ask for the derivatives of something that they're used to + // compute. What the helper returns is a `vector` of auto-differentiable + // numbers, but we can be sure that the zeroth element represents `x` + // and the first element `y`. Just to make completely sure that there's + // no ambiguity of what number type these variables are, we suffix all of + // the auto-differentiable variables with `ad`. + const std::vector independent_variables_ad = + ad_helper.get_sensitive_variables(); + const ADNumberType &x_ad = independent_variables_ad[0]; + const ADNumberType &y_ad = independent_variables_ad[1]; + + // We can immediately pass in our sensitive representation of the + // independent variables to our templated function that computes + // $f(x,y)$. + // This also returns an auto-differentiable number. + const ADNumberType f_ad = f(x_ad, y_ad); + + // So now the natural question to ask is what we have actually just + // computed by passing these special `x_ad` and `y_ad` variables to the + // function `f`, instead of the original `double` variables `x` and `y`? + // In other words, how is all of this related to the computation of the + // derivatives that we were wanting to determine? Or, more concisely: What + // is so special about this returned `ADNumberType` object that gives it + // the ability to magically return derivatives? + // + // In essence, how this *could* be done is the following: + // This special number can be viewed as a data structure that stores the + // function value, and the prescribed number of derivatives. For a + // once-differentiable number expecting two arguments, it might look like + // this: + // + // @code + // struct ADNumberType + // { + // double value; // The value of the object + // double derivatives[2]; // Array of derivatives of the object with + // respect + // // to x and y + // }; + // @endcode + // + // For our independent variable `x_ad`, the starting value of `x_ad.value` + // would simply be its assigned value (i.e., the real value of that this + // variable represents). The derivative `x_ad.derivatives[0]` would be + // initialized to `1`, since `x` is the zeroth independent variable and + // $\frac{d(x)}{dx} = 1$. The derivative `x.derivatives[1]` would be + // initialized to zero, since the first independent variable is `y` and + // $\frac{d(x)}{dy} = 0$. + // + // For the function derivatives to be meaningful, we must assume that not + // only is this function differentiable in an analytical sense, but that + // it is also differentiable at the evaluation point `x,y`. + // We can exploit both of these assumptions: when we use this number type + // in mathematical operations, the AD framework *could* + // overload the operations (e.g., `%operator+()`, `%operator*()` as well + // as `%sin()`, `%exp()`, etc.) such that the returned result has the + // expected value. At the same time, it would then compute the derivatives + // through the knowledge of exactly what function is being overloaded and + // rigorous application of the chain-rule. So, the `%sin()` function + // (with its argument `a` itself being a function of the independent + // variables `x` and `y`) *might* be defined as follows: + // + // @code + // ADNumberType sin(const ADNumberType &a) + // { + // ADNumberType output; + // + // // For the input argument "a", "a.value" is simply its value. + // output.value = sin(a.value); + // + // // We know that the derivative of sin(a) is cos(a), but we need + // // to also consider the chain rule and that the input argument + // // `a` is also differentiable with respect to the original + // // independent variables `x` and `y`. So `a.derivatives[0]` + // // and `a.derivatives[1]` respectively represent the partial + // // derivatives of `a` with respect to its inputs `x` and `y`. + // output.derivatives[0] = cos(a.value)*a.derivatives[0]; + // output.derivatives[1] = cos(a.value)*a.derivatives[1]; + // + // return output; + // }; + // @endcode + // + // All of that could of course also be done for second and even higher + // order derivatives. + // + // So it is now clear that with the above representation the + // `ADNumberType` is carrying around some extra data that represents the + // various derivatives of differentiable functions with respect to the + // original (sensitive) independent variables. It should therefore be + // noted that there is computational overhead associated with using them + // (as we compute extra functions when doing derivative computations) as + // well as memory overhead in storing these results. So the prescribed + // number of levels of differential operations should ideally be kept to a + // minimum to limit computational cost. We could, for instance, have + // computed the first derivatives ourself and then have used the + // Differentiation::AD::VectorFunction helper class to determine the + // gradient of the collection of dependent functions, which would be the + // second derivatives of the original scalar function. + // + // It is also worth noting that because the chain rule is indiscriminately + // applied and we only see the beginning and end-points of the calculation + // `{x,y}` $\rightarrow$ `f(x,y)`, we will only ever be able to query + // the total derivatives of `f`; the partial derivatives + // (`a.derivatives[0]` and `a.derivatives[1]` in the above example) are + // intermediate values and are hidden from us. + + // Okay, since we now at least have some idea as to exactly what `f_ad` + // represents and what is encoded within it, let's put all of that to + // some actual use. To gain access to those hidden derivative results, + // we register the final result with the helper class. After this point, + // we can no longer change the value of `f_ad` and have those changes + // reflected in the results returned by the helper class. + ad_helper.register_dependent_variable(f_ad); + + // The next step is to extract the derivatives (specifically, the function + // gradient and Hessian). To do so we first create some temporary data + // structures (with the result type `double`) to store the derivatives + // (noting that all derivatives are returned at once, and not + // individually)... + Vector Df(ad_helper.n_dependent_variables()); + FullMatrix D2f(ad_helper.n_dependent_variables(), + ad_helper.n_independent_variables()); + + // ... and we then request that the helper class compute these + // derivatives, and the function value itself. And that's it. We have + // everything that we were aiming to get. + const double computed_f = ad_helper.compute_value(); + ad_helper.compute_gradient(Df); + ad_helper.compute_hessian(D2f); + + // We can convince ourselves that the AD framework is + // correct by comparing it to the analytical solution. (Or, if you're + // like the author, you'll be doing the opposite and will rather verify + // that your implementation of the analytical solution is correct!) + AssertThrow(std::abs(f(x, y) - computed_f) < tol, + ExcMessage(std::string("Incorrect value computed for f. ") + + std::string("Hand-calculated value: ") + + Utilities::to_string(f(x, y)) + + std::string(" ; ") + + std::string("Value computed by AD: ") + + Utilities::to_string(computed_f))); + + // Because we know the ordering of the independent variables, we know + // which component of the gradient relates to which derivative... + const double computed_df_dx = Df[0]; + const double computed_df_dy = Df[1]; + + AssertThrow(std::abs(df_dx(x, y) - computed_df_dx) < tol, + ExcMessage( + std::string("Incorrect value computed for df/dx. ") + + std::string("Hand-calculated value: ") + + Utilities::to_string(df_dx(x, y)) + std::string(" ; ") + + std::string("Value computed by AD: ") + + Utilities::to_string(computed_df_dx))); + AssertThrow(std::abs(df_dy(x, y) - computed_df_dy) < tol, + ExcMessage( + std::string("Incorrect value computed for df/dy. ") + + std::string("Hand-calculated value: ") + + Utilities::to_string(df_dy(x, y)) + std::string(" ; ") + + std::string("Value computed by AD: ") + + Utilities::to_string(computed_df_dy))); + + // ... and similar for the Hessian. + const double computed_d2f_dx_dx = D2f[0][0]; + const double computed_d2f_dx_dy = D2f[0][1]; + const double computed_d2f_dy_dx = D2f[1][0]; + const double computed_d2f_dy_dy = D2f[1][1]; + + AssertThrow(std::abs(d2f_dx_dx(x, y) - computed_d2f_dx_dx) < tol, + ExcMessage( + std::string("Incorrect value computed for d2f/dx_dx. ") + + std::string("Hand-calculated value: ") + + Utilities::to_string(d2f_dx_dx(x, y)) + std::string(" ; ") + + std::string("Value computed by AD: ") + + Utilities::to_string(computed_d2f_dx_dx))); + AssertThrow(std::abs(d2f_dx_dy(x, y) - computed_d2f_dx_dy) < tol, + ExcMessage( + std::string("Incorrect value computed for d2f/dx_dy. ") + + std::string("Hand-calculated value: ") + + Utilities::to_string(d2f_dx_dy(x, y)) + std::string(" ; ") + + std::string("Value computed by AD: ") + + Utilities::to_string(computed_d2f_dx_dy))); + AssertThrow(std::abs(d2f_dy_dx(x, y) - computed_d2f_dy_dx) < tol, + ExcMessage( + std::string("Incorrect value computed for d2f/dy_dx. ") + + std::string("Hand-calculated value: ") + + Utilities::to_string(d2f_dy_dx(x, y)) + std::string(" ; ") + + std::string("Value computed by AD: ") + + Utilities::to_string(computed_d2f_dy_dx))); + AssertThrow(std::abs(d2f_dy_dy(x, y) - computed_d2f_dy_dy) < tol, + ExcMessage( + std::string("Incorrect value computed for d2f/dy_dy. ") + + std::string("Hand-calculated value: ") + + Utilities::to_string(d2f_dy_dy(x, y)) + std::string(" ; ") + + std::string("Value computed by AD: ") + + Utilities::to_string(computed_d2f_dy_dy))); + } + + // That's pretty great. There wasn't too much work involved in computing + // second-order derivatives of this trigonometric function. + + // @sect4{Hand-calculated derivatives of the analytical solution} + + // Since we now know how much "implementation effort" it takes to have the + // AD framework compute those derivatives for us, let's + // compare that to the same computed by hand and implemented in several + // stand-alone functions. + + // Here are the two first derivatives of $f(x,y) = + // \cos\left(\frac{y}{x}\right)$: + // + // $\frac{df(x,y)}{dx} = \frac{y}{x^2} \sin\left(\frac{y}{x}\right)$ + double df_dx(const double x, const double y) + { + Assert(x != 0.0, ExcDivideByZero()); + return y * std::sin(y / x) / (x * x); + } + + // $\frac{df(x,y)}{dx} = -\frac{1}{x} \sin\left(\frac{y}{x}\right)$ + double df_dy(const double x, const double y) + { + return -std::sin(y / x) / x; + } + + // And here are the four second derivatives of $f(x,y)$: + // + // $\frac{d^{2}f(x,y)}{dx^{2}} = -\frac{y}{x^4} (2x + // \sin\left(\frac{y}{x}\right) + y \cos\left(\frac{y}{x}\right))$ + double d2f_dx_dx(const double x, const double y) + { + return -y * (2 * x * std::sin(y / x) + y * std::cos(y / x)) / + (x * x * x * x); + } + + // $\frac{d^{2}f(x,y)}{dx dy} = \frac{1}{x^3} (x + // \sin\left(\frac{y}{x}\right) + y \cos\left(\frac{y}{x}\right))$ + double d2f_dx_dy(const double x, const double y) + { + return (x * std::sin(y / x) + y * std::cos(y / x)) / (x * x * x); + } + + // $\frac{d^{2}f(x,y)}{dy dx} = \frac{1}{x^3} (x + // \sin\left(\frac{y}{x}\right) + y \cos\left(\frac{y}{x}\right))$ (as + // expected, on the basis of [Schwarz's + // theorem](https://en.wikipedia.org/wiki/Symmetry_of_second_derivatives)) + double d2f_dy_dx(const double x, const double y) + { + return (x * std::sin(y / x) + y * std::cos(y / x)) / (x * x * x); + } + + // $\frac{d^{2}f(x,y)}{dy^{2}} = -\frac{1}{x^2} + // \cos\left(\frac{y}{x}\right)$ + double d2f_dy_dy(const double x, const double y) + { + return -(std::cos(y / x)) / (x * x); + } + + // Hmm... there's a lot of places in the above where we could have + // introduced an error in the above, especially when it comes to employing + // the chain rule. Although they're no silver bullet, at the very least + // these AD frameworks can serve as a verification tool to make sure that we + // haven't made any errors (either by calculation or by implementation) that + // would negatively affect our results. + + // The point of this example of course is that we might have + // chosen a relatively simple function $f(x,y)$ for which we can + // hand-verify that the derivatives the AD framework computed is + // correct. But the AD framework didn't care that the function was + // simple: It could have been a much much more convoluted + // expression, or could have depended on more than two variables, + // and it would still have been able to compute the derivatives -- + // the only difference would have been that *we* wouldn't have + // been able to come up with the derivatives any more to verify + // correctness of the AD framework. + + + + // @sect4{Computing derivatives using symbolic differentiation} + + // We'll now repeat the same exercise using symbolic differentiation. The + // term "symbolic differentiation" is a little bit misleading because + // differentiation is just one tool that the Computer Algebra System (CAS) + // (i.e., the symbolic framework) provides. Nevertheless, in the context + // of finite element modeling and applications it is the most common use + // of a CAS and will therefore be the one that we'll focus on. + // Once more, we'll supply the argument values `x` and `y` with which to + // evaluate our function $f(x,y) = \cos\left(\frac{y}{x}\right)$ and its + // derivatives, and a tolerance with which to test the correctness of the + // returned results. + void + run_and_verify_sd(const double x, const double y, const double tol = 1e-12) + { + // The first step that we need to take is to form the symbolic variables + // that represent the function arguments that we wish to differentiate + // with respect to. Again, these will be the independent variables for + // our problem and as such are, in some sense, primitive variables that + // have no dependencies on any other variable. We create these types of + // (independent) variables by initializing a symbolic type + // Differentiation::SD::Expression, which is a wrapper to a set of classes + // used by the symbolic framework, with a unique identifier. On this + // occasion it makes sense that this identifier, a `std::string`, be + // simply `"x"` for the $x$ argument, and likewise `"y"` for the $y$ + // argument to the dependent function. Like before, we'll suffix symbolic + // variable names with `sd` so that we can clearly see which variables are + // symbolic (as opposed to numeric) in nature. + const Differentiation::SD::Expression x_sd("x"); + const Differentiation::SD::Expression y_sd("y"); + + // Using the templated function that computes $f(x,y)$, we can pass + // these independent variables as arguments to the function. The returned + // result will be another symbolic type that represents the sequence of + // operations used to compute $\cos\left(\frac{y}{x}\right)$. + const Differentiation::SD::Expression f_sd = f(x_sd, y_sd); + + // At this point it is legitimate to print out the expression `f_sd`, and + // if we did so + // @code + // std::cout << "f(x,y) = " << f_sd << std::endl; + // @endcode + // we would see `f(x,y) = cos(y/x)` printed to the console. + // + // You might notice that we've constructed our symbolic function `f_sd` + // with no context as to how we might want to use it: In contrast to the + // AD approach shown above, what we were returned from calling + // `f(x_sd, y_sd)` is not the evaluation of the function `f` at some + // specific point, but is in fact a symbolic representation of the + // evaluation at a generic, as yet undetermined, point. This is one of the + // key points that makes symbolic frameworks (the CAS) different from + // automatic differentiation frameworks. Each of the variables `x_sd` and + // `y_sd`, and even the composite dependent function `f_sd`, are in some + // sense respectively "placeholders" for numerical values and a + // composition of operations. In fact, the individual components that are + // used to compose the function are also placeholders. The sequence of + // operations are encoded into in a tree-like data structure (conceptually + // simlar to an [abstract syntax + // tree](https://en.wikipedia.org/wiki/Abstract_syntax_tree)). + // + // Once we form these data structures we can defer any operations that we + // might want to do with them until some later time. Each of these + // placeholders represents something, but we have the opportunity to + // define or redefine what they represent at any convenient point in time. + // So for this particular problem it makes sense that we somehow want to + // associate "x" and "y" with *some* numerical value (with type yet to be + // determined), but we could conceptually (and if it made sense) assign + // the ratio "y/x" a value instead of the variables "x" and "y" + // individually. We could also associate with "x" or "y" some other + // symbolic function `g(a,b)`. Any of these operations involves + // manipulating the recorded tree of operations, and substituting the + // salient nodes on the tree (and that nodes' subtree) with something + // else. The key word here is "substitution", and indeed there are many + // functions in the Differentiation::SD namespace that have this word + // in their names. + // + // This capability makes the framework entirely generic. + // In the context of finite element simulations, the types of operations + // that we would typically perform with our symbolic types are + // function composition, differentiation, substitution (partial or + // complete), and evaluation (i.e., conversion of the symbolic type to its + // numerical counterpart). But should you need it, a CAS is often capable + // of more than just this: It could be forming anti-derivatives + // (integrals) of functions, perform simplifications on the expressions + // that form a function (e.g., replace $(\sin a)^2 + (\cos a)^2$ by + // $1$; or, more simply: if the function did an operation like `1+2`, a + // CAS could replace it by `3`), and so forth: The *expression* that a + // variable represents is obtained from how the function $f$ is + // implemented, but a CAS can do with it whatever its functionality + // happens to be. + // + // Specifically, to compute the symbolic representation of the first + // derivatives of the dependent function with respect to its individual + // independent variables, we use the + // Differentiation::SD::Expression::differentiate() function with the + // independent variable given as its argument. Each call will cause the + // CAS to go through the tree of operations that compose `f_sd` and + // differentiate each node of the expression tree with respect to the + // given symbolic argument. + const Differentiation::SD::Expression df_dx_sd = f_sd.differentiate(x_sd); + const Differentiation::SD::Expression df_dy_sd = f_sd.differentiate(y_sd); + + // To compute the symbolic representation of the second derivatives, we + // simply differentiate the first derivatives with respect to the + // independent variables. So to compute a higher order derivative, we + // first need to compute the lower order derivative. + // (As the return type of the call to `differentiate()` is an expression, + // we could in principal execute double differentiation directly from the + // scalar by chaining two calls together. But this is unnecessary in this + // particular case, since we have the intermediate results at hand.) + const Differentiation::SD::Expression d2f_dx_dx_sd = + df_dx_sd.differentiate(x_sd); + const Differentiation::SD::Expression d2f_dx_dy_sd = + df_dx_sd.differentiate(y_sd); + const Differentiation::SD::Expression d2f_dy_dx_sd = + df_dy_sd.differentiate(x_sd); + const Differentiation::SD::Expression d2f_dy_dy_sd = + df_dy_sd.differentiate(y_sd); + // Printing the expressions for the first and second derivatives, as + // computed by the CAS, using the statements + // @code + // std::cout << "df_dx_sd: " << df_dx_sd << std::endl; + // std::cout << "df_dy_sd: " << df_dy_sd << std::endl; + // std::cout << "d2f_dx_dx_sd: " << d2f_dx_dx_sd << std::endl; + // std::cout << "d2f_dx_dy_sd: " << d2f_dx_dy_sd << std::endl; + // std::cout << "d2f_dy_dx_sd: " << d2f_dy_dx_sd << std::endl; + // std::cout << "d2f_dy_dy_sd: " << d2f_dy_dy_sd << std::endl; + // @endcode + // renders the following output: + // @code{.sh} + // df_dx_sd: y*sin(y/x)/x**2 + // df_dy_sd: -sin(y/x)/x + // d2f_dx_dx_sd: -y**2*cos(y/x)/x**4 - 2*y*sin(y/x)/x**3 + // d2f_dx_dy_sd: sin(y/x)/x**2 + y*cos(y/x)/x**3 + // d2f_dy_dx_sd: sin(y/x)/x**2 + y*cos(y/x)/x**3 + // d2f_dy_dy_sd: -cos(y/x)/x**2 + // @endcode + // This compares favorably to the analytical expressions for these + // derivatives that were presented earlier. + + // Now that we have formed the symbolic expressions for the function and + // its derivatives, we want to evaluate them for the numeric values for + // the main function arguments `x` and `y`. To accomplish this, we + // construct a *substitution map*, which maps the symbolic values to their + // numerical counterparts. + const Differentiation::SD::types::substitution_map substitution_map = + Differentiation::SD::make_substitution_map( + std::pair{x_sd, x}, + std::pair{y_sd, y}); + + // The last step in the process is to convert all symbolic variables and + // operations into numerical values, and produce the numerical result of + // this operation. To do this we combine the substitution map with the + // symbolic variable in the step we have already mentioned above: + // "substitution". + // + // Once we pass this substitution map to the CAS, it will + // substitute each instance of the symbolic variable (or, more generally, + // sub-expression) with its numerical counterpart and then propagate these + // results up the operation tree, simplifying each node on the tree if + // possible. If the tree is reduced to a single value (i.e., we have + // substituted all of the independent variables with their numerical + // counterpart) then the evaluation is complete. + // + // Due to the strongly-typed nature of C++, we need to instruct the CAS to + // convert its representation of the result into an intrinsic data type + // (in this case a `double`). This is the "evaluation" step, and through + // the template type we define the return type of this process. + // Conveniently, these two steps can be done at once if we are certain + // that we've performed a full substitution. + const double computed_f = + f_sd.substitute_and_evaluate(substitution_map); + + AssertThrow(std::abs(f(x, y) - computed_f) < tol, + ExcMessage(std::string("Incorrect value computed for f. ") + + std::string("Hand-calculated value: ") + + Utilities::to_string(f(x, y)) + + std::string(" ; ") + + std::string("Value computed by AD: ") + + Utilities::to_string(computed_f))); + + // We can do the same for the first derivatives... + const double computed_df_dx = + df_dx_sd.substitute_and_evaluate(substitution_map); + const double computed_df_dy = + df_dy_sd.substitute_and_evaluate(substitution_map); + + AssertThrow(std::abs(df_dx(x, y) - computed_df_dx) < tol, + ExcMessage( + std::string("Incorrect value computed for df/dx. ") + + std::string("Hand-calculated value: ") + + Utilities::to_string(df_dx(x, y)) + std::string(" ; ") + + std::string("Value computed by AD: ") + + Utilities::to_string(computed_df_dx))); + AssertThrow(std::abs(df_dy(x, y) - computed_df_dy) < tol, + ExcMessage( + std::string("Incorrect value computed for df/dy. ") + + std::string("Hand-calculated value: ") + + Utilities::to_string(df_dy(x, y)) + std::string(" ; ") + + std::string("Value computed by AD: ") + + Utilities::to_string(computed_df_dy))); + + // ... and the second derivatives. + // Notice that we can reuse the same substitution map for each of these + // operations because we wish to evaluate all of these functions for the + // same values of `x` and `y`. Modifying the values in the substitution + // map renders the result of same symbolic expression evaluated with + // different values being assigned to the independent variables. + // We could also happily have each variable represent a real value in + // one pass, and a complex value in the next. + const double computed_d2f_dx_dx = + d2f_dx_dx_sd.substitute_and_evaluate(substitution_map); + const double computed_d2f_dx_dy = + d2f_dx_dy_sd.substitute_and_evaluate(substitution_map); + const double computed_d2f_dy_dx = + d2f_dy_dx_sd.substitute_and_evaluate(substitution_map); + const double computed_d2f_dy_dy = + d2f_dy_dy_sd.substitute_and_evaluate(substitution_map); + + AssertThrow(std::abs(d2f_dx_dx(x, y) - computed_d2f_dx_dx) < tol, + ExcMessage( + std::string("Incorrect value computed for d2f/dx_dx. ") + + std::string("Hand-calculated value: ") + + Utilities::to_string(d2f_dx_dx(x, y)) + std::string(" ; ") + + std::string("Value computed by SD: ") + + Utilities::to_string(computed_d2f_dx_dx))); + AssertThrow(std::abs(d2f_dx_dy(x, y) - computed_d2f_dx_dy) < tol, + ExcMessage( + std::string("Incorrect value computed for d2f/dx_dy. ") + + std::string("Hand-calculated value: ") + + Utilities::to_string(d2f_dx_dy(x, y)) + std::string(" ; ") + + std::string("Value computed by SD: ") + + Utilities::to_string(computed_d2f_dx_dy))); + AssertThrow(std::abs(d2f_dy_dx(x, y) - computed_d2f_dy_dx) < tol, + ExcMessage( + std::string("Incorrect value computed for d2f/dy_dx. ") + + std::string("Hand-calculated value: ") + + Utilities::to_string(d2f_dy_dx(x, y)) + std::string(" ; ") + + std::string("Value computed by SD: ") + + Utilities::to_string(computed_d2f_dy_dx))); + AssertThrow(std::abs(d2f_dy_dy(x, y) - computed_d2f_dy_dy) < tol, + ExcMessage( + std::string("Incorrect value computed for d2f/dy_dy. ") + + std::string("Hand-calculated value: ") + + Utilities::to_string(d2f_dy_dy(x, y)) + std::string(" ; ") + + std::string("Value computed by SD: ") + + Utilities::to_string(computed_d2f_dy_dy))); + } + + + // @sect4{The SimpleExample::run() function} + + // The function used to drive these initial examples is straightforward. + // We'll arbitrarily choose some values at which to evaluate the function + // (although knowing that `x = 0` is not permissible), and then pass these + // values to the functions that use the AD and SD frameworks. + void run() + { + const double x = 1.23; + const double y = 0.91; + + std::cout << "Simple example using automatic differentiation..." + << std::endl; + run_and_verify_ad(x, y); + std::cout << "... all calculations are correct!" << std::endl; + + std::cout << "Simple example using symbolic differentiation." + << std::endl; + run_and_verify_sd(x, y); + std::cout << "... all calculations are correct!" << std::endl; + } + + } // namespace SimpleExample + + + // @sect3{A more complex example: Using automatic and symbolic differentiation to compute derivatives at continuum points} + + // Now that we've introduced the principles behind automatic and symbolic + // differentiation, we'll put them into action by formulating two coupled + // magneto-mechanical constitutive laws: one that is rate-independent, and + // another that exhibits rate-dependent behavior. + // + // As you will recall from the introduction, the material + // constitutive laws we will consider are far more complicated than + // the simple example above. This is not just because of the form of + // the function $\psi_{0}$ that we will consider, but in particular + // because $\psi_{0}$ doesn't just depend on two scalar variables, but + // instead on a whole bunch of *tensors*, each with several + // components. In some cases, these are *symmetric* tensors, for + // which only a subset of components is in fact independent, and has + // to think about what it actually means to compute a derivative + // such as $\frac{\partial\psi_{0}}{\partial \mathbf{C}}$ where $\mathbf + // C$ is a symmetric tensor. How all of this will work will, + // hopefully, become clear below. It will also become clear that + // doing this by hand is going to be, at the very best, *exceedingly* + // *tedious* and, at worst, riddled with hard-to-find bugs. + namespace CoupledConstitutiveLaws + { + // @sect4{Constitutive parameters} + + // We start with a description of the various material parameters + // that appear in the description of the energy function $\psi_{0}$. + // + // The ConstitutiveParameters class is used to hold these values. + // Values for all parameters (both constitutive and rheological) are taken + // from @cite Pelteret2018a, and are given values that produce a + // constitutive response that is broadly representative of a real, + // laboratory-made magneto-active polymer, though the specific values used + // here are of no consequence to the purpose of this program of course. + // + // The first four constitutive parameters respectively represent + // - the elastic shear modulus $\mu_{e}$, + // - the elastic shear modulus at magnetic saturation $\mu_{e}^{\infty}$, + // - the saturation magnetic field strength for the elastic shear + // modulus $h_{e}^{\text{sat}}$, and + // - the Poisson ratio $\nu$. + class ConstitutiveParameters : public ParameterAcceptor + { + public: + ConstitutiveParameters(); + + double mu_e = 30.0e3; + double mu_e_inf = 250.0e3; + double mu_e_h_sat = 212.2e3; + double nu_e = 0.49; + + // The next four, which only pertain to the rate-dependent material, are + // parameters for + // - the viscoelastic shear modulus $\mu_{v}$, + // - the viscoelastic shear modulus at magnetic saturation + // $\mu_{v}^{\infty}$, + // - the saturation magnetic field strength for the viscoelastic + // shear modulus $h_{v}^{\text{sat}}$, and + // - the characteristic relaxation time $\tau$. + double mu_v = 20.0e3; + double mu_v_inf = 35.0e3; + double mu_v_h_sat = 92.84e3; + double tau_v = 0.6; + + // The last parameter is the relative magnetic permeability $\mu_{r}$. + double mu_r = 6.0; + + bool initialized = false; + }; + + // The parameters are initialized through the ParameterAcceptor + // framework, which is discussed in detail in step-60. + ConstitutiveParameters::ConstitutiveParameters() + : ParameterAcceptor("/Coupled Constitutive Laws/Constitutive Parameters/") + { + add_parameter("Elastic shear modulus", mu_e); + add_parameter("Elastic shear modulus at magnetic saturation", mu_e_inf); + add_parameter( + "Saturation magnetic field strength for elastic shear modulus", + mu_e_h_sat); + add_parameter("Poisson ratio", nu_e); + + add_parameter("Viscoelastic shear modulus", mu_v); + add_parameter("Viscoelastic shear modulus at magnetic saturation", + mu_v_inf); + add_parameter( + "Saturation magnetic field strength for viscoelastic shear modulus", + mu_v_h_sat); + add_parameter("Characteristic relaxation time", tau_v); + + add_parameter("Relative magnetic permeability", mu_r); + + parse_parameters_call_back.connect([&]() { initialized = true; }); + } + + + // @sect4{Constitutive laws: Base class} + + // Since we'll be formulating two constitutive laws for the same class of + // materials, it makes sense to define a base class that ensures a unified + // interface to them. + // + // The class declaration starts with the constructor that will + // accept the set of constitutive parameters that, in conjunction + // with the material law itself, dictate the material response. + template + class Coupled_Magnetomechanical_Constitutive_Law_Base + { + public: + Coupled_Magnetomechanical_Constitutive_Law_Base( + const ConstitutiveParameters &constitutive_parameters); + + // Instead of computing and returning the kinetic variables or their + // linearization at will, we'll calculate and store these values within a + // single method. These cached results will then be returned upon request. + // We'll defer the precise explanation as to why we'd want to do this to + // a later stage. What is important for now is to see that this function + // accepts all of the field variables, namely the magnetic field vector + // $\boldsymbol{\mathbb{H}}$ and right Cauchy-Green deformation tensor + // $\mathbf{C}$, as well as the time discretizer. These, in addition to + // the @p constitutive_parameters, are all the fundamental quantities that + // are required to compute the material response. + virtual void update_internal_data(const SymmetricTensor<2, dim> &C, + const Tensor<1, dim> & H, + const DiscreteTime &time) = 0; + + // The next few functions provide the interface to probe the material + // response due subject to the applied deformation and magnetic loading. + // + // Since the class of materials can be expressed in terms of a free energy + // $\psi_{0}$, we can compute that... + virtual double get_psi() const = 0; + + // ... as well as the two kinetic quantities: + // - the magnetic induction vector $\boldsymbol{\mathbb{B}}$, and + // - the total Piola-Kirchhoff stress tensor $\mathbf{S}^{\text{tot}}$ + virtual Tensor<1, dim> get_B() const = 0; + + virtual SymmetricTensor<2, dim> get_S() const = 0; + + // ... and the linearization of the kinetic quantities, which are: + // - the magnetostatic tangent tensor $\mathbb{D}$, + // - the total referential magnetoelastic coupling tensor + // $\mathfrak{P}^{\text{tot}}$, and + // - the total referential elastic tangent tensor + // $\mathcal{H}^{\text{tot}}$. + virtual SymmetricTensor<2, dim> get_DD() const = 0; + + virtual Tensor<3, dim> get_PP() const = 0; + + virtual SymmetricTensor<4, dim> get_HH() const = 0; + + // We'll also define a method that provides a mechanism for this class + // instance to do any additional tasks before moving on to the next + // timestep. Again, the reason for doing this will become clear a little + // later. + virtual void update_end_of_timestep(){}; + + // In the `protected` part of the class, + // we store a reference to an instance of the constitutive parameters + // that govern the material response. + // For convenience, we also define some functions that return + // various constitutive parameters (both explicitly defined, as well + // as calculated). + // + // The parameters related to the elastic response of the material are, + // in order: + // - the elastic shear modulus, + // - the elastic shear modulus at saturation magnetic field, + // - the saturation magnetic field strength for the elastic shear + // modulus, + // - the Poisson ratio, + // - the Lamé parameter, and + // - the bulk modulus. + protected: + const ConstitutiveParameters &constitutive_parameters; + + double get_mu_e() const; + + double get_mu_e_inf() const; + + double get_mu_e_h_sat() const; + + double get_nu_e() const; + + double get_lambda_e() const; + + double get_kappa_e() const; + + // The parameters related to the elastic response of the material are, + // in order: + // - the viscoelastic shear modulus, + // - the viscoelastic shear modulus at magnetic saturation, + // - the saturation magnetic field strength for the viscoelastic + // shear modulus, and + // - the characteristic relaxation time. + double get_mu_v() const; + + double get_mu_v_inf() const; + + double get_mu_v_h_sat() const; + + double get_tau_v() const; + + // The parameters related to the magnetic response of the material are, + // in order: + // - the relative magnetic permeability, and + // - the magnetic permeability constant $\mu_{0}$ (not really a material + // constant, + // but rather a universal constant that we'll group here for + // simplicity). + // + // We'll also implement a function that returns the + // timestep size from the time discretizion. + double get_mu_r() const; + + constexpr double get_mu_0() const; + double get_delta_t(const DiscreteTime &time) const; + }; + + + + // In the following, let us start by implementing the several + // relatively trivial member functions of the class just defined: + template + Coupled_Magnetomechanical_Constitutive_Law_Base:: + Coupled_Magnetomechanical_Constitutive_Law_Base( + const ConstitutiveParameters &constitutive_parameters) + : constitutive_parameters(constitutive_parameters) + { + Assert(get_kappa_e() > 0, ExcInternalError()); + } + + + template + double + Coupled_Magnetomechanical_Constitutive_Law_Base::get_mu_e() const + { + return constitutive_parameters.mu_e; + } + + + template + double + Coupled_Magnetomechanical_Constitutive_Law_Base::get_mu_e_inf() const + { + return constitutive_parameters.mu_e_inf; + } + + + template + double + Coupled_Magnetomechanical_Constitutive_Law_Base::get_mu_e_h_sat() const + { + return constitutive_parameters.mu_e_h_sat; + } + + + template + double + Coupled_Magnetomechanical_Constitutive_Law_Base::get_nu_e() const + { + return constitutive_parameters.nu_e; + } + + + template + double + Coupled_Magnetomechanical_Constitutive_Law_Base::get_lambda_e() const + { + return 2.0 * get_mu_e() * get_nu_e() / (1.0 - 2.0 * get_nu_e()); + } + + + template + double + Coupled_Magnetomechanical_Constitutive_Law_Base::get_kappa_e() const + { + return (2.0 * get_mu_e() * (1.0 + get_nu_e())) / + (3.0 * (1.0 - 2.0 * get_nu_e())); + } + + + template + double + Coupled_Magnetomechanical_Constitutive_Law_Base::get_mu_v() const + { + return constitutive_parameters.mu_v; + } + + + template + double + Coupled_Magnetomechanical_Constitutive_Law_Base::get_mu_v_inf() const + { + return constitutive_parameters.mu_v_inf; + } + + + template + double + Coupled_Magnetomechanical_Constitutive_Law_Base::get_mu_v_h_sat() const + { + return constitutive_parameters.mu_v_h_sat; + } + + + template + double + Coupled_Magnetomechanical_Constitutive_Law_Base::get_tau_v() const + { + return constitutive_parameters.tau_v; + } + + + template + double + Coupled_Magnetomechanical_Constitutive_Law_Base::get_mu_r() const + { + return constitutive_parameters.mu_r; + } + + + template + constexpr double + Coupled_Magnetomechanical_Constitutive_Law_Base::get_mu_0() const + { + return 4.0 * numbers::PI * 1e-7; + } + + + template + double Coupled_Magnetomechanical_Constitutive_Law_Base::get_delta_t( + const DiscreteTime &time) const + { + return time.get_previous_step_size(); + } + + + // @sect4{Magnetoelastic constitutive law (using automatic differentiation)} + + // We'll being by considering a non-dissipative material, namely one that + // is governed by a magneto-hyperelastic constitutive law that exhibits + // stiffening when immersed in a magnetic field. As described in + // the introduction, the stored energy density function for such a material + // might be given by + // @f[ + // \psi_{0} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) + // = \frac{1}{2} \mu_{e} f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} \right) + // \left[ \text{tr}(\mathbf{C}) - d - 2 \ln (\text{det}(\mathbf{F})) + // \right] + // + \lambda_{e} \ln^{2} \left(\text{det}(\mathbf{F}) \right) + // - \frac{1}{2} \mu_{0} \mu_{r} \text{det}(\mathbf{F}) + // \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}} \right] + // @f] + // with + // @f[ + // f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} \right) + // = 1 + \left[ \frac{\mu_{e}^{\infty}}{\mu_{e}} - 1 \right] + // \tanh \left( 2 \frac{\boldsymbol{\mathbb{H}} \cdot + // \boldsymbol{\mathbb{H}}} + // {\left(h_{e}^{\text{sat}}\right)^{2}} \right) . + // @f] + // + // Now on to the class that implements this behavior. + // Since we expect that this class fully describes a single material, we'll + // mark it as "final" so that the inheritance tree terminated here. + // At the top of the class, we define the helper type that we will use in + // the AD computations for our scalar energy density function. Note that we + // expect it to return values of type `double`. We also have to specify the + // number of spatial dimensions, `dim`, so that the link between vector, + // tensor and symmetric tensor fields and the number of components that they + // contain may be established. The concrete `ADTypeCode` used for the + // ADHelper class will be provided as a template argument at the point where + // this class is actually used. + template + class Magnetoelastic_Constitutive_Law_AD final + : public Coupled_Magnetomechanical_Constitutive_Law_Base + { + using ADHelper = + Differentiation::AD::ScalarFunction; + using ADNumberType = typename ADHelper::ad_type; + + public: + Magnetoelastic_Constitutive_Law_AD( + const ConstitutiveParameters &constitutive_parameters); + + // Since the public interface to the base class is pure-`virtual`, here + // we'll declare that this class will override all of these base class + // methods. + virtual void update_internal_data(const SymmetricTensor<2, dim> &C, + const Tensor<1, dim> & H, + const DiscreteTime &) override; + + virtual double get_psi() const override; + + virtual Tensor<1, dim> get_B() const override; + + virtual SymmetricTensor<2, dim> get_S() const override; + + virtual SymmetricTensor<2, dim> get_DD() const override; + + virtual Tensor<3, dim> get_PP() const override; + + virtual SymmetricTensor<4, dim> get_HH() const override; + + // In the `private` part of the class, + // we need to define some extractors that will help us set independent + // variables and later get the computed values related to the dependent + // variables. If this class were to be used in the context of a finite + // element problem, then each of these extractors is (most likely) related + // to the gradient of a component of the solution field (in this case, + // displacement and magnetic scalar potential). As you can probably infer + // by now, here "C" denotes the right Cauchy-Green tensor and "H" denotes + // the magnetic field vector. + private: + const FEValuesExtractors::Vector H_components; + const FEValuesExtractors::SymmetricTensor<2> C_components; + + // This is an instance of the automatic differentiation helper that + // we'll set up to do all of the differential calculations related to + // the constitutive law... + ADHelper ad_helper; + + // ... and the following three member variables will store the output from + // the + // @p ad_helper. The @p ad_helper returns the derivatives with respect + // to all field variables at once, so we'll retain the full gradient + // vector and Hessian matrix. From that, we'll extract the individual + // entries that we're actually interested in. + double psi; + Vector Dpsi; + FullMatrix D2psi; + }; + + // When setting up the field component extractors, its completely arbitrary + // as to how they are ordered. But it is important that the extractors do + // not have overlapping indices. The total number of components of these + // extractors defines the number of independent variables that the + // @p ad_helper needs to track, and with respect to which we'll be taking + // derivatives. The resulting data structures @p Dpsi and @p D2psi must also + // be sized accordingly. Once the @p ad_helper is configured (its input + // argument being the total number of components of $\mathbf{C}$ and + // $\boldsymbol{\mathbb{H}}$), we can directly interrogate it as to how many + // independent variables it uses. + template + Magnetoelastic_Constitutive_Law_AD:: + Magnetoelastic_Constitutive_Law_AD( + const ConstitutiveParameters &constitutive_parameters) + : Coupled_Magnetomechanical_Constitutive_Law_Base( + constitutive_parameters) + , H_components(0) + , C_components(Tensor<1, dim>::n_independent_components) + , ad_helper(Tensor<1, dim>::n_independent_components + + SymmetricTensor<2, dim>::n_independent_components) + , psi(0.0) + , Dpsi(ad_helper.n_independent_variables()) + , D2psi(ad_helper.n_independent_variables(), + ad_helper.n_independent_variables()) + {} + + // As stated before, due to the way that the automatic differentiation + // libraries + // work, the @p ad_helper will always returns the derivatives of the energy + // density function with respect to all field variables simultaneously. + // For this reason, it does not make sense to compute the derivatives in + // the functions `get_B()`, `get_S()`, etc. because we'd be doing a lot of + // extra computations that are then simply discarded. So, the best way to + // deal with that is to have a single function call that does all of the + // calculations up-front, and then we extract the stored data as its needed. + // That's what we'll do in the `update_internal_data()` method. As the + // material is rate-independent, we can ignore the DiscreteTime argument. + template + void + Magnetoelastic_Constitutive_Law_AD::update_internal_data( + const SymmetricTensor<2, dim> &C, + const Tensor<1, dim> & H, + const DiscreteTime &) + { + Assert(determinant(C) > 0, ExcInternalError()); + + // Since we reuse the @p ad_helper data structure at each time step, + // we need to clear it of all stale information before use. + ad_helper.reset(); + + // The next step is to set the values for all field components. + // These define the "point" around which we'll be computing the function + // gradients and their linearization. + // The extractors that we created before provide the association between + // the fields and the registry within the @p ad_helper -- they'll be used + // repeatedly to ensure that we have the correct interpretation of which + // variable corresponds to which component of `H` or `C`. + ad_helper.register_independent_variable(H, H_components); + ad_helper.register_independent_variable(C, C_components); + + // Now that we've done the initial setup, we can retrieve the AD + // counterparts of our fields. These are truly the independent variables + // for the energy function, and are "sensitive" to the calculations that + // are performed with them. Notice that the AD number are treated as a + // special number type, and can be used in many templated classes (in this + // example, as the scalar type for the Tensor and SymmetricTensor class). + const Tensor<1, dim, ADNumberType> H_ad = + ad_helper.get_sensitive_variables(H_components); + const SymmetricTensor<2, dim, ADNumberType> C_ad = + ad_helper.get_sensitive_variables(C_components); + + // We can also use them in many functions that are templated on the + // scalar type. So, for these intermediate values that we require, + // we can perform tensor operations and some mathematical functions. + // The resulting type will also be an automatically differentiable + // number, which encodes the operations performed in these functions. + const ADNumberType det_F_ad = std::sqrt(determinant(C_ad)); + const SymmetricTensor<2, dim, ADNumberType> C_inv_ad = invert(C_ad); + AssertThrow(det_F_ad > ADNumberType(0.0), + ExcMessage("Volumetric Jacobian must be positive.")); + + // Next we'll compute the scaling function that will cause the shear + // modulus to change (increase) under the influence of a magnetic field... + const ADNumberType f_mu_e_ad = + 1.0 + (this->get_mu_e_inf() / this->get_mu_e() - 1.0) * + std::tanh((2.0 * H_ad * H_ad) / + (this->get_mu_e_h_sat() * this->get_mu_e_h_sat())); + + // ... and then we can define the material stored energy density function. + // We'll see later that this example is sufficiently complex to warrant + // the use of AD to, at the very least, verify an unassisted + // implementation. + const ADNumberType psi_ad = + 0.5 * this->get_mu_e() * f_mu_e_ad * + (trace(C_ad) - dim - 2.0 * std::log(det_F_ad)) // + + this->get_lambda_e() * std::log(det_F_ad) * std::log(det_F_ad) // + - 0.5 * this->get_mu_0() * this->get_mu_r() * det_F_ad * + (H_ad * C_inv_ad * H_ad); // + + // The stored energy density function is, in fact, the dependent variable + // for this problem, so as a final step in the "configuration" phase, + // we register its definition with the @p ad_helper. + ad_helper.register_dependent_variable(psi_ad); + + // Finally, we can retrieve the resulting value of the stored energy + // density function, as well as its gradient and Hessian with respect + // to the input fields, and cache them. + psi = ad_helper.compute_value(); + ad_helper.compute_gradient(Dpsi); + ad_helper.compute_hessian(D2psi); + } + + // The following few functions then allow for querying the so-stored value + // of $\psi_{0}$, and to extract the desired components of the gradient + // vector and Hessian matrix. We again make use of the extractors to express + // which parts of the total gradient vector and Hessian matrix we wish to + // retrieve. They only return the derivatives of the energy function, so + // for our definitions of the kinetic variables and their linearization a + // few more manipulations are required to form the desired result. + template + double Magnetoelastic_Constitutive_Law_AD::get_psi() const + { + return psi; + } + + + template + Tensor<1, dim> + Magnetoelastic_Constitutive_Law_AD::get_B() const + { + const Tensor<1, dim> dpsi_dH = + ad_helper.extract_gradient_component(Dpsi, H_components); + return -dpsi_dH; + } + + + template + SymmetricTensor<2, dim> + Magnetoelastic_Constitutive_Law_AD::get_S() const + { + const SymmetricTensor<2, dim> dpsi_dC = + ad_helper.extract_gradient_component(Dpsi, C_components); + return 2.0 * dpsi_dC; + } + + + template + SymmetricTensor<2, dim> + Magnetoelastic_Constitutive_Law_AD::get_DD() const + { + const Tensor<2, dim> dpsi_dH_dH = + ad_helper.extract_hessian_component(D2psi, H_components, H_components); + return -symmetrize(dpsi_dH_dH); + } + + // Note that for coupled terms the order of the extractor + // arguments is especially important, as it dictates the order in which + // the directional derivatives are taken. So, if we'd reversed the order + // of the extractors in the call to `extract_hessian_component()` then we'd + // actually have been retrieving part of $\left[ \mathfrak{P}^{\text{tot}} + // \right]^{T}$. + template + Tensor<3, dim> + Magnetoelastic_Constitutive_Law_AD::get_PP() const + { + const Tensor<3, dim> dpsi_dC_dH = + ad_helper.extract_hessian_component(D2psi, C_components, H_components); + return -2.0 * dpsi_dC_dH; + } + + + template + SymmetricTensor<4, dim> + Magnetoelastic_Constitutive_Law_AD::get_HH() const + { + const SymmetricTensor<4, dim> dpsi_dC_dC = + ad_helper.extract_hessian_component(D2psi, C_components, C_components); + return 4.0 * dpsi_dC_dC; + } + + + // @sect4{Magneto-viscoelastic constitutive law (using symbolic algebra and differentiation)} + + // The second material law that we'll consider will be one that represents + // a magneto-viscoelastic material with a single dissipative mechanism. + // We'll consider the free energy density function for such a material to + // be defined as + // @f{align*}{ + // \psi_{0} \left( \mathbf{C}, \mathbf{C}_{v}, \boldsymbol{\mathbb{H}} + // \right) + // &= \psi_{0}^{ME} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) + // + \psi_{0}^{MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) + // \\ \psi_{0}^{ME} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) + // &= \frac{1}{2} \mu_{e} f_{\mu_{e}^{ME}} \left( \boldsymbol{\mathbb{H}} + // \right) + // \left[ \text{tr}(\mathbf{C}) - d - 2 \ln (\text{det}(\mathbf{F})) + // \right] + // + \lambda_{e} \ln^{2} \left(\text{det}(\mathbf{F}) \right) + // - \frac{1}{2} \mu_{0} \mu_{r} \text{det}(\mathbf{F}) + // \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}} \right] + // \\ \psi_{0}^{MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) + // &= \frac{1}{2} \mu_{v} f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} + // \right) + // \left[ \mathbf{C}_{v} : \left[ + // \left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + // \mathbf{C} \right] - d - \ln\left( + // \text{det}\left(\mathbf{C}_{v}\right) \right) \right] + // @f} + // with + // @f[ + // f_{\mu_{e}}^{ME} \left( \boldsymbol{\mathbb{H}} \right) + // = 1 + \left[ \frac{\mu_{e}^{\infty}}{\mu_{e}} - 1 \right] + // \tanh \left( 2 \frac{\boldsymbol{\mathbb{H}} \cdot + // \boldsymbol{\mathbb{H}}} + // {\left(h_{e}^{\text{sat}}\right)^{2}} \right) + // @f] + // @f[ + // f_{\mu_{v}}^{MVE} \left( \boldsymbol{\mathbb{H}} \right) + // = 1 + \left[ \frac{\mu_{v}^{\infty}}{\mu_{v}} - 1 \right] + // \tanh \left( 2 \frac{\boldsymbol{\mathbb{H}} \cdot + // \boldsymbol{\mathbb{H}}} + // {\left(h_{v}^{\text{sat}}\right)^{2}} \right), + // @f] + // in conjunction with the evolution law for the internal viscous variable + // @f[ + // \mathbf{C}_{v}^{(t)} + // = \frac{1}{1 + \frac{\Delta t}{\tau_{v}}} \left[ + // \mathbf{C}_{v}^{(t-1)} + // + \frac{\Delta t}{\tau_{v}} + // \left[\left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + // \mathbf{C} \right]^{-1} + // \right] + // @f] + // that was discretized using a first-order backward difference + // approximation. + // + // Again, let us see how this is implemented in a concrete class. Instead of + // the AD framework used in the previous class, we will now employ the SD + // approach. To support this, + // the class constructor accepts not only the @p constitutive_parameters, + // but also two additional variables that will be used to initialize + // a Differentiation::SD::BatchOptimizer. We'll give more context to this + // later. + template + class Magnetoviscoelastic_Constitutive_Law_SD final + : public Coupled_Magnetomechanical_Constitutive_Law_Base + { + public: + Magnetoviscoelastic_Constitutive_Law_SD( + const ConstitutiveParameters & constitutive_parameters, + const Differentiation::SD::OptimizerType optimizer_type, + const Differentiation::SD::OptimizationFlags optimization_flags); + + // Like for the automatic differentiation helper, the + // Differentiation::SD::BatchOptimizer will return a collection of + // results all at once. So, in order to do that just once, we'll employ + // a similar approach to before and do all of the expensive calculations + // within the `update_internal_data()` function, and cache the results + // for layer extraction. + virtual void update_internal_data(const SymmetricTensor<2, dim> &C, + const Tensor<1, dim> & H, + const DiscreteTime &time) override; + + virtual double get_psi() const override; + + virtual Tensor<1, dim> get_B() const override; + + virtual SymmetricTensor<2, dim> get_S() const override; + + virtual SymmetricTensor<2, dim> get_DD() const override; + + virtual Tensor<3, dim> get_PP() const override; + + virtual SymmetricTensor<4, dim> get_HH() const override; + + // Since we're dealing with a rate dependent material, we'll have to + // update the history variable at the appropriate time. That will be the + // purpose of this function. + virtual void update_end_of_timestep() override; + + // In the `private` part of the class, we will want to + // keep track of the internal viscous deformation, so the following + // two (real-valued, non-symbolic) member variables respectively hold + // - the value of internal variable time step (and, if embedded within a + // nonlinear solver framework, Newton step), and + // - the value of internal variable at the previous timestep. + // + // (We've labeled these variables "Q" so that they're easy to identify; + // in a sea of calculations it is not necessarily easy to distinguish + // `Cv` or `C_v` from `C`.) + private: + SymmetricTensor<2, dim> Q_t; + SymmetricTensor<2, dim> Q_t1; + + // As we'll be using symbolic types, we'll need to define some symbolic + // variables to use with the framework. (They are all suffixed with "SD" + // to make it easy to distinguish the symbolic types or expressions from + // real-valued types or scalars.) This can be done once up front + // (potentially even as `static` variables) to minimize the overhead + // associated with creating these variables. For the ultimate in generic + // programming, we can even describe the constitutive parameters + // symbolically, *potentially* allowing a single class instance to be + // reused with different inputs for these values too. + // + // These are the symbolic scalars that represent the elastic, viscous, and + // magnetic material parameters + // (defined mostly in the same order as they appear in the @p ConstitutiveParameters + // class). We also store a symbolic expression, @p delta_t_sd, that represents the + // time step size): + const Differentiation::SD::Expression mu_e_sd; + const Differentiation::SD::Expression mu_e_inf_sd; + const Differentiation::SD::Expression mu_e_h_sat_sd; + const Differentiation::SD::Expression lambda_e_sd; + const Differentiation::SD::Expression mu_v_sd; + const Differentiation::SD::Expression mu_v_inf_sd; + const Differentiation::SD::Expression mu_v_h_sat_sd; + const Differentiation::SD::Expression tau_v_sd; + const Differentiation::SD::Expression delta_t_sd; + const Differentiation::SD::Expression mu_r_sd; + + // Next we define some tensorial symbolic variables that represent the + // independent field variables, upon which the energy density function + // is parameterized: + const Tensor<1, dim, Differentiation::SD::Expression> H_sd; + const SymmetricTensor<2, dim, Differentiation::SD::Expression> C_sd; + + // And similarly we have the symbolic representation of the internal + // viscous variables (both its current value and its value at the + // previous timestep): + const SymmetricTensor<2, dim, Differentiation::SD::Expression> Q_t_sd; + const SymmetricTensor<2, dim, Differentiation::SD::Expression> Q_t1_sd; + + // We should also store the definitions of the dependent expressions: + // Although we'll only compute them once, we require them to retrieve + // data from the @p optimizer that is declared below. + // Furthermore, when serializing a material class like this one (not done + // as a part of this tutorial) we'd either need to serialize these + // expressions as well or we'd need to reconstruct them upon reloading. + Differentiation::SD::Expression psi_sd; + Tensor<1, dim, Differentiation::SD::Expression> B_sd; + SymmetricTensor<2, dim, Differentiation::SD::Expression> S_sd; + SymmetricTensor<2, dim, Differentiation::SD::Expression> BB_sd; + Tensor<3, dim, Differentiation::SD::Expression> PP_sd; + SymmetricTensor<4, dim, Differentiation::SD::Expression> HH_sd; + + // The next variable is then the optimizer that is used to evaluate the + // dependent functions. More specifically, it provides the possibility to + // accelerate the evaluation of the symbolic dependent expressions. This + // is a vital tool, because the native evaluation of lengthy expressions + // (using no method of acceleration, but rather direct evaluation directly + // of the symbolic expressions) can be very slow. The + // Differentiation::SD::BatchOptimizer class provides a mechanism by which + // to transform the symbolic expression tree into another code path that, + // for example, shares intermediate results between the various dependent + // expressions (meaning that these intermediate values only get calculated + // once per evaluation) and/or compiling the code using a just-in-time + // compiler (thereby retrieving near-native performance for the evaluation + // step). + // + // Performing this code transformation is very computationally expensive, + // so we store the optimizer so that it is done just once per class + // instance. This also further motivates the decision to make the + // constitutive parameters themselves symbolic. We could then reuse a + // single instance + // of this @p optimizer across several materials (with the same energy + // function, of course) and potentially multiple continuum points (if + // embedded within a finite element simulation). + // + // As specified by the template parameter, the numerical result will be of + // type double. + Differentiation::SD::BatchOptimizer optimizer; + + // During the evaluation phase, we must map the symbolic variables to + // their real-valued counterparts. The next method will provide this + // functionality. + // + // The final method of this class will configure the @p optimizer. + Differentiation::SD::types::substitution_map + make_substitution_map(const SymmetricTensor<2, dim> &C, + const Tensor<1, dim> & H, + const double delta_t) const; + + void initialize_optimizer(); + }; + + // As the resting deformation state is one at which the material is + // considered to be completely relaxed, the internal viscous variables are + // initialized with the identity tensor, i.e. $\mathbf{C}_{v} = \mathbf{I}$. + // The various symbolic variables representing the constitutive parameters, + // time step size, and field and internal variables all get a unique + // identifier. The optimizer is passed the two parameters that declare which + // optimization (acceleration) technique should be employed, as well as + // which additional steps should be taken by the CAS to help improve + // performance during evaluation. + template + Magnetoviscoelastic_Constitutive_Law_SD:: + Magnetoviscoelastic_Constitutive_Law_SD( + const ConstitutiveParameters & constitutive_parameters, + const Differentiation::SD::OptimizerType optimizer_type, + const Differentiation::SD::OptimizationFlags optimization_flags) + : Coupled_Magnetomechanical_Constitutive_Law_Base( + constitutive_parameters) + , Q_t(Physics::Elasticity::StandardTensors::I) + , Q_t1(Physics::Elasticity::StandardTensors::I) + , mu_e_sd("mu_e") + , mu_e_inf_sd("mu_e_inf") + , mu_e_h_sat_sd("mu_e_h_sat") + , lambda_e_sd("lambda_e") + , mu_v_sd("mu_v") + , mu_v_inf_sd("mu_v_inf") + , mu_v_h_sat_sd("mu_v_h_sat") + , tau_v_sd("tau_v") + , delta_t_sd("delta_t") + , mu_r_sd("mu_r") + , H_sd(Differentiation::SD::make_vector_of_symbols("H")) + , C_sd(Differentiation::SD::make_symmetric_tensor_of_symbols<2, dim>("C")) + , Q_t_sd( + Differentiation::SD::make_symmetric_tensor_of_symbols<2, dim>("Q_t")) + , Q_t1_sd( + Differentiation::SD::make_symmetric_tensor_of_symbols<2, dim>("Q_t1")) + , optimizer(optimizer_type, optimization_flags) + { + initialize_optimizer(); + } + + // The substitution map simply pairs all of the following data together: + // - the constitutive parameters (with values retrieved from the base + // class), + // - the time step size (with its value retrieved from the time + // discretizer), + // - the field values (with their values being prescribed by an external + // function that is calling into this @p Magnetoviscoelastic_Constitutive_Law_SD instance), and + // - the current and previous internal viscous deformation (with their + // values + // stored within this class instance). + template + Differentiation::SD::types::substitution_map + Magnetoviscoelastic_Constitutive_Law_SD::make_substitution_map( + const SymmetricTensor<2, dim> &C, + const Tensor<1, dim> & H, + const double delta_t) const + { + return Differentiation::SD::make_substitution_map( + std::make_pair(mu_e_sd, this->get_mu_e()), + std::make_pair(mu_e_inf_sd, this->get_mu_e_inf()), + std::make_pair(mu_e_h_sat_sd, this->get_mu_e_h_sat()), + std::make_pair(lambda_e_sd, this->get_lambda_e()), + std::make_pair(mu_v_sd, this->get_mu_v()), + std::make_pair(mu_v_inf_sd, this->get_mu_v_inf()), + std::make_pair(mu_v_h_sat_sd, this->get_mu_v_h_sat()), + std::make_pair(tau_v_sd, this->get_tau_v()), + std::make_pair(delta_t_sd, delta_t), + std::make_pair(mu_r_sd, this->get_mu_r()), + std::make_pair(H_sd, H), + std::make_pair(C_sd, C), + std::make_pair(Q_t_sd, Q_t), + std::make_pair(Q_t1_sd, Q_t1)); + } + + // Due to the "natural" use of the symbolic expressions, much of the + // procedure to configure the @p optimizer looks very similar to that which + // is used to construct the automatic differentiation helper. + // Nevertheless, we'll detail these steps again to highlight the differences + // that underlie the two frameworks. + // + // The function starts with expressions that symbolically encode the + // determinant of the deformation gradient (as expressed in terms of the + // right Cauchy-Green deformation tensor, our primary field variable), as + // well as the inverse of $\mathbf{C}$ itself: + template + void Magnetoviscoelastic_Constitutive_Law_SD::initialize_optimizer() + { + const Differentiation::SD::Expression det_F_sd = + std::sqrt(determinant(C_sd)); + const SymmetricTensor<2, dim, Differentiation::SD::Expression> C_inv_sd = + invert(C_sd); + + // Next is the symbolic representation of the saturation function for + // the elastic part of the free energy density function, followed by the + // magnetoelastic contribution to the free energy density function. + // This all has the same stucture as we'd seen previously. + const Differentiation::SD::Expression f_mu_e_sd = + 1.0 + + (mu_e_inf_sd / mu_e_sd - 1.0) * + std::tanh((2.0 * H_sd * H_sd) / (mu_e_h_sat_sd * mu_e_h_sat_sd)); + + const Differentiation::SD::Expression psi_ME_sd = + 0.5 * mu_e_sd * f_mu_e_sd * + (trace(C_sd) - dim - 2.0 * std::log(det_F_sd)) + + lambda_e_sd * std::log(det_F_sd) * std::log(det_F_sd) - + 0.5 * this->get_mu_0() * mu_r_sd * det_F_sd * (H_sd * C_inv_sd * H_sd); + + // In addition, we define the magneto-viscoelastic contribution to the + // free energy density function. The first component required to implement + // this is a scaling function that will cause the viscous shear modulus to + // change (increase) under the influence of a magnetic field (see + // @cite Pelteret2018a, equation 29). Thereafter we can compute the + // dissipative component of the energy density function; its expression + // is stated in @cite Pelteret2018a (equation 28), which is a + // straight-forward extension of an energy density function formulated by + // @cite Linder2011a (equation 46). + const Differentiation::SD::Expression f_mu_v_sd = + 1.0 + + (mu_v_inf_sd / mu_v_sd - 1.0) * + std::tanh((2.0 * H_sd * H_sd) / (mu_v_h_sat_sd * mu_v_h_sat_sd)); + + const Differentiation::SD::Expression psi_MVE_sd = + 0.5 * mu_v_sd * f_mu_v_sd * + (Q_t_sd * (std::pow(det_F_sd, -2.0 / dim) * C_sd) - dim - + std::log(determinant(Q_t_sd))); + + // From these building blocks, we can then define the material's total + // free energy density function: + psi_sd = psi_ME_sd + psi_MVE_sd; + + // As it stands, to the CAS the variable @p Q_t_sd appears + // to be independent of @p C_sd. Our tensorial symbolic expression + // @p Q_t_sd just has an identifier associated with it, and there is + // nothing that links it to the other tensorial symbolic expression + // @p C_sd. So any derivatives taken with respect to @p C_sd will ignore + // this inherent dependence which, as we can see from the evolution law, + // is in fact + // $\mathbf{C}_{v} = \mathbf{C}_{v} \left( \mathbf{C}, t \right)$. + // This means that deriving any function $f = f(\mathbf{C}, \mathbf{Q})$ + // with respect to $\mathbf{C}$ will return partial derivatives + // $\frac{\partial f(\mathbf{C}, \mathbf{Q})}{\partial \mathbf{C}} + // \Big\vert_{\mathbf{C}_{v}}$ as opposed to the total derivative + // $\frac{d f(\mathbf{C}, \mathbf{Q}(\mathbf{C}))}{d \mathbf{C}} = + // \frac{\partial f(\mathbf{C}, \mathbf{Q}(\mathbf{C}))}{\partial + // \mathbf{C}} \Big\vert_{\mathbf{C}_{v}} + \frac{\partial f(\mathbf{C}, + // \mathbf{Q}(\mathbf{C}))}{\partial \mathbf{C}_{v}} + // \Big\vert_{\mathbf{C}} : \frac{d \mathbf{Q}(\mathbf{C}))}{d + // \mathbf{C}}$. + // + // By contrast, with the current AD libraries the total derivative would + // always be returned. This implies that the computed kinetic variables + // would be incorrect for this class of material model, making AD the + // incorrect tool from which to derive (at the continuum point level) the + // constitutive law for this dissipative material from an energy density + // function. + // + // It is this specific level of control that characterizes a defining + // difference difference between the SD and AD frameworks. In a few lines + // we'll be manipulating the expression for the internal variable + // @p Q_t_sd such that it produces the correct linearization. + + // But, first, we'll compute the symbolic expressions for the kinetic + // variables, i.e., the magnetic induction vector and the Piola-Kirchhoff + // stress tensor. The code that performs the differentiation quite closely + // mimics the definition stated in the theory. + B_sd = -Differentiation::SD::differentiate(psi_sd, H_sd); + S_sd = 2.0 * Differentiation::SD::differentiate(psi_sd, C_sd); + + // Since the next step is to linearize the above, it is the appropriate + // time to inform the CAS of the explicit dependency of @p Q_t_sd on @p C_sd, + // i.e., state that $\mathbf{C}_{v} = \mathbf{C}_{v} \left( \mathbf{C}, t + // \right)$. This means that all future differential operations made with + // respect + // to @p C_sd will take into account this dependence (i.e., compute total derivatives). + // In other words, we will transform some expression such that their + // intrinsic parameterization changes from $f(\mathbf{C}, \mathbf{Q})$ + // to $f(\mathbf{C}, \mathbf{Q}(\mathbf{C}))$. + // + // To do this, we consider the time-discrete evolution law. + // From that, we have the explicit expression for the internal + // variable in terms of its history as well as the primary + // field variable. That is what it described in this expression: + const SymmetricTensor<2, dim, Differentiation::SD::Expression> + Q_t_sd_explicit = + (1.0 / (1.0 + delta_t_sd / tau_v_sd)) * + (Q_t1_sd + + (delta_t_sd / tau_v_sd * std::pow(det_F_sd, 2.0 / dim) * C_inv_sd)); + + // Next we produce an intermediate substitution map, which will take + // every instance of @p Q_t_sd (our identifier) found in an expression + // and replace it with the full expression held in @p Q_t_sd_explicit. + const Differentiation::SD::types::substitution_map + substitution_map_explicit = Differentiation::SD::make_substitution_map( + std::make_pair(Q_t_sd, Q_t_sd_explicit)); + + // We can the perform this substitution on the two kinetic variables + // and immediately differentiate the result that appears after that + // substitution with the field variables. (If you'd like, this could + // be split up into two steps with the intermediate results stored in + // a temporary variable.) Again, if you overlook the "complexity" + // generated by the substitution, these calls that linearize the kinetic + // variables and produce the three tangent tensors quite closely resembles + // what's stated in the theory. + BB_sd = symmetrize(Differentiation::SD::differentiate( + Differentiation::SD::substitute(B_sd, substitution_map_explicit), + H_sd)); + PP_sd = -Differentiation::SD::differentiate( + Differentiation::SD::substitute(S_sd, substitution_map_explicit), H_sd); + HH_sd = + 2.0 * + Differentiation::SD::differentiate( + Differentiation::SD::substitute(S_sd, substitution_map_explicit), + C_sd); + + // Now we need to tell the @p optimizer what entries we need to provide + // numerical values for in order for it to successfully perform its + // calculations. These essentially act as the input arguments to + // all dependent functions that the @p optimizer must evaluate. + // They are, collectively, the independent variables + // for the problem, the history variables, the time step sie and the + // constitutive parameters (since we've not hard encoded them in the + // energy density function). + // + // So what we really want is to provide it a collection of + // symbols, which one could accomplish in this way: + // @code + // optimizer.register_symbols(Differentiation::SD::make_symbol_map( + // mu_e_sd, mu_e_inf_sd, mu_e_h_sat_sd, lambda_e_sd, + // mu_v_sd, mu_v_inf_sd, mu_v_h_sat_sd, tau_v_sd, + // delta_t_sd, mu_r_sd, + // H_sd, C_sd, + // Q_t_sd, Q_t1_sd)); + // @endcode + // But this is all actually already encoded as the keys of the + // substitution map. Doing the above would also mean that we + // need to manage the symbols in two places (here and when constructing + // the substitution map), which is annoying and a potential source of + // error if this material class is modified or extended. + // Since we're not interested in the values at this point, + // it is alright if the substitution map is filled with invalid data + // for the values associated with each key entry. + // So we'll simply create a fake substitution map, and extract the + // symbols from that. Note that any substitution map passed to the + // @p optimizer will have to, at the very least, contain entries for + // these symbols. + optimizer.register_symbols( + Differentiation::SD::Utilities::extract_symbols( + make_substitution_map({}, {}, 0))); + + // We then inform the optimizer of what values we want calculated, which + // in our situation encompasses all of the dependent variables (namely + // the energy density function and its various derivatives). + optimizer.register_functions(psi_sd, B_sd, S_sd, BB_sd, PP_sd, HH_sd); + + // The last step is to finalize the optimizer. With this call it will + // determine an equivalent code path that will evaluate all of the + // dependent functions at once, but with less computational + // cost than when evaluating the symbolic expression directly. + // Note: This is an expensive call, so we want execute it as few times + // as possible. We've done it in the constructor of our class, which + // achieves the goal of being called only once per class instance. + optimizer.optimize(); + } + + // Since the configuration of the @p optimizer was done up front, there's + // very little to do each time we want to compute kinetic variables or + // their linearization (derivatives). + template + void Magnetoviscoelastic_Constitutive_Law_SD::update_internal_data( + const SymmetricTensor<2, dim> &C, + const Tensor<1, dim> & H, + const DiscreteTime & time) + { + // To update the internal history variable, we first need to compute + // a few fundamental quantities, which we've seen before. + // We can also ask the time discretizer for the time step size that + // was used to iterate from the previous time step to the current one. + const double delta_t = this->get_delta_t(time); + + const double det_F = std::sqrt(determinant(C)); + const SymmetricTensor<2, dim> C_inv = invert(C); + AssertThrow(det_F > 0.0, + ExcMessage("Volumetric Jacobian must be positive.")); + + // Now we can update the (real valued) internal viscous deformation + // tensor, as per the definition given by the evolution law in conjunction + // with the chosen time discretization scheme. + Q_t = (1.0 / (1.0 + delta_t / this->get_tau_v())) * + (Q_t1 + (delta_t / this->get_tau_v()) * std::pow(det_F, 2.0 / dim) * + C_inv); + + // Next we pass the optimizer the numeric values that we wish the + // independent variables, time step size and (implicit to this call), + // the constitutive parameters to represent. + const auto substitution_map = make_substitution_map(C, H, delta_t); + + // When making this next call, the call path used to (numerically) + // evaluate the dependent functions is quicker than dictionary + // substitution. + optimizer.substitute(substitution_map); + } + + // Having called `update_internal_data()`, it is then valid to + // extract data from the optimizer. + // When doing the evaluation, we need the exact symbolic expressions of + // the data to extracted from the optimizer. The implication of this + // is that we needed to store the symbolic expressions of all dependent + // variables for the lifetime of the optimizer (naturally, the same + // is implied for the input variables). + template + double Magnetoviscoelastic_Constitutive_Law_SD::get_psi() const + { + return optimizer.evaluate(psi_sd); + } + + + template + Tensor<1, dim> Magnetoviscoelastic_Constitutive_Law_SD::get_B() const + { + return optimizer.evaluate(B_sd); + } + + + template + SymmetricTensor<2, dim> + Magnetoviscoelastic_Constitutive_Law_SD::get_S() const + { + return optimizer.evaluate(S_sd); + } + + + template + SymmetricTensor<2, dim> + Magnetoviscoelastic_Constitutive_Law_SD::get_DD() const + { + return optimizer.evaluate(BB_sd); + } + + + template + Tensor<3, dim> Magnetoviscoelastic_Constitutive_Law_SD::get_PP() const + { + return optimizer.evaluate(PP_sd); + } + + + template + SymmetricTensor<4, dim> + Magnetoviscoelastic_Constitutive_Law_SD::get_HH() const + { + return optimizer.evaluate(HH_sd); + } + + // When moving forward in time, the "current" state of the internal variable + // instantaneously defines the state at the "previous" timestep. As such, we + // record value of history variable for use as the "past value" at the next + // time step. + template + void Magnetoviscoelastic_Constitutive_Law_SD::update_end_of_timestep() + { + Q_t1 = Q_t; + }; + + + // @sect3{A more complex example (continued): Parameters and hand-derived material classes} + + // Now that we've seen how the AD and SD frameworks can make light(er) work + // of defining these constitutive laws, we'll implement the equivalent + // classes by hand for the purpose of verification and to do some + // preliminary benchmarking of the frameworks versus a native + // implementation. + // + // At the expense of the author's sanity, what is documented below + // (hopefully accurately) are the full definitions for the kinetic variables + // and their tangents, as well as some intermediate computations. Since the + // structure and design of the constitutive law classes has been outlined + // earlier, we'll gloss over it and simply delineate between the various + // stages of calculations in the `update_internal_data()` method definition. + // It should be easy enough to link the derivative calculations (with their + // moderately expressive variable names) to their documented definitions + // that appear in the class descriptions. + // We will, however, take the opportunity to present two different paradigms + // for implementing constitutive law classes. The second will provide more + // flexibility than the first (thereby making it more easily extensible, + // in the author's opinion) at the expense of some performance. + + // @sect4{Magnetoelastic constitutive law (hand-derived)} + + // From the stored energy that, as mentioned earlier, is defined as + // @f[ + // \psi_{0} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) + // = \frac{1}{2} \mu_{e} f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} \right) + // \left[ \text{tr}(\mathbf{C}) - d - 2 \ln (\text{det}(\mathbf{F})) + // \right] + // + \lambda_{e} \ln^{2} \left(\text{det}(\mathbf{F}) \right) + // - \frac{1}{2} \mu_{0} \mu_{r} \text{det}(\mathbf{F}) + // \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}} \right] + // @f] + // with + // @f[ + // f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} \right) + // = 1 + \left[ \frac{\mu_{e}^{\infty}}{\mu_{e}} - 1 \right] + // \tanh \left( 2 \frac{\boldsymbol{\mathbb{H}} \cdot + // \boldsymbol{\mathbb{H}}} + // {\left(h_{e}^{\text{sat}}\right)^{2}} \right) , + // \\ \text{det}(\mathbf{F}) = \sqrt{\text{det}(\mathbf{C})} + // @f] + // for this magnetoelastic material, the first derivatives that correspond + // to the magnetic induction vector and total Piola-Kirchhoff stress + // tensor are + // @f[ + // \boldsymbol{\mathbb{B}} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} + // \right) + // \dealcoloneq - \frac{d \psi_{0}}{d \boldsymbol{\mathbb{H}}} + // = - \frac{1}{2} \mu_{e} \left[ \text{tr}(\mathbf{C}) - d - 2 \ln + // (\text{det}(\mathbf{F})) + // \right] \frac{d f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} + // \right)}{d \boldsymbol{\mathbb{H}}} + // + \mu_{0} \mu_{r} \text{det}(\mathbf{F}) \left[ \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}} + // \right] + // @f] + // @f{align} + // \mathbf{S}^{\text{tot}} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} + // \right) + // \dealcoloneq 2 \frac{d \psi_{0} \left( \mathbf{C}, + // \boldsymbol{\mathbb{H}} \right)}{d \mathbf{C}} + // &= \mu_{e} f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} \right) + // \left[ \frac{d\,\text{tr}(\mathbf{C})}{d \mathbf{C}} + // - 2 \frac{1}{\text{det}(\mathbf{F})} + // \frac{d\,\text{det}(\mathbf{F})}{d \mathbf{C}} \right] + // + 4 \lambda_{e} \ln \left(\text{det}(\mathbf{F}) \right) + // \frac{1}{\text{det}(\mathbf{F})} \frac{d\,\text{det}(\mathbf{F})}{d + // \mathbf{C}} + // - \mu_{0} \mu_{r} \left[ + // \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}} \right] \frac{d\,\text{det}(\mathbf{F})}{d + // \mathbf{C}} + \text{det}(\mathbf{F}) \frac{d \left[ + // \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}} + // \right]}{d \mathbf{C}} \right] + // \\ &= \mu_{e} f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} \right) + // \left[ \mathbf{I} - \mathbf{C}^{-1} \right] + // + 2 \lambda_{e} \ln \left(\text{det}(\mathbf{F}) \right) \mathbf{C}^{-1} + // - \mu_{0} \mu_{r} \left[ + // \frac{1}{2} \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} + // \cdot \boldsymbol{\mathbb{H}} \right] \text{det}(\mathbf{F}) + // \mathbf{C}^{-1} + // - \text{det}(\mathbf{F}) + // \left[ \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}} \right] \otimes + // \left[ \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}} \right] + // \right] + // @f} + // with + // @f[ + // \frac{d f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} \right)}{d + // \boldsymbol{\mathbb{H}}} + // = \left[ \frac{\mu_{e}^{\infty}}{\mu_{e}} - 1 \right] + // \text{sech}^{2} \left( 2 \frac{\boldsymbol{\mathbb{H}} \cdot + // \boldsymbol{\mathbb{H}}} + // {\left(h_{e}^{\text{sat}}\right)^{2}} \right) + // \left[ \frac{4} {\left(h_{e}^{\text{sat}}\right)^{2}} + // \boldsymbol{\mathbb{H}} \right] + // @f] + // @f[ + // \frac{d\,\text{tr}(\mathbf{C})}{d \mathbf{C}} + // = \mathbf{I} + // \quad \text{(the second-order identity tensor)} + // @f] + // @f[ + // \frac{d\,\text{det}(\mathbf{F})}{d \mathbf{C}} + // = \frac{1}{2} \text{det}(\mathbf{F}) \mathbf{C}^{-1} + // @f] + // @f[ + // \frac{d C^{-1}_{ab}}{d C_{cd}} + // = - \text{sym} \left( C^{-1}_{ac} C^{-1}_{bd} \right) + // = -\frac{1}{2} \left[ C^{-1}_{ac} C^{-1}_{bd} + C^{-1}_{ad} C^{-1}_{bc} + // \right] + // @f] + // @f[ + // \frac{d \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}} \right]}{d \mathbf{C}} + // = - \left[ \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}} \right] \otimes + // \left[ \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}} \right] + // @f] + // The use of the symmetry operator $\text{sym} \left( \bullet \right)$ in + // the one derivation above helps to ensure that the resulting rank-4 + // tensor, which holds minor symmetries due to the symmetry of $\mathbf{C}$, + // still maps rank-2 symmetric tensors to rank-2 symmetric tensors. See the + // SymmetricTensor class documentation and the introduction to step-44 and + // for further explanation as to what symmetry means in the context of + // fourth-order tensors. + // + // The linearization of each of the kinematic variables with respect to + // their arguments are + // @f[ + // \mathbb{D} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) + // = \frac{d \boldsymbol{\mathbb{B}}}{d \boldsymbol{\mathbb{H}}} + // = - \frac{1}{2} \mu_{e} \left[ \text{tr}(\mathbf{C}) - d - 2 \ln + // (\text{det}(\mathbf{F})) + // \right] \frac{d^{2} f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} + // \right)}{d \boldsymbol{\mathbb{H}} \otimes d \boldsymbol{\mathbb{H}}} + // + \mu_{0} \mu_{r} \text{det}(\mathbf{F}) \mathbf{C}^{-1} + // @f] + // @f{align} + // \mathfrak{P}^{\text{tot}} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} + // \right) = - \frac{d \mathbf{S}^{\text{tot}}}{d \boldsymbol{\mathbb{H}}} + // &= - \mu_{e} + // \left[ \frac{d\,\text{tr}(\mathbf{C})}{d \mathbf{C}} + // - 2 \frac{1}{\text{det}(\mathbf{F})} + // \frac{d\,\text{det}(\mathbf{F})}{d \mathbf{C}} \right] + // \otimes \frac{d f_{\mu_{e} \left( \boldsymbol{\mathbb{H}} + // \right)}}{d \boldsymbol{\mathbb{H}}} + // + \mu_{0} \mu_{r} \left[ + // \frac{d\,\text{det}(\mathbf{F})}{d \mathbf{C}} \otimes + // \frac{d \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}} + // \right]}{d \boldsymbol{\mathbb{H}}} \right] + // + \text{det}(\mathbf{F}) + // \frac{d^{2} \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} + // \cdot \boldsymbol{\mathbb{H}} + // \right]}{d \mathbf{C} \otimes d \boldsymbol{\mathbb{H}}} + // \\ &= - \mu_{e} + // \left[ \mathbf{I} - \mathbf{C}^{-1} \right] \otimes + // \frac{d f_{\mu_{e} \left( \boldsymbol{\mathbb{H}} \right)}}{d + // \boldsymbol{\mathbb{H}}} + // + \mu_{0} \mu_{r} \left[ + // \text{det}(\mathbf{F}) \mathbf{C}^{-1} \otimes + // \left[ \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}} \right] + // \right] + // + \text{det}(\mathbf{F}) + // \frac{d^{2} \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} + // \cdot \boldsymbol{\mathbb{H}} + // \right]}{d \mathbf{C} \otimes \mathbf{C} \boldsymbol{\mathbb{H}}} + // @f} + // @f{align} + // \mathcal{H}^{\text{tot}} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} + // \right) = 2 \frac{d \mathbf{S}^{\text{tot}}}{d \mathbf{C}} + // &= 2 \mu_{e} f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} \right) + // \left[ - \frac{d \mathbf{C}^{-1}}{d \mathbf{C}} \right] + // + 4 \lambda_{e} \left[ \mathbf{C}^{-1} \otimes \left[ + // \frac{1}{\text{det}(\mathbf{F})} \frac{d \, \text{det}(\mathbf{F})}{d + // \mathbf{C}} \right] + \ln \left(\text{det}(\mathbf{F}) \right) \frac{d + // \mathbf{C}^{-1}}{d \mathbf{C}} \right] + // \\ &- \mu_{0} \mu_{r} \left[ + // \text{det}(\mathbf{F}) \mathbf{C}^{-1} \otimes \frac{d \left[ + // \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}} \right]}{d \mathbf{C}} + // + \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}} \right] \mathbf{C}^{-1} \otimes \frac{d \, + // \text{det}(\mathbf{F})}{d \mathbf{C}} + // + \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}} \right] \text{det}(\mathbf{F}) \frac{d + // \mathbf{C}^{-1}}{d \mathbf{C}} + // \right] + // \\ &+ 2 \mu_{0} \mu_{r} \left[ \left[ + // \left[ \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}} \right] \otimes + // \left[ \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}} \right] + // \right] \otimes \frac{d \, \text{det}(\mathbf{F})}{d \mathbf{C}} + // - \text{det}(\mathbf{F}) + // \frac{d^{2} \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} + // \cdot \boldsymbol{\mathbb{H}}\right]}{d \mathbf{C} \otimes d + // \mathbf{C}} + // \right] + // \\ &= 2 \mu_{e} f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} \right) + // \left[ - \frac{d \mathbf{C}^{-1}}{d \mathbf{C}} \right] + // + 4 \lambda_{e} \left[ \frac{1}{2} \mathbf{C}^{-1} \otimes + // \mathbf{C}^{-1} + \ln \left(\text{det}(\mathbf{F}) \right) \frac{d + // \mathbf{C}^{-1}}{d \mathbf{C}} \right] + // \\ &- \mu_{0} \mu_{r} \left[ + // - \text{det}(\mathbf{F}) \mathbf{C}^{-1} \otimes \left[ \left[ + // \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}} \right] \otimes + // \left[ \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}} \right] \right] + // + \frac{1}{2} \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}} \right] \text{det}(\mathbf{F}) \mathbf{C}^{-1} + // \otimes \mathbf{C}^{-1} + // + \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}} \right] \text{det}(\mathbf{F}) \frac{d + // \mathbf{C}^{-1}}{d \mathbf{C}} + // \right] + // \\ &+ 2 \mu_{0} \mu_{r} \left[ \frac{1}{2} \text{det}(\mathbf{F}) \left[ + // \left[ \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}} \right] \otimes + // \left[ \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}} \right] + // \right] \otimes \mathbf{C}^{-1} + // - \text{det}(\mathbf{F}) + // \frac{d^{2} \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} + // \cdot \boldsymbol{\mathbb{H}}\right]}{d \mathbf{C} \otimes d + // \mathbf{C}} + // \right] + // @f} + // with + // @f[ + // \frac{d^{2} f_{\mu_{e}} \left( \boldsymbol{\mathbb{H}} \right)}{d + // \boldsymbol{\mathbb{H}} \otimes d \boldsymbol{\mathbb{H}}} + // = -2 \left[ \frac{\mu_{e}^{\infty}}{\mu_{e}} - 1 \right] + // \tanh \left( 2 \frac{\boldsymbol{\mathbb{H}} \cdot + // \boldsymbol{\mathbb{H}}} + // {\left(h_{e}^{\text{sat}}\right)^{2}} \right) + // \text{sech}^{2} \left( 2 \frac{\boldsymbol{\mathbb{H}} \cdot + // \boldsymbol{\mathbb{H}}} + // {\left(h_{e}^{\text{sat}}\right)^{2}} \right) + // \left[ \frac{4} {\left(h_{e}^{\text{sat}}\right)^{2}} \mathbf{I} + // \right] + // @f] + // @f[ + // \frac{d \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}} + // \right]}{d \boldsymbol{\mathbb{H}}} + // = 2 \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}} + // @f] + // @f[ + // \frac{d^{2} \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}}\right]}{d \mathbf{C} \otimes d + // \boldsymbol{\mathbb{H}}} \Rightarrow \frac{d^{2} \left[ \mathbb{H}_{e} + // C^{-1}_{ef} \mathbb{H}_{f} + // \right]}{d C_{ab} d \mathbb{H}_{c}} + // = - C^{-1}_{ac} C^{-1}_{be} \mathbb{H}_{e} - C^{-1}_{ae} \mathbb{H}_{e} + // C^{-1}_{bc} + // @f] + // @f{align} + // \frac{d^{2} \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}}\right]}{d \mathbf{C} \otimes d \mathbf{C}} + // &= -\frac{d \left[\left[ \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}} + // \right] \otimes + // \left[ \mathbf{C}^{-1} \cdot \boldsymbol{\mathbb{H}} + // \right]\right]}{d \mathbf{C}} + // \\ \Rightarrow + // \frac{d^{2} \left[ \mathbb{H}_{e} C^{-1}_{ef} \mathbb{H}_{f} + // \right]}{d C_{ab} d C_{cd}} + // &= \text{sym} \left( C^{-1}_{ae} \mathbb{H}_{e} C^{-1}_{cf} + // \mathbb{H}_{f} C^{-1}_{bd} + // + C^{-1}_{ce} \mathbb{H}_{e} C^{-1}_{bf} \mathbb{H}_{f} + // C^{-1}_{ad} \right) + // \\ &= \frac{1}{2} \left[ + // C^{-1}_{ae} \mathbb{H}_{e} C^{-1}_{cf} \mathbb{H}_{f} C^{-1}_{bd} + // + C^{-1}_{ae} \mathbb{H}_{e} C^{-1}_{df} \mathbb{H}_{f} C^{-1}_{bc} + // + C^{-1}_{ce} \mathbb{H}_{e} C^{-1}_{bf} \mathbb{H}_{f} C^{-1}_{ad} + // + C^{-1}_{be} \mathbb{H}_{e} C^{-1}_{df} \mathbb{H}_{f} C^{-1}_{ac} + // \right] + // @f} + // + // Well, that escalated quickly -- although the the definition of $\psi_{0}$ + // and $f_{\mu_e}$ might have given some hints that the calculating + // the kinetic fields and their linearization would take some effort, it is + // likely that there's a little more complexity to the final definitions + // that perhaps initially thought. + // Knowing what we now do, it's probably fair to say that we really do not + // want to compute first and second derivatives of these functions with + // respect to their arguments -- regardless of well we did in calculus + // classes, or how good a programmer we may be. + // + // In the class method definition where these are ultimately implemented, + // we've composed these calculations slightly differently. Some intermediate + // steps are also retained to give another perspective of how to + // systematically compute the derivatives. Additionally, some calculations + // are decomposed less or further to reuse some of the intermediate values + // and, hopefully, aid the reader to follow the derivative operations. + template + class Magnetoelastic_Constitutive_Law final + : public Coupled_Magnetomechanical_Constitutive_Law_Base + { + public: + Magnetoelastic_Constitutive_Law( + const ConstitutiveParameters &constitutive_parameters); + + virtual void update_internal_data(const SymmetricTensor<2, dim> &C, + const Tensor<1, dim> & H, + const DiscreteTime &) override; + + virtual double get_psi() const override; + + virtual Tensor<1, dim> get_B() const override; + + virtual SymmetricTensor<2, dim> get_S() const override; + + virtual SymmetricTensor<2, dim> get_DD() const override; + + virtual Tensor<3, dim> get_PP() const override; + + virtual SymmetricTensor<4, dim> get_HH() const override; + + private: + double psi; + Tensor<1, dim> B; + SymmetricTensor<2, dim> S; + SymmetricTensor<2, dim> BB; + Tensor<3, dim> PP; + SymmetricTensor<4, dim> HH; + }; + + + template + Magnetoelastic_Constitutive_Law::Magnetoelastic_Constitutive_Law( + const ConstitutiveParameters &constitutive_parameters) + : Coupled_Magnetomechanical_Constitutive_Law_Base( + constitutive_parameters) + , psi(0.0) + {} + + // For this class's update method, we'll simply precompute a collection of + // intermediate values (for function evaluations, derivative calculations, + // and the like) and "manually" arrange them in the order that's required + // to maximize their reuse. This means that we have to manage this + // ourselves, and decide what values must be compute before others, all + // while keeping some semblance of order or structure in the code itself. + // It's effective, but perhaps a little tedious. It also doesn't do too much + // to help future extension of the class, because all of these values remain + // local to this single method. + // + // Interestingly, this basic technique of precomputing intermediate + // expressions that are used in more than one place has a name: + // [common subexpression elimination + // (CSE)](https://en.wikipedia.org/wiki/Common_subexpression_elimination). + // It is a strategy used by Computer Algebra Systems to reduce the + // computational expense when they are tasked with evaluating similar + // expressions. + template + void Magnetoelastic_Constitutive_Law::update_internal_data( + const SymmetricTensor<2, dim> &C, + const Tensor<1, dim> & H, + const DiscreteTime &) + { + const double det_F = std::sqrt(determinant(C)); + const SymmetricTensor<2, dim> C_inv = invert(C); + AssertThrow(det_F > 0.0, + ExcMessage("Volumetric Jacobian must be positive.")); + + // The saturation function for the magneto-elastic energy. + const double two_h_dot_h_div_h_sat_squ = + (2.0 * H * H) / (this->get_mu_e_h_sat() * this->get_mu_e_h_sat()); + const double tanh_two_h_dot_h_div_h_sat_squ = + std::tanh(two_h_dot_h_div_h_sat_squ); + + const double f_mu_e = + 1.0 + (this->get_mu_e_inf() / this->get_mu_e() - 1.0) * + tanh_two_h_dot_h_div_h_sat_squ; + + // The first derivative of the saturation function, noting that + // $\frac{d \tanh(x)}{dx} = \text{sech}^{2}(x)$. + const double dtanh_two_h_dot_h_div_h_sat_squ = + std::pow(1.0 / std::cosh(two_h_dot_h_div_h_sat_squ), 2.0); + const Tensor<1, dim> dtwo_h_dot_h_div_h_sat_squ_dH = + 2.0 * 2.0 / (this->get_mu_e_h_sat() * this->get_mu_e_h_sat()) * H; + + const Tensor<1, dim> df_mu_e_dH = + (this->get_mu_e_inf() / this->get_mu_e() - 1.0) * + (dtanh_two_h_dot_h_div_h_sat_squ * dtwo_h_dot_h_div_h_sat_squ_dH); + + // The second derivative of saturation function, noting that + // $\frac{d \text{sech}^{2}(x)}{dx} = -2 \tanh(x) \text{sech}^{2}(x)$. + const double d2tanh_two_h_dot_h_div_h_sat_squ = + -2.0 * tanh_two_h_dot_h_div_h_sat_squ * dtanh_two_h_dot_h_div_h_sat_squ; + const SymmetricTensor<2, dim> d2two_h_dot_h_div_h_sat_squ_dH_dH = + 2.0 * 2.0 / (this->get_mu_e_h_sat() * this->get_mu_e_h_sat()) * + Physics::Elasticity::StandardTensors::I; + + const SymmetricTensor<2, dim> d2f_mu_e_dH_dH = + (this->get_mu_e_inf() / this->get_mu_e() - 1.0) * + (d2tanh_two_h_dot_h_div_h_sat_squ * + symmetrize(outer_product(dtwo_h_dot_h_div_h_sat_squ_dH, + dtwo_h_dot_h_div_h_sat_squ_dH)) + + dtanh_two_h_dot_h_div_h_sat_squ * d2two_h_dot_h_div_h_sat_squ_dH_dH); + + // Some intermediate quantities attained directly from the + // field / kinematic variables. + const double log_det_F = std::log(det_F); + const double tr_C = trace(C); + const Tensor<1, dim> C_inv_dot_H = C_inv * H; + const double H_dot_C_inv_dot_H = H * C_inv_dot_H; + + // First derivatives of the intermediate quantities. + const SymmetricTensor<2, dim> d_tr_C_dC = + Physics::Elasticity::StandardTensors::I; + const SymmetricTensor<2, dim> ddet_F_dC = 0.5 * det_F * C_inv; + const SymmetricTensor<2, dim> dlog_det_F_dC = 0.5 * C_inv; + + const Tensor<1, dim> dH_dot_C_inv_dot_H_dH = 2.0 * C_inv_dot_H; + + SymmetricTensor<4, dim> dC_inv_dC; + for (unsigned int A = 0; A < dim; ++A) + for (unsigned int B = A; B < dim; ++B) + for (unsigned int C = 0; C < dim; ++C) + for (unsigned int D = C; D < dim; ++D) + dC_inv_dC[A][B][C][D] -= // + 0.5 * (C_inv[A][C] * C_inv[B][D] // + + C_inv[A][D] * C_inv[B][C]); // + + const SymmetricTensor<2, dim> dH_dot_C_inv_dot_H_dC = + -symmetrize(outer_product(C_inv_dot_H, C_inv_dot_H)); + + // Second derivatives of the intermediate quantities. + const SymmetricTensor<4, dim> d2log_det_F_dC_dC = 0.5 * dC_inv_dC; + + const SymmetricTensor<4, dim> d2det_F_dC_dC = + 0.5 * (outer_product(C_inv, ddet_F_dC) + det_F * dC_inv_dC); + + const SymmetricTensor<2, dim> d2H_dot_C_inv_dot_H_dH_dH = 2.0 * C_inv; + + Tensor<3, dim> d2H_dot_C_inv_dot_H_dC_dH; + for (unsigned int A = 0; A < dim; ++A) + for (unsigned int B = 0; B < dim; ++B) + for (unsigned int C = 0; C < dim; ++C) + d2H_dot_C_inv_dot_H_dC_dH[A][B][C] -= + C_inv[A][C] * C_inv_dot_H[B] + // + C_inv_dot_H[A] * C_inv[B][C]; // + + SymmetricTensor<4, dim> d2H_dot_C_inv_dot_H_dC_dC; + for (unsigned int A = 0; A < dim; ++A) + for (unsigned int B = A; B < dim; ++B) + for (unsigned int C = 0; C < dim; ++C) + for (unsigned int D = C; D < dim; ++D) + d2H_dot_C_inv_dot_H_dC_dC[A][B][C][D] += + 0.5 * (C_inv_dot_H[A] * C_inv_dot_H[C] * C_inv[B][D] + + C_inv_dot_H[A] * C_inv_dot_H[D] * C_inv[B][C] + + C_inv_dot_H[B] * C_inv_dot_H[C] * C_inv[A][D] + + C_inv_dot_H[B] * C_inv_dot_H[D] * C_inv[A][C]); + + // The stored energy density function. + psi = + (0.5 * this->get_mu_e() * f_mu_e) * + (tr_C - dim - 2.0 * std::log(det_F)) + + this->get_lambda_e() * (std::log(det_F) * std::log(det_F)) - + (0.5 * this->get_mu_0() * this->get_mu_r()) * det_F * (H * C_inv * H); + + // The kinetic quantities. + B = -(0.5 * this->get_mu_e() * (tr_C - dim - 2.0 * log_det_F)) * + df_mu_e_dH // + + 0.5 * this->get_mu_0() * this->get_mu_r() * det_F * + dH_dot_C_inv_dot_H_dH; // + + S = 2.0 * (0.5 * this->get_mu_e() * f_mu_e) * // + (d_tr_C_dC - 2.0 * dlog_det_F_dC) // + + 2.0 * this->get_lambda_e() * (2.0 * log_det_F * dlog_det_F_dC) // + - 2.0 * (0.5 * this->get_mu_0() * this->get_mu_r()) * // + (H_dot_C_inv_dot_H * ddet_F_dC // + + det_F * dH_dot_C_inv_dot_H_dC); // + + // The linearization of the kinetic quantities. + BB = -(0.5 * this->get_mu_e() * (tr_C - dim - 2.0 * log_det_F)) * // + d2f_mu_e_dH_dH // + + 0.5 * this->get_mu_0() * this->get_mu_r() * det_F * + d2H_dot_C_inv_dot_H_dH_dH; // + + PP = -2.0 * (0.5 * this->get_mu_e()) * // + outer_product(Tensor<2, dim>(d_tr_C_dC - 2.0 * dlog_det_F_dC), // + df_mu_e_dH) // + + // + 2.0 * (0.5 * this->get_mu_0() * this->get_mu_r()) * // + (outer_product(Tensor<2, dim>(ddet_F_dC), dH_dot_C_inv_dot_H_dH) // + + det_F * d2H_dot_C_inv_dot_H_dC_dH); // + + HH = + 4.0 * (0.5 * this->get_mu_e() * f_mu_e) * (-2.0 * d2log_det_F_dC_dC) // + + 4.0 * this->get_lambda_e() * // + (2.0 * outer_product(dlog_det_F_dC, dlog_det_F_dC) // + + 2.0 * log_det_F * d2log_det_F_dC_dC) // + - 4.0 * (0.5 * this->get_mu_0() * this->get_mu_r()) * // + (H_dot_C_inv_dot_H * d2det_F_dC_dC // + + outer_product(ddet_F_dC, dH_dot_C_inv_dot_H_dC) // + + outer_product(dH_dot_C_inv_dot_H_dC, ddet_F_dC) // + + det_F * d2H_dot_C_inv_dot_H_dC_dC); // + } + + template + double Magnetoelastic_Constitutive_Law::get_psi() const + { + return psi; + } + + template + Tensor<1, dim> Magnetoelastic_Constitutive_Law::get_B() const + { + return B; + } + + template + SymmetricTensor<2, dim> Magnetoelastic_Constitutive_Law::get_S() const + { + return S; + } + + template + SymmetricTensor<2, dim> Magnetoelastic_Constitutive_Law::get_DD() const + { + return BB; + } + + template + Tensor<3, dim> Magnetoelastic_Constitutive_Law::get_PP() const + { + return PP; + } + + template + SymmetricTensor<4, dim> Magnetoelastic_Constitutive_Law::get_HH() const + { + return HH; + } + + + // @sect4{Magneto-viscoelastic constitutive law (hand-derived)} + + // As mentioned before, the free energy density function for the + // magneto-viscoelastic material with one dissipative mechanism that we'll + // be considering is defined as + // @f[ + // \psi_{0} \left( \mathbf{C}, \mathbf{C}_{v}, \boldsymbol{\mathbb{H}} + // \right) + // = \psi_{0}^{ME} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) + // + \psi_{0}^{MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) + // @f] + // @f[ + // \psi_{0}^{ME} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) + // = \frac{1}{2} \mu_{e} f_{\mu_{e}^{ME}} \left( \boldsymbol{\mathbb{H}} + // \right) + // \left[ \text{tr}(\mathbf{C}) - d - 2 \ln (\text{det}(\mathbf{F})) + // \right] + // + \lambda_{e} \ln^{2} \left(\text{det}(\mathbf{F}) \right) + // - \frac{1}{2} \mu_{0} \mu_{r} \text{det}(\mathbf{F}) + // \left[ \boldsymbol{\mathbb{H}} \cdot \mathbf{C}^{-1} \cdot + // \boldsymbol{\mathbb{H}} \right] + // @f] + // @f[ + // \psi_{0}^{MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) + // = \frac{1}{2} \mu_{v} f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} + // \right) + // \left[ \mathbf{C}_{v} : \left[ + // \left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + // \mathbf{C} \right] - d - \ln\left( + // \text{det}\left(\mathbf{C}_{v}\right) \right) \right] + // @f] + // with + // @f[ + // f_{\mu_{e}}^{ME} \left( \boldsymbol{\mathbb{H}} \right) + // = 1 + \left[ \frac{\mu_{e}^{\infty}}{\mu_{e}} - 1 \right] + // \tanh \left( 2 \frac{\boldsymbol{\mathbb{H}} \cdot + // \boldsymbol{\mathbb{H}}} + // {\left(h_{e}^{\text{sat}}\right)^{2}} \right) + // @f] + // @f[ + // f_{\mu_{v}}^{MVE} \left( \boldsymbol{\mathbb{H}} \right) + // = 1 + \left[ \frac{\mu_{v}^{\infty}}{\mu_{v}} - 1 \right] + // \tanh \left( 2 \frac{\boldsymbol{\mathbb{H}} \cdot + // \boldsymbol{\mathbb{H}}} + // {\left(h_{v}^{\text{sat}}\right)^{2}} \right) + // @f] + // and the evolution law + // @f[ + // \dot{\mathbf{C}}_{v} \left( \mathbf{C} \right) + // = \frac{1}{\tau} \left[ + // \left[\left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + // \mathbf{C}\right]^{-1} + // - \mathbf{C}_{v} \right] + // @f] + // that itself is parameterized in terms of $\mathbf{C}$. + // By design, the magnetoelastic part of the energy + // $\psi_{0}^{ME} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right)$ + // is identical to that of the magnetoelastic material presented earlier. + // So, for the derivatives of the various contributions stemming from this + // part of the energy, please refer to the previous section. We'll continue + // to highlight the specific contributions from those terms by + // superscripting the salient terms with $ME$, while contributions from the + // magneto-viscoelastic component are superscripted with $MVE$. + // Furthermore, the magnetic saturation function + // $f_{\mu_{v}}^{MVE} \left( \boldsymbol{\mathbb{H}} \right)$ + // for the damping term has the identical form as that of the elastic + // term (i.e., + // $f_{\mu_{e}}^{ME} \left( \boldsymbol{\mathbb{H}} \right)$ + // ), and so the structure of its derivatives are identical to that + // seen before; the only change is for the three constitutive parameters + // that are now associated with the viscous shear modulus $\mu_{v}$ rather + // than the elastic shear modulus $\mu_{e}$. + // + // For this magneto-viscoelastic material, the first derivatives that + // correspond to the magnetic induction vector and total Piola-Kirchhoff + // stress tensor are + // @f[ + // \boldsymbol{\mathbb{B}} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) + // \dealcoloneq - \frac{\partial \psi_{0} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right)}{\partial \boldsymbol{\mathbb{H}}} + // \Big\vert_{\mathbf{C}, \mathbf{C}_{v}} \equiv + // \boldsymbol{\mathbb{B}}^{ME} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} + // \right) + // + \boldsymbol{\mathbb{B}}^{MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) = - \frac{d \psi_{0}^{ME} \left( + // \mathbf{C}, \boldsymbol{\mathbb{H}} \right)}{d \boldsymbol{\mathbb{H}}} + // - \frac{\partial \psi_{0}^{MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right)}{\partial \boldsymbol{\mathbb{H}}} + // @f] + // @f[ + // \mathbf{S}^{\text{tot}} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) + // \dealcoloneq 2 \frac{\partial \psi_{0} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right)}{\partial \mathbf{C}} + // \Big\vert_{\mathbf{C}_{v}, \boldsymbol{\mathbb{H}}} \equiv + // \mathbf{S}^{\text{tot}, ME} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} + // \right) + // + \mathbf{S}^{\text{tot}, MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} + // \right) + // = 2 \frac{d \psi_{0}^{ME} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} + // \right)}{d \mathbf{C}} + // + 2 \frac{\partial \psi_{0}^{MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right)}{\partial \mathbf{C}} + // @f] + // with the viscous contributions being + // @f[ + // \boldsymbol{\mathbb{B}}^{MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) + // = - \frac{\partial \psi_{0}^{MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right)}{\partial \boldsymbol{\mathbb{H}}} + // \Big\vert_{\mathbf{C}, \mathbf{C}_{v}} = - \frac{1}{2} \mu_{v} + // \left[ \mathbf{C}_{v} : \left[ + // \left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + // \mathbf{C} \right] - d - \ln\left( + // \text{det}\left(\mathbf{C}_{v}\right) \right) \right] + // \frac{\partial f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} + // \right)}{\partial \boldsymbol{\mathbb{H}}} + // @f] + // @f[ + // \mathbf{S}^{\text{tot}, MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} + // \right) + // = 2 \frac{\partial \psi_{0}^{MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right)}{\partial \mathbf{C}} + // \Big\vert_{\mathbf{C}_{v}, \boldsymbol{\mathbb{H}}} = \mu_{v} + // f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} \right) + // \left[ \left[ \mathbf{C}_{v} : \mathbf{C} \right] \left[ - + // \frac{1}{d} + // \left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + // \mathbf{C}^{-1} \right] + // + \left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + // \mathbf{C}_{v} + // \right] + // @f] + // and with + // @f[ + // \frac{\partial f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} + // \right)}{\partial \boldsymbol{\mathbb{H}}} \equiv \frac{d + // f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} \right)}{d + // \boldsymbol{\mathbb{H}}} . + // @f] + // The time-discretized evolution law, + // @f[ + // \mathbf{C}_{v}^{(t)} \left( \mathbf{C} \right) + // = \frac{1}{1 + \frac{\Delta t}{\tau_{v}}} \left[ + // \mathbf{C}_{v}^{(t-1)} + // + \frac{\Delta t}{\tau_{v}} + // \left[\left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + // \mathbf{C} \right]^{-1} + // \right] + // @f] + // will also dictate how the linearization of the internal + // variable with respect to the field variables is composed. + // + // Observe that in order to attain the *correct* expressions for the + // magnetic induction vector and total Piola-Kirchhoff stress tensor for + // this dissipative material, we must adhere strictly to the outcome of + // applying the Coleman-Noll procedure: we must take *partial derivatives* + // of the free energy density function with respect to the field variables. + // (For out non-dissipative magnetoelastic material, taking either partial + // or total derivatives would have had the same result, so there was no + // need to draw your attention to this before.) + // The crucial part of the operation is to freeze the internal variable + // $\mathbf{C}_{v}^{(t)} \left( \mathbf{C} \right)$ while computing the + // derivatives of $\psi_{0}^{MVE} \left( \mathbf{C}, + // \mathbf{C}_{v} \left( \mathbf{C} \right), \boldsymbol{\mathbb{H}} + // \right)$ with respect to $\mathbf{C}$ -- the dependence of + // $\mathbf{C}_{v}^{(t)}$ on $\mathbf{C}$ is not to be taken into account. + // When deciding whether to use AD or SD to perform this task + // the choice is clear -- only the symbolic framework provides a mechanism + // to do this; as was mentioned before, AD can only return total derivatives + // so it is unsuitable for the task. + // + // To wrap things up, we'll present the material tangents for this + // rate-dependent coupled material. The linearization of both kinetic + // variables with respect to their arguments are + // @f[ + // \mathbb{D} \left( \mathbf{C}, \mathbf{C}_{v}, \boldsymbol{\mathbb{H}} + // \right) = \frac{d \boldsymbol{\mathbb{B}}}{d \boldsymbol{\mathbb{H}}} + // \equiv \mathbb{D}^{ME} \left( \mathbf{C}, \boldsymbol{\mathbb{H}} \right) + // + \mathbb{D}^{MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) = \frac{d \boldsymbol{\mathbb{B}}^{ME}}{d + // \boldsymbol{\mathbb{H}}} + // + \frac{d \boldsymbol{\mathbb{B}}^{MVE}}{d \boldsymbol{\mathbb{H}}} + // @f] + // @f[ + // \mathfrak{P}^{\text{tot}} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) = - \frac{d \mathbf{S}^{\text{tot}}}{d + // \boldsymbol{\mathbb{H}}} \equiv \mathfrak{P}^{\text{tot}, ME} \left( + // \mathbf{C}, \boldsymbol{\mathbb{H}} \right) + // + \mathfrak{P}^{\text{tot}, MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) = - \frac{d \mathbf{S}^{\text{tot}, + // ME}}{d \boldsymbol{\mathbb{H}}} + // - \frac{d \mathbf{S}^{\text{tot}, MVE}}{d \boldsymbol{\mathbb{H}}} + // @f] + // @f[ + // \mathcal{H}^{\text{tot}} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) = 2 \frac{d \mathbf{S}^{\text{tot}}}{d + // \mathbf{C}} \equiv \mathcal{H}^{\text{tot}, ME} \left( \mathbf{C}, + // \boldsymbol{\mathbb{H}} \right) + // + \mathcal{H}^{\text{tot}, MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) = 2 \frac{d \mathbf{S}^{\text{tot}, + // ME}}{d \mathbf{C}} + // + 2 \frac{d \mathbf{S}^{\text{tot}, MVE}}{d \mathbf{C}} + // @f] + // where the tangents for the viscous contributions are + // @f[ + // \mathbb{D}^{MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) = - \frac{1}{2} \mu_{v} + // \left[ \mathbf{C}_{v} : \left[ + // \left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + // \mathbf{C} \right] - d - \ln\left( + // \text{det}\left(\mathbf{C}_{v}\right) \right) \right] + // \frac{\partial^{2} f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} + // \right)}{\partial \boldsymbol{\mathbb{H}} \otimes + // d \boldsymbol{\mathbb{H}}} + // @f] + // @f[ + // \mathfrak{P}^{\text{tot}, MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) = - \mu_{v} + // \left[ \left[ \mathbf{C}_{v} : \mathbf{C} \right] \left[ - + // \frac{1}{d} + // \left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + // \mathbf{C}^{-1} \right] + // + \left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + // \mathbf{C}_{v} + // \right] \otimes \frac{d f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} + // \right)}{d \boldsymbol{\mathbb{H}}} + // @f] + // @f{align} + // \mathcal{H}^{\text{tot}, MVE} \left( \mathbf{C}, \mathbf{C}_{v}, + // \boldsymbol{\mathbb{H}} \right) + // &= 2 \mu_{v} f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} \right) + // \left[ - \frac{1}{d} + // \left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + // \mathbf{C}^{-1} \right] \otimes + // \left[ \mathbf{C}_{v} + \mathbf{C} : \frac{d \mathbf{C}_{v}}{d + // \mathbf{C}} \right] + // \\ &+ 2 \mu_{v} f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} \right) + // \left[ \mathbf{C}_{v} : \mathbf{C} \right] + // \left[ + // \frac{1}{d^{2}} + // \left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + // \mathbf{C}^{-1} \otimes \mathbf{C}^{-1} + // - \frac{1}{d} + // \left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} \frac{d + // \mathbf{C}^{-1}}{d \mathbf{C}} + // \right] + // \\ &+ 2 \mu_{v} f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} \right) + // \left[ + // -\frac{1}{d} + // \left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + // \mathbf{C}_{v} \otimes \mathbf{C}^{-1} + // + \left[\text{det}\left(\mathbf{F}\right)\right]^{-\frac{2}{d}} + // \frac{d \mathbf{C}_{v}}{d \mathbf{C}} + // \right] + // @f} + // with + // @f[ + // \frac{\partial^{2} f_{\mu_{v}^{MVE}} \left( \boldsymbol{\mathbb{H}} + // \right)}{\partial \boldsymbol{\mathbb{H}} \otimes + // d \boldsymbol{\mathbb{H}}} \equiv \frac{d^{2} f_{\mu_{v}^{MVE}} \left( + // \boldsymbol{\mathbb{H}} \right)}{d \boldsymbol{\mathbb{H}} \otimes d + // \boldsymbol{\mathbb{H}}} + // @f] + // and, from the evolution law, + // @f[ + // \frac{d \mathbf{C}_{v}}{d \mathbf{C}} + // \equiv \frac{d \mathbf{C}_{v}^{(t)}}{d \mathbf{C}} + // = \frac{\frac{\Delta t}{\tau_{v}} }{1 + \frac{\Delta t}{\tau_{v}}} + // \left[ + // \frac{1}{d} + // \left[\text{det}\left(\mathbf{F}\right)\right]^{\frac{2}{d}} + // \mathbf{C}^{-1} \otimes \mathbf{C}^{-1} + // + \left[\text{det}\left(\mathbf{F}\right)\right]^{\frac{2}{d}} \frac{d + // \mathbf{C}^{-1}}{d \mathbf{C}} + // \right] . + // @f] + // Notice that just the last term of $\mathcal{H}^{\text{tot}, MVE}$ + // contains the tangent of the internal variable. The linearization of this + // particular evolution law is linear. For an example of a nonlinear + // evolution law, for which this linearization must be solved for in an + // iterative manner, see @cite Koprowski-Theiss2011a. + template + class Magnetoviscoelastic_Constitutive_Law final + : public Coupled_Magnetomechanical_Constitutive_Law_Base + { + public: + Magnetoviscoelastic_Constitutive_Law( + const ConstitutiveParameters &constitutive_parameters); + + virtual void update_internal_data(const SymmetricTensor<2, dim> &C, + const Tensor<1, dim> & H, + const DiscreteTime &time) override; + + virtual double get_psi() const override; + + virtual Tensor<1, dim> get_B() const override; + + virtual SymmetricTensor<2, dim> get_S() const override; + + virtual SymmetricTensor<2, dim> get_DD() const override; + + virtual Tensor<3, dim> get_PP() const override; + + virtual SymmetricTensor<4, dim> get_HH() const override; + + virtual void update_end_of_timestep() override; + + private: + SymmetricTensor<2, dim> Q_t; + SymmetricTensor<2, dim> Q_t1; + + double psi; + Tensor<1, dim> B; + SymmetricTensor<2, dim> S; + SymmetricTensor<2, dim> BB; + Tensor<3, dim> PP; + SymmetricTensor<4, dim> HH; + + // A data structure that is used to store all intermediate calculations. + // We'll see shortly precisely how this can be leveraged to make the part + // of the code where we actually perform calculations clean and easy + // (well, at least easier) to follow and maintain. But for now, we can say + // that it will allow us to move the parts of the code where we compute + // the derivatives of intermediate quantities away from where they are + // used. + mutable GeneralDataStorage cache; + + // The next two functions are used to update the state of the field and + // internal variables, and will be called before we perform any + // detailed calculations. + void set_primary_variables(const SymmetricTensor<2, dim> &C, + const Tensor<1, dim> & H) const; + + void update_internal_variable(const DiscreteTime &time); + + // The remainder of the class interface is dedicated to methods that + // are used to compute the components required to calculate the free + // energy density function, and all of its derivatives: + + // The kinematic, or field, variables. + const Tensor<1, dim> &get_H() const; + + const SymmetricTensor<2, dim> &get_C() const; + + // A generalized formulation for the saturation function, with the + // required constitutive parameters passed as arguments to each function. + double get_two_h_dot_h_div_h_sat_squ(const double mu_h_sat) const; + + double get_tanh_two_h_dot_h_div_h_sat_squ(const double mu_h_sat) const; + + double get_f_mu(const double mu, + const double mu_inf, + const double mu_h_sat) const; + + // A generalized formulation for the first derivative of saturation + // function, with the required constitutive parameters passed as arguments + // to each function. + double get_dtanh_two_h_dot_h_div_h_sat_squ(const double mu_h_sat) const; + + Tensor<1, dim> + get_dtwo_h_dot_h_div_h_sat_squ_dH(const double mu_h_sat) const; + + Tensor<1, dim> get_df_mu_dH(const double mu, + const double mu_inf, + const double mu_h_sat) const; + + // A generalized formulation for the second derivative of saturation + // function, with the required constitutive parameters passed as arguments + // to each function. + double get_d2tanh_two_h_dot_h_div_h_sat_squ(const double mu_h_sat) const; + + SymmetricTensor<2, dim> + get_d2two_h_dot_h_div_h_sat_squ_dH_dH(const double mu_h_sat) const; + + SymmetricTensor<2, dim> get_d2f_mu_dH_dH(const double mu, + const double mu_inf, + const double mu_h_sat) const; + + // Intermediate quantities attained directly from the + // field / kinematic variables. + const double &get_det_F() const; + + const SymmetricTensor<2, dim> &get_C_inv() const; + + const double &get_log_det_F() const; + + const double &get_trace_C() const; + + const Tensor<1, dim> &get_C_inv_dot_H() const; + + const double &get_H_dot_C_inv_dot_H() const; + + // First derivatives of the intermediate quantities. + const SymmetricTensor<4, dim> &get_dC_inv_dC() const; + + const SymmetricTensor<2, dim> &get_d_tr_C_dC() const; + + const SymmetricTensor<2, dim> &get_ddet_F_dC() const; + + const SymmetricTensor<2, dim> &get_dlog_det_F_dC() const; + + const Tensor<1, dim> &get_dH_dot_C_inv_dot_H_dH() const; + + const SymmetricTensor<2, dim> &get_dH_dot_C_inv_dot_H_dC() const; + + // Derivative of internal variable with respect to field variables. + // Notice that we only need this one derivative of the internal variable, + // as this variable is only differentiated as part of the linearization + // of the kinetic variables. + const SymmetricTensor<4, dim> & + get_dQ_t_dC(const DiscreteTime &time) const; + + // Second derivatives of the intermediate quantities. + const SymmetricTensor<4, dim> &get_d2log_det_F_dC_dC() const; + + const SymmetricTensor<4, dim> &get_d2det_F_dC_dC() const; + + const SymmetricTensor<2, dim> &get_d2H_dot_C_inv_dot_H_dH_dH() const; + + const Tensor<3, dim> &get_d2H_dot_C_inv_dot_H_dC_dH() const; + + const SymmetricTensor<4, dim> &get_d2H_dot_C_inv_dot_H_dC_dC() const; + }; + + + template + Magnetoviscoelastic_Constitutive_Law< + dim>::Magnetoviscoelastic_Constitutive_Law(const ConstitutiveParameters + &constitutive_parameters) + : Coupled_Magnetomechanical_Constitutive_Law_Base( + constitutive_parameters) + , Q_t(Physics::Elasticity::StandardTensors::I) + , Q_t1(Physics::Elasticity::StandardTensors::I) + , psi(0.0) + {} + + + template + void Magnetoviscoelastic_Constitutive_Law::update_internal_data( + const SymmetricTensor<2, dim> &C, + const Tensor<1, dim> & H, + const DiscreteTime & time) + { + // Record the applied deformation state as well as the magnetic load. + // Thereafter, update internal (viscous) variable based on new deformation + // state. + set_primary_variables(C, H); + update_internal_variable(time); + + // Get the values for the elastic and viscous saturation function based + // on the current magnetic field... + const double f_mu_e = get_f_mu(this->get_mu_e(), + this->get_mu_e_inf(), + this->get_mu_e_h_sat()); + + const double f_mu_v = get_f_mu(this->get_mu_v(), + this->get_mu_v_inf(), + this->get_mu_v_h_sat()); + + // ... as well as their first derivatives... + const Tensor<1, dim> df_mu_e_dH = get_df_mu_dH(this->get_mu_e(), + this->get_mu_e_inf(), + this->get_mu_e_h_sat()); + + const Tensor<1, dim> df_mu_v_dH = get_df_mu_dH(this->get_mu_v(), + this->get_mu_v_inf(), + this->get_mu_v_h_sat()); + + + // ... and their second derivatives. + const SymmetricTensor<2, dim> d2f_mu_e_dH_dH = + get_d2f_mu_dH_dH(this->get_mu_e(), + this->get_mu_e_inf(), + this->get_mu_e_h_sat()); + + const SymmetricTensor<2, dim> d2f_mu_v_dH_dH = + get_d2f_mu_dH_dH(this->get_mu_v(), + this->get_mu_v_inf(), + this->get_mu_v_h_sat()); + + // Intermediate quantities. Note that, since we're fetching these values + // from a cache that has a lifetime that outlasts this function call, we + // can alias the result rather than copying the value from the cache. + const double & det_F = get_det_F(); + const SymmetricTensor<2, dim> &C_inv = get_C_inv(); + + const double &log_det_F = get_log_det_F(); + const double &tr_C = get_trace_C(); + const double &H_dot_C_inv_dot_H = get_H_dot_C_inv_dot_H(); + + // First derivatives of intermediate values, as well as the that of the + // internal variable with respect to the right Cauchy-Green deformation + // tensor. + const SymmetricTensor<2, dim> &d_tr_C_dC = get_d_tr_C_dC(); + const SymmetricTensor<2, dim> &ddet_F_dC = get_ddet_F_dC(); + const SymmetricTensor<2, dim> &dlog_det_F_dC = get_dlog_det_F_dC(); + + const SymmetricTensor<4, dim> &dQ_t_dC = get_dQ_t_dC(time); + + const Tensor<1, dim> &dH_dot_C_inv_dot_H_dH = get_dH_dot_C_inv_dot_H_dH(); + + const SymmetricTensor<2, dim> &dH_dot_C_inv_dot_H_dC = + get_dH_dot_C_inv_dot_H_dC(); + + // Second derivatives of intermediate values. + const SymmetricTensor<4, dim> &d2log_det_F_dC_dC = + get_d2log_det_F_dC_dC(); + + const SymmetricTensor<4, dim> &d2det_F_dC_dC = get_d2det_F_dC_dC(); + + const SymmetricTensor<2, dim> &d2H_dot_C_inv_dot_H_dH_dH = + get_d2H_dot_C_inv_dot_H_dH_dH(); + + const Tensor<3, dim> &d2H_dot_C_inv_dot_H_dC_dH = + get_d2H_dot_C_inv_dot_H_dC_dH(); + + const SymmetricTensor<4, dim> &d2H_dot_C_inv_dot_H_dC_dC = + get_d2H_dot_C_inv_dot_H_dC_dC(); + + // Since the definitions of the linearizations become particularly + // lengthy, we'll decompose the free energy density function into three + // additive components: + // - the "Neo-Hookean"-like term, + // - the rate-dependent term, and + // - the term that resembles that of the energy stored in the magnetic + // field. + // + // To remain consistent, each of these contributions will be individually + // added to the variables that we want to compute in that same order. + // + // So, first of all this is the energy density function itself: + psi = (0.5 * this->get_mu_e() * f_mu_e) * + (tr_C - dim - 2.0 * std::log(det_F)) + + this->get_lambda_e() * (std::log(det_F) * std::log(det_F)); + psi += (0.5 * this->get_mu_v() * f_mu_v) * + (Q_t * (std::pow(det_F, -2.0 / dim) * C) - dim - + std::log(determinant(Q_t))); + psi -= + (0.5 * this->get_mu_0() * this->get_mu_r()) * det_F * (H * C_inv * H); + + // ... followed by the magnetic induction vector and Piola-Kirchhoff + // stress: + B = + -(0.5 * this->get_mu_e() * (tr_C - dim - 2.0 * log_det_F)) * df_mu_e_dH; + B -= (0.5 * this->get_mu_v()) * + (Q_t * (std::pow(det_F, -2.0 / dim) * C) - dim - + std::log(determinant(Q_t))) * + df_mu_v_dH; + B += 0.5 * this->get_mu_0() * this->get_mu_r() * det_F * + dH_dot_C_inv_dot_H_dH; + + S = 2.0 * (0.5 * this->get_mu_e() * f_mu_e) * // + (d_tr_C_dC - 2.0 * dlog_det_F_dC) // + + 2.0 * this->get_lambda_e() * (2.0 * log_det_F * dlog_det_F_dC); // + S += 2.0 * (0.5 * this->get_mu_v() * f_mu_v) * + ((Q_t * C) * + ((-2.0 / dim) * std::pow(det_F, -2.0 / dim - 1.0) * ddet_F_dC) + + std::pow(det_F, -2.0 / dim) * Q_t); // dC/dC = II + S -= 2.0 * (0.5 * this->get_mu_0() * this->get_mu_r()) * // + (H_dot_C_inv_dot_H * ddet_F_dC // + + det_F * dH_dot_C_inv_dot_H_dC); // + + // ... and lastly the tangents due to the linearization of the kinetic + // variables. + BB = -(0.5 * this->get_mu_e() * (tr_C - dim - 2.0 * log_det_F)) * + d2f_mu_e_dH_dH; + BB -= (0.5 * this->get_mu_v()) * + (Q_t * (std::pow(det_F, -2.0 / dim) * C) - dim - + std::log(determinant(Q_t))) * + d2f_mu_v_dH_dH; + BB += 0.5 * this->get_mu_0() * this->get_mu_r() * det_F * + d2H_dot_C_inv_dot_H_dH_dH; + + PP = -2.0 * (0.5 * this->get_mu_e()) * + outer_product(Tensor<2, dim>(d_tr_C_dC - 2.0 * dlog_det_F_dC), + df_mu_e_dH); + PP -= 2.0 * (0.5 * this->get_mu_v()) * + outer_product(Tensor<2, dim>((Q_t * C) * + ((-2.0 / dim) * + std::pow(det_F, -2.0 / dim - 1.0) * + ddet_F_dC) + + std::pow(det_F, -2.0 / dim) * Q_t), + df_mu_v_dH); + PP += 2.0 * (0.5 * this->get_mu_0() * this->get_mu_r()) * + (outer_product(Tensor<2, dim>(ddet_F_dC), dH_dot_C_inv_dot_H_dH) + + det_F * d2H_dot_C_inv_dot_H_dC_dH); + + HH = + 4.0 * (0.5 * this->get_mu_e() * f_mu_e) * (-2.0 * d2log_det_F_dC_dC) // + + 4.0 * this->get_lambda_e() * // + (2.0 * outer_product(dlog_det_F_dC, dlog_det_F_dC) // + + 2.0 * log_det_F * d2log_det_F_dC_dC); // + HH += 4.0 * (0.5 * this->get_mu_v() * f_mu_v) * + (outer_product((-2.0 / dim) * std::pow(det_F, -2.0 / dim - 1.0) * + ddet_F_dC, + C * dQ_t_dC + Q_t) + + (Q_t * C) * + (outer_product(ddet_F_dC, + (-2.0 / dim) * (-2.0 / dim - 1.0) * + std::pow(det_F, -2.0 / dim - 2.0) * ddet_F_dC) + + ((-2.0 / dim) * std::pow(det_F, -2.0 / dim - 1.0) * + d2det_F_dC_dC)) + + outer_product(Q_t, + (-2.0 / dim) * std::pow(det_F, -2.0 / dim - 1.0) * + ddet_F_dC) + + std::pow(det_F, -2.0 / dim) * dQ_t_dC); + HH -= 4.0 * (0.5 * this->get_mu_0() * this->get_mu_r()) * // + (H_dot_C_inv_dot_H * d2det_F_dC_dC // + + outer_product(ddet_F_dC, dH_dot_C_inv_dot_H_dC) // + + outer_product(dH_dot_C_inv_dot_H_dC, ddet_F_dC) // + + det_F * d2H_dot_C_inv_dot_H_dC_dC); // + + + // Now that we're done using all of those temporary variables stored + // in our cache, we can clear it out to free up some memory. + cache.reset(); + } + + template + double Magnetoviscoelastic_Constitutive_Law::get_psi() const + { + return psi; + } + + + template + Tensor<1, dim> Magnetoviscoelastic_Constitutive_Law::get_B() const + { + return B; + } + + + template + SymmetricTensor<2, dim> + Magnetoviscoelastic_Constitutive_Law::get_S() const + { + return S; + } + + + template + SymmetricTensor<2, dim> + Magnetoviscoelastic_Constitutive_Law::get_DD() const + { + return BB; + } + + + template + Tensor<3, dim> Magnetoviscoelastic_Constitutive_Law::get_PP() const + { + return PP; + } + + + template + SymmetricTensor<4, dim> + Magnetoviscoelastic_Constitutive_Law::get_HH() const + { + return HH; + } + + + template + void Magnetoviscoelastic_Constitutive_Law::update_end_of_timestep() + { + Q_t1 = Q_t; + }; + + + template + void Magnetoviscoelastic_Constitutive_Law::update_internal_variable( + const DiscreteTime &time) + { + const double delta_t = this->get_delta_t(time); + + Q_t = (1.0 / (1.0 + delta_t / this->get_tau_v())) * + (Q_t1 + (delta_t / this->get_tau_v()) * + std::pow(get_det_F(), 2.0 / dim) * get_C_inv()); + } + + // The next few functions implement the generalized formulation for the + // saturation function, as well as its various derivatives. + template + double + Magnetoviscoelastic_Constitutive_Law::get_two_h_dot_h_div_h_sat_squ( + const double mu_h_sat) const + { + const Tensor<1, dim> &H = get_H(); + return (2.0 * H * H) / (mu_h_sat * mu_h_sat); + }; + + + template + double Magnetoviscoelastic_Constitutive_Law< + dim>::get_tanh_two_h_dot_h_div_h_sat_squ(const double mu_h_sat) const + { + return std::tanh(get_two_h_dot_h_div_h_sat_squ(mu_h_sat)); + }; + + // A scaling function that will cause the shear modulus + // to change (increase) under the influence of a magnetic + // field. + template + double Magnetoviscoelastic_Constitutive_Law::get_f_mu( + const double mu, + const double mu_inf, + const double mu_h_sat) const + { + return 1.0 + + (mu_inf / mu - 1.0) * get_tanh_two_h_dot_h_div_h_sat_squ(mu_h_sat); + }; + + // First derivative of scaling function + template + double Magnetoviscoelastic_Constitutive_Law< + dim>::get_dtanh_two_h_dot_h_div_h_sat_squ(const double mu_h_sat) const + { + return std::pow(1.0 / std::cosh(get_two_h_dot_h_div_h_sat_squ(mu_h_sat)), + 2.0); + }; + + + template + Tensor<1, dim> Magnetoviscoelastic_Constitutive_Law< + dim>::get_dtwo_h_dot_h_div_h_sat_squ_dH(const double mu_h_sat) const + { + return 2.0 * 2.0 / (mu_h_sat * mu_h_sat) * get_H(); + }; + + + template + Tensor<1, dim> Magnetoviscoelastic_Constitutive_Law::get_df_mu_dH( + const double mu, + const double mu_inf, + const double mu_h_sat) const + { + return (mu_inf / mu - 1.0) * + (get_dtanh_two_h_dot_h_div_h_sat_squ(mu_h_sat) * + get_dtwo_h_dot_h_div_h_sat_squ_dH(mu_h_sat)); + }; + + + template + double Magnetoviscoelastic_Constitutive_Law< + dim>::get_d2tanh_two_h_dot_h_div_h_sat_squ(const double mu_h_sat) const + { + return -2.0 * get_tanh_two_h_dot_h_div_h_sat_squ(mu_h_sat) * + get_dtanh_two_h_dot_h_div_h_sat_squ(mu_h_sat); + }; + + + template + SymmetricTensor<2, dim> Magnetoviscoelastic_Constitutive_Law< + dim>::get_d2two_h_dot_h_div_h_sat_squ_dH_dH(const double mu_h_sat) const + { + return 2.0 * 2.0 / (mu_h_sat * mu_h_sat) * + Physics::Elasticity::StandardTensors::I; + }; + + + template + SymmetricTensor<2, dim> + Magnetoviscoelastic_Constitutive_Law::get_d2f_mu_dH_dH( + const double mu, + const double mu_inf, + const double mu_h_sat) const + { + return (mu_inf / mu - 1.0) * + (get_d2tanh_two_h_dot_h_div_h_sat_squ(mu_h_sat) * + symmetrize( + outer_product(get_dtwo_h_dot_h_div_h_sat_squ_dH(mu_h_sat), + get_dtwo_h_dot_h_div_h_sat_squ_dH(mu_h_sat))) + + get_dtanh_two_h_dot_h_div_h_sat_squ(mu_h_sat) * + get_d2two_h_dot_h_div_h_sat_squ_dH_dH(mu_h_sat)); + }; + + // For the cached calculation approach that we've adopted for this material + // class, the root of all calculations are the field variables, and the + // immutable ancillary data such as the constitutive parameters and time + // step size. As such, we need to enter them into the cache in a different + // manner to the other variables, since they are inputs that are prescribed + // from outside the class itself. This function simply adds them to the + // cache directly from the input arguments, checking that there is no + // equivalent data there in the first place (we expect to call the + // `update_internal_data()` method only once per time step, or Newton + // iteration). + template + void Magnetoviscoelastic_Constitutive_Law::set_primary_variables( + const SymmetricTensor<2, dim> &C, + const Tensor<1, dim> & H) const + { + // Set value for $\boldsymbol{\mathbb{H}}$. + const std::string name_H("H"); + Assert(!cache.stores_object_with_name(name_H), + ExcMessage( + "The primary variable has already been added to the cache.")); + cache.add_unique_copy(name_H, H); + + // Set value for $\mathbf{C}$. + const std::string name_C("C"); + Assert(!cache.stores_object_with_name(name_C), + ExcMessage( + "The primary variable has already been added to the cache.")); + cache.add_unique_copy(name_C, C); + } + + // After that, we can fetch them from the cache at any point in time. + template + const Tensor<1, dim> & + Magnetoviscoelastic_Constitutive_Law::get_H() const + { + const std::string name("H"); + Assert(cache.stores_object_with_name(name), + ExcMessage("Primary variables must be added to the cache.")); + return cache.template get_object_with_name>(name); + } + + template + const SymmetricTensor<2, dim> & + Magnetoviscoelastic_Constitutive_Law::get_C() const + { + const std::string name("C"); + Assert(cache.stores_object_with_name(name), + ExcMessage("Primary variables must be added to the cache.")); + return cache.template get_object_with_name>(name); + } + + // With the primary variables guaranteed to be in the cache when we need + // them, we can not compute all intermediate values (either directly, or + // indirectly) from them. + // + // If the cache does not already store the value that we're looking for, + // then we quickly calculate it, store it in the cache and return the value + // just stored in the cache. That way we can return it as a reference and + // avoid copying the object. The same goes for any values that a compound + // function might depend on. Said another way, if there is a dependency + // chain of calculations that come before the one that we're currently + // interested in doing, then we're guaranteed to resolve the dependencies + // before we proceed with using any of those values. Although there is a + // cost to fetching data from the cache, the "resolved dependency" concept + // might be sufficiently convenient to make it worth looking past the extra + // cost. If these material laws are embedded within a finite element + // framework, then the added cost might not even be noticeable. + template + const double &Magnetoviscoelastic_Constitutive_Law::get_det_F() const + { + const std::string name("det_F"); + if (cache.stores_object_with_name(name) == false) + { + const double det_F = std::sqrt(determinant(get_C())); + AssertThrow(det_F > 0.0, + ExcMessage("Volumetric Jacobian must be positive.")); + cache.add_unique_copy(name, det_F); + } + + return cache.template get_object_with_name(name); + } + + + template + const SymmetricTensor<2, dim> & + Magnetoviscoelastic_Constitutive_Law::get_C_inv() const + { + const std::string name("C_inv"); + if (cache.stores_object_with_name(name) == false) + { + cache.add_unique_copy(name, invert(get_C())); + } + + return cache.template get_object_with_name>(name); + } + + + template + const double & + Magnetoviscoelastic_Constitutive_Law::get_log_det_F() const + { + const std::string name("log(det_F)"); + if (cache.stores_object_with_name(name) == false) + cache.add_unique_copy(name, std::log(get_det_F())); + + return cache.template get_object_with_name(name); + } + + + template + const double &Magnetoviscoelastic_Constitutive_Law::get_trace_C() const + { + const std::string name("trace(C)"); + if (cache.stores_object_with_name(name) == false) + cache.add_unique_copy(name, trace(get_C())); + + return cache.template get_object_with_name(name); + } + + + template + const Tensor<1, dim> & + Magnetoviscoelastic_Constitutive_Law::get_C_inv_dot_H() const + { + const std::string name("C_inv_dot_H"); + if (cache.stores_object_with_name(name) == false) + cache.add_unique_copy(name, get_C_inv() * get_H()); + + return cache.template get_object_with_name>(name); + } + + + template + const double & + Magnetoviscoelastic_Constitutive_Law::get_H_dot_C_inv_dot_H() const + { + const std::string name("H_dot_C_inv_dot_H"); + if (cache.stores_object_with_name(name) == false) + cache.add_unique_copy(name, get_H() * get_C_inv_dot_H()); + + return cache.template get_object_with_name(name); + } + + + template + const SymmetricTensor<4, dim> & + Magnetoviscoelastic_Constitutive_Law::get_dQ_t_dC( + const DiscreteTime &time) const + { + const std::string name("dQ_t_dC"); + if (cache.stores_object_with_name(name) == false) + { + const double delta_t = this->get_delta_t(time); + const double &det_F = get_det_F(); + + const SymmetricTensor<4, dim> dQ_t_dC = + (1.0 / (1.0 + delta_t / this->get_tau_v())) * + (delta_t / this->get_tau_v()) * + ((2.0 / dim) * std::pow(det_F, 2.0 / dim - 1.0) * + outer_product(get_C_inv(), get_ddet_F_dC()) + + std::pow(det_F, 2.0 / dim) * get_dC_inv_dC()); + + cache.add_unique_copy(name, dQ_t_dC); + } + + return cache.template get_object_with_name>(name); + } + + + template + const SymmetricTensor<4, dim> & + Magnetoviscoelastic_Constitutive_Law::get_dC_inv_dC() const + { + const std::string name("dC_inv_dC"); + if (cache.stores_object_with_name(name) == false) + { + const SymmetricTensor<2, dim> &C_inv = get_C_inv(); + SymmetricTensor<4, dim> dC_inv_dC; + + for (unsigned int A = 0; A < dim; ++A) + for (unsigned int B = A; B < dim; ++B) + for (unsigned int C = 0; C < dim; ++C) + for (unsigned int D = C; D < dim; ++D) + dC_inv_dC[A][B][C][D] -= // + 0.5 * (C_inv[A][C] * C_inv[B][D] // + + C_inv[A][D] * C_inv[B][C]); // + + cache.add_unique_copy(name, dC_inv_dC); + } + + return cache.template get_object_with_name>(name); + } + + + template + const SymmetricTensor<2, dim> & + Magnetoviscoelastic_Constitutive_Law::get_d_tr_C_dC() const + { + const std::string name("d_tr_C_dC"); + if (cache.stores_object_with_name(name) == false) + cache.add_unique_copy(name, + Physics::Elasticity::StandardTensors::I); + + return cache.template get_object_with_name>(name); + } + + + template + const SymmetricTensor<2, dim> & + Magnetoviscoelastic_Constitutive_Law::get_ddet_F_dC() const + { + const std::string name("ddet_F_dC"); + if (cache.stores_object_with_name(name) == false) + cache.add_unique_copy(name, 0.5 * get_det_F() * get_C_inv()); + + return cache.template get_object_with_name>(name); + } + + + template + const SymmetricTensor<2, dim> & + Magnetoviscoelastic_Constitutive_Law::get_dlog_det_F_dC() const + { + const std::string name("dlog_det_F_dC"); + if (cache.stores_object_with_name(name) == false) + cache.add_unique_copy(name, 0.5 * get_C_inv()); + + return cache.template get_object_with_name>(name); + } + + + template + const Tensor<1, dim> & + Magnetoviscoelastic_Constitutive_Law::get_dH_dot_C_inv_dot_H_dH() const + { + const std::string name("dH_dot_C_inv_dot_H_dH"); + if (cache.stores_object_with_name(name) == false) + cache.add_unique_copy(name, 2.0 * get_C_inv_dot_H()); + + return cache.template get_object_with_name>(name); + } + + + template + const SymmetricTensor<2, dim> & + Magnetoviscoelastic_Constitutive_Law::get_dH_dot_C_inv_dot_H_dC() const + { + const std::string name("dH_dot_C_inv_dot_H_dC"); + if (cache.stores_object_with_name(name) == false) + { + const Tensor<1, dim> C_inv_dot_H = get_C_inv_dot_H(); + cache.add_unique_copy( + name, -symmetrize(outer_product(C_inv_dot_H, C_inv_dot_H))); + } + + return cache.template get_object_with_name>(name); + } + + + template + const SymmetricTensor<4, dim> & + Magnetoviscoelastic_Constitutive_Law::get_d2log_det_F_dC_dC() const + { + const std::string name("d2log_det_F_dC_dC"); + if (cache.stores_object_with_name(name) == false) + cache.add_unique_copy(name, 0.5 * get_dC_inv_dC()); + + return cache.template get_object_with_name>(name); + } + + + template + const SymmetricTensor<4, dim> & + Magnetoviscoelastic_Constitutive_Law::get_d2det_F_dC_dC() const + { + const std::string name("d2det_F_dC_dC"); + if (cache.stores_object_with_name(name) == false) + cache.add_unique_copy(name, + 0.5 * + (outer_product(get_C_inv(), get_ddet_F_dC()) + + get_det_F() * get_dC_inv_dC())); + + return cache.template get_object_with_name>(name); + } + + + template + const SymmetricTensor<2, dim> & + Magnetoviscoelastic_Constitutive_Law::get_d2H_dot_C_inv_dot_H_dH_dH() + const + { + const std::string name("d2H_dot_C_inv_dot_H_dH_dH"); + if (cache.stores_object_with_name(name) == false) + cache.add_unique_copy(name, 2.0 * get_C_inv()); + + return cache.template get_object_with_name>(name); + } + + + template + const Tensor<3, dim> & + Magnetoviscoelastic_Constitutive_Law::get_d2H_dot_C_inv_dot_H_dC_dH() + const + { + const std::string name("d2H_dot_C_inv_dot_H_dC_dH"); + if (cache.stores_object_with_name(name) == false) + { + const Tensor<1, dim> & C_inv_dot_H = get_C_inv_dot_H(); + const SymmetricTensor<2, dim> &C_inv = get_C_inv(); + + Tensor<3, dim> d2H_dot_C_inv_dot_H_dC_dH; + for (unsigned int A = 0; A < dim; ++A) + for (unsigned int B = 0; B < dim; ++B) + for (unsigned int C = 0; C < dim; ++C) + d2H_dot_C_inv_dot_H_dC_dH[A][B][C] -= + C_inv[A][C] * C_inv_dot_H[B] + // + C_inv_dot_H[A] * C_inv[B][C]; // + + cache.add_unique_copy(name, d2H_dot_C_inv_dot_H_dC_dH); + } + + return cache.template get_object_with_name>(name); + } + + + template + const SymmetricTensor<4, dim> & + Magnetoviscoelastic_Constitutive_Law::get_d2H_dot_C_inv_dot_H_dC_dC() + const + { + const std::string name("d2H_dot_C_inv_dot_H_dC_dC"); + if (cache.stores_object_with_name(name) == false) + { + const Tensor<1, dim> & C_inv_dot_H = get_C_inv_dot_H(); + const SymmetricTensor<2, dim> &C_inv = get_C_inv(); + + SymmetricTensor<4, dim> d2H_dot_C_inv_dot_H_dC_dC; + for (unsigned int A = 0; A < dim; ++A) + for (unsigned int B = A; B < dim; ++B) + for (unsigned int C = 0; C < dim; ++C) + for (unsigned int D = C; D < dim; ++D) + d2H_dot_C_inv_dot_H_dC_dC[A][B][C][D] += + 0.5 * (C_inv_dot_H[A] * C_inv_dot_H[C] * C_inv[B][D] + + C_inv_dot_H[A] * C_inv_dot_H[D] * C_inv[B][C] + + C_inv_dot_H[B] * C_inv_dot_H[C] * C_inv[A][D] + + C_inv_dot_H[B] * C_inv_dot_H[D] * C_inv[A][C]); + + cache.add_unique_copy(name, d2H_dot_C_inv_dot_H_dC_dC); + } + + return cache.template get_object_with_name>(name); + } + + // @sect4{Rheological experiment parameters} + + // The @p RheologicalExperimentParameters class is used to drive the + // numerical experiments that are to be conducted on the coupled materials + // that we've implemented constitutive laws for. + class RheologicalExperimentParameters : public ParameterAcceptor + { + public: + RheologicalExperimentParameters(); + + // These are dimensions of the rheological specimen that is to be + // simulated. They, effectively, define the measurement point for our + // virtual experiment. + double sample_radius = 0.01; + double sample_height = 0.001; + + // The three steady-state loading parameters are respectively + // - the axial stretch, + // - the shear strain amplitude, and + // - the axial magnetic field strength. + double lambda_2 = 0.95; + double gamma_12 = 0.05; + double H_2 = 60.0e3; + + // Moreover, the parameters for the time-dependent rheological + // loading conditions are + // - the loading cycle frequency, + // - the number of load cycles, and + // - the number of discrete timesteps per cycle. + double frequency = 1.0 / (2.0 * numbers::PI); + unsigned int n_cycles = 5; + unsigned int n_steps_per_cycle = 2500; + + // We also declare some self-explanatory parameters related to output + // data generated for the experiments conducted with rate-dependent and + // rate-independent materials. + bool output_data_to_file = true; + std::string output_filename_rd = + "experimental_results-rate_dependent.csv"; + std::string output_filename_ri = + "experimental_results-rate_independent.csv"; + + // The next few functions compute time-related parameters for the + // experiment... + double start_time() const; + + double end_time() const; + + double delta_t() const; + + // ... while the following two prescribe the mechanical and magnetic + // loading at any given time... + Tensor<1, 3> get_H(const double time) const; + + Tensor<2, 3> get_F(const double time) const; + + // ... and this last one outputs the status of the experiment to the + // console. + bool print_status(const int step_number) const; + + bool initialized = false; + }; + + + + RheologicalExperimentParameters::RheologicalExperimentParameters() + : ParameterAcceptor("/Coupled Constitutive Laws/Rheological Experiment/") + { + add_parameter("Experimental sample radius", sample_radius); + add_parameter("Experimental sample radius", sample_height); + + add_parameter("Axial stretch", lambda_2); + add_parameter("Shear strain amplitude", gamma_12); + add_parameter("Axial magnetic field strength", H_2); + + add_parameter("Frequency", frequency); + add_parameter("Number of loading cycles", n_cycles); + add_parameter("Discretisation for each cycle", n_steps_per_cycle); + + add_parameter("Output experimental results to file", output_data_to_file); + add_parameter("Output file name (rate dependent constitutive law)", + output_filename_rd); + add_parameter("Output file name (rate independent constitutive law)", + output_filename_ri); + + parse_parameters_call_back.connect([&]() -> void { initialized = true; }); + } + + + double RheologicalExperimentParameters::start_time() const + { + return 0.0; + } + + + double RheologicalExperimentParameters::end_time() const + { + return n_cycles / frequency; + } + + + double RheologicalExperimentParameters::delta_t() const + { + return (end_time() - start_time()) / (n_steps_per_cycle * n_cycles); + } + + + bool + RheologicalExperimentParameters::print_status(const int step_number) const + { + return (step_number % (n_cycles * n_steps_per_cycle / 100)) == 0; + } + + // The applied magnetic field is always aligned with the axis of rotation + // of the rheometer's rotor. + Tensor<1, 3> RheologicalExperimentParameters::get_H(const double) const + { + return Tensor<1, 3>({0.0, 0.0, H_2}); + } + + // The applied deformation (gradient) is computed based on the geometry + // of the rheometer and the sample, the sampling point, and the experimental + // parameters. From the displacement profile documented in the introduction, + // the deformation gradient may be expressed in Cartesian coordinates as + // @f[ + // \mathbf{F} = \begin{bmatrix} + // \frac{\cos\left(\alpha\right)}{\sqrt{\lambda_{3}}} + // & -\frac{\sin\left(\alpha\right)}{\sqrt{\lambda_{3}}} + // & -\tau R \sqrt{\lambda_{3}} \sin\left(\Theta + \alpha\right) + // \\ \frac{\sin\left(\alpha\right)}{\sqrt{\lambda_{3}}} + // & \frac{\cos\left(\alpha\right)}{\sqrt{\lambda_{3}}} + // & -\tau R \sqrt{\lambda_{3}} \cos\left(\Theta + \alpha\right) + // \\ 0 & 0 & \lambda_{3} + // \end{bmatrix} + // @f] + Tensor<2, 3> RheologicalExperimentParameters::get_F(const double time) const + { + AssertThrow((sample_radius > 0.0 && sample_height > 0.0), + ExcMessage("Non-physical sample dimensions")); + AssertThrow(lambda_2 > 0.0, + ExcMessage("Non-physical applied axial stretch")); + + const double sqrt_lambda_2 = std::sqrt(lambda_2); + const double inv_sqrt_lambda_2 = 1.0 / sqrt_lambda_2; + + const double alpha_max = + std::atan(std::tan(gamma_12) * sample_height / + sample_radius); // Small strain approximation + const double A = sample_radius * alpha_max; + const double w = 2.0 * numbers::PI * frequency; // in rad /s + const double gamma_t = A * std::sin(w * time); + const double tau_t = + gamma_t / + (sample_radius * sample_height); // Torsion angle per unit length + const double alpha_t = tau_t * lambda_2 * sample_height; + + Tensor<2, 3> F; + F[0][0] = inv_sqrt_lambda_2 * std::cos(alpha_t); + F[0][1] = -inv_sqrt_lambda_2 * std::sin(alpha_t); + F[0][2] = -tau_t * sample_radius * sqrt_lambda_2 * std::sin(alpha_t); + F[1][0] = inv_sqrt_lambda_2 * std::sin(alpha_t); + F[1][1] = inv_sqrt_lambda_2 * std::cos(alpha_t); + F[1][2] = tau_t * sample_radius * sqrt_lambda_2 * std::cos(alpha_t); + F[2][0] = 0.0; + F[2][1] = 0.0; + F[2][2] = lambda_2; + + AssertThrow((F[0][0] > 0) && (F[1][1] > 0) && (F[2][2] > 0), + ExcMessage("Non-physical deformation gradient component.")); + AssertThrow(std::abs(determinant(F) - 1.0) < 1e-6, + ExcMessage("Volumetric Jacobian is not equal to unity.")); + + return F; + } + + // @sect4{Rheological experiment: Parallel plate rotational rheometer} + + // This is the function that will drive the numerical experiments. + template + void run_rheological_experiment( + const RheologicalExperimentParameters &experimental_parameters, + Coupled_Magnetomechanical_Constitutive_Law_Base + &material_hand_calculated, + Coupled_Magnetomechanical_Constitutive_Law_Base + & material_assisted_computation, + TimerOutput & timer, + const std::string filename) + { + // We can take the hand-implemented constitutive law and compare the + // results that we attain with it to those that we get using AD or SD. + // In this way, we can verify that they produce identical results (which + // indicates that either both implementations have a high probability of + // being correct, or that they're incorrect with identical flaws being + // present in both). Either way, it is a decent sanity check for the + // fully self-implemented variants and can certainly be used as a + // debugging strategy when differences between the results are + // detected). + const auto check_material_class_results = + []( + const Coupled_Magnetomechanical_Constitutive_Law_Base &to_verify, + const Coupled_Magnetomechanical_Constitutive_Law_Base &blessed, + const double tol = 1e-6) { + (void)to_verify; + (void)blessed; + (void)tol; + + Assert(std::abs(blessed.get_psi() - to_verify.get_psi()) < tol, + ExcMessage("No match for psi. Error: " + + Utilities::to_string(std::abs( + blessed.get_psi() - to_verify.get_psi())))); + + Assert((blessed.get_B() - to_verify.get_B()).norm() < tol, + ExcMessage("No match for B. Error: " + + Utilities::to_string( + (blessed.get_B() - to_verify.get_B()).norm()))); + Assert((blessed.get_S() - to_verify.get_S()).norm() < tol, + ExcMessage("No match for S. Error: " + + Utilities::to_string( + (blessed.get_S() - to_verify.get_S()).norm()))); + + Assert((blessed.get_DD() - to_verify.get_DD()).norm() < tol, + ExcMessage("No match for BB. Error: " + + Utilities::to_string( + (blessed.get_DD() - to_verify.get_DD()).norm()))); + Assert((blessed.get_PP() - to_verify.get_PP()).norm() < tol, + ExcMessage("No match for PP. Error: " + + Utilities::to_string( + (blessed.get_PP() - to_verify.get_PP()).norm()))); + Assert((blessed.get_HH() - to_verify.get_HH()).norm() < tol, + ExcMessage("No match for HH. Error: " + + Utilities::to_string( + (blessed.get_HH() - to_verify.get_HH()).norm()))); + }; + + // We'll be outputting the constitutive response of the material to file + // for post-processing, so here we declare a `stream` that will act as + // a buffer for this output. We'll use a simple CSV format for the + // outputted results. + std::ostringstream stream; + stream + << "Time;Axial magnetic field strength [A/m];Axial magnetic induction [T];Shear strain [%];Shear stress [Pa]\n"; + + // Using the DiscreteTime class, we iterate through each timestep using + // a fixed time step size. + for (DiscreteTime time(experimental_parameters.start_time(), + experimental_parameters.end_time() + + experimental_parameters.delta_t(), + experimental_parameters.delta_t()); + time.is_at_end() == false; + time.advance_time()) + { + if (experimental_parameters.print_status(time.get_step_number())) + std::cout << "Timestep = " << time.get_step_number() + << " @ time = " << time.get_current_time() << "s." + << std::endl; + + // We fetch and compute the loading to be applied to the material + // at this time step... + const Tensor<1, dim> H = + experimental_parameters.get_H(time.get_current_time()); + const Tensor<2, dim> F = + experimental_parameters.get_F(time.get_current_time()); + const SymmetricTensor<2, dim> C = + Physics::Elasticity::Kinematics::C(F); + + // ... then we update the state of the materials... + { + TimerOutput::Scope timer_section(timer, "Hand calculated"); + material_hand_calculated.update_internal_data(C, H, time); + material_hand_calculated.update_end_of_timestep(); + } + + { + TimerOutput::Scope timer_section(timer, "Assisted computation"); + material_assisted_computation.update_internal_data(C, H, time); + material_assisted_computation.update_end_of_timestep(); + } + + // ... and test for discrepencies between the two. + check_material_class_results(material_hand_calculated, + material_assisted_computation); + + if (experimental_parameters.output_data_to_file) + { + // The next thing that we will do is collect some results to + // post-process. All quantities are in the "current configuration" + // (rather than the "reference configuration", in which all + // quantities computed by the constitutive laws are framed). + const Tensor<1, dim> h = + Physics::Transformations::Covariant::push_forward(H, F); + const Tensor<1, dim> b = + Physics::Transformations::Piola::push_forward( + material_hand_calculated.get_B(), F); + const SymmetricTensor<2, dim> sigma = + Physics::Transformations::Piola::push_forward( + material_hand_calculated.get_S(), F); + stream << time.get_current_time() << ";" << h[2] << ";" << b[2] + << ";" << F[1][2] * 100.0 << ";" << sigma[1][2] << "\n"; + } + } + + // Finally, we output the strain-stress and magnetic loading history to + // file. + if (experimental_parameters.output_data_to_file) + { + std::ofstream output(filename); + output << stream.str(); + } + }; + + // @sect4{The CoupledConstitutiveLaws::run() function} + + // The purpose of this driver function is to read in all of the parameters + // from file and, based off of that, create a representative instance of + // each constitutive law and invoke the function that conducts a rheological + // experiment with it. + void run(int argc, char *argv[]) + { + using namespace dealii; + + constexpr unsigned int dim = 3; + + const ConstitutiveParameters constitutive_parameters; + const RheologicalExperimentParameters experimental_parameters; + + std::string parameter_file; + if (argc > 1) + parameter_file = argv[1]; + else + parameter_file = "parameters.prm"; + ParameterAcceptor::initialize(parameter_file, "used_parameters.prm"); + + // We start the actual work by configuring and running the experiment + // using our rate-independent constitutive law. The automatically + // differentiable number type is hard-coded here, but with some clever + // templating it is possible to select which framework to use at run time + // (e.g., as selected through the parameter file). We'll simultaneously + // perform the experiments with the counterpary material law that was + // fully implemented by hand, and check what it computes against our + // assisted implementation. + { + TimerOutput timer(std::cout, + TimerOutput::summary, + TimerOutput::wall_times); + std::cout + << "Coupled magnetoelastic constitutive law using automatic differentiation." + << std::endl; + + constexpr Differentiation::AD::NumberTypes ADTypeCode = + Differentiation::AD::NumberTypes::sacado_dfad_dfad; + + Magnetoelastic_Constitutive_Law material(constitutive_parameters); + Magnetoelastic_Constitutive_Law_AD material_ad( + constitutive_parameters); + + run_rheological_experiment(experimental_parameters, + material, + material_ad, + timer, + experimental_parameters.output_filename_ri); + + std::cout << "... all calculations are correct!" << std::endl; + } + + // Next we do the same for the rate-dependent constitutive law. + // The highest performance option is selected as default if SymEngine + // is set up to use the LLVM just-in-time compiler which (in conjunction + // with some aggressive compilation flags) produces the fastest code + // evaluation path of all of the available option. As a fall-back, the + // so called "lambda" optimizer (which only requires a C++11 compliant + // compiler) will be selected. At the same time, we'll ask the CAS to + // perform common subexpression elimination to minimize the number of + // intermediate calculations used during evaluation. + // We'll record how long it takes to execute the "initialization" step + // inside the constructor for the SD implementation, as this is where the + // abovementioned transformations occur. + { + TimerOutput timer(std::cout, + TimerOutput::summary, + TimerOutput::wall_times); + std::cout + << "Coupled magneto-viscoelastic constitutive law using symbolic differentiation." + << std::endl; + +#ifdef DEAL_II_SYMENGINE_WITH_LLVM + std::cout << "Using LLVM optimizer." << std::endl; + constexpr Differentiation::SD::OptimizerType optimizer_type = + Differentiation::SD::OptimizerType::llvm; + constexpr Differentiation::SD::OptimizationFlags optimization_flags = + Differentiation::SD::OptimizationFlags::optimize_all; +#else + std::cout << "Using lambda optimizer." << std::endl; + constexpr Differentiation::SD::OptimizerType optimizer_type = + Differentiation::SD::OptimizerType::lambda; + constexpr Differentiation::SD::OptimizationFlags optimization_flags = + Differentiation::SD::OptimizationFlags::optimize_cse; +#endif + + Magnetoviscoelastic_Constitutive_Law material( + constitutive_parameters); + + timer.enter_subsection("Initialize symbolic CL"); + Magnetoviscoelastic_Constitutive_Law_SD material_sd( + constitutive_parameters, optimizer_type, optimization_flags); + timer.leave_subsection(); + + run_rheological_experiment(experimental_parameters, + material, + material_sd, + timer, + experimental_parameters.output_filename_rd); + + std::cout << "... all calculations are correct!" << std::endl; + } + } + + } // namespace CoupledConstitutiveLaws + +} // namespace Step71 + + +// @sect3{The main() function} + +// The main function only calls the driver functions for the two sets of +// examples that are to be executed. +int main(int argc, char *argv[]) +{ + Step71::SimpleExample::run(); + Step71::CoupledConstitutiveLaws::run(argc, argv); + + return 0; +}