Merge pull request #19981 from xueyang94/next

A revised KKS solving method

modules/phase_field/doc/content/source/kernels/NestedKKSACBulkC.md

@@ -0,0 +1,65 @@
+# NestedKKSACBulkC
+
+!syntax description /Kernels/NestedKKSACBulkC
+
+Kim-Kim-Suzuki (KKS) nested solve kernel (1 of 3). An Allen-Cahn kernel for the terms with a direct composition dependence. This kernel can be used for one or multiple species.
+
+### Residual
+
+If a model has N species:
+
+\begin{equation}
+\begin{aligned}
+R&=&-\frac{dh}{d\eta}\left[-\sum_{c=1}^N\frac{dF_a}{dc_a}(c_a-c_b)\right]\\
+&=&\frac{dh}{d\eta}\left[\sum_{c=1}^N\frac{dF_a}{dc_a}(c_a-c_b)\right]
+\end{aligned}
+\end{equation}
+
+### Jacobian
+
+#### On-diagonal
+
+We are looking for the $\frac \partial{\partial u_j}$ derivative of $R$, where
+$u\equiv\eta$. We need to apply the chain rule and will only keep terms
+with the $\frac{\partial \eta}{\partial u_j}\frac{\partial}{\partial \eta} = \phi_j \frac{\partial}{\partial \eta}$ derivative.
+
+\begin{equation}
+\begin{aligned}
+J &=& \phi_j\frac{\partial}{\partial\eta}\left[\frac{dh}{d\eta}\sum_{c=1}^N\left(\frac{\partial F_a}{\partial c_a}(c_a-c_b)\right)\right]    \\
+&=&\phi_j\left[\frac{d^2h}{d\eta^2}\sum_{c=1}^N\left(\frac{\partial F_a}{\partial c_a}(c_a-c_b)\right) + \frac{dh}{d\eta}\frac{\partial}{\partial\eta}\left(\sum_{c=1}^N\frac{\partial F_a}{\partial c_a}(c_a-c_b)\right)  \right] \\
+&=&\phi_j\left[\frac{dh}{d\eta^2}\sum_{c=1}^N\left(\frac{\partial F_a}{\partial c_a}(c_a-c_b)\right) + \frac{dh}{d\eta}\sum_{c=1}^N\left(\sum_{b=1}^N(\frac{\partial^2F_a}{\partial c_a \partial b_a}\frac{\partial b_a}{\partial\eta})(c_a-b_a) + \frac{\partial F_a}{\partial c_a}(\frac{\partial c_a}{\partial\eta} - \frac{\partial c_b}{\partial\eta})\right)\right]
+\end{aligned}
+\end{equation}
+
+#### Off-diagonal
+
+Since the partial derivative of phase concentrations $c_a$ w.r.t global concentration $c$ is hidden when computing the $\frac \partial{\partial u_j}$ derivative of $R$, where $u\equiv c$, the derivative w.r.t $c$ must be calculated separetely from any other variable dependence. We need to
+apply the chain rule and will again only keep terms with the
+$\frac{\partial c}{\partial u_j}\frac{\partial}{\partial c} = \phi_j \frac{\partial}{\partial c}$
+derivative.
+
+\begin{equation}
+\begin{aligned}
+J &=& \phi_j\frac{dh}{d\eta}\frac{\partial}{\partial c}\left[\sum_{c=1}^N\left(\frac{\partial F_a}{\partial c_a}(c_a-c_b)\right)\right] \\
+&=& \phi_j\frac{dh}{d\eta}\sum_{c=1}^N\left(\sum_{b=1}^N(\frac{\partial^2F_a}{\partial c_a\partial b_a}\frac{\partial b_a}{\partial c})(c_a-c_b) + \frac{\partial F_a}{\partial c_a}(\frac{\partial c_a}{\partial c} - \frac{\partial c_b}{\partial c})\right)
+\end{aligned}
+\end{equation}
+
+If $F_a$ contains any other *explicit* variables, for example temperature $T$:
+
+\begin{equation}
+\begin{aligned}
+J &=& \phi_j \frac{\partial}{\partial T} \left( \frac{dh}{d\eta}\frac{dF_a}{dc_a}(c_a-c_b) \right)\\
+&=& \frac{dh}{d\eta} \phi_j  \frac{\partial}{\partial T}\left(\frac{d F_a}{d c_a}\right) (c_a - c_b)\\
+\end{aligned}
+\end{equation}
+
+The off-diagonal Jacobian contributions are multiplied by the Allen-Cahn
+mobility $L$ at each point for consistency with the other terms in the Allen-Cahn
+equation.
+
+!syntax parameters /Kernels/NestedKKSACBulkC
+
+!syntax inputs /Kernels/NestedKKSACBulkC
+
+!syntax children /Kernels/NestedKKSACBulkC

