one_d_womersley.cc File Reference

#include <cstring>
#include "generic.h"
#include "navier_stokes.h"
#include "meshes/collapsible_channel_mesh.h"
#include "womersley.h"
#include "multi_physics.h"

Classes
class	OscillatingWall
class	CollapsibleChannelProblem< ELEMENT >
	Problem class. More...

Namespaces
	BL_Squash
	Global_Physical_Variables
	Global variables.

Functions
double	BL_Squash::squash_fct (const double &s)
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...
int	main (int argc, char *argv[])

Function Documentation

main()

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

Driver code for an unsteady adaptive collapsible channel problem with prescribed wall motion. Presence of command line arguments indicates validation run with coarse resolution and small number of timesteps.

{
 // Store command line arguments
 CommandLineArgs::setup(argc,argv);
 
 //Uncomment for parallel stuff
 /*#ifdef OOMPH_HAS_MPI
 // Initialise MPI
 MPI_Helpers::init(argc,argv);
 
 // Swtich off output modifier
 oomph_info.output_modifier_pt() = &default_output_modifier;
 
 // switch off oomph_info output for all processors but rank 0
 if (MPI_Helpers::My_rank!=0)
  {
   oomph_info.stream_pt() = &oomph_nullstream;
   OomphLibWarning::set_stream_pt(&oomph_nullstream);
   OomphLibError::set_stream_pt(&oomph_nullstream);
  }
 else
  {
   oomph_info << "\n\n=====================================================\n";
   oomph_info << "Number of processors: " 
              << MPI_Helpers::Nproc << "\n";
  }
  #endif */
 
 // Set default values
 string directory_for_data = "RESLT";
 unsigned outflow = 0;
 Global_Physical_Variables::Re = 10;
 unsigned solver=1;
 unsigned p_solver=2;
 unsigned f_solver=2;
 unsigned refinement=0;
 unsigned nsteps=200;
 double amplitude=1.0e-2;
 double period=0.45;
 bool validation_run=false;
 
 //Set Trilinos default value
#ifdef OOMPH_HAS_TRILINOS 
 unsigned f_ml_settings=1;
#endif
 
 // Set Hypre default values
#ifdef OOMPH_HAS_HYPRE
 unsigned f_bamg_smoother=1;
 double f_bamg_damping=1.0;
 double f_bamg_strength=0.25;
#endif
 
 
 
 // Parse command line
 int arg_index = 0;
 unsigned print_help = 0;
 
 while (arg_index < argc)
  {
   if ( strcmp(argv[arg_index], "-dir") == 0 )
    {
     arg_index++;
     string directory_number = argv[arg_index++];
     directory_for_data += directory_number;
    }
   else if ( strcmp(argv[arg_index], "-outflow") == 0 )
    {
     arg_index++;
     outflow = atoi(argv[arg_index++]);
    }
   else if ( strcmp(argv[arg_index], "-Re") == 0 )
    {
     arg_index++;
     Global_Physical_Variables::Re = atof(argv[arg_index++]);
    }
   else if ( strcmp(argv[arg_index], "-solver") == 0 )
    {
     arg_index++;
     solver = atoi(argv[arg_index++]);
    }
   else if ( strcmp(argv[arg_index], "-p_solver") == 0 )
    {
     arg_index++;
     p_solver = atoi(argv[arg_index++]);
    }
   else if ( strcmp(argv[arg_index], "-f_solver") == 0 )
    {
     arg_index++;
     f_solver = atoi(argv[arg_index++]);
    }
#ifdef OOMPH_HAS_HYPRE
   else if ( strcmp(argv[arg_index], "-f_bamg_smoother") == 0 )
    {
     arg_index++;
     f_bamg_smoother = atoi(argv[arg_index++]);
    }
   else if ( strcmp(argv[arg_index], "-f_bamg_damping") == 0 )
    {
     arg_index++;
     f_bamg_damping = atof(argv[arg_index++]);
    }  
   else if ( strcmp(argv[arg_index], "-f_bamg_strength") == 0 )
    {
     arg_index++;
     f_bamg_strength = atof(argv[arg_index++]);
    }
#endif
#ifdef OOMPH_HAS_TRILINOS
   else if ( strcmp(argv[arg_index], "-f_ml_settings") == 0 )
    {
     arg_index++;
     f_ml_settings = atoi(argv[arg_index++]);
    }
#endif
   else if ( strcmp(argv[arg_index], "-refine") == 0 )
    {
     arg_index++;
     refinement = atoi(argv[arg_index++]);
    }
   else if ( strcmp(argv[arg_index], "-amplitude") == 0 )
    {
     arg_index++;
     amplitude = atof(argv[arg_index++]);
    }
   else if ( strcmp(argv[arg_index], "-nsteps") == 0 )
    {
     arg_index++;
     nsteps = atoi(argv[arg_index++]);
    }
   else if ( strcmp(argv[arg_index], "-period") == 0 )
    {
     arg_index++;
     period = atof(argv[arg_index++]);
    }
   else if ( strcmp(argv[arg_index], "-validation_run") == 0 )
    {
     arg_index++;
     validation_run=true;
    }
   else if ( strcmp(argv[arg_index], "-help") == 0 )
    {
     print_help = 1;
     break;
    }
   else
    {
     arg_index++;
    }
  }
 
 if (print_help)
  {
   if (MPI_Helpers::communicator_pt()->my_rank()==0)
    {
     oomph_info << "\n\nOption flags:\n";
     oomph_info << "-dir <n>                  Data saved to /RESLTn\n";
     oomph_info << "-outflow <ID>             Outflow type:\n"
                << "                          ID=0 flux control\n"
                << "                          ID=1 pressure driven with impedance"
                << " with LSC preconditioner\n"
                << "                          ID=2 pressure driven with impedance"
                << " with block triangluar and LSC preconditioner\n"
                << "                          ID=3 prescribed outflow\n"
               
                << "                          otherwise pressure driven flow\n"; 
     oomph_info << "-Re <value>               Sets Re=value\n";
     oomph_info << "-nsteps <value>           Number of timesteps\n";
     oomph_info << "-solver <ID>              Linear solver type:\n"
                << "                          ID=0 GMRES\n"
                << "                          otherwise SuperLU is used\n";
     oomph_info << "-p_solver <ID>            P matrix solve type:\n"
                << "                          ID=0 BoomerAMG\n"
                << "                          ID=1 ML SA-AMG\n"
                << "                          otherwise SuperLU is used\n";
     oomph_info << "-f_solver <ID>            F matrix solve type\n"
                << "                          ID=0 BoomerAMG\n"
                << "                          ID=1 ML SA-AMG\n"
                << "                          otherwise SuperLU is used\n";
     oomph_info << "-f_bamg_smoother <ID>     BoomerAMG smoother for F matrix\n"
                << "                          ID=0 damped Jacobi\n"
                << "                          ID=1 Gauss-Seidel\n";
     oomph_info << "-f_bamg_damping <value>   BoomerAMG Jacobi damping for F matrix\n";
     oomph_info << "-f_bamg_strength <value>  BoomerAMG strength paramater F matrix\n";
     oomph_info << "-f_ml_settings <ID>       Settings for ML on F matrix\n"
                << "                          ID=0 NSSA\n"
                << "                          otherwise SA\n";
     oomph_info << "-refine <n>               Refines mesh n times\n";
     oomph_info << "-period <value>           Set period to value\n";
     oomph_info << "-amplitude <value>        Set amplitude to value\n";
     oomph_info << "-validation_run           Generate validation data for 1D Womersley\n";
    }
//#ifdef OOMPH_HAS_MPI 
//   // finalize MPI
//  MPI_Helpers::finalize();
//#endif
   return (0);
  }
 
 // Set validation values
 if (validation_run)
  {
   if (outflow==0)
    {
     directory_for_data = "RESLT_flux_control";
    }
   else if (outflow==1)
    {
     directory_for_data = "RESLT_impedance_tube";
    }
   else if (outflow==2)
    {
     directory_for_data = "RESLT_impedance_tube_with_flux_control";
    }
   Global_Physical_Variables::Re = 100;
   Global_Physical_Variables::ReSt = 100;
   solver=0;
   refinement=0;
   amplitude=0.5;
   period=10.0;
  }
 
 // Let St=1
 Global_Physical_Variables::ReSt = Global_Physical_Variables::Re;
 
 // Output settings
 oomph_info << "Storing data in /"
            << directory_for_data << "\n";
 oomph_info << "Outflow type: "
            << outflow << "\n";
 oomph_info << "Re="
            <<  Global_Physical_Variables::Re
            << "\n";
 if (solver==0)
  {
   oomph_info << "Using GMRES solver\n";
   if (p_solver==0)
    {
#ifdef OOMPH_HAS_HYPRE
     oomph_info << "Using BoomerAMG on P matrix\n";
#else
     oomph_info 
       << "Warning: Hypre not available. Using SuperLU on P matrix\n";
#endif
    }
   else if (p_solver==1)
    {
#ifdef OOMPH_HAS_TRILINOS
     oomph_info << "Using ML SA-AMG on P matrix\n";
#else
     oomph_info 
       << "Warning: Trilinos not available. Using SuperLU on P matrix\n";
#endif
    }
   else
    {
     oomph_info << "Using SuperLU on P matrix\n";
    }
   if (f_solver==0)
    {
#ifdef OOMPH_HAS_HYPRE
     oomph_info << "Using BoomerAMG on F matrix\n"
                << "smoother: " << f_bamg_smoother
                << "\ndamping: " << f_bamg_damping
                << "\nstrength: " << f_bamg_strength 
                << "\n";
#else
     oomph_info 
       << "Warning: Hypre not available. Using SuperLU on F matrix\n";
#endif
    }
   else if (f_solver==1)
    {
#ifdef OOMPH_HAS_TRILINOS
     oomph_info << "Using ML SA-AMG on F matrix ";
     if (f_ml_settings==0)
      {
       oomph_info << "with NSSA settings\n";
      }
     else
      {
       oomph_info << "with SA settings\n";
      }
#else
     oomph_info 
       << "Warning: Trilinos not available. Using SuperLU on F matrix\n";
#endif
    }
   else
    {
     oomph_info << "Using SuperLU on F matrix\n";
    }
  }
 else
  {
   oomph_info << "Using SuperLU solver\n";
  }
 
 oomph_info << "Refining mesh " << refinement << " times\n";
 
 // Number of elements in the domain
 unsigned nup=20;
 unsigned ncollapsible=8;
 unsigned ndown=40;
 unsigned ny=16;
 
 // Length of the domain
 double lup=5.0;
 double lcollapsible=1.0;
 double ldown=10.0;
 double ly=1.0;
 double limpedance=0.0;
 
 // Modify values if using impedance tube
 if ((outflow==1) || (outflow==2))
  {
   limpedance = 5.0;
   ldown-=limpedance;
   ndown=20;
  }
 
 // Modify values for validation run
 if (validation_run)
  {
   limpedance = 8.0;
   ldown=2.0;
   nup=10;
   ncollapsible=4;
   ny=8;
   ndown=4;
  }
 
 oomph_info << "Amplitude=" << amplitude
            << "\nPeriod=" << period << "\n";
 
 // Solver stuff
 GMRES<CRDoubleMatrix>* iterative_solver_pt=0;
 NavierStokesSchurComplementPreconditioner* ns_preconditioner_pt=0;
 
 // Don't build the problem with Crouzeix Raviart Elements if using the 
 // SchurComplement preconditioner!
 CollapsibleChannelProblem<RefineableQTaylorHoodElement<2> > 
  problem(nup, ncollapsible, ndown, ny, 
          lup, lcollapsible, ldown, ly, limpedance,
          amplitude, period, outflow);
 
 //problem.switch_off_wall_oscillations();
 
 // Refine problem
 for (unsigned i=0;i<refinement;i++)
  {
   problem.refine_uniformly();
  }
 
 if (solver==0)
  {
   // Set up the solver and pass to the problem
   iterative_solver_pt = new GMRES<CRDoubleMatrix>;
   //  iterative_solver_pt->set_preconditioner_RHS();
   iterative_solver_pt->max_iter() = 200;
   iterative_solver_pt->tolerance() = 1.0e-8;
   problem.linear_solver_pt() = iterative_solver_pt;
 
   // Set up the preconditioner
   NavierStokesSchurComplementPreconditioner* ns_preconditioner_pt = 0;
   FSIPreconditioner* fsi_preconditioner_pt = 0;
 
   if (outflow==2)
    {
     // Construct the FSI preconditioner
     fsi_preconditioner_pt = new FSIPreconditioner(&problem);
     fsi_preconditioner_pt->enable_doc_time();
 
     // Get a pointer to the preconditioner
     ns_preconditioner_pt = 
      dynamic_cast<NavierStokesSchurComplementPreconditioner*>
      (fsi_preconditioner_pt->navier_stokes_preconditioner_pt());
     
     // Set solid mesh with "true" to tolerate multiple element types.
     fsi_preconditioner_pt->
      set_wall_mesh(problem.outflow_impedance_master_mesh_pt(),true);
 
     // Set fluid mesh
     fsi_preconditioner_pt->
      set_navier_stokes_mesh
      (problem.create_mesh_for_navier_stokes_preconditioner());
 
     // Pass the FSI preconditioner to the solver
     iterative_solver_pt->preconditioner_pt() = fsi_preconditioner_pt;
    }
   else
    {
     // Construct the preconditioner and pass in the problem pointer.
     ns_preconditioner_pt = new NavierStokesSchurComplementPreconditioner
      (&problem);
 
     // Setup the fluid mesh
     ns_preconditioner_pt->set_navier_stokes_mesh
      (problem.create_mesh_for_navier_stokes_preconditioner());
     
     // Pass the preconditioner to the solver
     iterative_solver_pt->preconditioner_pt() = ns_preconditioner_pt;
    }
 
   // Settings for the LCS preconditioner
   ns_preconditioner_pt->set_navier_stokes_mesh(
    problem.create_mesh_for_navier_stokes_preconditioner());
   ns_preconditioner_pt->enable_doc_time();
   
   // Set up the P sub block preconditioner
#ifdef OOMPH_HAS_HYPRE
   if (p_solver==0)
    {
#ifdef OOMPH_HAS_HYPRE
     Preconditioner* p_preconditioner_pt = new HyprePreconditioner;
     HyprePreconditioner* hypre_preconditioner_pt = 
      static_cast<HyprePreconditioner*>(p_preconditioner_pt);
     Hypre_default_settings::
      set_defaults_for_2D_poisson_problem(hypre_preconditioner_pt);
     ns_preconditioner_pt->set_p_preconditioner(p_preconditioner_pt);
#endif
    }
#endif
 
#ifdef OOMPH_HAS_TRILINOS
   if (p_solver==1)
    {
#ifdef OOMPH_HAS_TRILINOS
     Preconditioner* p_preconditioner_pt = new TrilinosMLPreconditioner;
     ns_preconditioner_pt->set_p_preconditioner(p_preconditioner_pt);
#endif
    }
#endif
    
   // Set up the F sub block precondioner
#ifdef OOMPH_HAS_HYPRE
   if (f_solver==0)
    {
#ifdef OOMPH_HAS_HYPRE
     Preconditioner* f_preconditioner_pt = new HyprePreconditioner;
     HyprePreconditioner* hypre_preconditioner_pt = 
      static_cast<HyprePreconditioner*>(f_preconditioner_pt);
     hypre_preconditioner_pt->use_BoomerAMG();
     hypre_preconditioner_pt->set_amg_iterations(1);
     hypre_preconditioner_pt->amg_using_simple_smoothing();
     hypre_preconditioner_pt->amg_simple_smoother()=f_bamg_smoother;
     hypre_preconditioner_pt->amg_damping()=f_bamg_damping;
     hypre_preconditioner_pt->amg_strength()=f_bamg_strength;
     ns_preconditioner_pt->set_f_preconditioner(f_preconditioner_pt);
#endif 
   }
#endif
 
#ifdef OOMPH_HAS_TRILINOS    
   if (f_solver==1)
    {
#ifdef OOMPH_HAS_TRILINOS
     Preconditioner* f_preconditioner_pt = new TrilinosMLPreconditioner;
     TrilinosMLPreconditioner* trilinos_preconditioner_pt = 
      static_cast<TrilinosMLPreconditioner*>(f_preconditioner_pt);
     if (f_ml_settings==1)
      {
       trilinos_preconditioner_pt->set_NSSA_default_values();
      }
     else
      {
       trilinos_preconditioner_pt->set_SA_default_values();
      }
     ns_preconditioner_pt->set_f_preconditioner(f_preconditioner_pt);
#endif
    }
#endif
  }
 
 // run the problem
 problem.unsteady_run(directory_for_data, nsteps, validation_run);
 
 // delete stuff - remembering P and F matrix preconditioners are 
 // deleted in the Schur complement preconditioner
 delete iterative_solver_pt;
 iterative_solver_pt=0;
 delete ns_preconditioner_pt;
 ns_preconditioner_pt=0;
 
 
//#ifdef OOMPH_HAS_MPI 
// // finalize MPI
// MPI_Helpers::finalize();
//#endif
 
} //end of driver code

