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 "BroydenSolver.h"
10 : #include "TensorProblem.h"
11 : #include "DomainAction.h"
12 :
13 : registerMooseObject("SwiftApp", BroydenSolver);
14 :
15 : InputParameters
16 0 : BroydenSolver::validParams()
17 : {
18 0 : InputParameters params = SplitOperatorBase::validParams();
19 0 : params.addClassDescription("Implicit secant solver time integration.");
20 0 : params.addParam<unsigned int>("substeps", 1, "secant solver substeps per time step.");
21 0 : params.addParam<unsigned int>("max_iterations", 5, "Maximum number of secant solver iteration.");
22 0 : params.addParam<Real>("relative_tolerance", 1e-9, "Convergence tolerance.");
23 0 : params.addParam<Real>("absolute_tolerance", 1e-9, "Convergence tolerance.");
24 0 : params.addParam<Real>("damping", 1.0, "Damping factor for the update step.");
25 0 : params.addParam<Real>(
26 0 : "initial_jacobian_guess", 1.0, "Factor for the initial inverse jacobian guess.");
27 0 : params.addParam<Real>(
28 0 : "dt_epsilon", 1e-4, "Semi-implicit stable timestep to bootstrap secant solve.");
29 0 : params.set<unsigned int>("substeps") = 0;
30 0 : params.addParam<bool>("verbose", false, "Show convergence history.");
31 0 : return params;
32 0 : }
33 :
34 0 : BroydenSolver::BroydenSolver(const InputParameters & parameters)
35 : : SplitOperatorBase(parameters),
36 : IterativeTensorSolverInterface(),
37 0 : _substeps(getParam<unsigned int>("substeps")),
38 0 : _max_iterations(getParam<unsigned int>("max_iterations")),
39 0 : _relative_tolerance(getParam<Real>("relative_tolerance")),
40 0 : _absolute_tolerance(getParam<Real>("absolute_tolerance")),
41 0 : _verbose(getParam<bool>("verbose")),
42 0 : _damping(getParam<Real>("damping")),
43 0 : _eye_factor(getParam<Real>("initial_jacobian_guess")),
44 0 : _dim(_domain.getDim()),
45 0 : _options(MooseTensor::complexFloatTensorOptions())
46 : {
47 : // no history required
48 0 : getVariables(0);
49 :
50 : // Jacobian dimensions
51 : const auto n = _variables.size();
52 0 : const auto & s = _domain.getReciprocalShape();
53 0 : std::vector<int64_t> v(s.begin(), s.end());
54 0 : v.push_back(n);
55 0 : v.push_back(n);
56 :
57 : // Create a 3x3 identity matrix and expand to all grid points
58 0 : _M = torch::eye(n, _options) * _eye_factor;
59 0 : for (const auto i : make_range(_dim))
60 : {
61 : libmesh_ignore(i);
62 : _M.unsqueeze_(0);
63 : }
64 0 : _M = _M.expand(v);
65 0 : }
66 :
67 : void
68 0 : BroydenSolver::computeBuffer()
69 : {
70 0 : for (_substep = 0; _substep < _substeps; ++_substep)
71 0 : broydenSolve();
72 0 : }
73 :
74 : void
75 0 : BroydenSolver::broydenSolve()
76 : {
77 0 : const auto dt = _dt / _substeps;
78 0 : const auto n = _variables.size();
79 :
80 : // stack u_old
81 0 : std::vector<torch::Tensor> u_old_v(n);
82 0 : for (const auto i : make_range(n))
83 0 : if (_variables[i]._reciprocal_buffer.defined())
84 : u_old_v[i] = _variables[i]._reciprocal_buffer;
85 : else
86 0 : u_old_v[i] = _domain.fft(_variables[i]._buffer);
87 0 : const auto u_old = torch::stack(u_old_v, -1);
88 :
89 0 : auto stackVariables = [&]()
90 : {
91 0 : std::vector<torch::Tensor> u(n);
92 0 : std::vector<torch::Tensor> N(n);
93 0 : std::vector<torch::Tensor> L(n);
94 0 : for (const auto i : make_range(n))
95 : {
96 0 : u[i] = _variables[i]._reciprocal_buffer;
97 0 : N[i] = _variables[i]._nonlinear_reciprocal;
98 : L[i] =
99 0 : _variables[i]._linear_reciprocal ? *(_variables[i]._linear_reciprocal) : torch::tensor(0);
100 : }
101 0 : return std::make_tuple(torch::stack(u, -1), torch::stack(N, -1), torch::stack(L, -1));
102 0 : };
103 :
104 : // initial residual
105 0 : _compute->computeBuffer();
106 0 : forwardBuffers();
107 :
108 0 : const auto [u0, N, L] = stackVariables();
109 : torch::Tensor u = u0;
110 0 : torch::Tensor R = (N + L * u) * dt;
111 :
112 : // initial residual norm
113 0 : const auto R0norm = torch::norm(R).item<double>();
114 :
115 : // secant iterations
116 0 : for (_iterations = 0; _iterations < _max_iterations; ++_iterations)
117 : {
118 : // check for convergence
119 0 : const auto Rnorm = torch::norm(R).item<double>();
120 :
121 : // NaN check
122 0 : if (std::isnan(Rnorm))
123 0 : mooseError("NAN!");
124 :
125 : // residual divergence check
126 0 : if (_iterations > 4 && Rnorm * 10.0 / _iterations > R0norm)
127 0 : mooseWarning("Diverging residual ", Rnorm, " ", Rnorm * 10.0 / _iterations, ' ', R0norm);
128 :
129 0 : if (Rnorm < _absolute_tolerance || Rnorm / R0norm < _relative_tolerance)
130 : {
131 0 : std::cout << "Broyden solve converged after " << _iterations << " iterations. |R|=" << Rnorm
132 0 : << " |R|/|R0|=" << Rnorm / R0norm << '\n';
133 0 : _is_converged = true;
134 0 : return;
135 : }
136 0 : else if (_verbose)
137 0 : std::cout << _iterations << " |R|=" << Rnorm << std::endl;
138 :
139 : // update step dx
140 0 : const auto sk = -torch::matmul(_M, R.unsqueeze(-1)); // column vector
141 0 : const auto skT = sk.squeeze(-1).unsqueeze(-2); // row vector
142 :
143 : // update u
144 0 : const auto u_out_v = torch::unbind(u + sk.squeeze(-1) * 0.5, -1);
145 0 : for (const auto i : make_range(n))
146 : {
147 : // look at min max here and maybe apply bounds?
148 0 : _variables[i]._buffer = _domain.ifft(u_out_v[i]);
149 : }
150 :
151 : // update residual
152 0 : _compute->computeBuffer();
153 0 : forwardBuffers();
154 :
155 0 : const auto [u0, N, L] = stackVariables();
156 : u = u0;
157 0 : const auto Rnew = (N + L * u) * dt + u_old - u;
158 :
159 : // residual change
160 0 : const auto yk = (Rnew - R).unsqueeze(-1);
161 0 : const auto ykT = yk.squeeze(-1).unsqueeze(-2);
162 :
163 : const auto denom = torch::matmul(skT, yk);
164 0 : if (_verbose)
165 : {
166 0 : const auto denom_zero = torch::sum(torch::where(torch::abs(denom) == 0, 1, 0)).item<double>();
167 0 : if (denom_zero)
168 0 : std::cout << "Matrix update denominator is zero in " << denom_zero << " entries.\n";
169 : }
170 :
171 0 : _M = _M + torch::where(torch::abs(denom) > 1e-12,
172 0 : torch::matmul((sk - torch::matmul(_M, yk)), skT) / denom,
173 : 0.0);
174 : // _M = _M + torch::matmul((sk - torch::matmul(_M, yk)), skT) /
175 : // torch::where(torch::abs(denom) > 1e-12, denom, 1.0);
176 :
177 : // _M = _M + torch::matmul((sk - torch::matmul(_M, yk)), skT) / torch::matmul(skT, yk);
178 :
179 : R = Rnew;
180 0 : }
181 :
182 0 : std::cerr << "Broyden solve did not converge within the maximum number of iterations.\n";
183 0 : _is_converged = false;
184 0 : }
|