modules/phase_field/doc/content/source/kernels/NestedKKSACBulkF.md

@@ -0,0 +1,54 @@
+# NestedKKSACBulkF
+
+!syntax description /Kernels/NestedKKSACBulkF
+
+Kim-Kim-Suzuki (KKS) nested solve kernel (2 of 3). An Allen-Cahn kernel for the terms without a direct composition dependence.
+
+### Residual
+
+\begin{equation}
+R = -\frac{dh}{d\eta}(F_a-F_b) + w\frac{dg}{d\eta}
+\end{equation}
+
+### Jacobian
+
+#### On-diagonal
+
+We are looking for the $\frac \partial{\partial u_j}$ derivative of $R$, where
+$u\equiv\eta$. We need to apply the chain rule and will only keep terms
+with the $\frac{\partial \eta}{\partial u_j}\frac{\partial}{\partial \eta}=\phi_j \frac{\partial}{\partial\eta}$
+derivative.
+
+\begin{equation}
+\begin{aligned}
+J &=& -\phi_j\frac{\partial}{\partial\eta}\left(\frac{dh}{d\eta}(F_a-F_b)\right) + w\phi_j\frac{\partial}{\partial\eta}\frac{dg}{dh}    \\
+&=& -\phi_j\left(\frac{d^2h}{d\eta^2}(F_a-F_b) + \frac{dh}{d\eta}\sum_{c=1}^N(\frac{\partial F_a}{\partial c_a}\frac{\partial c_a}{\partial\eta} - \frac{\partial F_b}{\partial c_b}\frac{\partial c_b}{\partial\eta})\right) + w\phi_j\frac{d^2g}{d\eta^2}
+\end{aligned}
+\end{equation}
+
+#### Off-diagonal
+
+Since the partial derivative of phase concentrations $c_a$ w.r.t global concentration $c$ is hidden when computing the $\frac \partial{\partial u_j}$ derivative of $R$, where $u\equiv c$, the derivative w.r.t $c$ must be calculated separetely from any other variable dependence. We need to apply the chain rule and will again only keep terms
+with the $\frac{\partial c}{\partial u_j}\frac{\partial}{\partial c}=\phi_j \frac{\partial}{\partial c}$
+derivative.
+
+\begin{equation}
+\begin{aligned}
+J &=& \phi_j\frac{\partial}{\partial c}\left(-\frac{dh}{d\eta}(F_a-F_b) + w\frac{dg}{d\eta}\right)   \\
+&=& -\phi_j\frac{dh}{d\eta}\sum_{c=1}^N(\frac{\partial F_a}{\partial c_a}\frac{\partial c_a}{\partial c} - \frac{\partial F_b}{\partial c_b}\frac{\partial c_b}{\partial c})
+\end{aligned}
+\end{equation}
+
+If $F_a$ and $F_b$ contain any other *explicit* variables, for example temperature $T$:
+
+\begin{equation}
+J = -\phi_j\frac{dh}{d\eta}\left( \frac{\partial F_a}{\partial T} - \frac{\partial F_b}{\partial T}\right)
+\end{equation}
+
+The off-diagonal Jacobian contribution is multiplied by the Allen-Cahn mobility $L$ at each point for consistency with the other terms in the Allen-Cahn equation.
+
+!syntax parameters /Kernels/NestedKKSACBulkF
+
+!syntax inputs /Kernels/NestedKKSACBulkF
+
+!syntax children /Kernels/NestedKKSACBulkF

modules/phase_field/doc/content/source/kernels/NestedKKSSplitCHCRes.md

