segregated_fsi_solver.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 
27 #include <algorithm>
28 
29 #include "segregated_fsi_solver.h"
30 #include <algorithm>
31 
32 namespace oomph
33 {
34  //============================================================================
35  /// Extrapolate solid data and update fluid mesh. Extrapolation is based
36  /// on previous history values and is therefore only called in
37  /// time-dependent runs.
38  //============================================================================
40  {
41  // Number of solid Data items:
42  unsigned n_data = Solid_data_pt.size();
43 
44  // Loop over all solid Data items:
45  for (unsigned i = 0; i < n_data; i++)
46  {
47  // Number of values stored at this Data item:
48  unsigned n_value = Solid_data_pt[i]->nvalue();
49 
50  // Number of values stored at this Data item:
51  unsigned n_prev = Solid_data_pt[i]->time_stepper_pt()->nprev_values();
52 
53  // Can't use this extrapolation if we don't have at least two
54  // history values; may be able to do better for higher order
55  // timesteppers, though.
56  if (n_prev >= 2)
57  {
58  // Extrapolate all values
59  for (unsigned k = 0; k < n_value; k++)
60  {
61  // Linear extrapolation based on previous two values,
62  // assuming constant timestep.
63  double new_value =
64  2.0 * Solid_data_pt[i]->value(1, k) - Solid_data_pt[i]->value(2, k);
65  *(Solid_data_pt[i]->value_pt(0, k)) = new_value;
66  }
67  }
68  }
69 
70  // Update fluid node position in response to changes in wall shape
71  // (If fluid mesh was specified via non-NULL Fluid_mesh_pt
72  // this node update will only act on the fluid nodes).
73  if (Fluid_mesh_pt != 0)
74  {
76  }
77  else
78  {
79  mesh_pt()->node_update();
80  }
81  }
82 
83 
84  //============================================================================
85  /// Rebuild global mesh for monolithic problem
86  //============================================================================
88  {
89  // Get rid of the previous submeshes
91 
92  // Add original submeshes
93  unsigned orig_n_sub_mesh = Orig_sub_mesh_pt.size();
94  for (unsigned i = 0; i < orig_n_sub_mesh; i++)
95  {
97  }
98 
99  // Rebuild global mesh
101  }
102 
103 
104  //============================================================================
105  /// Only include solid elements in the Problem's mesh. This is
106  /// called before the segregated solid solve. The solid elements are
107  /// identified by the user via the solid_mesh_pt argument
108  /// in the pure virtual function identify_fluid_and_solid_dofs(...).
109  /// If a NULL pointer is returned by this function (i.e. if the user
110  /// hasn't bothered to identify the solids elements in a submesh, then
111  /// no stripping of non-solid elements is performed. The result
112  /// of the computation will be correct but
113  /// it won't be very efficient.
114  //============================================================================
116  {
117  // Wipe global mesh and rebuild it with just the solid elements in it
118  if (Solid_mesh_pt != 0)
119  {
123  }
124  }
125 
126 
127  //============================================================================
128  /// Only include fluid elements in the Problem's mesh. This is
129  /// called before the segregated fluid solve. The fluid elements are
130  /// identified by the user via the fluid_mesh_pt argument
131  /// in the pure virtual function identify_fluid_and_solid_dofs(...).
132  /// If a NULL pointer is returned by this function (i.e. if the user
133  /// hasn't bothered to identify the fluids elements in a submesh, then
134  /// no stripping of non-fluid elements is performed. The result
135  /// of the computation will be correct but
136  /// it won't be very efficient.
137  //============================================================================
139  {
140  // Wipe global mesh and rebuild it with just the fluid elements in it
141  if (Fluid_mesh_pt != 0)
142  {
146  }
147  }
148 
149 
150  //============================================================================
151  /// Pin all fluid dofs
152  //============================================================================
154  {
155  // Number of fluid Data items:
156  unsigned n_data = Fluid_data_pt.size();
157 
158  // Loop over all fluid Data items:
159  for (unsigned i = 0; i < n_data; i++)
160  {
161  // Number of values stored at this Data item:
162  unsigned n_value = Fluid_data_pt[i]->nvalue();
163 
164  // Pin all values
165  for (unsigned k = 0; k < n_value; k++)
166  {
167  Fluid_data_pt[i]->pin(k);
168  }
169  }
170  }
171 
172 
173  //============================================================================
174  /// Restore the pinned status of all fluid dofs
175  //============================================================================
177  {
178  // Number of fluid Data items:
179  unsigned n_data = Fluid_data_pt.size();
180 
181  // Loop over all fluid Data items:
182  for (unsigned i = 0; i < n_data; i++)
183  {
184  // Number of values stored at this Data item:
185  unsigned n_value = Fluid_data_pt[i]->nvalue();
186 
187  // Pin all values
188  for (unsigned k = 0; k < n_value; k++)
189  {
190  if (Fluid_value_is_pinned[i][k])
191  {
192  Fluid_data_pt[i]->pin(k);
193  }
194  else
195  {
196  Fluid_data_pt[i]->unpin(k);
197  }
198  }
199  }
200  }
201 
202 
203  //============================================================================
204  /// Pin all solid dofs
205  //============================================================================
207  {
208  // Number of solid Data items:
209  unsigned n_data = Solid_data_pt.size();
210 
211  // Loop over all solid Data items:
212  for (unsigned i = 0; i < n_data; i++)
213  {
214  // Number of values stored at this Data item:
215  unsigned n_value = Solid_data_pt[i]->nvalue();
216 
217  // Pin all values
218  for (unsigned k = 0; k < n_value; k++)
219  {
220  Solid_data_pt[i]->pin(k);
221  }
222  }
223  }
224 
225 
226  //============================================================================
227  /// Restore the pinned status of all solid dofs
228  //============================================================================
230  {
231  // Number of solid Data items:
232  unsigned n_data = Solid_data_pt.size();
233 
234  // Loop over all solid Data items:
235  for (unsigned i = 0; i < n_data; i++)
236  {
237  // Number of values stored at this Data item:
238  unsigned n_value = Solid_data_pt[i]->nvalue();
239 
240  // Pin all values
241  for (unsigned k = 0; k < n_value; k++)
242  {
243  if (Solid_value_is_pinned[i][k])
244  {
245  Solid_data_pt[i]->pin(k);
246  }
247  else
248  {
249  Solid_data_pt[i]->unpin(k);
250  }
251  }
252  }
253  }
254 
255 
256  //============================================================================
257  /// Store the current solid values as reference values for future
258  /// convergence check. Also add another entry to pointwise Aitken history.
259  //============================================================================
261  {
262  // Counter for the number of values:
263  unsigned value_count = 0;
264 
265  // Number of solid Data items:
266  unsigned n_data = Solid_data_pt.size();
267 
268  // Loop over all solid Data items:
269  for (unsigned i = 0; i < n_data; i++)
270  {
271  // Number of values stored at this Data item:
272  unsigned n_value = Solid_data_pt[i]->nvalue();
273 
274  // Loop over all values
275  for (unsigned k = 0; k < n_value; k++)
276  {
277  // Make backup
278  Previous_solid_value[value_count] = Solid_data_pt[i]->value(k);
279 
280  // Store in pointwise Aitken history
282  {
284  Solid_data_pt[i]->value(k);
285  }
286 
287  // Increment counter
288  value_count++;
289  }
290  }
291 
292  // We stored another level of Aitken history values
294  }
295 
296 
297  //============================================================================
298  /// Get rms of change in the solid dofs; the max. change of the
299  /// solid dofs and the rms norm of the solid dofs themselves.
300  /// Change is computed relative to the reference values stored when
301  /// store_solid_dofs() was last called.
302  //============================================================================
304  double& max_change,
305  double& rms_norm)
306  {
307  // Initialise
308  rms_change = 0.0;
309  max_change = 0.0;
310  rms_norm = 0.0;
311 
312  // Counter for the number of values:
313  unsigned value_count = 0;
314 
315  // Number of solid Data items:
316  unsigned n_data = Solid_data_pt.size();
317 
318  // Loop over all solid Data items:
319  for (unsigned i = 0; i < n_data; i++)
320  {
321  // Number of values stored at this Data item:
322  unsigned n_value = Solid_data_pt[i]->nvalue();
323 
324  // Loop over all values
325  for (unsigned k = 0; k < n_value; k++)
326  {
327  // Change
328  double change =
329  Previous_solid_value[value_count] - Solid_data_pt[i]->value(k);
330 
331  // Max change?
332  if (std::fabs(change) > max_change) max_change = std::fabs(change);
333 
334  // Add square of change relative to previous value
335  rms_change += pow(change, 2);
336 
337  // Add square of previous value
338  rms_norm += pow(Previous_solid_value[value_count], 2);
339 
340  // Increment counter
341  value_count++;
342  }
343  }
344 
345  // Turn into rms:
346  rms_change = sqrt(rms_change / double(value_count));
347  rms_norm = sqrt(rms_norm / double(value_count));
348  }
349 
350 
351  //============================================================================
352  /// Under-relax solid, either by classical relaxation or by Irons & Tuck.
353  //============================================================================
355  {
356  // Irons and Tuck extrapolation/relaxation; an extension of Aitken's method
357  //-------------------------------------------------------------------------
359  {
360  double top = 0.0;
361  double den = 0.0;
362  double crit = 0.0;
363 
364  // Counter for the number of values:
365  unsigned value_count = 0;
366 
367  // Number of solid Data items:
368  unsigned n_data = Solid_data_pt.size();
369 
370  // Loop over all solid Data items:
371  for (unsigned i = 0; i < n_data; i++)
372  {
373  // Number of values stored at this Data item:
374  unsigned n_value = Solid_data_pt[i]->nvalue();
375 
376  // Loop over all values
377  for (unsigned k = 0; k < n_value; k++)
378  {
379  // Prediction from solid solver
380  double new_solid_value = Solid_data_pt[i]->value(k);
381 
382  // Previus value
383  double old_solid_value = Previous_solid_value[value_count];
384 
385  // Change
386  double change = old_solid_value - new_solid_value;
387 
388  // Change of change
389  double del2 = Del_irons_and_tuck[value_count] - change;
390 
391  // Update change
392  Del_irons_and_tuck[value_count] = change;
393 
394  // Update top
395  top += del2 * change;
396 
397  // Update denominator
398  den += del2 * del2;
399 
400  // Update convergence criterion
401  crit += std::fabs(change);
402 
403  // Increment counter
404  value_count++;
405  }
406  }
407 
408  // Update relaxation factor. The if buffers the case in which
409  // we haven't realised that we've converged (so that den=0).
410  // This can happen, e.g. if the convergence assessment is based on the
411  // global residual or during validation. In that case we
412  // obviously don't want any changes to the iterates.
413  if (den != 0.0)
414  {
415  double new_r = R_irons_and_tuck + (R_irons_and_tuck - 1.0) * top / den;
416  R_irons_and_tuck = new_r;
417  }
418  else
419  {
420  R_irons_and_tuck = 0.0;
421  }
422 
423  // Loop over all solid Data items for update
424  value_count = 0;
425  for (unsigned i = 0; i < n_data; i++)
426  {
427  // Number of values stored at this Data item:
428  unsigned n_value = Solid_data_pt[i]->nvalue();
429 
430  // Loop over all values
431  for (unsigned k = 0; k < n_value; k++)
432  {
433  // Compute relaxed/extrapolated value
434  double new_value = Solid_data_pt[i]->value(k) +
435  R_irons_and_tuck * Del_irons_and_tuck[value_count];
436 
437  // Assign
438  Solid_data_pt[i]->set_value(k, new_value);
439 
440  // Increment counter
441  value_count++;
442  }
443  }
444  return;
445  }
446 
447  // Standard relaxation
448  //--------------------
449  else
450  {
451  // No relaxation: Can return immediately
452  if (Omega_relax == 1.0) return;
453 
454  // Counter for the number of values:
455  unsigned value_count = 0;
456 
457  // Number of solid Data items:
458  unsigned n_data = Solid_data_pt.size();
459 
460  // Loop over all solid Data items:
461  for (unsigned i = 0; i < n_data; i++)
462  {
463  // Number of values stored at this Data item:
464  unsigned n_value = Solid_data_pt[i]->nvalue();
465 
466  // Loop over all values
467  for (unsigned k = 0; k < n_value; k++)
468  {
469  // Prediction from solid solver
470  double new_solid_value = Solid_data_pt[i]->value(k);
471 
472  // Previus value
473  double old_solid_value = Previous_solid_value[value_count];
474 
475  // Relax
476  Solid_data_pt[i]->set_value(k,
477  Omega_relax * new_solid_value +
478  (1.0 - Omega_relax) * old_solid_value);
479 
480  // Increment counter
481  value_count++;
482  }
483  }
484  }
485  }
486 
487  //============================================================================
488  /// Pointwise Aitken extrapolation for solid variables
489  //============================================================================
491  {
492  // Counter for the number of values:
493  unsigned value_count = 0;
494 
495  // Number of solid Data items:
496  unsigned n_data = Solid_data_pt.size();
497 
498  // Loop over all solid Data items:
499  for (unsigned i = 0; i < n_data; i++)
500  {
501  // Number of values stored at this Data item:
502  unsigned n_value = Solid_data_pt[i]->nvalue();
503 
504  // Loop over all values
505  for (unsigned k = 0; k < n_value; k++)
506  {
507  // Shorthand
508  double v0 = Pointwise_aitken_solid_value[value_count][0];
509  double v1 = Pointwise_aitken_solid_value[value_count][1];
510  double v2 = Pointwise_aitken_solid_value[value_count][2];
511 
512  double new_value = v2;
513 
514  double max_diff = std::max(std::fabs(v1 - v0), std::fabs(v2 - v1));
515  if (max_diff > 1.0e-10)
516  {
517  new_value = v2 - std::pow((v2 - v1), int(2)) / (v2 - 2.0 * v1 + v0);
518  }
519  Solid_data_pt[i]->set_value(k, new_value);
520 
521  // Increment counter
522  value_count++;
523  }
524  }
525 
526  // Reset the counter for the Aitken convergence check
527  // (setting counter to -1 forces three new genuine
528  // iterates to be computed).
530  }
531 
532 
533  //============================================================================
534  /// Segregated fixed point solver with optional pointwise Aitken acceleration
535  /// on selected solid dofs. Returns PicardConvergenceData object
536  /// that contains the vital stats of the iteration.
537  //============================================================================
539  {
540  // Initialise timer for essential bits of code
541  reset_timer();
542  // Start timer for overall solve
543  clock_t t_total_start = clock();
544 
545  // If convergence is based on max. global residual we may as well
546  // document it...
547  bool get_max_global_residual = Doc_max_global_residual;
548  if (Convergence_criterion ==
550  {
551  get_max_global_residual = true;
552  }
553 
554  // Create object to doc convergence stats
555  PicardConvergenceData conv_data;
556 
557  // Initialise total time for computation of global residuals
558  double cpu_for_global_residual = 0.0;
559 
560  // Update anything that needs updating
562 
563  // Set flags to values that are appropriate if Picard iteration
564  // does not converge with Max_picard iterations
565  bool converged = false;
566  unsigned iter_taken = Max_picard;
567 
568  // This one will be overwritten during the convergence checks
569  double tol_achieved = 0.0;
570 
571  // Store the current values of the solid dofs as reference values
572  // and for the pointwise Aitken extrapolation
574 
575  // Loop over Picard iterations
576  //----------------------------
577  for (unsigned iter = 1; iter <= Max_picard; iter++)
578  {
579  // Calculate the initial residuals
580  if (iter == 1)
581  {
582  // Problem is always non-linear?
583  // Perform any actions before the convergence check
585 
586  double max_res = 0.0;
587  if (get_max_global_residual)
588  {
589  clock_t t_start = clock();
594  // This is now a full problem
596  DoubleVector residual;
597  get_residuals(residual);
598  // Get maximum residuals using our own abscmp function
599  max_res = residual.max();
600  clock_t t_end = clock();
601  cpu_for_global_residual += double(t_end - t_start) / CLOCKS_PER_SEC;
602  }
603 
604  oomph_info << "==================================================\n";
605  oomph_info << "Initial iteration : " << 0 << std::endl;
606  oomph_info << "RMS change : " << 0 << std::endl;
607  oomph_info << "Max. change : " << 0 << std::endl;
608  oomph_info << "RMS norm : " << 0 << std::endl;
609  if (get_max_global_residual)
610  {
611  oomph_info << "Max. global residual : " << max_res
612  << std::endl;
613  }
614  oomph_info << "==================================================\n\n";
615 
616  // Check for convergence, but this only makes sense
617  // for the maximum (rather than relative case)
618  if ((Convergence_criterion ==
620  (max_res < Convergence_tolerance))
621  {
622  oomph_info
623  << "\n\n\n////////////////////////////////////////////////////////"
624  << "\nPicard iteration converged after " << 0 << " steps!"
625  << std::endl
626  << "Convergence was based on max. residual of coupled eqns \n"
627  << "being less than " << Convergence_tolerance << "." << std::endl
628  << "////////////////////////////////////////////////////////\n\n\n"
629  << std::endl;
630 
631  // Converged, so bail out
632  // Number of iterations taken
633  iter_taken = 0;
634 
635  // We have converged (overwrites default of false)
636  conv_data.set_solver_converged();
637 
638  // Break loop using a GO TO! This is the first (and hopefully
639  // the last) one in the whole of oomph-lib. Here it's
640  // it's actually the cleanest way to exit from these
641  // nested control structures
642  goto jump_out_of_picard;
643  }
644  }
645 
646  // Stage 1: Freeze wall and solve single-physics fluid problem
647  //------------------------------------------------------------
648  // Pin the solid dofs, restore the pinned status of the fluid dofs
649  // and re-assign the equation numbers
651  pin_solid_dofs();
654 
655  // Solve the fluid problem for the current wall geometry
656  oomph_info << "\n\nDOING FLUID SOLVE\n\n";
657  // This is a fluid solve at the moment done by a newton solve
659  newton_solve();
660 
661  // Stage 2: Freeze fluid variables and update wall shape
662  //------------------------------------------------------
663 
664  // Now restore the pinned status of the wall and pin the fluid
665  // dofs in anticipation of the wall update:
667  pin_fluid_dofs();
670 
671  // Solve the solid problem for the wall solution
672  oomph_info << "\n\nDOING SOLID SOLVE\n\n";
673  // This is a solid solve, at the moment done by a newton method
675  newton_solve();
676 
677  // Under-relax, either by classical relaxation or by Irons & Tuck
678  // Note that we are under-relaxing before the convergence check
679  // because under-relaxtion may be required to acheieve any
680  // kind of convergence. If the convergence check is on the RELATIVE
681  // change of the residuals, however, then a small under-relaxation
682  // parameter will cause a false convergence. You have been warned!
684 
685  // Stage 3: Convergence check (possibly again after pointwise Aitken
686  //------------------------------------------------------------------
687  // extrapolation)
688  //---------------
689  // Assume that we are not doing Aitken acceleration
691  do
692  {
693  // Perform any actions before the convergence check
695 
696  // Get the change in the solid variables
697  double rms_change;
698  double max_change;
699  double rms_norm;
700  double max_res = 0.0;
701  get_solid_change(rms_change, max_change, rms_norm);
702 
703 
704  // If we are computing the maximum global residual, do so
705  if (get_max_global_residual)
706  {
707  clock_t t_start = clock();
708 
709  // Free all dofs so the residuals of the global equations are computed
714 
715  // We are now treating this as a full solve
717 
718  // Get the residuals
719  DoubleVector residual;
720  get_residuals(residual);
721 
722  // Get maximum residuals, using our own abscmp function
723  max_res = residual.max();
724 
725  clock_t t_end = clock();
726  cpu_for_global_residual += double(t_end - t_start) / CLOCKS_PER_SEC;
727  }
728 
729  oomph_info << "==================================================\n";
730  oomph_info << "Iteration : " << iter << std::endl;
731  oomph_info << "RMS change : " << rms_change
732  << std::endl;
733  oomph_info << "Max. change : " << max_change
734  << std::endl;
735  oomph_info << "RMS norm : " << rms_norm << std::endl;
736  if (get_max_global_residual)
737  {
738  oomph_info << "Max. global residual : " << max_res
739  << std::endl;
740  }
741  oomph_info << "==================================================\n\n";
742 
743  // Check for convergence
744  switch (Convergence_criterion)
745  {
747  tol_achieved = max_change;
748  if (tol_achieved < Convergence_tolerance)
749  {
750  oomph_info
751  << "\n\n\n/////////////////////////////////////////////////////"
752  "///"
753  << "\nPicard iteration converged after " << iter << " steps!"
754  << std::endl
755  << "Convergence was based on absolute change in solid dofs \n"
756  << "being less than " << Convergence_tolerance << "."
757  << std::endl
758  << "////////////////////////////////////////////////////////"
759  "\n\n\n"
760  << std::endl;
761  converged = true;
762  }
763  break;
764 
765 
767  tol_achieved = std::fabs(rms_change / rms_norm);
768  if (tol_achieved < Convergence_tolerance)
769  {
770  oomph_info
771  << "\n\n\n/////////////////////////////////////////////////////"
772  "//"
773  << "\nPicard iteration converged after " << iter << " steps!"
774  << std::endl
775  << "Convergence was based on relative change in solid dofs \n"
776  << "being less than " << Convergence_tolerance << "."
777  << std::endl
778  << "////////////////////////////////////////////////////////"
779  "\n\n\n"
780  << std::endl;
781  converged = true;
782  }
783  break;
784 
786  tol_achieved = max_res;
787  if (tol_achieved < Convergence_tolerance)
788  {
789  oomph_info
790  << "\n\n\n/////////////////////////////////////////////////////"
791  "///"
792  << "\nPicard iteration converged after " << iter << " steps!"
793  << std::endl
794  << "Convergence was based on max. residual of coupled eqns \n"
795  << "being less than " << Convergence_tolerance << "."
796  << std::endl
797  << "////////////////////////////////////////////////////////"
798  "\n\n\n"
799  << std::endl;
800  converged = true;
801  }
802  break;
803  }
804 
805  // If converged bail out
806  if (converged)
807  {
808  // Number of iterations taken
809  iter_taken = iter;
810 
811  // We have converged (overwrites default of false)
812  conv_data.set_solver_converged();
813 
814  // Break loop using a GO TO! This is the first (and hopefully
815  // the last) one in the whole of oomph-lib. Here it's
816  // it's actually the cleanest way to exit from these
817  // nested control structures
818  goto jump_out_of_picard;
819  }
820 
821  // Store the current values of the solid dofs as reference values
822  // and for the pointwise Aitken extrapolation
824 
825  // Correct wall shape by pointwise Aitken extrapolation if possible
826  //-----------------------------------------------------------------
827  // and desired
828  //------------
829  // This is an acceleration method for the (possibly under-relaxed)
830  // sequence of iterates.
832  (iter > Pointwise_aitken_start))
833  {
835 
836  // Repeat the convergence check
838  }
839  else
840  {
841  // Didn't change anything: Don't repeat the convergence check
843  }
844  }
845  // Repeat convergence while we are doing aitken extrapolation
847 
848  } // End of loop over iterations
849 
850 
851  // Jump address for converged or diverged iteration
852  jump_out_of_picard:
853 
854  // Reset everything
859  // This is a full solve again
861 
862  // Do any updates that are required
864 
865  // Number of iterations (either this is still Max_iter from
866  // the initialisation or it's been overwritten on convergence)
867  conv_data.niter() = iter_taken;
868 
869  // Total cpu time
870  clock_t t_total_end = clock();
871  conv_data.cpu_total() =
872  double(t_total_end - t_total_start) / CLOCKS_PER_SEC;
873 
874  // Total essential cpu time (exluding doc etc)
876 
877  // cpu time for check/doc of global residual
878  conv_data.cpu_for_global_residual() = cpu_for_global_residual;
879 
880  // Final tolerance achieved by the iteration
881  conv_data.tol_achieved() = tol_achieved;
882 
883  // Doc non-convergence
884  if (!converged)
885  {
886  switch (Convergence_criterion)
887  {
889  oomph_info
890  << "\n\n\n////////////////////////////////////////////////////////"
891  << "\nPicard iteration did not converge after " << iter_taken
892  << " steps!" << std::endl
893  << "Convergence was based on absolute change in solid dofs \n"
894  << "being less than " << Convergence_tolerance << " " << std::endl
895  << "but we achieved only " << tol_achieved << "." << std::endl
896  << "////////////////////////////////////////////////////////\n\n\n"
897  << std::endl;
898  // Throw an error indicating if we ran out of iterations
899  if (iter_taken == Max_picard)
900  {
901  throw SegregatedSolverError(true);
902  }
903  else
904  {
905  throw SegregatedSolverError(false);
906  }
907  break;
908 
909 
911  oomph_info
912  << "\n\n\n///////////////////////////////////////////////////////"
913  << "\nPicard iteration did not converge after " << iter_taken
914  << " steps!" << std::endl
915  << "Convergence was based on relative change in solid dofs \n"
916  << "being less than " << Convergence_tolerance << " " << std::endl
917  << "but we achieved only " << tol_achieved << "." << std::endl
918  << "////////////////////////////////////////////////////////\n\n\n"
919  << std::endl;
920  // Throw an error indicating if we ran out of iterations
921  if (iter_taken == Max_picard)
922  {
923  throw SegregatedSolverError(true);
924  }
925  else
926  {
927  throw SegregatedSolverError(false);
928  }
929  break;
930 
932  oomph_info
933  << "\n\n\n////////////////////////////////////////////////////////"
934  << "\nPicard iteration did not converge after " << iter_taken
935  << " steps!" << std::endl
936  << "Convergence was based on max. residual of coupled eqns \n"
937  << "being less than " << Convergence_tolerance << " " << std::endl
938  << "but we achieved only " << tol_achieved << "." << std::endl
939 
940  << "////////////////////////////////////////////////////////\n\n\n"
941  << std::endl;
942  // Throw an error indicating if we ran out of iterations
943  if (iter_taken == Max_picard)
944  {
945  throw SegregatedSolverError(true);
946  }
947  else
948  {
949  throw SegregatedSolverError(false);
950  }
951  break;
952  }
953  }
954 
955  return conv_data;
956  }
957 
958 
959  //============================================================================
960  /// Segregated fixed point solver with optional pointwise Aitken acceleration
961  /// on selected solid dofs. Returns PicardConvergenceData object
962  /// that contains the vital stats of the iteration.
963  //============================================================================
965  {
966  // Find out how many timesteppers there are
967  unsigned n_time_steppers = ntime_stepper();
968 
969  // Vector of bools to store the is_steady status of the various
970  // timesteppers when we came in here
971  std::vector<bool> was_steady(n_time_steppers);
972 
973  // Loop over them all and make them (temporarily) static
974  for (unsigned i = 0; i < n_time_steppers; i++)
975  {
976  was_steady[i] = time_stepper_pt(i)->is_steady();
978  }
979 
980  // Create object to doc convergence stats
981  PicardConvergenceData conv_data;
982 
983  // Solve the non-linear problem by the segregated solver
984  try
985  {
986  conv_data = segregated_solve();
987  }
988  // Catch any exceptions thrown in the segregated solver
989  catch (SegregatedSolverError& error)
990  {
991  if (!error.Ran_out_of_iterations)
992  {
993  std::ostringstream error_stream;
994  error_stream << "Error occured in Segregated solver. " << std::endl;
995  throw OomphLibError(
996  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
997  }
998  else
999  {
1000  oomph_info << "Note: Ran out of iterations but continuing anyway"
1001  << std::endl;
1002  }
1003  }
1004 
1005  // Reset the is_steady status of all timesteppers that
1006  // weren't already steady when we came in here and reset their
1007  // weights
1008  for (unsigned i = 0; i < n_time_steppers; i++)
1009  {
1010  if (!was_steady[i])
1011  {
1013  }
1014  }
1015 
1016  // Since we performed a steady solve, the history values
1017  // now have to be set as if we had performed an impulsive start from
1018  // the current solution. This ensures that the time-derivatives
1019  // evaluate to zero even now that the timesteppers have been
1020  // reactivated.
1022 
1023  // Return the convergence data
1024  return conv_data;
1025  }
1026 
1027 
1028  //============================================================================
1029  /// Segregated fixed point solver with optional pointwise Aitken acceleration
1030  /// on selected solid dofs. Returns PicardConvergenceData object
1031  /// that contains the vital stats of the iteration.
1032  //============================================================================
1034  const double& dt)
1035  {
1036  // We shift the values, so shift_values is true
1037  return unsteady_segregated_solve(dt, true);
1038  }
1039 
1040  //============================================================================
1041  /// Segregated fixed point solver with optional pointwise Aitken acceleration
1042  /// on selected solid dofs. Returns PicardConvergenceData object
1043  /// that contains the vital stats of the iteration.
1044  //============================================================================
1046  const double& dt, const bool& shift_values)
1047  {
1048  // Shift the time values and the dts according to the control flag
1049  if (shift_values)
1050  {
1052  }
1053 
1054  // Advance global time and set current value of dt
1055  time_pt()->time() += dt;
1056  time_pt()->dt() = dt;
1057 
1058  // Find out how many timesteppers there are
1059  unsigned n_time_steppers = ntime_stepper();
1060 
1061  // Loop over them all and set the weights
1062  for (unsigned i = 0; i < n_time_steppers; i++)
1063  {
1065  }
1066 
1067  // Now update anything that needs updating before the timestep
1068  // This could be time-dependent boundary conditions, for example.
1070 
1071  // Extrapolate the solid data and then update fluid mesh during unsteady run
1073 
1074  // Create object to doc convergence stats
1075  PicardConvergenceData conv_data;
1076 
1077  try
1078  {
1079  // Solve the non-linear problem for this timestep with Newton's method
1080  conv_data = segregated_solve();
1081  }
1082  // Catch any exceptions thrown in the segregated solver
1083  catch (SegregatedSolverError& error)
1084  {
1085  if (!error.Ran_out_of_iterations)
1086  {
1087  std::ostringstream error_stream;
1088  error_stream << "Error occured in Segregated solver. " << std::endl;
1089  throw OomphLibError(
1090  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1091  }
1092  else
1093  {
1094  oomph_info << "Note: Ran out of iterations but continuing anyway"
1095  << std::endl;
1096  }
1097  }
1098 
1099  // Now update anything that needs updating after the timestep
1101 
1102  return conv_data;
1103  }
1104 
1105 
1106  //============================================================================
1107  /// Setup segregated solver, using the information provided by the user
1108  /// in his/her implementation of the pure virtual function
1109  /// identify_fluid_and_solid_dofs(...).
1110  //============================================================================
1112  const bool& full_setup_of_fluid_and_solid_dofs)
1113  {
1114  // If we are doing a full setup
1115  if (full_setup_of_fluid_and_solid_dofs)
1116  {
1117  // Identify the fluid and solid Data
1118  Vector<Data*> fluid_data_pt;
1119  Vector<Data*> solid_data_pt;
1121  fluid_data_pt, solid_data_pt, Fluid_mesh_pt, Solid_mesh_pt);
1122 
1123  if (Fluid_mesh_pt == 0)
1124  {
1125  oomph_info
1126  << std::endl
1127  << std::endl
1128  << "Warning: Your implementation of the pure virtual\n"
1129  << " function identify_fluid_and_solid_dofs(...)\n"
1130  << " returned a NULL pointer for Fluid_mesh_pt.\n"
1131  << " --> The fluid elements will remain activated\n"
1132  << " during the solid solve. This is inefficient!\n"
1133  << " You should combine all fluid elements into a combined\n"
1134  << " mesh and specify this mesh in your\n"
1135  << " implementation of \n\n"
1136  << " "
1137  "SegregatableFSIProblem::identify_fluid_and_solid_dofs(...)"
1138  << std::endl
1139  << std::endl;
1140  }
1141 
1142 
1143  if (Solid_mesh_pt == 0)
1144  {
1145  oomph_info
1146  << std::endl
1147  << std::endl
1148  << "Warning: Your implementation of the pure virtual\n"
1149  << " function identify_fluid_and_solid_dofs(...)\n"
1150  << " returned a NULL pointer for Solid_mesh_pt.\n"
1151  << " --> The solid elements will remain activated\n"
1152  << " during the fluid solve. This is inefficient!\n"
1153  << " You should combine all solid elements into a combined\n"
1154  << " mesh and specify this mesh in your\n"
1155  << " implementation of \n\n"
1156  << " "
1157  "SegregatableFSIProblem::identify_fluid_and_solid_dofs(...)"
1158  << std::endl
1159  << std::endl;
1160  }
1161 
1162  // Back up the pointers to the submeshes in the original problem
1163  // so we can restore the problem when we're done
1164  unsigned orig_n_sub_mesh = nsub_mesh();
1165  Orig_sub_mesh_pt.resize(orig_n_sub_mesh);
1166  for (unsigned i = 0; i < orig_n_sub_mesh; i++)
1167  {
1169  }
1170 
1171  // Fluid
1172  //------
1173 
1174  // Reset
1175  Fluid_data_pt.clear();
1176  Fluid_value_is_pinned.clear();
1177 
1178  // Make space for fluid data items
1179  unsigned n_fluid_data = fluid_data_pt.size();
1180  Fluid_data_pt.resize(n_fluid_data);
1181  Fluid_value_is_pinned.resize(n_fluid_data);
1182 
1183  // Counter for number of fluid values
1184  unsigned n_fluid_values = 0;
1185 
1186  // Loop over fluid data
1187  for (unsigned i = 0; i < n_fluid_data; i++)
1188  {
1189  // Copy data
1190  Fluid_data_pt[i] = fluid_data_pt[i];
1191 
1192  // Number of values stored at this Data item:
1193  unsigned n_value = fluid_data_pt[i]->nvalue();
1194  Fluid_value_is_pinned[i].resize(n_value);
1195 
1196  // Copy pinned status for all values
1197  for (unsigned k = 0; k < n_value; k++)
1198  {
1199  Fluid_value_is_pinned[i][k] = fluid_data_pt[i]->is_pinned(k);
1200 
1201  // Increment counter for number of fluid values
1202  n_fluid_values++;
1203  }
1204  }
1205 
1206 
1207  // Solid
1208  //------
1209 
1210  // Reset
1211  Solid_data_pt.clear();
1212  Solid_value_is_pinned.clear();
1213  Previous_solid_value.clear();
1215  Del_irons_and_tuck.clear();
1216 
1217 
1218  unsigned n_solid_data = solid_data_pt.size();
1219 
1220  // Make space for solid data items
1221  unsigned nsolid_data = solid_data_pt.size();
1222  Solid_data_pt.resize(nsolid_data);
1223  Solid_value_is_pinned.resize(nsolid_data);
1224 
1225  // Counter for number of solid values
1226  unsigned n_solid_values = 0;
1227 
1228  // Loop over solid data
1229  for (unsigned i = 0; i < n_solid_data; i++)
1230  {
1231  // Copy data
1232  Solid_data_pt[i] = solid_data_pt[i];
1233 
1234  // Number of values stored at this Data item:
1235  unsigned n_value = solid_data_pt[i]->nvalue();
1236  Solid_value_is_pinned[i].resize(n_value);
1237 
1238  // Copy pinned status for all values
1239  for (unsigned k = 0; k < n_value; k++)
1240  {
1241  Solid_value_is_pinned[i][k] = solid_data_pt[i]->is_pinned(k);
1242 
1243  // Increment counter for number of solid values
1244  n_solid_values++;
1245  }
1246  }
1247 
1248  // Make space for previous solid values
1249  Previous_solid_value.resize(n_solid_values);
1250 
1251  // Allocate storage and initialise Irons and Tuck extrapolation
1253  {
1254  Del_irons_and_tuck.resize(n_solid_values);
1255  }
1256 
1257  // Make space for pointwise Aitken extrapolation
1259  {
1260  Pointwise_aitken_solid_value.resize(n_solid_values);
1261  for (unsigned i = 0; i < n_solid_values; i++)
1262  {
1263  Pointwise_aitken_solid_value[i].resize(3);
1264  }
1265  }
1266 
1267  } // End of full setup
1268 
1269 
1270  // Initialise Irons and Tuck relaxation factor
1271  R_irons_and_tuck = 1.0 - Omega_relax;
1272 
1273  // Initialise Irons and Tuck delta values
1274  unsigned n_del = Del_irons_and_tuck.size();
1275  for (unsigned i = 0; i < n_del; i++)
1276  {
1277  Del_irons_and_tuck[i] = 1.0e20;
1278  }
1279 
1280  // Initialise counter for the number of pointwise Aitken values stored
1282  }
1283 } // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
double max() const
returns the maximum coefficient
virtual void node_update(const bool &update_all_solid_nodes=false)
Update nodal positions in response to changes in the domain shape. Uses the FiniteElement::get_x(....
Definition: mesh.cc:287
An OomphLibError object which should be thrown when an run-time error is encountered....
Object that collates convergence data of Picard iteration.
void set_solver_converged()
Set the flag to indicate that the solver has converged.
unsigned & niter()
Number of iterations performed.
double & cpu_for_global_residual()
CPU time for computation of global residual vectors Note: This time is contained in Total_CPU and is ...
double & tol_achieved()
Final tolerance achieved by the iteration.
double & essential_cpu_total()
Total essential CPU time for segregated solve (excluding any actions that merely doc the progress of ...
double & cpu_total()
Total CPU time for segregated solve.
virtual void actions_after_implicit_timestep()
Actions that should be performed after each implicit time step. This is needed when one wants to solv...
Definition: problem.h:1070
unsigned add_sub_mesh(Mesh *const &mesh_pt)
Add a submesh to the problem and return its number, i, by which it can be accessed via mesh_pt(i).
Definition: problem.h:1330
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Assign all equation numbers for problem: Deals with global data (= data that isn't attached to any el...
Definition: problem.cc:2075
void assign_initial_values_impulsive()
Initialise data and nodal positions to simulate impulsive start from initial configuration/solution.
Definition: problem.cc:11692
void flush_sub_meshes()
Flush the problem's collection of sub-meshes. Must be followed by call to rebuild_global_mesh().
Definition: problem.h:1339
Time *& time_pt()
Return a pointer to the global time object.
Definition: problem.h:1504
virtual void actions_before_implicit_timestep()
Actions that should be performed before each implicit time step. This is needed when one wants to sol...
Definition: problem.h:1064
void newton_solve()
Use Newton method to solve the problem.
Definition: problem.cc:8976
virtual void get_residuals(DoubleVector &residuals)
Return the fully-assembled residuals Vector for the problem: Virtual so it can be overloaded in for m...
Definition: problem.cc:3800
virtual void shift_time_values()
Shift all values along to prepare for next timestep.
Definition: problem.cc:11827
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
void rebuild_global_mesh()
If one of the submeshes has changed (e.g. by mesh adaptation) we need to update the global mesh....
Definition: problem.cc:1619
unsigned nsub_mesh() const
Return number of submeshes.
Definition: problem.h:1323
unsigned ntime_stepper() const
Return the number of time steppers.
Definition: problem.h:1684
TimeStepper *& time_stepper_pt()
Access function for the pointer to the first (presumably only) timestepper.
Definition: problem.h:1524
bool Recheck_convergence_after_pointwise_aitken
Have we just done a pointwise Aitken step.
void extrapolate_solid_data()
Extrapolate solid data and update fluid mesh during unsteady run.
Vector< std::vector< bool > > Solid_value_is_pinned
Vector of vectors that store the pinned status of solid Data values.
void setup_segregated_solver(const bool &full_setup_of_fluid_and_solid_dofs=true)
Setup the segregated solver: Backup the pinned status of the fluid and solid dofs and allocate the in...
void use_only_fluid_elements()
Only include fluid elements in the Problem's mesh. This is called before the segregated fluid solve....
Vector< Mesh * > Orig_sub_mesh_pt
Backup for the pointers to the submeshes in the original problem.
bool Use_pointwise_aitken
Use pointwise Aitken extrapolation?
Vector< Data * > Fluid_data_pt
Vector storing the Data objects associated with the fluid problem: Tyically the nodal and internal da...
bool Use_irons_and_tuck_extrapolation
Boolean flag to indicate use of Irons and Tuck's extrapolation for solid values.
virtual void actions_after_segregated_solve()
This function is called once at the end of each segregated solve.
void use_only_solid_elements()
Only include solid elements in the Problem's mesh. This is called before the segregated solid solve....
Mesh * Solid_mesh_pt
Mesh containing only solid elements – the elements in this mesh will be excluded from the assembly pr...
unsigned Pointwise_aitken_start
Start pointwise Aitken extrpolation after specified number of Picard iterations.
void rebuild_monolithic_mesh()
Rebuild global mesh for monolithic discretisation.
virtual void identify_fluid_and_solid_dofs(Vector< Data * > &fluid_data_pt, Vector< Data * > &solid_data_pt, Mesh *&fluid_mesh_pt, Mesh *&solid_mesh_pt)=0
Identify the fluid and solid Data. This is a pure virtual function that MUST be implemented for every...
void restore_solid_dofs()
Restore pinned status of solid dofs.
Vector< double > Previous_solid_value
Vector storing the previous solid values – used for convergence check.
double R_irons_and_tuck
Irons and Tuck relaxation factor.
PicardConvergenceData segregated_solve()
Segregated solver. Peform a segregated step from the present state of the system. Returns PicardConve...
Vector< double > Del_irons_and_tuck
Vector of changes in Irons and Tuck under-relaxation.
PicardConvergenceData steady_segregated_solve()
Steady version of segregated solver. Makes all timesteppers steady before solving....
int Convergence_criterion
Convergence criterion (enumerated flag)
void pointwise_aitken_extrapolate()
Do pointwise Aitken extrapolation for solid.
void under_relax_solid()
Under-relax the most recently computed solid variables, either by classical relaxation or by Irons & ...
Vector< Vector< double > > Pointwise_aitken_solid_value
Vector of Vectors containing up to three previous iterates for the solid dofs; used for pointwise Ait...
unsigned Max_picard
Max. number of Picard iterations.
Vector< std::vector< bool > > Fluid_value_is_pinned
Vector of vectors that store the pinned status of fluid Data values.
bool Doc_max_global_residual
Doc maximum global residual during iteration? (default: false)
PicardConvergenceData unsteady_segregated_solve(const double &dt)
Unsteady segregated solver, advance time by dt and solve by the segregated solver....
void restore_fluid_dofs()
Restore pinned status of fluid dofs.
virtual void actions_before_segregated_convergence_check()
This function is to be filled with actions that take place before the check for convergence of the en...
Vector< Data * > Solid_data_pt
Vector storing the Data objects associated with the solid problem: Typically the positional data of s...
void get_solid_change(double &rms_change, double &max_change, double &rms_norm)
Get rms of change in the solid dofs; the max. change of the solid dofs and the rms norm of the solid ...
int Pointwise_aitken_counter
Number of Aitken histories available (int because after extrapolation it's re-initialised to -1 to fo...
double Convergence_tolerance
Convergence tolerance for Picard iteration.
Mesh * Fluid_mesh_pt
Mesh containing only fluid elements – the elements in this Mesh will be excluded from the assembly pr...
double t_spent_on_actual_solve()
Total elapsed time since start of solve.
double Omega_relax
Under-relaxation parameter. (1.0: no under-relaxation; 0.0: Freeze wall shape)
virtual void actions_before_segregated_solve()
This function is called once at the start of each segregated solve.
int Solve_type
Solve that is taking place (enumerated flag)
void store_solid_dofs()
Store the current solid values as reference values for future convergence check. Also add another ent...
////////////////////////////////////////////////////////////////// //////////////////////////////////...
bool Ran_out_of_iterations
Boolean indicating if the error occured because we ran out of iterations.
virtual void set_weights()=0
Function to set the weights for present timestep (don't need to pass present timestep or previous tim...
void make_steady()
Function to make the time stepper temporarily steady. This is trivially achieved by setting all the w...
Definition: timesteppers.h:374
virtual void undo_make_steady()
Reset the is_steady status of a specific TimeStepper to its default and re-assign the weights.
Definition: timesteppers.h:482
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:389
double & dt(const unsigned &t=0)
Return the value of the t-th stored timestep (t=0: present; t>0: previous).
Definition: timesteppers.h:136
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...