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 "NestedKKSMultiACBulkC.h"
      11             : 
      12             : registerMooseObject("PhaseFieldApp", NestedKKSMultiACBulkC);
      13             : 
      14             : InputParameters
      15          90 : NestedKKSMultiACBulkC::validParams()
      16             : {
      17          90 :   InputParameters params = KKSMultiACBulkBase::validParams();
      18          90 :   params.addClassDescription("Multi-phase KKS model kernel (part 2 of 2) for the Bulk Allen-Cahn. "
      19             :                              "This includes all terms dependent on chemical potential.");
      20         180 :   params.addRequiredCoupledVar("global_cs", "Global concentrations, for example, c, b.");
      21         180 :   params.addRequiredCoupledVar("all_etas", "Order parameters.");
      22         180 :   params.addRequiredParam<std::vector<MaterialPropertyName>>(
      23             :       "ci_names",
      24             :       "Phase concentrations. They must have the same order as Fj_names and global_cs, for "
      25             :       "example, c1, c2, b1, b2.");
      26          90 :   return params;
      27           0 : }
      28             : 
      29          48 : NestedKKSMultiACBulkC::NestedKKSMultiACBulkC(const InputParameters & parameters)
      30             :   : KKSMultiACBulkBase(parameters),
      31          48 :     _c_names(coupledNames("global_cs")),
      32          48 :     _c_map(getParameterJvarMap("global_cs")),
      33          48 :     _num_c(coupledComponents("global_cs")),
      34          48 :     _eta_names(coupledNames("all_etas")),
      35          48 :     _eta_map(getParameterJvarMap("all_etas")),
      36          48 :     _k(-1),
      37          96 :     _ci_names(getParam<std::vector<MaterialPropertyName>>("ci_names")),
      38          48 :     _ci_name_matrix(_num_c),
      39          48 :     _prop_ci(_num_c),
      40          48 :     _dcidetaj(_num_c),
      41          48 :     _dcidb(_num_c),
      42          48 :     _prop_d2hjdetaidetap(_num_j),
      43          48 :     _dF1dc1(_num_c),
      44          48 :     _d2F1dc1db1(_num_c),
      45          96 :     _d2F1dc1darg(_n_args)
      46             : {
      47         192 :   for (unsigned int i = 0; i < _num_j; ++i)
      48             :   {
      49             :     // get order parameter names and variable indices
      50         288 :     _eta_names[i] = getVar("all_etas", i)->name();
      51             : 
      52             :     // Set _k to the position of the nonlinear variable in the list of etaj's
      53         144 :     if (coupled("all_etas", i) == _var.number())
      54          48 :       _k = i;
      55             :   }
      56             : 
      57             :   // initialize _ci_name_matrix and _prop_ci for easy reference
      58          96 :   for (unsigned int m = 0; m < _num_c; ++m)
      59             :   {
      60          48 :     _ci_name_matrix[m].resize(_num_j);
      61          48 :     _prop_ci[m].resize(_num_j);
      62             : 
      63         192 :     for (const auto n : make_range(_num_j))
      64             :     {
      65         144 :       _ci_name_matrix[m][n] = _ci_names[m * _num_j + n];
      66         144 :       _prop_ci[m][n] = &getMaterialPropertyByName<Real>(_ci_name_matrix[m][n]);
      67             :     }
      68             :   }
      69             : 
      70             :   // _dcidetaj and _dcidb are computed in KKSPhaseConcentrationMultiPhaseDerivatives
      71          96 :   for (const auto m : make_range(_num_c))
      72             :   {
      73          48 :     _dcidetaj[m].resize(_num_j);
      74          48 :     _dcidb[m].resize(_num_j);
      75             : 
      76         192 :     for (const auto n : make_range(_num_j))
      77             :     {
      78         144 :       _dcidetaj[m][n].resize(_num_j);
      79         144 :       _dcidb[m][n].resize(_num_c);
      80             : 
      81         576 :       for (const auto l : make_range(_num_j))
      82         432 :         _dcidetaj[m][n][l] =
      83         864 :             &getMaterialPropertyDerivative<Real>(_ci_names[n + m * _num_j], _eta_names[l]);
      84             : 
      85         288 :       for (const auto l : make_range(_num_c))
      86         144 :         _dcidb[m][n][l] =
      87         288 :             &getMaterialPropertyDerivative<Real>(_ci_names[n + m * _num_j], _c_names[l]);
      88             :     }
      89             :   }
      90             : 
      91         192 :   for (const auto m : make_range(_num_j))
      92             :   {
      93         144 :     _prop_d2hjdetaidetap[m].resize(_num_j);
      94             : 
      95         576 :     for (const auto n : make_range(_num_j))
      96         432 :       _prop_d2hjdetaidetap[m][n] =
      97         864 :           &getMaterialPropertyDerivative<Real>(_hj_names[m], _eta_names[_k], _eta_names[n]);
      98             :   }
      99             : 
     100             :   // _dF1dc1 and _d2F1dc1db1 are computed in KKSPhaseConcentrationMultiPhaseMaterial
     101          96 :   for (const auto m : make_range(_num_c))
     102             :   {
     103          48 :     _dF1dc1[m] = &getMaterialPropertyDerivative<Real>(_Fj_names[0], _ci_names[m * _num_j]);
     104          48 :     _d2F1dc1db1[m].resize(_num_c);
     105             : 
     106          96 :     for (const auto n : make_range(_num_c))
     107          48 :       _d2F1dc1db1[m][n] = &getMaterialPropertyDerivative<Real>(
     108          48 :           _Fj_names[0], _ci_name_matrix[m][0], _ci_name_matrix[n][0]);
     109             :   }
     110             : 
     111             :   // _d2F1dc1darg are computed in KKSPhaseConcentrationMultiPhaseMaterial
     112          96 :   for (const auto m : make_range(_num_c))
     113             :   {
     114          48 :     _d2F1dc1darg[m].resize(_n_args);
     115             : 
     116         432 :     for (const auto n : make_range(_n_args))
     117         384 :       _d2F1dc1darg[m][n] =
     118         384 :           &getMaterialPropertyDerivative<Real>(_Fj_names[0], _ci_name_matrix[m][0], n);
     119             :   }
     120          48 : }
     121             : 
     122             : Real
     123    26500800 : NestedKKSMultiACBulkC::computeDFDOP(PFFunctionType type)
     124             : {
     125             :   Real sum = 0.0;
     126             : 
     127    26500800 :   switch (type)
     128             :   {
     129    24811200 :     case Residual:
     130             : 
     131    49622400 :       for (const auto m : make_range(_num_c))
     132             :       {
     133             :         Real sum1 = 0.0;
     134             : 
     135    99244800 :         for (const auto n : make_range(_num_j))
     136    74433600 :           sum1 += (*_prop_dhjdetai[n])[_qp] * (*_prop_ci[m][n])[_qp];
     137             : 
     138    24811200 :         sum += (*_dF1dc1[m])[_qp] * sum1;
     139             :       }
     140             : 
     141    24811200 :       return -sum;
     142             : 
     143     1689600 :     case Jacobian:
     144             :       /** For when this kernel is used in the Lagrange multiplier equation
     145             :           In that case the Lagrange multiplier is the nonlinear variable */
     146     1689600 :       if (_etai_var != _var.number())
     147             :         return 0.0;
     148             : 
     149     3379200 :       for (const auto m : make_range(_num_c))
     150             :       {
     151             :         Real sum1 = 0.0;
     152             :         Real sum2 = 0.0;
     153             :         Real sum3 = 0.0;
     154             : 
     155     3379200 :         for (const auto n : make_range(_num_c))
     156     1689600 :           sum1 += (*_d2F1dc1db1[m][n])[_qp] * (*_dcidetaj[n][0][_k])[_qp];
     157             : 
     158     6758400 :         for (const auto l : make_range(_num_j))
     159             :         {
     160     5068800 :           sum2 += (*_prop_dhjdetai[l])[_qp] * (*_prop_ci[m][l])[_qp];
     161             : 
     162     5068800 :           sum3 += (*_prop_d2hjdetaidetap[l][_k])[_qp] * (*_prop_ci[m][l])[_qp] +
     163     5068800 :                   (*_prop_dhjdetai[l])[_qp] * (*_dcidetaj[m][l][_k])[_qp];
     164             :         }
     165             : 
     166     1689600 :         sum += sum1 * sum2 + (*_dF1dc1[m])[_qp] * sum3;
     167             :       }
     168             : 
     169     1689600 :       return -sum * _phi[_j][_qp];
     170             :   }
     171             : 
     172           0 :   mooseError("Invalid type passed in");
     173             : }
     174             : 
     175             : Real
     176    20275200 : NestedKKSMultiACBulkC::computeQpOffDiagJacobian(unsigned int jvar)
     177             : {
     178             :   // first get dependence of mobility _L on other variables using parent class member function Real
     179    20275200 :   Real res = ACBulk<Real>::computeQpOffDiagJacobian(jvar);
     180             : 
     181             :   Real sum = 0.0;
     182             : 
     183             :   // Then add dependence of KKSACBulkF on other variables if other cs are the coupled variables
     184    20275200 :   auto compvar = mapJvarToCvar(jvar, _c_map);
     185    20275200 :   if (compvar >= 0)
     186             :   {
     187    13516800 :     for (const auto m : make_range(_num_c))
     188             :     {
     189             :       Real sum1 = 0.0;
     190             :       Real sum2 = 0.0;
     191             :       Real sum3 = 0.0;
     192             : 
     193    13516800 :       for (const auto n : make_range(_num_c))
     194     6758400 :         sum1 += (*_d2F1dc1db1[m][n])[_qp] * (*_dcidb[n][0][compvar])[_qp];
     195             : 
     196    27033600 :       for (const auto l : make_range(_num_j))
     197             :       {
     198    20275200 :         sum2 += (*_prop_dhjdetai[l])[_qp] * (*_prop_ci[m][l])[_qp];
     199    20275200 :         sum3 += (*_prop_dhjdetai[l])[_qp] * (*_dcidb[m][l][compvar])[_qp];
     200             :       }
     201             : 
     202     6758400 :       sum += sum1 * sum2 + (*_dF1dc1[m])[_qp] * sum3;
     203             :     }
     204             : 
     205     6758400 :     res += -_L[_qp] * sum * _phi[_j][_qp] * _test[_i][_qp];
     206             : 
     207     6758400 :     return res;
     208             :   }
     209             : 
     210             :   // if order parameters are the coupled variables
     211    13516800 :   auto etavar = mapJvarToCvar(jvar, _eta_map);
     212    13516800 :   if (etavar >= 0)
     213             :   {
     214    27033600 :     for (const auto m : make_range(_num_c))
     215             :     {
     216             :       Real sum1 = 0.0;
     217             :       Real sum2 = 0.0;
     218             :       Real sum3 = 0.0;
     219             : 
     220    27033600 :       for (const auto n : make_range(_num_c))
     221    13516800 :         sum1 += (*_d2F1dc1db1[m][n])[_qp] * (*_dcidetaj[n][0][etavar])[_qp];
     222             : 
     223    54067200 :       for (const auto l : make_range(_num_j))
     224             :       {
     225    40550400 :         sum2 += (*_prop_dhjdetai[l])[_qp] * (*_prop_ci[m][l])[_qp];
     226             : 
     227    40550400 :         sum3 += (*_prop_d2hjdetaidetap[l][etavar])[_qp] * (*_prop_ci[m][l])[_qp] +
     228    40550400 :                 (*_prop_dhjdetai[l])[_qp] * (*_dcidetaj[m][l][etavar])[_qp];
     229             :       }
     230             : 
     231    13516800 :       sum += sum1 * sum2 + (*_dF1dc1[m])[_qp] * sum3;
     232             :     }
     233             : 
     234    13516800 :     res += -_L[_qp] * sum * _phi[_j][_qp] * _test[_i][_qp];
     235             :   }
     236             : 
     237             :   // for all other vars get the coupled variable jvar is referring to
     238             :   const unsigned int cvar = mapJvarToCvar(jvar);
     239    27033600 :   for (unsigned int m = 0; m < _num_c; ++m)
     240             :   {
     241             :     Real sum1 = 0.0;
     242             : 
     243    54067200 :     for (const auto n : make_range(_num_j))
     244    40550400 :       sum1 += (*_d2F1dc1darg[m][cvar])[_qp] * (*_prop_dhjdetai[n])[_qp] * (*_prop_ci[m][n])[_qp];
     245             : 
     246    13516800 :     res += -_L[_qp] * sum1 * _phi[_j][_qp] * _test[_i][_qp];
     247             :   }
     248             : 
     249             :   return res;
     250             : }