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-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26
27#include <algorithm>
28
30#include <algorithm>
31
32namespace 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 {
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) +
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;
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)
620 (max_res < Convergence_tolerance))
621 {
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
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:
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 {
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 {
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 {
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 {
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
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
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 {
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 {
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 {
1264 }
1265 }
1266
1267 } // End of full setup
1268
1269
1270 // Initialise Irons and Tuck relaxation factor
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.
double & cpu_total()
Total CPU time for segregated solve.
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 ...
unsigned & niter()
Number of iterations performed.
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:1989
void assign_initial_values_impulsive()
Initialise data and nodal positions to simulate impulsive start from initial configuration/solution.
Definition: problem.cc:11499
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
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:8783
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:3714
TimeStepper *& time_stepper_pt()
Access function for the pointer to the first (presumably only) timestepper.
Definition: problem.h:1524
virtual void shift_time_values()
Shift all values along to prepare for next timestep.
Definition: problem.cc:11634
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:1533
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:1680
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
Time *& time_pt()
Return a pointer to the global time object.
Definition: problem.h:1504
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 & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
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
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...