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