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 62 : KKSPhaseConcentrationMultiPhaseDerivatives::validParams()
17 : {
18 62 : InputParameters params = DerivativeMaterialInterface<Material>::validParams();
19 62 : 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 124 : params.addRequiredCoupledVar("global_cs", "The interpolated concentrations c, b, etc");
24 124 : params.addRequiredCoupledVar("all_etas", "Order parameters.");
25 124 : 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 124 : params.addRequiredParam<std::vector<MaterialName>>(
30 : "Fj_names", "Free energy material objects in the same order as all_etas.");
31 124 : params.addRequiredParam<std::vector<MaterialPropertyName>>(
32 : "hj_names", "witching functions in the same order as all_etas.");
33 62 : return params;
34 0 : }
35 :
36 48 : KKSPhaseConcentrationMultiPhaseDerivatives::KKSPhaseConcentrationMultiPhaseDerivatives(
37 48 : const InputParameters & parameters)
38 : : DerivativeMaterialInterface<Material>(parameters),
39 48 : _num_c(coupledComponents("global_cs")),
40 48 : _c_names(coupledNames("global_cs")),
41 48 : _eta_names(coupledNames("all_etas")),
42 48 : _num_j(coupledComponents("all_etas")),
43 48 : _prop_ci(_num_c * _num_j),
44 96 : _ci_names(getParam<std::vector<MaterialPropertyName>>("ci_names")),
45 48 : _dcidetaj(_num_c),
46 48 : _dcidb(_num_c),
47 96 : _Fj_names(getParam<std::vector<MaterialName>>("Fj_names")),
48 48 : _d2Fidcidbi(_num_j),
49 96 : _hj_names(getParam<std::vector<MaterialPropertyName>>("hj_names")),
50 48 : _prop_hj(_num_j),
51 96 : _dhjdetai(_num_j)
52 :
53 : {
54 192 : for (const auto m : make_range(_num_c * _num_j))
55 144 : _prop_ci[m] = &getMaterialPropertyByName<Real>(_ci_names[m]);
56 :
57 96 : for (const auto m : make_range(_num_c))
58 : {
59 48 : _dcidb[m].resize(_num_j);
60 48 : _dcidetaj[m].resize(_num_j);
61 :
62 192 : for (const auto n : make_range(_num_j))
63 : {
64 144 : _dcidb[m][n].resize(_num_c);
65 144 : _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 288 : for (const auto l : make_range(_num_c))
70 144 : _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 576 : for (const auto l : make_range(_num_j))
75 432 : _dcidetaj[m][n][l] =
76 864 : &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 192 : for (const auto m : make_range(_num_j))
84 : {
85 144 : _d2Fidcidbi[m].resize(_num_c);
86 :
87 288 : for (const auto n : make_range(_num_c))
88 : {
89 144 : _d2Fidcidbi[m][n].resize(_num_c);
90 :
91 288 : for (const auto l : make_range(_num_c))
92 144 : _d2Fidcidbi[m][n][l] = &getMaterialPropertyDerivative<Real>(
93 144 : _Fj_names[m], _ci_names[m + n * _num_j], _ci_names[m + l * _num_j]);
94 : }
95 : }
96 :
97 192 : for (const auto m : make_range(_num_j))
98 : {
99 144 : _prop_hj[m] = &getMaterialPropertyByName<Real>(_hj_names[m]);
100 :
101 144 : _dhjdetai[m].resize(_num_j);
102 :
103 576 : for (const auto n : make_range(_num_j))
104 432 : _dhjdetai[m][n] = &getMaterialPropertyDerivative<Real>(_hj_names[m], _eta_names[n]);
105 : }
106 48 : }
107 :
108 : void
109 2073600 : KKSPhaseConcentrationMultiPhaseDerivatives::computeQpProperties()
110 : {
111 : // declare Jacobian matrix A
112 2073600 : 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 4147200 : for (const auto m : make_range(_num_c))
119 : {
120 : // equal chemical potential derivative equations
121 6220800 : for (const auto n : make_range(_num_j - 1))
122 : {
123 8294400 : for (const auto l : make_range(_num_c))
124 : {
125 4147200 : A(m * _num_j + n, n + l * _num_j) = (*_d2Fidcidbi[n][m][l])[_qp];
126 4147200 : 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 8294400 : for (const auto n : make_range(_num_j))
132 6220800 : A((m + 1) * _num_j - 1, m * _num_j + n) = (*_prop_hj[n])[_qp];
133 : }
134 :
135 2073600 : 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 4147200 : for (const auto i : make_range(_num_c))
140 : {
141 2073600 : std::vector<Real> k_c(_num_j * _num_c);
142 2073600 : std::vector<Real> x_c(_num_j * _num_c);
143 :
144 : // assign non-zero elements in k_c
145 2073600 : k_c[i * _num_j + _num_j - 1] = 1;
146 :
147 : // compute x_c
148 8294400 : for (const auto m : make_range(_num_j * _num_c))
149 : {
150 24883200 : for (const auto n : make_range(_num_j * _num_c))
151 18662400 : x_c[m] += A(m, n) * k_c[n];
152 : }
153 :
154 : // assign the values in x_c to _dcidb
155 4147200 : for (const auto m : make_range(_num_c))
156 : {
157 8294400 : for (const auto n : make_range(_num_j))
158 6220800 : (*_dcidb[m][n][i])[_qp] = x_c[m * _num_j + n];
159 : }
160 2073600 : }
161 :
162 : // solve linear system of constraint derivatives wrt eta for computing dcidetaj use the same
163 : // linear matrix as computing dcidb
164 8294400 : for (const auto i : make_range(_num_j))
165 : {
166 6220800 : std::vector<Real> k_eta(_num_j * _num_c);
167 6220800 : std::vector<Real> x_eta(_num_j * _num_c);
168 :
169 : // assign non-zero elements in k_eta
170 12441600 : for (const auto m : make_range(_num_c))
171 : {
172 : Real sum = 0.0;
173 :
174 24883200 : for (const auto n : make_range(_num_j))
175 18662400 : sum += (*_dhjdetai[n][i])[_qp] * (*_prop_ci[m * _num_j + n])[_qp];
176 :
177 6220800 : k_eta[m * _num_j + _num_j - 1] = -sum;
178 : }
179 :
180 : // compute x_eta
181 24883200 : for (const auto m : make_range(_num_j * _num_c))
182 : {
183 74649600 : for (const auto n : make_range(_num_j * _num_c))
184 55987200 : x_eta[m] += A(m, n) * k_eta[n];
185 : }
186 :
187 : // assign the values in x_eta to _dcidetaj
188 12441600 : for (const auto m : make_range(_num_c))
189 : {
190 24883200 : for (const auto n : make_range(_num_j))
191 18662400 : (*_dcidetaj[m][n][i])[_qp] = x_eta[m * _num_j + n];
192 : }
193 6220800 : }
194 2073600 : }
|
