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 "ACMultiInterface.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "MooseVariable.h"
      14             : #include "NonlinearSystem.h"
      15             : 
      16             : registerMooseObject("PhaseFieldApp", ACMultiInterface);
      17             : 
      18             : InputParameters
      19         166 : ACMultiInterface::validParams()
      20             : {
      21         166 :   InputParameters params = Kernel::validParams();
      22         166 :   params.addClassDescription("Gradient energy Allen-Cahn Kernel with cross terms");
      23         332 :   params.addRequiredCoupledVar("etas", "All eta_i order parameters of the multiphase problem");
      24         332 :   params.addRequiredParam<std::vector<MaterialPropertyName>>("kappa_names",
      25             :                                                              "The kappa used with the kernel");
      26         332 :   params.addParam<MaterialPropertyName>("mob_name", "L", "The mobility used with the kernel");
      27         166 :   return params;
      28           0 : }
      29             : 
      30          87 : ACMultiInterface::ACMultiInterface(const InputParameters & parameters)
      31             :   : Kernel(parameters),
      32          87 :     _num_etas(coupledComponents("etas")),
      33          87 :     _eta(coupledValues("etas")),
      34          87 :     _grad_eta(coupledGradients("etas")),
      35          87 :     _eta_vars(_sys.nVariables(), -1),
      36         174 :     _kappa_names(getParam<std::vector<MaterialPropertyName>>("kappa_names")),
      37          87 :     _kappa(_num_etas),
      38         261 :     _L(getMaterialProperty<Real>("mob_name"))
      39             : {
      40          87 :   if (_num_etas != _kappa_names.size())
      41           0 :     paramError("kappa_names", "Supply the same number of etas and kappa_names.");
      42             : 
      43          87 :   unsigned int nvariables = _sys.nVariables();
      44             : 
      45             :   int a = -1;
      46         348 :   for (unsigned int i = 0; i < _num_etas; ++i)
      47             :   {
      48             :     // populate lookup table form jvar to _eta index
      49         261 :     unsigned int var = coupled("etas", i);
      50         261 :     if (var < nvariables)
      51         237 :       _eta_vars[var] = i;
      52             : 
      53             :     // get the index of the variable the kernel is operating on
      54         261 :     if (coupled("etas", i) == _var.number())
      55          87 :       a = i;
      56             : 
      57             :     // get gradient prefactors
      58         261 :     _kappa[i] = &getMaterialPropertyByName<Real>(_kappa_names[i]);
      59             :   }
      60             : 
      61          87 :   if (a < 0)
      62           0 :     paramError(
      63             :         "etas", "Kernel variable must be listed in etas for ACMultiInterface kernel ", name());
      64             :   else
      65          87 :     _a = a;
      66          87 : }
      67             : 
      68             : Real
      69     7728000 : ACMultiInterface::computeQpResidual()
      70             : {
      71     7728000 :   const VariableValue & _eta_a = _u;
      72     7728000 :   const VariableGradient & _grad_eta_a = _grad_u;
      73             : 
      74             :   Real sum = 0.0;
      75    30912000 :   for (unsigned int b = 0; b < _num_etas; ++b)
      76             :   {
      77             :     // skip the diagonal term (does not contribute)
      78    23184000 :     if (b == _a)
      79     7728000 :       continue;
      80             : 
      81    15456000 :     sum += (*_kappa[b])[_qp] *
      82             :            (
      83             :                // order 1 terms
      84    15456000 :                2.0 * _test[_i][_qp] *
      85    15456000 :                    (_eta_a[_qp] * (*_grad_eta[b])[_qp] - (*_eta[b])[_qp] * _grad_eta_a[_qp]) *
      86    15456000 :                    (*_grad_eta[b])[_qp]
      87             :                // volume terms
      88    15456000 :                + (-(_eta_a[_qp] * (*_eta[b])[_qp] * _grad_test[_i][_qp] +
      89    15456000 :                     _test[_i][_qp] * (*_eta[b])[_qp] * _grad_eta_a[_qp] +
      90    15456000 :                     _test[_i][_qp] * _eta_a[_qp] * (*_grad_eta[b])[_qp]) *
      91    15456000 :                   (*_grad_eta[b])[_qp]) -
      92    15456000 :                (-((*_eta[b])[_qp] * (*_eta[b])[_qp] * _grad_test[_i][_qp] +
      93    15456000 :                   2.0 * _test[_i][_qp] * (*_eta[b])[_qp] * (*_grad_eta[b])[_qp]) *
      94             :                 _grad_eta_a[_qp]));
      95             :   }
      96             : 
      97     7728000 :   return _L[_qp] * sum;
      98             : }
      99             : 
     100             : Real
     101     3507200 : ACMultiInterface::computeQpJacobian()
     102             : {
     103             :   Real sum = 0.0;
     104    14028800 :   for (unsigned int b = 0; b < _num_etas; ++b)
     105             :   {
     106             :     // skip the diagonal term (does not contribute)
     107    10521600 :     if (b == _a)
     108     3507200 :       continue;
     109             : 
     110     7014400 :     sum += (*_kappa[b])[_qp] *
     111     7014400 :            (2.0 * _test[_i][_qp] *
     112     7014400 :                 ((_phi[_j][_qp] * (*_grad_eta[b])[_qp] - (*_eta[b])[_qp] * _grad_phi[_j][_qp]) *
     113     7014400 :                  (*_grad_eta[b])[_qp]) +
     114     7014400 :             (-(_phi[_j][_qp] * (*_eta[b])[_qp] * _grad_test[_i][_qp] +
     115     7014400 :                _test[_i][_qp] * (*_eta[b])[_qp] * _grad_phi[_j][_qp] +
     116     7014400 :                _test[_i][_qp] * _phi[_j][_qp] * (*_grad_eta[b])[_qp]) *
     117     7014400 :              (*_grad_eta[b])[_qp]) -
     118     7014400 :             (-((*_eta[b])[_qp] * (*_eta[b])[_qp] * _grad_test[_i][_qp] +
     119     7014400 :                2.0 * _test[_i][_qp] * (*_eta[b])[_qp] * (*_grad_eta[b])[_qp]) *
     120             :              _grad_phi[_j][_qp]));
     121             :   }
     122             : 
     123     3507200 :   return _L[_qp] * sum;
     124             : }
     125             : 
     126             : Real
     127     9011200 : ACMultiInterface::computeQpOffDiagJacobian(unsigned int jvar)
     128             : {
     129     9011200 :   const VariableValue & _eta_a = _u;
     130     9011200 :   const VariableGradient & _grad_eta_a = _grad_u;
     131             : 
     132     9011200 :   const int b = _eta_vars[jvar];
     133     9011200 :   if (b < 0)
     134             :     return 0.0;
     135             : 
     136     5504000 :   return _L[_qp] * (*_kappa[b])[_qp] *
     137     5504000 :          (2.0 * _test[_i][_qp] *
     138     5504000 :               ((_eta_a[_qp] * _grad_phi[_j][_qp] - _phi[_j][_qp] * _grad_eta_a[_qp]) *
     139     5504000 :                    (*_grad_eta[b])[_qp] +
     140     5504000 :                (_eta_a[_qp] * (*_grad_eta[b])[_qp] - (*_eta[b])[_qp] * _grad_eta_a[_qp]) *
     141     5504000 :                    _grad_phi[_j][_qp]) +
     142     5504000 :           (-(_eta_a[_qp] * _phi[_j][_qp] * _grad_test[_i][_qp] +
     143     5504000 :              _test[_i][_qp] * _phi[_j][_qp] * _grad_eta_a[_qp] +
     144     5504000 :              _test[_i][_qp] * _eta_a[_qp] * _grad_phi[_j][_qp]) *
     145     5504000 :                (*_grad_eta[b])[_qp] -
     146     5504000 :            (_eta_a[_qp] * (*_eta[b])[_qp] * _grad_test[_i][_qp] +
     147     5504000 :             _test[_i][_qp] * (*_eta[b])[_qp] * _grad_eta_a[_qp] +
     148             :             _test[_i][_qp] * _eta_a[_qp] * (*_grad_eta[b])[_qp]) *
     149     5504000 :                _grad_phi[_j][_qp]) -
     150     5504000 :           (-(2.0 * (*_eta[b])[_qp] * _phi[_j][_qp] * _grad_test[_i][_qp] +
     151             :              2.0 * _test[_i][_qp] *
     152             :                  (_phi[_j][_qp] * (*_grad_eta[b])[_qp] + (*_eta[b])[_qp] * _grad_phi[_j][_qp])) *
     153     5504000 :            _grad_eta_a[_qp]));
     154             : }