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 340 : CrackGeometricFunction::validParams() 11 : { 12 340 : InputParameters params = CustomParsedFunctionBase::validParams(); 13 340 : 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 340 : params.set<std::string>("property_name") = "alpha"; 19 680 : params.addRequiredCoupledVar("phase_field", "The phase-field variable"); 20 : 21 340 : params.set<unsigned int>("derivative_order") = 1; 22 340 : params.suppressParameter<unsigned int>("derivative_order"); 23 : 24 680 : 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 680 : 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 680 : params.addParam<Real>( 35 680 : "tolerance", 1e-8, "Tolerance of the numerically computed normalization constant"); 36 680 : params.addParam<unsigned int>("maximum_iterations", 37 680 : 1e9, 38 : "Maximum number of iterations allowed for the numerical " 39 : "computation of the normalization constant"); 40 340 : return params; 41 0 : } 42 : 43 153 : CrackGeometricFunction::CrackGeometricFunction(const InputParameters & parameters) 44 : : CustomParsedFunctionBase(parameters), 45 306 : _d_idx(argIndex(coupled("phase_field"))), 46 306 : _xi(declareADProperty<Real>(getParam<MaterialPropertyName>("initial_derivative"))), 47 306 : _c0(declareADProperty<Real>(getParam<MaterialPropertyName>("normalization_constant"))), 48 306 : _tolerance(getParam<Real>("tolerance")), 49 459 : _max_its(getParam<unsigned int>("maximum_iterations")) 50 : { 51 : // set d to 0 and evaluate the first derivative of the crack geometric function 52 153 : _func_params[_d_idx] = 0; 53 153 : _xi_0 = evaluate(_derivatives[0]._F, _name); 54 153 : _c0_0 = computeNormalizationConstant(); 55 153 : } 56 : 57 : ADReal 58 153 : CrackGeometricFunction::computeNormalizationConstant() 59 : { 60 : // We use an adaptive trapezoidal rule to integrate the normalization constant to a target 61 : // precision 62 153 : std::stack<std::pair<Real, Real>> S; 63 153 : S.emplace(0, 1); 64 153 : ADReal I = 0; 65 : unsigned int its = 0; 66 2457078 : while (!S.empty()) 67 : { 68 2456925 : auto interval = S.top(); 69 : S.pop(); 70 : 71 2456925 : Real a = interval.first; 72 2456925 : Real b = interval.second; 73 2456925 : Real m = (a + b) / 2; 74 : 75 2456925 : Real I1 = (normalizationIntegrand(a) + normalizationIntegrand(b)) * (b - a) / 2; 76 : Real I2 = 77 2456925 : (normalizationIntegrand(a) + 2 * normalizationIntegrand(m) + normalizationIntegrand(b)) * 78 2456925 : (b - a) / 4; 79 : 80 2456925 : 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 1228386 : its++; 87 : } 88 : 89 2456925 : if (its >= _max_its) 90 0 : mooseError("Maximum number of iterations reached, but the crack geometric function still " 91 : "hasn't converge."); 92 : } 93 153 : return I; 94 : } 95 : 96 : Real 97 12284625 : CrackGeometricFunction::normalizationIntegrand(const ADReal & d) 98 : { 99 12284625 : _func_params[_d_idx] = d; 100 12284625 : return 4.0 * std::sqrt(raw_value(evaluate(_func_F, _name))); 101 : } 102 : 103 : void 104 19387953 : CrackGeometricFunction::computeQpProperties() 105 : { 106 19387953 : _xi[_qp] = _xi_0; 107 19387953 : _c0[_qp] = _c0_0; 108 : 109 19387953 : CustomParsedFunctionBase::computeQpProperties(); 110 19387953 : }