navier_stokes_preconditioners.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-2023 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 #ifndef OOMPH_NAVIER_STOKES_PRECONDITIONERS_HEADER
27 #define OOMPH_NAVIER_STOKES_PRECONDITIONERS_HEADER
28 
29 
30 // Config header generated by autoconfig
31 #ifdef HAVE_CONFIG_H
32 #include <oomph-lib-config.h>
33 #endif
34 
35 
36 // oomphlib headers
37 #include "../generic/matrices.h"
38 #include "../generic/assembly_handler.h"
39 #include "../generic/problem.h"
40 #include "../generic/block_preconditioner.h"
41 #include "../generic/preconditioner.h"
42 #include "../generic/SuperLU_preconditioner.h"
43 #include "../generic/matrix_vector_product.h"
44 #include "navier_stokes_elements.h"
46 
47 
48 namespace oomph
49 {
50  /// ////////////////////////////////////////////////////////////////////
51  /// ////////////////////////////////////////////////////////////////////
52  /// ////////////////////////////////////////////////////////////////////
53 
54 
55  //======start_of_namespace============================================
56  /// Namespace for exact solution for pressure advection diffusion
57  /// problem
58  //====================================================================
59  namespace PressureAdvectionDiffusionValidation
60  {
61  /// Flag for solution
62  extern unsigned Flag;
63 
64  /// Peclet number -- overwrite with actual Reynolds number
65  extern double Peclet;
66 
67  /// Wind
68  extern void wind_function(const Vector<double>& x, Vector<double>& wind);
69 
70  /// Exact solution as a Vector
71  extern void get_exact_u(const Vector<double>& x, Vector<double>& u);
72 
73  /// Exact solution as a scalar
74  extern void get_exact_u(const Vector<double>& x, double& u);
75 
76  /// Source function required to make the solution above an exact solution
77  extern double source_function(const Vector<double>& x_vect);
78 
79  } // namespace PressureAdvectionDiffusionValidation
80 
81 
82  /// /////////////////////////////////////////////////////////////
83  /// /////////////////////////////////////////////////////////////
84  /// /////////////////////////////////////////////////////////////
85 
86 
87  //=============================================================
88  /// A class that is used to define the functions used to
89  /// assemble the elemental contributions to the pressure
90  /// advection diffusion problem used by the Fp preconditioner.
91  //===============================================================
93  {
94  public:
95  /// Constructor. Pass spatial dimension
96  FpPreconditionerAssemblyHandler(const unsigned& ndim) // : Ndim(ndim)
97  {
98  }
99 
100  /// Empty virtual destructor
102 
103  /// Return the contribution to the residuals of the element elem_pt
104  void get_residuals(GeneralisedElement* const& elem_pt,
105  Vector<double>& residuals)
106  {
107  unsigned n_dof = elem_pt->ndof();
108  for (unsigned i = 0; i < n_dof; i++)
109  {
110  residuals[i] = 0.0;
111  }
112 
113  dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(elem_pt)
114  ->fill_in_pressure_advection_diffusion_residuals(residuals);
115  }
116 
117  /// Calculate the elemental Jacobian matrix "d equation
118  /// / d variable" for elem_pt.
119  void get_jacobian(GeneralisedElement* const& elem_pt,
120  Vector<double>& residuals,
121  DenseMatrix<double>& jacobian)
122  {
123  // Initialise
124  unsigned n_dof = elem_pt->ndof();
125  for (unsigned i = 0; i < n_dof; i++)
126  {
127  residuals[i] = 0.0;
128  for (unsigned j = 0; j < n_dof; j++)
129  {
130  jacobian(i, j) = 0.0;
131  }
132  }
133 
134  dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(elem_pt)
135  ->fill_in_pressure_advection_diffusion_jacobian(residuals, jacobian);
136  }
137 
138  private:
139  /// Spatial dimension of problem
140  // unsigned Ndim; // cgj: Ndim is never used.
141  };
142 
143 
144  /// ////////////////////////////////////////////////////////////////
145  /// ////////////////////////////////////////////////////////////////
146  /// ////////////////////////////////////////////////////////////////
147 
148 
149  //===========================================================================
150  /// Auxiliary Problem that can be used to assemble the
151  /// pressure advection diffusion matrix needed by the FpPreconditoner
152  //===========================================================================
153  template<class ELEMENT>
155  {
156  public:
157  /// Constructor: Pass Navier-Stokes mesh and pointer to orig problem
159  Problem* orig_problem_pt)
160  {
161  // Pass across mesh -- boundary conditions have already
162  // been applied and the equations have been enumerated.
163  mesh_pt() = navier_stokes_mesh_pt;
164 
165  // Store pointer to orig problem so we can re-assign the
166  // orig eqn numbers when we're done
167  Orig_problem_pt = orig_problem_pt;
168 
169  // Get the spatial dimension
170  Ndim = mesh_pt()->finite_element_pt(0)->dim();
171 
172  // Set new assembly handler
174  }
175 
176 
177  /// Get the pressure advection diffusion Jacobian
179  {
180  // Pin all non-pressure dofs
182 
183  // Re-assign the equation numbers for all elements in Navier-Stokes mesh
184  unsigned n_dof = assign_eqn_numbers();
185  oomph_info << "Eqn numbers after pinning veloc dofs: " << n_dof
186  << std::endl;
187 
188  // Get "Jacobian" of the modified system
189  DoubleVector dummy_residuals;
190  this->get_jacobian(dummy_residuals, jacobian);
191 
192  // Reset pin status
194 
195  // Reassign equation numbering of original problem
196  oomph_info << "Eqn numbers in orig problem after re-assignment: "
197  << Orig_problem_pt->assign_eqn_numbers() << std::endl;
198  }
199 
200 
201  /// Reset pin status of all values
203  {
204  // Reset pin status for nodes
205  unsigned nnod = mesh_pt()->nnode();
206  for (unsigned j = 0; j < nnod; j++)
207  {
208  Node* nod_pt = mesh_pt()->node_pt(j);
209  unsigned nval = nod_pt->nvalue();
210  for (unsigned i = 0; i < nval; i++)
211  {
212  nod_pt->eqn_number(i) = Eqn_number_backup[nod_pt][i];
213  }
214  }
215 
216  // ... and internal data
217  unsigned nelem = mesh_pt()->nelement();
218  for (unsigned e = 0; e < nelem; e++)
219  {
220  // Get actual element
221  ELEMENT* el_pt =
222  dynamic_cast<ELEMENT*>(mesh_pt()->finite_element_pt(e));
223 
224 
225 #ifdef PARANOID
226  if (el_pt == 0)
227  {
228  std::ostringstream error_message;
229  error_message << "Navier Stokes mesh must contain only Navier Stokes"
230  << "bulk elements\n";
231  throw OomphLibError(error_message.str(),
232  OOMPH_CURRENT_FUNCTION,
233  OOMPH_EXCEPTION_LOCATION);
234  }
235 #endif
236 
237  unsigned nint = el_pt->ninternal_data();
238  for (unsigned j = 0; j < nint; j++)
239  {
240  Data* data_pt = el_pt->internal_data_pt(j);
241  unsigned nvalue = data_pt->nvalue();
242  for (unsigned i = 0; i < nvalue; i++)
243  {
244  data_pt->eqn_number(i) = Eqn_number_backup[data_pt][i];
245  }
246  }
247  }
248 
249  // Free up storage
250  Eqn_number_backup.clear();
251  }
252 
253 
254  /// Pin all non-pressure dofs and (if boolean flag is set to true
255  /// also all pressure dofs along domain boundaries -- only used
256  /// for validation)
258  const bool& set_pressure_bc_for_validation = false)
259  {
260  // Backup pin status and then pin all non-pressure degrees of freedom
261  unsigned nelem = mesh_pt()->nelement();
262  for (unsigned e = 0; e < nelem; e++)
263  {
264  // Get actual element (needed because different elements
265  // store nodal pressure in different places)
266  ELEMENT* el_pt =
267  dynamic_cast<ELEMENT*>(mesh_pt()->finite_element_pt(e));
268 
269 
270 #ifdef PARANOID
271  if (el_pt == 0)
272  {
273  std::ostringstream error_message;
274  error_message << "Navier Stokes mesh must contain only Navier Stokes"
275  << "bulk elements\n";
276  throw OomphLibError(error_message.str(),
277  OOMPH_CURRENT_FUNCTION,
278  OOMPH_EXCEPTION_LOCATION);
279  }
280 #endif
281 
282  // Check if element has internal pressure representation
283  // usually discontinuous -- preconditioner doesn't work for that case!
284  if (el_pt->p_nodal_index_nst() < 0)
285  {
286  std::ostringstream error_message;
287  error_message << "Cannot use Fp preconditioner with discontinuous\n"
288  << "pressures.\n";
289  throw OomphLibError(error_message.str(),
290  OOMPH_CURRENT_FUNCTION,
291  OOMPH_EXCEPTION_LOCATION);
292  }
293 
294  // Loop over internal data and pin the values (having established that
295  // pressure dofs aren't amongst those)
296  unsigned nint = el_pt->ninternal_data();
297  for (unsigned j = 0; j < nint; j++)
298  {
299  Data* data_pt = el_pt->internal_data_pt(j);
300  if (Eqn_number_backup[data_pt].size() == 0)
301  {
302  unsigned nvalue = data_pt->nvalue();
303  Eqn_number_backup[data_pt].resize(nvalue);
304  for (unsigned i = 0; i < nvalue; i++)
305  {
306  // Backup
307  Eqn_number_backup[data_pt][i] = data_pt->eqn_number(i);
308 
309  // Pin everything
310  data_pt->pin(i);
311  }
312  }
313  }
314 
315  // Now deal with nodal values
316  unsigned nnod = el_pt->nnode();
317  for (unsigned j = 0; j < nnod; j++)
318  {
319  Node* nod_pt = el_pt->node_pt(j);
320  if (Eqn_number_backup[nod_pt].size() == 0)
321  {
322  unsigned nvalue = nod_pt->nvalue();
323  Eqn_number_backup[nod_pt].resize(nvalue);
324  for (unsigned i = 0; i < nvalue; i++)
325  {
326  // Pin everything apart from the nodal pressure
327  // value
328  if (int(i) != el_pt->p_nodal_index_nst())
329  {
330  // Backup and pin
331  Eqn_number_backup[nod_pt][i] = nod_pt->eqn_number(i);
332  nod_pt->pin(i);
333  }
334  // Else it's a pressure value
335  else
336  {
337  // Exclude non-nodal pressure based elements
338  if (el_pt->p_nodal_index_nst() >= 0)
339  {
340  // Backup
341  Eqn_number_backup[nod_pt][i] = nod_pt->eqn_number(i);
342  }
343  }
344  }
345  }
346 
347  // Set wind
348  if (set_pressure_bc_for_validation)
349  {
350  Vector<double> x(2);
351  x[0] = nod_pt->x(0);
352  x[1] = nod_pt->x(1);
353  Vector<double> u(2);
355  nod_pt->set_value(el_pt->u_index_nst(0), u[0]);
356  nod_pt->set_value(el_pt->u_index_nst(1), u[1]);
357  }
358  }
359  }
360  }
361 
362 
363  /// Validate pressure advection diffusion problem and doc exact
364  /// and computed solution in directory specified in DocInfo object
365  void validate(DocInfo& doc_info)
366  {
367  oomph_info << "\n\n==============================================\n\n";
368  oomph_info << "Doing validation for pressure adv diff problem\n";
369 
370  // Choose exact solution
372 
373  // Pin all non-pressure dofs and pressure dofs along boundaries for
374  // validation
375  bool set_pressure_bc_for_validation = true;
376  pin_all_non_pressure_dofs(set_pressure_bc_for_validation);
377 
378 
379  // Set Peclet number
381  dynamic_cast<NavierStokesEquations<2>*>(mesh_pt()->element_pt(0))->re();
382 
383  // Loop over all elements and set source function pointer
384  unsigned nel = mesh_pt()->nelement();
385  for (unsigned e = 0; e < nel; e++)
386  {
387  dynamic_cast<NavierStokesEquations<2>*>(mesh_pt()->element_pt(e))
388  ->source_fct_for_pressure_adv_diff() =
390  }
391 
392  // Re-assign the equation numbers for all elements in Navier-Stokes mesh
393  oomph_info << "Eqn numbers after pinning veloc dofs: "
394  << assign_eqn_numbers() << std::endl;
395 
396  // Attach Robin BC elements
397  unsigned nbound = mesh_pt()->nboundary();
399  {
400  // Loop over all boundaries of Navier Stokes mesh
401  for (unsigned b = 0; b < nbound; b++)
402  {
403  // How many bulk elements are adjacent to boundary b?
404  unsigned n_element = mesh_pt()->nboundary_element(b);
405 
406  // Loop over the bulk elements adjacent to boundary b?
407  for (unsigned e = 0; e < n_element; e++)
408  {
411  mesh_pt()->boundary_element_pt(b, e));
412 
413  // What is the index of the face of the bulk element e on bondary b
414  int face_index = mesh_pt()->face_index_at_boundary(b, e);
415 
416  // Build face element
417  bulk_elem_pt->build_fp_press_adv_diff_robin_bc_element(face_index);
418 
419  } // end of loop over bulk elements adjacent to boundary b
420  }
421  }
422 
423 
424  // Loop over all elements and set source function pointer
425  std::ofstream outfile;
426  outfile.open("robin_elements.dat");
427  for (unsigned e = 0; e < nel; e++)
428  {
429  dynamic_cast<NavierStokesEquations<2>*>(mesh_pt()->element_pt(e))
430  ->output_pressure_advection_diffusion_robin_elements(outfile);
431  }
432  outfile.close();
433 
434  // Solve it
435  newton_solve();
436 
437  // And output the solution...
438  doc_solution(doc_info);
439 
440  // Kill Robin BC elements
442  {
443  // Loop over all boundaries of Navier Stokes mesh
444  for (unsigned b = 0; b < nbound; b++)
445  {
446  // How many bulk elements are adjacent to boundary b?
447  unsigned n_element = mesh_pt()->nboundary_element(b);
448 
449  // Loop over the bulk elements adjacent to boundary b?
450  for (unsigned e = 0; e < n_element; e++)
451  {
454  mesh_pt()->boundary_element_pt(b, e));
455 
456  // Kill associated Robin elements
458 
459  } // end of loop over bulk elements adjacent to boundary b
460  }
461  }
462 
463  // Reset pin status
465 
466  // Reassign equation numbering of original problem
467  oomph_info << "Eqn numbers in orig problem after re-assignment: "
468  << Orig_problem_pt->assign_eqn_numbers() << std::endl;
469 
470 
471  oomph_info << "Done validation for pressure adv diff problem\n";
472  oomph_info << "\n\n==============================================\n\n";
473  }
474 
475 
476  /// Doc solution (only needed during development -- kept alive
477  /// in case validation is required during code maintenance)
478  void doc_solution(DocInfo& doc_info)
479  {
480  std::ofstream some_file;
481  std::ostringstream filename;
482 
483  // Number of plot points
484  unsigned npts;
485  npts = 5;
486 
487  // Check value of FE solution in first element
488  Vector<double> s(Ndim, 0.0);
489  Vector<double> x(Ndim);
490  double p = dynamic_cast<ELEMENT*>(mesh_pt()->finite_element_pt(0))
491  ->interpolated_p_nst(s);
492  dynamic_cast<ELEMENT*>(mesh_pt()->finite_element_pt(0))
493  ->interpolated_x(s, x);
494 
495  // Get offset-free exact solution
496  double u_exact = 0;
498 
499  // Adjust offset
500  unsigned nnode = mesh_pt()->nnode();
501  for (unsigned j = 0; j < nnode; j++)
502  {
503  if (mesh_pt()->node_pt(j)->nvalue() == 3)
504  {
505  *(mesh_pt()->node_pt(j)->value_pt(2)) -= p - u_exact;
506  }
507  }
508 
509  // Output solution
510  filename << doc_info.directory() << "/fp_soln" << doc_info.number()
511  << ".dat";
512  some_file.open(filename.str().c_str());
513  some_file.precision(20);
514  mesh_pt()->output(some_file, npts);
515  some_file.close();
516 
517  filename.str("");
518  filename << doc_info.directory() << "/fp_exact_soln" << doc_info.number()
519  << ".dat";
520  some_file.open(filename.str().c_str());
521  some_file.precision(20);
522  mesh_pt()->output_fct(
524  some_file.close();
525  }
526 
527  private:
528  /// The spatial dimension
529  unsigned Ndim;
530 
531  /// Pointer to orig problem (required so we can re-assign
532  /// eqn numbers)
534 
535  /// Map to store original equation numbers
536  std::map<Data*, std::vector<int>> Eqn_number_backup;
537  };
538 
539 
540  /// /////////////////////////////////////////////////////////////////////////
541  /// /////////////////////////////////////////////////////////////////////////
542  /// /////////////////////////////////////////////////////////////////////////
543 
544 
545  //===========================================================================
546  /// The least-squares commutator (LSC; formerly BFBT) Navier Stokes
547  /// preconditioner. It uses blocks corresponding to the velocity
548  /// and pressure unknowns, i.e. there are a total of 2x2 blocks,
549  /// and all velocity components are treated as a single block of unknowns.
550  ///
551  /// Here are the details: An "ideal" Navier-Stokes preconditioner
552  /// would solve the system
553  /// \f[ \left( \begin{array}{cc} {\bf F} & {\bf G} \\ {\bf D} & {\bf 0} \end{array} \right) \left( \begin{array}{c} {\bf z}_u \\ {\bf z}_p \end{array} \right) = \left( \begin{array}{c} {\bf r}_u \\ {\bf r}_p \end{array} \right) \f]
554  /// where \f$ {\bf F}\f$ is the momentum block, \f$ {\bf G} \f$ the
555  /// discrete gradient operator, and \f$ {\bf D}\f$ the discrete
556  /// divergence operator. (For unstabilised elements, we have
557  /// \f$ {\bf D} = {\bf G}^T \f$ and in much of the literature
558  /// the divergence matrix is denoted by \f$ {\bf B} \f$ .)
559  /// The use of this preconditioner would ensure the convergence
560  /// of any iterative linear solver in a single iteration but its
561  /// application is, of course, exactly as expensive as a direct solve.
562  /// The LSC/BFBT preconditioner replaces the exact Jacobian by
563  /// a block-triangular approximation
564  /// \f[ \left( \begin{array}{cc} {\bf F} & {\bf G} \\ {\bf 0} & -{\bf M}_s \end{array} \right) \left( \begin{array}{c} {\bf z}_u \\ {\bf z}_p \end{array} \right) = \left( \begin{array}{c} {\bf r}_u \\ {\bf r}_p \end{array} \right), \f]
565  /// where \f${\bf M}_s\f$ is an approximation to the pressure
566  /// Schur-complement \f$ {\bf S} = {\bf D} {\bf F}^{-1}{\bf G}. \f$
567  /// This system can be solved in two steps:
568  /// -# Solve the second row for \f$ {\bf z}_p\f$ via
569  ///
570  /// \f[ {\bf z}_p = - {\bf M}_s^{-1} {\bf r}_p \f]
571  /// -# Given \f$ {\bf z}_p \f$ , solve the first row for \f$ {\bf z}_u\f$ via
572  ///
573  /// \f[ {\bf z}_u = {\bf F}^{-1} \big( {\bf r}_u - {\bf G} {\bf z}_p \big) \f]
574  /// .
575  /// In the LSC/BFBT preconditioner, the action of the inverse pressure
576  /// Schur complement
577  /// \f[ {\bf z}_p = - {\bf M}_s^{-1} {\bf r}_p \f]
578  /// is approximated by
579  /// \f[ {\bf z}_p = - \big({\bf D} \widehat{\bf Q}^{-1}{\bf G} \big)^{-1} \big({\bf D} \widehat{\bf Q}^{-1}{\bf F} \widehat{\bf Q}^{-1}{\bf G}\big) \big({\bf D} \widehat{\bf Q}^{-1}{\bf G} \big)^{-1} {\bf r}_p, \f]
580  /// where \f$ \widehat{\bf Q} \f$ is the diagonal of the velocity
581  /// mass matrix. The evaluation of this expression involves
582  /// two linear solves involving the matrix
583  /// \f[ {\bf P} = \big({\bf D} \widehat{\bf Q}^{-1}{\bf G} \big) \f]
584  /// which has the character of a matrix arising from the discretisation
585  /// of a Poisson problem on the pressure space. We also have
586  /// to evaluate matrix-vector products with the matrix
587  /// \f[ {\bf E}={\bf D}\widehat{\bf Q}^{-1}{\bf F}\widehat{\bf Q}^{-1}{\bf G} \f]
588  /// Details of the theory can be found in "Finite Elements and
589  /// Fast Iterative Solvers with Applications in Incompressible Fluid
590  /// Dynamics" by Howard C. Elman, David J. Silvester, and Andrew J. Wathen,
591  /// published by Oxford University Press, 2006.
592  ///
593  /// In our implementation of the preconditioner, the linear systems
594  /// can either be solved "exactly", using SuperLU (in its incarnation
595  /// as an exact preconditioner; this is the default) or by any
596  /// other Preconditioner (inexact solver) specified via the access functions
597  /// \code
598  /// NavierStokesSchurComplementPreconditioner::set_f_preconditioner(...)
599  /// \endcode
600  /// or
601  /// \code
602  /// NavierStokesSchurComplementPreconditioner::set_p_preconditioner(...)
603  /// \endcode
604  //===========================================================================
606  : public BlockPreconditioner<CRDoubleMatrix>
607  {
608  public:
609  /// Constructor - sets defaults for control flags
612  {
613  // Insist that all elements are of type
614  // NavierStokesElementWithDiagonalMassMatrices and issue a warning
615  // if one is not.
617 
618  // Use Robin BC elements for Fp preconditioner -- yes by default
619  Use_robin_for_fp = true;
620 
621  // Flag to indicate that the preconditioner has been setup
622  // previously -- if setup() is called again, data can
623  // be wiped.
625 
626  // By default we use the LSC version
627  Use_LSC = true;
628 
629  // By default we use SuperLU for both p and f blocks
632 
633  // Pin pressure dof in press adv diff problem for Fp precond
635 
637 
638  // Set default preconditioners (inexact solvers) -- they are
639  // members of this class!
642 
643  // set Doc_time to false
644  Doc_time = false;
645 
646  // null the off diagonal Block matrix pt
647  Bt_mat_vec_pt = 0;
648 
649  // null the F matrix vector product helper
650  F_mat_vec_pt = 0;
651 
652  // null the QBt matrix vector product pt
653  QBt_mat_vec_pt = 0;
654 
655  // null the E matrix vector product helper in Fp
656  E_mat_vec_pt = 0;
657 
658  // Initially assume that there are no multiple element types in the
659  // Navier-Stokes mesh
661  }
662 
663  /// Destructor
665  {
666  clean_up_memory();
667  }
668 
669  /// Broken copy constructor
672 
673  /// Broken assignment operator
674  // Commented out broken assignment operator because this can lead to a
675  // conflict warning when used in the virtual inheritence hierarchy.
676  // Essentially the compiler doesn't realise that two separate
677  // implementations of the broken function are the same and so, quite
678  // rightly, it shouts.
679  /*void operator=(const NavierStokesSchurComplementPreconditioner&) =
680  * delete;*/
681 
682  /// Set the problem pointer (non-const pointer, the problem WILL be
683  /// changed) for use in get_pressure_advection_diffusion_matrix().
685  {
687  }
688 
690  {
691 #ifdef PARANOID
692  if (Problem_pt == 0)
693  {
694  std::ostringstream error_msg;
695  error_msg << "Problem pointer is null, did you set it yet?";
696  throw OomphLibError(
697  error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
698  }
699 #endif
700  return Problem_pt;
701  }
702 
703  /// Accept presence of elements that are not of type
704  /// NavierStokesElementWithDiagonalMassMatrices without issuing a warning.
706  {
708  }
709 
710 
711  /// Do not accept presence of elements that are not of type
712  /// NavierStokesElementWithDiagonalMassMatrices without issuing a warning.
714  {
716  }
717 
718 
719  /// Setup the preconditioner
720  void setup();
721 
722  /// for some reason we have to remind the compiler that there is a
723  /// setup() function in Preconditioner base class.
724  using Preconditioner::setup;
725 
726  /// Apply preconditioner to Vector r
728 
729  /// Specify the mesh containing the block-preconditionable Navier-Stokes
730  /// elements. The optional argument indicates if there are multiple types
731  /// of elements in the same mesh.
733  Mesh* mesh_pt,
734  const bool& allow_multiple_element_type_in_navier_stokes_mesh = false)
735  {
736  // Store the mesh pointer.
738 
739  // Are there multiple element types in this mesh?
741  allow_multiple_element_type_in_navier_stokes_mesh;
742  }
743 
744  /// Function to set a new pressure matrix preconditioner (inexact solver)
745  void set_p_preconditioner(Preconditioner* new_p_preconditioner_pt)
746  {
747  // If the default preconditioner has been used
748  // clean it up now...
750  {
751  delete P_preconditioner_pt;
752  }
753  P_preconditioner_pt = new_p_preconditioner_pt;
755  }
756 
757  /// Function to (re-)set pressure matrix preconditioner (inexact
758  /// solver) to SuperLU
760  {
762  {
765  }
766  }
767 
768  /// Function to set a new momentum matrix preconditioner (inexact solver)
769  void set_f_preconditioner(Preconditioner* new_f_preconditioner_pt)
770  {
771  // If the default preconditioner has been used
772  // clean it up now...
774  {
775  delete F_preconditioner_pt;
776  }
777  F_preconditioner_pt = new_f_preconditioner_pt;
779  }
780 
781  /// Use LSC version of the preconditioner
782  void use_lsc()
783  {
784  Use_LSC = true;
785  }
786 
787  /// Use Fp version of the preconditioner
788  void use_fp()
789  {
790  Use_LSC = false;
791  }
792 
793  /// Function to (re-)set momentum matrix preconditioner (inexact
794  /// solver) to SuperLU
796  {
798  {
801  }
802  }
803 
804  /// Enable documentation of time
806  {
807  Doc_time = true;
808  }
809 
810  /// Disable documentation of time
812  {
813  Doc_time = false;
814  }
815 
816  /// Helper function to delete preconditioner data.
817  void clean_up_memory();
818 
819  /// Use Robin BC elements for the Fp preconditioner
821  {
822  Use_robin_for_fp = true;
823  }
824 
825  /// Don't use Robin BC elements for the Fp preconditioenr
827  {
828  Use_robin_for_fp = false;
829  }
830 
831  /// Set boolean indicating that we want to pin first pressure
832  /// dof in Navier Stokes mesh when
833  /// assembling the pressure advection diffusion system for
834  /// Fp preconditoner -- needed at zero Reynolds number
835  /// for non-enclosed flows but seems harmless in any case
837  {
839  }
840 
841  /// Set boolean indicating that we do not want to pin first pressure
842  /// dof in Navier Stokes mesh when
843  /// assembling the pressure advection diffusion system for
844  /// Fp preconditoner -- needed at zero Reynolds number
845  /// for non-enclosed flows but seems harmless in any case
847  {
849  }
850 
851  /// Validate auxiliary pressure advection diffusion problem
852  /// in 2D
853  template<class ELEMENT>
854  void validate(DocInfo& doc_info, Problem* orig_problem_pt)
855  {
857  Navier_stokes_mesh_pt, orig_problem_pt);
858  aux_problem.validate(doc_info);
859  }
860 
861  /// Pin all non-pressure dofs
863  {
864  // Backup pin status and then pin all non-pressure degrees of freedom
865  unsigned nelem = Navier_stokes_mesh_pt->nelement();
866  for (unsigned e = 0; e < nelem; e++)
867  {
868  // Get pointer to the bulk element that is adjacent to boundary b
872 
873  // Pin
875  }
876  }
877 
878 
879  /// Get the pressure advection diffusion matrix
881  {
882  // Pin all non-pressure dofs
884 
885  // Get the spatial dimension
886  unsigned ndim = Navier_stokes_mesh_pt->finite_element_pt(0)->dim();
887 
888  // Backup assembly handler and set new one
889  AssemblyHandler* backed_up_assembly_handler_pt =
893 
894  // Backup submeshes (if any)
895  unsigned n_sub_mesh = problem_pt()->nsub_mesh();
896  Vector<Mesh*> backed_up_sub_mesh_pt(n_sub_mesh);
897  for (unsigned i = 0; i < n_sub_mesh; i++)
898  {
899  backed_up_sub_mesh_pt[i] = problem_pt()->mesh_pt(i);
900  }
901  // Flush sub-meshes but don't call rebuild_global_mesh()
902  // so we can simply re-instate the sub-meshes below.
904 
905  // Backup the problem's own mesh pointer
906  Mesh* backed_up_mesh_pt = problem_pt()->mesh_pt();
908 
909 #ifdef OOMPH_HAS_MPI
910 
911  // Backup the start and end elements for the distributed
912  // assembly process
913  Vector<unsigned> backed_up_first_el_for_assembly;
914  Vector<unsigned> backed_up_last_el_for_assembly;
916  backed_up_first_el_for_assembly, backed_up_last_el_for_assembly);
917 
918  // Now re-assign
920 
921 #endif
922 
923  // Identify pinned pressure dof
924  int pinned_pressure_eqn = -2;
926  {
927  // Loop over all Navier Stokes elements to find first non-pinned
928  // pressure dof
929  unsigned nel = Navier_stokes_mesh_pt->nelement();
930  for (unsigned e = 0; e < nel; e++)
931  {
935  int local_eqn = bulk_elem_pt->p_local_eqn(0);
936  if (local_eqn >= 0)
937  {
938  pinned_pressure_eqn = bulk_elem_pt->eqn_number(local_eqn);
939  break;
940  }
941  }
942 
943 
944 #ifdef OOMPH_HAS_MPI
945  if (problem_pt()->distributed())
946  {
947  int pinned_pressure_eqn_local = pinned_pressure_eqn;
948  int pinned_pressure_eqn_global = pinned_pressure_eqn;
949  MPI_Allreduce(&pinned_pressure_eqn_local,
950  &pinned_pressure_eqn_global,
951  1,
952  MPI_INT,
953  MPI_MIN,
954  this->comm_pt()->mpi_comm());
955  pinned_pressure_eqn = pinned_pressure_eqn_global;
956  }
957 
958 #endif
959  }
960 
961 
962  // Loop over all Navier Stokes elements
963  unsigned nel = Navier_stokes_mesh_pt->nelement();
964  for (unsigned e = 0; e < nel; e++)
965  {
969 
970  // Set pinned pressure equation
971  bulk_elem_pt->pinned_fp_pressure_eqn() = pinned_pressure_eqn;
972  }
973 
974 
975  // Attach Robin BC elements
976  unsigned nbound = Navier_stokes_mesh_pt->nboundary();
977  if (Use_robin_for_fp)
978  {
979  // Loop over all boundaries of Navier Stokes mesh
980  for (unsigned b = 0; b < nbound; b++)
981  {
982  // How many bulk elements are adjacent to boundary b?
983  unsigned n_element = Navier_stokes_mesh_pt->nboundary_element(b);
984 
985  // Loop over the bulk elements adjacent to boundary b?
986  for (unsigned e = 0; e < n_element; e++)
987  {
991 
992  // What is the index of the face of the bulk element e on bondary b
993  int face_index =
995 
996  // Build face element
997  bulk_elem_pt->build_fp_press_adv_diff_robin_bc_element(face_index);
998 
999  } // end of loop over bulk elements adjacent to boundary b
1000  }
1001  }
1002 
1003  // Get "Jacobian" of the modified system
1004  DoubleVector dummy_residuals;
1005  problem_pt()->get_jacobian(dummy_residuals, fp_matrix);
1006 
1007  // Kill Robin BC elements
1008  if (Use_robin_for_fp)
1009  {
1010  // Loop over all boundaries of Navier Stokes mesh
1011  for (unsigned b = 0; b < nbound; b++)
1012  {
1013  // How many bulk elements are adjacent to boundary b?
1014  unsigned n_element = Navier_stokes_mesh_pt->nboundary_element(b);
1015 
1016  // Loop over the bulk elements adjacent to boundary b?
1017  for (unsigned e = 0; e < n_element; e++)
1018  {
1022 
1023  // Kill associated Robin elements
1025 
1026  } // end of loop over bulk elements adjacent to boundary b
1027  }
1028  }
1029 
1030  // Reset pin status
1031  reset_pin_status();
1032 
1033 #ifdef OOMPH_HAS_MPI
1034 
1035  // Reset start and end elements for the distributed
1036  // assembly process
1038  backed_up_first_el_for_assembly, backed_up_last_el_for_assembly);
1039 
1040 #endif
1041 
1042  // Cleanup and reset assembly handler
1043  delete problem_pt()->assembly_handler_pt();
1044  problem_pt()->assembly_handler_pt() = backed_up_assembly_handler_pt;
1045 
1046  // Re-instate submeshes. (No need to call rebuild_global_mesh()
1047  // as it was never unbuilt).
1048  for (unsigned i = 0; i < n_sub_mesh; i++)
1049  {
1050  problem_pt()->add_sub_mesh(backed_up_sub_mesh_pt[i]);
1051  }
1052 
1053 
1054  // Reset the problem's mesh pointer
1055  problem_pt()->mesh_pt() = backed_up_mesh_pt;
1056  }
1057 
1058 
1059  /// Reset pin status of all values
1061  {
1062  // Reset pin status for nodes
1063  unsigned nnod = Navier_stokes_mesh_pt->nnode();
1064  for (unsigned j = 0; j < nnod; j++)
1065  {
1066  Node* nod_pt = Navier_stokes_mesh_pt->node_pt(j);
1067  unsigned nval = nod_pt->nvalue();
1068  for (unsigned i = 0; i < nval; i++)
1069  {
1070  nod_pt->eqn_number(i) = Eqn_number_backup[nod_pt][i];
1071  }
1072 
1073  // If it's a solid node deal with its positional data too
1074  SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
1075  if (solid_nod_pt != 0)
1076  {
1077  Data* solid_posn_data_pt = solid_nod_pt->variable_position_pt();
1078  unsigned nval = solid_posn_data_pt->nvalue();
1079  for (unsigned i = 0; i < nval; i++)
1080  {
1081  solid_posn_data_pt->eqn_number(i) =
1082  Eqn_number_backup[solid_posn_data_pt][i];
1083  }
1084  }
1085  }
1086 
1087  // ... and internal data
1088  unsigned nelem = Navier_stokes_mesh_pt->nelement();
1089  for (unsigned e = 0; e < nelem; e++)
1090  {
1091  // Pointer to element
1093 
1094  // Pin/unpin internal data
1095  unsigned nint = el_pt->ninternal_data();
1096  for (unsigned j = 0; j < nint; j++)
1097  {
1098  Data* data_pt = el_pt->internal_data_pt(j);
1099  unsigned nvalue = data_pt->nvalue();
1100  for (unsigned i = 0; i < nvalue; i++)
1101  {
1102  data_pt->eqn_number(i) = Eqn_number_backup[data_pt][i];
1103  }
1104  }
1105  }
1106 
1107  // Free up storage
1108  Eqn_number_backup.clear();
1109  }
1110 
1111  private:
1112  // oomph-lib objects
1113  // -----------------
1114 
1115  // Pointers to preconditioner (=inexact solver) objects
1116  // -----------------------------------------------------
1117  /// Pointer to the 'preconditioner' for the pressure matrix
1119 
1120  /// Pointer to the 'preconditioner' for the F matrix
1122 
1123  /// flag indicating whether the default F preconditioner is used
1125 
1126  /// flag indicating whether the default P preconditioner is used
1128 
1129  /// Control flag is true if the preconditioner has been setup
1130  /// (used so we can wipe the data when the preconditioner is
1131  /// called again)
1133 
1134  /// Boolean to indicate that presence of elements that are not of
1135  /// type NavierStokesElementWithDiagonalMassMatrices is acceptable
1136  /// (suppresses warning that issued otherwise).
1138 
1139  /// Helper function to assemble the inverse diagonals of the pressure
1140  /// and velocity mass matrices from the elemental contributions defined in
1141  /// NavierStokesEquations<DIM>.
1142  /// If do_both=true, both are computed, otherwise only the velocity
1143  /// mass matrix (the LSC version of the preconditioner only needs
1144  /// that one)
1146  CRDoubleMatrix*& inv_p_mass_pt,
1147  CRDoubleMatrix*& inv_v_mass_pt,
1148  const bool& do_both);
1149 
1150  /// Boolean indicating whether the momentum system preconditioner
1151  /// is a block preconditioner
1153 
1154  /// Set Doc_time to true for outputting results of timings
1155  bool Doc_time;
1156 
1157  /// MatrixVectorProduct operator for Qv^{-1} Bt
1159 
1160  /// MatrixVectorProduct operator for Bt
1162 
1163  /// MatrixVectorProduct operator for F
1165 
1166  /// MatrixVectorProduct operator for E = Fp Qp^{-1} (only for Fp variant)
1168 
1169  /// the pointer to the mesh of block preconditionable Navier
1170  /// Stokes elements.
1172 
1173  /// Flag to indicate if there are multiple element types in the
1174  /// Navier-Stokes mesh.
1176 
1177  /// Boolean to indicate use of LSC (true) or Fp (false) variant
1178  bool Use_LSC;
1179 
1180  /// Use Robin BC elements for Fp preconditoner?
1182 
1183  /// Map to store original eqn numbers of various Data values when
1184  /// assembling pressure advection diffusion matrix
1185  std::map<Data*, std::vector<int>> Eqn_number_backup;
1186 
1187  /// Boolean indicating if we want to pin first pressure
1188  /// dof in Navier Stokes mesh when
1189  /// assembling the pressure advection diffusion system for
1190  /// Fp preconditoner -- needed at zero Reynolds number
1191  /// for non-enclosed flows but seems harmless in any case
1193 
1194  /// Storage for the (non-const!) problem pointer for use in
1195  /// get_pressure_advection_diffusion_matrix().
1197  };
1198 
1199 
1200  /// ////////////////////////////////////////////////////////////////////////////
1201  /// ////////////////////////////////////////////////////////////////////////////
1202  /// ////////////////////////////////////////////////////////////////////////////
1203 
1204 
1205  //============================================================================
1206  /// The exact Navier Stokes preconditioner. This extracts 2x2 blocks
1207  /// (corresponding to the velocity and pressure unknowns) and uses these to
1208  /// build a single preconditioner matrix for testing purposes.
1209  /// Iterative solvers should converge in a single step if this is used.
1210  /// If it doesn't something is wrong in the setup of the block matrices.
1211  //=============================================================================
1212  template<typename MATRIX>
1214  {
1215  public:
1216  /// Constructor - do nothing
1218 
1219 
1220  /// Destructor - do nothing
1222 
1223 
1224  /// Broken copy constructor
1226  delete;
1227 
1228  /// Broken assignment operator
1229  /*void operator=(const NavierStokesExactPreconditioner&) = delete;*/
1230 
1231  /// Setup the preconditioner
1232  void setup();
1233 
1234 
1235  /// Apply preconditioner to r
1237 
1238  protected:
1239  /// Preconditioner matrix
1240  MATRIX P_matrix;
1241  };
1242 
1243 } // namespace oomph
1244 #endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A class that is used to define the functions used to assemble the elemental contributions to the resi...
Block Preconditioner base class. The block structure of the overall problem is determined from the Me...
const Mesh * mesh_pt(const unsigned &i) const
Access to i-th mesh (of the various meshes that contain block preconditionable elements of the same n...
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
double * value_pt(const unsigned &i) const
Return the pointer to the i-the stored value. Typically this is required when direct access to the st...
Definition: nodes.h:324
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:271
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
Information for documentation of results: Directory and file number to enable output in the form RESL...
std::string directory() const
Output directory.
unsigned & number()
Number used (e.g.) for labeling output files.
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
///////////////////////////////////////////////////////////// ///////////////////////////////////////...
virtual ~FpPreconditionerAssemblyHandler()
Empty virtual destructor.
void get_jacobian(GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable" for elem_pt.
FpPreconditionerAssemblyHandler(const unsigned &ndim)
Constructor. Pass spatial dimension.
void get_residuals(GeneralisedElement *const &elem_pt, Vector< double > &residuals)
Return the contribution to the residuals of the element elem_pt.
//////////////////////////////////////////////////////////////// ////////////////////////////////////...
void get_pressure_advection_diffusion_jacobian(CRDoubleMatrix &jacobian)
Get the pressure advection diffusion Jacobian.
FpPressureAdvectionDiffusionProblem(Mesh *navier_stokes_mesh_pt, Problem *orig_problem_pt)
Constructor: Pass Navier-Stokes mesh and pointer to orig problem.
void pin_all_non_pressure_dofs(const bool &set_pressure_bc_for_validation=false)
Pin all non-pressure dofs and (if boolean flag is set to true also all pressure dofs along domain bou...
Problem * Orig_problem_pt
Pointer to orig problem (required so we can re-assign eqn numbers)
void doc_solution(DocInfo &doc_info)
Doc solution (only needed during development – kept alive in case validation is required during code ...
std::map< Data *, std::vector< int > > Eqn_number_backup
Map to store original equation numbers.
void reset_pin_status()
Reset pin status of all values.
void validate(DocInfo &doc_info)
Validate pressure advection diffusion problem and doc exact and computed solution in directory specif...
A Generalised Element class.
Definition: elements.h:73
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number.
Definition: elements.h:704
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
unsigned ninternal_data() const
Return the number of internal data objects.
Definition: elements.h:823
Matrix vector product helper class - primarily a wrapper to Trilinos's Epetra matrix vector product m...
A general mesh class.
Definition: mesh.h:67
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
int face_index_at_boundary(const unsigned &b, const unsigned &e) const
For the e-th finite element on boundary b, return int to indicate the face_index of the face adjacent...
Definition: mesh.h:896
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Definition: mesh.h:878
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt)
Output a given Vector function at f(n_plot) points in each element.
Definition: mesh.cc:2199
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
Definition: mesh.h:840
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
void output(std::ostream &outfile)
Output for all elements.
Definition: mesh.cc:2027
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
void setup()
Broken assignment operator.
NavierStokesExactPreconditioner(const NavierStokesExactPreconditioner &)=delete
Broken copy constructor.
void preconditioner_solve(const Vector< double > &r, Vector< double > &z)
Apply preconditioner to r.
///////////////////////////////////////////////////////////////////////// ///////////////////////////...
void set_p_superlu_preconditioner()
Function to (re-)set pressure matrix preconditioner (inexact solver) to SuperLU.
bool Doc_time
Set Doc_time to true for outputting results of timings.
Preconditioner * F_preconditioner_pt
Pointer to the 'preconditioner' for the F matrix.
void validate(DocInfo &doc_info, Problem *orig_problem_pt)
Validate auxiliary pressure advection diffusion problem in 2D.
void use_fp()
Use Fp version of the preconditioner.
bool Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements
Boolean to indicate that presence of elements that are not of type NavierStokesElementWithDiagonalMas...
Problem * Problem_pt
Storage for the (non-const!) problem pointer for use in get_pressure_advection_diffusion_matrix().
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 Use_robin_for_fp
Use Robin BC elements for Fp preconditoner?
bool Pin_first_pressure_dof_in_press_adv_diff
Boolean indicating if we want to pin first pressure dof in Navier Stokes mesh when assembling the pre...
void set_problem_pt(Problem *problem_pt)
Broken assignment operator.
void disable_accept_non_NavierStokesElementWithDiagonalMassMatrices_elements()
Do not accept presence of elements that are not of type NavierStokesElementWithDiagonalMassMatrices w...
void set_p_preconditioner(Preconditioner *new_p_preconditioner_pt)
Function to set a new pressure matrix preconditioner (inexact solver)
void unpin_first_pressure_dof_in_press_adv_diff()
Set boolean indicating that we do not want to pin first pressure dof in Navier Stokes mesh when assem...
void clean_up_memory()
Helper function to delete preconditioner data.
Mesh * Navier_stokes_mesh_pt
the pointer to the mesh of block preconditionable Navier Stokes elements.
MatrixVectorProduct * F_mat_vec_pt
MatrixVectorProduct operator for F.
MatrixVectorProduct * E_mat_vec_pt
MatrixVectorProduct operator for E = Fp Qp^{-1} (only for Fp variant)
bool Use_LSC
Boolean to indicate use of LSC (true) or Fp (false) variant.
bool Using_default_p_preconditioner
flag indicating whether the default P preconditioner is used
MatrixVectorProduct * QBt_mat_vec_pt
MatrixVectorProduct operator for Qv^{-1} Bt.
bool Using_default_f_preconditioner
flag indicating whether the default F preconditioner is used
bool Allow_multiple_element_type_in_navier_stokes_mesh
Flag to indicate if there are multiple element types in the Navier-Stokes mesh.
Preconditioner * P_preconditioner_pt
Pointer to the 'preconditioner' for the pressure matrix.
void set_f_preconditioner(Preconditioner *new_f_preconditioner_pt)
Function to set a new momentum matrix preconditioner (inexact solver)
void assemble_inv_press_and_veloc_mass_matrix_diagonal(CRDoubleMatrix *&inv_p_mass_pt, CRDoubleMatrix *&inv_v_mass_pt, const bool &do_both)
Helper function to assemble the inverse diagonals of the pressure and velocity mass matrices from the...
NavierStokesSchurComplementPreconditioner(const NavierStokesSchurComplementPreconditioner &)=delete
Broken copy constructor.
void set_navier_stokes_mesh(Mesh *mesh_pt, const bool &allow_multiple_element_type_in_navier_stokes_mesh=false)
Specify the mesh containing the block-preconditionable Navier-Stokes elements. The optional argument ...
void get_pressure_advection_diffusion_matrix(CRDoubleMatrix &fp_matrix)
Get the pressure advection diffusion matrix.
NavierStokesSchurComplementPreconditioner(Problem *problem_pt)
Constructor - sets defaults for control flags.
void enable_robin_for_fp()
Use Robin BC elements for the Fp preconditioner.
bool F_preconditioner_is_block_preconditioner
Boolean indicating whether the momentum system preconditioner is a block preconditioner.
MatrixVectorProduct * Bt_mat_vec_pt
MatrixVectorProduct operator for Bt.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to Vector r.
void disable_robin_for_fp()
Don't use Robin BC elements for the Fp preconditioenr.
void use_lsc()
Use LSC version of the preconditioner.
void set_f_superlu_preconditioner()
Function to (re-)set momentum matrix preconditioner (inexact solver) to SuperLU.
std::map< Data *, std::vector< int > > Eqn_number_backup
Map to store original eqn numbers of various Data values when assembling pressure advection diffusion...
void enable_accept_non_NavierStokesElementWithDiagonalMassMatrices_elements()
Accept presence of elements that are not of type NavierStokesElementWithDiagonalMassMatrices without ...
void pin_first_pressure_dof_in_press_adv_diff()
Set boolean indicating that we want to pin first pressure dof in Navier Stokes mesh when assembling t...
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
An OomphLibError object which should be thrown when an run-time error is encountered....
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
virtual const OomphCommunicator * comm_pt() const
Get function for comm pointer.
virtual void setup()=0
Setup the preconditioner. Pure virtual generic interface function.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
unsigned add_sub_mesh(Mesh *const &mesh_pt)
Add a submesh to the problem and return its number, i, by which it can be accessed via mesh_pt(i).
Definition: problem.h:1330
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Assign all equation numbers for problem: Deals with global data (= data that isn't attached to any el...
Definition: problem.cc:2075
void get_first_and_last_element_for_assembly(Vector< unsigned > &first_el_for_assembly, Vector< unsigned > &last_el_for_assembly) const
Get first and last elements for parallel assembly of non-distributed problem.
Definition: problem.h:1004
AssemblyHandler *& assembly_handler_pt()
Return a pointer to the assembly handler object.
Definition: problem.h:1570
void flush_sub_meshes()
Flush the problem's collection of sub-meshes. Must be followed by call to rebuild_global_mesh().
Definition: problem.h:1339
virtual void get_jacobian(DoubleVector &residuals, DenseDoubleMatrix &jacobian)
Return the fully-assembled Jacobian and residuals for the problem Interface for the case when the Jac...
Definition: problem.cc:3921
void newton_solve()
Use Newton method to solve the problem.
Definition: problem.cc:8976
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
void set_default_first_and_last_element_for_assembly()
Set default first and last elements for parallel assembly of non-distributed problem.
Definition: problem.cc:1699
void set_first_and_last_element_for_assembly(Vector< unsigned > &first_el_for_assembly, Vector< unsigned > &last_el_for_assembly)
Manually set first and last elements for parallel assembly of non-distributed problem.
Definition: problem.h:988
unsigned nsub_mesh() const
Return number of submeshes.
Definition: problem.h:1323
bool distributed() const
If we have MPI return the "problem has been distributed" flag, otherwise it can't be distributed so r...
Definition: problem.h:808
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
Definition: nodes.h:1686
Data *const & variable_position_pt() const
Pointer to variable_position data (const version)
Definition: nodes.h:1765
An interface to allow SuperLU to be used as an (exact) Preconditioner.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
virtual void pin_all_non_pressure_dofs(std::map< Data *, std::vector< int >> &eqn_number_backup)=0
Pin all non-pressure dofs and backup eqn numbers of all Data.
virtual int p_local_eqn(const unsigned &n) const =0
Access function for the local equation number information for the pressure. p_local_eqn[n] = local eq...
virtual int & pinned_fp_pressure_eqn()=0
Global eqn number of pressure dof that's pinned in pressure adv diff problem.
virtual void build_fp_press_adv_diff_robin_bc_element(const unsigned &face_index)=0
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
virtual void delete_pressure_advection_diffusion_robin_elements()=0
Delete the FaceElements that apply the Robin boundary condition to the pressure advection diffusion p...
void get_exact_u(const Vector< double > &x, Vector< double > &u)
Exact solution as a Vector.
double source_function(const Vector< double > &x_vect)
Source function required to make the solution above an exact solution.
void wind_function(const Vector< double > &x, Vector< double > &wind)
Wind.
double Peclet
Peclet number – overwrite with actual Reynolds number.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...