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