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 : }