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 "ACBarrierFunction.h"
      11             : #include "libmesh/utility.h"
      12             : 
      13             : registerMooseObject("PhaseFieldApp", ACBarrierFunction);
      14             : registerMooseObject("PhaseFieldApp", ADACBarrierFunction);
      15             : 
      16             : template <bool is_ad>
      17             : InputParameters
      18         133 : ACBarrierFunctionTempl<is_ad>::validParams()
      19             : {
      20         133 :   InputParameters params = ACBarrierFunctionBase<is_ad>::validParams();
      21         266 :   params.addRequiredParam<MaterialPropertyName>(
      22             :       "gamma", "The interface profile coefficient to use with the kernel");
      23         133 :   params.addClassDescription("Allen-Cahn kernel used when 'mu' is a function of variables");
      24         133 :   return params;
      25           0 : }
      26             : 
      27             : template <bool is_ad>
      28          70 : ACBarrierFunctionTempl<is_ad>::ACBarrierFunctionTempl(const InputParameters & parameters)
      29             :   : ACBarrierFunctionBase<is_ad>(parameters),
      30         140 :     _uname(this->template getParam<NonlinearVariableName>("variable")),
      31          70 :     _gamma_name(this->template getParam<MaterialPropertyName>("gamma")),
      32          70 :     _gamma(this->template getGenericMaterialProperty<Real, is_ad>(_gamma_name)),
      33         140 :     _dmudvar(this->template getGenericMaterialProperty<Real, is_ad>(
      34          70 :         this->derivativePropertyNameFirst("mu", _uname)))
      35             : {
      36          70 : }
      37             : 
      38          36 : ACBarrierFunction::ACBarrierFunction(const InputParameters & parameters)
      39             :   : ACBarrierFunctionTempl<false>(parameters),
      40          36 :     _n_eta(_vals.size()),
      41          36 :     _d2mudvar2(getMaterialPropertyDerivative<Real>("mu", _uname, _uname)),
      42          72 :     _vname(getParam<std::vector<VariableName>>("v")),
      43          36 :     _d2mudvardeta(_n_eta),
      44          72 :     _vmap(getParameterJvarMap("v"))
      45             : {
      46          84 :   for (unsigned int i = 0; i < _n_eta; ++i)
      47          48 :     _d2mudvardeta[i] = &getMaterialPropertyDerivative<Real>("mu", _uname, _vname[i]);
      48          36 : }
      49             : 
      50             : Real
      51     2766400 : ACBarrierFunction::computeDFDOP(PFFunctionType type)
      52             : {
      53     2766400 :   switch (type)
      54             :   {
      55     1627200 :     case Residual:
      56     1627200 :       return _dmudvar[_qp] * calculateF0();
      57             : 
      58             :     case Jacobian:
      59             :     {
      60             :       Real df0dvar = 0.0;
      61     2393600 :       for (unsigned int i = 0; i < _n_eta; ++i)
      62     1254400 :         df0dvar += (*_vals[i])[_qp] * (*_vals[i])[_qp];
      63             : 
      64     1139200 :       df0dvar *= 2.0 * _gamma[_qp] * _u[_qp];
      65     1139200 :       df0dvar += _u[_qp] * _u[_qp] * _u[_qp] - _u[_qp];
      66             : 
      67     1139200 :       return (_d2mudvar2[_qp] * calculateF0() + _dmudvar[_qp] * df0dvar) * _phi[_j][_qp];
      68             :     }
      69             : 
      70           0 :     default:
      71           0 :       mooseError("Invalid type passed in");
      72             :   }
      73             : }
      74             : 
      75             : ADReal
      76    34579200 : ADACBarrierFunction::computeDFDOP()
      77             : {
      78   103737600 :   ADReal f0 = 0.25 + 0.25 * _u[_qp] * _u[_qp] * _u[_qp] * _u[_qp] - 0.5 * _u[_qp] * _u[_qp];
      79    34579200 :   ADReal sum_etaj = _u[_qp] * _u[_qp];
      80             : 
      81    69158400 :   for (unsigned int i = 0; i < _op_num; ++i)
      82             :   {
      83    34579200 :     for (unsigned int j = i + 1; j < _op_num; ++j)
      84           0 :       sum_etaj += (*_vals[j])[_qp] * (*_vals[j])[_qp];
      85             : 
      86   103737600 :     f0 += 0.25 * Utility::pow<4>((*_vals[i])[_qp]) - 0.5 * Utility::pow<2>((*_vals[i])[_qp]);
      87    69158400 :     f0 += sum_etaj * Utility::pow<2>((*_vals[i])[_qp]) * _gamma[_qp];
      88             :   }
      89    34579200 :   return _dmudvar[_qp] * f0;
      90             : }
      91             : 
      92             : Real
      93     4454400 : ACBarrierFunction::computeQpOffDiagJacobian(unsigned int jvar)
      94             : {
      95             :   const unsigned int j = mapJvarToCvar(jvar);
      96             :   Real sum_etai2 = 0.0;
      97             :   Real df0deta_base = 0.0;
      98             :   Real df0deta = 0.0;
      99             : 
     100     9267200 :   for (unsigned int i = 0; i < _n_eta; ++i)
     101     4812800 :     if (i != j)
     102      358400 :       sum_etai2 += (*_vals[i])[_qp] * (*_vals[i])[_qp];
     103             : 
     104     4454400 :   auto etavar = mapJvarToCvar(jvar, _vmap);
     105     4454400 :   if (etavar >= 0)
     106             :   {
     107     4454400 :     df0deta_base = (*_vals[etavar])[_qp] * (*_vals[etavar])[_qp] - 1.0 +
     108     4454400 :                    2.0 * _gamma[_qp] * (_u[_qp] * _u[_qp] + sum_etai2);
     109     4454400 :     df0deta = (*_vals[etavar])[_qp] * df0deta_base;
     110             : 
     111     4454400 :     return ((*_d2mudvardeta[etavar])[_qp] * calculateF0() + _dmudvar[_qp] * df0deta) *
     112     4454400 :            _phi[_j][_qp] * _test[_i][_qp] * _L[_qp];
     113             :   }
     114             :   return 0.0;
     115             : }
     116             : 
     117             : Real
     118     7220800 : ACBarrierFunction::calculateF0()
     119             : {
     120     7220800 :   Real var_phase = 0.25 * _u[_qp] * _u[_qp] * _u[_qp] * _u[_qp] - 0.5 * _u[_qp] * _u[_qp];
     121             :   Real eta_phase = 0.0;
     122             :   Real eta_interface = 0.0;
     123             :   Real sum_eta_j2;
     124             : 
     125    15201600 :   for (unsigned int i = 0; i < _n_eta; ++i)
     126             :   {
     127             :     sum_eta_j2 = 0.0;
     128     8740800 :     for (unsigned int j = i + 1; j < _n_eta; ++j)
     129      760000 :       sum_eta_j2 += (*_vals[j])[_qp] * (*_vals[j])[_qp];
     130             :     // eta_phase += 0.25 * (*_vals[i])[_qp] * (*_vals[i])[_qp] * (*_vals[i])[_qp] * (*_vals[i])[_qp]
     131             :     // -
     132             :     //             0.5 * (*_vals[i])[_qp] * (*_vals[i])[_qp];
     133     7980800 :     eta_phase += 0.25 * Utility::pow<4>((*_vals[i])[_qp]) - 0.5 * Utility::pow<2>((*_vals[i])[_qp]);
     134     7980800 :     eta_interface +=
     135     7980800 :         (_u[_qp] * _u[_qp] + sum_eta_j2) * (*_vals[i])[_qp] * (*_vals[i])[_qp] * _gamma[_qp];
     136             :   }
     137     7220800 :   return 0.25 + var_phase + eta_phase + eta_interface;
     138             : }
     139             : 
     140             : template class ACBarrierFunctionTempl<false>;
     141             : template class ACBarrierFunctionTempl<true>;