segregated_fsi_solver.cc

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//====================================================================

#include<algorithm>

#include "segregated_fsi_solver.h"
#include <algorithm>

namespace oomph
{

//============================================================================
/// Extrapolate solid data and update fluid mesh. Extrapolation is based
/// on previous history values and is therefore only called in 
/// time-dependent runs.
//============================================================================
 void SegregatableFSIProblem::extrapolate_solid_data()
 {
  // Number of solid Data items:
  unsigned n_data=Solid_data_pt.size();
   
  // Loop over all solid Data items:
  for (unsigned i=0;i<n_data;i++)
   {
    // Number of values stored at this Data item:
    unsigned n_value=Solid_data_pt[i]->nvalue();
     
    // Number of values stored at this Data item:
    unsigned n_prev=Solid_data_pt[i]->time_stepper_pt()->nprev_values();
      
    // Can't use this extrapolation if we don't have at least two
    // history values; may be able to do better for higher order
    // timesteppers, though.
    if (n_prev>=2)
     {
      // Extrapolate all values
      for (unsigned k=0;k<n_value;k++)
       {
        // Linear extrapolation based on previous two values,
        // assuming constant timestep.
        double new_value=2.0*Solid_data_pt[i]->value(1,k)-
         Solid_data_pt[i]->value(2,k);
        *(Solid_data_pt[i]->value_pt(0,k))=new_value;
       }
     }
   }

  // Update fluid node position in response to changes in wall shape
  // (If fluid mesh was specified via non-NULL Fluid_mesh_pt
  // this node update will only act on the fluid nodes).
  if (Fluid_mesh_pt!=0)
   {
    Fluid_mesh_pt->node_update();
   }
  else
   {
    mesh_pt()->node_update();
   }

 }


//============================================================================
/// Rebuild global mesh for monolithic problem
//============================================================================
 void SegregatableFSIProblem::rebuild_monolithic_mesh()
 {
  // Get rid of the previous submeshes
  flush_sub_meshes();

  // Add original submeshes
  unsigned orig_n_sub_mesh=Orig_sub_mesh_pt.size();
  for (unsigned i=0;i<orig_n_sub_mesh;i++)
   {
    add_sub_mesh(Orig_sub_mesh_pt[i]);
   }

  // Rebuild global mesh
  rebuild_global_mesh();
 }



//============================================================================
/// 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 solids 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 SegregatableFSIProblem::use_only_solid_elements()
 {

  // Wipe global mesh and rebuild it with just the solid elements in it
  if (Solid_mesh_pt!=0)
   {
    flush_sub_meshes();
    add_sub_mesh(Solid_mesh_pt);
    rebuild_global_mesh();
   }

 }



//============================================================================
/// 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 SegregatableFSIProblem::use_only_fluid_elements()
 {

  // Wipe global mesh and rebuild it with just the fluid elements in it
  if (Fluid_mesh_pt!=0)
   {
    flush_sub_meshes();
    add_sub_mesh(Fluid_mesh_pt);
    rebuild_global_mesh();
   }
 }



//============================================================================
/// Pin all fluid dofs
//============================================================================
 void SegregatableFSIProblem::pin_fluid_dofs()
 {

  // Number of fluid Data items:
  unsigned n_data=Fluid_data_pt.size();

  // Loop over all fluid Data items:
  for (unsigned i=0;i<n_data;i++)
   {
    // Number of values stored at this Data item:
    unsigned n_value=Fluid_data_pt[i]->nvalue();
   
    // Pin all values
    for (unsigned k=0;k<n_value;k++)
     {
      Fluid_data_pt[i]->pin(k);
     }
   }
 }



//============================================================================
/// Restore the pinned status of all fluid dofs
//============================================================================
 void SegregatableFSIProblem::restore_fluid_dofs()
 {

  // Number of fluid Data items:
  unsigned n_data=Fluid_data_pt.size();

  // Loop over all fluid Data items:
  for (unsigned i=0;i<n_data;i++)
   {
    // Number of values stored at this Data item:
    unsigned n_value=Fluid_data_pt[i]->nvalue();
   
    // Pin all values
    for (unsigned k=0;k<n_value;k++)
     {
      if (Fluid_value_is_pinned[i][k])
       {
        Fluid_data_pt[i]->pin(k);
       }
      else
       {
        Fluid_data_pt[i]->unpin(k);
       }
     }
   }

 }



//============================================================================
/// Pin all solid dofs
//============================================================================
 void SegregatableFSIProblem::pin_solid_dofs()
 {

  // Number of solid Data items:
  unsigned n_data=Solid_data_pt.size();

  // Loop over all solid Data items:
  for (unsigned i=0;i<n_data;i++)
   {
    // Number of values stored at this Data item:
    unsigned n_value=Solid_data_pt[i]->nvalue();
   
    // Pin all values
    for (unsigned k=0;k<n_value;k++)
     {
      Solid_data_pt[i]->pin(k);
     }
   }

 }



//============================================================================
/// Restore the pinned status of all solid dofs
//============================================================================
 void SegregatableFSIProblem::restore_solid_dofs()
 {

  // Number of solid Data items:
  unsigned n_data=Solid_data_pt.size();

  // Loop over all solid Data items: 
  for (unsigned i=0;i<n_data;i++)
   {
    // Number of values stored at this Data item:
    unsigned n_value=Solid_data_pt[i]->nvalue();
   
    // Pin all values
    for (unsigned k=0;k<n_value;k++)
     {
      if (Solid_value_is_pinned[i][k])
       {
        Solid_data_pt[i]->pin(k);
       }
      else
       {
        Solid_data_pt[i]->unpin(k);
       }
     }
   }
 }





//============================================================================
/// Store the current solid values as reference values for future
/// convergence check. Also add another entry to pointwise Aitken history. 
//============================================================================
 void SegregatableFSIProblem::store_solid_dofs()
 {

  // Counter for the number of values:
  unsigned value_count=0;

  // Number of solid Data items:
  unsigned n_data=Solid_data_pt.size();

  // Loop over all solid Data items:
  for (unsigned i=0;i<n_data;i++)
   {
    // Number of values stored at this Data item:
    unsigned n_value=Solid_data_pt[i]->nvalue();
   
    // Loop over all values
    for (unsigned k=0;k<n_value;k++)
     {
      // Make backup
      Previous_solid_value[value_count]=Solid_data_pt[i]->value(k);
     
      // Store in pointwise Aitken history
      if (Use_pointwise_aitken&&(Pointwise_aitken_counter>=0))
       {
        Pointwise_aitken_solid_value[value_count][Pointwise_aitken_counter]=
         Solid_data_pt[i]->value(k);
       }

      // Increment counter
      value_count++;
     }
   }

  // We stored another level of Aitken history values
  Pointwise_aitken_counter++;
 }





//============================================================================
/// 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 SegregatableFSIProblem::get_solid_change(
  double& rms_change,  double& max_change, double& rms_norm)
 {

  // Initialise
  rms_change=0.0;
  max_change=0.0;
  rms_norm=0.0;

  // Counter for the number of values:
  unsigned value_count=0;

  // Number of solid Data items:
  unsigned n_data=Solid_data_pt.size();

  // Loop over all solid Data items:
  for (unsigned i=0;i<n_data;i++)
   {
    // Number of values stored at this Data item:
    unsigned n_value=Solid_data_pt[i]->nvalue();
   
    // Loop over all values
    for (unsigned k=0;k<n_value;k++)
     {

      // Change
      double change=Previous_solid_value[value_count]-
       Solid_data_pt[i]->value(k);

      // Max change?
      if (std::fabs(change)>max_change) max_change=std::fabs(change);

      //Add square of change relative to previous value
      rms_change+=pow(change,2);

      //Add square of previous value
      rms_norm+=pow(Previous_solid_value[value_count],2);

      // Increment counter
      value_count++;
     }
   }

  // Turn into rms:
  rms_change=sqrt(rms_change/double(value_count));
  rms_norm=sqrt(rms_norm/double(value_count));
 }




