fsi_jacobian_approximation.cc File Reference

#include <iostream>
#include <cstring>
#include "generic.h"
#include "navier_stokes.h"
#include "beam.h"
#include "meshes/one_d_lagrangian_mesh.h"
#include "meshes/collapsible_channel_mesh.h"

Classes
class	UndeformedWall
class	FSICollapsibleChannelProblem< ELEMENT >
	Problem class. More...

Namespaces
	BL_Squash
	Control_Flags
	Namespace for phyical parameters.
	Global_Physical_Variables
	Global variables.

Functions
double	BL_Squash::squash_fct (const double &s)
double &	Global_Physical_Variables::external_pressure ()
	Access function to value of external pressure. More...
void	Global_Physical_Variables::prescribed_traction (const double &t, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
	Traction applied on the fluid at the left (inflow) boundary. More...
void	Global_Physical_Variables::load (const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
	Load function: Apply a constant external pressure to the beam. More...
int	main (int argc, char *argv[])

Variables
bool	Control_Flags::Validation_run =false
	Normal run or validation run? More...
bool	Control_Flags::Steady_run =false
	Steady run or unsteady run? More...

Function Documentation

main()

int main	(	int argc	,
		char *	argv[]
	)

Driver code for a collapsible channel problem with FSI.

Run with the -help flag for more information on command line arguments to run with an approximation to the Jacobian, to run steady calculations and to perform validation runs with coarse resolution and small number of timesteps.

Disclaimer: this code is intended to demonstrate how to ask for various approximations to be made to the Jacobian matrix which can result in reduced calculations per Newton step (and to test this).However these approximations to the Jacobian may increase the overall number of Newton steps, or even cause the Newton solve to diverge. Therefore for any particular FSI problem some experiementation is required to determine whether approximating the Jacobian is beneficial. Note: In tests using this code, approximating the Jacobain increased the overall time per Newton solve, however for a 3D equivalent to this problem significant speed ups were found.

{
 
 // Store command line arguments
 CommandLineArgs::setup(argc,argv);
 
 // Parse command line
 int arg_index = 0;
 bool print_help = false;
 bool approx_jacobian = false;
 
 while (arg_index < argc)
  {
   if ( strcmp(argv[arg_index], "-validation_run") == 0 )
    {
     arg_index++;
     Control_Flags::Validation_run=true;
    }
   else if ( strcmp(argv[arg_index], "-steady_run") == 0 )
    {
     arg_index++;
     Control_Flags::Steady_run=true;
    }
   else if ( strcmp(argv[arg_index], "-approx_jacobian") == 0 )
    {
     arg_index++;
     approx_jacobian=true;
    }
   else if ( strcmp(argv[arg_index], "-help") == 0 )
    {
     print_help = true;
     break;
    }
   else
    {
     arg_index++;
    }
  }
 
 if (print_help)
  {
   oomph_info << "\n\nOption flags:\n\n";
   oomph_info << "-validation run   Perform validation run with coarse mesh\n"
              << "                  and fewer Newton solves\n\n";
   oomph_info << "-steady_run       Solve series of steady problems with\n"
              << "                  different wall positions, otherwise\n"
              << "                  solve unsteady problem\n\n";
   oomph_info << "-approx_jacobian  Use Jacobian approximation\n\n";
   return (0);
  }
 
 // Set flags controlling approximations to the Jacobian
 bool wall_jacobian_ignores_fluid_shear_stress_data = false;
 bool wall_jacobian_ignores_geometric_data = false;
 bool fluid_jacobian_ignores_geometric_data = false;
 if (approx_jacobian)
  {
   wall_jacobian_ignores_fluid_shear_stress_data = true;
   wall_jacobian_ignores_geometric_data = true;
   
   if (Control_Flags::Steady_run)
    {
     // This can speed up the calculation of the Jacobian but causes the
     // Newton method to fail for unsteady problems
     fluid_jacobian_ignores_geometric_data = true;
    }
  }
 
 // Reduction in resolution for validation run?
 unsigned coarsening_factor=1;
 if (Control_Flags::Validation_run)
  {
   coarsening_factor=4;
  }
 
 // Number of elements in the domain
 unsigned nup=20/coarsening_factor;
 unsigned ncollapsible=40/coarsening_factor;
 unsigned ndown=40/coarsening_factor;
 unsigned ny=16/coarsening_factor;
 
 // Length of the domain
 double lup=5.0;
 double lcollapsible=10.0;
 double ldown=10.0;
 double ly=1.0;
 
 // Build the problem with QTaylorHoodElements
 FSICollapsibleChannelProblem<AlgebraicElement<QTaylorHoodElement<2> > > 
  problem(nup,
          ncollapsible,
          ndown,
          ny, 
          lup,
          lcollapsible,
          ldown,
          ly,
          wall_jacobian_ignores_fluid_shear_stress_data,
          wall_jacobian_ignores_geometric_data,
          fluid_jacobian_ignores_geometric_data);
 
 // Set external pressure (on the wall stiffness scale). 
 Global_Physical_Variables::external_pressure() = 1.0e-1;
 
 // Pressure on the left boundary: This is consistent with steady
 // Poiseuille flow
 Global_Physical_Variables::P_up=12.0*(lup+lcollapsible+ldown);
 
 // Set initial wall displacement
 Global_Physical_Variables::Prescribed_y = 1.0;
 
 // Set Re and ReSt - it appears small values are required in this 2D problem, 
 // with the Jacobian approximations otherwise the Newton method converges, this
 // is not the experience found with the 3D equivalent problem.
 Global_Physical_Variables::Re=1;
 Global_Physical_Variables::ReSt=Global_Physical_Variables::Re;
 
 // Timestep. Note: Preliminary runs indicate that the period of
 // the oscillation is about 1 so this gives us 40 steps per period.
 double dt=1.0/40.0; 
 
 // Initial time for the simulation
 double t_min=0.0;
 
 // Maximum time for simulation
 double t_max=3.5; 
 
 // Initialise timestep 
 problem.time_pt()->time()=t_min;
 problem.initialise_dt(dt);
 
 // Apply initial condition
 problem.set_initial_condition();
 
 //Set output directory
 DocInfo doc_info;
 if (Control_Flags::Validation_run)
  {
   if (approx_jacobian)
    {
     oomph_info << "Using approximate Jacobian\n";
     if (Control_Flags::Steady_run)
      doc_info.set_directory("RESLT_APPROX_STEADY");
     else
      doc_info.set_directory("RESLT_APPROX_UNSTEADY");
    }
   else
    {
     oomph_info << "Using exact Jacobian\n";
     if (Control_Flags::Steady_run)
      doc_info.set_directory("RESLT_EXACT_STEADY");
     else
      doc_info.set_directory("RESLT_EXACT_UNSTEADY");
    }
  }
 else
  {
   doc_info.set_directory("RESLT");
  }
 
 // Open a trace file 
 ofstream trace_file;
 char filename[100];   
 sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
 trace_file.open(filename);
   
 // Output the initial solution
 problem.doc_solution(doc_info, trace_file);
 
 // Increment step number
 doc_info.number()++;
 
 if (Control_Flags::Steady_run)
  {
   // Number of solves (reduced for validation)
   unsigned nstep = 20;
   if (Control_Flags::Validation_run)
    {
     nstep=3;
    }
 
   // Loop over different wall displacement
   for (unsigned istep=0;istep<nstep;istep++)
    {
     oomph_info << "\nSteady Newton solve " << istep
                << " for wall displacement of " 
                << Global_Physical_Variables::Prescribed_y << "\n"
                << "====================================================\n";
     
     // Solve the problem
     problem.newton_solve();
 
     cout << "External pressure="
          << Global_Physical_Variables::external_pressure()
          << "\n";
     
     // Output the solution
     problem.doc_solution(doc_info, trace_file);
     
     // Step number
     doc_info.number()++;
 
     // Change wall position
     Global_Physical_Variables::Prescribed_y -= 0.01;
    }
  }
 else
  {
   // Find number of timesteps (reduced for validation)
   unsigned nstep = unsigned((t_max-t_min)/dt);
   if (Control_Flags::Validation_run)
    {
     nstep=10;
    }
   
   // Timestepping loop
   for (unsigned istep=0;istep<nstep;istep++)
    {
     oomph_info << "\nNewton solve " << istep << "\n"
                << "================\n";
     
     // Solve the problem
     problem.unsteady_newton_solve(dt);
     
     // Output the solution
     problem.doc_solution(doc_info, trace_file);
     
     // Step number
     doc_info.number()++;
    }
   
  }
 
 // Close trace file.
 trace_file.close();
 
}//end of main

References oomph::DocInfo::directory(), Global_Physical_Variables::external_pressure(), MergeRestartFiles::filename, Mesh_Parameters::ly, oomph::DocInfo::number(), Mesh_Parameters::ny, oomph::oomph_info, Global_Physical_Variables::P_up, Global_Physical_Variables::Prescribed_y, problem, Global_Physical_Variables::Re, Global_Physical_Variables::ReSt, oomph::DocInfo::set_directory(), Flag_definition::setup(), Control_Flags::Steady_run, and Control_Flags::Validation_run.