29 #ifdef OOMPH_HAS_HYPRE
30 #include "../../generic/hypre_solver.h"
39 #ifdef OOMPH_HAS_HYPRE
44 namespace HypreSubsidiaryPreconditionerHelper
64 bool doc_block_matrices =
false;
85 std::ostringstream error_message;
88 error_message <<
"NavierStokesSchurComplementPreconditioner only works "
89 <<
"with CRDoubleMatrix matrices" << std::endl;
93 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
107 unsigned n_dof_types = 0;
116 oomph_info <<
"Number of dof types: " << n_dof_types << std::endl;
119 if (n_dof_types != 3)
122 std::ostringstream error_message;
125 error_message <<
"Should only be 3 dof types! You have " << n_dof_types
126 <<
" dof types!" << std::endl;
130 OOMPH_CURRENT_FUNCTION,
131 OOMPH_EXCEPTION_LOCATION);
137 throw OomphLibError(
"Can only be used as a subsidiary preconditioner!",
138 OOMPH_CURRENT_FUNCTION,
139 OOMPH_EXCEPTION_LOCATION);
146 dof_to_block_map[n_dof_types - 1] = 1;
155 double block_setup_time = t_block_setup_finish - t_block_setup_start;
158 oomph_info <<
"Time for block_setup(...) [sec]: " << block_setup_time
165 oomph_info <<
"Number of block types: " << n_block_type << std::endl;
189 d_pt->
multiply(*g_pt, *p_matrix_pt);
192 if (doc_block_matrices)
199 n_block_type, n_block_type, 0);
202 for (
unsigned i = 0;
i < n_block_type;
i++)
205 for (
unsigned j = 0; j < n_block_type; j++)
225 oomph_info <<
"Done output of j_" + suffix +
".csv" << std::endl;
228 for (
unsigned i = 0;
i < n_block_type;
i++)
231 for (
unsigned j = 0; j < n_block_type; j++)
234 delete block_matrix_pt(
i, j);
237 block_matrix_pt(
i, j) = 0;
251 oomph_info <<
"Done output of f_matrix" + suffix +
".dat" << std::endl;
257 oomph_info <<
"Done output of d_matrix" + suffix +
".dat" << std::endl;
263 oomph_info <<
"Done output of g_matrix" + suffix +
".dat" << std::endl;
270 oomph_info <<
"Done output of p_matrix" + suffix +
".dat" << std::endl;
301 unsigned n_row_f = f_pt->
nrow();
304 unsigned n_nnz_f = f_pt->
nnz();
309 (n_nnz_f * (
sizeof(int) +
sizeof(
double))));
312 unsigned n_row_g = g_pt->
nrow();
315 unsigned n_nnz_g = g_pt->
nnz();
319 (n_nnz_g * (
sizeof(int) +
sizeof(
double))));
328 #ifdef OOMPH_HAS_HYPRE
393 ->get_total_memory_needed_for_superlu();
399 unsigned n_row_p = p_matrix_pt->
nrow();
402 unsigned n_nnz_p = p_matrix_pt->
nnz();
407 (n_nnz_p * (
sizeof(int) +
sizeof(
double))));
439 ->get_total_memory_needed_for_superlu();
446 unsigned n_row_f = f_pt->
nrow();
449 unsigned n_nnz_f = f_pt->
nnz();
454 (n_nnz_f * (
sizeof(int) +
sizeof(
double))));
516 another_temp_vec.
clear();
528 another_temp_vec.
clear();
539 temp_vec -= another_temp_vec;
556 another_temp_vec.
clear();
562 another_temp_vec += temp_vec;
601 unsigned n_dof_types = 0;
610 oomph_info <<
"Number of dof types: " << n_dof_types << std::endl;
613 if (n_dof_types != 3)
616 std::ostringstream error_message;
619 error_message <<
"Should only be 3 dof types! You have " << n_dof_types
620 <<
" dof types!" << std::endl;
624 OOMPH_CURRENT_FUNCTION,
625 OOMPH_EXCEPTION_LOCATION);
632 throw OomphLibError(
"Currently only used as a subsidiary preconditioner!",
633 OOMPH_CURRENT_FUNCTION,
634 OOMPH_EXCEPTION_LOCATION);
648 double block_setup_time = t_block_setup_finish - t_block_setup_start;
651 oomph_info <<
"Time for block_setup(...) [sec]: " << block_setup_time
658 oomph_info <<
"Number of block types: " << n_block_type << std::endl;
682 (n_nnz * (
sizeof(int) +
sizeof(
double))));
696 for (
unsigned i = 0;
i < n_dof_types;
i++)
736 if (n_block_types != 1)
740 OOMPH_CURRENT_FUNCTION,
741 OOMPH_EXCEPTION_LOCATION);
760 unsigned n_dof = block_r[0].nrow();
778 double t_prec_application_total = 0.0;
803 r_updated, block_z_with_size_of_full_z);
812 t_prec_application_total +=
813 (t_prec_application_end - t_prec_application_start);
823 double normb = r.
norm();
836 resid = beta / normb;
864 oomph_info <<
"GMRES block preconditioner converged immediately. "
865 <<
"Normalised residual norm: " << resid << std::endl;
881 oomph_info <<
"Time for all preconditioner applications [sec]: "
882 << t_prec_application_total
883 <<
"\nTotal time for solve with GMRES block preconditioner "
913 unsigned iter_restart;
916 for (iter_restart = 0; iter_restart < n_dof && iter <=
Max_iter;
917 iter_restart++, iter++)
920 H[iter_restart].resize(iter_restart + 2);
949 r_updated, block_z_with_size_of_full_z);
959 t_prec_application_total +=
960 (t_prec_application_end - t_prec_application_start);
980 r_updated, block_z_with_size_of_full_z);
990 t_prec_application_total +=
991 (t_prec_application_end - t_prec_application_start);
1002 for (
unsigned k = 0; k <= iter_restart; k++)
1005 H[iter_restart][k] = 0.0;
1008 double* vk_pt = v[k].values_pt();
1011 H[iter_restart][k] += w.
dot(v[k]);
1014 for (
unsigned i = 0;
i < n_dof;
i++)
1017 w_pt[
i] -= H[iter_restart][k] * vk_pt[
i];
1024 H[iter_restart][iter_restart + 1] = w.
norm();
1027 v[iter_restart + 1] = w;
1030 v[iter_restart + 1] /= H[iter_restart][iter_restart + 1];
1033 for (
unsigned k = 0; k < iter_restart; k++)
1037 H[iter_restart][k], H[iter_restart][k + 1], cs[k], sn[k]);
1042 H[iter_restart][iter_restart + 1],
1049 H[iter_restart][iter_restart + 1],
1055 s[iter_restart + 1],
1060 beta = std::fabs(
s[iter_restart + 1]);
1063 resid = beta / normb;
1072 oomph_info << iter <<
" " << resid << std::endl;
1098 iter_restart, H,
s, v, block_z_with_size_of_full_z, block_z[0]);
1109 <<
"\nGMRES block preconditioner converged (1). Normalised "
1110 <<
"residual norm: " << resid
1111 <<
"\nNumber of iterations to convergence: " << iter <<
"\n"
1130 <<
"Time for all preconditioner applications [sec]: "
1131 << t_prec_application_total
1132 <<
"\nTotal time for solve with GMRES block preconditioner "
1142 if (iter_restart > 0)
1153 (iter_restart - 1), H,
s, v, block_z_with_size_of_full_z, block_z[0]);
1186 r_updated, block_z_with_size_of_full_z);
1196 t_prec_application_total +=
1197 (t_prec_application_end - t_prec_application_start);
1205 resid = beta / normb;
1214 <<
"\nGMRES block preconditioner converged (2). Normalised "
1215 <<
"residual norm: " << resid
1216 <<
"\nNumber of iterations to convergence: " << iter <<
"\n"
1229 <<
"Time for all preconditioner applications [sec]: "
1230 << t_prec_application_total
1231 <<
"\nTotal time for solve with GMRES block preconditioner "
1240 oomph_info <<
"\nGMRES block preconditioner did not converge to required "
1241 <<
"tolerance! \nReturning with normalised residual norm: "
1242 << resid <<
"\nafter " <<
Max_iter <<
" iterations.\n"
1247 std::ostringstream error_message_stream;
1248 error_message_stream <<
"Solver failed to converge and you requested an "
1249 <<
"error on convergence failures.";
1251 OOMPH_EXCEPTION_LOCATION,
1252 OOMPH_CURRENT_FUNCTION);
void return_block_vector(const unsigned &n, const DoubleVector &b, DoubleVector &v) const
Takes the n-th block ordered vector, b, and copies its entries to the appropriate entries in the natu...
void get_block(const unsigned &i, const unsigned &j, CRDoubleMatrix &output_matrix, const bool &ignore_replacement_block=false) const
Put block (i,j) into output_matrix. This block accounts for any coarsening of dof types and any repla...
void get_block_vector(const unsigned &n, const DoubleVector &v, DoubleVector &b) const
Takes the naturally ordered vector, v and returns the n-th block vector, b. Here n is the block numbe...
unsigned nblock_types() const
Return the number of block types.
CRDoubleMatrix * matrix_pt() const
Access function to matrix_pt. If this is the master then cast the matrix pointer to MATRIX*,...
bool is_subsidiary_block_preconditioner() const
Return true if this preconditioner is a subsidiary preconditioner.
unsigned ndof_types() const
Return the total number of DOF types.
void turn_into_subsidiary_block_preconditioner(BlockPreconditioner< MATRIX > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)
Function to turn this preconditioner into a subsidiary preconditioner that operates within a bigger "...
void get_block_vectors(const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const
Takes the naturally ordered vector and rearranges it into a vector of sub vectors corresponding to th...
void setup_matrix_vector_product(MatrixVectorProduct *matvec_prod_pt, CRDoubleMatrix *block_pt, const Vector< unsigned > &block_col_indices)
Setup a matrix vector product. matvec_prod_pt is a pointer to the MatrixVectorProduct,...
virtual void block_setup()
Determine the size of the matrix blocks and setup the lookup schemes relating the global degrees of f...
A class for compressed row matrices. This is a distributable object.
void sparse_indexed_output_with_offset(std::string filename)
Indexed output function to print a matrix to a file as i,j,a(i,j) for a(i,j)!=0 only....
void multiply(const DoubleVector &x, DoubleVector &soln) const
Multiply the matrix by the vector x: soln=Ax.
unsigned long nnz() const
Return the number of nonzero entries (the local nnz)
unsigned long nrow() const
Return the number of rows of the matrix.
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
virtual void residual(const DoubleVector &x, const DoubleVector &b, DoubleVector &residual_)
Find the residual, i.e. r=b-Ax the residual.
A vector in the mathematical sense, initially developed for linear algebra type applications....
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
double * values_pt()
access function to the underlying values
double norm() const
compute the 2 norm of this vector
double dot(const DoubleVector &vec) const
compute the dot product of this vector with the vector vec.
void clear()
wipes the DoubleVector
double Memory_usage_in_bytes
Storage for the memory usage of the solver if the flag above is set to true (in bytes)
void apply_plane_rotation(double &dx, double &dy, double &cs, double &sn)
Helper function: Apply plane rotation. This is done using the update:
virtual void clean_up_memory()
Clean up the memory (empty for now...)
SpaceTimeNavierStokesSubsidiaryPreconditioner * Navier_stokes_subsidiary_preconditioner_pt
Pointer to the preconditioner for the block matrix.
bool Compute_memory_statistics
Flag to indicate whether or not to record the memory statistics this preconditioner.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to r.
bool Preconditioner_has_been_setup
Control flag is true if the preconditioner has been setup (used so we can wipe the data when the prec...
unsigned Iterations
Number of iterations taken.
bool Preconditioner_LHS
boolean indicating use of left hand preconditioning (if true) or right hand preconditioning (if false...
CRDoubleMatrix * Matrix_pt
Pointer to matrix.
void setup()
Setup the preconditioner.
void update(const unsigned &k, const Vector< Vector< double >> &H, const Vector< double > &s, const Vector< DoubleVector > &v, const DoubleVector &block_x_with_size_of_full_x, DoubleVector &x)
Helper function to update the result vector using the result, x=x_0+V_m*y.
void generate_plane_rotation(double &dx, double &dy, double &cs, double &sn)
Helper function: Generate a plane rotation. This is done by finding the values of (i....
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
double Tolerance
Convergence tolerance.
bool Throw_error_after_max_iter
Should we throw an error instead of just returning when we hit the max iterations?
unsigned Max_iter
Maximum number of iterations.
double Solution_time
linear solver solution time
std::ofstream Output_file_stream
Output file stream for convergence history.
bool Doc_convergence_history
Flag indicating if the convergence history is to be documented.
bool built() const
if the communicator_pt is null then the distribution is not setup then false is returned,...
bool Doc_time
Boolean flag that indicates whether the time taken.
Matrix vector product helper class - primarily a wrapper to Trilinos's Epetra matrix vector product m...
void multiply_transpose(const DoubleVector &x, DoubleVector &y) const
Apply the transpose of the operator to the vector x and return the result in the vector y.
void multiply(const DoubleVector &x, DoubleVector &y) const
Apply the operator to the vector x and return the result in the vector y.
std::ostream *& stream_pt()
Access function for the stream pointer.
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...
void setup(DoubleMatrixBase *matrix_pt)
Setup the preconditioner: store the matrix pointer and the communicator pointer then call preconditio...
std::ostream * Stream_pt
Pointer to the output stream – defaults to std::cout.
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
bool Silent_preconditioner_setup
Boolean to indicate whether or not the build should be done silently.
General purpose block triangular preconditioner. By default this is Upper triangular....
Preconditioner * P_preconditioner_pt
Pointer to the 'preconditioner' for the pressure matrix.
double Memory_usage_in_bytes
Storage for the memory usage of the solver if the flag above is set to true (in bytes)
void setup()
Setup the preconditioner.
bool Preconditioner_has_been_setup
Control flag is true if the preconditioner has been setup (used so we can wipe the data when the prec...
bool Compute_memory_statistics
Flag to indicate whether or not to record the memory statistics this preconditioner.
MatrixVectorProduct * G_mat_vec_pt
MatrixVectorProduct operator for G.
MatrixVectorProduct * F_mat_vec_pt
MatrixVectorProduct operator for F.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to r.
bool Using_default_f_preconditioner
Flag indicating whether the default F preconditioner is used.
virtual void clean_up_memory()
Clean up the memory.
void enable_doc_memory_usage()
Document the memory usage.
bool Using_default_p_preconditioner
Flag indicating whether the default P preconditioner is used.
Preconditioner * F_preconditioner_pt
Pointer to the 'preconditioner' for the F matrix.
An interface to allow SuperLU to be used as an (exact) Preconditioner.
void initialise(const _Tp &__value)
Iterate over all values and set to the desired value.
void concatenate(const DenseMatrix< CRDoubleMatrix * > &matrix_pt, CRDoubleMatrix &result_matrix)
Concatenate CRDoubleMatrix matrices. The in matrices are concatenated such that the block structure o...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
Preconditioner * set_hypre_preconditioner()
Assign the Hypre preconditioner pointer.
void set_defaults_for_2D_poisson_problem(HyprePreconditioner *hypre_preconditioner_pt)
Set default parameters for use as preconditioner in 2D Poisson-type problem.
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Nullstream oomph_nullstream
///////////////////////////////////////////////////////////////////////// ///////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...