//============================================================================
/// Under-relax solid, either by classical relaxation or by Irons & Tuck.
//============================================================================
 void SegregatableFSIProblem::under_relax_solid()
 {

  // Irons and Tuck extrapolation/relaxation; an extension of Aitken's method
  //-------------------------------------------------------------------------
  if (Use_irons_and_tuck_extrapolation)
   {
    double top=0.0;
    double den=0.0;
    double crit=0.0;

    // Counter for the number of values:
    unsigned value_count=0;
     
    // Number of solid Data items:
    unsigned n_data=Solid_data_pt.size();
     
    // Loop over all solid Data items:
    for (unsigned i=0;i<n_data;i++)
     {
      // Number of values stored at this Data item:
      unsigned n_value=Solid_data_pt[i]->nvalue();
       
      // Loop over all values
      for (unsigned k=0;k<n_value;k++)
       {         
        // Prediction from solid solver
        double new_solid_value=Solid_data_pt[i]->value(k);
         
        // Previus value
        double old_solid_value=Previous_solid_value[value_count];
         
        // Change
        double change=old_solid_value-new_solid_value;

        // Change of change
        double del2=Del_irons_and_tuck[value_count]-change;

        // Update change
        Del_irons_and_tuck[value_count]=change;

        // Update top
        top+=del2*change;
         
        // Update denominator
        den+=del2*del2;

        // Update convergence criterion
        crit+=std::fabs(change);

        // Increment counter
        value_count++;
       }
     }

    // Update relaxation factor. The if buffers the case in which
    // we haven't realised that we've converged (so that den=0).
    // This can happen, e.g. if the convergence assessment is based on the
    // global residual or during validation. In that case we 
    // obviously don't want any changes to the iterates.
    if (den!=0.0)
     {
      double new_r=R_irons_and_tuck+(R_irons_and_tuck-1.0)*top/den;
      R_irons_and_tuck=new_r;
     }
    else
     {
      R_irons_and_tuck=0.0;
     }

    // Loop over all solid Data items for update
    value_count=0;
    for (unsigned i=0;i<n_data;i++)
     {
      // Number of values stored at this Data item:
      unsigned n_value=Solid_data_pt[i]->nvalue();
       
      // Loop over all values
      for (unsigned k=0;k<n_value;k++)
       {
        // Compute relaxed/extrapolated value
        double new_value=Solid_data_pt[i]->value(k)+
         R_irons_and_tuck*Del_irons_and_tuck[value_count];

        // Assign 
        Solid_data_pt[i]->set_value(k,new_value);
         
        // Increment counter
        value_count++;
       }
     }
    return;
   }

  // Standard relaxation
  //--------------------
  else
   {
     
     
    // No relaxation: Can return immediately
    if (Omega_relax==1.0) return;
     
    // Counter for the number of values:
    unsigned value_count=0;
     
    // Number of solid Data items:
    unsigned n_data=Solid_data_pt.size();
     
    // Loop over all solid Data items:
    for (unsigned i=0;i<n_data;i++)
     {
      // Number of values stored at this Data item:
      unsigned n_value=Solid_data_pt[i]->nvalue();
       
      // Loop over all values
      for (unsigned k=0;k<n_value;k++)
       {
         
        // Prediction from solid solver
        double new_solid_value=Solid_data_pt[i]->value(k);
         
        // Previus value
        double old_solid_value=Previous_solid_value[value_count];
         
        // Relax
        Solid_data_pt[i]->set_value(k,Omega_relax*new_solid_value+
                                    (1.0-Omega_relax)*old_solid_value);
                  
        // Increment counter
        value_count++;
       }
     }
   }
 }

