Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/libmesh_config.h"
22 :
23 : // This class is only enabled in complex numbers builds
24 : #if defined(LIBMESH_USE_COMPLEX_NUMBERS)
25 :
26 : // Local includes
27 : #include "libmesh/frequency_system.h"
28 :
29 : #include "libmesh/equation_systems.h"
30 : #include "libmesh/int_range.h"
31 : #include "libmesh/libmesh_logging.h"
32 : #include "libmesh/linear_solver.h"
33 : #include "libmesh/numeric_vector.h"
34 :
35 : namespace libMesh
36 : {
37 :
38 1 : FrequencySystem::FrequencySystem (EquationSystems & es,
39 : const std::string & name_in,
40 1 : const unsigned int number_in) :
41 : LinearImplicitSystem (es, name_in, number_in),
42 1 : solve_system (nullptr),
43 1 : _finished_set_frequencies (false),
44 1 : _keep_solution_duplicates (true),
45 1 : _finished_init (false),
46 1 : _finished_assemble (false)
47 : {
48 : // default value for wave speed & fluid density
49 : //_equation_systems.parameters.set<Real>("wave speed") = 340.;
50 : //_equation_systems.parameters.set<Real>("rho") = 1.225;
51 1 : }
52 :
53 :
54 :
55 2 : FrequencySystem::~FrequencySystem () = default;
56 :
57 :
58 :
59 0 : void FrequencySystem::clear ()
60 : {
61 0 : LinearImplicitSystem::clear();
62 :
63 0 : _finished_set_frequencies = false;
64 0 : _keep_solution_duplicates = true;
65 0 : _finished_init = false;
66 0 : _finished_assemble = false;
67 :
68 : // We have to distinguish between the
69 : // simple straightforward "clear()"
70 : // and the clear that also touches the
71 : // EquationSystems parameters "current frequency" etc.
72 : // Namely, when reading from file (through equation_systems_io.C
73 : // methods), the param's are read in, then the systems.
74 : // Prior to reading a system, this system gets cleared...
75 : // And there, all the previously loaded frequency parameters
76 : // would get lost...
77 0 : }
78 :
79 :
80 :
81 0 : void FrequencySystem::clear_all ()
82 : {
83 0 : this->clear ();
84 :
85 : EquationSystems & es = this->get_equation_systems();
86 :
87 : // clear frequencies in the parameters section of the
88 : // EquationSystems object
89 0 : if (es.parameters.have_parameter<unsigned int> ("n_frequencies"))
90 : {
91 0 : unsigned int n_freq = es.parameters.get<unsigned int>("n_frequencies");
92 0 : for (unsigned int n=0; n < n_freq; n++)
93 0 : es.parameters.remove(this->form_freq_param_name(n));
94 0 : es.parameters.remove("current frequency");
95 : }
96 0 : }
97 :
98 :
99 :
100 :
101 1 : void FrequencySystem::init_data ()
102 : {
103 : // initialize parent data and additional solution vectors
104 1 : LinearImplicitSystem::init_data();
105 :
106 : // Log how long initializing the system takes
107 : LOG_SCOPE("init()", "FrequencySystem");
108 :
109 : EquationSystems & es =
110 : this->get_equation_systems();
111 :
112 : // make sure we have frequencies to solve for
113 1 : if (!_finished_set_frequencies)
114 : {
115 : // not supported for now
116 0 : if (parameters.have_parameter<unsigned int> ("n_frequencies"))
117 0 : libmesh_not_implemented_msg("ERROR: Setting the n_frequencies parameter on the system is not supported");
118 :
119 : // when this system was read from file, check
120 : // if this has a "n_frequencies" parameter,
121 : // and initialize us with these.
122 0 : if (es.parameters.have_parameter<unsigned int> ("n_frequencies"))
123 : {
124 : #ifndef NDEBUG
125 : const unsigned int n_freq =
126 : es.parameters.get<unsigned int>("n_frequencies");
127 :
128 : libmesh_assert_greater (n_freq, 0);
129 : #endif
130 0 : _finished_set_frequencies = true;
131 :
132 0 : this->set_current_frequency(0.);
133 : }
134 : else
135 0 : libmesh_error_msg("ERROR: Need to set frequencies before calling init().");
136 : }
137 :
138 1 : _finished_init = true;
139 1 : }
140 :
141 :
142 :
143 1 : void FrequencySystem::assemble ()
144 : {
145 : libmesh_assert (_finished_init);
146 :
147 1 : libmesh_error_msg_if(_finished_assemble, "ERROR: Matrices already assembled.");
148 :
149 : // Log how long assemble() takes
150 : LOG_SCOPE("assemble()", "FrequencySystem");
151 :
152 : // prepare matrix with the help of the _dof_map,
153 : // fill with sparsity pattern, initialize the
154 : // additional matrices
155 : LinearImplicitSystem::assemble();
156 :
157 : //matrix.print ();
158 : //rhs.print ();
159 :
160 1 : _finished_assemble = true;
161 1 : }
162 :
163 :
164 :
165 1 : void FrequencySystem::set_frequencies_by_steps (const Number base_freq,
166 : const Number freq_step,
167 : const unsigned int n_freq,
168 : const bool allocate_solution_duplicates)
169 : {
170 1 : this->_keep_solution_duplicates = allocate_solution_duplicates;
171 :
172 : // sanity check
173 1 : libmesh_error_msg_if(_finished_set_frequencies, "ERROR: frequencies already initialized.");
174 :
175 : EquationSystems & es =
176 : this->get_equation_systems();
177 :
178 : // store number of frequencies as parameter
179 1 : es.parameters.set<unsigned int>("n_frequencies") = n_freq;
180 :
181 4 : for (unsigned int n=0; n<n_freq; n++)
182 : {
183 : // remember frequencies as parameters
184 3 : es.parameters.set<Number>(this->form_freq_param_name(n)) =
185 3 : base_freq + Number(n) * freq_step;
186 :
187 : // build storage for solution vector, if wanted
188 3 : if (this->_keep_solution_duplicates)
189 6 : this->add_vector(this->form_solu_vec_name(n));
190 : }
191 :
192 1 : _finished_set_frequencies = true;
193 :
194 : // set the current frequency
195 1 : this->set_current_frequency(0);
196 1 : }
197 :
198 :
199 :
200 0 : void FrequencySystem::set_frequencies_by_range (const Number min_freq,
201 : const Number max_freq,
202 : const unsigned int n_freq,
203 : const bool allocate_solution_duplicates)
204 : {
205 0 : this->_keep_solution_duplicates = allocate_solution_duplicates;
206 :
207 : // sanity checks. Only look at the real part.
208 : libmesh_assert_greater_equal (std::real(max_freq), std::real(min_freq));
209 : libmesh_assert_greater (n_freq, 0);
210 0 : libmesh_error_msg_if(_finished_set_frequencies, "ERROR: frequencies already initialized.");
211 :
212 : EquationSystems & es =
213 : this->get_equation_systems();
214 :
215 : // store number of frequencies as parameter
216 0 : es.parameters.set<unsigned int>("n_frequencies") = n_freq;
217 :
218 : // set frequencies, build solution storage
219 0 : for (unsigned int n=0; n<n_freq; n++)
220 : {
221 : // remember frequencies as parameters
222 0 : es.parameters.set<Number>(this->form_freq_param_name(n)) =
223 0 : min_freq + static_cast<Number>(n)*(max_freq-min_freq)/static_cast<Number>(n_freq-1);
224 :
225 : // build storage for solution vector, if wanted
226 0 : if (this->_keep_solution_duplicates)
227 0 : System::add_vector(this->form_solu_vec_name(n));
228 : }
229 :
230 0 : _finished_set_frequencies = true;
231 :
232 : // set the current frequency
233 0 : this->set_current_frequency(0);
234 0 : }
235 :
236 :
237 :
238 0 : void FrequencySystem::set_frequencies (const std::vector<Real> & frequencies,
239 : const bool allocate_solution_duplicates)
240 : {
241 : libmesh_deprecated();
242 0 : this->_keep_solution_duplicates = allocate_solution_duplicates;
243 :
244 : // sanity checks
245 : libmesh_assert(!frequencies.empty());
246 0 : libmesh_error_msg_if(_finished_set_frequencies, "ERROR: frequencies already initialized.");
247 :
248 : EquationSystems & es =
249 : this->get_equation_systems();
250 :
251 : // store number of frequencies as parameter
252 0 : es.parameters.set<unsigned int>("n_frequencies") = frequencies.size();
253 :
254 : // set frequencies, build solution storage
255 0 : for (auto n : index_range(frequencies))
256 : {
257 : // remember frequencies as parameters
258 0 : es.parameters.set<Number>(this->form_freq_param_name(n)) = static_cast<Number>(frequencies[n]);
259 :
260 : // build storage for solution vector, if wanted
261 0 : if (this->_keep_solution_duplicates)
262 0 : System::add_vector(this->form_solu_vec_name(n));
263 : }
264 :
265 0 : _finished_set_frequencies = true;
266 :
267 : // set the current frequency
268 0 : this->set_current_frequency(0);
269 0 : }
270 :
271 :
272 0 : void FrequencySystem::set_frequencies (const std::vector<Number> & frequencies,
273 : const bool allocate_solution_duplicates)
274 : {
275 : libmesh_deprecated();
276 0 : this->_keep_solution_duplicates = allocate_solution_duplicates;
277 :
278 : // sanity checks
279 : libmesh_assert(!frequencies.empty());
280 0 : libmesh_error_msg_if(_finished_set_frequencies, "ERROR: frequencies already initialized.");
281 :
282 : EquationSystems & es =
283 : this->get_equation_systems();
284 :
285 : // store number of frequencies as parameter
286 0 : es.parameters.set<unsigned int>("n_frequencies") = frequencies.size();
287 :
288 : // set frequencies, build solution storage
289 0 : for (auto n : index_range(frequencies))
290 : {
291 : // remember frequencies as parameters
292 0 : es.parameters.set<Number>(this->form_freq_param_name(n)) = frequencies[n];
293 :
294 : // build storage for solution vector, if wanted
295 0 : if (this->_keep_solution_duplicates)
296 0 : System::add_vector(this->form_solu_vec_name(n));
297 : }
298 :
299 0 : _finished_set_frequencies = true;
300 :
301 : // set the current frequency
302 0 : this->set_current_frequency(0);
303 0 : }
304 :
305 :
306 :
307 :
308 0 : unsigned int FrequencySystem::n_frequencies () const
309 : {
310 : libmesh_assert(_finished_set_frequencies);
311 0 : return this->get_equation_systems().parameters.get<unsigned int>("n_frequencies");
312 : }
313 :
314 :
315 :
316 0 : void FrequencySystem::solve ()
317 : {
318 : libmesh_assert_greater (this->n_frequencies(), 0);
319 :
320 : // Solve for all the specified frequencies
321 0 : this->solve (0, this->n_frequencies()-1);
322 0 : }
323 :
324 :
325 :
326 3 : void FrequencySystem::solve (const unsigned int n_start,
327 : const unsigned int n_stop)
328 : {
329 : // Assemble the linear system, if not already done
330 3 : if (!_finished_assemble)
331 1 : this->assemble ();
332 :
333 : // the user-supplied solve method _has_ to be provided by the user
334 : libmesh_assert(solve_system);
335 :
336 : // existence & range checks
337 : libmesh_assert_greater (this->n_frequencies(), 0);
338 : libmesh_assert_less (n_stop, this->n_frequencies());
339 :
340 : EquationSystems & es =
341 : this->get_equation_systems();
342 :
343 : // Get the user-specified linear solver tolerance,
344 : // the user-specified maximum # of linear solver iterations,
345 : // the user-specified wave speed
346 3 : const auto [maxits, tol] = this->get_linear_solve_parameters();
347 :
348 : // start solver loop
349 6 : for (unsigned int n=n_start; n<= n_stop; n++)
350 : {
351 : // set the current frequency
352 3 : this->set_current_frequency(n);
353 :
354 : // Call the user-supplied pre-solve method
355 3 : LOG_CALL("user_pre_solve()", "FrequencySystem", this->solve_system(es, this->name()));
356 :
357 : // Solve the linear system for this specific frequency
358 : const std::pair<unsigned int, Real> rval =
359 3 : linear_solver->solve (*matrix, *solution, *rhs, tol, maxits);
360 :
361 3 : std::tie(_n_linear_iterations, _final_linear_residual) = rval;
362 3 : vec_rval.push_back(rval);
363 :
364 : /**
365 : * store the current solution in the additional vector
366 : */
367 3 : if (this->_keep_solution_duplicates)
368 6 : this->get_vector(this->form_solu_vec_name(n)) = *solution;
369 : }
370 :
371 : // sanity check
372 : //libmesh_assert_equal_to (vec_rval.size(), (n_stop-n_start+1));
373 3 : }
374 :
375 :
376 :
377 1 : void FrequencySystem::attach_solve_function(void fptr(EquationSystems & es,
378 : const std::string & name))
379 : {
380 : libmesh_assert(fptr);
381 :
382 1 : solve_system = fptr;
383 1 : }
384 :
385 :
386 :
387 4 : void FrequencySystem::set_current_frequency(unsigned int n)
388 : {
389 : libmesh_assert_less (n, n_frequencies());
390 :
391 : EquationSystems & es =
392 : this->get_equation_systems();
393 :
394 4 : es.parameters.set<Number>("current frequency") =
395 4 : es.parameters.get<Number>(this->form_freq_param_name(n));
396 4 : }
397 :
398 :
399 :
400 7 : std::string FrequencySystem::form_freq_param_name(const unsigned int n) const
401 : {
402 : libmesh_assert_less (n, 9999);
403 7 : std::string nstr = std::to_string(n);
404 7 : constexpr std::size_t digits = 4;
405 7 : const int zeros = digits - std::min(digits,nstr.size());
406 7 : nstr.insert(0, zeros, '0');
407 14 : return "frequency "+nstr;
408 : }
409 :
410 :
411 :
412 6 : std::string FrequencySystem::form_solu_vec_name(const unsigned int n) const
413 : {
414 : libmesh_assert_less (n, 9999);
415 6 : std::string nstr = std::to_string(n);
416 6 : constexpr std::size_t digits = 4;
417 6 : const int zeros = digits - std::min(digits,nstr.size());
418 6 : nstr.insert(0, zeros, '0');
419 12 : return "solution "+nstr;
420 : }
421 :
422 : } // namespace libMesh
423 :
424 :
425 : #endif // if defined(LIBMESH_USE_COMPLEX_NUMBERS)
|
