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         138 : NestedKKSMultiACBulkC::validParams()
      16             : {
      17         138 :   InputParameters params = KKSMultiACBulkBase::validParams();
      18         138 :   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         276 :   params.addRequiredCoupledVar("global_cs", "Global concentrations, for example, c, b.");
      21         276 :   params.addRequiredCoupledVar("all_etas", "Order parameters.");
      22         276 :   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         138 :   return params;
      27           0 : }
      28             : 
      29          72 : NestedKKSMultiACBulkC::NestedKKSMultiACBulkC(const InputParameters & parameters)
      30             :   : KKSMultiACBulkBase(parameters),
      31          72 :     _c_names(coupledNames("global_cs")),
      32          72 :     _c_map(getParameterJvarMap("global_cs")),
      33          72 :     _num_c(coupledComponents("global_cs")),
      34          72 :     _eta_names(coupledNames("all_etas")),
      35          72 :     _eta_map(getParameterJvarMap("all_etas")),
      36          72 :     _k(-1),
      37         144 :     _ci_names(getParam<std::vector<MaterialPropertyName>>("ci_names")),
      38          72 :     _ci_name_matrix(_num_c),
      39          72 :     _prop_ci(_num_c),
      40          72 :     _dcidetaj(_num_c),
      41          72 :     _dcidb(_num_c),
      42          72 :     _prop_d2hjdetaidetap(_num_j),
      43          72 :     _dF1dc1(_num_c),
      44          72 :     _d2F1dc1db1(_num_c),
      45         144 :     _d2F1dc1darg(_n_args)
      46             : {
      47         288 :   for (unsigned int i = 0; i < _num_j; ++i)
      48             :   {
      49             :     // get order parameter names and variable indices
      50         432 :     _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         216 :     if (coupled("all_etas", i) == _var.number())
      54          72 :       _k = i;
      55             :   }
      56             : 
      57             :   // initialize _ci_name_matrix and _prop_ci for easy reference
      58         144 :   for (unsigned int m = 0; m < _num_c; ++m)
      59             :   {
      60          72 :     _ci_name_matrix[m].resize(_num_j);
      61          72 :     _prop_ci[m].resize(_num_j);
      62             : 
      63         288 :     for (const auto n : make_range(_num_j))
      64             :     {
      65         216 :       _ci_name_matrix[m][n] = _ci_names[m * _num_j + n];
      66         216 :       _prop_ci[m][n] = &getMaterialPropertyByName<Real>(_ci_name_matrix[m][n]);
      67             :     }
      68             :   }
      69             : 
      70             :   // _dcidetaj and _dcidb are computed in KKSPhaseConcentrationMultiPhaseDerivatives
      71         144 :   for (const auto m : make_range(_num_c))
      72             :   {
      73          72 :     _dcidetaj[m].resize(_num_j);
      74          72 :     _dcidb[m].resize(_num_j);
      75             : 
      76         288 :     for (const auto n : make_range(_num_j))
      77             :     {
      78         216 :       _dcidetaj[m][n].resize(_num_j);
      79         216 :       _dcidb[m][n].resize(_num_c);
      80             : 
      81         864 :       for (const auto l : make_range(_num_j))
      82         648 :         _dcidetaj[m][n][l] =
      83        1296 :             &getMaterialPropertyDerivative<Real>(_ci_names[n + m * _num_j], _eta_names[l]);
      84             : 
      85         432 :       for (const auto l : make_range(_num_c))
      86         216 :         _dcidb[m][n][l] =
      87         432 :             &getMaterialPropertyDerivative<Real>(_ci_names[n + m * _num_j], _c_names[l]);
      88             :     }
      89             :   }
      90             : 
      91         288 :   for (const auto m : make_range(_num_j))
      92             :   {
      93         216 :     _prop_d2hjdetaidetap[m].resize(_num_j);
      94             : 
      95         864 :     for (const auto n : make_range(_num_j))
      96         648 :       _prop_d2hjdetaidetap[m][n] =
      97        1296 :           &getMaterialPropertyDerivative<Real>(_hj_names[m], _eta_names[_k], _eta_names[n]);
      98             :   }
      99             : 
     100             :   // _dF1dc1 and _d2F1dc1db1 are computed in KKSPhaseConcentrationMultiPhaseMaterial
     101         144 :   for (const auto m : make_range(_num_c))
     102             :   {
     103          72 :     _dF1dc1[m] = &getMaterialPropertyDerivative<Real>(_Fj_names[0], _ci_names[m * _num_j]);
     104          72 :     _d2F1dc1db1[m].resize(_num_c);
     105             : 
     106         144 :     for (const auto n : make_range(_num_c))
     107          72 :       _d2F1dc1db1[m][n] = &getMaterialPropertyDerivative<Real>(
     108          72 :           _Fj_names[0], _ci_name_matrix[m][0], _ci_name_matrix[n][0]);
     109             :   }
     110             : 
     111             :   // _d2F1dc1darg are computed in KKSPhaseConcentrationMultiPhaseMaterial
     112         144 :   for (const auto m : make_range(_num_c))
     113             :   {
     114          72 :     _d2F1dc1darg[m].resize(_n_args);
     115             : 
     116         648 :     for (const auto n : make_range(_n_args))
     117         576 :       _d2F1dc1darg[m][n] =
     118         576 :           &getMaterialPropertyDerivative<Real>(_Fj_names[0], _ci_name_matrix[m][0], n);
     119             :   }
     120          72 : }
     121             : 
     122             : Real
     123    31324800 : NestedKKSMultiACBulkC::computeDFDOP(PFFunctionType type)
     124             : {
     125             :   Real sum = 0.0;
     126             : 
     127    31324800 :   switch (type)
     128             :   {
     129    29328000 :     case Residual:
     130             : 
     131    58656000 :       for (const auto m : make_range(_num_c))
     132             :       {
     133             :         Real sum1 = 0.0;
     134             : 
     135   117312000 :         for (const auto n : make_range(_num_j))
     136    87984000 :           sum1 += (*_prop_dhjdetai[n])[_qp] * (*_prop_ci[m][n])[_qp];
     137             : 
     138    29328000 :         sum += (*_dF1dc1[m])[_qp] * sum1;
     139             :       }
     140             : 
     141    29328000 :       return -sum;
     142             : 
     143     1996800 :     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     1996800 :       if (_etai_var != _var.number())
     147             :         return 0.0;
     148             : 
     149     3993600 :       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     3993600 :         for (const auto n : make_range(_num_c))
     156     1996800 :           sum1 += (*_d2F1dc1db1[m][n])[_qp] * (*_dcidetaj[n][0][_k])[_qp];
     157             : 
     158     7987200 :         for (const auto l : make_range(_num_j))
     159             :         {
     160     5990400 :           sum2 += (*_prop_dhjdetai[l])[_qp] * (*_prop_ci[m][l])[_qp];
     161             : 
     162     5990400 :           sum3 += (*_prop_d2hjdetaidetap[l][_k])[_qp] * (*_prop_ci[m][l])[_qp] +
     163     5990400 :                   (*_prop_dhjdetai[l])[_qp] * (*_dcidetaj[m][l][_k])[_qp];
     164             :         }
     165             : 
     166     1996800 :         sum += sum1 * sum2 + (*_dF1dc1[m])[_qp] * sum3;
     167             :       }
     168             : 
     169     1996800 :       return -sum * _phi[_j][_qp];
     170             :   }
     171             : 
     172           0 :   mooseError("Invalid type passed in");
     173             : }
     174             : 
     175             : Real
     176    23961600 : NestedKKSMultiACBulkC::computeQpOffDiagJacobian(unsigned int jvar)
     177             : {
     178             :   // first get dependence of mobility _L on other variables using parent class member function Real
     179    23961600 :   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    23961600 :   auto compvar = mapJvarToCvar(jvar, _c_map);
     185    23961600 :   if (compvar >= 0)
     186             :   {
     187    15974400 :     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    15974400 :       for (const auto n : make_range(_num_c))
     194     7987200 :         sum1 += (*_d2F1dc1db1[m][n])[_qp] * (*_dcidb[n][0][compvar])[_qp];
     195             : 
     196    31948800 :       for (const auto l : make_range(_num_j))
     197             :       {
     198    23961600 :         sum2 += (*_prop_dhjdetai[l])[_qp] * (*_prop_ci[m][l])[_qp];
     199    23961600 :         sum3 += (*_prop_dhjdetai[l])[_qp] * (*_dcidb[m][l][compvar])[_qp];
     200             :       }
     201             : 
     202     7987200 :       sum += sum1 * sum2 + (*_dF1dc1[m])[_qp] * sum3;
     203             :     }
     204             : 
     205     7987200 :     res += -_L[_qp] * sum * _phi[_j][_qp] * _test[_i][_qp];
     206             : 
     207     7987200 :     return res;
     208             :   }
     209             : 
     210             :   // if order parameters are the coupled variables
     211    15974400 :   auto etavar = mapJvarToCvar(jvar, _eta_map);
     212    15974400 :   if (etavar >= 0)
     213             :   {
     214    31948800 :     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    31948800 :       for (const auto n : make_range(_num_c))
     221    15974400 :         sum1 += (*_d2F1dc1db1[m][n])[_qp] * (*_dcidetaj[n][0][etavar])[_qp];
     222             : 
     223    63897600 :       for (const auto l : make_range(_num_j))
     224             :       {
     225    47923200 :         sum2 += (*_prop_dhjdetai[l])[_qp] * (*_prop_ci[m][l])[_qp];
     226             : 
     227    47923200 :         sum3 += (*_prop_d2hjdetaidetap[l][etavar])[_qp] * (*_prop_ci[m][l])[_qp] +
     228    47923200 :                 (*_prop_dhjdetai[l])[_qp] * (*_dcidetaj[m][l][etavar])[_qp];
     229             :       }
     230             : 
     231    15974400 :       sum += sum1 * sum2 + (*_dF1dc1[m])[_qp] * sum3;
     232             :     }
     233             : 
     234    15974400 :     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    31948800 :   for (unsigned int m = 0; m < _num_c; ++m)
     240             :   {
     241             :     Real sum1 = 0.0;
     242             : 
     243    63897600 :     for (const auto n : make_range(_num_j))
     244    47923200 :       sum1 += (*_d2F1dc1darg[m][cvar])[_qp] * (*_prop_dhjdetai[n])[_qp] * (*_prop_ci[m][n])[_qp];
     245             : 
     246    15974400 :     res += -_L[_qp] * sum1 * _phi[_j][_qp] * _test[_i][_qp];
     247             :   }
     248             : 
     249             :   return res;
     250             : }