Line data    Source code

       1             : //* This file is part of the RACCOON application
       2             : //* being developed at Dolbow lab at Duke University
       3             : //* http://dolbow.pratt.duke.edu
       4             : 
       5             : #include "CrackGeometricFunction.h"
       6             : 
       7             : registerMooseObject("raccoonApp", CrackGeometricFunction);
       8             : 
       9             : InputParameters
      10         330 : CrackGeometricFunction::validParams()
      11             : {
      12         330 :   InputParameters params = CustomParsedFunctionBase::validParams();
      13         330 :   params.addClassDescription(
      14             :       "This is a wrapper of ADDerivativeParsedMaterial to conveniently define a crack geometric "
      15             :       "function. The initial derivative as well as the normalization constant are automatically "
      16             :       "populated given the function definition.");
      17             : 
      18         330 :   params.set<std::string>("property_name") = "alpha";
      19         660 :   params.addRequiredCoupledVar("phase_field", "The phase-field variable");
      20             : 
      21         330 :   params.set<unsigned int>("derivative_order") = 1;
      22         330 :   params.suppressParameter<unsigned int>("derivative_order");
      23             : 
      24         660 :   params.addParam<MaterialPropertyName>("initial_derivative",
      25             :                                         "xi",
      26             :                                         "Name of the material to store the initial slope of the "
      27             :                                         "crack geometric function, $\\alpha(d=0)$");
      28             : 
      29         660 :   params.addParam<MaterialPropertyName>("normalization_constant",
      30             :                                         "c0",
      31             :                                         "Name of the material to store the normalization constant, "
      32             :                                         "$4\\int_0^1 \\sqrt{\\alpha(s)} \\diff{s}$");
      33             : 
      34         660 :   params.addParam<Real>(
      35         660 :       "tolerance", 1e-8, "Tolerance of the numerically computed normalization constant");
      36         660 :   params.addParam<unsigned int>("maximum_iterations",
      37         660 :                                 1e9,
      38             :                                 "Maximum number of iterations allowed for the numerical "
      39             :                                 "computation of the normalization constant");
      40         330 :   return params;
      41           0 : }
      42             : 
      43         144 : CrackGeometricFunction::CrackGeometricFunction(const InputParameters & parameters)
      44             :   : CustomParsedFunctionBase(parameters),
      45         288 :     _d_idx(argIndex(coupled("phase_field"))),
      46         288 :     _xi(declareADProperty<Real>(getParam<MaterialPropertyName>("initial_derivative"))),
      47         288 :     _c0(declareADProperty<Real>(getParam<MaterialPropertyName>("normalization_constant"))),
      48         288 :     _tolerance(getParam<Real>("tolerance")),
      49         432 :     _max_its(getParam<unsigned int>("maximum_iterations"))
      50             : {
      51             :   // set d to 0 and evaluate the first derivative of the crack geometric function
      52         144 :   _func_params[_d_idx] = 0;
      53         144 :   _xi_0 = evaluate(_derivatives[0]._F, _name);
      54         144 :   _c0_0 = computeNormalizationConstant();
      55         144 : }
      56             : 
      57             : ADReal
      58         144 : CrackGeometricFunction::computeNormalizationConstant()
      59             : {
      60             :   // We use an adaptive trapezoidal rule to integrate the normalization constant to a target
      61             :   // precision
      62         144 :   std::stack<std::pair<Real, Real>> S;
      63         144 :   S.emplace(0, 1);
      64         144 :   ADReal I = 0;
      65             :   unsigned int its = 0;
      66     2308128 :   while (!S.empty())
      67             :   {
      68     2307984 :     auto interval = S.top();
      69             :     S.pop();
      70             : 
      71     2307984 :     Real a = interval.first;
      72     2307984 :     Real b = interval.second;
      73     2307984 :     Real m = (a + b) / 2;
      74             : 
      75     2307984 :     Real I1 = (normalizationIntegrand(a) + normalizationIntegrand(b)) * (b - a) / 2;
      76             :     Real I2 =
      77     2307984 :         (normalizationIntegrand(a) + 2 * normalizationIntegrand(m) + normalizationIntegrand(b)) *
      78     2307984 :         (b - a) / 4;
      79             : 
      80     2307984 :     if (std::abs(I1 - I2) < 3 * (b - a) * _tolerance)
      81             :       I += I2;
      82             :     else
      83             :     {
      84             :       S.emplace(a, m);
      85             :       S.emplace(m, b);
      86     1153920 :       its++;
      87             :     }
      88             : 
      89     2307984 :     if (its >= _max_its)
      90           0 :       mooseError("Maximum number of iterations reached, but the crack geometric function still "
      91             :                  "hasn't converge.");
      92             :   }
      93         144 :   return I;
      94             : }
      95             : 
      96             : Real
      97    11539920 : CrackGeometricFunction::normalizationIntegrand(const ADReal & d)
      98             : {
      99    11539920 :   _func_params[_d_idx] = d;
     100    11539920 :   return 4.0 * std::sqrt(raw_value(evaluate(_func_F, _name)));
     101             : }
     102             : 
     103             : void
     104    18861177 : CrackGeometricFunction::computeQpProperties()
     105             : {
     106    18861177 :   _xi[_qp] = _xi_0;
     107    18861177 :   _c0[_qp] = _c0_0;
     108             : 
     109    18861177 :   CustomParsedFunctionBase::computeQpProperties();
     110    18861177 : }