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 "NestedKKSMultiACBulkF.h"
      11             : 
      12             : registerMooseObject("PhaseFieldApp", NestedKKSMultiACBulkF);
      13             : 
      14             : InputParameters
      15         138 : NestedKKSMultiACBulkF::validParams()
      16             : {
      17         138 :   InputParameters params = KKSMultiACBulkBase::validParams();
      18         138 :   params.addClassDescription("KKS model kernel (part 1 of 2) for the Bulk Allen-Cahn. This "
      19             :                              "includes all terms NOT dependent on chemical potential.");
      20         276 :   params.addRequiredCoupledVar("global_cs", "The interpolated concentrations c, b, etc.");
      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         276 :   params.addRequiredParam<Real>("wi", "Double well height parameter.");
      27         276 :   params.addRequiredParam<MaterialPropertyName>(
      28             :       "gi_name", "Base name for the double well function g_i(eta_i).");
      29         138 :   return params;
      30           0 : }
      31             : 
      32          72 : NestedKKSMultiACBulkF::NestedKKSMultiACBulkF(const InputParameters & parameters)
      33             :   : KKSMultiACBulkBase(parameters),
      34          72 :     _c_names(coupledNames("global_cs")),
      35          72 :     _c_map(getParameterJvarMap("global_cs")),
      36          72 :     _num_c(coupledComponents("global_cs")),
      37          72 :     _eta_names(coupledNames("all_etas")),
      38          72 :     _eta_map(getParameterJvarMap("all_etas")),
      39          72 :     _k(-1),
      40         144 :     _ci_names(getParam<std::vector<MaterialPropertyName>>("ci_names")),
      41          72 :     _dcidetaj(_num_c),
      42          72 :     _dcidb(_num_c),
      43         144 :     _wi(getParam<Real>("wi")),
      44          72 :     _gi_name(getParam<MaterialPropertyName>("gi_name")),
      45          72 :     _dgi(getMaterialPropertyDerivative<Real>("gi_name", _etai_name)),
      46          72 :     _d2gi(getMaterialPropertyDerivative<Real>("gi_name", _etai_name, _etai_name)),
      47          72 :     _d2hjdetaidetap(_num_j),
      48          72 :     _dF1dc1(_num_c),
      49         144 :     _dFidarg(_num_j)
      50             : {
      51         288 :   for (const auto i : make_range(_num_j))
      52             :   {
      53             :     // get order parameter names and variable indices
      54         432 :     _eta_names[i] = getVar("all_etas", i)->name();
      55             : 
      56             :     // Set _k to the position of the nonlinear variable in the list of etaj's
      57         216 :     if (coupled("all_etas", i) == _var.number())
      58          72 :       _k = i;
      59             :   }
      60             : 
      61             :   // _dcideta and _dcidb are computed in KKSPhaseConcentrationDerivatives
      62         144 :   for (const auto m : make_range(_num_c))
      63             :   {
      64          72 :     _dcidetaj[m].resize(_num_j);
      65          72 :     _dcidb[m].resize(_num_j);
      66             : 
      67         288 :     for (const auto n : make_range(_num_j))
      68             :     {
      69         216 :       _dcidetaj[m][n].resize(_num_j);
      70         216 :       _dcidb[m][n].resize(_num_c);
      71             : 
      72         864 :       for (const auto l : make_range(_num_j))
      73         648 :         _dcidetaj[m][n][l] =
      74        1296 :             &getMaterialPropertyDerivative<Real>(_ci_names[n + m * _num_j], _eta_names[l]);
      75             : 
      76         432 :       for (const auto l : make_range(_num_c))
      77         216 :         _dcidb[m][n][l] =
      78         432 :             &getMaterialPropertyDerivative<Real>(_ci_names[n + m * _num_j], _c_names[l]);
      79             :     }
      80             :   }
      81             : 
      82             :   // _dF1dc1 is computed in KKSPhaseConcentrationMaterial
      83         144 :   for (const auto m : make_range(_num_c))
      84          72 :     _dF1dc1[m] = &getMaterialPropertyDerivative<Real>(_Fj_names[0], _ci_names[m * _num_j]);
      85             : 
      86         288 :   for (const auto m : make_range(_num_j))
      87             :   {
      88         216 :     _d2hjdetaidetap[m].resize(_num_j);
      89             : 
      90         864 :     for (const auto n : make_range(_num_j))
      91         648 :       _d2hjdetaidetap[m][n] =
      92        1296 :           &getMaterialPropertyDerivative<Real>(_hj_names[m], _eta_names[_k], _eta_names[n]);
      93             :   }
      94             : 
      95             :   // _dFidarg is computed in KKSPhaseConcentrationMaterial
      96         288 :   for (const auto m : make_range(_num_j))
      97             :   {
      98         216 :     _dFidarg[m].resize(_n_args);
      99             : 
     100        1944 :     for (const auto n : make_range(_n_args))
     101        1728 :       _dFidarg[m][n] = &getMaterialPropertyDerivative<Real>(_Fj_names[m], m);
     102             :   }
     103          72 : }
     104             : 
     105             : Real
     106    31324800 : NestedKKSMultiACBulkF::computeDFDOP(PFFunctionType type)
     107             : {
     108             :   Real sum = 0.0;
     109             : 
     110    31324800 :   switch (type)
     111             :   {
     112    29328000 :     case Residual:
     113   117312000 :       for (const auto m : make_range(_num_j))
     114    87984000 :         sum += (*_prop_dhjdetai[m])[_qp] * (*_prop_Fj[m])[_qp];
     115             : 
     116    29328000 :       return sum + _wi * _dgi[_qp];
     117             : 
     118     1996800 :     case Jacobian:
     119             :       // For when this kernel is used in the Lagrange multiplier equation. In that case the
     120             :       // Lagrange multiplier is the nonlinear variable
     121     1996800 :       if (_etai_var != _var.number())
     122             :         return 0.0;
     123             : 
     124             :       // if eta_i is the nonlinear variable
     125     7987200 :       for (const auto m : make_range(_num_j))
     126             :       {
     127             :         Real sum1 = 0.0;
     128             : 
     129    11980800 :         for (const auto n : make_range(_num_c))
     130     5990400 :           sum1 += (*_dF1dc1[n])[_qp] * (*_dcidetaj[n][m][_k])[_qp];
     131             : 
     132     5990400 :         sum +=
     133     5990400 :             (*_d2hjdetaidetap[m][_k])[_qp] * (*_prop_Fj[m])[_qp] + (*_prop_dhjdetai[m])[_qp] * sum1;
     134             :       }
     135             : 
     136     1996800 :       return _phi[_j][_qp] * (sum + _wi * _d2gi[_qp]);
     137             :   }
     138             : 
     139           0 :   mooseError("Invalid type passed in");
     140             : }
     141             : 
     142             : Real
     143    23961600 : NestedKKSMultiACBulkF::computeQpOffDiagJacobian(unsigned int jvar)
     144             : {
     145             :   // first get dependence of mobility _L on other variables using parent class member function Real
     146    23961600 :   Real res = ACBulk<Real>::computeQpOffDiagJacobian(jvar);
     147             : 
     148             :   Real sum = 0.0;
     149             : 
     150             :   // if concentrations are the coupled variables
     151    23961600 :   auto compvar = mapJvarToCvar(jvar, _c_map);
     152    23961600 :   if (compvar >= 0)
     153             :   {
     154    31948800 :     for (const auto m : make_range(_num_j))
     155             :     {
     156             :       Real sum1 = 0.0;
     157             : 
     158    47923200 :       for (const auto n : make_range(_num_c))
     159    23961600 :         sum1 += (*_dF1dc1[n])[_qp] * (*_dcidb[n][m][compvar])[_qp];
     160             : 
     161    23961600 :       sum += (*_prop_dhjdetai[m])[_qp] * sum1;
     162             :     }
     163             : 
     164     7987200 :     res += _L[_qp] * sum * _phi[_j][_qp] * _test[_i][_qp];
     165             : 
     166     7987200 :     return res;
     167             :   }
     168             : 
     169             :   // if order parameters are the coupled variables
     170    15974400 :   auto etavar = mapJvarToCvar(jvar, _eta_map);
     171    15974400 :   if (etavar >= 0)
     172             :   {
     173    63897600 :     for (const auto m : make_range(_num_j))
     174             :     {
     175             :       Real sum1 = 0.0;
     176             : 
     177    95846400 :       for (const auto n : make_range(_num_c))
     178    47923200 :         sum1 += (*_dF1dc1[n])[_qp] * (*_dcidetaj[n][m][etavar])[_qp];
     179             : 
     180    47923200 :       sum += (*_d2hjdetaidetap[m][etavar])[_qp] * (*_prop_Fj[m])[_qp] +
     181    47923200 :              (*_prop_dhjdetai[m])[_qp] * sum1;
     182             :     }
     183             : 
     184    15974400 :     res += _L[_qp] * sum * _phi[_j][_qp] * _test[_i][_qp];
     185             : 
     186    15974400 :     return res;
     187             :   }
     188             : 
     189             :   // Handle the case when this kernel is used in the Lagrange multiplier equation. In this case
     190             :   // the second derivative of the barrier function contributes to the off-diagonal Jacobian
     191           0 :   if (jvar == _etai_var)
     192             :   {
     193           0 :     sum += _wi * _d2gi[_qp];
     194             : 
     195           0 :     res += _L[_qp] * sum * _phi[_j][_qp] * _test[_i][_qp];
     196             :   }
     197             : 
     198             :   // for all other vars get the coupled variable jvar is referring to
     199             :   const unsigned int cvar = mapJvarToCvar(jvar);
     200           0 :   for (const auto m : make_range(_num_j))
     201           0 :     res -= _L[_qp] * (*_prop_dhjdetai[m])[_qp] * (*_dFidarg[m][cvar])[_qp] * _phi[_j][_qp] *
     202           0 :            _test[_i][_qp];
     203             : 
     204             :   return res;
     205             : }