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 "AdamsBashforthMoulton.h"
10 : #include "TensorProblem.h"
11 : #include "Conversion.h"
12 : #include "DomainAction.h"
13 : #include <array>
14 :
15 : registerMooseObject("SwiftApp", AdamsBashforthMoulton);
16 : registerMooseObjectRenamed("SwiftApp",
17 : SemiImplicitSolver,
18 : "10/01/2025 00:01",
19 : AdamsBashforthMoulton);
20 :
21 : namespace
22 : {
23 : // Max order supported (up to ABM5)
24 : constexpr std::size_t max_order = 5;
25 : }
26 :
27 : InputParameters
28 158 : AdamsBashforthMoulton::validParams()
29 : {
30 158 : InputParameters params = SplitOperatorBase::validParams();
31 158 : params.addClassDescription("Adams-Bashforth-Moulton semi-implicit/explicit time integration "
32 : "solver with optional implicit corrector.");
33 316 : params.addParam<unsigned int>("substeps", 1, "semi-implicit substeps per time step.");
34 316 : params.addRangeCheckedParam<std::size_t>("predictor_order",
35 316 : 2,
36 158 : "predictor_order > 0 & predictor_order <= " +
37 158 : Moose::stringify(max_order),
38 : "Order of the Adams-Bashforth predictor.");
39 316 : params.addRangeCheckedParam<std::size_t>("corrector_order",
40 316 : 2,
41 158 : "corrector_order > 0 & corrector_order <= " +
42 158 : Moose::stringify(max_order),
43 : "Order of the Adams-Moulton corrector.");
44 316 : params.addParam<std::size_t>(
45 : "corrector_steps",
46 316 : 0,
47 : "Number the Adams-Moulton corrector steps to take (one is usually sufficient).");
48 158 : return params;
49 0 : }
50 :
51 78 : AdamsBashforthMoulton::AdamsBashforthMoulton(const InputParameters & parameters)
52 : : SplitOperatorBase(parameters),
53 78 : _substeps(getParam<unsigned int>("substeps")),
54 156 : _predictor_order(getParam<std::size_t>("predictor_order") - 1),
55 156 : _corrector_order(getParam<std::size_t>("corrector_order") - 1),
56 156 : _corrector_steps(getParam<std::size_t>("corrector_steps")),
57 78 : _sub_dt(_tensor_problem.subDt()),
58 156 : _sub_time(_tensor_problem.subTime())
59 : {
60 78 : getVariables(_predictor_order);
61 78 : }
62 :
63 : void
64 1004 : AdamsBashforthMoulton::computeBuffer()
65 : {
66 : // Adams–Bashforth coefficients (zero-padded)
67 1004 : constexpr std::array<std::array<double, max_order>, max_order> beta = {{
68 : {1.0, 0.0, 0.0, 0.0, 0.0}, // AB1
69 : {3.0 / 2.0, -1.0 / 2.0, 0.0, 0.0, 0.0}, // AB2
70 : {23.0 / 12.0, -16.0 / 12.0, 5.0 / 12.0, 0.0, 0.0}, // AB3
71 : {55.0 / 24.0, -59.0 / 24.0, 37.0 / 24.0, -9.0 / 24.0, 0.0}, // AB4
72 : {190.0 / 720.0, -2774.0 / 720.0, 2616.0 / 720.0, -1274.0 / 720.0, 251.0 / 720.0}, // AB5
73 : }};
74 :
75 : // Adams–Moulton coefficients (zero-padded)
76 1004 : constexpr std::array<std::array<double, max_order>, max_order> alpha = {{
77 : {1.0, 0.0, 0.0, 0.0, 0.0}, // AM1
78 : {0.5, 0.5, 0.0, 0.0, 0.0}, // AM2
79 : {5.0 / 12.0, 8.0 / 12.0, -1.0 / 12.0, 0.0, 0.0}, // AM3
80 : {9.0 / 24.0, 19.0 / 24.0, -5.0 / 24.0, 1.0 / 24.0, 0.0}, // AM4
81 : {251.0 / 720.0, 646.0 / 720.0, -264.0 / 720.0, 106.0 / 720.0, -19.0 / 720.0}, // AM5
82 : }};
83 :
84 1004 : const bool dt_changed = (_dt != _dt_old);
85 :
86 : torch::Tensor ubar;
87 1004 : _sub_dt = _dt / _substeps;
88 :
89 : // subcycles
90 50708 : for (const auto substep : make_range(_substeps))
91 : {
92 : // re-evaluate the solve compute
93 49704 : _compute->computeBuffer();
94 49704 : forwardBuffers();
95 :
96 : // Adams-Bashforth predictor on all variables
97 49704 : for (auto & [u,
98 : reciprocal_buffer,
99 : linear_reciprocal,
100 : nonlinear_reciprocal,
101 147408 : old_nonlinear_reciprocal] : _variables)
102 : {
103 97704 : const auto n_old = old_nonlinear_reciprocal.size();
104 :
105 : // Order is what the user requested, or what the available history allows for.
106 : // If dt changes between steps, we start at first order again
107 : const auto order =
108 195408 : std::min(substep < _predictor_order && dt_changed ? 0 : n_old, _predictor_order);
109 :
110 : // Adams-Bashforth
111 195408 : ubar = reciprocal_buffer + (_sub_dt * beta[order][0]) * nonlinear_reciprocal;
112 262382 : for (const auto i : make_range(order))
113 494034 : ubar += (_sub_dt * beta[order][i + 1]) * old_nonlinear_reciprocal[i];
114 :
115 97704 : if (linear_reciprocal)
116 390816 : ubar /= (1.0 - _sub_dt * *linear_reciprocal);
117 :
118 195408 : u = _domain.ifft(ubar);
119 : }
120 :
121 : // AB: y[n+1] = y[n] + dt * f(y[n])
122 : // AM: y[n+1] = y[n] + dt * f(y[n+1])
123 :
124 : // increment substep time
125 49704 : _sub_time += _sub_dt;
126 :
127 : // Adams-Moulton corrector
128 49704 : if (_corrector_steps)
129 : {
130 : // we need to keep the previous time step reciprocal_buffer if we run the AM corrector
131 3000 : std::vector<torch::Tensor> ubar_n(_variables.size());
132 9000 : for (const auto k : index_range(_variables))
133 6000 : ubar_n[k] = _variables[k]._reciprocal_buffer;
134 :
135 : // if the corrector order is AM2 or higher we also need the f calculated by AB prior to the
136 : // update to the current step values
137 : std::vector<torch::Tensor> N_n;
138 3000 : if (_corrector_order > 0)
139 : {
140 1000 : N_n.resize(_variables.size());
141 3000 : for (const auto k : index_range(_variables))
142 2000 : N_n[k] = _variables[k]._nonlinear_reciprocal;
143 : }
144 :
145 : // apply multiple corrector steps, going forward we probably want to allow users to fixpoint
146 : // iterate until a given convergence criterion is fulfilled.
147 8000 : for (std::size_t j = 0; j < _corrector_steps; ++j)
148 : {
149 : // re-evaluate the solve compute with the predicted variable values
150 5000 : _compute->computeBuffer();
151 5000 : forwardBuffers();
152 :
153 15000 : for (const auto k : index_range(_variables))
154 : {
155 10000 : auto & u = _variables[k]._buffer;
156 10000 : const auto * linear_reciprocal = _variables[k]._linear_reciprocal;
157 10000 : const auto & nonlinear_reciprocal_pred = _variables[k]._nonlinear_reciprocal;
158 10000 : const auto & old_nonlinear_reciprocal = _variables[k]._old_nonlinear_reciprocal;
159 :
160 : const auto n_old = old_nonlinear_reciprocal.size();
161 : const auto order =
162 10000 : std::min(substep < _corrector_order && dt_changed ? 1 : n_old + 1, _corrector_order);
163 10000 : if (order == 0)
164 6000 : continue;
165 :
166 12000 : ubar = ubar_n[k] + (_sub_dt * alpha[order][0]) * nonlinear_reciprocal_pred;
167 :
168 : if (order > 0)
169 : {
170 8000 : ubar += (_sub_dt * alpha[order][1]) * N_n[k];
171 4000 : for (const auto i : make_range(order - 1))
172 0 : ubar += (_sub_dt * alpha[order][i + 2]) * old_nonlinear_reciprocal[i];
173 : }
174 :
175 4000 : if (linear_reciprocal)
176 16000 : ubar /= (1.0 - _sub_dt * *linear_reciprocal);
177 :
178 8000 : u = _domain.ifft(ubar);
179 : }
180 : }
181 3000 : }
182 :
183 : // we skip the advanceState on the last substep because MOOSE will call that automatically
184 49704 : if (substep < _substeps - 1)
185 48700 : _tensor_problem.advanceState();
186 : }
187 1004 : }
|
