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 "SwitchingFunctionPenalty.h"
      11             : 
      12             : registerMooseObject("PhaseFieldApp", SwitchingFunctionPenalty);
      13             : 
      14             : InputParameters
      15          69 : SwitchingFunctionPenalty::validParams()
      16             : {
      17          69 :   InputParameters params = Kernel::validParams();
      18          69 :   params.addClassDescription(
      19             :       "Penalty kernel to constrain the sum of all switching functions in a multiphase system.");
      20         138 :   params.addParam<std::vector<MaterialPropertyName>>(
      21             :       "h_names", "Switching Function Materials that provide h(eta_i)");
      22         138 :   params.addRequiredCoupledVar("etas", "eta_i order parameters, one for each h");
      23         138 :   params.addParam<Real>("penalty", 1.0, "Penalty scaling factor");
      24          69 :   return params;
      25           0 : }
      26             : 
      27          37 : SwitchingFunctionPenalty::SwitchingFunctionPenalty(const InputParameters & parameters)
      28             :   : DerivativeMaterialInterface<Kernel>(parameters),
      29          74 :     _h_names(getParam<std::vector<MaterialPropertyName>>("h_names")),
      30          37 :     _num_h(_h_names.size()),
      31          37 :     _h(_num_h),
      32          37 :     _dh(_num_h),
      33          74 :     _penalty(getParam<Real>("penalty")),
      34          37 :     _number_of_nl_variables(_sys.nVariables()),
      35          37 :     _j_eta(_number_of_nl_variables, -1),
      36          37 :     _a(-1)
      37             : {
      38             :   // parameter check. We need exactly one eta per h
      39          74 :   if (_num_h != coupledComponents("etas"))
      40           0 :     paramError("h_names", "Need to pass in as many h_names as etas");
      41             : 
      42             :   // fetch switching functions (for the residual) and h derivatives (for the Jacobian)
      43         132 :   for (unsigned int i = 0; i < _num_h; ++i)
      44             :   {
      45          95 :     _h[i] = &getMaterialPropertyByName<Real>(_h_names[i]);
      46         190 :     _dh[i] = &getMaterialPropertyDerivative<Real>(_h_names[i], coupledName("etas", i));
      47             : 
      48             :     // generate the lookup table from j_var -> eta index
      49          95 :     unsigned int num = coupled("etas", i);
      50          95 :     if (num < _number_of_nl_variables)
      51          95 :       _j_eta[num] = i;
      52             : 
      53             :     // figure out which variable this kernel is acting on
      54          95 :     if (num == _var.number())
      55          37 :       _a = i;
      56             :   }
      57             : 
      58          37 :   if (_a < 0)
      59           0 :     paramError("etas", "Kernel variable must be listed in etas");
      60             : 
      61          37 :   _d2h = &getMaterialPropertyDerivative<Real>(_h_names[_a], _var.name(), _var.name());
      62          37 : }
      63             : 
      64             : Real
      65     2356480 : SwitchingFunctionPenalty::computeQpResidual()
      66             : {
      67             :   Real g = -1.0;
      68     7069440 :   for (unsigned int i = 0; i < _num_h; ++i)
      69     4712960 :     g += (*_h[i])[_qp];
      70             : 
      71     2356480 :   return _test[_i][_qp] * _penalty * 2.0 * g * (*_dh[_a])[_qp];
      72             : }
      73             : 
      74             : Real
      75     1505280 : SwitchingFunctionPenalty::computeQpJacobian()
      76             : {
      77             :   Real g = -1.0;
      78     4515840 :   for (unsigned int i = 0; i < _num_h; ++i)
      79     3010560 :     g += (*_h[i])[_qp];
      80             : 
      81     1505280 :   return _test[_i][_qp] * _penalty * _phi[_j][_qp] * 2.0 *
      82     1505280 :          ((*_dh[_a])[_qp] * (*_dh[_a])[_qp] + g * (*_d2h)[_qp]);
      83             : }
      84             : 
      85             : Real
      86     4515840 : SwitchingFunctionPenalty::computeQpOffDiagJacobian(unsigned int j_var)
      87             : {
      88     4515840 :   const int eta = _j_eta[j_var];
      89             : 
      90     4515840 :   if (eta >= 0)
      91     1505280 :     return _test[_i][_qp] * _penalty * _phi[_j][_qp] * 2.0 * (*_dh[eta])[_qp] * (*_dh[_a])[_qp];
      92             :   else
      93             :     return 0.0;
      94             : }