Line data Source code
1 : /**********************************************************************/
2 : /* DO NOT MODIFY THIS HEADER */
3 : /* Swift, a Fourier spectral solver for MOOSE */
4 : /* */
5 : /* Copyright 2024 Battelle Energy Alliance, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /**********************************************************************/
8 :
9 : #include "FFTMechanics.h"
10 : #include "DomainAction.h"
11 : #include "MooseError.h"
12 : #include "SwiftUtils.h"
13 : #include <ATen/TensorIndexing.h>
14 : #include <ATen/core/TensorBody.h>
15 : #include <ATen/ops/unsqueeze_ops.h>
16 : #include <util/Optional.h>
17 :
18 : registerMooseObject("SwiftApp", FFTMechanics);
19 :
20 : InputParameters
21 0 : FFTMechanics::validParams()
22 : {
23 0 : InputParameters params = TensorOperator<>::validParams();
24 0 : params.addClassDescription("deGeus variational mechanics solve. Updates the coupled buffer "
25 : "holding the deformation gradient tensor.");
26 : // params.addParam<TensorInputBufferName>("C0", "Stiffness tensor estimate");
27 : // params.addParam<TensorInputBufferName>("C", "Stiffness tensor");
28 0 : params.addRequiredParam<TensorInputBufferName>("K", "Bulk modulus");
29 0 : params.addParam<TensorInputBufferName>("mu", "Shear modulus");
30 0 : params.addParam<Real>("l_tol", 1e-2, "Linear congugate gradient solve tolerance");
31 0 : params.addParam<unsigned int>("l_max_its",
32 : "Maximum number of congugate gradient solve iterations");
33 0 : params.addParam<Real>("nl_rel_tol", 1e-5, "Nonlinear solve absolute tolerance");
34 0 : params.addParam<Real>("nl_abs_tol", 1e-8, "Nonlinear solve relative tolerance");
35 0 : params.addParam<unsigned int>("nl_max_its", 100, "Maximum number of nonlinear solve iterations");
36 0 : params.addParam<TensorInputBufferName>("stress", "stress", "Computed stress");
37 0 : params.addParam<TensorInputBufferName>("tangent_operator", "dstressdstrain", "Tangent operator");
38 0 : params.addRequiredParam<TensorComputeName>("constitutive_model",
39 : "Tensor compute for the constitutive model (computes "
40 : "stress from displacement gradeint tensor)");
41 0 : params.addParam<TensorInputBufferName>("applied_macroscopic_strain",
42 : "Applied macroscopic strain");
43 0 : params.addParam<TensorInputBufferName>("F", "F", "Deformation gradient tensor.");
44 0 : params.addParam<bool>("verbose", false, "Print non-linear residuals.");
45 0 : return params;
46 0 : }
47 :
48 0 : FFTMechanics::FFTMechanics(const InputParameters & parameters)
49 : : TensorOperator<>(parameters),
50 0 : _ti(torch::eye(_dim, MooseTensor::floatTensorOptions())),
51 0 : _tI(MooseTensor::unsqueeze0(_ti, _dim)),
52 0 : _tI4(MooseTensor::unsqueeze0(torch::einsum("il,jk", {_ti, _ti}), _dim)),
53 0 : _tI4rt(MooseTensor::unsqueeze0(torch::einsum("ik,jl", {_ti, _ti}), _dim)),
54 0 : _tI4s((_tI4 + _tI4rt) / 2.0),
55 0 : _tII(MooseTensor::dyad22(_tI, _tI)),
56 0 : _tK(getInputBuffer("K")),
57 0 : _tmu(getInputBuffer("mu")),
58 0 : _r2_shape(_domain.getValueShape({_dim, _dim})),
59 0 : _tF(getInputBuffer<>("F")),
60 0 : _tP(getInputBuffer("stress")),
61 0 : _tK4(getInputBuffer("tangent_operator")),
62 0 : _l_tol(getParam<Real>("l_tol")),
63 0 : _l_max_its(isParamValid("l_max_its") ? getParam<unsigned int>("l_max_its")
64 0 : : _domain.getNumberOfCells()),
65 0 : _nl_rel_tol(getParam<Real>("nl_rel_tol")),
66 0 : _nl_abs_tol(getParam<Real>("nl_abs_tol")),
67 0 : _nl_max_its(getParam<unsigned int>("nl_max_its")),
68 0 : _constitutive_model(getCompute("constitutive_model")),
69 0 : _applied_macroscopic_strain(isParamValid("applied_macroscopic_strain")
70 0 : ? &getInputBuffer("applied_macroscopic_strain")
71 : : nullptr),
72 0 : _verbose(getParam<bool>("verbose"))
73 : {
74 : // Build projection tensor once
75 0 : const auto & q = _domain.getKGrid();
76 0 : const auto Q = _domain.getKSquare().unsqueeze(-1).unsqueeze(-1);
77 :
78 0 : auto M = torch::where(Q == 0, 0.0, q.unsqueeze(-2) * q.unsqueeze(-1) / Q);
79 :
80 0 : M = M.unsqueeze(-3).unsqueeze(-1);
81 :
82 0 : const auto delta_im = _ti.unsqueeze(1).unsqueeze(1).expand({_dim, _dim, _dim, _dim});
83 :
84 0 : _Ghat4 = (M * delta_im).to(MooseTensor::complexFloatTensorOptions());
85 0 : }
86 :
87 : void
88 0 : FFTMechanics::check()
89 : {
90 0 : const auto & stress_name = getParam<TensorOutputBufferName>("stress");
91 0 : if (!_constitutive_model.getSuppliedItems().count(stress_name))
92 0 : paramError("constitutive_model", "does not provide stress tensor '", stress_name, "'.");
93 0 : }
94 :
95 : void
96 0 : FFTMechanics::computeBuffer()
97 : {
98 : using namespace MooseTensor;
99 :
100 0 : const auto shape = _domain.getShape();
101 0 : std::vector<int64_t> _r2shape(shape.begin(), shape.end());
102 0 : _r2shape.push_back(_dim);
103 0 : _r2shape.push_back(_dim);
104 :
105 0 : const auto G = [&](const torch::Tensor & A2)
106 0 : { return _domain.ifft(ddot42(_Ghat4, _domain.fft(A2))).reshape(-1); };
107 0 : const auto K_dF = [&](const torch::Tensor & dFm)
108 0 : { return trans2(ddot42(_tK4, trans2(dFm.reshape(_r2shape)))); };
109 0 : const auto G_K_dF = [&](const torch::Tensor & dFm) { return G(K_dF(dFm)); };
110 :
111 : // initialize deformation gradient, and stress/stiffness [grid of tensors]
112 0 : _u = _tF;
113 0 : _constitutive_model.computeBuffer();
114 :
115 : // initial residual: distribute "barF" over grid using "K4"
116 0 : auto b = _applied_macroscopic_strain ? -G_K_dF(_applied_macroscopic_strain->expand(_r2_shape))
117 0 : : -G_K_dF(torch::zeros_like(_tF));
118 :
119 : // _u += DbarF;
120 0 : if (_applied_macroscopic_strain)
121 0 : _u = _u + _applied_macroscopic_strain->expand(_r2_shape);
122 :
123 : const auto Fn =
124 0 : at::linalg_norm(_u, c10::nullopt, c10::nullopt, false, c10::nullopt).cpu().item<double>();
125 :
126 : unsigned int iiter = 0;
127 0 : auto dFm = torch::zeros_like(b);
128 :
129 : // iterate as long as the iterative update does not vanish
130 : while (true)
131 : {
132 : const auto [dFm_new, iterations, lnorm] =
133 0 : conjugateGradientSolve(G_K_dF, b, dFm, _l_tol, _l_max_its);
134 : dFm = dFm_new;
135 :
136 : // update DOFs (array -> tens.grid)
137 0 : _u = _u + dFm.reshape(_r2_shape);
138 :
139 : // new residual stress and tangent
140 0 : _constitutive_model.computeBuffer();
141 :
142 : // convert res.stress to residual
143 0 : b = -G(_tP);
144 :
145 : const auto anorm =
146 0 : at::linalg_norm(dFm, c10::nullopt, c10::nullopt, false, c10::nullopt).cpu().item<double>();
147 0 : const auto rnorm = anorm / Fn;
148 :
149 : // print nonlinear residual to the screen
150 0 : if (_verbose)
151 0 : _console << "|R|=" << anorm << "\t|R/R0|=" << rnorm << '\n';
152 :
153 : // check convergence
154 0 : if ((rnorm < _nl_rel_tol || anorm < _nl_abs_tol) && iiter > 0)
155 : break;
156 :
157 0 : iiter++;
158 :
159 0 : if (iiter > _nl_max_its)
160 0 : paramError("nl_max_its",
161 : "Exceeded the maximum number of nonlinear iterations without converging.");
162 : }
163 0 : }
|
