// LIC// ==================================================================== // LIC// This file forms part of oomph-lib, the object-oriented, // LIC// multi-physics finite-element library, available // LIC// at http://www.oomph-lib.org. // LIC// // LIC// Copyright (C) 2006-2026 Matthias Heil and Andrew Hazel // LIC// // LIC// This library is free software; you can redistribute it and/or // LIC// modify it under the terms of the GNU Lesser General Public // LIC// License as published by the Free Software Foundation; either // LIC// version 2.1 of the License, or (at your option) any later version. // LIC// // LIC// This library is distributed in the hope that it will be useful, // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU // LIC// Lesser General Public License for more details. // LIC// // LIC// You should have received a copy of the GNU Lesser General Public // LIC// License along with this library; if not, write to the Free Software // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA // LIC// 02110-1301 USA. // LIC// // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk. // LIC// // LIC//==================================================================== #include "pseudo_elastic_preconditioner.h" namespace oomph { //============================================================================= /// Functions to create instances of optimal subsidiary operators for /// the PseudoElasticPreconditioner //============================================================================= namespace Pseudo_Elastic_Preconditioner_Subsidiary_Operator_Helper { #ifdef OOMPH_HAS_HYPRE /// AMG w/ GS smoothing for the augmented elastic subsidiary linear /// systems Preconditioner* get_elastic_preconditioner_hypre() { HyprePreconditioner* hypre_preconditioner_pt = new HyprePreconditioner("Hypre for diagonal blocks in pseudo-solid"); hypre_preconditioner_pt->set_amg_iterations(2); hypre_preconditioner_pt->amg_using_simple_smoothing(); if (MPI_Helpers::communicator_pt()->nproc() > 1) { // Jacobi in parallel hypre_preconditioner_pt->amg_simple_smoother() = 0; } else { // Gauss Seidel in serial (was 3 actually...) hypre_preconditioner_pt->amg_simple_smoother() = 1; } hypre_preconditioner_pt->hypre_method() = HyprePreconditioner::BoomerAMG; hypre_preconditioner_pt->amg_damping() = 1.0; hypre_preconditioner_pt->amg_coarsening() = 6; return hypre_preconditioner_pt; } /// AMG w/ GS smoothing for the augmented elastic subsidiary linear /// systems -- calls Hypre version to stay consistent with previous default Preconditioner* get_elastic_preconditioner() { return get_elastic_preconditioner_hypre(); } #endif #ifdef OOMPH_HAS_TRILINOS /// TrilinosML smoothing for the augmented elastic /// subsidiary linear systems Preconditioner* get_elastic_preconditioner_trilinos_ml() { TrilinosMLPreconditioner* trilinos_prec_pt = new TrilinosMLPreconditioner; return trilinos_prec_pt; } /// CG with diagonal preconditioner for the lagrange multiplier /// subsidiary linear systems. Preconditioner* get_lagrange_multiplier_preconditioner() { InnerIterationPreconditioner* prec_pt = new InnerIterationPreconditioner; // Note: This makes CG a proper "inner iteration" for // which GMRES (may) no longer converge. We should really // use FGMRES or GMRESR for this. However, here the solver // is so good that it'll converge very quickly anyway // so there isn't much to be gained by limiting the number // of iterations... // prec_pt->max_iter() = 4; prec_pt->solver_pt()->solver_type() = TrilinosAztecOOSolver::CG; prec_pt->solver_pt()->disable_doc_time(); return prec_pt; } #endif } // namespace Pseudo_Elastic_Preconditioner_Subsidiary_Operator_Helper /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// //============================================================================= // Setup method for the PseudoElasticPreconditioner. //============================================================================= void PseudoElasticPreconditioner::setup() { // clean this->clean_up_memory(); #ifdef PARANOID // paranoid check that meshes have been set if (Elastic_mesh_pt == 0) { std::ostringstream error_message; error_message << "The elastic mesh must be set."; throw OomphLibError( error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION); } if (Lagrange_multiplier_mesh_pt == 0) { std::ostringstream error_message; error_message << "The Lagrange multiplier mesh must be set."; throw OomphLibError( error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION); } #endif // set the meshes this->set_nmesh(2); this->set_mesh(0, Elastic_mesh_pt); this->set_mesh(1, Lagrange_multiplier_mesh_pt); // Figure out number of dof types unsigned n_solid_dof_types = 0; unsigned n_dof_types = 0; if (this->is_master_block_preconditioner()) { // get the number of solid dof types from the first element n_solid_dof_types = this->ndof_types_in_mesh(0); // get the total number of dof types n_dof_types = n_solid_dof_types + this->ndof_types_in_mesh(1); } else { n_dof_types = this->ndof_types(); n_solid_dof_types = n_dof_types - (n_dof_types / 3); } #ifdef PARANOID if (n_dof_types % 3 != 0) { std::ostringstream error_message; error_message << "This preconditioner requires DIM*3 types of DOF"; throw OomphLibError( error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION); } #endif // determine the dimension Dim = n_dof_types / 3; #ifdef PARANOID // Recast Jacobian matrix to CRDoubleMatrix CRDoubleMatrix* cr_matrix_pt = dynamic_cast(matrix_pt()); if (cr_matrix_pt == 0) { std::ostringstream error_message; error_message << "FSIPreconditioner only works with" << " CRDoubleMatrix matrices" << std::endl; throw OomphLibError( error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION); } #endif // Setup the dof_list scheme for block_setup. // The natural DOF types order is (in 3D): // 0 1 2 3 4 5 6 7 8 <- natural DOF type index // x y z xc yc zc lx ly lz <- DOf type // // We need to group the directional displacements together: // 0 1 2 3 4 5 6 7 8 <- block index // x xc y yc xz zc lx ly lz <- DOF type // // The mapping required is: // 0 2 4 1 3 5 6 7 8 // // In general we want: // // dof_to_block_map[DOF type] = block type Vector dof_list_for_block_setup(n_dof_types); for (unsigned dim_i = 0; dim_i < Dim; dim_i++) { // bulk solid dof types dof_list_for_block_setup[dim_i] = 2 * dim_i; // constrained solid dof types dof_list_for_block_setup[dim_i + Dim] = 2 * dim_i + 1; // lagr dof types dof_list_for_block_setup[dim_i + 2 * Dim] = 2 * Dim + dim_i; } // Setup the block ordering. We have the following: // Block types [0, 2*Dim) are the solid blocks. // Block types [2*Dim, 3*Dim) are the Lagrange multiplier blocks. // // For d = 0, 1, 2, // Bulk solid doftypes: 2*d // Constrained solid doftypes: 2*d+1 // Lgr doftyoes: 2*Dim + d this->block_setup(dof_list_for_block_setup); // Create dof_list for subsidiary preconditioner. This will be used later // in turn_into_subsidiary_block_preconditioner(...) Vector dof_list_for_subsidiary_prec(n_solid_dof_types); // The dof types are currently in the order (in 3D): // 0 1 2 3 4 5 6 7 8 <- DOF index // x y z xc yc zc lx ly lz <- DOF type // // We need to group the directional displacements together: // x xc y yc xz zc // // The mapping required is: // 0 3 1 4 2 5 // // In general we want: // dof_list[subsidiary DOF type] = master DOF type for (unsigned dim_i = 0; dim_i < Dim; dim_i++) { // bulk solid dof dof_list_for_subsidiary_prec[2 * dim_i] = dim_i; // constrained solid dof dof_list_for_subsidiary_prec[2 * dim_i + 1] = dim_i + Dim; } // Get the solid blocks DenseMatrix solid_matrix_pt( n_solid_dof_types, n_solid_dof_types, 0); for (unsigned row_i = 0; row_i < n_solid_dof_types; row_i++) { for (unsigned col_i = 0; col_i < n_solid_dof_types; col_i++) { solid_matrix_pt(row_i, col_i) = new CRDoubleMatrix; this->get_block(row_i, col_i, *solid_matrix_pt(row_i, col_i)); } } // compute the scaling (if required) if (Use_inf_norm_of_s_scaling) { Scaling = CRDoubleMatrixHelpers::inf_norm(solid_matrix_pt); } else { Scaling = 1.0; } // Add the scaled identify matrix to the constrained solid blocks. for (unsigned d = 0; d < Dim; d++) { // Where is the constrained block located? unsigned block_i = 2 * d + 1; // Data from the constrained block. double* s_values = solid_matrix_pt(block_i, block_i)->value(); int* s_column_index = solid_matrix_pt(block_i, block_i)->column_index(); int* s_row_start = solid_matrix_pt(block_i, block_i)->row_start(); int s_nrow_local = solid_matrix_pt(block_i, block_i)->nrow_local(); int s_first_row = solid_matrix_pt(block_i, block_i)->first_row(); // Loop through the diagonal entries and add the scaling. for (int i = 0; i < s_nrow_local; i++) { bool found = false; for (int j = s_row_start[i]; j < s_row_start[i + 1] && !found; j++) { if (s_column_index[j] == i + s_first_row) { s_values[j] += Scaling; found = true; } // if } // for row start // Check if the diagonal entry is found. if (!found) { std::ostringstream error_message; error_message << "The diagonal entry for the constained block(" << block_i << "," << block_i << ")\n" << "on local row " << i << " does not exist." << std::endl; throw OomphLibError(error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION); } } // for nrow_local } // for Dim // setup the solid subsidiary preconditioner // //////////////////////////////// this preconditioner uses the full S // matrix if (E_preconditioner_type == Exact_block_preconditioner) { ExactBlockPreconditioner* s_prec_pt = new ExactBlockPreconditioner; // Add mesh (not actually used since this only acts as a subsidiary // preconditioner... // s_prec_pt->add_mesh(Elastic_mesh_pt); Vector> doftype_to_doftype_map(n_solid_dof_types, Vector(1, 0)); for (unsigned i = 0; i < n_solid_dof_types; i++) { doftype_to_doftype_map[i][0] = i; } s_prec_pt->turn_into_subsidiary_block_preconditioner( this, dof_list_for_subsidiary_prec, doftype_to_doftype_map); if (Elastic_subsidiary_preconditioner_function_pt != 0) { s_prec_pt->set_subsidiary_preconditioner_function( Elastic_subsidiary_preconditioner_function_pt); } // Set the replacement blocks. for (unsigned d = 0; d < Dim; d++) { // The dof-block located in the block-ordering. unsigned block_i = 2 * d + 1; // The dof-block located in the dof-ordering. unsigned dof_block_i = Dim + d; this->set_replacement_dof_block( dof_block_i, dof_block_i, solid_matrix_pt(block_i, block_i)); } s_prec_pt->Preconditioner::setup(matrix_pt()); Elastic_preconditioner_pt = s_prec_pt; } // otherwise it is a block based preconditioner else { GeneralPurposeBlockPreconditioner* s_prec_pt = 0; // set the block preconditioning method switch (E_preconditioner_type) { case Block_diagonal_preconditioner: { s_prec_pt = new BlockDiagonalPreconditioner; } break; case Block_upper_triangular_preconditioner: { BlockTriangularPreconditioner* block_triangular_prec_pt = new BlockTriangularPreconditioner; block_triangular_prec_pt->upper_triangular(); s_prec_pt = block_triangular_prec_pt; } break; case Block_lower_triangular_preconditioner: { BlockTriangularPreconditioner* block_triangular_prec_pt = new BlockTriangularPreconditioner; block_triangular_prec_pt->lower_triangular(); s_prec_pt = block_triangular_prec_pt; } break; default: { std::ostringstream error_msg; error_msg << "There is no such block based preconditioner.\n" << "Candidates are:\n" << "PseudoElasticPreconditioner::Block_diagonal_preconditioner\n" << "PseudoElasticPreconditioner::Block_upper_triangular_" "preconditioner\n" << "PseudoElasticPreconditioner::Block_lower_triangular_" "preconditioner\n"; throw OomphLibError( error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION); } break; } // Add mesh (not actually used since this only acts as a subsidiary // preconditioner... // s_prec_pt->add_mesh(Elastic_mesh_pt); // The block to block map Vector> doftype_to_doftype_map(Dim, Vector(2, 0)); Vector s_prec_dof_to_block_map(Dim, 0); unsigned tmp_index = 0; for (unsigned d = 0; d < Dim; d++) { s_prec_dof_to_block_map[d] = d; doftype_to_doftype_map[d][0] = tmp_index++; doftype_to_doftype_map[d][1] = tmp_index++; } s_prec_pt->turn_into_subsidiary_block_preconditioner( this, dof_list_for_subsidiary_prec, doftype_to_doftype_map); if (Elastic_subsidiary_preconditioner_function_pt != 0) { s_prec_pt->set_subsidiary_preconditioner_function( Elastic_subsidiary_preconditioner_function_pt); } // Set the replacement blocks. for (unsigned d = 0; d < Dim; d++) { // The dof-block located in the block-ordering. unsigned block_i = 2 * d + 1; // Then dof-block located in the dof-ordering. unsigned dof_block_i = Dim + d; this->set_replacement_dof_block( dof_block_i, dof_block_i, solid_matrix_pt(block_i, block_i)); } s_prec_pt->set_dof_to_block_map(s_prec_dof_to_block_map); s_prec_pt->Preconditioner::setup(matrix_pt()); Elastic_preconditioner_pt = s_prec_pt; } // No longer require the solid blocks for (unsigned row_i = 0; row_i < n_solid_dof_types; row_i++) { for (unsigned col_i = 0; col_i < n_solid_dof_types; col_i++) { delete solid_matrix_pt(row_i, col_i); solid_matrix_pt(row_i, col_i) = 0; } } // next setup the lagrange multiplier preconditioners Lagrange_multiplier_preconditioner_pt.resize(Dim); for (unsigned d = 0; d < Dim; d++) { CRDoubleMatrix* b_pt = new CRDoubleMatrix; this->get_block(2 * Dim + d, 2 * d + 1, *b_pt); // if a non default preconditioner is specified create // the preconditioners if (Lagrange_multiplier_subsidiary_preconditioner_function_pt != 0) { Lagrange_multiplier_preconditioner_pt[d] = (*Lagrange_multiplier_subsidiary_preconditioner_function_pt)(); } // else use default superlu preconditioner else { Lagrange_multiplier_preconditioner_pt[d] = ExactPreconditionerFactory::create_exact_preconditioner(); } // and setup Lagrange_multiplier_preconditioner_pt[d]->setup(b_pt); delete b_pt; b_pt = 0; } } //============================================================================= /// Apply the elastic subsidiary preconditioner. //============================================================================= void PseudoElasticPreconditioner::elastic_preconditioner_solve( const DoubleVector& r, DoubleVector& z) { // apply the solid preconditioner Elastic_preconditioner_pt->preconditioner_solve(r, z); } //============================================================================= /// Apply the lagrange multiplier subsidiary preconditioner. //============================================================================= void PseudoElasticPreconditioner::lagrange_multiplier_preconditioner_solve( const DoubleVector& r, DoubleVector& z) { // apply the lagrange multiplier preconditioner for (unsigned d = 0; d < Dim; d++) { DoubleVector x; this->get_block_vector(Dim * 2 + d, r, x); DoubleVector y; Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(x, y); Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(y, x); unsigned nrow_local = x.nrow_local(); double* x_pt = x.values_pt(); for (unsigned i = 0; i < nrow_local; i++) { x_pt[i] = x_pt[i] * Scaling; } this->return_block_vector(Dim * 2 + d, x, z); } } //============================================================================= /// Clears the memory. //============================================================================= void PseudoElasticPreconditioner::clean_up_memory() { // clean the block preconditioner base class memory this->clear_block_preconditioner_base(); // delete the solid preconditioner delete Elastic_preconditioner_pt; Elastic_preconditioner_pt = 0; // delete the lagrange multiplier preconditioner pt unsigned sz = Lagrange_multiplier_preconditioner_pt.size(); for (unsigned i = 0; i < sz; i++) { delete Lagrange_multiplier_preconditioner_pt[i]; Lagrange_multiplier_preconditioner_pt[i] = 0; } } /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// //============================================================================= // Setup method for the PseudoElasticPreconditioner. //============================================================================= void PseudoElasticPreconditionerOld::setup() { // clean this->clean_up_memory(); #ifdef PARANOID // paranoid check that meshes have been set if (Elastic_mesh_pt == 0) { std::ostringstream error_message; error_message << "The elastic mesh must be set."; throw OomphLibError( error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION); } if (Lagrange_multiplier_mesh_pt == 0) { std::ostringstream error_message; error_message << "The Lagrange multiplier mesh must be set."; throw OomphLibError( error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION); } #endif // set the mesh unsigned n_solid_dof_types = 0; unsigned n_dof_types = 0; this->set_mesh(0, Elastic_mesh_pt); this->set_mesh(1, Lagrange_multiplier_mesh_pt); if (this->is_master_block_preconditioner()) { // get the number of solid dof types from the first element n_solid_dof_types = this->ndof_types_in_mesh(0); // get the total number of dof types n_dof_types = n_solid_dof_types + this->ndof_types_in_mesh(1); } else { n_dof_types = this->ndof_types(); n_solid_dof_types = n_dof_types - (n_dof_types / 3); } #ifdef PARANOID if (n_dof_types % 3 != 0) { std::ostringstream error_message; error_message << "This preconditioner requires DIM*3 types of DOF"; throw OomphLibError( error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION); } #endif // determine the dimension Dim = n_dof_types / 3; // Recast Jacobian matrix to CRDoubleMatrix CRDoubleMatrix* cr_matrix_pt = dynamic_cast(matrix_pt()); #ifdef PARANOID if (cr_matrix_pt == 0) { std::ostringstream error_message; error_message << "FSIPreconditioner only works with" << " CRDoubleMatrix matrices" << std::endl; throw OomphLibError( error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION); } #endif // Setup the block ordering. this->block_setup(); // Create dof_list for subsidiary preconditioner. This will be used later // in turn_into_subsidiary_block_preconditioner(...) Vector dof_list_for_subsidiary_prec(n_solid_dof_types); for (unsigned i = 0; i < n_solid_dof_types; i++) { dof_list_for_subsidiary_prec[i] = i; } // compute the scaling (if required) if (Use_inf_norm_of_s_scaling) { Vector dof_list(n_solid_dof_types); for (unsigned i = 0; i < n_solid_dof_types; i++) { dof_list[i] = i; } PseudoElasticPreconditionerScalingHelperOld* helper_pt = new PseudoElasticPreconditionerScalingHelperOld( this, cr_matrix_pt, dof_list, Elastic_mesh_pt, comm_pt()); Scaling = helper_pt->s_inf_norm(); delete helper_pt; helper_pt = 0; } else { Scaling = 1.0; } // setup the solid subsidiary preconditioner // //////////////////////////////// this preconditioner uses the full S // matrix if (E_preconditioner_type == Exact_block_preconditioner) { PseudoElasticPreconditionerSubsidiaryPreconditionerOld* s_prec_pt = new PseudoElasticPreconditionerSubsidiaryPreconditionerOld; Vector dof_list(n_solid_dof_types); for (unsigned i = 0; i < n_solid_dof_types; i++) { dof_list[i] = i; } s_prec_pt->turn_into_subsidiary_block_preconditioner(this, dof_list); if (Elastic_subsidiary_preconditioner_function_pt != 0) { s_prec_pt->set_subsidiary_preconditioner_function( Elastic_subsidiary_preconditioner_function_pt); } s_prec_pt->scaling() = Scaling; s_prec_pt->Preconditioner::setup(matrix_pt()); Elastic_preconditioner_pt = s_prec_pt; } // otherwise it is a block based preconditioner else { // create the preconditioner PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld* s_prec_pt = new PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld; Vector dof_list(n_solid_dof_types); for (unsigned i = 0; i < n_solid_dof_types; i++) { dof_list[i] = i; } s_prec_pt->turn_into_subsidiary_block_preconditioner(this, dof_list); // set the subsidiary solve method if (Elastic_subsidiary_preconditioner_function_pt != 0) { s_prec_pt->set_subsidiary_preconditioner_function( Elastic_subsidiary_preconditioner_function_pt); } // set the scaling s_prec_pt->scaling() = Scaling; // set the block preconditioning method switch (E_preconditioner_type) { case Block_diagonal_preconditioner: s_prec_pt->use_block_diagonal_approximation(); break; case Block_upper_triangular_preconditioner: s_prec_pt->use_upper_triangular_approximation(); break; case Block_lower_triangular_preconditioner: s_prec_pt->use_lower_triangular_approximation(); break; default: break; } // setup s_prec_pt->Preconditioner::setup(matrix_pt()); Elastic_preconditioner_pt = s_prec_pt; } // next setup the lagrange multiplier preconditioners Lagrange_multiplier_preconditioner_pt.resize(Dim); for (unsigned d = 0; d < Dim; d++) { CRDoubleMatrix* b_pt = new CRDoubleMatrix; this->get_block(2 * Dim + d, Dim + d, *b_pt); // if a non default preconditioner is specified create // the preconditioners if (Lagrange_multiplier_subsidiary_preconditioner_function_pt != 0) { Lagrange_multiplier_preconditioner_pt[d] = (*Lagrange_multiplier_subsidiary_preconditioner_function_pt)(); } // else use default superlu preconditioner else { Lagrange_multiplier_preconditioner_pt[d] = ExactPreconditionerFactory::create_exact_preconditioner(); } // and setup Lagrange_multiplier_preconditioner_pt[d]->setup(b_pt); delete b_pt; b_pt = 0; } } //============================================================================= /// Apply the elastic subsidiary preconditioner. //============================================================================= void PseudoElasticPreconditionerOld::elastic_preconditioner_solve( const DoubleVector& r, DoubleVector& z) { // apply the solid preconditioner Elastic_preconditioner_pt->preconditioner_solve(r, z); } //============================================================================= /// Apply the lagrange multiplier subsidiary preconditioner. //============================================================================= void PseudoElasticPreconditionerOld::lagrange_multiplier_preconditioner_solve( const DoubleVector& r, DoubleVector& z) { // apply the lagrange multiplier preconditioner for (unsigned d = 0; d < Dim; d++) { DoubleVector x; this->get_block_vector(Dim * 2 + d, r, x); DoubleVector y; Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(x, y); Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(y, x); unsigned nrow_local = x.nrow_local(); double* x_pt = x.values_pt(); for (unsigned i = 0; i < nrow_local; i++) { x_pt[i] = x_pt[i] * Scaling; } this->return_block_vector(Dim * 2 + d, x, z); } } //============================================================================= /// Clears the memory. //============================================================================= void PseudoElasticPreconditionerOld::clean_up_memory() { // clean the block preconditioner base class memory this->clear_block_preconditioner_base(); // delete the solid preconditioner delete Elastic_preconditioner_pt; Elastic_preconditioner_pt = 0; // delete the lagrange multiplier preconditioner pt unsigned sz = Lagrange_multiplier_preconditioner_pt.size(); for (unsigned i = 0; i < sz; i++) { delete Lagrange_multiplier_preconditioner_pt[i]; Lagrange_multiplier_preconditioner_pt[i] = 0; } } /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// //============================================================================= /// Setup the preconditioner //============================================================================= void PseudoElasticPreconditionerSubsidiaryPreconditionerOld::setup() { // clean memory this->clean_up_memory(); #ifdef PARANOID // paranoid check that this preconditioner has an even number of DOF types if (this->ndof_types() % 2 != 0) { std::ostringstream error_message; error_message << "This SUBSIDIARY preconditioner requires an even number of " << "types of DOF"; throw OomphLibError( error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION); } #endif // assemble dof_to_block_map unsigned ndof_types = this->ndof_types(); Vector dof_to_block_map(ndof_types, 0); for (unsigned i = ndof_types / 2; i < ndof_types; i++) { dof_to_block_map[i] = 1; } this->block_setup(dof_to_block_map); // get block 11 CRDoubleMatrix* s11_pt = new CRDoubleMatrix; this->get_block(1, 1, *s11_pt); // add the scaled identity matrix to block 11 double* s11_values = s11_pt->value(); int* s11_column_index = s11_pt->column_index(); int* s11_row_start = s11_pt->row_start(); int s11_nrow_local = s11_pt->nrow_local(); int s11_first_row = s11_pt->first_row(); for (int i = 0; i < s11_nrow_local; i++) { bool found = false; for (int j = s11_row_start[i]; j < s11_row_start[i + 1] && !found; j++) { if (s11_column_index[j] == i + s11_first_row) { s11_values[j] += Scaling; found = true; } } } VectorMatrix required_blocks(2, 2); const bool want_block = true; for (unsigned b_i = 0; b_i < 2; b_i++) { for (unsigned b_j = 0; b_j < 2; b_j++) { required_blocks[b_i][b_j].select_block(b_i, b_j, want_block); } } required_blocks[1][1].set_replacement_block_pt(s11_pt); CRDoubleMatrix s_prec_pt = this->get_concatenated_block(required_blocks); delete s11_pt; s11_pt = 0; // setup the preconditioner if (Subsidiary_preconditioner_function_pt != 0) { Preconditioner_pt = (*Subsidiary_preconditioner_function_pt)(); } else { Preconditioner_pt = ExactPreconditionerFactory::create_exact_preconditioner(); } Preconditioner_pt->setup(&s_prec_pt); } //============================================================================= /// Apply the preconditioner. //============================================================================= void PseudoElasticPreconditionerSubsidiaryPreconditionerOld:: preconditioner_solve(const DoubleVector& r, DoubleVector& z) { DoubleVector x; this->get_block_ordered_preconditioner_vector(r, x); DoubleVector y; Preconditioner_pt->preconditioner_solve(x, y); this->return_block_ordered_preconditioner_vector(y, z); } /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// //============================================================================= /// clean up the memory //============================================================================= void PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld:: clean_up_memory() { // number of block types unsigned n_block = Diagonal_block_preconditioner_pt.size(); // delete diagonal blocks for (unsigned i = 0; i < n_block; i++) { delete Diagonal_block_preconditioner_pt[i]; Diagonal_block_preconditioner_pt[i] = 0; if (Method == 1) { for (unsigned j = i + 1; j < n_block; j++) { delete Off_diagonal_matrix_vector_products(i, j); Off_diagonal_matrix_vector_products(i, j) = 0; } } else if (Method == 2) { for (unsigned j = 0; j < i; j++) { delete Off_diagonal_matrix_vector_products(i, j); Off_diagonal_matrix_vector_products(i, j) = 0; } } } // clean up the block preconditioner this->clear_block_preconditioner_base(); } //============================================================================= /// Setup the preconditioner. //============================================================================= void PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld::setup() { // clean the memory this->clean_up_memory(); // determine the number of DOF types unsigned n_dof_types = this->ndof_types(); #ifdef PARANOID // must be Dim*2 dof types if (n_dof_types % 2 != 0) { std::ostringstream error_message; error_message << "This preconditioner requires DIM*3 types of DOF"; throw OomphLibError( error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION); } #endif // store the dimension of the problem unsigned dim = n_dof_types / 2; // assemble the dof to block lookup scheme Vector dof_to_block_map(n_dof_types, 0); for (unsigned d = 0; d < dim; d++) { dof_to_block_map[d] = d; dof_to_block_map[d + dim] = d; } // setup the blocks look up schemes this->block_setup(dof_to_block_map); // Storage for the diagonal block preconditioners Diagonal_block_preconditioner_pt.resize(dim); // storage for the off diagonal matrix vector products Off_diagonal_matrix_vector_products.resize(dim, dim, 0); // setup the subsidiary preconditioners for (unsigned d = 0; d < dim; d++) { Vector dof_list(2); dof_list[0] = d; dof_list[1] = d + dim; Diagonal_block_preconditioner_pt[d] = new PseudoElasticPreconditionerSubsidiaryPreconditionerOld; Diagonal_block_preconditioner_pt[d] ->turn_into_subsidiary_block_preconditioner(this, dof_list); if (Subsidiary_preconditioner_function_pt != 0) { Diagonal_block_preconditioner_pt[d] ->set_subsidiary_preconditioner_function( Subsidiary_preconditioner_function_pt); } Diagonal_block_preconditioner_pt[d]->scaling() = Scaling; Diagonal_block_preconditioner_pt[d]->Preconditioner::setup(matrix_pt()); // the preconditioning method.\n // 0 - block diagonal\n // 1 - upper triangular\n // 2 - lower triangular\n // next setup the off diagonal mat vec operators if required if (Method == 1 || Method == 2) { unsigned l = d + 1; unsigned u = dim; if (Method == 2) { l = 0; u = d; } for (unsigned j = l; j < u; j++) { CRDoubleMatrix* block_matrix_pt = new CRDoubleMatrix; this->get_block(d, j, *block_matrix_pt); Off_diagonal_matrix_vector_products(d, j) = new MatrixVectorProduct(); // Off_diagonal_matrix_vector_products(d,j)->setup(block_matrix_pt); this->setup_matrix_vector_product( Off_diagonal_matrix_vector_products(d, j), block_matrix_pt, j); delete block_matrix_pt; block_matrix_pt = 0; } } } // setup the subsidiary preconditioner. } // preconditioner setup //============================================================================= /// Apply preconditioner to r //============================================================================= void PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld:: preconditioner_solve(const DoubleVector& res, DoubleVector& z) { // copy r DoubleVector r(res); unsigned n_block; // Cache umber of block types (also the spatial DIM) n_block = this->nblock_types(); // loop parameters int start = n_block - 1; int end = -1; int step = -1; if (Method != 1) { start = 0; end = n_block; step = 1; } // the preconditioning method. // 0 - block diagonal // 1 - upper triangular // 2 - lower triangular // // loop over the DIM // // For Method = 0 or 2 (diagonal, lower) // start = 2, end = -1, step = -1 // i = 2,1,0 // // For Method = 1 (upper) // start = 0, end = 3 step = 1 // i = 0, 1, 2 for (int i = start; i != end; i += step) { // solve Diagonal_block_preconditioner_pt[i]->preconditioner_solve(r, z); // if upper or lower triangular if (Method != 0) { // substitute // for (int j = i + step; j != end; j += step) { DoubleVector x; this->get_block_vector(i, z, x); DoubleVector y; Off_diagonal_matrix_vector_products(j, i)->multiply(x, y); x.clear(); this->get_block_vector(j, r, x); x -= y; this->return_block_vector(j, x, r); } // substitute } // if upper or lower } // for loop over DIM } // Block preconditioner solve } // namespace oomph