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 "KKSPhaseConcentrationMultiPhaseDerivatives.h"
11 : #include "MatrixTools.h"
12 :
13 : registerMooseObject("PhaseFieldApp", KKSPhaseConcentrationMultiPhaseDerivatives);
14 :
15 : InputParameters
16 94 : KKSPhaseConcentrationMultiPhaseDerivatives::validParams()
17 : {
18 94 : InputParameters params = DerivativeMaterialInterface<Material>::validParams();
19 94 : params.addClassDescription(
20 : "Computes the KKS phase concentration derivatives wrt global concentrations and order "
21 : "parameters, which are used for the chain rule in the KKS kernels. This class is intended to "
22 : "be used with KKSPhaseConcentrationMultiPhaseMaterial.");
23 188 : params.addRequiredCoupledVar("global_cs", "The interpolated concentrations c, b, etc");
24 188 : params.addRequiredCoupledVar("all_etas", "Order parameters.");
25 188 : params.addRequiredParam<std::vector<MaterialPropertyName>>(
26 : "ci_names",
27 : "Phase concentrations. They must have the same order as Fj_names and global_cs, for "
28 : "example, c1, c2, b1, b2.");
29 188 : params.addRequiredParam<std::vector<MaterialName>>(
30 : "Fj_names", "Free energy material objects in the same order as all_etas.");
31 188 : params.addRequiredParam<std::vector<MaterialPropertyName>>(
32 : "hj_names", "witching functions in the same order as all_etas.");
33 94 : return params;
34 0 : }
35 :
36 72 : KKSPhaseConcentrationMultiPhaseDerivatives::KKSPhaseConcentrationMultiPhaseDerivatives(
37 72 : const InputParameters & parameters)
38 : : DerivativeMaterialInterface<Material>(parameters),
39 72 : _num_c(coupledComponents("global_cs")),
40 72 : _c_names(coupledNames("global_cs")),
41 72 : _eta_names(coupledNames("all_etas")),
42 72 : _num_j(coupledComponents("all_etas")),
43 72 : _prop_ci(_num_c * _num_j),
44 144 : _ci_names(getParam<std::vector<MaterialPropertyName>>("ci_names")),
45 72 : _dcidetaj(_num_c),
46 72 : _dcidb(_num_c),
47 144 : _Fj_names(getParam<std::vector<MaterialName>>("Fj_names")),
48 72 : _d2Fidcidbi(_num_j),
49 144 : _hj_names(getParam<std::vector<MaterialPropertyName>>("hj_names")),
50 72 : _prop_hj(_num_j),
51 144 : _dhjdetai(_num_j)
52 :
53 : {
54 288 : for (const auto m : make_range(_num_c * _num_j))
55 216 : _prop_ci[m] = &getMaterialPropertyByName<Real>(_ci_names[m]);
56 :
57 144 : for (const auto m : make_range(_num_c))
58 : {
59 72 : _dcidb[m].resize(_num_j);
60 72 : _dcidetaj[m].resize(_num_j);
61 :
62 288 : for (const auto n : make_range(_num_j))
63 : {
64 216 : _dcidb[m][n].resize(_num_c);
65 216 : _dcidetaj[m][n].resize(_num_j);
66 :
67 : // Derivative of phase concentration wrt global concentration. In _dcidb[m][n][l], m is the
68 : // species index of ci, n is the phase index of ci, and l is the species index of b
69 432 : for (const auto l : make_range(_num_c))
70 216 : _dcidb[m][n][l] = &declarePropertyDerivative<Real>(_ci_names[n + m * _num_j], _c_names[l]);
71 :
72 : // Derivative of phase concentration wrt eta. In _dcidetaj[m][n][l], m is the species index
73 : // of ci, n is the phase index of ci, and l is the phase of etaj
74 864 : for (const auto l : make_range(_num_j))
75 648 : _dcidetaj[m][n][l] =
76 1296 : &declarePropertyDerivative<Real>(_ci_names[n + m * _num_j], _eta_names[l]);
77 : }
78 : }
79 :
80 : // Second derivative of free energy wrt phase concentrations for use in this material. In
81 : // _d2Fidcidbi[m][n][l], m is phase index of Fi, n is the species index of ci, l is the species
82 : // index of bi.
83 288 : for (const auto m : make_range(_num_j))
84 : {
85 216 : _d2Fidcidbi[m].resize(_num_c);
86 :
87 432 : for (const auto n : make_range(_num_c))
88 : {
89 216 : _d2Fidcidbi[m][n].resize(_num_c);
90 :
91 432 : for (const auto l : make_range(_num_c))
92 216 : _d2Fidcidbi[m][n][l] = &getMaterialPropertyDerivative<Real>(
93 216 : _Fj_names[m], _ci_names[m + n * _num_j], _ci_names[m + l * _num_j]);
94 : }
95 : }
96 :
97 288 : for (const auto m : make_range(_num_j))
98 : {
99 216 : _prop_hj[m] = &getMaterialPropertyByName<Real>(_hj_names[m]);
100 :
101 216 : _dhjdetai[m].resize(_num_j);
102 :
103 864 : for (const auto n : make_range(_num_j))
104 648 : _dhjdetai[m][n] = &getMaterialPropertyDerivative<Real>(_hj_names[m], _eta_names[n]);
105 : }
106 72 : }
107 :
108 : void
109 2454400 : KKSPhaseConcentrationMultiPhaseDerivatives::computeQpProperties()
110 : {
111 : // declare Jacobian matrix A
112 2454400 : Eigen::MatrixXd A(_num_c * _num_j, _num_c * _num_j);
113 :
114 : // initialize all elements in A to be zero
115 : A.setZero();
116 :
117 : // fill in the non-zero elements in A
118 4908800 : for (const auto m : make_range(_num_c))
119 : {
120 : // equal chemical potential derivative equations
121 7363200 : for (const auto n : make_range(_num_j - 1))
122 : {
123 9817600 : for (const auto l : make_range(_num_c))
124 : {
125 4908800 : A(m * _num_j + n, n + l * _num_j) = (*_d2Fidcidbi[n][m][l])[_qp];
126 4908800 : A(m * _num_j + n, n + l * _num_j + 1) = -(*_d2Fidcidbi[n + 1][m][l])[_qp];
127 : }
128 : }
129 :
130 : // concentration conservation derivative equations
131 9817600 : for (const auto n : make_range(_num_j))
132 7363200 : A((m + 1) * _num_j - 1, m * _num_j + n) = (*_prop_hj[n])[_qp];
133 : }
134 :
135 2454400 : A = A.inverse();
136 :
137 : // solve linear system of constraint derivatives wrt b for computing dcidb loop through
138 : // derivatives wrt the ith component; they have the same A, but different k_c
139 4908800 : for (const auto i : make_range(_num_c))
140 : {
141 2454400 : std::vector<Real> k_c(_num_j * _num_c);
142 2454400 : std::vector<Real> x_c(_num_j * _num_c);
143 :
144 : // assign non-zero elements in k_c
145 2454400 : k_c[i * _num_j + _num_j - 1] = 1;
146 :
147 : // compute x_c
148 9817600 : for (const auto m : make_range(_num_j * _num_c))
149 : {
150 29452800 : for (const auto n : make_range(_num_j * _num_c))
151 22089600 : x_c[m] += A(m, n) * k_c[n];
152 : }
153 :
154 : // assign the values in x_c to _dcidb
155 4908800 : for (const auto m : make_range(_num_c))
156 : {
157 9817600 : for (const auto n : make_range(_num_j))
158 7363200 : (*_dcidb[m][n][i])[_qp] = x_c[m * _num_j + n];
159 : }
160 2454400 : }
161 :
162 : // solve linear system of constraint derivatives wrt eta for computing dcidetaj use the same
163 : // linear matrix as computing dcidb
164 9817600 : for (const auto i : make_range(_num_j))
165 : {
166 7363200 : std::vector<Real> k_eta(_num_j * _num_c);
167 7363200 : std::vector<Real> x_eta(_num_j * _num_c);
168 :
169 : // assign non-zero elements in k_eta
170 14726400 : for (const auto m : make_range(_num_c))
171 : {
172 : Real sum = 0.0;
173 :
174 29452800 : for (const auto n : make_range(_num_j))
175 22089600 : sum += (*_dhjdetai[n][i])[_qp] * (*_prop_ci[m * _num_j + n])[_qp];
176 :
177 7363200 : k_eta[m * _num_j + _num_j - 1] = -sum;
178 : }
179 :
180 : // compute x_eta
181 29452800 : for (const auto m : make_range(_num_j * _num_c))
182 : {
183 88358400 : for (const auto n : make_range(_num_j * _num_c))
184 66268800 : x_eta[m] += A(m, n) * k_eta[n];
185 : }
186 :
187 : // assign the values in x_eta to _dcidetaj
188 14726400 : for (const auto m : make_range(_num_c))
189 : {
190 29452800 : for (const auto n : make_range(_num_j))
191 22089600 : (*_dcidetaj[m][n][i])[_qp] = x_eta[m * _num_j + n];
192 : }
193 7363200 : }
194 2454400 : }
|
