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// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 
28 
29 namespace oomph
30 {
31  //=============================================================================
32  /// clean up memory method
33  //=============================================================================
35  {
36  // wipe the preconditioner
38  Navier_stokes_preconditioner_pt->clean_up_memory();
40  Solid_preconditioner_pt->clean_up_memory();
41 
42  // clean the subsidiary matvec operators
43  Fluid_pseudo_elastic_matvec_pt->clean_up_memory();
44  Solid_fluid_matvec_pt->clean_up_memory();
45  Solid_pseudo_elastic_matvec_pt->clean_up_memory();
46  Lagrange_solid_matvec_pt->clean_up_memory();
47  }
48 
49  //=============================================================================
50  /// Setup the precoonditioner.
51  //=============================================================================
53  {
54  // clean the memory
55  this->clean_up_memory();
56 
57 #ifdef PARANOID
58  // paranoid check that the meshes have been set
60  {
61  std::ostringstream error_message;
62  error_message << "The fluid and pseudo elastic mesh must be set.";
63  throw OomphLibError(
64  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
65  }
66  if (Solid_mesh_pt == 0)
67  {
68  std::ostringstream error_message;
69  error_message << "The solid mesh must be set.";
70  throw OomphLibError(
71  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
72  }
74  {
75  std::ostringstream error_message;
76  error_message << "The Lagrange multiplier mesh must be set.";
77  throw OomphLibError(
78  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
79  }
80 #endif
81 
82  // add the meshes
83  this->set_mesh(0, Fluid_and_pseudo_elastic_mesh_pt);
84  this->set_mesh(1, Solid_mesh_pt);
85  this->set_mesh(2, Lagrange_multiplier_mesh_pt);
86 
87  // determine the number of fluid dofs
88  unsigned nfluid_dof = Dim + 1;
89 
90  // determine the number of pseudo solid dofs
91  unsigned npseudo_elastic_dof = this->ndof_types_in_mesh(0) - nfluid_dof;
92 
93  // determine the number of solid dofs
94  unsigned nsolid_dof = this->ndof_types_in_mesh(1);
95 
96  // determine the number of lagrange multiplier dofs
97  unsigned nlagr_mult_dof = this->ndof_types_in_mesh(2);
98 
99  // total number of dof types
100  unsigned ntotal_dof =
101  nfluid_dof + npseudo_elastic_dof + nsolid_dof + nlagr_mult_dof;
102 
103  // setup the block lookup scheme
104  // block 0 - fluid
105  // 1 - solid
106  // 2 - pseudo solid inc. lagrange mult
107  Vector<unsigned> dof_to_block_map(ntotal_dof, 0);
108  int c = nfluid_dof;
109  for (unsigned i = 0; i < npseudo_elastic_dof; i++)
110  {
111  dof_to_block_map[c] = 2;
112  c++;
113  }
114  for (unsigned i = 0; i < nsolid_dof; i++)
115  {
116  dof_to_block_map[c] = 1;
117  c++;
118  }
119  for (unsigned i = 0; i < nlagr_mult_dof; i++)
120  {
121  dof_to_block_map[c] = 3;
122  c++;
123  }
124 
125  // Recast Jacobian matrix to CRDoubleMatrix
126  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
127 #ifdef PARANOID
128  if (cr_matrix_pt == 0)
129  {
130  std::ostringstream error_message;
131  error_message << "FSIPreconditioner only works with"
132  << " CRDoubleMatrix matrices" << std::endl;
133  throw OomphLibError(
134  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
135  }
136 #endif
137 
138  // Call block setup for this preconditioner
139  this->block_setup(dof_to_block_map);
140 
141  // SETUP THE PRECONDITIONERS
142  // =========================
143 
144  // setup the navier stokes preconditioner
146  {
147  Navier_stokes_schur_complement_preconditioner_pt->set_navier_stokes_mesh(
149 
150  Vector<unsigned> ns_dof_list(nfluid_dof, 0);
151  for (unsigned i = 0; i < nfluid_dof; i++)
152  {
153  ns_dof_list[i] = i;
154  }
155 
157  ->turn_into_subsidiary_block_preconditioner(this, ns_dof_list);
158 
160  matrix_pt());
161  }
162  else
163  {
164  CRDoubleMatrix* ns_matrix_pt = new CRDoubleMatrix;
165  this->get_block(0, 0, *ns_matrix_pt);
166 
167  Navier_stokes_preconditioner_pt->setup(ns_matrix_pt);
168  delete ns_matrix_pt;
169  ns_matrix_pt = 0;
170  }
171 
172  // next the solid preconditioner
173  if (dynamic_cast<BlockPreconditioner<CRDoubleMatrix>*>(
175  {
177  GeneralPurposeBlockPreconditioner<CRDoubleMatrix>*
178  solid_block_preconditioner_pt =
179  dynamic_cast<GeneralPurposeBlockPreconditioner<CRDoubleMatrix>*>(
181 
182  if (solid_block_preconditioner_pt != 0)
183  {
184  unsigned offset = nfluid_dof + npseudo_elastic_dof;
185  Vector<unsigned> solid_prec_dof_list(nsolid_dof);
186  for (unsigned i = 0; i < nsolid_dof; i++)
187  {
188  solid_prec_dof_list[i] = offset + i;
189  }
190  solid_block_preconditioner_pt
191  ->turn_into_subsidiary_block_preconditioner(this,
192  solid_prec_dof_list);
193  solid_block_preconditioner_pt->setup(cr_matrix_pt);
194  }
195  else
196  {
197  std::ostringstream error_message;
198  error_message << "If the (real) solid preconditioner is a "
199  << "BlockPreconditioner then is must be a "
200  << "GeneralPurposeBlockPreconditioner";
201  throw OomphLibError(error_message.str(),
202  OOMPH_CURRENT_FUNCTION,
203  OOMPH_EXCEPTION_LOCATION);
204  }
205  }
206  // otherwise it is a solid preconditioner
207  else
208  {
210  CRDoubleMatrix* s_matrix_pt = new CRDoubleMatrix;
211  this->get_block(1, 1, *s_matrix_pt);
212  Solid_preconditioner_pt->setup(s_matrix_pt);
213  delete s_matrix_pt;
214  s_matrix_pt = 0;
215  }
216 
217  // next the pseudo solid preconditioner
218  unsigned ndof_for_pseudo_elastic_prec = Dim * 3;
219  Vector<unsigned> pseudo_elastic_prec_dof_list(ndof_for_pseudo_elastic_prec,
220  0);
221  for (unsigned i = 0; i < Dim * 2; i++)
222  {
223  pseudo_elastic_prec_dof_list[i] = nfluid_dof + i;
224  }
225  for (unsigned i = 0; i < Dim; i++)
226  {
227  pseudo_elastic_prec_dof_list[i + Dim * 2] =
228  nfluid_dof + npseudo_elastic_dof + nsolid_dof + i;
229  }
230  Pseudo_elastic_preconditioner_pt->turn_into_subsidiary_block_preconditioner(
231  this, pseudo_elastic_prec_dof_list);
236  Pseudo_elastic_preconditioner_pt->Preconditioner::setup(matrix_pt());
237 
238  // SETUP THE MATRIX VECTOR PRODUCT OPERATORS
239  // =========================================
240 
241  // setup the fluid pseudo-solid matvec operator
242  CRDoubleMatrix* fp_matrix_pt = new CRDoubleMatrix;
243  get_block(0, 2, *fp_matrix_pt);
244  // Fluid_pseudo_elastic_matvec_pt->setup(fp_matrix_pt);
245  this->setup_matrix_vector_product(
246  Fluid_pseudo_elastic_matvec_pt, fp_matrix_pt, 2);
247  delete fp_matrix_pt;
248  fp_matrix_pt = 0;
249 
250  // setup the solid fluid matvec operator
251  CRDoubleMatrix* sf_matrix_pt = new CRDoubleMatrix;
252  get_block(1, 0, *sf_matrix_pt);
253  // Solid_fluid_matvec_pt->setup(sf_matrix_pt);
254  this->setup_matrix_vector_product(Solid_fluid_matvec_pt, sf_matrix_pt, 0);
255  delete sf_matrix_pt;
256  sf_matrix_pt = 0;
257 
258  // setup the solid pseudo-solid matvec operator
259  CRDoubleMatrix* sp_matrix_pt = new CRDoubleMatrix;
260  get_block(1, 2, *sp_matrix_pt);
261  // Solid_pseudo_elastic_matvec_pt->setup(sp_matrix_pt);
262  this->setup_matrix_vector_product(
263  Solid_pseudo_elastic_matvec_pt, sp_matrix_pt, 2);
264  delete sp_matrix_pt;
265  sp_matrix_pt = 0;
266 
267  // build the lagrange solid matvec operator
268  CRDoubleMatrix* ls_matrix_pt = new CRDoubleMatrix;
269  get_block(3, 1, *ls_matrix_pt);
270  // Lagrange_solid_matvec_pt->setup(ls_matrix_pt);
271  this->setup_matrix_vector_product(
272  Lagrange_solid_matvec_pt, ls_matrix_pt, 1);
273  delete ls_matrix_pt;
274  ls_matrix_pt = 0;
275  }
276 
277  //=============================================================================
278  /// Apply the preconditioner
279  //=============================================================================
281  const DoubleVector& r, DoubleVector& z)
282  {
283  // apply the "pseudo solid" component of the pseudo solid preconditioner
285 
286  // apply the fluid on pseudo solid matrix vector product operator
287  DoubleVector x;
288  this->get_block_vector(2, z, x);
289  DoubleVector y;
290  Fluid_pseudo_elastic_matvec_pt->multiply(x, y);
291  DoubleVector w;
292  Solid_pseudo_elastic_matvec_pt->multiply(x, w);
293  x.clear();
294  this->get_block_vector(0, r, x);
295  x -= y;
296  y.clear();
297 
298  // storage for a copy of z
299  DoubleVector z_copy;
300 
301  // apply the ns preconditioner
303  {
304  z_copy.build(z);
305  this->return_block_vector(0, x, z_copy);
306  x.clear();
308  z_copy, z);
309  z_copy.clear();
310  }
311  else
312  {
313  Navier_stokes_preconditioner_pt->preconditioner_solve(x, y);
314  x.clear();
315  this->return_block_vector(0, y, z);
316  y.clear();
317  }
318 
319  // apply the solid onto fluid matrix vector product operator
320  this->get_block_vector(0, z, x);
321  Solid_fluid_matvec_pt->multiply(x, y);
322  x.clear();
323  this->get_block_vector(1, r, x);
324  x -= y;
325  y.clear();
326 
327  // apply the result of the solid onto pseudo solid matrix vector product
328  x -= w;
329  w.clear();
330 
331  // apply the solid preconditioner
333  {
334  DoubleVector z_copy(z);
335  this->return_block_vector(1, x, z_copy);
336  x.clear();
337  (dynamic_cast<GeneralPurposeBlockPreconditioner<CRDoubleMatrix>*>(
339  ->preconditioner_solve(z_copy, z);
340  this->get_block_vector(1, z, y);
341  }
342  else
343  {
344  Solid_preconditioner_pt->preconditioner_solve(x, y);
345  x.clear();
346  this->return_block_vector(1, y, z);
347  }
348 
349  // apply the lagrange multiplier solid matrix vector product operator
350  Lagrange_solid_matvec_pt->multiply(y, x);
351  y.clear();
352  this->get_block_vector(3, r, y);
353  y -= x;
354  x.clear();
355  z_copy.build(z);
356  this->return_block_vector(3, y, z_copy);
357 
358  // apply the lagrange multiplier compenent of the pseudo solid
359  // preconditioner
361  z_copy, z);
362  z_copy.clear();
363  }
364 } // namespace oomph
MatrixVectorProduct * Solid_pseudo_elastic_matvec_pt
solid onto pseudo solid matrix vector operatio
Mesh * Lagrange_multiplier_mesh_pt
Mesh containing the lagrange multiplier elements.
Preconditioner * Solid_preconditioner_pt
pointer to the solid preconditioner
NavierStokesSchurComplementPreconditioner * Navier_stokes_schur_complement_preconditioner_pt
Navier Stokes Schur complement preconditioner.
Preconditioner * Navier_stokes_preconditioner_pt
pointer to the navier stokes precondtioner
MatrixVectorProduct * Fluid_pseudo_elastic_matvec_pt
fluid onto pseudosolid matrix vector operator
MatrixVectorProduct * Solid_fluid_matvec_pt
solid onto fluid matrix vector operatio
Mesh * Fluid_and_pseudo_elastic_mesh_pt
Mesh containing the combined fluid and pseudo solid element.
bool Solid_preconditioner_is_block_preconditioner
boolean flag to indicate whether the Solid preconditioner is a block preconditioner
bool Use_navier_stokes_schur_complement_preconditioner
If true the Navier Stokes Schur complement preconditioner is used. Otherwise SuperLUPreconditioner is...
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply the preconditioner.
PseudoElasticPreconditioner * Pseudo_elastic_preconditioner_pt
pointer to the pseudo solid preconditioner
Mesh * Solid_mesh_pt
Mesh containing the solid elements.
void lagrange_multiplier_preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply the lagrange multiplier subsidiary preconditioner.
void set_lagrange_multiplier_mesh(Mesh *mesh_pt)
Access function to mesh containing the block-preconditionable lagrange multiplier elements.
void elastic_preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply the elastic subsidiary preconditioner.
void set_elastic_mesh(Mesh *mesh_pt)
Access function to mesh containing the block-preconditionable elastic elements.