Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "KKSPhaseConcentrationDerivatives.h"
      11             : #include "MatrixTools.h"
      12             : 
      13             : registerMooseObject("PhaseFieldApp", KKSPhaseConcentrationDerivatives);
      14             : 
      15             : InputParameters
      16          62 : KKSPhaseConcentrationDerivatives::validParams()
      17             : {
      18          62 :   InputParameters params = DerivativeMaterialInterface<Material>::validParams();
      19          62 :   params.addClassDescription(
      20             :       "Computes the KKS phase concentration derivatives wrt global concentrations and order "
      21             :       "parameters, which are used in the chain rules in the KKS kernels. This class is intended to "
      22             :       "be used with KKSPhaseConcentrationMaterial.");
      23         124 :   params.addRequiredCoupledVar("global_cs", "The interpolated concentrations c, b, etc");
      24         124 :   params.addRequiredCoupledVar("eta", "Order parameter.");
      25         124 :   params.addRequiredParam<std::vector<MaterialPropertyName>>(
      26             :       "ci_names",
      27             :       "Phase concentrations. The order must match Fa, Fb, and global_cs, for example, c1, "
      28             :       "c2, b1, b2, etc");
      29         124 :   params.addRequiredParam<MaterialName>("fa_name", "Fa material object.");
      30         124 :   params.addRequiredParam<MaterialName>("fb_name", "Fb material object.");
      31         124 :   params.addParam<MaterialPropertyName>("h_name", "h", "Switching function h(eta).");
      32          62 :   return params;
      33           0 : }
      34             : 
      35          48 : KKSPhaseConcentrationDerivatives::KKSPhaseConcentrationDerivatives(
      36          48 :     const InputParameters & parameters)
      37             :   : DerivativeMaterialInterface<Material>(parameters),
      38          48 :     _num_c(coupledComponents("global_cs")),
      39          48 :     _c_names(coupledNames("global_cs")),
      40          48 :     _eta_name(getVar("eta", 0)->name()),
      41          96 :     _ci_names(getParam<std::vector<MaterialPropertyName>>("ci_names")),
      42          48 :     _prop_ci(_num_c * 2),
      43          48 :     _dcidb(_num_c),
      44          48 :     _dcideta(_num_c),
      45          48 :     _Fa_name(getParam<MaterialName>("fa_name")),
      46          48 :     _Fb_name(getParam<MaterialName>("fb_name")),
      47          48 :     _d2Fidcidbi(2),
      48          96 :     _prop_h(getMaterialProperty<Real>("h_name")),
      49          96 :     _prop_dh(getMaterialPropertyDerivative<Real>("h_name", _eta_name))
      50             : {
      51         144 :   for (const auto m : make_range(_num_c * 2))
      52          96 :     _prop_ci[m] = &getMaterialPropertyByName<Real>(_ci_names[m]);
      53             : 
      54          96 :   for (const auto m : make_range(_num_c))
      55             :   {
      56          48 :     _dcideta[m].resize(2);
      57          48 :     _dcidb[m].resize(2);
      58         144 :     for (const auto n : make_range(2))
      59             :     {
      60             :       // Derivative of phase concentration wrt eta. In _dcideta[m][n], m is the species index of
      61             :       // ci, n is the phase index of ci
      62          96 :       _dcideta[m][n] = &declarePropertyDerivative<Real>(_ci_names[m * 2 + n], _eta_name);
      63          96 :       _dcidb[m][n].resize(_num_c);
      64             : 
      65             :       // Derivative of phase concentration wrt global concentration. In _dcidb[m][n][l], m is the
      66             :       //  species index of ci, n is the phase index of ci, l is the species index of b
      67         192 :       for (const auto l : make_range(_num_c))
      68          96 :         _dcidb[m][n][l] = &declarePropertyDerivative<Real>(_ci_names[m * 2 + n], _c_names[l]);
      69             :     }
      70             :   }
      71             : 
      72             :   /** Second derivative of free energy wrt phase concentrations for use in this material. In
      73             :       _d2Fidcidbi[m][n][l], m is phase index of Fi, n is the species index of ci, l is the species
      74             :       index of bi.
      75             :   */
      76         144 :   for (const auto m : make_range(2))
      77             :   {
      78          96 :     _d2Fidcidbi[m].resize(_num_c);
      79         192 :     for (const auto n : make_range(_num_c))
      80             :     {
      81          96 :       _d2Fidcidbi[m][n].resize(_num_c);
      82         192 :       for (const auto l : make_range(_num_c))
      83             :       {
      84          96 :         if (m == 0)
      85          48 :           _d2Fidcidbi[0][n][l] =
      86          96 :               &getMaterialPropertyDerivative<Real>(_Fa_name, _ci_names[n * 2], _ci_names[l * 2]);
      87             : 
      88             :         else
      89          48 :           _d2Fidcidbi[1][n][l] = &getMaterialPropertyDerivative<Real>(
      90          48 :               _Fb_name, _ci_names[n * 2 + 1], _ci_names[l * 2 + 1]);
      91             :       }
      92             :     }
      93             :   }
      94          48 : }
      95             : 
      96             : void
      97     1501200 : KKSPhaseConcentrationDerivatives::computeQpProperties()
      98             : {
      99             :   // declare Jacobian matrix A
     100     1501200 :   Eigen::MatrixXd A(_num_c * 2, _num_c * 2);
     101             : 
     102             :   A.setZero();
     103             : 
     104             :   // fill in the non-zero elements in A
     105     3002400 :   for (const auto m : make_range(_num_c))
     106             :   {
     107     3002400 :     for (const auto n : make_range(_num_c))
     108             :     {
     109             :       // equal chemical potential derivative equations
     110     1501200 :       A(m * 2, n * 2) = (*_d2Fidcidbi[0][m][n])[_qp];
     111     1501200 :       A(m * 2, n * 2 + 1) = -(*_d2Fidcidbi[1][m][n])[_qp];
     112             :     }
     113             : 
     114             :     // concentration conservation derivative equations
     115     1501200 :     A(m * 2 + 1, m * 2) = 1 - _prop_h[_qp];
     116     1501200 :     A(m * 2 + 1, m * 2 + 1) = _prop_h[_qp];
     117             :   }
     118             : 
     119     1501200 :   A = A.inverse();
     120             : 
     121             :   // solve linear system of constraint derivatives wrt b for computing dcidb loop through
     122             :   // derivatives wrt the ith component; they have the same A, but different k_c
     123     3002400 :   for (const auto i : make_range(_num_c))
     124             :   {
     125     1501200 :     std::vector<Real> k_c(_num_c * 2);
     126     1501200 :     std::vector<Real> x_c(_num_c * 2);
     127             : 
     128             :     // assign the non-zero elements in k_c
     129     1501200 :     k_c[i * 2 + 1] = 1;
     130             : 
     131             :     // compute x_c
     132     4503600 :     for (const auto m : make_range(_num_c * 2))
     133             :     {
     134     9007200 :       for (const auto n : make_range(_num_c * 2))
     135     6004800 :         x_c[m] += A(m, n) * k_c[n];
     136             :     }
     137             : 
     138             :     // assign the values in x_c to _dcidb
     139     3002400 :     for (const auto m : make_range(_num_c))
     140             :     {
     141     4503600 :       for (const auto n : make_range(2))
     142     3002400 :         (*_dcidb[m][n][i])[_qp] = x_c[m * 2 + n];
     143             :     }
     144     1501200 :   }
     145             : 
     146             :   // solve linear system of constraint derivatives wrt eta for computing dcideta use the same
     147             :   // linear matrix as computing dcidb
     148     1501200 :   std::vector<Real> k_eta(_num_c * 2);
     149     1501200 :   std::vector<Real> x_eta(_num_c * 2);
     150             : 
     151             :   // fill in k_eta
     152     3002400 :   for (const auto m : make_range(_num_c))
     153             :   {
     154     1501200 :     k_eta[m * 2] = 0;
     155     1501200 :     k_eta[m * 2 + 1] = _prop_dh[_qp] * ((*_prop_ci[m * 2])[_qp] - (*_prop_ci[m * 2 + 1])[_qp]);
     156             :   }
     157             : 
     158             :   // compute x_eta
     159     4503600 :   for (const auto m : make_range(_num_c * 2))
     160             :   {
     161     9007200 :     for (const auto n : make_range(_num_c * 2))
     162     6004800 :       x_eta[m] += A(m, n) * k_eta[n];
     163             :   }
     164             : 
     165             :   // assign the values in x_eta to _dcideta
     166     3002400 :   for (const auto m : make_range(_num_c))
     167             :   {
     168     4503600 :     for (const auto n : make_range(2))
     169     3002400 :       (*_dcideta[m][n])[_qp] = x_eta[m * 2 + n];
     170             :   }
     171     1501200 : }