segregated_fsi_solver.h

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
#ifndef OOMPH_SEGREGATED_FSI_SOLVER
#define OOMPH_SEGREGATED_FSI_SOLVER


#include "../generic/problem.h"
#include "../generic/geom_objects.h"
#include "../generic/mesh.h"

namespace oomph
{


//===============================================================
/// Object that collates convergence data of Picard iteration
//===============================================================
class PicardConvergenceData
 {
 
   public:

  /// Constructor initialises all data
  PicardConvergenceData() : Niter(0), CPU_total(0.0), 
   Essential_cpu_total(0.0), CPU_for_global_residual(0.0), 
   Tol_achieved(0.0), Has_converged(false) 
   {}

  /// Empty destructor 
  ~PicardConvergenceData(){};
  
  /// Number of iterations performed
  unsigned& niter(){return Niter;}
  
  /// Total CPU time for segregated solve
  double& cpu_total(){return CPU_total;}
  
  /// \short Total essential CPU time for segregated solve
  /// (excluding any actions that merely doc the progress 
  /// of the iteration, etc.)
  double& essential_cpu_total(){return Essential_cpu_total;}

  /// \short CPU time for computation of global residual vectors
  /// Note: This time is contained in Total_CPU and is
  /// only used if convergence is based on the residual
  /// of the fully coupled system.
  double& cpu_for_global_residual(){return CPU_for_global_residual;}
  
  /// Final tolerance achieved by the iteration
  double& tol_achieved(){return Tol_achieved;}
  
  /// Flag to indicate if the solver has converged
  bool has_converged() const {return Has_converged;}
  
  /// Set the flag to indicate that the solver has converged
  void set_solver_converged() {Has_converged=true;}

  /// Set the flag to indicate that the solver has not converged
  void set_solver_not_converged() {Has_converged=false;}

   private:
  
  /// Number of iterations performed
  unsigned Niter;
  
  /// Total CPU time for segregated solve
  double CPU_total;
  
  /// \short Total essential CPU time for segregated solve
  /// (excluding any actions that merely doc the progress 
  /// of the iteration, etc.)
  double Essential_cpu_total;
  
  /// \short  CPU time for computation of global residual vectors
  /// Note: This time is contained in Total_CPU and is
  /// only used if convergence is based on the residual
  /// of the fully coupled system
 double CPU_for_global_residual;

 /// Final tolerance achieved by the iteration
 double Tol_achieved;

 /// Flag to indicate if the solver has converged
 bool Has_converged;

 };


/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////



//=======================================================================
/// A class to handle errors in the Segregated solver
//=======================================================================
class SegregatedSolverError
{
  public:

 /// Default constructor, does nothing
 SegregatedSolverError(const bool& ran_out_of_iterations=false)
  {
   Ran_out_of_iterations=ran_out_of_iterations;
  }

 /// Boolean indicating if the error occured because
 /// we ran out of iterations
 bool Ran_out_of_iterations;

};



/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////



//===============================================================
/// Base class for problems that can be solved by segregated 
/// FSI solver
//===============================================================
 class SegregatableFSIProblem : public virtual Problem
  {
    protected:

   /// \short This function is called once at the start of each
   /// segregated solve.
   virtual void actions_before_segregated_solve() {}

   /// \short This function is called once at the end of each
   /// segregated solve.
   virtual void actions_after_segregated_solve() {}

   /// \short This function is to be filled with actions that take place
   /// before the check for convergence of the entire segregated solve
   virtual void actions_before_segregated_convergence_check() {}

    public:

   /// \short Constructor. Set default values for solver parameters:
   /// - Don't use pointwise Aitken extrapolation but if it's used at 
   ///   all, it's used immediately.
   /// - No under-relaxation at all (neither classical nor Irons&Tuck)
   /// - Base convergence on max. residual of coupled system of eqns
   /// - Convergence tolerance = tolerance for Newton solver
   ///   defined in Problem base class
   /// - Max. 50 Picard iterations
   SegregatableFSIProblem() 
    {
     // Use pointwise Aitken extrapolation?
     Use_pointwise_aitken=false;

     // Default: No under-relaxation
     Omega_relax=1.0;
 
     // Don't use of Irons and Tuck's extrapolation for solid values
     Use_irons_and_tuck_extrapolation=false;

     // Start using pointwise Aitken immediately
     Pointwise_aitken_start=0;

     // By default we don't recheck convergence
     Recheck_convergence_after_pointwise_aitken=false;

     // Default solve type is full solve
     Solve_type = Full_solve;

     // Convergence criterion
     Convergence_criterion=Assess_convergence_based_on_max_global_residual;

     // Convergence tolerance (as in global Newton solver)
     Convergence_tolerance=Problem::Newton_solver_tolerance;

     // Doc max. global residual during iteration?
     Doc_max_global_residual=false;

     // Max. number of Picard iterations
     Max_picard=50;

     // Pointer to Mesh containing only fluid elements -- the elements in this
     // Mesh will be excluded from the assembly process when
     // the solid problem is solved
     Fluid_mesh_pt=0;

     // Pointer to Mesh containing only solid elements -- the elements in this
     // mesh will be excluded from the assembly process when
     // the fluid problem is solved
     Solid_mesh_pt=0;

     // Initialise timer that allows doc of iteration/cpu time
     T_ref = clock();
     T_spent_on_actual_solve=0.0;

     /// \short boolean flag to indicate if timer has been halted
     Timer_has_been_halted=false;
    }

