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 "KKSPhaseConcentrationMultiPhaseMaterial.h"
11 : #include "MatrixTools.h"
12 :
13 : registerMooseObject("PhaseFieldApp", KKSPhaseConcentrationMultiPhaseMaterial);
14 :
15 : InputParameters
16 94 : KKSPhaseConcentrationMultiPhaseMaterial::validParams()
17 : {
18 94 : InputParameters params = DerivativeMaterialInterface<Material>::validParams();
19 94 : params.addClassDescription(
20 : "Computes the KKS phase concentrations by using a nested Newton iteration "
21 : "to solve the equal chemical potential and concentration conservation equations for "
22 : "multiphase systems. This class is intented to be used with "
23 : "KKSPhaseConcentrationMultiPhaseDerivatives.");
24 188 : params.addRequiredCoupledVar("global_cs", "The interpolated concentrations c, b, etc.");
25 188 : params.addRequiredCoupledVar("all_etas", "Order parameters.");
26 188 : params.addRequiredParam<std::vector<MaterialPropertyName>>(
27 : "hj_names", "Switching functions in the same order as all_etas.");
28 188 : params.addRequiredParam<std::vector<MaterialName>>(
29 : "Fj_names", "Free energy material objects in the same order as all_etas.");
30 188 : params.addRequiredParam<std::vector<MaterialPropertyName>>(
31 : "ci_names",
32 : "Phase concentrations. They must have the same order as Fj_names and global_cs, for "
33 : "example, c1, c2, b1, b2.");
34 188 : params.addRequiredParam<std::vector<Real>>("ci_IC",
35 : "Initial values of ci in the same order of ci_names");
36 188 : params.addParam<MaterialPropertyName>(
37 : "nested_iterations",
38 : "The output number of nested Newton iterations at each quadrature point.");
39 188 : params.addCoupledVar("args", "The coupled variables of free energies.");
40 188 : params.addParam<bool>(
41 188 : "damped_Newton", false, "Whether or not to use the damped Newton's method.");
42 188 : params.addParam<MaterialName>("conditions",
43 : "C",
44 : "Material property that checks bounds and conditions on the "
45 : "material properties being solved for.");
46 94 : params += NestedSolve::validParams();
47 94 : return params;
48 0 : }
49 :
50 72 : KKSPhaseConcentrationMultiPhaseMaterial::KKSPhaseConcentrationMultiPhaseMaterial(
51 72 : const InputParameters & parameters)
52 : : DerivativeMaterialInterface<Material>(parameters),
53 72 : _prop_c(coupledValues("global_cs")),
54 72 : _num_c(coupledComponents("global_cs")),
55 72 : _num_j(coupledComponents("all_etas")),
56 144 : _hj_names(getParam<std::vector<MaterialPropertyName>>("hj_names")),
57 72 : _prop_hj(_num_j),
58 144 : _Fj_names(getParam<std::vector<MaterialName>>("Fj_names")),
59 72 : _prop_Fi(_num_j),
60 144 : _ci_names(getParam<std::vector<MaterialPropertyName>>("ci_names")),
61 72 : _prop_ci(_num_c * _num_j),
62 72 : _ci_old(_num_c * _num_j),
63 144 : _ci_IC(getParam<std::vector<Real>>("ci_IC")),
64 72 : _dFidci(_num_j),
65 72 : _d2Fidcidbi(_num_j),
66 72 : _args_names(coupledNames("args")),
67 72 : _n_args(coupledComponents("args")),
68 72 : _dFidarg(_num_j),
69 72 : _d2F1dc1darg(_num_c),
70 72 : _iter(declareProperty<Real>("nested_iterations")),
71 144 : _abs_tol(getParam<Real>("absolute_tolerance")),
72 144 : _rel_tol(getParam<Real>("relative_tolerance")),
73 144 : _damped_newton(getParam<bool>("damped_Newton")),
74 72 : _condition_name(getParam<MaterialName>("conditions")),
75 144 : _nested_solve(NestedSolve(parameters))
76 :
77 : {
78 : // phase concentrations
79 288 : for (const auto m : make_range(_num_c * _num_j))
80 : {
81 216 : _ci_old[m] = &getMaterialPropertyOld<Real>(_ci_names[m]);
82 216 : _prop_ci[m] = &declareProperty<Real>(_ci_names[m]);
83 : }
84 :
85 : // free energies
86 288 : for (const auto m : make_range(_num_j))
87 216 : _prop_Fi[m] = &getMaterialPropertyByName<Real>(_Fj_names[m]);
88 :
89 : // free energy derivatives w.r.t. phase concentrations
90 288 : for (const auto m : make_range(_num_j))
91 : {
92 216 : _prop_hj[m] = &getMaterialPropertyByName<Real>(_hj_names[m]);
93 :
94 216 : _dFidci[m].resize(_num_c);
95 216 : _d2Fidcidbi[m].resize(_num_c);
96 :
97 432 : for (const auto n : make_range(_num_c))
98 : {
99 216 : _dFidci[m][n] = &getMaterialPropertyDerivative<Real>(_Fj_names[m], _ci_names[m + n * _num_j]);
100 :
101 216 : _d2Fidcidbi[m][n].resize(_num_c);
102 :
103 432 : for (unsigned int l = 0; l < _num_c; ++l)
104 : {
105 216 : _d2Fidcidbi[m][n][l] = &getMaterialPropertyDerivative<Real>(
106 216 : _Fj_names[m], _ci_names[m + n * _num_j], _ci_names[m + l * _num_j]);
107 : }
108 : }
109 : }
110 :
111 : // derivative of free energies wrt coupled variables
112 288 : for (const auto m : make_range(_num_j))
113 : {
114 216 : _dFidarg[m].resize(_n_args);
115 :
116 216 : for (const auto n : make_range(_n_args))
117 0 : _dFidarg[m][n] = &getMaterialPropertyDerivative<Real>(_Fj_names[m], _args_names[n]);
118 : }
119 :
120 : // second derivatives of F1 wrt c1 and other coupled variables
121 144 : for (const auto m : make_range(_num_c))
122 : {
123 72 : _d2F1dc1darg[m].resize(_n_args);
124 :
125 72 : for (const auto n : make_range(_n_args))
126 : {
127 0 : _d2F1dc1darg[m][n] =
128 0 : &getMaterialPropertyDerivative<Real>(_Fj_names[0], _ci_names[m * _num_j], _args_names[n]);
129 : }
130 : }
131 :
132 72 : if (_damped_newton)
133 36 : _C = &getMaterialPropertyByName<Real>(_condition_name);
134 : else
135 36 : _C = nullptr;
136 72 : }
137 :
138 : void
139 38400 : KKSPhaseConcentrationMultiPhaseMaterial::initQpStatefulProperties()
140 : {
141 153600 : for (const auto m : make_range(_num_c * _num_j))
142 115200 : (*_prop_ci[m])[_qp] = _ci_IC[m];
143 38400 : }
144 :
145 : void
146 72 : KKSPhaseConcentrationMultiPhaseMaterial::initialSetup()
147 : {
148 72 : _Fj_mat.resize(_num_j);
149 :
150 288 : for (unsigned int m = 0; m < _num_j; ++m)
151 216 : _Fj_mat[m] = &getMaterialByName(_Fj_names[m]);
152 72 : if (_damped_newton)
153 36 : _condition = &getMaterialByName(_condition_name);
154 72 : }
155 :
156 : void
157 2454400 : KKSPhaseConcentrationMultiPhaseMaterial::computeQpProperties()
158 : {
159 : // parameters for nested Newton iteration
160 2454400 : NestedSolve::Value<> solution(_num_c * _num_j);
161 :
162 : // initialize first guess using the solution from previous step
163 9817600 : for (const auto m : make_range(_num_c * _num_j))
164 7363200 : solution(m) = (*_ci_old[m])[_qp];
165 :
166 2454400 : _nested_solve.setAbsoluteTolerance(_abs_tol);
167 2454400 : _nested_solve.setRelativeTolerance(_rel_tol);
168 :
169 8180389 : auto compute = [&](const NestedSolve::Value<> & guess,
170 : NestedSolve::Value<> & residual,
171 : NestedSolve::Jacobian<> & jacobian)
172 : {
173 32721556 : for (const auto m : make_range(_num_c * _num_j))
174 24541167 : (*_prop_ci[m])[_qp] = guess(m);
175 :
176 32721556 : for (const auto m : make_range(_num_j))
177 24541167 : _Fj_mat[m]->computePropertiesAtQp(_qp);
178 :
179 : // assign residual functions
180 16360778 : for (const auto m : make_range(_num_c))
181 : {
182 24541167 : for (const auto n : make_range(_num_j - 1))
183 16360778 : residual(m * _num_j + n) = (*_dFidci[n][m])[_qp] - (*_dFidci[n + 1][m])[_qp];
184 :
185 8180389 : residual((m + 1) * _num_j - 1) = -(*_prop_c[m])[_qp];
186 :
187 32721556 : for (const auto l : make_range(_num_j))
188 24541167 : residual((m + 1) * _num_j - 1) += (*_prop_hj[l])[_qp] * (*_prop_ci[m * _num_j + l])[_qp];
189 : }
190 :
191 : jacobian.setZero();
192 :
193 : // fill in the non-zero terms in jacobian
194 16360778 : for (const auto m : make_range(_num_c))
195 : {
196 : // equal chemical potential derivative equations
197 24541167 : for (const auto n : make_range(_num_j - 1))
198 : {
199 32721556 : for (const auto l : make_range(_num_c))
200 : {
201 16360778 : jacobian(m * _num_j + n, n + l * _num_j) = (*_d2Fidcidbi[n][m][l])[_qp];
202 16360778 : jacobian(m * _num_j + n, n + l * _num_j + 1) = -(*_d2Fidcidbi[n + 1][m][l])[_qp];
203 : }
204 : }
205 :
206 : // concentration conservation derivative equations
207 32721556 : for (const auto n : make_range(_num_j))
208 24541167 : jacobian((m + 1) * _num_j - 1, m * _num_j + n) = (*_prop_hj[n])[_qp];
209 : }
210 10634789 : };
211 :
212 22344276 : auto computeCondition = [&](const NestedSolve::Value<> & guess) -> Real
213 : {
214 89377104 : for (const auto m : make_range(_num_c * _num_j))
215 67032828 : (*_prop_ci[m])[_qp] = guess(m);
216 22344276 : _condition->computePropertiesAtQp(_qp);
217 22344276 : return ((*_C)[_qp]);
218 2454400 : };
219 :
220 : // choose between Newton or damped Newton's method
221 2454400 : if (!_damped_newton)
222 1329600 : _nested_solve.nonlinear(solution, compute);
223 : else
224 1124800 : _nested_solve.nonlinearDamped(solution, compute, computeCondition);
225 :
226 2454400 : _iter[_qp] = _nested_solve.getIterations();
227 :
228 2454400 : if (_nested_solve.getState() == NestedSolve::State::NOT_CONVERGED)
229 0 : mooseException("Nested Newton iteration did not converge.");
230 :
231 : // assign solution to ci
232 9817600 : for (const auto m : make_range(_num_c * _num_j))
233 7363200 : (*_prop_ci[m])[_qp] = solution[m];
234 2454400 : }
|