//============================================================================
/// Pointwise Aitken extrapolation for solid variables
//============================================================================
 void SegregatableFSIProblem::pointwise_aitken_extrapolate()
 {

  // Counter for the number of values:
  unsigned value_count=0;

  // Number of solid Data items:
  unsigned n_data=Solid_data_pt.size();

  // Loop over all solid Data items:
  for (unsigned i=0;i<n_data;i++)
   {
    // Number of values stored at this Data item:
    unsigned n_value=Solid_data_pt[i]->nvalue();
   
    // Loop over all values
    for (unsigned k=0;k<n_value;k++)
     {
      // Shorthand
      double v0=Pointwise_aitken_solid_value[value_count][0];
      double v1=Pointwise_aitken_solid_value[value_count][1];
      double v2=Pointwise_aitken_solid_value[value_count][2];

      double new_value=v2;

      double max_diff=std::max(std::fabs(v1-v0),std::fabs(v2-v1));
      if (max_diff>1.0e-10)
       {
        new_value=v2-std::pow((v2-v1),int(2))/(v2-2.0*v1+v0);
       }
      Solid_data_pt[i]->set_value(k,new_value);
      
      // Increment counter
      value_count++;
     }
   }

  // Reset the counter for the Aitken convergence check
  // (setting counter to -1 forces three new genuine
  // iterates to be computed). 
  Pointwise_aitken_counter=-1;
 
 }


