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 "SwitchingFunctionConstraintLagrange.h"
      11             : 
      12             : registerMooseObject("PhaseFieldApp", SwitchingFunctionConstraintLagrange);
      13             : 
      14             : InputParameters
      15         155 : SwitchingFunctionConstraintLagrange::validParams()
      16             : {
      17         155 :   InputParameters params = Kernel::validParams();
      18         155 :   params.addClassDescription("Lagrange multiplier kernel to constrain the sum of all switching "
      19             :                              "functions in a multiphase system. This kernel acts on the Lagrange "
      20             :                              "multiplier variable.");
      21         310 :   params.addParam<std::vector<MaterialPropertyName>>("h_names", "Switching function materials");
      22         310 :   params.addRequiredCoupledVar("etas", "eta order parameters");
      23         310 :   params.addParam<Real>("epsilon", 1e-9, "Shift factor to avoid a zero pivot");
      24         155 :   return params;
      25           0 : }
      26             : 
      27          81 : SwitchingFunctionConstraintLagrange::SwitchingFunctionConstraintLagrange(
      28          81 :     const InputParameters & parameters)
      29             :   : DerivativeMaterialInterface<JvarMapKernelInterface<Kernel>>(parameters),
      30         162 :     _h_names(getParam<std::vector<MaterialPropertyName>>("h_names")),
      31          81 :     _num_h(_h_names.size()),
      32          81 :     _h(_num_h),
      33          81 :     _dh(_num_h),
      34          81 :     _eta_map(getParameterJvarMap("etas")),
      35         243 :     _epsilon(getParam<Real>("epsilon"))
      36             : {
      37             :   // parameter check. We need exactly one eta per h
      38         162 :   if (_num_h != coupledComponents("etas"))
      39           0 :     paramError("etas", "Need to pass in as many etas as h_names");
      40             : 
      41             :   // fetch switching functions (for the residual) and h derivatives (for the Jacobian)
      42         291 :   for (std::size_t i = 0; i < _num_h; ++i)
      43             :   {
      44         210 :     _h[i] = &getMaterialPropertyByName<Real>(_h_names[i]);
      45             : 
      46         210 :     _dh[i].resize(_num_h);
      47         774 :     for (std::size_t j = 0; j < _num_h; ++j)
      48        1128 :       _dh[i][j] = &getMaterialPropertyDerivative<Real>(_h_names[i], coupledName("etas", j));
      49             :   }
      50          81 : }
      51             : 
      52             : Real
      53     9491080 : SwitchingFunctionConstraintLagrange::computeQpResidual()
      54             : {
      55     9491080 :   Real g = -_epsilon * _u[_qp] - 1.0;
      56    36092440 :   for (std::size_t i = 0; i < _num_h; ++i)
      57    26601360 :     g += (*_h[i])[_qp];
      58             : 
      59     9491080 :   return _test[_i][_qp] * g;
      60             : }
      61             : 
      62             : Real
      63     4975520 : SwitchingFunctionConstraintLagrange::computeQpJacobian()
      64             : {
      65     4975520 :   return _test[_i][_qp] * -_epsilon * _phi[_j][_qp];
      66             : }
      67             : 
      68             : Real
      69    13279040 : SwitchingFunctionConstraintLagrange::computeQpOffDiagJacobian(unsigned int jvar)
      70             : {
      71    13279040 :   auto eta = mapJvarToCvar(jvar, _eta_map);
      72    13279040 :   if (eta >= 0)
      73             :   {
      74             :     Real g = 0.0;
      75    51331520 :     for (std::size_t i = 0; i < _num_h; ++i)
      76    38052480 :       g += (*_dh[i][eta])[_qp] * _phi[_j][_qp];
      77    13279040 :     return g * _test[_i][_qp];
      78             :   }
      79             : 
      80             :   return 0.0;
      81             : }