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-2022 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
29namespace oomph
30{
31 //=============================================================================
32 /// clean up memory method
33 //=============================================================================
35 {
36 // wipe the preconditioner
41
42 // clean the subsidiary matvec operators
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
84 this->set_mesh(1, Solid_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 {
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
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
168 delete ns_matrix_pt;
169 ns_matrix_pt = 0;
170 }
171
172 // next the solid preconditioner
173 if (dynamic_cast<BlockPreconditioner<CRDoubleMatrix>*>(
175 {
178 solid_block_preconditioner_pt =
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
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 }
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);
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);
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);
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);
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;
291 DoubleVector 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 {
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);
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();
339 ->preconditioner_solve(z_copy, z);
340 this->get_block_vector(1, z, y);
341 }
342 else
343 {
345 x.clear();
346 this->return_block_vector(1, y, z);
347 }
348
349 // apply the lagrange multiplier solid matrix vector product operator
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
cstr elem_len * i
Definition: cfortran.h:603
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...
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 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...
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 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,...
CRDoubleMatrix * matrix_pt() const
Access function to matrix_pt. If this is the master then cast the matrix pointer to MATRIX*,...
virtual void block_setup()
Determine the size of the matrix blocks and setup the lookup schemes relating the global degrees of f...
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...
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
void clear()
wipes the DoubleVector
Base class for general purpose block preconditioners. Deals with setting subsidiary preconditioners a...
void clean_up_memory()
clear the memory
void multiply(const DoubleVector &x, DoubleVector &y) const
Apply the operator to the vector x and return the result in the vector y.
void clean_up_memory()
Helper function to delete preconditioner data.
void set_navier_stokes_mesh(Mesh *mesh_pt, const bool &allow_multiple_element_type_in_navier_stokes_mesh=false)
Specify the mesh containing the block-preconditionable Navier-Stokes elements. The optional argument ...
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to Vector r.
An OomphLibError object which should be thrown when an run-time error is encountered....
void setup(DoubleMatrixBase *matrix_pt)
Setup the preconditioner: store the matrix pointer and the communicator pointer then call preconditio...
virtual void clean_up_memory()
Clean up memory (empty). Generic interface function.
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
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.
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...