spacetime_navier_stokes_block_preconditioner.h
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2024 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Header file for SpaceTimeNavierStokes elements
27 #ifndef OOMPH_SPACETIME_NAVIER_STOKES_BLOCK_PRECONDITIONERS_HEADER
28 #define OOMPH_SPACETIME_NAVIER_STOKES_BLOCK_PRECONDITIONERS_HEADER
29 
30 // Config header generated by autoconfig
31 #ifdef HAVE_CONFIG_H
32 #include <oomph-lib-config.h>
33 #endif
34 
35 // Oomph-lib headers
36 #include "generic.h"
37 
38 // Space-time block preconditioning machinery
39 #include "../../SpaceTimeBlockPreconditioner/spacetime_block_preconditioner.cc"
40 
41 /// /////////////////////////////////////////////////////////////////////////
42 /// /////////////////////////////////////////////////////////////////////////
43 /// /////////////////////////////////////////////////////////////////////////
44 
45 namespace oomph
46 {
47  //=============================================================================
48  /// General purpose block triangular preconditioner. By default this is
49  /// Upper triangular. Also, by default SuperLUPreconditioner (or
50  /// SuperLUDistPreconditioner) is used to solve the subsidiary systems, but
51  /// other preconditioners can be used by setting them using passing a pointer
52  /// to a function of type SubsidiaryPreconditionerFctPt to the method
53  /// subsidiary_preconditioner_function_pt().
54  //=============================================================================
55  class SpaceTimeNavierStokesSubsidiaryPreconditioner
56  : public BlockPreconditioner<CRDoubleMatrix>
57  {
58  public:
59  /// Constructor. (By default this preconditioner is upper triangular).
62  {
63  // By default, don't store the memory statistics of this preconditioner
65 
66  // Initialise the value of Memory_usage_in_bytes
68 
69  // Flag to indicate that the preconditioner has been setup
70  // previously -- if setup() is called again, data can
71  // be wiped.
73 
74  // By default we use SuperLU for both p and f blocks
77 
78  // Set default preconditioners (inexact solvers) -- they are
79  // members of this class!
82 
83  // Null the momentum block matrix-vector product helper
84  F_mat_vec_pt = 0;
85 
86  // Null the gradient block matrix-vector product helper
87  G_mat_vec_pt = 0;
88  }
89 
90  /// Destructor - delete the preconditioner matrices
92  {
93  // Call the auxiliary clean up function
94  this->clean_up_memory();
95  } // End of ~SpaceTimeNavierStokesSubsidiaryPreconditioner
96 
97  /// Clean up the memory
98  virtual void clean_up_memory()
99  {
100  // If we've actually set the preconditioner up
102  {
103  // Delete matvecs
104  delete F_mat_vec_pt;
105  F_mat_vec_pt = 0;
106 
107  delete G_mat_vec_pt;
108  G_mat_vec_pt = 0;
109 
110  // Delete stuff from velocity solve
112  {
113  delete F_preconditioner_pt;
115  }
116 
117  // Delete stuff from Schur complement approx
119  {
120  delete P_preconditioner_pt;
122  }
123  } // if (Preconditioner_has_been_setup)
124  } // End of clean_up_memory
125 
126  /// Broken copy constructor
129 
130  /// Broken assignment operator
132  delete;
133 
134  /// For some reason we need to remind the compiler that there is
135  /// also a function named setup in the base class.
136  using Preconditioner::setup;
137 
138  /// Setup the preconditioner
139  void setup();
140 
141  /// Apply preconditioner to r
143 
144 
145  /// Document the memory usage
147  {
148  /// Set the appropriate flag to true
150  } // End of enable_doc_memory_usage
151 
152 
153  /// Don't document the memory usage!
155  {
156  /// Set the appropriate flag to false
158  } // End of disable_doc_memory_usage
159 
160 
161  /// Get the memory statistics
163  {
164  // Has the preconditioner even been set up yet?
166  {
167  // Were we meant to compute the statistics?
169  {
170  // Return the appropriate variable value
171  return Memory_usage_in_bytes;
172  }
173  else
174  {
175  // Allocate storage for an output stream
176  std::ostringstream warning_message_stream;
177 
178  // Create a warning message
179  warning_message_stream
180  << "The memory statistics have not been calculated "
181  << "so I'm returning\nthe value zero." << std::endl;
182 
183  // Give the user a warning
184  OomphLibWarning(warning_message_stream.str(),
185  OOMPH_CURRENT_FUNCTION,
186  OOMPH_EXCEPTION_LOCATION);
187 
188  // Return the value zero
189  return 0.0;
190  }
191  }
192  // If the preconditioner hasn't been set up yet
193  else
194  {
195  // Allocate storage for an output stream
196  std::ostringstream warning_message_stream;
197 
198  // Create a warning message
199  warning_message_stream
200  << "The preconditioner hasn't even been set up yet "
201  << "so I'm returning\nthe value zero." << std::endl;
202 
203  // Give the user a warning
204  OomphLibWarning(warning_message_stream.str(),
205  OOMPH_CURRENT_FUNCTION,
206  OOMPH_EXCEPTION_LOCATION);
207 
208  // Return the value zero
209  return 0.0;
210  } // if (Preconditioner_has_been_setup)
211  } // End of get_memory_usage_in_bytes
212 
213  private:
214  /// Pointer to the 'preconditioner' for the F matrix
216 
217  /// Pointer to the 'preconditioner' for the pressure matrix
219 
220  /// Flag indicating whether the default F preconditioner is used
222 
223  /// Flag indicating whether the default P preconditioner is used
225 
226  /// Control flag is true if the preconditioner has been setup
227  /// (used so we can wipe the data when the preconditioner is called again)
229 
230  /// Flag to indicate whether or not to record the memory statistics
231  /// this preconditioner
233 
234  /// Storage for the memory usage of the solver if the flag above
235  /// is set to true (in bytes)
236  double Memory_usage_in_bytes;
237 
238  /// MatrixVectorProduct operator for F
240 
241  /// MatrixVectorProduct operator for G
243  };
244 
245 
246  //=============================================================================
247  /// The block preconditioner form of GMRES. This version extracts
248  /// the blocks from the global systems and assembles the system by
249  /// concatenating all the matrices together
250  //=============================================================================
251  class GMRESBlockPreconditioner
252  : public IterativeLinearSolver,
253  public virtual BlockPreconditioner<CRDoubleMatrix>
254  {
255  public:
256  /// Constructor (empty)
259  Matrix_pt(0),
261  Iterations(0),
263  Preconditioner_LHS(false)
264  {
265  // By default, don't store the memory statistics of this preconditioner
267 
268  // Initialise the value of Memory_usage_in_bytes
269  Memory_usage_in_bytes = 0.0;
270  } // End of GMRESBlockPreconditioner
271 
272  /// Destructor
274  {
275  // Call the auxiliary clean up function
276  this->clean_up_memory();
277  } // End of ~GMRESBlockPreconditioner
278 
279  /// Clean up the memory (empty for now...)
280  virtual void clean_up_memory()
281  {
282  // If the block preconditioner has been set up previously
284  {
285  // Delete the subsidiary preconditioner
287 
288  // Make it a null pointer
290 
291  // Delete the matrix!
292  delete Matrix_pt;
293 
294  // Make it a null pointer
295  Matrix_pt = 0;
296  }
297  } // End of clean_up_memory
298 
299  /// Broken copy constructor
301 
302  /// Broken assignment operator
303  void operator=(const GMRESBlockPreconditioner&) = delete;
304 
305  /// For some reason we need to remind the compiler that there is
306  /// also a function named setup in the base class.
307  using Preconditioner::setup;
308 
309  /// Setup the preconditioner
310  void setup();
311 
312  /// Apply preconditioner to r
314 
315  /// Solver: Takes pointer to problem and returns the results vector
316  /// which contains the solution of the linear system defined by
317  /// the problem's fully assembled Jacobian and residual vector.
318  void solve(Problem* const& problem_pt, DoubleVector& result)
319  {
320  // Broken
321  throw OomphLibError("Can't use this interface!",
322  OOMPH_CURRENT_FUNCTION,
323  OOMPH_EXCEPTION_LOCATION);
324  } // End of solve
325 
326  /// Handle to the number of iterations taken
327  unsigned iterations() const
328  {
329  // Return the number of iterations
330  return Iterations;
331  } // End of iterations
332 
333  /// Set left preconditioning (the default)
335  {
336  Preconditioner_LHS = true;
337  }
338 
339  /// Enable right preconditioning
341  {
342  Preconditioner_LHS = false;
343  }
344 
345 
346  /// Document the memory usage
348  {
349  /// Set the appropriate flag to true
351  } // End of enable_doc_memory_statistics
352 
353 
354  /// Don't document the memory usage!
356  {
357  /// Set the appropriate flag to false
359  } // End of disable_doc_memory_statistics
360 
361 
362  /// Get the memory statistics
364  {
365  // Has the preconditioner even been set up yet?
367  {
368  // Were we meant to compute the statistics?
370  {
371  // Return the appropriate variable value
372  return Memory_usage_in_bytes;
373  }
374  else
375  {
376  // Allocate storage for an output stream
377  std::ostringstream warning_message_stream;
378 
379  // Create a warning message
380  warning_message_stream
381  << "The memory statistics have not been calculated "
382  << "so I'm returning\nthe value zero." << std::endl;
383 
384  // Give the user a warning
385  OomphLibWarning(warning_message_stream.str(),
386  OOMPH_CURRENT_FUNCTION,
387  OOMPH_EXCEPTION_LOCATION);
388 
389  // Return the value zero
390  return 0.0;
391  }
392  }
393  // If the preconditioner hasn't been set up yet
394  else
395  {
396  // Allocate storage for an output stream
397  std::ostringstream warning_message_stream;
398 
399  // Create a warning message
400  warning_message_stream
401  << "The preconditioner hasn't even been set up yet "
402  << "so I'm returning\nthe value zero." << std::endl;
403 
404  // Give the user a warning
405  OomphLibWarning(warning_message_stream.str(),
406  OOMPH_CURRENT_FUNCTION,
407  OOMPH_EXCEPTION_LOCATION);
408 
409  // Return the value zero
410  return 0.0;
411  } // if (Preconditioner_has_been_setup)
412  } // End of get_memory_usage_in_bytes
413 
414 
415  /// Handle to the Navier-Stokes subsidiary block preconditioner
416  /// DRAIG: Make sure the desired const-ness is correct later...
418  const
419  {
420  // Return a pointer to the appropriate member data
422  } // End of navier_stokes_subsidiary_preconditioner_pt
423 
424  private:
425  /// Helper function to update the result vector using the result,
426  /// x=x_0+V_m*y
427  void update(const unsigned& k,
428  const Vector<Vector<double>>& H,
429  const Vector<double>& s,
430  const Vector<DoubleVector>& v,
431  const DoubleVector& block_x_with_size_of_full_x,
432  DoubleVector& x)
433  {
434  // Make a local copy of s
435  Vector<double> y(s);
436 
437  // Backsolve:
438  for (int i = int(k); i >= 0; i--)
439  {
440  // Divide the i-th entry of y by the i-th diagonal entry of H
441  y[i] /= H[i][i];
442 
443  // Loop over the previous values of y and update
444  for (int j = i - 1; j >= 0; j--)
445  {
446  // Update the j-th entry of y
447  y[j] -= H[i][j] * y[i];
448  }
449  } // for (int i=int(k);i>=0;i--)
450 
451  // Store the number of rows in the result vector
452  unsigned n_x = x.nrow();
453 
454  // Build a temporary vector with entries initialised to 0.0
455  DoubleVector temp(x.distribution_pt(), 0.0);
456 
457  // Build a temporary vector with entries initialised to 0.0
458  DoubleVector z(x.distribution_pt(), 0.0);
459 
460  // Get access to the underlying values
461  double* temp_pt = temp.values_pt();
462 
463  // Calculate x=Vy
464  for (unsigned j = 0; j <= k; j++)
465  {
466  // Get access to j-th column of v
467  const double* vj_pt = v[j].values_pt();
468 
469  // Loop over the entries of the vector, temp
470  for (unsigned i = 0; i < n_x; i++)
471  {
472  temp_pt[i] += vj_pt[i] * y[j];
473  }
474  } // for (unsigned j=0;j<=k;j++)
475 
476  // If we're using LHS preconditioning
477  if (Preconditioner_LHS)
478  {
479  // Since we're using LHS preconditioning the preconditioner is applied
480  // to the matrix and RHS vector so we simply update the value of x
481  x += temp;
482  }
483  // If we're using RHS preconditioning
484  else
485  {
486  // This is pretty inefficient but there is no other way to do this with
487  // block sub preconditioners as far as I can tell because they demand
488  // to have the full r and z vectors...
489  DoubleVector block_z_with_size_of_full_z(
490  block_x_with_size_of_full_x.distribution_pt());
491  DoubleVector temp_with_size_of_full_rhs(
492  block_x_with_size_of_full_x.distribution_pt());
493 
494  // Reorder the entries of temp and put them into the global-sized vector
495  this->return_block_vector(0, temp, temp_with_size_of_full_rhs);
496 
497  // Solve on the global-sized vectors
499  temp_with_size_of_full_rhs, block_z_with_size_of_full_z);
500 
501  // Now reorder the entries and put them into z
502  this->get_block_vector(0, block_z_with_size_of_full_z, z);
503 
504  // Since we're using RHS preconditioning the preconditioner is applied
505  // to the solution vector
506  // preconditioner_pt()->preconditioner_solve(temp,z);
507 
508  // Use the update: x_m=x_0+inv(M)Vy [see Saad Y,"Iterative methods for
509  // sparse linear systems", p.284]
510  x += z;
511  }
512  } // End of update
513 
514  /// Helper function: Generate a plane rotation. This is done by
515  /// finding the values of \f$ \cos(\theta) \f$ (i.e. cs) and \sin(\theta)
516  /// (i.e. sn) such that:
517  /// \f[ \begin{bmatrix} \cos\theta & \sin\theta \newline -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \newline dy \end{bmatrix} = \begin{bmatrix} r \newline 0 \end{bmatrix}, \f]
518  /// where \f$ r=\sqrt{pow(dx,2)+pow(dy,2)} \f$. The values of a and b are
519  /// given by:
520  /// \f[ \cos\theta=\dfrac{dx}{\sqrt{pow(dx,2)+pow(dy,2)}}, \f]
521  /// and
522  /// \f[ \sin\theta=\dfrac{dy}{\sqrt{pow(dx,2)+pow(dy,2)}}. \f]
523  /// Taken from: Saad Y."Iterative methods for sparse linear systems", p.192
524  void generate_plane_rotation(double& dx, double& dy, double& cs, double& sn)
525  {
526  // If dy=0 then we do not need to apply a rotation
527  if (dy == 0.0)
528  {
529  // Using theta=0 gives cos(theta)=1
530  cs = 1.0;
531 
532  // Using theta=0 gives sin(theta)=0
533  sn = 0.0;
534  }
535  // If dx or dy is large using the normal form of calculting cs and sn
536  // is naive since this may overflow or underflow so instead we calculate
537  // r=sqrt(pow(dx,2)+pow(dy,2)) by using r=|dy|sqrt(1+pow(dx/dy,2)) if
538  // |dy|>|dx| [see <A
539  // HREF=https://en.wikipedia.org/wiki/Hypot">Hypot</A>.].
540  else if (fabs(dy) > fabs(dx))
541  {
542  // Since |dy|>|dx| calculate the ratio dx/dy
543  double temp = dx / dy;
544 
545  // Calculate sin(theta)=dy/sqrt(pow(dx,2)+pow(dy,2))
546  sn = 1.0 / sqrt(1.0 + temp * temp);
547 
548  // Calculate cos(theta)=dx/sqrt(pow(dx,2)+pow(dy,2))=(dx/dy)*sin(theta)
549  cs = temp * sn;
550  }
551  // Otherwise, we have |dx|>=|dy| so to, again, avoid overflow or underflow
552  // calculate the values of cs and sn using the method above
553  else
554  {
555  // Since |dx|>=|dy| calculate the ratio dy/dx
556  double temp = dy / dx;
557 
558  // Calculate cos(theta)=dx/sqrt(pow(dx,2)+pow(dy,2))
559  cs = 1.0 / sqrt(1.0 + temp * temp);
560 
561  // Calculate sin(theta)=dy/sqrt(pow(dx,2)+pow(dy,2))=(dy/dx)*cos(theta)
562  sn = temp * cs;
563  }
564  } // End of generate_plane_rotation
565 
566  /// Helper function: Apply plane rotation. This is done using the
567  /// update:
568  /// \f[ \begin{bmatrix} dx \newline dy \end{bmatrix} \leftarrow \begin{bmatrix} \cos\theta & \sin\theta \newline -\sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} dx \newline dy \end{bmatrix}. \f]
569  void apply_plane_rotation(double& dx, double& dy, double& cs, double& sn)
570  {
571  // Calculate the value of dx but don't update it yet
572  double temp = cs * dx + sn * dy;
573 
574  // Set the value of dy
575  dy = -sn * dx + cs * dy;
576 
577  // Set the value of dx using the correct values of dx and dy
578  dx = temp;
579  } // End of apply_plane_rotation
580 
581  /// Pointer to matrix
583 
584  /// Pointer to the preconditioner for the block matrix
587 
588  /// Number of iterations taken
589  unsigned Iterations;
590 
591  /// Flag to indicate whether or not to record the memory statistics
592  /// this preconditioner
594 
595  /// Storage for the memory usage of the solver if the flag above
596  /// is set to true (in bytes)
597  double Memory_usage_in_bytes;
598 
599  /// Control flag is true if the preconditioner has been setup (used
600  /// so we can wipe the data when the preconditioner is called again)
602 
603  /// boolean indicating use of left hand preconditioning (if true)
604  /// or right hand preconditioning (if false)
605  bool Preconditioner_LHS;
606  };
607 } // End of namespace oomph
608 #endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
Block Preconditioner base class. The block structure of the overall problem is determined from the Me...
void return_block_vector(const unsigned &n, const DoubleVector &b, DoubleVector &v) const
Takes the n-th block ordered vector, b, and copies its entries to the appropriate entries in the natu...
void get_block_vector(const unsigned &n, const DoubleVector &v, DoubleVector &b) const
Takes the naturally ordered vector, v and returns the n-th block vector, b. Here n is the block numbe...
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
unsigned nrow() const
access function to the number of global rows.
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
double * values_pt()
access function to the underlying values
The block preconditioner form of GMRES. This version extracts the blocks from the global systems and ...
void operator=(const GMRESBlockPreconditioner &)=delete
Broken assignment operator.
double Memory_usage_in_bytes
Storage for the memory usage of the solver if the flag above is set to true (in bytes)
void apply_plane_rotation(double &dx, double &dy, double &cs, double &sn)
Helper function: Apply plane rotation. This is done using the update:
unsigned iterations() const
Handle to the number of iterations taken.
GMRESBlockPreconditioner(const GMRESBlockPreconditioner &)=delete
Broken copy constructor.
void disable_doc_memory_statistics()
Don't document the memory usage!
virtual void clean_up_memory()
Clean up the memory (empty for now...)
SpaceTimeNavierStokesSubsidiaryPreconditioner * Navier_stokes_subsidiary_preconditioner_pt
Pointer to the preconditioner for the block matrix.
bool Compute_memory_statistics
Flag to indicate whether or not to record the memory statistics this preconditioner.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to r.
bool Preconditioner_has_been_setup
Control flag is true if the preconditioner has been setup (used so we can wipe the data when the prec...
SpaceTimeNavierStokesSubsidiaryPreconditioner * navier_stokes_subsidiary_preconditioner_pt() const
Handle to the Navier-Stokes subsidiary block preconditioner DRAIG: Make sure the desired const-ness i...
void set_preconditioner_RHS()
Enable right preconditioning.
bool Preconditioner_LHS
boolean indicating use of left hand preconditioning (if true) or right hand preconditioning (if false...
void solve(Problem *const &problem_pt, DoubleVector &result)
Solver: Takes pointer to problem and returns the results vector which contains the solution of the li...
double get_memory_usage_in_bytes()
Get the memory statistics.
void set_preconditioner_LHS()
Set left preconditioning (the default)
void setup()
Setup the preconditioner.
void update(const unsigned &k, const Vector< Vector< double >> &H, const Vector< double > &s, const Vector< DoubleVector > &v, const DoubleVector &block_x_with_size_of_full_x, DoubleVector &x)
Helper function to update the result vector using the result, x=x_0+V_m*y.
void generate_plane_rotation(double &dx, double &dy, double &cs, double &sn)
Helper function: Generate a plane rotation. This is done by finding the values of (i....
Matrix vector product helper class - primarily a wrapper to Trilinos's Epetra matrix vector product m...
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
virtual void setup()=0
Setup the preconditioner. Pure virtual generic interface function.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:154
General purpose block triangular preconditioner. By default this is Upper triangular....
virtual ~SpaceTimeNavierStokesSubsidiaryPreconditioner()
Destructor - delete the preconditioner matrices.
Preconditioner * P_preconditioner_pt
Pointer to the 'preconditioner' for the pressure matrix.
double Memory_usage_in_bytes
Storage for the memory usage of the solver if the flag above is set to true (in bytes)
bool Preconditioner_has_been_setup
Control flag is true if the preconditioner has been setup (used so we can wipe the data when the prec...
bool Compute_memory_statistics
Flag to indicate whether or not to record the memory statistics this preconditioner.
void operator=(const SpaceTimeNavierStokesSubsidiaryPreconditioner &)=delete
Broken assignment operator.
SpaceTimeNavierStokesSubsidiaryPreconditioner()
Constructor. (By default this preconditioner is upper triangular).
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to r.
bool Using_default_f_preconditioner
Flag indicating whether the default F preconditioner is used.
bool Using_default_p_preconditioner
Flag indicating whether the default P preconditioner is used.
Preconditioner * F_preconditioner_pt
Pointer to the 'preconditioner' for the F matrix.
SpaceTimeNavierStokesSubsidiaryPreconditioner(const SpaceTimeNavierStokesSubsidiaryPreconditioner &)=delete
Broken copy constructor.
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...