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