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 "KKSPhaseConcentrationMaterial.h"
11 : #include "MatrixTools.h"
12 :
13 : registerMooseObject("PhaseFieldApp", KKSPhaseConcentrationMaterial);
14 :
15 : InputParameters
16 94 : KKSPhaseConcentrationMaterial::validParams()
17 : {
18 94 : InputParameters params = DerivativeMaterialInterface<Material>::validParams();
19 94 : params.addClassDescription("Computes the KKS phase concentrations by using nested Newton "
20 : "iteration to solve the equal chemical "
21 : "potential and concentration conservation equations. This class is "
22 : "intended to be used with KKSPhaseConcentrationDerivatives.");
23 188 : params.addRequiredCoupledVar("global_cs", "The interpolated concentrations c, b, etc.");
24 188 : params.addRequiredParam<MaterialPropertyName>("h_name", "Switching function h(eta).");
25 188 : params.addRequiredParam<std::vector<MaterialPropertyName>>(
26 : "ci_names",
27 : "Phase concentrations. The order must match Fa, Fb, and global_cs, for example, c1, c2, b1, "
28 : "b2, etc.");
29 188 : params.addRequiredParam<std::vector<Real>>("ci_IC",
30 : "Initial values of ci in the same order as ci_names.");
31 188 : params.addRequiredParam<MaterialName>("fa_name", "Fa material object.");
32 188 : params.addRequiredParam<MaterialName>("fb_name", "Fb material object.");
33 188 : params.addParam<MaterialPropertyName>(
34 : "nested_iterations",
35 : "The output number of nested Newton iterations at each quadrature point.");
36 188 : params.addCoupledVar("args", "The coupled variables of Fa and Fb.");
37 188 : params.addParam<bool>(
38 188 : "damped_Newton", false, "Whether or not to use the damped Newton's method.");
39 188 : params.addParam<MaterialName>("conditions",
40 : "C",
41 : "Material property that checks bounds and conditions on the "
42 : "material properties being solved for.");
43 94 : params += NestedSolve::validParams();
44 94 : return params;
45 0 : }
46 :
47 72 : KKSPhaseConcentrationMaterial::KKSPhaseConcentrationMaterial(const InputParameters & parameters)
48 : : DerivativeMaterialInterface<Material>(parameters),
49 72 : _prop_c(coupledValues("global_cs")),
50 72 : _num_c(coupledComponents("global_cs")),
51 144 : _prop_h(getMaterialProperty<Real>("h_name")),
52 144 : _ci_names(getParam<std::vector<MaterialPropertyName>>("ci_names")),
53 72 : _prop_ci(_num_c * 2),
54 72 : _ci_old(_num_c * 2),
55 144 : _ci_IC(getParam<std::vector<Real>>("ci_IC")),
56 72 : _Fa_name(getParam<MaterialName>("fa_name")),
57 72 : _Fb_name(getParam<MaterialName>("fb_name")),
58 72 : _prop_Fi(2),
59 72 : _Fi_copy(2),
60 72 : _dFidci(_num_c * 2),
61 72 : _dFidci_copy(_num_c * 2),
62 72 : _d2Fidcidbi(2),
63 72 : _d2Fadc1db1_copy(_num_c),
64 72 : _args_names(coupledNames("args")),
65 72 : _n_args(coupledComponents("args")),
66 72 : _dFadarg(_n_args),
67 72 : _dFadarg_copy(_n_args),
68 72 : _dFbdarg(_n_args),
69 72 : _dFbdarg_copy(_n_args),
70 72 : _d2Fadcadarg(_n_args),
71 72 : _d2Fadcadarg_copy(_n_args),
72 72 : _iter(declareProperty<Real>("nested_iterations")),
73 144 : _abs_tol(getParam<Real>("absolute_tolerance")),
74 144 : _rel_tol(getParam<Real>("relative_tolerance")),
75 144 : _damped_newton(getParam<bool>("damped_Newton")),
76 72 : _condition_name(getParam<MaterialName>("conditions")),
77 144 : _nested_solve(NestedSolve(parameters))
78 :
79 : {
80 : // phase concentrations
81 216 : for (const auto m : make_range(_num_c * 2))
82 : {
83 144 : _prop_ci[m] = &declareProperty<Real>(_ci_names[m]);
84 144 : _ci_old[m] = &getMaterialPropertyOld<Real>(_ci_names[m]);
85 : }
86 :
87 : // free energies
88 72 : _prop_Fi[0] = &getMaterialPropertyByName<Real>(_Fa_name);
89 72 : _prop_Fi[1] = &getMaterialPropertyByName<Real>(_Fb_name);
90 :
91 : // declare _fi_copy to be passed to the kernels
92 72 : _Fi_copy[0] = &declareProperty<Real>("cp" + _Fa_name);
93 72 : _Fi_copy[1] = &declareProperty<Real>("cp" + _Fb_name);
94 :
95 : // derivative of free energy wrt phase concentrations
96 144 : for (const auto m : make_range(_num_c))
97 : {
98 72 : _dFidci[m * 2] = &getMaterialPropertyDerivative<Real>(_Fa_name, _ci_names[m * 2]);
99 72 : _dFidci[m * 2 + 1] = &getMaterialPropertyDerivative<Real>(_Fb_name, _ci_names[m * 2 + 1]);
100 :
101 : // declare _dFidci_copy to be passed to the kernels
102 72 : _dFidci_copy[m * 2] = &declarePropertyDerivative<Real>("cp" + _Fa_name, _ci_names[m * 2]);
103 72 : _dFidci_copy[m * 2 + 1] =
104 144 : &declarePropertyDerivative<Real>("cp" + _Fb_name, _ci_names[m * 2 + 1]);
105 : }
106 :
107 : // Second derivative of free energy wrt phase concentrations for use in this material. In
108 : // _d2Fidcidbi[m][n][l], m is phase index of Fi, n is the species index of ci, l is the species
109 : // index of bi.
110 216 : for (const auto m : make_range(2))
111 : {
112 144 : _d2Fidcidbi[m].resize(_num_c);
113 288 : for (const auto n : make_range(_num_c))
114 : {
115 144 : _d2Fidcidbi[m][n].resize(_num_c);
116 288 : for (const auto l : make_range(_num_c))
117 : {
118 144 : if (m == 0)
119 72 : _d2Fidcidbi[0][n][l] =
120 144 : &getMaterialPropertyDerivative<Real>(_Fa_name, _ci_names[n * 2], _ci_names[l * 2]);
121 :
122 : else
123 72 : _d2Fidcidbi[1][n][l] = &getMaterialPropertyDerivative<Real>(
124 72 : _Fb_name, _ci_names[n * 2 + 1], _ci_names[l * 2 + 1]);
125 : }
126 : }
127 : }
128 :
129 : // _d2Fadc1db1_copy (2D symmetric matrix), to be passed to kernels
130 144 : for (const auto m : make_range(_num_c))
131 : {
132 72 : _d2Fadc1db1_copy[m].resize(_num_c);
133 144 : for (const auto n : make_range(_num_c))
134 72 : _d2Fadc1db1_copy[m][n] =
135 144 : &declarePropertyDerivative<Real>("cp" + _Fa_name, _ci_names[m * 2], _ci_names[n * 2]);
136 : }
137 :
138 : // partial derivative of Fa and Fb wrt coupled variables, to be passed to kernels
139 72 : for (const auto m : make_range(_n_args))
140 : {
141 0 : _dFadarg[m] = &getMaterialPropertyDerivative<Real>(_Fa_name, _args_names[m]);
142 0 : _dFadarg_copy[m] = &declarePropertyDerivative<Real>("cp" + _Fa_name, _args_names[m]);
143 0 : _dFbdarg[m] = &getMaterialPropertyDerivative<Real>(_Fb_name, _args_names[m]);
144 0 : _dFbdarg_copy[m] = &declarePropertyDerivative<Real>("cp" + _Fb_name, _args_names[m]);
145 : }
146 :
147 : // second partial derivatives of Fa wrt ca and another coupled variable, to be passed to
148 : // kernels
149 72 : for (const auto m : make_range(_n_args))
150 : {
151 0 : _d2Fadcadarg[m].resize(_num_c);
152 0 : _d2Fadcadarg_copy[m].resize(_num_c);
153 0 : for (const auto n : make_range(_num_c))
154 : {
155 0 : _d2Fadcadarg[m][n] =
156 0 : &getMaterialPropertyDerivative<Real>(_Fa_name, _ci_names[n * 2], _args_names[m]);
157 0 : _d2Fadcadarg_copy[m][n] =
158 0 : &declarePropertyDerivative<Real>("cp" + _Fa_name, _ci_names[n * 2], _args_names[m]);
159 : }
160 : }
161 :
162 72 : if (_damped_newton)
163 36 : _C = &getMaterialPropertyByName<Real>(_condition_name);
164 : else
165 36 : _C = nullptr;
166 72 : }
167 :
168 : void
169 21600 : KKSPhaseConcentrationMaterial::initQpStatefulProperties()
170 : {
171 64800 : for (const auto m : make_range(_num_c * 2))
172 43200 : (*_prop_ci[m])[_qp] = _ci_IC[m];
173 21600 : }
174 :
175 : void
176 72 : KKSPhaseConcentrationMaterial::initialSetup()
177 : {
178 72 : _Fa = &getMaterial("fa_name");
179 72 : _Fb = &getMaterial("fb_name");
180 72 : if (_damped_newton)
181 36 : _condition = &getMaterialByName(_condition_name);
182 72 : }
183 :
184 : void
185 1884600 : KKSPhaseConcentrationMaterial::computeQpProperties()
186 : {
187 : // parameters for nested Newton iteration
188 1884600 : NestedSolve::Value<> solution(_num_c * 2);
189 :
190 5653800 : for (unsigned int m = 0; m < _num_c * 2; ++m)
191 3769200 : solution(m) = (*_ci_old[m])[_qp];
192 :
193 1884600 : _nested_solve.setAbsoluteTolerance(_abs_tol);
194 1884600 : _nested_solve.setRelativeTolerance(_rel_tol);
195 :
196 6047394 : auto compute = [&](const NestedSolve::Value<> & guess,
197 : NestedSolve::Value<> & residual,
198 : NestedSolve::Jacobian<> & jacobian)
199 : {
200 18142182 : for (const auto m : make_range(_num_c * 2))
201 12094788 : (*_prop_ci[m])[_qp] = guess(m);
202 :
203 6047394 : _Fa->computePropertiesAtQp(_qp);
204 6047394 : _Fb->computePropertiesAtQp(_qp);
205 :
206 : // assign residual functions
207 12094788 : for (const auto m : make_range(_num_c))
208 : {
209 6047394 : residual(m * 2) = (*_dFidci[m * 2])[_qp] - (*_dFidci[m * 2 + 1])[_qp];
210 6047394 : residual(m * 2 + 1) = (1 - _prop_h[_qp]) * (*_prop_ci[m * 2])[_qp] +
211 6047394 : _prop_h[_qp] * (*_prop_ci[m * 2 + 1])[_qp] - (*_prop_c[m])[_qp];
212 : }
213 :
214 : jacobian.setZero();
215 :
216 : // fill in the non-zero elements in jacobian
217 12094788 : for (const auto m : make_range(_num_c))
218 : {
219 12094788 : for (const auto n : make_range(_num_c))
220 : {
221 : // equal chemical potential derivative equations
222 6047394 : jacobian(m * 2, n * 2) = (*_d2Fidcidbi[0][m][n])[_qp];
223 6047394 : jacobian(m * 2, n * 2 + 1) = -(*_d2Fidcidbi[1][m][n])[_qp];
224 : }
225 : // concentration conservation derivative equations
226 6047394 : jacobian(m * 2 + 1, m * 2) = 1 - _prop_h[_qp];
227 6047394 : jacobian(m * 2 + 1, m * 2 + 1) = _prop_h[_qp];
228 : }
229 7931994 : };
230 23620316 : auto computeCondition = [&](const NestedSolve::Value<> & guess) -> Real
231 : {
232 70860948 : for (const auto m : make_range(_num_c * 2))
233 47240632 : (*_prop_ci[m])[_qp] = guess(m);
234 23620316 : _condition->computePropertiesAtQp(_qp);
235 23620316 : return ((*_C)[_qp]);
236 1884600 : };
237 : // choose between Newton or damped Newton's method
238 1884600 : if (!_damped_newton)
239 934200 : _nested_solve.nonlinear(solution, compute);
240 : else
241 950400 : _nested_solve.nonlinearDamped(solution, compute, computeCondition);
242 1884600 : _iter[_qp] = _nested_solve.getIterations();
243 :
244 1884600 : if (_nested_solve.getState() == NestedSolve::State::NOT_CONVERGED)
245 0 : mooseException("Nested Newton iteration did not converge.");
246 :
247 : // assign solution to ci
248 5653800 : for (const auto m : make_range(_num_c * 2))
249 3769200 : (*_prop_ci[m])[_qp] = solution[m];
250 :
251 : // assign to the copied parameters to be used in kernels
252 5653800 : for (const auto m : make_range(2))
253 3769200 : (*_Fi_copy[m])[_qp] = (*_prop_Fi[m])[_qp];
254 :
255 5653800 : for (const auto m : make_range(_num_c * 2))
256 3769200 : (*_dFidci_copy[m])[_qp] = (*_dFidci[m])[_qp];
257 :
258 3769200 : for (const auto m : make_range(_num_c))
259 : {
260 3769200 : for (const auto n : make_range(_num_c))
261 1884600 : (*_d2Fadc1db1_copy[m][n])[_qp] = (*_d2Fidcidbi[0][m][n])[_qp];
262 : }
263 :
264 1884600 : for (const auto m : make_range(_n_args))
265 : {
266 0 : (*_dFadarg_copy[m])[_qp] = (*_dFadarg[m])[_qp];
267 0 : (*_dFbdarg_copy[m])[_qp] = (*_dFbdarg[m])[_qp];
268 : }
269 :
270 1884600 : for (const auto m : make_range(_n_args))
271 : {
272 0 : for (const auto n : make_range(_num_c))
273 0 : (*_d2Fadcadarg_copy[m][n])[_qp] = (*_d2Fadcadarg[m][n])[_qp];
274 : }
275 1884600 : }
|