@@ -0,0 +1,56 @@
+# NestedKKSSplitCHCRes
+
+!syntax description /Kernels/NestedKKSSplitCHCRes
+
+Kim-Kim-Suzuki (KKS) nested solve kernel (3 of 3) for two-phase or multiphase models. A kernel for the split Cahn-Hilliard term. This kernel operates on the global concentration $c_i$ as the non-linear variable.
+
+## Residual
+
+\begin{equation}
+R = \frac{\partial F_1}{\partial c_{i,1}} - \mu,
+\end{equation}
+
+where $c_{i,1}$ is the phase concentration of species $i$ in the first phase.
+
+### Jacobian
+
+#### On-Diagonal
+
+We need to apply the chain rule and will only keep terms
+with the $\frac{\partial c_i}{\partial u_j}\frac{\partial}{\partial c_i}=\phi_j \frac{\partial}{\partial c_i}$
+derivative. If a system has C components:
+
+\begin{equation}
+\begin{aligned}
+J &=& \phi_j\frac{\partial}{\partial c_i}(\frac{\partial F_1}{\partial c_{i,1}} - \mu)   \\
+&=& \phi_j\sum_{j=1}^C(\frac{\partial^2F_1}{\partial c_{i,1}\partial c_{j,1}}\frac{\partial c_{j,1}}{\partial c_i})
+\end{aligned}
+\end{equation}
+
+#### Off-diagonal
+
+Since the partial derivative of phase concentrations $c_{i,1}$ w.r.t phase parameter $eta_p$ is hidden when computing the $\frac \partial{\partial u_j}$ derivative of $R$, where $u\equiv \eta_p$, the derivative w.r.t $\eta_p$ must be calculated separately from any other variable dependence. We need to apply the chain rule and will again only keep terms
+with the $\frac{\partial \eta_p}{\partial u_j}\frac{\partial}{\partial \eta_p}=\phi_j \frac{\partial}{\partial \eta_p}$ derivative.
+
+\begin{equation}
+\begin{aligned}
+J &=& \phi_j\frac{\partial}{\partial\eta_p}(\frac{\partial F_1}{\partial c_{i,1}} - \mu)   \\
+&=& \phi_j\sum_{j=1}^C(\frac{\partial^2 F_1}{\partial c_{i,1}\partial c{j,1}}\frac{\partial c_{j,1}}{\partial\eta_p})
+\end{aligned}
+\end{equation}
+
+
+If $F_1$ contains any other *explicit* variables, for example temperature $T$:
+
+\begin{equation}
+\begin{aligned}
+J &=& \phi_j\frac{\partial}{\partial T}(\frac{\partial F_1}{\partial c_{i,1}} - \mu)   \\
+&=& \phi_j\frac{\partial^2 F_1}{\partial c_{i,1}\partial T }
+\end{aligned}
+\end{equation}
+
+!syntax parameters /Kernels/NestedKKSSplitCHCRes
+
+!syntax inputs /Kernels/NestedKKSSplitCHCRes
+
+!syntax children /Kernels/NestedKKSSplitCHCRes

modules/phase_field/doc/content/source/materials/KKSPhaseConcentrationDerivatives.md

@@ -0,0 +1,17 @@
+# KKSPhaseConcentrationDerivatives
+
+Kim-Kim-Suzuki (KKS) nested solve material (part 2 of 2). KKSPhaseConcentrationDerivatives computes the partial derivatives of phase concentrations w.r.t global concentrations and phase parameters, for example, $\frac{\partial c_a}{\partial c}$, where $c$ is the global concentration and $c_a$ is the phase concentration of species $c$ in phase $a$. Another example is $\frac{\partial c_b}{\eta}$ where $\eta$ is one of the two phase parameters in the model and $c_b$ is the phase concentration of species $c$ in phase $b$. These partial derivatives are used in KKS kernels as chain rule terms in the residual and Jacobian.
+
+## Example input:
+
+!listing kks_example_nested.i block=Materials
+
+## Class Description
+
+!syntax description /Materials/KKSPhaseConcentrationDerivatives
+
+!syntax parameters /Materials/KKSPhaseConcentrationDerivatives
+
+!syntax inputs /Materials/KKSPhaseConcentrationDerivatives
+
+!syntax children /Materials/KKSPhaseConcentrationDerivatives

modules/phase_field/doc/content/source/materials/KKSPhaseConcentrationMaterial.md

@@ -0,0 +1,29 @@
+# KKSPhaseConcentrationMaterial
+
+Kim-Kim-Suzuki (KKS) nested solve material (part 1 of 2). KKSPhaseConcentrationMaterial implements a nested Newton iteration to solve the KKS constraint equations for the phase concentrations $c_a$ and $c_b$ as material properties (instead of non-linear variables as in the traditional solve in MOOSE). The constraint equations are the mass conservation equation for each global concentration ($c$):
+
+\begin{equation}
+c=\left(1-h(\eta)\right)c_a + h(\eta)c_b,
+\end{equation}
+
+and the pointwise equality of the phase chemical potentials:
+
+\begin{equation}
+\frac{\partial f_a}{\partial c_a} = \frac{\partial f_b}{\partial c_b}.
+\end{equation}
+
+The parameters Fa_material and Fb_material must have *compute=false*. This material also passes the phase free energies and their partial derivatives w.r.t phase concentrations to the KKS kernels (NestKKSACBulkC, NestKKSACBulkF, NestKKSSplitCHCRes).
+
+## Example input:
+
+!listing kks_example_nested.i block=Materials
+
+## Class Description
+
+!syntax description /Materials/KKSPhaseConcentrationMaterial
+
+!syntax parameters /Materials/KKSPhaseConcentrationMaterial
+
+!syntax inputs /Materials/KKSPhaseConcentrationMaterial
+
+!syntax children /Materials/KKSPhaseConcentrationMaterial

modules/phase_field/include/kernels/NestedKKSACBulkC.h

