pseudo_elastic_fsi_preconditioner.cc
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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 
32 
33 namespace oomph {
34 
35  //=============================================================================
36  /// clean up memory method
37  //=============================================================================
39  {
40  // wipe the preconditioner
45 
46  // clean the subsidiary matvec operators
51  }
52 
53  //=============================================================================
54  /// \short Setup the precoonditioner.
55  //=============================================================================
57  {
58  // clean the memory
59  this->clean_up_memory();
60 
61 #ifdef PARANOID
62  // paranoid check that the meshes have been set
64  {
65  std::ostringstream error_message;
66  error_message << "The fluid and pseudo elastic mesh must be set.";
67  throw OomphLibError(error_message.str(),
68  OOMPH_CURRENT_FUNCTION,
69  OOMPH_EXCEPTION_LOCATION);
70 
71  }
72  if (Solid_mesh_pt==0)
73  {
74  std::ostringstream error_message;
75  error_message << "The solid mesh must be set.";
76  throw OomphLibError(error_message.str(),
77  OOMPH_CURRENT_FUNCTION,
78  OOMPH_EXCEPTION_LOCATION);
79  }
81  {
82  std::ostringstream error_message;
83  error_message << "The Lagrange multiplier mesh must be set.";
84  throw OomphLibError(error_message.str(),
85  OOMPH_CURRENT_FUNCTION,
86  OOMPH_EXCEPTION_LOCATION);
87  }
88 #endif
89 
90  // add the meshes
92  this->set_mesh(1,Solid_mesh_pt);
94 
95  // determine the number of fluid dofs
96  unsigned nfluid_dof = Dim + 1;
97 
98  // determine the number of pseudo solid dofs
99  unsigned npseudo_elastic_dof = this->ndof_types_in_mesh(0)-nfluid_dof;
100 
101  // determine the number of solid dofs
102  unsigned nsolid_dof = this->ndof_types_in_mesh(1);
103 
104  // determine the number of lagrange multiplier dofs
105  unsigned nlagr_mult_dof = this->ndof_types_in_mesh(2);
106 
107  // total number of dof types
108  unsigned ntotal_dof = nfluid_dof + npseudo_elastic_dof + nsolid_dof
109  + nlagr_mult_dof;
110 
111  // setup the block lookup scheme
112  // block 0 - fluid
113  // 1 - solid
114  // 2 - pseudo solid inc. lagrange mult
115  Vector<unsigned> dof_to_block_map(ntotal_dof,0);
116  int c = nfluid_dof;
117  for (unsigned i = 0; i < npseudo_elastic_dof; i++)
118  {
119  dof_to_block_map[c] = 2;
120  c++;
121  }
122  for (unsigned i = 0; i < nsolid_dof; i++)
123  {
124  dof_to_block_map[c] = 1;
125  c++;
126  }
127  for (unsigned i = 0; i < nlagr_mult_dof; i++)
128  {
129  dof_to_block_map[c] = 3;
130  c++;
131  }
132 
133  // Recast Jacobian matrix to CRDoubleMatrix
134  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
135 #ifdef PARANOID
136  if (cr_matrix_pt==0)
137  {
138  std::ostringstream error_message;
139  error_message << "FSIPreconditioner only works with"
140  << " CRDoubleMatrix matrices" << std::endl;
141  throw OomphLibError(error_message.str(),
142  OOMPH_CURRENT_FUNCTION,
143  OOMPH_EXCEPTION_LOCATION);
144  }
145 #endif
146 
147  // Call block setup for this preconditioner
148  this->block_setup(dof_to_block_map);
149 
150  // SETUP THE PRECONDITIONERS
151  // =========================
152 
153  // setup the navier stokes preconditioner
155  {
157  set_navier_stokes_mesh(Fluid_and_pseudo_elastic_mesh_pt);
158 
159  Vector<unsigned> ns_dof_list(nfluid_dof,0);
160  for (unsigned i = 0; i < nfluid_dof; i++)
161  {
162  ns_dof_list[i] = i;
163  }
164 
167 
170  }
171  else
172  {
173  CRDoubleMatrix* ns_matrix_pt = new CRDoubleMatrix;
174  this->get_block(0,0,*ns_matrix_pt);
175 
176  Navier_stokes_preconditioner_pt->setup(ns_matrix_pt);
177  delete ns_matrix_pt; ns_matrix_pt = 0;
178  }
179 
180  // next the solid preconditioner
181  if (dynamic_cast<BlockPreconditioner<CRDoubleMatrix>* >
183  {
186  solid_block_preconditioner_pt =
189 
190  if (solid_block_preconditioner_pt != 0)
191  {
192  unsigned offset = nfluid_dof+npseudo_elastic_dof;
193  Vector<unsigned> solid_prec_dof_list(nsolid_dof);
194  for (unsigned i = 0; i < nsolid_dof; i++)
195  {
196  solid_prec_dof_list[i]=offset+i;
197  }
198  solid_block_preconditioner_pt
200  solid_prec_dof_list);
201  solid_block_preconditioner_pt->setup(cr_matrix_pt);
202  }
203  else
204  {
205  std::ostringstream error_message;
206  error_message << "If the (real) solid preconditioner is a "
207  << "BlockPreconditioner then is must be a "
208  << "GeneralPurposeBlockPreconditioner";
209  throw OomphLibError(error_message.str(),
210  OOMPH_CURRENT_FUNCTION,
211  OOMPH_EXCEPTION_LOCATION);
212  }
213  }
214  // otherwise it is a solid preconditioner
215  else
216  {
218  CRDoubleMatrix* s_matrix_pt = new CRDoubleMatrix;
219  this->get_block(1,1,*s_matrix_pt);
220  Solid_preconditioner_pt->setup(s_matrix_pt);
221  delete s_matrix_pt; s_matrix_pt = 0;
222  }
223 
224  // next the pseudo solid preconditioner
225  unsigned ndof_for_pseudo_elastic_prec = Dim*3;
226  Vector<unsigned> pseudo_elastic_prec_dof_list(ndof_for_pseudo_elastic_prec,0);
227  for (unsigned i = 0; i < Dim*2; i++)
228  {
229  pseudo_elastic_prec_dof_list[i] = nfluid_dof+i;
230  }
231  for (unsigned i = 0; i < Dim; i++)
232  {
233  pseudo_elastic_prec_dof_list[i+Dim*2] = nfluid_dof
234  + npseudo_elastic_dof + nsolid_dof + i;
235  }
238  pseudo_elastic_prec_dof_list);
240  set_elastic_mesh(this->Fluid_and_pseudo_elastic_mesh_pt);
242  set_lagrange_multiplier_mesh(this->Lagrange_multiplier_mesh_pt);
245 
246  // SETUP THE MATRIX VECTOR PRODUCT OPERATORS
247  // =========================================
248 
249  // setup the fluid pseudo-solid matvec operator
250  CRDoubleMatrix* fp_matrix_pt = new CRDoubleMatrix;
251  get_block(0,2,*fp_matrix_pt);
252 // Fluid_pseudo_elastic_matvec_pt->setup(fp_matrix_pt);
254  fp_matrix_pt,2);
255  delete fp_matrix_pt; fp_matrix_pt = 0;
256 
257  // setup the solid fluid matvec operator
258  CRDoubleMatrix* sf_matrix_pt = new CRDoubleMatrix;
259  get_block(1,0,*sf_matrix_pt);
260 // Solid_fluid_matvec_pt->setup(sf_matrix_pt);
262  delete sf_matrix_pt; sf_matrix_pt = 0;
263 
264  // setup the solid pseudo-solid matvec operator
265  CRDoubleMatrix* sp_matrix_pt = new CRDoubleMatrix;
266  get_block(1,2,*sp_matrix_pt);
267 // Solid_pseudo_elastic_matvec_pt->setup(sp_matrix_pt);
269  sp_matrix_pt,2);
270  delete sp_matrix_pt; sp_matrix_pt = 0;
271 
272  // build the lagrange solid matvec operator
273  CRDoubleMatrix* ls_matrix_pt = new CRDoubleMatrix;
274  get_block(3,1,*ls_matrix_pt);
275 // Lagrange_solid_matvec_pt->setup(ls_matrix_pt);
277  ls_matrix_pt,1);
278  delete ls_matrix_pt; ls_matrix_pt = 0;
279 
280  }
281 
282  //=============================================================================
283  /// \short Apply the preconditioner
284  //=============================================================================
287  {
288  // apply the "pseudo solid" component of the pseudo solid preconditioner
289  Pseudo_elastic_preconditioner_pt->elastic_preconditioner_solve(r,z);
290 
291  // apply the fluid on pseudo solid matrix vector product operator
292  DoubleVector x;
293  this->get_block_vector(2,z,x);
294  DoubleVector y;
295  Fluid_pseudo_elastic_matvec_pt->multiply(x,y);
296  DoubleVector w;
297  Solid_pseudo_elastic_matvec_pt->multiply(x,w);
298  x.clear();
299  this->get_block_vector(0,r,x);
300  x -= y;
301  y.clear();
302 
303  // storage for a copy of z
304  DoubleVector z_copy;
305 
306  // apply the ns preconditioner
307  if (Use_navier_stokes_schur_complement_preconditioner)
308  {
309  z_copy.build(z);
310  this->return_block_vector(0,x,z_copy);
311  x.clear();
312  Navier_stokes_schur_complement_preconditioner_pt->
313  preconditioner_solve(z_copy,z);
314  z_copy.clear();
315  }
316  else
317  {
318  Navier_stokes_preconditioner_pt->preconditioner_solve(x,y);
319  x.clear();
320  this->return_block_vector(0,y,z);
321  y.clear();
322  }
323 
324  // apply the solid onto fluid matrix vector product operator
325  this->get_block_vector(0,z,x);
326  Solid_fluid_matvec_pt->multiply(x,y);
327  x.clear();
328  this->get_block_vector(1,r,x);
329  x -= y;
330  y.clear();
331 
332  // apply the result of the solid onto pseudo solid matrix vector product
333  x -= w;
334  w.clear();
335 
336  // apply the solid preconditioner
337  if (Solid_preconditioner_is_block_preconditioner)
338  {
339  DoubleVector z_copy(z);
340  this->return_block_vector(1,x,z_copy);
341  x.clear();
343  (Solid_preconditioner_pt))->preconditioner_solve(z_copy,z);
344  this->get_block_vector(1,z,y);
345  }
346  else
347  {
348  Solid_preconditioner_pt->preconditioner_solve(x,y);
349  x.clear();
350  this->return_block_vector(1,y,z);
351  }
352 
353  // apply the lagrange multiplier solid matrix vector product operator
354  Lagrange_solid_matvec_pt->multiply(y,x);
355  y.clear();
356  this->get_block_vector(3,r,y);
357  y -= x;
358  x.clear();
359  z_copy.build(z);
360  this->return_block_vector(3,y,z_copy);
361 
362  // apply the lagrange multiplier compenent of the pseudo solid
363  // preconditioner
364  Pseudo_elastic_preconditioner_pt->
365  lagrange_multiplier_preconditioner_solve(z_copy,z);
366  z_copy.clear();
367  }
368 }
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply the preconditioner.
MatrixVectorProduct * Solid_fluid_matvec_pt
solid onto fluid matrix vector operatio
void clear()
wipes the DoubleVector
cstr elem_len * i
Definition: cfortran.h:607
void clean_up_memory()
Broken assignment operator.
void clean_up_memory()
Helper function to delete preconditioner data.
unsigned ndof_types_in_mesh(const unsigned &i) const
Return the number of DOF types in mesh i. WARNING: This should only be used by the upper-most master ...
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, block_pt is a pointer to the block matrix, block_col_indices is a vector indicating which block indices does the RHS vector we want to multiply the matrix by.
Preconditioner * Navier_stokes_preconditioner_pt
pointer to the navier stokes precondtioner
NavierStokesSchurComplementPreconditioner * Navier_stokes_schur_complement_preconditioner_pt
Navier Stokes Schur complement preconditioner.
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
Mesh * Lagrange_multiplier_mesh_pt
Mesh containing the lagrange multiplier elements.
virtual void block_setup()
Determine the size of the matrix blocks and setup the lookup schemes relating the global degrees of f...
MatrixVectorProduct * Fluid_pseudo_elastic_matvec_pt
fluid onto pseudosolid matrix vector operator
bool Solid_preconditioner_is_block_preconditioner
boolean flag to indicate whether the Solid preconditioner is a block preconditioner ...
void set_mesh(const unsigned &i, const Mesh *const mesh_pt, const bool &allow_multiple_element_type_in_mesh=false)
Set the i-th mesh for this block preconditioner. Note: The method set_nmesh(...) must be called befor...
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...
Mesh * Solid_mesh_pt
Mesh containing the solid elements.
void setup(DoubleMatrixBase *matrix_pt)
Setup the preconditioner: store the matrix pointer and the communicator pointer then call preconditio...
MatrixVectorProduct * Solid_pseudo_elastic_matvec_pt
solid onto pseudo solid matrix vector operatio
Preconditioner * Solid_preconditioner_pt
pointer to the solid preconditioner
bool Use_navier_stokes_schur_complement_preconditioner
If true the Navier Stokes Schur complement preconditioner is used. Otherwise SuperLUPreconditioner is...
PseudoElasticPreconditioner * Pseudo_elastic_preconditioner_pt
pointer to the pseudo solid preconditioner
Mesh * Fluid_and_pseudo_elastic_mesh_pt
Mesh containing the combined fluid and pseudo solid element.
void turn_into_subsidiary_block_preconditioner(BlockPreconditioner< CRDoubleMatrix > *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 "...
virtual void setup()=0
Setup the preconditioner. Pure virtual generic interface function.
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
CRDoubleMatrix * matrix_pt() const
Access function to matrix_pt. If this is the master then cast the matrix pointer to MATRIX*...
void clean_up_memory()
clear the memory
virtual void clean_up_memory()
Clean up memory (empty). Generic interface function.
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872