Line data Source code
1 : /**********************************************************************/
2 : /* DO NOT MODIFY THIS HEADER */
3 : /* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
4 : /* */
5 : /* Copyright 2017 Battelle Energy Alliance, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /**********************************************************************/
8 :
9 : #include "SpectralExecutionerLinearElastic.h"
10 : #include "InputParameters.h"
11 : #include "RankFourTensor.h"
12 :
13 : registerMooseObject("MagpieApp", SpectralExecutionerLinearElastic);
14 :
15 : InputParameters
16 6 : SpectralExecutionerLinearElastic::validParams()
17 : {
18 6 : InputParameters params = SpectralExecutionerBase::validParams();
19 6 : params.addClassDescription("Executioner for spectral solve of diffusion equation");
20 12 : params.addParam<Real>("time_step", 1.0, "Time step for ODE integration");
21 12 : params.addParam<unsigned int>("number_iterations", 0, "Maximum number of iterations for solver");
22 12 : params.addParam<Real>("young_modulus", 1.0, "First parameter for isotropic materials");
23 12 : params.addParam<Real>("poisson_ratio", 1.0, "Second parameter for isotropic materials");
24 12 : params.addParam<Real>("average_material_factor", 1.0, "Homogeneized factor for multiphase");
25 12 : params.addRequiredParam<std::vector<Real>>(
26 : "global_strain_tensor",
27 : "Vector of values defining the constant applied global strain "
28 : "to add. Components are XX, YY, ZZ, YZ, XZ, XY");
29 12 : params.addParam<Real>("solver_error", 1.0e-4, "Error for fixed iteration solver");
30 :
31 6 : return params;
32 0 : }
33 :
34 3 : SpectralExecutionerLinearElastic::SpectralExecutionerLinearElastic(
35 3 : const InputParameters & parameters)
36 : : SpectralExecutionerBase(parameters),
37 3 : _dt(getParam<Real>("time_step")),
38 6 : _nsteps(getParam<unsigned int>("number_iterations")),
39 6 : _young_modulus(getParam<Real>("young_modulus")),
40 6 : _poisson_ratio(getParam<Real>("poisson_ratio")),
41 6 : _average_factor(getParam<Real>("average_material_factor")),
42 6 : _solver_error(getParam<Real>("solver_error")),
43 6 : _t_current(0.0)
44 : {
45 6 : _initial_strain_tensor.fillFromInputVector(getParam<std::vector<Real>>("global_strain_tensor"));
46 3 : }
47 :
48 : void
49 3 : SpectralExecutionerLinearElastic::fillOutEpsilonBuffer(
50 : FFTBufferBase<RankTwoTensor> & epsilon_buffer)
51 : {
52 3 : epsilon_buffer.realSpace() = _initial_strain_tensor;
53 3 : }
54 :
55 : void
56 3 : SpectralExecutionerLinearElastic::getGreensFunction(FFTBufferBase<RankFourTensor> & gamma_hat,
57 : const RankFourTensor & elasticity_tensor)
58 : {
59 : const int ndim = 3;
60 : std::size_t index = 0;
61 :
62 : const Complex I(0.0, 1.0);
63 : auto & gamma_hat_F = gamma_hat.reciprocalSpace();
64 :
65 : // Fill fourth-order Green operator based on homogeneous material properties
66 99 : for (auto ivec : gamma_hat.kTable(0))
67 3168 : for (auto jvec : gamma_hat.kTable(1))
68 55296 : for (auto kvec : gamma_hat.kTable(2))
69 : {
70 52224 : const std::array<Complex, 3> freq{ivec * I, jvec * I, kvec * I};
71 :
72 52224 : Real lambda0 = _young_modulus * _average_factor * _poisson_ratio /
73 52224 : ((1 + _poisson_ratio) * (1 - 2 * _poisson_ratio));
74 52224 : Real nu0 = _young_modulus * _average_factor / (2 * (1 + _poisson_ratio));
75 52224 : Real constant = (lambda0 + nu0) / (nu0 * (lambda0 + 2.0 * nu0));
76 :
77 208896 : for (int i = 0; i < ndim; i++)
78 626688 : for (int j = 0; j < ndim; j++)
79 1880064 : for (int k = 0; k < ndim; k++)
80 5640192 : for (int l = 0; l < ndim; l++)
81 : {
82 : Complex q_square = freq[0] * freq[0] + freq[1] * freq[1] + freq[2] * freq[2];
83 4230144 : if (std::abs(q_square) > 1.0e-12)
84 : {
85 4229901 : gamma_hat_F[index](i, j, k, l) =
86 4229901 : -1.0 * constant * (freq[i] * freq[j] * freq[k] * freq[l]) /
87 : (q_square * q_square) +
88 7049835 : (Real(i == k) * freq[j] * freq[l] + Real(j == k) * freq[i] * freq[l] +
89 7049835 : Real(i == l) * freq[j] * freq[k] + Real(j == l) * freq[i] * freq[k]) /
90 4229901 : (4.0 * nu0 * q_square);
91 : }
92 : else
93 243 : gamma_hat_F[index] = elasticity_tensor.invSymm();
94 : }
95 :
96 52224 : index++;
97 : }
98 3 : }
99 :
100 : FFTBufferBase<RankTwoTensor> &
101 3 : SpectralExecutionerLinearElastic::getInitialStress(
102 : FFTBufferBase<RankTwoTensor> & epsilon_buffer,
103 : const FFTBufferBase<RankFourTensor> & elastic_tensor)
104 : {
105 3 : auto & stress_buffer = getFFTBuffer<RankTwoTensor>("stress");
106 :
107 : // Homogeneous initial state of strain
108 3 : fillOutEpsilonBuffer(epsilon_buffer);
109 :
110 : // Set real stress buffer to product E * epsilon
111 3 : stress_buffer.realSpace().setToProductRealSpace(elastic_tensor.realSpace(),
112 : epsilon_buffer.realSpace());
113 3 : return stress_buffer;
114 : }
115 :
116 : void
117 15 : SpectralExecutionerLinearElastic::advanceReciprocalEpsilon(
118 : FFTBufferBase<RankTwoTensor> & epsilon_buffer,
119 : const FFTBufferBase<RankTwoTensor> & stress_buffer,
120 : const FFTBufferBase<RankFourTensor> & gamma_hat)
121 : {
122 15 : Complex I(1.0, 0.0);
123 :
124 15 : epsilon_buffer.reciprocalSpace().applyLambdaReciprocalSpace(
125 261120 : [&gamma_hat, &stress_buffer, &epsilon_buffer](std::size_t index)
126 : {
127 : return epsilon_buffer.reciprocalSpace()[index] -
128 261120 : gamma_hat.reciprocalSpace()[index] * stress_buffer.reciprocalSpace()[index];
129 : });
130 :
131 : // Avoid divide by zero
132 15 : epsilon_buffer.reciprocalSpace()[0] = _initial_strain_tensor * I;
133 15 : }
134 :
135 : void
136 15 : SpectralExecutionerLinearElastic::updateRealSigma(
137 : const FFTBufferBase<RankTwoTensor> & epsilon_buffer,
138 : FFTBufferBase<RankTwoTensor> & stress_buffer,
139 : const FFTBufferBase<RankFourTensor> & elastic_tensor)
140 : {
141 : // Set real stress buffer to product E * epsilon
142 15 : stress_buffer.realSpace().setToProductRealSpace(elastic_tensor.realSpace(),
143 : epsilon_buffer.realSpace());
144 15 : }
145 :
146 : void
147 3 : SpectralExecutionerLinearElastic::filloutElasticTensor(
148 : const FFTBufferBase<Real> & ratio_buffer,
149 : FFTBufferBase<Real> & index_buffer,
150 : FFTBufferBase<RankFourTensor> & elastic_tensor_buffer)
151 : {
152 : const auto & grid = elastic_tensor_buffer.grid();
153 :
154 3 : int ni = grid[0];
155 3 : int nj = grid[1];
156 3 : int nk = grid[2];
157 :
158 : size_t index = 0;
159 3 : std::vector<Real> iso_const(2);
160 :
161 99 : for (int freq_x = 0; freq_x < ni; ++freq_x)
162 3168 : for (int freq_y = 0; freq_y < nj; ++freq_y)
163 101376 : for (int freq_z = 0; freq_z < nk; ++freq_z)
164 : {
165 : // Define elastic tensor
166 98304 : iso_const[0] = _young_modulus * ratio_buffer.realSpace()[index];
167 98304 : iso_const[1] = _poisson_ratio;
168 :
169 98304 : elastic_tensor_buffer.realSpace()[index].fillFromInputVector(
170 : iso_const, RankFourTensor::symmetric_isotropic_E_nu);
171 98304 : index_buffer.realSpace()[index] = index;
172 98304 : index++;
173 : }
174 3 : }
175 :
176 : bool
177 15 : SpectralExecutionerLinearElastic::hasStressConverged(const FFTBufferBase<RankTwoTensor> & stress)
178 : {
179 :
180 : const auto & grid = stress.grid();
181 :
182 15 : const int ni = grid[0];
183 15 : const int nj = grid[1];
184 15 : const int nk = (grid[2] >> 1) + 1;
185 :
186 : std::size_t index = 0;
187 :
188 : const auto & ivec = stress.kTable(0);
189 : const auto & jvec = stress.kTable(1);
190 : const auto & kvec = stress.kTable(2);
191 :
192 15 : const std::vector<int> grid_vector = stress.grid();
193 :
194 : const Complex I(0.0, 1.0);
195 :
196 : Complex error_n(0.0, 0.0);
197 : Complex error_0(0.0, 0.0);
198 :
199 : // Error: Ensure overall equilibrium
200 495 : for (int freq_x = 0; freq_x < ni; ++freq_x)
201 15840 : for (int freq_y = 0; freq_y < nj; ++freq_y)
202 276480 : for (int freq_z = 0; freq_z < nk; ++freq_z)
203 : {
204 261120 : const ComplexVectorValue freq{ivec[freq_x] * I, jvec[freq_y] * I, kvec[freq_z] * I};
205 :
206 : ComplexVectorValue kvector_stress;
207 261120 : kvector_stress = stress.reciprocalSpace()[index] * freq;
208 :
209 : error_n += kvector_stress(0) * kvector_stress(0) + kvector_stress(1) * kvector_stress(1) +
210 : kvector_stress(2) * kvector_stress(2);
211 :
212 261120 : if (freq_x == 0 && freq_y == 0 && freq_z == 0)
213 15 : error_0 = stress.reciprocalSpace()[0].L2norm() * stress.reciprocalSpace()[0].L2norm();
214 :
215 261120 : index++;
216 : }
217 :
218 15 : Real iteration_error = std::sqrt(std::norm(error_n)) / std::sqrt(std::norm(error_0));
219 : Moose::out << "Iteration error: " << iteration_error << "\n";
220 :
221 15 : if (iteration_error > _solver_error)
222 : return false;
223 : else
224 3 : return true;
225 : }
226 :
227 : void
228 3 : SpectralExecutionerLinearElastic::execute()
229 : {
230 : unsigned int thisStep = 0;
231 :
232 3 : _time_step = thisStep;
233 3 : _time = _time_step;
234 3 : _fe_problem.outputStep(EXEC_INITIAL);
235 3 : _fe_problem.advanceState();
236 :
237 3 : auto & epsilon_buffer = getFFTBuffer<RankTwoTensor>("epsilon");
238 3 : if (epsilon_buffer.dim() != 3)
239 0 : mooseError("Error: Problem dimension not implemented in SpectralExecutionerLinearElastic.");
240 :
241 3 : auto & ratio_buffer = getFFTBuffer<Real>("stiffness_ratio");
242 3 : auto & elastic_tensor_buffer = getFFTBuffer<RankFourTensor>("elastic");
243 3 : auto & index_buffer = getFFTBuffer<Real>("index_buffer");
244 :
245 : // Fill out space-varying elasticity tensor
246 3 : filloutElasticTensor(ratio_buffer, index_buffer, elastic_tensor_buffer);
247 :
248 : // Get corresponding initial stress (also fill out epsilon_buffer with initial strain)
249 3 : auto & stress_buffer = getInitialStress(epsilon_buffer, elastic_tensor_buffer);
250 :
251 : // Get specific Green's function
252 3 : auto & gamma_hat = getFFTBuffer<RankFourTensor>("gamma");
253 :
254 : thisStep++;
255 3 : _t_current += _dt;
256 3 : _time_step = thisStep;
257 :
258 3 : _fe_problem.execute(EXEC_FINAL);
259 3 : _time = _t_current;
260 :
261 : Moose::out << "_t_current: " << _t_current << ". \n";
262 :
263 3 : _fe_problem.outputStep(EXEC_FINAL);
264 :
265 : // Fill out a homogeneous elasticity tensor with some average properties
266 3 : std::vector<Real> iso_const(2);
267 3 : iso_const[0] = _young_modulus * _average_factor;
268 3 : iso_const[1] = _poisson_ratio;
269 3 : RankFourTensor elasticity_homo;
270 3 : elasticity_homo.fillFromInputVector(iso_const, RankFourTensor::symmetric_isotropic_E_nu);
271 3 : getGreensFunction(gamma_hat, elasticity_homo);
272 :
273 : // Our FFTW "many" plans do not preserve the input, so explicit copies are made
274 : FFTData<RankTwoTensor> stress_buffer_backup_real = stress_buffer.realSpace();
275 3 : epsilon_buffer.forward();
276 :
277 : FFTData<ComplexRankTwoTensor> epsilon_buffer_backup_reciprocal = epsilon_buffer.reciprocalSpace();
278 :
279 : bool is_converged = false;
280 :
281 15 : for (unsigned int step_no = 0; step_no < _nsteps; step_no++)
282 : {
283 : // Update sigma in the real space
284 15 : updateRealSigma(epsilon_buffer, stress_buffer, elastic_tensor_buffer);
285 :
286 15 : stress_buffer_backup_real = stress_buffer.realSpace();
287 15 : stress_buffer.forwardRaw();
288 15 : stress_buffer.reciprocalSpace() *= stress_buffer.backwardScale();
289 15 : stress_buffer.realSpace() = stress_buffer_backup_real;
290 :
291 : // Convergence check: Ensure global equilibrium
292 15 : is_converged = hasStressConverged(stress_buffer);
293 :
294 : // Compute new strain tensor in Fourier space
295 15 : epsilon_buffer.reciprocalSpace() = epsilon_buffer_backup_reciprocal;
296 15 : advanceReciprocalEpsilon(epsilon_buffer, stress_buffer, gamma_hat);
297 :
298 : // Cache reciprocal epsilon to avoid being overwritten upon backward (inverse) operation
299 15 : epsilon_buffer_backup_reciprocal = epsilon_buffer.reciprocalSpace();
300 15 : epsilon_buffer.backwardRaw();
301 :
302 15 : thisStep++;
303 15 : _t_current += _dt;
304 15 : _time_step = thisStep;
305 :
306 15 : _fe_problem.execute(EXEC_FINAL);
307 15 : _time = _t_current;
308 :
309 : Moose::out << "_t_current: " << _t_current << ". \n";
310 :
311 15 : _fe_problem.outputStep(EXEC_FINAL);
312 :
313 15 : if (is_converged)
314 : break;
315 :
316 12 : if (step_no != _nsteps - 1)
317 12 : _fe_problem.advanceState();
318 : }
319 3 : }
|