//============================================================================
/// Segregated fixed point solver with optional pointwise Aitken acceleration
/// on selected solid dofs. Returns PicardConvergenceData object 
/// that contains the vital stats of the iteration.
//============================================================================
 PicardConvergenceData SegregatableFSIProblem::segregated_solve()
 {
  // Initialise timer for essential bits of code
  reset_timer();
  // Start timer for overall solve
  clock_t t_total_start = clock();

  // If convergence is based on max. global residual we may as well
  // document it...
  bool get_max_global_residual=Doc_max_global_residual;
  if(Convergence_criterion==Assess_convergence_based_on_max_global_residual)
   {
    get_max_global_residual=true;
   }

  // Create object to doc convergence stats
  PicardConvergenceData conv_data;
   
  // Initialise total time for computation of global residuals
  double cpu_for_global_residual=0.0;
   
  //Update anything that needs updating
  actions_before_segregated_solve();
   
  // Set flags to values that are appropriate if Picard iteration
  // does not converge with Max_picard iterations
  bool converged=false;
  unsigned iter_taken=Max_picard;

  // This one will be overwritten during the convergence checks
  double tol_achieved=0.0;
         
  // Store the current values of the solid dofs as reference values
  // and for the pointwise Aitken extrapolation 
  store_solid_dofs();

  // Loop over Picard iterations
  //----------------------------
  for (unsigned iter=1;iter<=Max_picard;iter++)
   {

    //Calculate the initial residuals
    if(iter==1)
     {
      //Problem is always non-linear?
      //Perform any actions before the convergence check
      actions_before_segregated_convergence_check();
       
      double max_res = 0.0;
      if (get_max_global_residual)
       {
        clock_t t_start = clock(); 
        restore_solid_dofs();
        restore_fluid_dofs();
        rebuild_monolithic_mesh();
        assign_eqn_numbers();
        //This is now a full problem
        Solve_type = Full_solve;
        DoubleVector residual;
        get_residuals(residual);
        //Get maximum residuals using our own abscmp function
        max_res = residual.max();
        clock_t t_end = clock();
        cpu_for_global_residual+=double(t_end-t_start)/CLOCKS_PER_SEC;
       }
       
      oomph_info << "==================================================\n";
      oomph_info <<   "Initial iteration             : " 
                 << 0 << std::endl;
      oomph_info <<   "RMS  change           :       "
                 << 0 << std::endl;
      oomph_info <<   "Max. change           :       "
                 << 0 << std::endl;
      oomph_info <<   "RMS norm              :       "
                 << 0   << std::endl;
      if (get_max_global_residual)
       {
        oomph_info << "Max. global residual  :       "
                   << max_res   << std::endl;
       }
      oomph_info << "==================================================\n\n";
       
      // Check for convergence, but this only makes sense
      // for the maximum (rather than relative case)
      if((Convergence_criterion==
          Assess_convergence_based_on_max_global_residual) &&
         (max_res < Convergence_tolerance))
       {
        oomph_info
         << "\n\n\n////////////////////////////////////////////////////////"
         << "\nPicard iteration converged after " 
         << 0 << " steps!" << std::endl
         << "Convergence was based on max. residual of coupled eqns \n"
         << "being less than " << Convergence_tolerance << "." << std::endl
         << "////////////////////////////////////////////////////////\n\n\n"
         << std::endl;

        //Converged, so bail out
        // Number of iterations taken
        iter_taken=0;
           
        // We have converged (overwrites default of false)
        conv_data.set_solver_converged();
           
        // Break loop using a GO TO! This is the first (and hopefully 
        // the last) one in the whole of oomph-lib. Here it's 
        // it's actually the cleanest way to exit from these
        // nested control structures
        goto jump_out_of_picard;
       }
     }
     
    // Stage 1: Freeze wall and solve single-physics fluid problem
    //------------------------------------------------------------
    // Pin the solid dofs, restore the pinned status of the fluid dofs
    // and re-assign the equation numbers
    restore_fluid_dofs();
    pin_solid_dofs();
    use_only_fluid_elements();
    assign_eqn_numbers();
     
    // Solve the fluid problem for the current wall geometry
    oomph_info <<"\n\nDOING FLUID SOLVE\n\n";
    //This is a fluid solve at the moment done by a newton solve
    Solve_type=Fluid_solve;
    newton_solve();
     
    // Stage 2: Freeze fluid variables and update wall shape
    //------------------------------------------------------
     
    // Now restore the pinned status of the wall and pin the fluid
    // dofs in anticipation of the wall update:
    restore_solid_dofs();
    pin_fluid_dofs();
    use_only_solid_elements(); 
    assign_eqn_numbers();
     
    // Solve the solid problem for the wall solution
    oomph_info <<"\n\nDOING SOLID SOLVE\n\n";
    //This is a solid solve, at the moment done by a newton method
    Solve_type = Solid_solve;
    newton_solve();

    // Under-relax, either by classical relaxation or by Irons & Tuck
    // Note that we are under-relaxing before the convergence check
    // because under-relaxtion may be required to acheieve any
    // kind of convergence. If the convergence check is on the RELATIVE
    // change of the residuals, however, then a small under-relaxation
    // parameter will cause a false convergence. You have been warned!
    under_relax_solid();
     
    // Stage 3: Convergence check (possibly again after pointwise Aitken 
    //------------------------------------------------------------------
    // extrapolation)
    //---------------
    //Assume that we are not doing Aitken acceleration
    Recheck_convergence_after_pointwise_aitken=false;
    do
     {
      //Perform any actions before the convergence check
      actions_before_segregated_convergence_check();
       
      //Get the change in the solid variables
      double rms_change;
      double max_change;
      double rms_norm;
      double max_res=0.0;
      get_solid_change(rms_change,max_change,rms_norm);
       
       
      //If we are computing the maximum global residual, do so
      if (get_max_global_residual)
       {
        clock_t t_start = clock();

        // Free all dofs so the residuals of the global equations are computed
        restore_solid_dofs();
        restore_fluid_dofs(); 
        rebuild_monolithic_mesh();
        assign_eqn_numbers();

        //We are now treating this as a full solve
        Solve_type = Full_solve;

        //Get the residuals
        DoubleVector residual;
        get_residuals(residual);
         
        //Get maximum residuals, using our own abscmp function
        max_res =  residual.max();
         
        clock_t t_end = clock();
        cpu_for_global_residual+=double(t_end-t_start)/CLOCKS_PER_SEC;
       }
       
      oomph_info << "==================================================\n";
      oomph_info <<   "Iteration             : " 
                 << iter << std::endl;
      oomph_info <<   "RMS  change           :       "
                 << rms_change << std::endl;
      oomph_info <<   "Max. change           :       "
                 << max_change << std::endl;
      oomph_info <<   "RMS norm              :       "
                 << rms_norm   << std::endl;
      if (get_max_global_residual)
       {
        oomph_info << "Max. global residual  :       "
                   << max_res   << std::endl;
       }
      oomph_info << "==================================================\n\n";
     
      // Check for convergence
      switch (Convergence_criterion)
       {
       
       case Assess_convergence_based_on_absolute_solid_change:
        tol_achieved=max_change;
        if (tol_achieved<Convergence_tolerance)
         {
          oomph_info 
           << "\n\n\n////////////////////////////////////////////////////////"
           << "\nPicard iteration converged after " 
           << iter << " steps!" << std::endl
           << "Convergence was based on absolute change in solid dofs \n"
           << "being less than " << Convergence_tolerance << "." << std::endl
           << "////////////////////////////////////////////////////////\n\n\n"
           << std::endl;
          converged=true;
         }
        break;
         
         
       case Assess_convergence_based_on_relative_solid_change:
        tol_achieved=std::fabs(rms_change/rms_norm);
        if (tol_achieved<Convergence_tolerance)
         {
          oomph_info
           << "\n\n\n///////////////////////////////////////////////////////" 
           << "\nPicard iteration converged after " 
           << iter << " steps!" << std::endl
           << "Convergence was based on relative change in solid dofs \n"
           << "being less than " << Convergence_tolerance << "." << std::endl
           << "////////////////////////////////////////////////////////\n\n\n"
           << std::endl;
          converged=true;
         }
        break;
         
       case Assess_convergence_based_on_max_global_residual:
        tol_achieved=max_res;
        if (tol_achieved<Convergence_tolerance)
         {
          oomph_info
           << "\n\n\n////////////////////////////////////////////////////////"
           << "\nPicard iteration converged after " 
           << iter << " steps!" << std::endl
           << "Convergence was based on max. residual of coupled eqns \n"
           << "being less than " << Convergence_tolerance << "." << std::endl
           << "////////////////////////////////////////////////////////\n\n\n"
           << std::endl;
          converged=true;
         }
        break;

       }

      // If converged bail out
      if (converged)
       {
        // Number of iterations taken
        iter_taken=iter;

        // We have converged (overwrites default of false)
        conv_data.set_solver_converged();

        // Break loop using a GO TO! This is the first (and hopefully 
        // the last) one in the whole of oomph-lib. Here it's 
        // it's actually the cleanest way to exit from these
        // nested control structures
        goto jump_out_of_picard;
       }
      
      // Store the current values of the solid dofs as reference values
      // and for the pointwise Aitken extrapolation 
      store_solid_dofs();
     
      // Correct wall shape by pointwise Aitken extrapolation if possible
      //-----------------------------------------------------------------
      // and desired
      //------------
      //This is an acceleration method for the (possibly under-relaxed)
      //sequence of iterates.
      if ((3==Pointwise_aitken_counter)&&(Use_pointwise_aitken)&&
          (iter>Pointwise_aitken_start))
       {
        pointwise_aitken_extrapolate();
         
        // Repeat the convergence check
        Recheck_convergence_after_pointwise_aitken=true;
       }
      else
       {
        // Didn't change anything: Don't repeat the convergence check
        Recheck_convergence_after_pointwise_aitken=false;
       } 
     }
    //Repeat convergence while we are doing aitken extrapolation
    while(Recheck_convergence_after_pointwise_aitken);
   
   } // End of loop over iterations
   
   

  // Jump address for converged or diverged iteration
   jump_out_of_picard:

  // Reset everything
  restore_solid_dofs();
  restore_fluid_dofs();
  rebuild_monolithic_mesh();
  assign_eqn_numbers();
  //This is a full solve again
  Solve_type = Full_solve;

  // Do any updates that are required 
  actions_after_segregated_solve();
   
  // Number of iterations (either this is still Max_iter from 
  // the initialisation or it's been overwritten on convergence)
  conv_data.niter()=iter_taken;
   
  // Total cpu time
  clock_t t_total_end = clock();
  conv_data.cpu_total()=double(t_total_end-t_total_start)/
   CLOCKS_PER_SEC;
   
  // Total essential cpu time (exluding doc etc)
  conv_data.essential_cpu_total()=t_spent_on_actual_solve();

  // cpu time for check/doc of global residual
  conv_data.cpu_for_global_residual()=cpu_for_global_residual;
   
  // Final tolerance achieved by the iteration
  conv_data.tol_achieved()=tol_achieved;

  // Doc non-convergence
  if (!converged)
   {
    switch (Convergence_criterion)
     {
       
     case Assess_convergence_based_on_absolute_solid_change:
      oomph_info 
       << "\n\n\n////////////////////////////////////////////////////////"
       << "\nPicard iteration did not converge after " 
       << iter_taken << " steps!" << std::endl
       << "Convergence was based on absolute change in solid dofs \n"
       << "being less than " << Convergence_tolerance << " " << std::endl
       << "but we achieved only " << tol_achieved << "." << std::endl
       << "////////////////////////////////////////////////////////\n\n\n"
       << std::endl;
      //Throw an error indicating if we ran out of iterations
      if (iter_taken==Max_picard)
       {
        throw SegregatedSolverError(true);
       }
      else 
       {
        throw SegregatedSolverError(false);
       }
      break;
       
       
     case Assess_convergence_based_on_relative_solid_change:
      oomph_info
       << "\n\n\n///////////////////////////////////////////////////////" 
       << "\nPicard iteration did not converge after " 
       << iter_taken << " steps!" << std::endl
       << "Convergence was based on relative change in solid dofs \n"
       << "being less than " << Convergence_tolerance << " " << std::endl
       << "but we achieved only " << tol_achieved << "." << std::endl
       << "////////////////////////////////////////////////////////\n\n\n"
       << std::endl;
      //Throw an error indicating if we ran out of iterations
      if (iter_taken==Max_picard)
       {
        throw SegregatedSolverError(true);
       }
      else 
       {
        throw SegregatedSolverError(false);
       }
      break;
       
     case Assess_convergence_based_on_max_global_residual:
      oomph_info
       << "\n\n\n////////////////////////////////////////////////////////"
       << "\nPicard iteration did not converge after " 
       << iter_taken << " steps!" << std::endl
       << "Convergence was based on max. residual of coupled eqns \n"
       << "being less than " << Convergence_tolerance << " " << std::endl
       << "but we achieved only " << tol_achieved << "." << std::endl
        
       << "////////////////////////////////////////////////////////\n\n\n"
       << std::endl;
      //Throw an error indicating if we ran out of iterations
      if (iter_taken==Max_picard)
       {
        throw SegregatedSolverError(true);
       }
      else 
       {
        throw SegregatedSolverError(false);
       }
      break;
       
     }
   }

  return conv_data;
 }



