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