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-2025 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:
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
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 {
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 {
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
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 =
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
339
340 // Increment counter
341 value_count++;
342 }
343 }
344
345 // Turn into rms:
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
384
385 // Change
387
388 // Change of change
390
391 // Update 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;
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
474
475 // Relax
476 Solid_data_pt[i]->set_value(k,
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
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
544
545 // If convergence is based on max. global residual we may as well
546 // document it...
550 {
552 }
553
554 // Create object to doc convergence stats
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;
588 {
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;
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)
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
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;
702
703
704 // If we are computing the maximum global residual, do so
706 {
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;
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
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
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.
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
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
871 conv_data.cpu_total() =
873
874 // Total essential cpu time (exluding doc etc)
875 conv_data.essential_cpu_total() = t_spent_on_actual_solve();
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 {
902 "Error occured in Segregated solver. \n",
905 }
906 else
907 {
909 "Error occured in Segregated solver. \n",
912 }
913 break;
914
915
918 << "\n\n\n///////////////////////////////////////////////////////"
919 << "\nPicard iteration did not converge after " << iter_taken
920 << " steps!" << std::endl
921 << "Convergence was based on relative change in solid dofs \n"
922 << "being less than " << Convergence_tolerance << " " << std::endl
923 << "but we achieved only " << tol_achieved << "." << std::endl
924 << "////////////////////////////////////////////////////////\n\n\n"
925 << std::endl;
926 // Throw an error indicating if we ran out of iterations
927 if (iter_taken == Max_picard)
928 {
930 "Error occured in Segregated solver. \n",
933 }
934 else
935 {
937 "Error occured in Segregated solver. \n",
940 }
941 break;
942
945 << "\n\n\n////////////////////////////////////////////////////////"
946 << "\nPicard iteration did not converge after " << iter_taken
947 << " steps!" << std::endl
948 << "Convergence was based on max. residual of coupled eqns \n"
949 << "being less than " << Convergence_tolerance << " " << std::endl
950 << "but we achieved only " << tol_achieved << "." << std::endl
951
952 << "////////////////////////////////////////////////////////\n\n\n"
953 << std::endl;
954 // Throw an error indicating if we ran out of iterations
955 if (iter_taken == Max_picard)
956 {
958 "Error occured in Segregated solver. \n",
961 }
962 else
963 {
965 "Error occured in Segregated solver. \n",
968 }
969 break;
970 }
971 }
972
973 return conv_data;
974 }
975
976
977 //============================================================================
978 /// Segregated fixed point solver with optional pointwise Aitken acceleration
979 /// on selected solid dofs. Returns PicardConvergenceData object
980 /// that contains the vital stats of the iteration.
981 //============================================================================
983 {
984 // Find out how many timesteppers there are
985 unsigned n_time_steppers = ntime_stepper();
986
987 // Vector of bools to store the is_steady status of the various
988 // timesteppers when we came in here
989 std::vector<bool> was_steady(n_time_steppers);
990
991 // Loop over them all and make them (temporarily) static
992 for (unsigned i = 0; i < n_time_steppers; i++)
993 {
996 }
997
998 // Create object to doc convergence stats
1000
1001 // Solve the non-linear problem by the segregated solver
1002 try
1003 {
1005 }
1006 // Catch any exceptions thrown in the segregated solver
1008 {
1009 // Continue, but output note
1010 oomph_info << "Note: Ran out of iterations but continuing anyway"
1011 << std::endl;
1012 }
1013
1014 // Reset the is_steady status of all timesteppers that
1015 // weren't already steady when we came in here and reset their
1016 // weights
1017 for (unsigned i = 0; i < n_time_steppers; i++)
1018 {
1019 if (!was_steady[i])
1020 {
1022 }
1023 }
1024
1025 // Since we performed a steady solve, the history values
1026 // now have to be set as if we had performed an impulsive start from
1027 // the current solution. This ensures that the time-derivatives
1028 // evaluate to zero even now that the timesteppers have been
1029 // reactivated.
1031
1032 // Return the convergence data
1033 return conv_data;
1034 }
1035
1036
1037 //============================================================================
1038 /// Segregated fixed point solver with optional pointwise Aitken acceleration
1039 /// on selected solid dofs. Returns PicardConvergenceData object
1040 /// that contains the vital stats of the iteration.
1041 //============================================================================
1043 const double& dt)
1044 {
1045 // We shift the values, so shift_values is true
1046 return unsteady_segregated_solve(dt, true);
1047 }
1048
1049 //============================================================================
1050 /// Segregated fixed point solver with optional pointwise Aitken acceleration
1051 /// on selected solid dofs. Returns PicardConvergenceData object
1052 /// that contains the vital stats of the iteration.
1053 //============================================================================
1055 const double& dt, const bool& shift_values)
1056 {
1057 // Shift the time values and the dts according to the control flag
1058 if (shift_values)
1059 {
1061 }
1062
1063 // Advance global time and set current value of dt
1064 time_pt()->time() += dt;
1065 time_pt()->dt() = dt;
1066
1067 // Find out how many timesteppers there are
1068 unsigned n_time_steppers = ntime_stepper();
1069
1070 // Loop over them all and set the weights
1071 for (unsigned i = 0; i < n_time_steppers; i++)
1072 {
1074 }
1075
1076 // Now update anything that needs updating before the timestep
1077 // This could be time-dependent boundary conditions, for example.
1079
1080 // Extrapolate the solid data and then update fluid mesh during unsteady run
1082
1083 // Create object to doc convergence stats
1085
1086 try
1087 {
1088 // Solve the non-linear problem for this timestep with Newton's method
1090 }
1091 // Catch any exceptions thrown in the segregated solver
1093 {
1094 // Continue, but output note
1095 oomph_info << "Note: Ran out of iterations but continuing anyway"
1096 << std::endl;
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 //============================================================================
1113 {
1114 // If we are doing a full setup
1116 {
1117 // Identify the fluid and solid Data
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();
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();
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
1191
1192 // Number of values stored at this Data item:
1193 unsigned n_value = fluid_data_pt[i]->nvalue();
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
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);
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
1233
1234 // Number of values stored at this Data item:
1235 unsigned n_value = solid_data_pt[i]->nvalue();
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
1245 }
1246 }
1247
1248 // Make space for previous solid values
1250
1251 // Allocate storage and initialise Irons and Tuck extrapolation
1253 {
1255 }
1256
1257 // Make space for pointwise Aitken extrapolation
1259 {
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....
double max() const
returns the maximum coefficient
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
Definition elements.cc:4320
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
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
Object that collates convergence data of Picard iteration.
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:1090
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:1350
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:2085
void assign_initial_values_impulsive()
Initialise data and nodal positions to simulate impulsive start from initial configuration/solution.
Definition problem.cc:11575
void flush_sub_meshes()
Flush the problem's collection of sub-meshes. Must be followed by call to rebuild_global_mesh().
Definition problem.h:1359
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:1084
void newton_solve()
Use Newton method to solve the problem.
Definition problem.cc:8816
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:3810
TimeStepper *& time_stepper_pt()
Access function for the pointer to the first (presumably only) timestepper.
Definition problem.h:1544
virtual void shift_time_values()
Shift all values along to prepare for next timestep.
Definition problem.cc:11710
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:1629
unsigned nsub_mesh() const
Return number of submeshes.
Definition problem.h:1343
unsigned ntime_stepper() const
Return the number of time steppers.
Definition problem.h:1710
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition problem.h:1300
Time *& time_pt()
Return a pointer to the global time object.
Definition problem.h:1524
////////////////////////////////////////////////////////////////// //////////////////////////////////...
bool Recheck_convergence_after_pointwise_aitken
Have we just done a pointwise Aitken step.
void extrapolate_solid_data()
Extrapolate solid data and update fluid mesh during unsteady run.
Vector< std::vector< bool > > Solid_value_is_pinned
Vector of vectors that store the pinned status of solid Data values.
void setup_segregated_solver(const bool &full_setup_of_fluid_and_solid_dofs=true)
Setup the segregated solver: Backup the pinned status of the fluid and solid dofs and allocate the in...
void use_only_fluid_elements()
Only include fluid elements in the Problem's mesh. This is called before the segregated fluid solve....
Vector< Mesh * > Orig_sub_mesh_pt
Backup for the pointers to the submeshes in the original problem.
bool Use_pointwise_aitken
Use pointwise Aitken extrapolation?
Vector< Data * > Fluid_data_pt
Vector storing the Data objects associated with the fluid problem: Tyically the nodal and internal da...
bool Use_irons_and_tuck_extrapolation
Boolean flag to indicate use of Irons and Tuck's extrapolation for solid values.
virtual void actions_after_segregated_solve()
This function is called once at the end of each segregated solve.
void use_only_solid_elements()
Only include solid elements in the Problem's mesh. This is called before the segregated solid solve....
Mesh * Solid_mesh_pt
Mesh containing only solid elements – the elements in this mesh will be excluded from the assembly pr...
unsigned Pointwise_aitken_start
Start pointwise Aitken extrpolation after specified number of Picard iterations.
void rebuild_monolithic_mesh()
Rebuild global mesh for monolithic discretisation.
virtual void identify_fluid_and_solid_dofs(Vector< Data * > &fluid_data_pt, Vector< Data * > &solid_data_pt, Mesh *&fluid_mesh_pt, Mesh *&solid_mesh_pt)=0
Identify the fluid and solid Data. This is a pure virtual function that MUST be implemented for every...
void restore_solid_dofs()
Restore pinned status of solid dofs.
Vector< double > Previous_solid_value
Vector storing the previous solid values – used for convergence check.
double R_irons_and_tuck
Irons and Tuck relaxation factor.
PicardConvergenceData segregated_solve()
Segregated solver. Peform a segregated step from the present state of the system. Returns PicardConve...
Vector< double > Del_irons_and_tuck
Vector of changes in Irons and Tuck under-relaxation.
PicardConvergenceData steady_segregated_solve()
Steady version of segregated solver. Makes all timesteppers steady before solving....
int Convergence_criterion
Convergence criterion (enumerated flag)
void pointwise_aitken_extrapolate()
Do pointwise Aitken extrapolation for solid.
void under_relax_solid()
Under-relax the most recently computed solid variables, either by classical relaxation or by Irons & ...
Vector< Vector< double > > Pointwise_aitken_solid_value
Vector of Vectors containing up to three previous iterates for the solid dofs; used for pointwise Ait...
unsigned Max_picard
Max. number of Picard iterations.
Vector< std::vector< bool > > Fluid_value_is_pinned
Vector of vectors that store the pinned status of fluid Data values.
bool Doc_max_global_residual
Doc maximum global residual during iteration? (default: false)
PicardConvergenceData unsteady_segregated_solve(const double &dt)
Unsteady segregated solver, advance time by dt and solve by the segregated solver....
void restore_fluid_dofs()
Restore pinned status of fluid dofs.
virtual void actions_before_segregated_convergence_check()
This function is to be filled with actions that take place before the check for convergence of the en...
Vector< Data * > Solid_data_pt
Vector storing the Data objects associated with the solid problem: Typically the positional data of s...
void get_solid_change(double &rms_change, double &max_change, double &rms_norm)
Get rms of change in the solid dofs; the max. change of the solid dofs and the rms norm of the solid ...
int Pointwise_aitken_counter
Number of Aitken histories available (int because after extrapolation it's re-initialised to -1 to fo...
double Convergence_tolerance
Convergence tolerance for Picard iteration.
Mesh * Fluid_mesh_pt
Mesh containing only fluid elements – the elements in this Mesh will be excluded from the assembly pr...
double t_spent_on_actual_solve()
Total elapsed time since start of solve.
double Omega_relax
Under-relaxation parameter. (1.0: no under-relaxation; 0.0: Freeze wall shape)
virtual void actions_before_segregated_solve()
This function is called once at the start of each segregated solve.
int Solve_type
Solve that is taking place (enumerated flag)
void store_solid_dofs()
Store the current solid values as reference values for future convergence check. Also add another ent...
A class to handle errors in the Segregated solver.
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
virtual void set_weights()=0
Function to set the weights for present timestep (don't need to pass present timestep or previous tim...
virtual unsigned nprev_values() const =0
Number of previous values available: 0 for static, 1 for BDF<1>,...
void make_steady()
Function to make the time stepper temporarily steady. This is trivially achieved by setting all the w...
virtual void undo_make_steady()
Reset the is_steady status of a specific TimeStepper to its default and re-assign the weights.
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
double & time()
Return the current value of the continuous time.
double & dt(const unsigned &t=0)
Return the value of the t-th stored timestep (t=0: present; t>0: previous).
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...