40 //==========================================================================
41 namespace HypreSubsidiaryPreconditionerHelper
44 Preconditioner* set_hypre_preconditioner()
46 return new HyprePreconditioner;
48 } // End of HypreSubsidiaryPreconditionerHelper
58 bool doc_block_matrices =
false;
76 if (
dynamic_cast<CRDoubleMatrix*
>(this->
matrix_pt()) == 0)
79 std::ostringstream error_message;
82 error_message <<
"NavierStokesSchurComplementPreconditioner only works "
83 <<
"with CRDoubleMatrix matrices" << std::endl;
87 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
101 unsigned n_dof_types = 0;
110 oomph_info <<
"Number of dof types: " << n_dof_types << std::endl;
113 if (n_dof_types != 3)
116 std::ostringstream error_message;
119 error_message <<
"Should only be 3 dof types! You have " << n_dof_types
120 <<
" dof types!" << std::endl;
123 throw OomphLibError(error_message.str(),
124 OOMPH_CURRENT_FUNCTION,
125 OOMPH_EXCEPTION_LOCATION);
131 throw OomphLibError(
"Can only be used as a subsidiary preconditioner!",
132 OOMPH_CURRENT_FUNCTION,
133 OOMPH_EXCEPTION_LOCATION);
137 Vector<unsigned> dof_to_block_map(n_dof_types);
140 dof_to_block_map[n_dof_types - 1] = 1;
149 double block_setup_time = t_block_setup_finish - t_block_setup_start;
152 oomph_info <<
"Time for block_setup(...) [sec]: " << block_setup_time
159 oomph_info <<
"Number of block types: " << n_block_type << std::endl;
162 CRDoubleMatrix* f_pt =
new CRDoubleMatrix;
165 CRDoubleMatrix* d_pt =
new CRDoubleMatrix;
168 CRDoubleMatrix* g_pt =
new CRDoubleMatrix;
171 CRDoubleMatrix* p_matrix_pt =
new CRDoubleMatrix;
183 d_pt->multiply(*g_pt, *p_matrix_pt);
186 if (doc_block_matrices)
192 DenseMatrix<CRDoubleMatrix*> block_matrix_pt(
193 n_block_type, n_block_type, 0);
196 for (
unsigned i = 0;
i < n_block_type;
i++)
199 for (
unsigned j = 0; j < n_block_type; j++)
202 block_matrix_pt(
i, j) =
new CRDoubleMatrix;
210 CRDoubleMatrix* j_pt =
new CRDoubleMatrix;
216 j_pt->sparse_indexed_output_with_offset(
"j_" + suffix +
".csv");
219 oomph_info <<
"Done output of j_" + suffix +
".csv" << std::endl;
222 for (
unsigned i = 0;
i < n_block_type;
i++)
225 for (
unsigned j = 0; j < n_block_type; j++)
228 delete block_matrix_pt(
i, j);
231 block_matrix_pt(
i, j) = 0;
242 f_pt->sparse_indexed_output_with_offset(
"f_matrix" + suffix +
".dat");
245 oomph_info <<
"Done output of f_matrix" + suffix +
".dat" << std::endl;
248 d_pt->sparse_indexed_output_with_offset(
"d_matrix" + suffix +
".dat");
251 oomph_info <<
"Done output of d_matrix" + suffix +
".dat" << std::endl;
254 g_pt->sparse_indexed_output_with_offset(
"g_matrix" + suffix +
".dat");
257 oomph_info <<
"Done output of g_matrix" + suffix +
".dat" << std::endl;
260 p_matrix_pt->sparse_indexed_output_with_offset(
"p_matrix" + suffix +
264 oomph_info <<
"Done output of p_matrix" + suffix +
".dat" << std::endl;
295 unsigned n_row_f = f_pt->nrow();
298 unsigned n_nnz_f = f_pt->nnz();
303 (n_nnz_f * (
sizeof(int) +
sizeof(
double))));
306 unsigned n_row_g = g_pt->nrow();
309 unsigned n_nnz_g = g_pt->nnz();
313 (n_nnz_g * (
sizeof(int) +
sizeof(
double))));
337 ->get_total_memory_needed_for_superlu();
365 ->get_total_memory_needed_for_superlu();
397 const DoubleVector& r, DoubleVector& z)
400 if (!z.distribution_pt()->built())
403 z.build(r.distribution_pt(), 0.0);
410 DoubleVector temp_vec;
413 DoubleVector another_temp_vec;
430 another_temp_vec.clear();
442 another_temp_vec.clear();
450 temp_vec.build(another_temp_vec.distribution_pt(), 0.0);
453 temp_vec -= another_temp_vec;
470 another_temp_vec.clear();
476 another_temp_vec += temp_vec;
515 unsigned n_dof_types = 0;
524 oomph_info <<
"Number of dof types: " << n_dof_types << std::endl;
527 if (n_dof_types != 3)
530 std::ostringstream error_message;
533 error_message <<
"Should only be 3 dof types! You have " << n_dof_types
534 <<
" dof types!" << std::endl;
537 throw OomphLibError(error_message.str(),
538 OOMPH_CURRENT_FUNCTION,
539 OOMPH_EXCEPTION_LOCATION);
546 throw OomphLibError(
"Currently only used as a subsidiary preconditioner!",
547 OOMPH_CURRENT_FUNCTION,
548 OOMPH_EXCEPTION_LOCATION);
553 Vector<unsigned> dof_to_block_map(n_dof_types, 0);
562 double block_setup_time = t_block_setup_finish - t_block_setup_start;
565 oomph_info <<
"Time for block_setup(...) [sec]: " << block_setup_time
572 oomph_info <<
"Number of block types: " << n_block_type << std::endl;
582 new SpaceTimeNavierStokesSubsidiaryPreconditioner;
596 (n_nnz * (
sizeof(int) +
sizeof(
double))));
606 Vector<unsigned> dof_map(n_dof_types);
610 for (
unsigned i = 0;
i < n_dof_types;
i++)
644 DoubleVector& solution)
650 if (n_block_types != 1)
653 throw OomphLibError(
"Can only deal with one dof type!",
654 OOMPH_CURRENT_FUNCTION,
655 OOMPH_EXCEPTION_LOCATION);
659 Vector<DoubleVector> block_z(n_block_types);
662 Vector<DoubleVector> block_r(n_block_types);
671 block_z[0].initialise(0.0);
674 unsigned n_dof = block_r[0].nrow();
686 Vector<double>
s(n_dof + 1, 0);
687 Vector<double> cs(n_dof + 1);
688 Vector<double> sn(n_dof + 1);
692 double t_prec_application_total = 0.0;
706 DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt(), 0.0);
707 DoubleVector r_updated(rhs.distribution_pt(), 0.0);
717 r_updated, block_z_with_size_of_full_z);
726 t_prec_application_total +=
727 (t_prec_application_end - t_prec_application_start);
737 double normb = r.norm();
750 resid = beta / normb;
778 oomph_info <<
"GMRES block preconditioner converged immediately. "
779 <<
"Normalised residual norm: " << resid << std::endl;
795 oomph_info <<
"Time for all preconditioner applications [sec]: "
796 << t_prec_application_total
797 <<
"\nTotal time for solve with GMRES block preconditioner "
806 Vector<DoubleVector> v(n_dof + 1);
811 Vector<Vector<double>> H(n_dof + 1);
827 unsigned iter_restart;
830 for (iter_restart = 0; iter_restart < n_dof && iter <=
Max_iter;
831 iter_restart++, iter++)
834 H[iter_restart].resize(iter_restart + 2);
853 DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt(),
855 DoubleVector r_updated(rhs.distribution_pt(), 0.0);
863 r_updated, block_z_with_size_of_full_z);
873 t_prec_application_total +=
874 (t_prec_application_end - t_prec_application_start);
885 DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt());
886 DoubleVector r_updated(rhs.distribution_pt());
894 r_updated, block_z_with_size_of_full_z);
904 t_prec_application_total +=
905 (t_prec_application_end - t_prec_application_start);
913 double* w_pt = w.values_pt();
916 for (
unsigned k = 0; k <= iter_restart; k++)
919 H[iter_restart][k] = 0.0;
922 double* vk_pt = v[k].values_pt();
925 H[iter_restart][k] += w.dot(v[k]);
928 for (
unsigned i = 0;
i < n_dof;
i++)
931 w_pt[
i] -= H[iter_restart][k] * vk_pt[
i];
938 H[iter_restart][iter_restart + 1] = w.norm();
941 v[iter_restart + 1] = w;
944 v[iter_restart + 1] /= H[iter_restart][iter_restart + 1];
947 for (
unsigned k = 0; k < iter_restart; k++)
951 H[iter_restart][k], H[iter_restart][k + 1], cs[k], sn[k]);
956 H[iter_restart][iter_restart + 1],
963 H[iter_restart][iter_restart + 1],
974 beta = std::fabs(
s[iter_restart + 1]);
977 resid = beta / normb;
986 oomph_info << iter <<
" " << resid << std::endl;
1002 DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt(), 0.0);
1012 iter_restart, H,
s, v, block_z_with_size_of_full_z, block_z[0]);
1023 <<
"\nGMRES block preconditioner converged (1). Normalised "
1024 <<
"residual norm: " << resid
1025 <<
"\nNumber of iterations to convergence: " << iter <<
"\n"
1044 <<
"Time for all preconditioner applications [sec]: "
1045 << t_prec_application_total
1046 <<
"\nTotal time for solve with GMRES block preconditioner "
1056 if (iter_restart > 0)
1059 DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt(), 0.0);
1067 (iter_restart - 1), H,
s, v, block_z_with_size_of_full_z, block_z[0]);
1091 DoubleVector block_z_with_size_of_full_z(rhs.distribution_pt());
1092 DoubleVector r_updated(rhs.distribution_pt());
1100 r_updated, block_z_with_size_of_full_z);
1110 t_prec_application_total +=
1111 (t_prec_application_end - t_prec_application_start);
1119 resid = beta / normb;
1128 <<
"\nGMRES block preconditioner converged (2). Normalised "
1129 <<
"residual norm: " << resid
1130 <<
"\nNumber of iterations to convergence: " << iter <<
"\n"
1143 <<
"Time for all preconditioner applications [sec]: "
1144 << t_prec_application_total
1145 <<
"\nTotal time for solve with GMRES block preconditioner "
1154 oomph_info <<
"\nGMRES block preconditioner did not converge to required "
1155 <<
"tolerance! \nReturning with normalised residual norm: "
1156 << resid <<
"\nafter " <<
Max_iter <<
" iterations.\n"
1161 std::ostringstream error_message_stream;
1162 error_message_stream <<
"Solver failed to converge and you requested an "
1163 <<
"error on convergence failures.";
1164 throw OomphLibError(error_message_stream.str(),
1165 OOMPH_EXCEPTION_LOCATION,
1166 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...
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.
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 Doc_time
Boolean flag that indicates whether the time taken.
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.
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.
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.
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.
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...