libMesh
frequency_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 /*
24  * Require complex arithmetic
25  */
26 #if defined(LIBMESH_USE_COMPLEX_NUMBERS)
27 
28 
29 // C++ includes
30 #include <cstdio> // for sprintf
31 
32 // Local includes
33 #include "libmesh/frequency_system.h"
34 #include "libmesh/equation_systems.h"
35 #include "libmesh/libmesh_logging.h"
36 #include "libmesh/linear_solver.h"
37 #include "libmesh/numeric_vector.h"
38 
39 namespace libMesh
40 {
41 
42 
43 
44 // ------------------------------------------------------------
45 // FrequencySystem implementation
47  const std::string & name_in,
48  const unsigned int number_in) :
49  LinearImplicitSystem (es, name_in, number_in),
50  solve_system (libmesh_nullptr),
51  _finished_set_frequencies (false),
52  _keep_solution_duplicates (true),
53  _finished_init (false),
54  _finished_assemble (false)
55 {
56  // default value for wave speed & fluid density
57  //_equation_systems.parameters.set<Real>("wave speed") = 340.;
58  //_equation_systems.parameters.set<Real>("rho") = 1.225;
59 }
60 
61 
62 
64 {
65  this->clear ();
66 
67  // the additional matrices and vectors are cleared and zero'ed in System
68 }
69 
70 
71 
72 
74 {
76 
79  _finished_init = false;
80  _finished_assemble = false;
81 
82  /*
83  * We have to distinguish between the
84  * simple straightforward "clear()"
85  * and the clear that also touches the
86  * EquationSystems parameters "current frequency" etc.
87  * Namely, when reading from file (through equation_systems_io.C
88  * methods), the param's are read in, then the systems.
89  * Prior to reading a system, this system gets cleared...
90  * And there, all the previously loaded frequency parameters
91  * would get lost...
92  */
93 }
94 
95 
96 
98 {
99  this->clear ();
100 
101  EquationSystems & es = this->get_equation_systems();
102 
103  // clear frequencies in the parameters section of the
104  // EquationSystems object
105  if (es.parameters.have_parameter<unsigned int> ("n_frequencies"))
106  {
107  unsigned int n_freq = es.parameters.get<unsigned int>("n_frequencies");
108  for (unsigned int n=0; n < n_freq; n++)
109  es.parameters.remove(this->form_freq_param_name(n));
110  es.parameters.remove("current frequency");
111  }
112 }
113 
114 
115 
116 
118 {
119  // initialize parent data and additional solution vectors
121 
122  // Log how long initializing the system takes
123  LOG_SCOPE("init()", "FrequencySystem");
124 
125  EquationSystems & es =
126  this->get_equation_systems();
127 
128  // make sure we have frequencies to solve for
130  {
131  /*
132  * when this system was read from file, check
133  * if this has a "n_frequencies" parameter,
134  * and initialize us with these.
135  */
136  if (es.parameters.have_parameter<unsigned int> ("n_frequencies"))
137  {
138 #ifndef NDEBUG
139  const unsigned int n_freq =
140  es.parameters.get<unsigned int>("n_frequencies");
141 
142  libmesh_assert_greater (n_freq, 0);
143 #endif
145 
146  this->set_current_frequency(0.);
147  }
148  else
149  libmesh_error_msg("ERROR: Need to set frequencies before calling init().");
150  }
151 
152  _finished_init = true;
153 }
154 
155 
156 
158 {
160 
161  if (_finished_assemble)
162  libmesh_error_msg("ERROR: Matrices already assembled.");
163 
164  // Log how long assemble() takes
165  LOG_SCOPE("assemble()", "FrequencySystem");
166 
167  // prepare matrix with the help of the _dof_map,
168  // fill with sparsity pattern, initialize the
169  // additional matrices
171 
172  //matrix.print ();
173  //rhs.print ();
174 
175  _finished_assemble = true;
176 }
177 
178 
179 
181  const Number freq_step,
182  const unsigned int n_freq,
183  const bool allocate_solution_duplicates)
184 {
185  this->_keep_solution_duplicates = allocate_solution_duplicates;
186 
187  // sanity check
189  libmesh_error_msg("ERROR: frequencies already initialized.");
190 
191  EquationSystems & es =
192  this->get_equation_systems();
193 
194  // store number of frequencies as parameter
195  es.parameters.set<unsigned int>("n_frequencies") = n_freq;
196 
197  for (unsigned int n=0; n<n_freq; n++)
198  {
199  // remember frequencies as parameters
200  es.parameters.set<Number>(this->form_freq_param_name(n)) =
201  base_freq + Number(n) * freq_step;
202 
203  // build storage for solution vector, if wanted
204  if (this->_keep_solution_duplicates)
205  this->add_vector(this->form_solu_vec_name(n));
206  }
207 
209 
210  // set the current frequency
211  this->set_current_frequency(0);
212 }
213 
214 
215 
217  const Number max_freq,
218  const unsigned int n_freq,
219  const bool allocate_solution_duplicates)
220 {
221  this->_keep_solution_duplicates = allocate_solution_duplicates;
222 
223  // sanity checks. Only look at the real part.
224  libmesh_assert_greater_equal (std::real(max_freq), std::real(min_freq));
225  libmesh_assert_greater (n_freq, 0);
226 
228  libmesh_error_msg("ERROR: frequencies already initialized.");
229 
230  EquationSystems & es =
231  this->get_equation_systems();
232 
233  // store number of frequencies as parameter
234  es.parameters.set<unsigned int>("n_frequencies") = n_freq;
235 
236  // set frequencies, build solution storage
237  for (unsigned int n=0; n<n_freq; n++)
238  {
239  // remember frequencies as parameters
240  es.parameters.set<Number>(this->form_freq_param_name(n)) =
241  min_freq + static_cast<Number>(n)*(max_freq-min_freq)/static_cast<Number>(n_freq-1);
242 
243  // build storage for solution vector, if wanted
244  if (this->_keep_solution_duplicates)
246  }
247 
249 
250  // set the current frequency
251  this->set_current_frequency(0);
252 }
253 
254 
255 
256 void FrequencySystem::set_frequencies (const std::vector<Real> & frequencies,
257  const bool allocate_solution_duplicates)
258 {
259  libmesh_deprecated();
260  this->_keep_solution_duplicates = allocate_solution_duplicates;
261 
262  // sanity checks
263  libmesh_assert(!frequencies.empty());
264 
266  libmesh_error_msg("ERROR: frequencies already initialized.");
267 
268  EquationSystems & es =
269  this->get_equation_systems();
270 
271  // store number of frequencies as parameter
272  es.parameters.set<unsigned int>("n_frequencies") = frequencies.size();
273 
274  // set frequencies, build solution storage
275  for (std::size_t n=0; n<frequencies.size(); n++)
276  {
277  // remember frequencies as parameters
278  es.parameters.set<Number>(this->form_freq_param_name(n)) = static_cast<Number>(frequencies[n]);
279 
280  // build storage for solution vector, if wanted
281  if (this->_keep_solution_duplicates)
283  }
284 
286 
287  // set the current frequency
288  this->set_current_frequency(0);
289 }
290 
291 
292 void FrequencySystem::set_frequencies (const std::vector<Number> & frequencies,
293  const bool allocate_solution_duplicates)
294 {
295  libmesh_deprecated();
296  this->_keep_solution_duplicates = allocate_solution_duplicates;
297 
298  // sanity checks
299  libmesh_assert(!frequencies.empty());
300 
302  libmesh_error_msg("ERROR: frequencies already initialized.");
303 
304  EquationSystems & es =
305  this->get_equation_systems();
306 
307  // store number of frequencies as parameter
308  es.parameters.set<unsigned int>("n_frequencies") = frequencies.size();
309 
310  // set frequencies, build solution storage
311  for (std::size_t n=0; n<frequencies.size(); n++)
312  {
313  // remember frequencies as parameters
314  es.parameters.set<Number>(this->form_freq_param_name(n)) = frequencies[n];
315 
316  // build storage for solution vector, if wanted
317  if (this->_keep_solution_duplicates)
319  }
320 
322 
323  // set the current frequency
324  this->set_current_frequency(0);
325 }
326 
327 
328 
329 
330 unsigned int FrequencySystem::n_frequencies () const
331 {
333  return this->get_equation_systems().parameters.get<unsigned int>("n_frequencies");
334 }
335 
336 
337 
339 {
340  libmesh_assert_greater (this->n_frequencies(), 0);
341 
342  // Solve for all the specified frequencies
343  this->solve (0, this->n_frequencies()-1);
344 }
345 
346 
347 
348 void FrequencySystem::solve (const unsigned int n_start,
349  const unsigned int n_stop)
350 {
351  // Assemble the linear system, if not already done
352  if (!_finished_assemble)
353  this->assemble ();
354 
355  // the user-supplied solve method _has_ to be provided by the user
357 
358  // existence & range checks
359  libmesh_assert_greater (this->n_frequencies(), 0);
360  libmesh_assert_less (n_stop, this->n_frequencies());
361 
362  EquationSystems & es =
363  this->get_equation_systems();
364 
365  // Get the user-specified linear solver tolerance,
366  // the user-specified maximum # of linear solver iterations,
367  // the user-specified wave speed
368  const Real tol =
369  es.parameters.get<Real>("linear solver tolerance");
370  const unsigned int maxits =
371  es.parameters.get<unsigned int>("linear solver maximum iterations");
372 
373 
374  // // return values
375  // std::vector<std::pair<unsigned int, Real>> vec_rval;
376 
377  // start solver loop
378  for (unsigned int n=n_start; n<= n_stop; n++)
379  {
380  // set the current frequency
381  this->set_current_frequency(n);
382 
383  // Call the user-supplied pre-solve method
384  START_LOG("user_pre_solve()", "FrequencySystem");
385 
386  this->solve_system (es, this->name());
387 
388  STOP_LOG("user_pre_solve()", "FrequencySystem");
389 
390 
391  // Solve the linear system for this specific frequency
392  const std::pair<unsigned int, Real> rval =
393  linear_solver->solve (*matrix, *solution, *rhs, tol, maxits);
394 
395  _n_linear_iterations = rval.first;
396  _final_linear_residual = rval.second;
397 
398  vec_rval.push_back(rval);
399 
403  if (this->_keep_solution_duplicates)
404  this->get_vector(this->form_solu_vec_name(n)) = *solution;
405  }
406 
407  // sanity check
408  //libmesh_assert_equal_to (vec_rval.size(), (n_stop-n_start+1));
409 
410  //return vec_rval;
411 }
412 
413 
414 
416  const std::string & name))
417 {
419 
420  solve_system = fptr;
421 }
422 
423 
424 
426 {
427  libmesh_assert_less (n, n_frequencies());
428 
429  EquationSystems & es =
430  this->get_equation_systems();
431 
432  es.parameters.set<Number>("current frequency") =
433  es.parameters.get<Number>(this->form_freq_param_name(n));
434 }
435 
436 
437 
438 std::string FrequencySystem::form_freq_param_name(const unsigned int n) const
439 {
440  libmesh_assert_less (n, 9999);
441  char buf[15];
442  sprintf(buf, "frequency %04u", n);
443  return (buf);
444 }
445 
446 
447 
448 std::string FrequencySystem::form_solu_vec_name(const unsigned int n) const
449 {
450  libmesh_assert_less (n, 9999);
451  char buf[15];
452  sprintf(buf, "solution %04u", n);
453  return (buf);
454 }
455 
456 } // namespace libMesh
457 
458 
459 #endif // if defined(LIBMESH_USE_COMPLEX_NUMBERS)
std::string form_solu_vec_name(const unsigned int n) const
This is the EquationSystems class.
void clear_all()
The full clear method also clears the frequencies (stored as parameters of the EquationSystems object...
void(* solve_system)(EquationSystems &es, const std::string &name)
Function that computes frequency-dependent data of the system.
virtual void solve() libmesh_override
Solves the system for all frequencies.
unsigned int _n_linear_iterations
The number of linear iterations required to solve the linear system Ax=b.
This class provides a specific system class.
bool _keep_solution_duplicates
when the solution for each frequency should be stored in an additional vector, then this bool is true...
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
void set_frequencies_by_range(const Number min_freq, const Number max_freq, const unsigned int n_freq, const bool allocate_solution_duplicates=true)
Set the frequency range for which the system should be solved.
const T & get(const std::string &) const
Definition: parameters.h:430
void set_frequencies(const std::vector< Real > &frequencies, const bool allocate_solution_duplicates=true)
Set the frequency range by simply copying the values from frequencies.
virtual void clear() libmesh_override
Clear all the data structures associated with the system, but leave the frequencies untouched...
The libMesh namespace provides an interface to certain functionality in the library.
FrequencySystem(EquationSystems &es, const std::string &name_in, const unsigned int number_in)
Constructor.
std::string form_freq_param_name(const unsigned int n) const
const std::string & name() const
Definition: system.h:1998
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
Number fptr(const Point &p, const Parameters &, const std::string &libmesh_dbg_var(sys_name), const std::string &unknown_name)
Definition: projection.C:78
libmesh_assert(j)
void attach_solve_function(void fptr(EquationSystems &es, const std::string &name))
Register a required user function to use in assembling/solving the system.
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
bool _finished_assemble
true when we have finished the assemble() phase.
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:681
~FrequencySystem()
Destructor.
virtual void assemble() libmesh_override
Prepares matrix and _dof_map for matrix assembly.
bool _finished_init
true when we have finished the init() phase.
void remove(const std::string &)
Removes the specified parameter from the list, if it exists.
Definition: parameters.h:480
void set_frequencies_by_steps(const Number base_freq, const Number freq_step=0., const unsigned int n_freq=1, const bool allocate_solution_duplicates=true)
Set the frequency range for which the system should be solved.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:794
unsigned int n_frequencies() const
virtual void assemble() libmesh_override
Assemble the linear system.
Real _final_linear_residual
The final residual for the linear system Ax=b.
virtual void init_data() libmesh_override
Initializes new data members of the system.
UniquePtr< LinearSolver< Number > > linear_solver
The LinearSolver defines the interface used to solve the linear_implicit system.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
T & set(const std::string &)
Definition: parameters.h:469
SparseMatrix< Number > * matrix
The system matrix.
const EquationSystems & get_equation_systems() const
Definition: system.h:712
bool _finished_set_frequencies
true when we have frequencies to solve for.
Parameters parameters
Data structure holding arbitrary parameters.
void set_current_frequency(unsigned int n)
Sets the current frequency to the n-th entry in the vector _frequencies.
bool have_parameter(const std::string &) const
Definition: parameters.h:411
std::vector< std::pair< unsigned int, Real > > vec_rval
The number of iterations and the final residual when the Ax=b is solved for multiple frequencies...