   /// Empty virtual destructor
   virtual ~SegregatableFSIProblem(){}

   /// \short Identify the fluid and solid Data. This is a pure virtual
   /// function that MUST be implemented for every specific problem that
   /// is to be solved by the segregated solver.
   /// The two mesh pointers identify meshes that contain
   /// elements and nodes used during the fluid 
   /// or solid solves respectively. Elements that feature during
   /// both phases of the segretated solution must be included in both.
   /// These pointers may be set to NULL. In this case, all elements
   /// in the global mesh (set up during the monolithic discretisation 
   /// of the problem) contribute to the global Jacobian matrix
   /// and the residual vector, even if their contributions only contain
   /// zero entries. This can be costly, though the code will
   /// still compute the correct results.
   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;

   /// \short Setup the segregated solver: Backup the pinned status of
   /// the fluid and solid dofs and allocate the internal storage
   /// based on the input provided by identify_fluid_and_solid_dofs(...)
   /// In addition, reset storage associated with convergence acceleration
   /// techniques.
   /// If the problem and degrees of freedom has not changed between
   /// calls to the solver then it is wasteful to call 
   /// identify_fluid_and_solid_dofs(...) again and again. If the optional
   /// boolean flag is set to false then the storage for convergence
   /// acceleration techniques is reset, but the fluid and solid dofs
   /// are not altered.
   void setup_segregated_solver(
    const bool &full_setup_of_fluid_and_solid_dofs=true);
   
   /// \short Segregated solver. Peform a segregated step from
   /// the present state of the system.
   /// Returns PicardConvergenceData object that contains the vital
   /// stats of the iteration
   PicardConvergenceData segregated_solve();
   
   /// \short Steady version of segregated solver. Makes all
   /// timesteppers steady before solving.
   /// Returns PicardConvergenceData object that contains the
   /// vital stats of the iteration. 
   PicardConvergenceData steady_segregated_solve();
   

   /// \short 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. 
   PicardConvergenceData unsteady_segregated_solve(const double& dt);


   /// \short Unsteady segregated solver. Advance time by dt and solve
   /// the system by a segregated method. The boolean flag is used to
   /// control whether the time values should be shifted. If it is true 
   /// the current data values will be shifted (stored as previous 
   /// timesteps) before the solution step.
   /// Returns PicardConvergenceData object that contains the
   /// vital stats of the iteration. 
   PicardConvergenceData unsteady_segregated_solve(const double& dt,
                                                   const bool &shift_values);
   

   /// \short Assess convergence based on max. residual of coupled system of 
   /// eqns. The argument specifies the convergence tolerance.
   void assess_convergence_based_on_max_global_residual(
    const double& tol)
    {
     Convergence_criterion=Assess_convergence_based_on_max_global_residual;
     Convergence_tolerance=tol;
    }

   /// \short Assess convergence based on max. residuals of coupled
   /// system of eqns. This interface has no argument
   /// and the default convergence tolerance
   /// for the Newton solver, Problem::Newton_solver_tolerance is used.
   void assess_convergence_based_on_max_global_residual()
    {
     assess_convergence_based_on_max_global_residual(
      Problem::Newton_solver_tolerance);
    }

   /// \short Assess convergence based on max. absolute change of solid
   /// dofs. The argument specifies the convergence tolerance.
   void assess_convergence_based_on_absolute_solid_change(
    const double& tol)
    {
     Convergence_criterion=Assess_convergence_based_on_absolute_solid_change;
     Convergence_tolerance=tol;
    }

   /// \short Assess convergence based on max. absolute change of solid
   /// dofs. This interface has no argument and the default 
   /// convergence tolerance
   /// for the Newton solver, Problem::Newton_solver_tolerance is used.
   void assess_convergence_based_on_absolute_solid_change()
    {
     assess_convergence_based_on_absolute_solid_change(
      Problem::Newton_solver_tolerance);
    }
   
   /// \short Assess convergence based on max. relative change of solid
   /// dofs. The argument specifies the convergence tolerance.
   void assess_convergence_based_on_relative_solid_change(
    const double& tol)
    {
     Convergence_criterion=Assess_convergence_based_on_relative_solid_change;
     Convergence_tolerance=tol;
    }
      
   /// \short Assess convergence based on max. relative change of solid
   /// dofs. This interface has no argument and the default 
   /// convergence tolerance
   /// for the Newton solver, Problem::Newton_solver_tolerance is used.
   void assess_convergence_based_on_relative_solid_change()
    {
     assess_convergence_based_on_relative_solid_change(
      Problem::Newton_solver_tolerance);
    }


   /// \short Use pointwise Aitken extrapolation. The argument is used to 
   /// specify the Picard iteration after which pointwise Aitken extrapolation 
   /// is to be used for the first time. 
   void enable_pointwise_aitken(const unsigned& pointwise_aitken_start)
    {
     Pointwise_aitken_start = pointwise_aitken_start;
     Use_pointwise_aitken = true;
    }

   /// \short Use pointwise Aitken extrapolation. This interface has
   /// no argument and the current value of Pointwise_aitken_start will
   /// be used. The default is zero, extrapolation starts immediately
   void enable_pointwise_aitken() {Use_pointwise_aitken = true;}

   /// \short Disable the use of pointwise Aitken extrapolation
   void disable_pointwise_aitken() {Use_pointwise_aitken = false;}

   ///\short Use under-relaxation and (optionally) specify under-relaxation 
   /// parameter. Default: omega=1.0 (i.e. no actual under-relaxation; 
   /// Other extreme: omega=0.0 (freeze wall shape). Under-relaxation
   /// parameter can also be computed dynamically by setting
   /// use_irons_and_tuck_extrapolation()
   void enable_under_relaxation(const double& omega=1.0)
    {Omega_relax=omega;}

   ///\short Use Irons and Tuck extrapolation for solid dofs
   void enable_irons_and_tuck_extrapolation()
   {Use_irons_and_tuck_extrapolation = true;}

   ///\short Do not use Irons and Tuck extrapolation for solid dofs
   void disable_irons_and_tuck_extrapolation()
   {Use_irons_and_tuck_extrapolation = false;}
   
   /// Enumerated flags for convergence criteria
   enum convergence_criteria{Assess_convergence_based_on_absolute_solid_change,
                             Assess_convergence_based_on_relative_solid_change,
                             Assess_convergence_based_on_max_global_residual};


   /// Enumerated flags to indicate which solve is taking place
   enum solve_type{Full_solve,Fluid_solve,Solid_solve};

   /// \short Get rms of change in the solid dofs; the max. change of the
   /// solid dofs and the rms norm of the solid dofs themselves.
   /// Change is computed relative to the reference values stored when
   /// store_solid_dofs() was last called.
   void get_solid_change(double& rms_change, 
                         double& max_change,
                         double& rms_norm);

   /// \short Store the current solid values as reference values for 
   /// future convergence check. Also add another entry to pointwise 
   /// Aitken history if required.
   void store_solid_dofs();

   /// \short Reset timer
   void reset_timer()
    {
     T_spent_on_actual_solve=0.0;
     T_ref=clock();
     Timer_has_been_halted=false;
    }


   /// \short (Re-)start timer (e.g. after completing non-essential
   /// parts of the code such as documentation of the iteration's
   /// progress)
   void restart_timer()
    {
     T_ref=clock();
     Timer_has_been_halted=false;
    }


   /// \short Halt timer (e.g. before performing non-essential
   /// parts of the code such as documentation of the iteration's
   /// progress)
   void halt_timer()
    {
     if (!Timer_has_been_halted)
      {
       T_spent_on_actual_solve+=double(clock()-T_ref)/CLOCKS_PER_SEC;
       Timer_has_been_halted=true;
      }
    }


   /// \short Total elapsed time since start of solve
   double t_spent_on_actual_solve()
    {
     halt_timer();
     double time=T_spent_on_actual_solve;
     restart_timer();
     return time;
    }


    protected:

   /// Rebuild global mesh for monolithic discretisation
   void rebuild_monolithic_mesh();

   /// \short Number of Aitken histories available (int because after
   /// extrapolation it's re-initialised to -1 to force the computation
   /// of three new genuine iterates).
   int Pointwise_aitken_counter;

   /// Use pointwise Aitken extrapolation?
   bool Use_pointwise_aitken;

   /// \short Start pointwise Aitken extrpolation after specified number 
   /// of Picard iterations
   unsigned Pointwise_aitken_start;

   /// Solve that is taking place (enumerated flag)
   int Solve_type;
   
   /// Convergence tolerance for Picard iteration
   double Convergence_tolerance;

   /// Max. number of Picard iterations
   unsigned Max_picard;

   /// Doc maximum global residual during iteration? (default: false)
   bool Doc_max_global_residual;

   /// Restore pinned status of fluid dofs
   void restore_fluid_dofs();

   /// Pin solid dofs
   void pin_solid_dofs();

   /// Restore pinned status of solid dofs
   void restore_solid_dofs();
 

   /// Do pointwise Aitken extrapolation for solid
   void pointwise_aitken_extrapolate();
 
   /// \short Vector storing the Data objects associated with the fluid
   /// problem: Tyically the nodal and internal data of the elements in the
   /// fluid bulk mesh 
   Vector<Data*> Fluid_data_pt;

   /// \short Vector of vectors that store the pinned status of
   /// fluid Data values
   Vector<std::vector<bool> > Fluid_value_is_pinned;

   /// \short Vector storing the Data objects associated with the solid
   /// problem: Typically the positional data of solid nodes and
   /// any quantities associated with displacement control, say.
   Vector<Data*> Solid_data_pt;

   /// \short Vector of vectors that store the pinned status of
   /// solid Data values
   Vector<std::vector<bool> > Solid_value_is_pinned;

   /// \short Vector storing the previous solid values -- used for
   /// convergence check
   Vector<double> Previous_solid_value;

   /// \short Mesh containing only fluid elements -- the elements in this
   /// Mesh will be excluded from the assembly process when
   /// the solid problem is solved
   Mesh* Fluid_mesh_pt;

   /// \short Mesh containing only solid elements -- the elements in this
   /// mesh will be excluded from the assembly process when
   /// the fluid problem is solved
   Mesh* Solid_mesh_pt;
   
   /// \short Backup for the pointers to the submeshes in the original problem
   Vector<Mesh*> Orig_sub_mesh_pt;

   /// Vector of changes in Irons and Tuck under-relaxation
   Vector<double> Del_irons_and_tuck;

   /// Irons and Tuck relaxation factor
   double R_irons_and_tuck;

   /// \short Vector of Vectors containing up to three previous
   /// iterates for the solid dofs; used for pointwise Aitken extrapolation
   Vector<Vector<double> > Pointwise_aitken_solid_value;

   /// Have we just done a pointwise Aitken step
   bool Recheck_convergence_after_pointwise_aitken;

    private:

   /// Extrapolate solid data and update fluid mesh during unsteady run
   void extrapolate_solid_data();

   /// \short Under-relax the most recently computed solid variables, either
   /// by classical relaxation or by Irons & Tuck
   void under_relax_solid();

   /// \short Only include fluid elements in the Problem's mesh. This is
   /// called before the segregated fluid solve. The fluid elements are 
   /// identified by the user via the fluid_mesh_pt argument
   /// in the pure virtual function identify_fluid_and_solid_dofs(...). 
   /// If a NULL pointer is returned by this function (i.e. if the user 
   /// hasn't bothered to identify the fluids elements in a submesh, then 
   /// no stripping of non-fluid elements is performed. The result
   /// of the computation will be correct but
   /// it won't be very efficient. 
   void use_only_fluid_elements();

   /// \short Only include solid elements in the Problem's mesh. This is
   /// called before the segregated solid solve. The solid elements are 
   /// identified by the user via the solid_mesh_pt argument
   /// in the pure virtual function identify_fluid_and_solid_dofs(...). 
   /// If a NULL pointer is returned by this function (i.e. if the user 
   /// hasn't bothered to identify the solid elements in a submesh, then 
   /// no stripping of non-solid elements is performed. The result
   /// of the computation will be correct but
   /// it won't be very efficient. 
   void use_only_solid_elements();

   /// Pin fluid dofs
   void pin_fluid_dofs();

   /// \short Under-relaxation parameter. (1.0: no under-relaxation; 
   /// 0.0: Freeze wall shape)
   double Omega_relax; 

   /// \short Boolean flag to indicate use of Irons and Tuck's extrapolation 
   /// for solid values
   bool Use_irons_and_tuck_extrapolation;

   /// Convergence criterion (enumerated flag)
   int Convergence_criterion;
  
   /// \short Reference time for segregated solve. Can be 
   /// re-initialised whenever total elapsed time has been stored
   /// (before entering non-essential doc sections of the code)
   clock_t T_ref;

   /// \short Total elapsed time since start of solve, can be
   /// accumulated by adding bits of time spent in relevant parts of
   /// code (bypassing sections that only document the progress)
   double T_spent_on_actual_solve;

   /// \short boolean flag to indicate if timer has been halted
   bool Timer_has_been_halted;

  };

}



#endif