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 62 : KKSPhaseConcentrationMaterial::validParams()
17 : {
18 62 : InputParameters params = DerivativeMaterialInterface<Material>::validParams();
19 62 : 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 124 : params.addRequiredCoupledVar("global_cs", "The interpolated concentrations c, b, etc.");
24 124 : params.addRequiredParam<MaterialPropertyName>("h_name", "Switching function h(eta).");
25 124 : 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 124 : params.addRequiredParam<std::vector<Real>>("ci_IC",
30 : "Initial values of ci in the same order as ci_names.");
31 124 : params.addRequiredParam<MaterialName>("fa_name", "Fa material object.");
32 124 : params.addRequiredParam<MaterialName>("fb_name", "Fb material object.");
33 124 : params.addParam<MaterialPropertyName>(
34 : "nested_iterations",
35 : "The output number of nested Newton iterations at each quadrature point.");
36 124 : params.addCoupledVar("args", "The coupled variables of Fa and Fb.");
37 124 : params.addParam<bool>(
38 124 : "damped_Newton", false, "Whether or not to use the damped Newton's method.");
39 124 : params.addParam<MaterialName>("conditions",
40 : "C",
41 : "Material property that checks bounds and conditions on the "
42 : "material properties being solved for.");
43 62 : params += NestedSolve::validParams();
44 62 : return params;
45 0 : }
46 :
47 48 : KKSPhaseConcentrationMaterial::KKSPhaseConcentrationMaterial(const InputParameters & parameters)
48 : : DerivativeMaterialInterface<Material>(parameters),
49 48 : _prop_c(coupledValues("global_cs")),
50 48 : _num_c(coupledComponents("global_cs")),
51 96 : _prop_h(getMaterialProperty<Real>("h_name")),
52 96 : _ci_names(getParam<std::vector<MaterialPropertyName>>("ci_names")),
53 48 : _prop_ci(_num_c * 2),
54 48 : _ci_old(_num_c * 2),
55 96 : _ci_IC(getParam<std::vector<Real>>("ci_IC")),
56 48 : _Fa_name(getParam<MaterialName>("fa_name")),
57 48 : _Fb_name(getParam<MaterialName>("fb_name")),
58 48 : _prop_Fi(2),
59 48 : _Fi_copy(2),
60 48 : _dFidci(_num_c * 2),
61 48 : _dFidci_copy(_num_c * 2),
62 48 : _d2Fidcidbi(2),
63 48 : _d2Fadc1db1_copy(_num_c),
64 48 : _args_names(coupledNames("args")),
65 48 : _n_args(coupledComponents("args")),
66 48 : _dFadarg(_n_args),
67 48 : _dFadarg_copy(_n_args),
68 48 : _dFbdarg(_n_args),
69 48 : _dFbdarg_copy(_n_args),
70 48 : _d2Fadcadarg(_n_args),
71 48 : _d2Fadcadarg_copy(_n_args),
72 48 : _iter(declareProperty<Real>("nested_iterations")),
73 96 : _abs_tol(getParam<Real>("absolute_tolerance")),
74 96 : _rel_tol(getParam<Real>("relative_tolerance")),
75 96 : _damped_newton(getParam<bool>("damped_Newton")),
76 48 : _condition_name(getParam<MaterialName>("conditions")),
77 96 : _nested_solve(NestedSolve(parameters))
78 :
79 : {
80 : // phase concentrations
81 144 : for (const auto m : make_range(_num_c * 2))
82 : {
83 96 : _prop_ci[m] = &declareProperty<Real>(_ci_names[m]);
84 96 : _ci_old[m] = &getMaterialPropertyOld<Real>(_ci_names[m]);
85 : }
86 :
87 : // free energies
88 48 : _prop_Fi[0] = &getMaterialPropertyByName<Real>(_Fa_name);
89 48 : _prop_Fi[1] = &getMaterialPropertyByName<Real>(_Fb_name);
90 :
91 : // declare _fi_copy to be passed to the kernels
92 48 : _Fi_copy[0] = &declareProperty<Real>("cp" + _Fa_name);
93 48 : _Fi_copy[1] = &declareProperty<Real>("cp" + _Fb_name);
94 :
95 : // derivative of free energy wrt phase concentrations
96 96 : for (const auto m : make_range(_num_c))
97 : {
98 48 : _dFidci[m * 2] = &getMaterialPropertyDerivative<Real>(_Fa_name, _ci_names[m * 2]);
99 48 : _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 48 : _dFidci_copy[m * 2] = &declarePropertyDerivative<Real>("cp" + _Fa_name, _ci_names[m * 2]);
103 48 : _dFidci_copy[m * 2 + 1] =
104 96 : &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 144 : for (const auto m : make_range(2))
111 : {
112 96 : _d2Fidcidbi[m].resize(_num_c);
113 192 : for (const auto n : make_range(_num_c))
114 : {
115 96 : _d2Fidcidbi[m][n].resize(_num_c);
116 192 : for (const auto l : make_range(_num_c))
117 : {
118 96 : if (m == 0)
119 48 : _d2Fidcidbi[0][n][l] =
120 96 : &getMaterialPropertyDerivative<Real>(_Fa_name, _ci_names[n * 2], _ci_names[l * 2]);
121 :
122 : else
123 48 : _d2Fidcidbi[1][n][l] = &getMaterialPropertyDerivative<Real>(
124 48 : _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 96 : for (const auto m : make_range(_num_c))
131 : {
132 48 : _d2Fadc1db1_copy[m].resize(_num_c);
133 96 : for (const auto n : make_range(_num_c))
134 48 : _d2Fadc1db1_copy[m][n] =
135 96 : &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 48 : 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 48 : 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 48 : if (_damped_newton)
163 24 : _C = &getMaterialPropertyByName<Real>(_condition_name);
164 : else
165 24 : _C = nullptr;
166 48 : }
167 :
168 : void
169 18000 : KKSPhaseConcentrationMaterial::initQpStatefulProperties()
170 : {
171 54000 : for (const auto m : make_range(_num_c * 2))
172 36000 : (*_prop_ci[m])[_qp] = _ci_IC[m];
173 18000 : }
174 :
175 : void
176 48 : KKSPhaseConcentrationMaterial::initialSetup()
177 : {
178 48 : _Fa = &getMaterial("fa_name");
179 48 : _Fb = &getMaterial("fb_name");
180 48 : if (_damped_newton)
181 24 : _condition = &getMaterialByName(_condition_name);
182 48 : }
183 :
184 : void
185 1501200 : KKSPhaseConcentrationMaterial::computeQpProperties()
186 : {
187 : // parameters for nested Newton iteration
188 1501200 : NestedSolve::Value<> solution(_num_c * 2);
189 :
190 4503600 : for (unsigned int m = 0; m < _num_c * 2; ++m)
191 3002400 : solution(m) = (*_ci_old[m])[_qp];
192 :
193 1501200 : _nested_solve.setAbsoluteTolerance(_abs_tol);
194 1501200 : _nested_solve.setRelativeTolerance(_rel_tol);
195 :
196 4778004 : auto compute = [&](const NestedSolve::Value<> & guess,
197 : NestedSolve::Value<> & residual,
198 : NestedSolve::Jacobian<> & jacobian)
199 : {
200 14334012 : for (const auto m : make_range(_num_c * 2))
201 9556008 : (*_prop_ci[m])[_qp] = guess(m);
202 :
203 4778004 : _Fa->computePropertiesAtQp(_qp);
204 4778004 : _Fb->computePropertiesAtQp(_qp);
205 :
206 : // assign residual functions
207 9556008 : for (const auto m : make_range(_num_c))
208 : {
209 4778004 : residual(m * 2) = (*_dFidci[m * 2])[_qp] - (*_dFidci[m * 2 + 1])[_qp];
210 4778004 : residual(m * 2 + 1) = (1 - _prop_h[_qp]) * (*_prop_ci[m * 2])[_qp] +
211 4778004 : _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 9556008 : for (const auto m : make_range(_num_c))
218 : {
219 9556008 : for (const auto n : make_range(_num_c))
220 : {
221 : // equal chemical potential derivative equations
222 4778004 : jacobian(m * 2, n * 2) = (*_d2Fidcidbi[0][m][n])[_qp];
223 4778004 : jacobian(m * 2, n * 2 + 1) = -(*_d2Fidcidbi[1][m][n])[_qp];
224 : }
225 : // concentration conservation derivative equations
226 4778004 : jacobian(m * 2 + 1, m * 2) = 1 - _prop_h[_qp];
227 4778004 : jacobian(m * 2 + 1, m * 2 + 1) = _prop_h[_qp];
228 : }
229 6279204 : };
230 18467212 : auto computeCondition = [&](const NestedSolve::Value<> & guess) -> Real
231 : {
232 55401636 : for (const auto m : make_range(_num_c * 2))
233 36934424 : (*_prop_ci[m])[_qp] = guess(m);
234 18467212 : _condition->computePropertiesAtQp(_qp);
235 18467212 : return ((*_C)[_qp]);
236 1501200 : };
237 : // choose between Newton or damped Newton's method
238 1501200 : if (!_damped_newton)
239 753300 : _nested_solve.nonlinear(solution, compute);
240 : else
241 747900 : _nested_solve.nonlinearDamped(solution, compute, computeCondition);
242 1501200 : _iter[_qp] = _nested_solve.getIterations();
243 :
244 1501200 : if (_nested_solve.getState() == NestedSolve::State::NOT_CONVERGED)
245 0 : mooseException("Nested Newton iteration did not converge.");
246 :
247 : // assign solution to ci
248 4503600 : for (const auto m : make_range(_num_c * 2))
249 3002400 : (*_prop_ci[m])[_qp] = solution[m];
250 :
251 : // assign to the copied parameters to be used in kernels
252 4503600 : for (const auto m : make_range(2))
253 3002400 : (*_Fi_copy[m])[_qp] = (*_prop_Fi[m])[_qp];
254 :
255 4503600 : for (const auto m : make_range(_num_c * 2))
256 3002400 : (*_dFidci_copy[m])[_qp] = (*_dFidci[m])[_qp];
257 :
258 3002400 : for (const auto m : make_range(_num_c))
259 : {
260 3002400 : for (const auto n : make_range(_num_c))
261 1501200 : (*_d2Fadc1db1_copy[m][n])[_qp] = (*_d2Fidcidbi[0][m][n])[_qp];
262 : }
263 :
264 1501200 : 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 1501200 : 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 1501200 : }
|
