frequency_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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
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,
48  const unsigned int number) :
49  LinearImplicitSystem (es, name, number),
50  solve_system (NULL),
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 =
102  this->get_equation_systems();
103 
104  // clear frequencies in the parameters section of the
105  // EquationSystems object
106  if (es.parameters.have_parameter<unsigned int> ("n_frequencies"))
107  {
108  unsigned int n_frequencies = es.parameters.get<unsigned int>("n_frequencies");
109  for (unsigned int n=0; n < n_frequencies; n++)
110  es.parameters.remove(this->form_freq_param_name(n));
111  es.parameters.remove("current frequency");
112  }
113 }
114 
115 
116 
117 
119 {
120  // initialize parent data and additional solution vectors
122 
123  // Log how long initializing the system takes
124  START_LOG("init()", "FrequencySystem");
125 
126  EquationSystems& es =
127  this->get_equation_systems();
128 
129  // make sure we have frequencies to solve for
131  {
132  /*
133  * when this system was read from file, check
134  * if this has a "n_frequencies" parameter,
135  * and initialize us with these.
136  */
137  if (es.parameters.have_parameter<unsigned int> ("n_frequencies"))
138  {
139  const unsigned int n_freq =
140  es.parameters.get<unsigned int>("n_frequencies");
141 
142  libmesh_assert_greater (n_freq, 0);
143 
145 
146  this->set_current_frequency(0);
147  }
148  else
149  {
150  libMesh::err << "ERROR: Need to set frequencies before calling init(). " << std::endl;
151  libmesh_error();
152  }
153  }
154 
155  _finished_init = true;
156 
157  // Stop logging init()
158  STOP_LOG("init()", "FrequencySystem");
159 }
160 
161 
162 
164 {
166 
167  if (_finished_assemble)
168  {
169  libMesh::err << "ERROR: Matrices already assembled." << std::endl;
170  libmesh_error ();
171  }
172 
173  // Log how long assemble() takes
174  START_LOG("assemble()", "FrequencySystem");
175 
176  // prepare matrix with the help of the _dof_map,
177  // fill with sparsity pattern, initialize the
178  // additional matrices
180 
181  //matrix.print ();
182  //rhs.print ();
183 
184  _finished_assemble = true;
185 
186  // Log how long assemble() takes
187  STOP_LOG("assemble()", "FrequencySystem");
188 }
189 
190 
191 
193  const Real freq_step,
194  const unsigned int n_freq,
195  const bool allocate_solution_duplicates)
196 {
197  this->_keep_solution_duplicates = allocate_solution_duplicates;
198 
199  // sanity check
201  {
202  libMesh::err << "ERROR: frequencies already initialized."
203  << std::endl;
204  libmesh_error();
205  }
206 
207  EquationSystems& es =
208  this->get_equation_systems();
209 
210  // store number of frequencies as parameter
211  es.parameters.set<unsigned int>("n_frequencies") = n_freq;
212 
213  for (unsigned int n=0; n<n_freq; n++)
214  {
215  // remember frequencies as parameters
216  es.parameters.set<Real>(this->form_freq_param_name(n)) =
217  base_freq + n * freq_step;
218 
219  // build storage for solution vector, if wanted
220  if (this->_keep_solution_duplicates)
221  this->add_vector(this->form_solu_vec_name(n));
222  }
223 
225 
226  // set the current frequency
227  this->set_current_frequency(0);
228 }
229 
230 
231 
233  const Real max_freq,
234  const unsigned int n_freq,
235  const bool allocate_solution_duplicates)
236 {
237  this->_keep_solution_duplicates = allocate_solution_duplicates;
238 
239  // sanity checks
240  libmesh_assert_greater (max_freq, min_freq);
241  libmesh_assert_greater (n_freq, 0);
242 
244  {
245  libMesh::err << "ERROR: frequencies already initialized. " << std::endl;
246  libmesh_error();
247  }
248 
249  EquationSystems& es =
250  this->get_equation_systems();
251 
252  // store number of frequencies as parameter
253  es.parameters.set<unsigned int>("n_frequencies") = n_freq;
254 
255  // set frequencies, build solution storage
256  for (unsigned int n=0; n<n_freq; n++)
257  {
258  // remember frequencies as parameters
259  es.parameters.set<Real>(this->form_freq_param_name(n)) =
260  min_freq + n*(max_freq-min_freq)/(n_freq-1);
261 
262  // build storage for solution vector, if wanted
263  if (this->_keep_solution_duplicates)
265  }
266 
268 
269  // set the current frequency
270  this->set_current_frequency(0);
271 }
272 
273 
274 
275 void FrequencySystem::set_frequencies (const std::vector<Real>& frequencies,
276  const bool allocate_solution_duplicates)
277 {
278  this->_keep_solution_duplicates = allocate_solution_duplicates;
279 
280  // sanity checks
281  libmesh_assert(!frequencies.empty());
282 
284  {
285  libMesh::err << "ERROR: frequencies already initialized. " << std::endl;
286  libmesh_error();
287  }
288 
289  EquationSystems& es =
290  this->get_equation_systems();
291 
292  // store number of frequencies as parameter
293  es.parameters.set<unsigned int>("n_frequencies") = frequencies.size();
294 
295  // set frequencies, build solution storage
296  for (unsigned int n=0; n<frequencies.size(); n++)
297  {
298  // remember frequencies as parameters
299  es.parameters.set<Real>(this->form_freq_param_name(n)) = frequencies[n];
300 
301  // build storage for solution vector, if wanted
302  if (this->_keep_solution_duplicates)
304  }
305 
307 
308  // set the current frequency
309  this->set_current_frequency(0);
310 }
311 
312 
313 
314 
315 unsigned int FrequencySystem::n_frequencies () const
316 {
318  return this->get_equation_systems().parameters.get<unsigned int>("n_frequencies");
319 }
320 
321 
322 
324 {
326 
327  // Solve for all the specified frequencies
328  this->solve (0, this->n_frequencies()-1);
329 }
330 
331 
332 
333 void FrequencySystem::solve (const unsigned int n_start,
334  const unsigned int n_stop)
335 {
336  // Assemble the linear system, if not already done
337  if (!_finished_assemble)
338  this->assemble ();
339 
340  // the user-supplied solve method _has_ to be provided by the user
342 
343  // existence & range checks
345  libmesh_assert_less (n_stop, this->n_frequencies());
346 
347  EquationSystems& es =
348  this->get_equation_systems();
349 
350  // Get the user-specified linear solver tolerance,
351  // the user-specified maximum # of linear solver iterations,
352  // the user-specified wave speed
353  const Real tol =
354  es.parameters.get<Real>("linear solver tolerance");
355  const unsigned int maxits =
356  es.parameters.get<unsigned int>("linear solver maximum iterations");
357 
358 
359 // // return values
360 // std::vector< std::pair<unsigned int, Real> > vec_rval;
361 
362  // start solver loop
363  for (unsigned int n=n_start; n<= n_stop; n++)
364  {
365  // set the current frequency
366  this->set_current_frequency(n);
367 
368  // Call the user-supplied pre-solve method
369  START_LOG("user_pre_solve()", "FrequencySystem");
370 
371  this->solve_system (es, this->name());
372 
373  STOP_LOG("user_pre_solve()", "FrequencySystem");
374 
375 
376  // Solve the linear system for this specific frequency
377  const std::pair<unsigned int, Real> rval =
378  linear_solver->solve (*matrix, *solution, *rhs, tol, maxits);
379 
380  _n_linear_iterations = rval.first;
381  _final_linear_residual = rval.second;
382 
383  vec_rval.push_back(rval);
384 
388  if (this->_keep_solution_duplicates)
389  this->get_vector(this->form_solu_vec_name(n)) = *solution;
390  }
391 
392  // sanity check
393  //libmesh_assert_equal_to (vec_rval.size(), (n_stop-n_start+1));
394 
395  //return vec_rval;
396 }
397 
398 
399 
401  const std::string& name))
402 {
403  libmesh_assert(fptr);
404 
405  solve_system = fptr;
406 }
407 
408 
409 
411 {
412  libmesh_assert_less (n, n_frequencies());
413 
414  EquationSystems& es =
415  this->get_equation_systems();
416 
417  es.parameters.set<Real>("current frequency") =
418  es.parameters.get<Real>(this->form_freq_param_name(n));
419 }
420 
421 
422 
423 std::string FrequencySystem::form_freq_param_name(const unsigned int n) const
424 {
425  libmesh_assert_less (n, 9999);
426  char buf[15];
427  sprintf(buf, "frequency %04d", n);
428  return (buf);
429 }
430 
431 
432 
433 std::string FrequencySystem::form_solu_vec_name(const unsigned int n) const
434 {
435  libmesh_assert_less (n, 9999);
436  char buf[15];
437  sprintf(buf, "solution %04d", n);
438  return (buf);
439 }
440 
441 } // namespace libMesh
442 
443 
444 #endif // if defined(LIBMESH_USE_COMPLEX_NUMBERS)

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo