26 #ifndef OOMPH_NAVIER_STOKES_PRECONDITIONERS_HEADER
27 #define OOMPH_NAVIER_STOKES_PRECONDITIONERS_HEADER
32 #include <oomph-lib-config.h>
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"
59 namespace PressureAdvectionDiffusionValidation
68 extern void wind_function(
const Vector<double>& x, Vector<double>& wind);
71 extern void get_exact_u(
const Vector<double>& x, Vector<double>& u);
74 extern void get_exact_u(
const Vector<double>& x,
double& u);
107 unsigned n_dof = elem_pt->
ndof();
108 for (
unsigned i = 0;
i < n_dof;
i++)
114 ->fill_in_pressure_advection_diffusion_residuals(residuals);
124 unsigned n_dof = elem_pt->
ndof();
125 for (
unsigned i = 0;
i < n_dof;
i++)
128 for (
unsigned j = 0; j < n_dof; j++)
130 jacobian(
i, j) = 0.0;
135 ->fill_in_pressure_advection_diffusion_jacobian(residuals, jacobian);
153 template<
class ELEMENT>
163 mesh_pt() = navier_stokes_mesh_pt;
185 oomph_info <<
"Eqn numbers after pinning veloc dofs: " << n_dof
196 oomph_info <<
"Eqn numbers in orig problem after re-assignment: "
206 for (
unsigned j = 0; j < nnod; j++)
209 unsigned nval = nod_pt->
nvalue();
210 for (
unsigned i = 0;
i < nval;
i++)
218 for (
unsigned e = 0;
e < nelem;
e++)
228 std::ostringstream error_message;
229 error_message <<
"Navier Stokes mesh must contain only Navier Stokes"
230 <<
"bulk elements\n";
232 OOMPH_CURRENT_FUNCTION,
233 OOMPH_EXCEPTION_LOCATION);
237 unsigned nint = el_pt->ninternal_data();
238 for (
unsigned j = 0; j < nint; j++)
240 Data* data_pt = el_pt->internal_data_pt(j);
241 unsigned nvalue = data_pt->
nvalue();
242 for (
unsigned i = 0;
i < nvalue;
i++)
258 const bool& set_pressure_bc_for_validation =
false)
262 for (
unsigned e = 0;
e < nelem;
e++)
273 std::ostringstream error_message;
274 error_message <<
"Navier Stokes mesh must contain only Navier Stokes"
275 <<
"bulk elements\n";
277 OOMPH_CURRENT_FUNCTION,
278 OOMPH_EXCEPTION_LOCATION);
284 if (el_pt->p_nodal_index_nst() < 0)
286 std::ostringstream error_message;
287 error_message <<
"Cannot use Fp preconditioner with discontinuous\n"
290 OOMPH_CURRENT_FUNCTION,
291 OOMPH_EXCEPTION_LOCATION);
296 unsigned nint = el_pt->ninternal_data();
297 for (
unsigned j = 0; j < nint; j++)
299 Data* data_pt = el_pt->internal_data_pt(j);
302 unsigned nvalue = data_pt->
nvalue();
304 for (
unsigned i = 0;
i < nvalue;
i++)
316 unsigned nnod = el_pt->nnode();
317 for (
unsigned j = 0; j < nnod; j++)
319 Node* nod_pt = el_pt->node_pt(j);
322 unsigned nvalue = nod_pt->
nvalue();
324 for (
unsigned i = 0;
i < nvalue;
i++)
328 if (
int(
i) != el_pt->p_nodal_index_nst())
338 if (el_pt->p_nodal_index_nst() >= 0)
348 if (set_pressure_bc_for_validation)
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]);
367 oomph_info <<
"\n\n==============================================\n\n";
368 oomph_info <<
"Doing validation for pressure adv diff problem\n";
375 bool set_pressure_bc_for_validation =
true;
385 for (
unsigned e = 0;
e < nel;
e++)
388 ->source_fct_for_pressure_adv_diff() =
393 oomph_info <<
"Eqn numbers after pinning veloc dofs: "
401 for (
unsigned b = 0; b < nbound; b++)
407 for (
unsigned e = 0;
e < n_element;
e++)
425 std::ofstream outfile;
426 outfile.open(
"robin_elements.dat");
427 for (
unsigned e = 0;
e < nel;
e++)
430 ->output_pressure_advection_diffusion_robin_elements(outfile);
444 for (
unsigned b = 0; b < nbound; b++)
450 for (
unsigned e = 0;
e < n_element;
e++)
467 oomph_info <<
"Eqn numbers in orig problem after re-assignment: "
471 oomph_info <<
"Done validation for pressure adv diff problem\n";
472 oomph_info <<
"\n\n==============================================\n\n";
480 std::ofstream some_file;
481 std::ostringstream filename;
491 ->interpolated_p_nst(
s);
493 ->interpolated_x(
s, x);
501 for (
unsigned j = 0; j < nnode; j++)
503 if (
mesh_pt()->node_pt(j)->nvalue() == 3)
512 some_file.open(filename.str().c_str());
513 some_file.precision(20);
518 filename << doc_info.
directory() <<
"/fp_exact_soln" << doc_info.
number()
520 some_file.open(filename.str().c_str());
521 some_file.precision(20);
694 std::ostringstream error_msg;
695 error_msg <<
"Problem pointer is null, did you set it yet?";
697 error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
734 const bool& allow_multiple_element_type_in_navier_stokes_mesh =
false)
741 allow_multiple_element_type_in_navier_stokes_mesh;
853 template<
class ELEMENT>
866 for (
unsigned e = 0;
e < nelem;
e++)
897 for (
unsigned i = 0;
i < n_sub_mesh;
i++)
916 backed_up_first_el_for_assembly, backed_up_last_el_for_assembly);
924 int pinned_pressure_eqn = -2;
930 for (
unsigned e = 0;
e < nel;
e++)
938 pinned_pressure_eqn = bulk_elem_pt->
eqn_number(local_eqn);
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,
955 pinned_pressure_eqn = pinned_pressure_eqn_global;
964 for (
unsigned e = 0;
e < nel;
e++)
980 for (
unsigned b = 0; b < nbound; b++)
986 for (
unsigned e = 0;
e < n_element;
e++)
1011 for (
unsigned b = 0; b < nbound; b++)
1017 for (
unsigned e = 0;
e < n_element;
e++)
1033 #ifdef OOMPH_HAS_MPI
1038 backed_up_first_el_for_assembly, backed_up_last_el_for_assembly);
1048 for (
unsigned i = 0;
i < n_sub_mesh;
i++)
1064 for (
unsigned j = 0; j < nnod; j++)
1067 unsigned nval = nod_pt->
nvalue();
1068 for (
unsigned i = 0;
i < nval;
i++)
1075 if (solid_nod_pt != 0)
1078 unsigned nval = solid_posn_data_pt->
nvalue();
1079 for (
unsigned i = 0;
i < nval;
i++)
1089 for (
unsigned e = 0;
e < nelem;
e++)
1096 for (
unsigned j = 0; j < nint; j++)
1099 unsigned nvalue = data_pt->
nvalue();
1100 for (
unsigned i = 0;
i < nvalue;
i++)
1148 const bool& do_both);
1212 template<
typename MATRIX>
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.
A class that represents a collection of data; each Data object may contain many different individual ...
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
void pin(const unsigned &i)
Pin the i-th stored variable.
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...
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...
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
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....
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
///////////////////////////////////////////////////////////// ///////////////////////////////////////...
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.
//////////////////////////////////////////////////////////////// ////////////////////////////////////...
unsigned Ndim
The spatial dimension.
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.
unsigned ndof() const
Return the number of equations/dofs in the element.
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number.
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
unsigned ninternal_data() const
Return the number of internal data objects.
Matrix vector product helper class - primarily a wrapper to Trilinos's Epetra matrix vector product m...
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
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...
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
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.
unsigned nboundary() const
Return number of boundaries.
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
unsigned long nnode() const
Return number of nodes in the mesh.
void output(std::ostream &outfile)
Output for all elements.
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
unsigned long nelement() const
Return number of elements in the mesh.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
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.
~NavierStokesExactPreconditioner()
Destructor - do nothing.
MATRIX P_matrix
Preconditioner matrix.
NavierStokesExactPreconditioner()
Constructor - do nothing.
///////////////////////////////////////////////////////////////////////// ///////////////////////////...
void pin_all_non_pressure_dofs()
Pin all non-pressure dofs.
void setup()
Setup the preconditioner.
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.
virtual ~NavierStokesSchurComplementPreconditioner()
Destructor.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to Vector r.
void disable_doc_time()
Disable documentation of time.
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 enable_doc_time()
Enable documentation of time.
void reset_pin_status()
Reset pin status of all values.
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...
Problem * problem_pt() const
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
double & x(const unsigned &i)
Return the i-th nodal coordinate.
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.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
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).
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...
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.
AssemblyHandler *& assembly_handler_pt()
Return a pointer to the assembly handler object.
void flush_sub_meshes()
Flush the problem's collection of sub-meshes. Must be followed by call to rebuild_global_mesh().
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...
void newton_solve()
Use Newton method to solve the problem.
Mesh *& mesh_pt()
Return a pointer to the global mesh.
void set_default_first_and_last_element_for_assembly()
Set default first and last elements for parallel assembly of non-distributed problem.
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.
unsigned nsub_mesh() const
Return number of submeshes.
bool distributed() const
If we have MPI return the "problem has been distributed" flag, otherwise it can't be distributed so r...
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
Data *const & variable_position_pt() const
Pointer to variable_position data (const version)
An interface to allow SuperLU to be used as an (exact) Preconditioner.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
virtual void pin_all_non_pressure_dofs(std::map< Data *, std::vector< int >> &eqn_number_backup)=0
Pin all non-pressure dofs and backup eqn numbers of all Data.
virtual int p_local_eqn(const unsigned &n) const =0
Access function for the local equation number information for the pressure. p_local_eqn[n] = local eq...
virtual int & pinned_fp_pressure_eqn()=0
Global eqn number of pressure dof that's pinned in pressure adv diff problem.
virtual void build_fp_press_adv_diff_robin_bc_element(const unsigned &face_index)=0
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
virtual void delete_pressure_advection_diffusion_robin_elements()=0
Delete the FaceElements that apply the Robin boundary condition to the pressure advection diffusion p...
void get_exact_u(const Vector< double > &x, Vector< double > &u)
Exact solution as a Vector.
double source_function(const Vector< double > &x_vect)
Source function required to make the solution above an exact solution.
void wind_function(const Vector< double > &x, Vector< double > &wind)
Wind.
unsigned Flag
Flag for solution.
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...