//============================================================================
/// Segregated fixed point solver with optional pointwise Aitken acceleration
/// on selected solid dofs. Returns PicardConvergenceData object 
/// that contains the vital stats of the iteration.
//============================================================================
 PicardConvergenceData SegregatableFSIProblem::steady_segregated_solve()
 {
  //Find out how many timesteppers there are
  unsigned n_time_steppers = ntime_stepper();
   
  // Vector of bools to store the is_steady status of the various
  // timesteppers when we came in here
  std::vector<bool> was_steady(n_time_steppers);

  //Loop over them all and make them (temporarily) static
  for(unsigned i=0;i<n_time_steppers;i++)
   {
    was_steady[i]=time_stepper_pt(i)->is_steady();
    time_stepper_pt(i)->make_steady();
   }

  // Create object to doc convergence stats
  PicardConvergenceData conv_data;
   
  //Solve the non-linear problem by the segregated solver
  try
   {
    conv_data = segregated_solve();
   }
  //Catch any exceptions thrown in the segregated solver
  catch(SegregatedSolverError &error)
   {
    if (!error.Ran_out_of_iterations)
     {
      std::ostringstream error_stream;
      error_stream << "Error occured in Segregated solver. "
                   << std::endl;
      throw OomphLibError(error_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
    else
     {
      oomph_info << "Note: Ran out of iterations but continuing anyway"
                 << std::endl;
     }
   }
   
  // Reset the is_steady status of all timesteppers that
  // weren't already steady when we came in here and reset their
  // weights
  for(unsigned i=0;i<n_time_steppers;i++)
   {
    if (!was_steady[i])
     {
      time_stepper_pt(i)->undo_make_steady();
     }
   }
   
  // Since we performed a steady solve, the history values 
  // now have to be set as if we had performed an impulsive start from
  // the current solution. This ensures that the time-derivatives
  // evaluate to zero even now that the timesteppers have been
  // reactivated. 
  assign_initial_values_impulsive();
   
  //Return the convergence data
  return conv_data;
 }


//============================================================================
/// Segregated fixed point solver with optional pointwise Aitken acceleration
/// on selected solid dofs. Returns PicardConvergenceData object 
/// that contains the vital stats of the iteration.
//============================================================================
 PicardConvergenceData SegregatableFSIProblem::unsteady_segregated_solve(
  const double& dt)
 {
  //We shift the values, so shift_values is true
  return unsteady_segregated_solve(dt,true);
 }

//============================================================================
/// Segregated fixed point solver with optional pointwise Aitken acceleration
/// on selected solid dofs. Returns PicardConvergenceData object 
/// that contains the vital stats of the iteration.
//============================================================================
 PicardConvergenceData SegregatableFSIProblem::unsteady_segregated_solve(
  const double& dt, const bool &shift_values)
 {
  //Shift the time values and the dts according to the control flag
  if(shift_values) {shift_time_values();}

  // Advance global time and set current value of dt 
  time_pt()->time()+=dt;
  time_pt()->dt()=dt;
   
  //Find out how many timesteppers there are
  unsigned n_time_steppers = ntime_stepper();
   
  //Loop over them all and set the weights
  for(unsigned i=0;i<n_time_steppers;i++)
   {
    time_stepper_pt(i)->set_weights();
   }
   
  //Now update anything that needs updating before the timestep
  //This could be time-dependent boundary conditions, for example.
  actions_before_implicit_timestep();
   
  // Extrapolate the solid data and then update fluid mesh during unsteady run
  extrapolate_solid_data();
   
  // Create object to doc convergence stats
  PicardConvergenceData conv_data;
   
  try
   {
    //Solve the non-linear problem for this timestep with Newton's method
    conv_data = segregated_solve();
   }
  //Catch any exceptions thrown in the segregated solver
  catch(SegregatedSolverError &error)
   {
    if (!error.Ran_out_of_iterations)
     {
      std::ostringstream error_stream;
      error_stream << "Error occured in Segregated solver. "
                   << std::endl;
      throw OomphLibError(error_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
    else
     {
      oomph_info << "Note: Ran out of iterations but continuing anyway"
                 << std::endl;
     }
   }
   
  //Now update anything that needs updating after the timestep
  actions_after_implicit_timestep();

  return conv_data;
 }


//============================================================================
/// Setup segregated solver, using the information provided by the user
/// in his/her implementation of the pure virtual function 
/// identify_fluid_and_solid_dofs(...).
//============================================================================
 void SegregatableFSIProblem::setup_segregated_solver(
  const bool &full_setup_of_fluid_and_solid_dofs)
 {
  //If we are doing a full setup
  if(full_setup_of_fluid_and_solid_dofs)
   {
    // Identify the fluid and solid Data
    Vector<Data*> fluid_data_pt;
    Vector<Data*> solid_data_pt;
    identify_fluid_and_solid_dofs(fluid_data_pt,solid_data_pt,
                                  Fluid_mesh_pt,Solid_mesh_pt);
    
    if (Fluid_mesh_pt==0)
     {
      oomph_info 
       << std::endl << std::endl 
       << "Warning: Your implementation of the pure virtual\n"
       << "         function identify_fluid_and_solid_dofs(...)\n"
       << "         returned a NULL pointer for Fluid_mesh_pt.\n"
       << "         --> The fluid elements will remain activated\n"
       << "         during the solid solve. This is inefficient!\n"
       << "         You should combine all fluid elements into a combined\n"
       << "         mesh and specify this mesh in your\n"
       << "         implementation of \n\n"
       << "         SegregatableFSIProblem::identify_fluid_and_solid_dofs(...)"
       << std::endl << std::endl;
     }
    
    
    if (Solid_mesh_pt==0)
     {
      oomph_info 
       << std::endl << std::endl 
       << "Warning: Your implementation of the pure virtual\n"
       << "         function identify_fluid_and_solid_dofs(...)\n"
       << "         returned a NULL pointer for Solid_mesh_pt.\n"
       << "         --> The solid elements will remain activated\n"
       << "         during the fluid solve. This is inefficient!\n"
       << "         You should combine all solid elements into a combined\n"
       << "         mesh and specify this mesh in your\n"
       << "         implementation of \n\n"
       << "         SegregatableFSIProblem::identify_fluid_and_solid_dofs(...)"
       << std::endl << std::endl;
     }
    
    // Back up the pointers to the submeshes in the original problem
    // so we can restore the problem when we're done
    unsigned orig_n_sub_mesh=nsub_mesh();
    Orig_sub_mesh_pt.resize(orig_n_sub_mesh);
    for (unsigned i=0;i<orig_n_sub_mesh;i++)
     {
      Orig_sub_mesh_pt[i]=mesh_pt(i);
     }
    
    // Fluid
    //------
    
    // Reset
    Fluid_data_pt.clear();
    Fluid_value_is_pinned.clear();
    
    // Make space for fluid data items
    unsigned n_fluid_data=fluid_data_pt.size();
    Fluid_data_pt.resize(n_fluid_data);
    Fluid_value_is_pinned.resize(n_fluid_data);
    
    // Counter for number of fluid values
    unsigned n_fluid_values=0;
    
    // Loop over fluid data
    for (unsigned i=0;i<n_fluid_data;i++)
     {
      // Copy data
      Fluid_data_pt[i]=fluid_data_pt[i];
      
      // Number of values stored at this Data item:
      unsigned n_value=fluid_data_pt[i]->nvalue();
      Fluid_value_is_pinned[i].resize(n_value);
      
      // Copy pinned status for all values
      for (unsigned k=0;k<n_value;k++)
       {
        Fluid_value_is_pinned[i][k]=
         fluid_data_pt[i]->is_pinned(k);
        
        // Increment counter for number of fluid values
        n_fluid_values++;
       }
     }
    
    
    // Solid
    //------
    
    // Reset
    Solid_data_pt.clear();
    Solid_value_is_pinned.clear();
    Previous_solid_value.clear();
    Pointwise_aitken_solid_value.clear();
    Del_irons_and_tuck.clear();
    
    
    unsigned n_solid_data=solid_data_pt.size();
    
    // Make space for solid data items
    unsigned nsolid_data=solid_data_pt.size();
    Solid_data_pt.resize(nsolid_data);
    Solid_value_is_pinned.resize(nsolid_data);
    
    // Counter for number of solid values
    unsigned n_solid_values=0;
    
    // Loop over solid data
    for (unsigned i=0;i<n_solid_data;i++)
     {
      // Copy data
      Solid_data_pt[i]=solid_data_pt[i];
      
      // Number of values stored at this Data item:
      unsigned n_value=solid_data_pt[i]->nvalue();
      Solid_value_is_pinned[i].resize(n_value);
      
      // Copy pinned status for all values
      for (unsigned k=0;k<n_value;k++)
       {
        Solid_value_is_pinned[i][k]=
         solid_data_pt[i]->is_pinned(k);
        
        // Increment counter for number of solid values
        n_solid_values++;
       }
     }
    
    // Make space for previous solid values
    Previous_solid_value.resize(n_solid_values);
    
    // Allocate storage and initialise Irons and Tuck extrapolation
    if (Use_irons_and_tuck_extrapolation)
     {
      Del_irons_and_tuck.resize(n_solid_values);
     }
    
    // Make space for pointwise Aitken extrapolation
    if (Use_pointwise_aitken)
     {
      Pointwise_aitken_solid_value.resize(n_solid_values);
      for (unsigned i=0;i<n_solid_values;i++)
       {
        Pointwise_aitken_solid_value[i].resize(3);
       }
     }

   } //End of full setup
  

      

  // Initialise Irons and Tuck relaxation factor
  R_irons_and_tuck=1.0-Omega_relax;

  //Initialise Irons and Tuck delta values
  unsigned n_del = Del_irons_and_tuck.size();
  for (unsigned i=0;i<n_del;i++) {Del_irons_and_tuck[i]=1.0e20;}

  // Initialise counter for the number of pointwise Aitken values stored
  Pointwise_aitken_counter=0;
  
 }
}