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 "SLKKSMultiACBulkC.h"
      11             : 
      12             : registerMooseObject("PhaseFieldApp", SLKKSMultiACBulkC);
      13             : 
      14             : InputParameters
      15          56 : SLKKSMultiACBulkC::validParams()
      16             : {
      17          56 :   auto params = SLKKSMultiPhaseBase::validParams();
      18          56 :   params.addClassDescription("Multi-phase SLKKS model kernel for the bulk Allen-Cahn. "
      19             :                              "This includes all terms dependent on chemical potential.");
      20         112 :   params.addRequiredParam<MaterialPropertyName>(
      21             :       "F", "Phase free energy function that is a function of 'c'");
      22         112 :   params.addRequiredCoupledVar("c", "Concentration variable F depends on");
      23         112 :   params.addCoupledVar("eta_i",
      24             :                        "Order parameter that derivatives are taken with respect to (kernel "
      25             :                        "variable is used if this is not specified)");
      26         112 :   params.addParam<MaterialPropertyName>("mob_name", "L", "The mobility used with the kernel");
      27          56 :   return params;
      28           0 : }
      29             : 
      30          30 : SLKKSMultiACBulkC::SLKKSMultiACBulkC(const InputParameters & parameters)
      31             :   : SLKKSMultiPhaseBase(parameters),
      32          30 :     _c_name(coupledName("c", 0)),
      33          30 :     _lagrange(isCoupled("eta_i")),
      34          60 :     _etai_name(_lagrange ? coupledName("eta_i", 0) : _var.name()),
      35          30 :     _etai_var(_lagrange ? coupled("eta_i") : _var.number()),
      36          30 :     _cs_names(_ncs),
      37          30 :     _prop_dFdc(getMaterialPropertyDerivative<Real>("F", _c_name)),
      38          30 :     _prop_d2Fdcdcs(_ncs),
      39          30 :     _prop_dhdni(_nh),
      40          30 :     _prop_d2hdnidn(_nh),
      41          30 :     _l_cs(-1),
      42          30 :     _l_etai(-1),
      43          90 :     _mob(getMaterialProperty<Real>("mob_name"))
      44             : {
      45             :   // Determine position of the selected concentration variable
      46         150 :   for (std::size_t i = 0; i < _ncs; ++i)
      47             :   {
      48         120 :     _cs_names[i] = coupledName("cs", i);
      49         120 :     if (coupled("cs", i) == _c_var)
      50          30 :       _l_cs = i;
      51             :   }
      52             : 
      53             :   // Check to make sure the nonlinear variable is in the cs list
      54          30 :   if (_l_cs < 0)
      55           0 :     paramError("cs", "One of the listed variables must be the 'c' variable");
      56             : 
      57             :   // Determine position of the selected concentration variable
      58          90 :   for (std::size_t i = 0; i < _neta; ++i)
      59          60 :     if (coupled("eta", i) == _etai_var)
      60          30 :       _l_etai = i;
      61             : 
      62             :   // Check to make sure the nonlinear variable or eta_i (if supplied) is in the eta list
      63          30 :   if (_l_etai < 0)
      64           0 :     paramError("eta_i",
      65             :                "Either eta_i or the kernel variable must be one of the listed 'eta' variables");
      66             : 
      67             :   // Iterate over all sublattice concentrations
      68         150 :   for (std::size_t i = 0; i < _ncs; ++i)
      69             :     // get second partial derivatives w.r.t. c and cs (only if c and cs are in the same phase)
      70         240 :     _prop_d2Fdcdcs[i] = (_phase[i] == _phase[_l_cs])
      71         180 :                             ? &getMaterialPropertyDerivative<Real>("F", _c_name, _cs_names[i])
      72             :                             : nullptr;
      73             : 
      74             :   // fetch all switching function derivatives
      75          90 :   for (std::size_t i = 0; i < _nh; ++i)
      76             :   {
      77          60 :     _prop_dhdni[i] = &getMaterialPropertyDerivativeByName<Real>(_h_names[i], _etai_name);
      78             : 
      79          60 :     _prop_d2hdnidn[i].resize(_neta);
      80         180 :     for (std::size_t j = 0; j < _neta; ++j)
      81         120 :       _prop_d2hdnidn[i][j] =
      82         240 :           &getMaterialPropertyDerivativeByName<Real>(_h_names[i], _etai_name, _eta_names[j]);
      83             :   }
      84          30 : }
      85             : 
      86             : Real
      87       70320 : SLKKSMultiACBulkC::precomputeQpResidual()
      88             : { // sum over phases
      89             :   std::size_t k = 0;
      90             :   Real sum = 0.0;
      91      210960 :   for (std::size_t i = 0; i < _nh; ++i)
      92             :   {
      93             :     // sum sublattice concentrations
      94             :     Real csum = 0.0;
      95      421920 :     for (unsigned int j = 0; j < _ns[i]; ++j)
      96             :     {
      97      281280 :       csum += (*_cs[k])[_qp] * _a_cs[k];
      98      281280 :       k++;
      99             :     }
     100      140640 :     sum += (*_prop_dhdni[i])[_qp] * csum;
     101             :   }
     102             : 
     103       70320 :   return -_mob[_qp] * _prop_dFdc[_qp] / _a_cs[_l_cs] * sum;
     104             : }
     105             : 
     106             : Real
     107      116400 : SLKKSMultiACBulkC::precomputeQpJacobian()
     108             : {
     109             :   // For when this kernel is used in the Lagrange multiplier equation
     110      116400 :   if (_lagrange)
     111             :     return 0.0;
     112             : 
     113             :   // For when eta_i is the nonlinear variable
     114             :   std::size_t k = 0;
     115             :   Real sum = 0.0;
     116      349200 :   for (std::size_t i = 0; i < _nh; ++i)
     117             :   {
     118             :     // sum sublattice concentrations
     119             :     Real csum = 0.0;
     120      698400 :     for (unsigned int j = 0; j < _ns[i]; ++j)
     121             :     {
     122      465600 :       csum += (*_cs[k])[_qp] * _a_cs[k];
     123      465600 :       k++;
     124             :     }
     125      232800 :     sum += (*_prop_d2hdnidn[i][_l_etai])[_qp] * csum;
     126             :   }
     127             : 
     128      116400 :   return -_mob[_qp] * _phi[_j][_qp] * _prop_dFdc[_qp] / _a_cs[_l_cs] * sum;
     129             : }
     130             : 
     131             : Real
     132     1164000 : SLKKSMultiACBulkC::computeQpOffDiagJacobian(unsigned int jvar)
     133             : {
     134             :   // concentration variables
     135     1164000 :   auto csvar = mapJvarToCvar(jvar, _cs_map);
     136     1164000 :   if (csvar >= 0)
     137             :   {
     138             :     // does F depend on the csvar? if so we need to apply the product rule
     139      931200 :     if (_phase[csvar] == _phase[_l_cs])
     140             :     {
     141             :       // sum over phases
     142             :       std::size_t k = 0;
     143             :       Real sum = 0.0;
     144      698400 :       for (std::size_t i = 0; i < _nh; ++i)
     145             :       {
     146             :         // sum sublattice concentrations
     147             :         Real csum = 0.0;
     148     1396800 :         for (unsigned int j = 0; j < _ns[i]; ++j)
     149             :         {
     150      931200 :           csum += (*_cs[k])[_qp] * _a_cs[k];
     151      931200 :           k++;
     152             :         }
     153      465600 :         sum += (*_prop_dhdni[i])[_qp] * csum;
     154             :       }
     155             : 
     156      232800 :       return -_mob[_qp] *
     157      232800 :              ((*_prop_d2Fdcdcs[csvar])[_qp] * sum +
     158      232800 :               _prop_dFdc[_qp] * (*_prop_dhdni[_phase[csvar]])[_qp] * _a_cs[csvar]) /
     159      232800 :              _a_cs[_l_cs] * _phi[_j][_qp] * _test[_i][_qp];
     160             :     }
     161             : 
     162      698400 :     return -_mob[_qp] * _prop_dFdc[_qp] / _a_cs[_l_cs] * (*_prop_dhdni[_phase[csvar]])[_qp] *
     163      698400 :            _a_cs[csvar] * _phi[_j][_qp] * _test[_i][_qp];
     164             :   }
     165             : 
     166             :   // order parameters
     167      232800 :   auto etavar = mapJvarToCvar(jvar, _eta_map);
     168      232800 :   if (etavar >= 0)
     169             :   {
     170             :     // sum over phases
     171             :     std::size_t k = 0;
     172             :     Real sum = 0.0;
     173      698400 :     for (std::size_t i = 0; i < _nh; ++i)
     174             :     {
     175             :       // sum sublattice concentrations
     176             :       Real csum = 0.0;
     177     1396800 :       for (unsigned int j = 0; j < _ns[i]; ++j)
     178             :       {
     179      931200 :         csum += (*_cs[k])[_qp] * _a_cs[k];
     180      931200 :         k++;
     181             :       }
     182      465600 :       sum += (*_prop_d2hdnidn[i][etavar])[_qp] * csum;
     183             :     }
     184             : 
     185      232800 :     return -_mob[_qp] * _prop_dFdc[_qp] / _a_cs[_l_cs] * sum * _phi[_j][_qp] * _test[_i][_qp];
     186             :   }
     187             : 
     188             :   return 0.0;
     189             : }