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-2022 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"
46
47
48namespace 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
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 {
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 {
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);
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);
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 {
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.
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();
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
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
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
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
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
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
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
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
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
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
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
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
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
void output(std::ostream &outfile)
Output for all elements.
Definition: mesh.cc:2027
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:1989
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
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:3890
void newton_solve()
Use Newton method to solve the problem.
Definition: problem.cc:8783
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:1613
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
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
AssemblyHandler *& assembly_handler_pt()
Return a pointer to the assembly handler object.
Definition: problem.h:1570
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 int & pinned_fp_pressure_eqn()=0
Global eqn number of pressure dof that's pinned in pressure adv diff problem.
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 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...