References oomph::HyprePreconditioner::amg_damping(), oomph::HyprePreconditioner::amg_simple_smoother(), oomph::HyprePreconditioner::amg_strength(), oomph::HyprePreconditioner::amg_using_simple_smoothing(), oomph::FSIPreconditioner::enable_doc_time(), oomph::NavierStokesSchurComplementPreconditioner::enable_doc_time(), i, Mesh_Parameters::ly, oomph::IterativeLinearSolver::max_iter(), oomph::FSIPreconditioner::navier_stokes_preconditioner_pt(), Mesh_Parameters::ny, oomph::oomph_info, oomph::IterativeLinearSolver::preconditioner_pt(), problem, Global_Physical_Variables::Re, Global_Physical_Variables::ReSt, oomph::HyprePreconditioner::set_amg_iterations(), oomph::Hypre_default_settings::set_defaults_for_2D_poisson_problem(), oomph::NavierStokesSchurComplementPreconditioner::set_f_preconditioner(), oomph::NavierStokesSchurComplementPreconditioner::set_navier_stokes_mesh(), oomph::TrilinosMLPreconditioner::set_NSSA_default_values(), oomph::NavierStokesSchurComplementPreconditioner::set_p_preconditioner(), oomph::TrilinosMLPreconditioner::set_SA_default_values(), Flag_definition::setup(), solver, oomph::IterativeLinearSolver::tolerance(), and oomph::HyprePreconditioner::use_BoomerAMG().