@@ -0,0 +1,61 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "KKSACBulkBase.h"
+
+/**
+ * KKSACBulkBase child class for the phase concentration difference term
+ * \f$ \frac{dh}{d\eta}\frac{\partial F_a}{\partial c_a}(c_a-c_b) \f$
+ * in the the Allen-Cahn bulk residual.
+ */
+class NestedKKSACBulkC : public KKSACBulkBase
+{
+public:
+  static InputParameters validParams();
+
+  NestedKKSACBulkC(const InputParameters & parameters);
+
+protected:
+  virtual Real computeDFDOP(PFFunctionType type);
+  virtual Real computeQpOffDiagJacobian(unsigned int jvar);
+
+  /// Global concentrations
+  std::vector<VariableName> _c_names;
+  const JvarMap & _c_map;
+
+  /// Number of global concentrations
+  const unsigned int _num_c;
+
+  /// Phase concentrations
+  const std::vector<MaterialPropertyName> _ci_names;
+
+  /// Phase concentration properties
+  std::vector<std::vector<const MaterialProperty<Real> *>> _prop_ci;
+
+  /// Derivative of phase concentrations wrt eta \f$ \frac d{d{eta}} c_i \f$
+  std::vector<std::vector<const MaterialProperty<Real> *>> _dcideta;
+
+  /// Derivative of phase concentrations wrt global concentrations \f$ \frac d{db} c_i \f$
+  std::vector<std::vector<std::vector<const MaterialProperty<Real> *>>> _dcidb;
+
+  /// Free energy of phase a
+  const MaterialPropertyName _Fa_name;
+
+  /// Derivative of the free energy function \f$ \frac d{dc_a} F_a \f$
+  std::vector<const MaterialProperty<Real> *> _dFadca;
+
+  /// Second derivative of the free energy function \f$ \frac {d^2}{dc_a db_a} F_a \f$
+  std::vector<std::vector<const MaterialProperty<Real> *>> _d2Fadcadba;
+
+  /// Mixed partial derivatives of the free energy function wrt c and
+  /// any other coupled variables \f$ \frac {d^2}{dc_a dq} F_a \f$
+  std::vector<std::vector<const MaterialProperty<Real> *>> _d2Fadcadarg;
+};

modules/phase_field/include/kernels/NestedKKSACBulkF.h

@@ -0,0 +1,76 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "KKSACBulkBase.h"
+
+/**
+ * KKSACBulkBase child class for the free energy difference term
+ * \f$ -\frac{dh}{d\eta}(F_a-F_b)+w\frac{dg}{d\eta} \f$
+ * in the the Allen-Cahn bulk residual.
+ */
+class NestedKKSACBulkF : public KKSACBulkBase
+{
+public:
+  static InputParameters validParams();
+
+  NestedKKSACBulkF(const InputParameters & parameters);
+
+protected:
+  virtual Real computeDFDOP(PFFunctionType type);
+  virtual Real computeQpOffDiagJacobian(unsigned int jvar);
+
+  const JvarMap & _c_map;
+
+  /// Number of global concentrations
+  const unsigned int _num_c;
+
+  /// Global concentrations
+  const std::vector<VariableName> _c_names;
+
+  /// Phase concentrations
+  const std::vector<MaterialPropertyName> _ci_names;
+
+  /// Free energy of phase a
+  const MaterialPropertyName _Fa_name;
+
+  /// Derivative of the free energy function \f$ \frac d{dc_a} F_a \f$
+  std::vector<const MaterialProperty<Real> *> _dFadca;
+
+  /// Free energy of phase b
+  const MaterialPropertyName _Fb_name;
+
+  /// Derivative of the free energy function \f$ \frac d{dc_b} F_b \f$
+  std::vector<const MaterialProperty<Real> *> _dFbdcb;
+
+  /// Derivative of barrier function g
+  const MaterialProperty<Real> & _prop_dg;
+
+  /// Second derivative of barrier function g
+  const MaterialProperty<Real> & _prop_d2g;
+
+  /// Double well height parameter
+  const Real _w;
+
+  /// Free energy properties
+  std::vector<const MaterialProperty<Real> *> _prop_Fi;
+
+  /// Derivative of phase concentration wrt eta \f$ \frac d{d{eta}} c_i \f$
+  std::vector<std::vector<const MaterialProperty<Real> *>> _dcideta;
+
+  /// Derivative of phase concentration wrt global concentration \f$ \frac d{b} c_i \f$
+  std::vector<std::vector<std::vector<const MaterialProperty<Real> *>>> _dcidb;
+
+  /// Partial derivative of the free energy function Fa wrt coupled variables \f$ \frac d{dq} F_a \f$
+  std::vector<const MaterialProperty<Real> *> _dFadarg;
+
+  /// Partial derivative of the free energy function Fb wrt coupled variables \f$ \frac d{dq} F_b \f$
+  std::vector<const MaterialProperty<Real> *> _dFbdarg;
+};