time_periodic_taylor_couette.cc File Reference

#include "generic.h"
#include "navier_stokes.h"
#include "linearised_axisym_navier_stokes.h"
#include "meshes/rectangular_quadmesh.h"
#include "../self_starting_BDF2_timestepper.h"

Classes
class	BaseStateProblem< BASE_ELEMENT >
	Base state problem class for Tuckerman counter-rotating lids problem. More...
class	PerturbedStateProblem< BASE_ELEMENT, PERTURBED_ELEMENT >
class	StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >

Namespaces
	GlobalPhysicalVariables
	Namespace for physical parameters.

Functions
int	main (int argc, char *argv[])
	Driver for time-periodic Taylor–Couette problem. More...

Variables
double	GlobalPhysicalVariables::MMC_Re_current
double	GlobalPhysicalVariables::MMC_Re_lower_limit = 20.93
	Loop information. More...
double	GlobalPhysicalVariables::MMC_Re_upper_limit = 20.94
double	GlobalPhysicalVariables::AngularFrequency = 10.58
	Angular frequency of outer cylinder rotation (omega = 2*gamma^2) More...
double	GlobalPhysicalVariables::Epsilon = 0.9
	Dimensionless modulation amplitude (epsilon) More...
double	GlobalPhysicalVariables::Eta = 0.88
	Radius ratio (radius of inner cylinder / radius of outer cylinder) More...
double	GlobalPhysicalVariables::a = 3.2
	Wavenumber of azimuthal waves e^iaz. More...
Vector< double >	GlobalPhysicalVariables::G (3)
	Direction of gravity. More...

Function Documentation

main()

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

Driver for time-periodic Taylor–Couette problem.

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

{
 // Store command line arguments
 CommandLineArgs::setup(argc,argv);
 
 // Set number of elements in radial (r) direction
 const unsigned base_n_r = 16;
 const unsigned perturbed_n_r = 16;
 
 // Set number of elements in axial (z) direction
 const unsigned base_n_z = 16;
 const unsigned perturbed_n_z = 16;
 
 // Cache radius ratio
 const double eta = GlobalPhysicalVariables::Eta;
 
 // Determine tadius of inner cylinder
 const double radius_inner_cylinder = eta/(1.0-eta);
 
 // Determine radius of outer cylinder
 const double radius_outer_cylinder = 1.0/(1.0-eta);
 
 // Length in axial (z) direction (set to fit in half a wavelength)
 const double l_z = MathematicalConstants::Pi/GlobalPhysicalVariables::a;
 
 // Set direction of gravity (vertically downwards)
 GlobalPhysicalVariables::G[0] = 0.0;
 GlobalPhysicalVariables::G[1] = -1.0;
 GlobalPhysicalVariables::G[2] = 0.0;
 
 // Determine period of the base state problem
 const double T = 
  (2*MathematicalConstants::Pi)/GlobalPhysicalVariables::AngularFrequency;
 
 // Set number of timesteps per period
 unsigned nstep = 40;
 
 // Set maximum number of iterations of the power method
 unsigned max_iter = 1000;
 
 // Set tolerance for power method convergence (this is usually machine
 // precision according to Bai, Demmel, Dongarra, Ruhe & Van der Vorst)
 const double power_method_tolerance = 1e-8;
 
 // Set tolerance for L2 norm of difference in dof vectors from one
 // period to the next for defining the point at which the base flow
 // has settled down to a periodic state.
 double base_flow_norm_tolerance = 1e-6;
 
 // Cache upper and lower limits and number of steps for loop over MMC_Re
 double lower = GlobalPhysicalVariables::MMC_Re_lower_limit;
 double upper = GlobalPhysicalVariables::MMC_Re_upper_limit;
 double n_param_steps = GlobalPhysicalVariables::Nparam_steps;
 
 // If we are doing a validation run, only compute one value of the
 // Reynolds number, use a lower number of timesteps per period,
 // accept a lower base flow convergence tolerance and insist that
 // only one iteration of the power method is performed
 if(CommandLineArgs::Argc>1)
  {
   n_param_steps = 1;
   nstep = 20;
   base_flow_norm_tolerance = 1e-2;
   max_iter = 1;
  }
 
 // Check that nstep is exactly divisible by both 2 and 5
 // (This is necessary for the self_starting_BDF2 timestepper)
 if(nstep%2 != 0 || nstep%5 !=0)
  {
   cout << "\nThe self-starting BDF2 timestepper requires that the number"
        << "\nof timesteps be exactly divisible by both 2 and 5. This is"
        << "\nnot the case for nstep = " << nstep << ". Exiting...\n"
        << endl;
   exit(2);
  }
 
 // Calculate length of timestep
 const double dt = T/nstep;
 
 // -----------------------------------------
 // LinearisedAxisymmetricQTaylorHoodElements
 // -----------------------------------------
 {
  cout << "\nDoing LinearisedAxisymmetricQTaylorHoodElements\n" << endl;
 
  // Build stability problem (this creates base and perturbed state problems)
  StabilityProblem<AxisymmetricQTaylorHoodElement,
   LinearisedAxisymmetricQTaylorHoodMultiDomainElement>
   problem(base_n_r,base_n_z,perturbed_n_r,perturbed_n_z,
           radius_inner_cylinder,radius_outer_cylinder,l_z);
  
  // Set up labels for output
  DocInfo doc_info;
 
  // Output directory
  doc_info.set_directory("RESLT_TH");
  
  // Initialise timestep (also sets weights for all timesteppers)
  problem.initialise_dt(dt);
  
  // Set initial conditions
  problem.set_initial_condition();
  
  // Set up storage for dominant eigenvector and eigenvalue
  DoubleVector eigenvector;
  double eigenvalue = 0.0;
  
  // Get initial guess for eigenvector (these are the initial conditions
  // set up in PerturbedStateProblem::set_initial_conditions())
  problem.get_perturbed_state_problem_dofs(eigenvector);
  
  // Create and initialise global trace file
  ofstream global_trace;
  char filename[256];
  sprintf(filename,"%s/global_trace_epsilon%2.1f.dat",
          doc_info.directory().c_str(),
          GlobalPhysicalVariables::Epsilon);
  global_trace.open(filename);
  global_trace << "MMC_Re, dominating eigenvalue, "
               << "n_periods_to_set_up_periodic_base_flow, "
               << "n_power_method_iterations" << std::endl;
  
  // Begin loop over Reynolds number
  for(unsigned param=0;param<n_param_steps;param++)
   {
    // Set MMC_Re_current to the correct value
    if(lower==upper || n_param_steps<=1)
     {
      GlobalPhysicalVariables::MMC_Re_current = lower;
     }
    else
     {
      GlobalPhysicalVariables::MMC_Re_current =
       lower + (((upper-lower)/(n_param_steps-1))*param);
     }
    
    cout << "\n====================================================" << endl;
    cout << "Beginning parameter run " << param+1 << " of "
         << n_param_steps << ": MMC_Re = "
         << GlobalPhysicalVariables::MMC_Re_current << endl;
    cout << "====================================================" << endl;
    
    // Reset the global time (of both problems!) to zero
    problem.reset_global_time_to_zero();
    
    // Initialise counter for solutions
    doc_info.number()=0;
    
    // Create trace files
    problem.create_trace_files(doc_info);
    
    // Initialise trace files
    problem.initialise_trace_files();
    
    // Output initial conditions
    problem.doc_solution(doc_info,false,false);
    
    // Increment counter for solutions 
    doc_info.number()++;
    
    // Set up periodic base flow and record how many period it takes
    const unsigned n_periods_to_set_up_periodic_base_flow
     = problem.set_up_periodic_base_flow(dt,nstep,doc_info,
                                         base_flow_norm_tolerance);
    
    // Set up storage for base state dofs over one period
    const unsigned n_entries_required = (27*(nstep/10))-3;
    Vector<DoubleVector> base_state_dofs_over_period(n_entries_required);
    
    // Perform power method and store number of iterations
    const unsigned n_power_method_iterations =
     problem.perform_power_method(dt,nstep,doc_info,power_method_tolerance,
                                  max_iter,eigenvalue,eigenvector,
                                  base_state_dofs_over_period);
    
    cout << "\nDominating eigenvalue is " << eigenvalue << endl;
    
    // The corresponding dominating eigenvector is "input"
    
    // Document in the global trace file
    global_trace << GlobalPhysicalVariables::MMC_Re_current << " "
                 << eigenvalue << " "
                 << n_periods_to_set_up_periodic_base_flow << " "
                 << n_power_method_iterations << std::endl;
    
    // The initial guess for the next power method will be the eigenvector
    // calculated by this power method
    
    // Clear and close trace files
    problem.close_trace_files();
    
   } // End loop over Reynolds number
 }
 
 // ----------------------------------------------
 // LinearisedAxisymmetricQCrouzeixRaviartElements
 // ----------------------------------------------
 {
  cout << "\nDoing LinearisedAxisymmetricQCrouzeixRaviartElements\n" << endl;
 
  // Build stability problem (this creates base and perturbed state problems)
  StabilityProblem<AxisymmetricQCrouzeixRaviartElement,
   LinearisedAxisymmetricQCrouzeixRaviartMultiDomainElement>
   problem(base_n_r,base_n_z,perturbed_n_r,perturbed_n_z,
           radius_inner_cylinder,radius_outer_cylinder,l_z);
  
  // Set up labels for output
  DocInfo doc_info;
 
  // Output directory
  doc_info.set_directory("RESLT_CR");
  
  // Initialise timestep (also sets weights for all timesteppers)
  problem.initialise_dt(dt);
  
  // Set initial conditions
  problem.set_initial_condition();
  
  // Set up storage for dominant eigenvector and eigenvalue
  DoubleVector eigenvector;
  double eigenvalue = 0.0;
  
  // Get initial guess for eigenvector (these are the initial conditions
  // set up in PerturbedStateProblem::set_initial_conditions())
  problem.get_perturbed_state_problem_dofs(eigenvector);
  
  // Create and initialise global trace file
  ofstream global_trace;
  char filename[256];
  sprintf(filename,"%s/global_trace_epsilon%2.1f.dat",
          doc_info.directory().c_str(),
          GlobalPhysicalVariables::Epsilon);
  global_trace.open(filename);
  global_trace << "MMC_Re, dominating eigenvalue, "
               << "n_periods_to_set_up_periodic_base_flow, "
               << "n_power_method_iterations" << std::endl;
  
  // Begin loop over Reynolds number
  for(unsigned param=0;param<n_param_steps;param++)
   {
    // Set MMC_Re_current to the correct value
    if(lower==upper || n_param_steps<=1)
     {
      GlobalPhysicalVariables::MMC_Re_current = lower;
     }
    else
     {
      GlobalPhysicalVariables::MMC_Re_current =
       lower + (((upper-lower)/(n_param_steps-1))*param);
     }
    
    cout << "\n====================================================" << endl;
    cout << "Beginning parameter run " << param+1 << " of "
         << n_param_steps << ": MMC_Re = "
         << GlobalPhysicalVariables::MMC_Re_current << endl;
    cout << "====================================================" << endl;
    
    // Reset the global time (of both problems!) to zero
    problem.reset_global_time_to_zero();
    
    // Initialise counter for solutions
    doc_info.number()=0;
    
    // Create trace files
    problem.create_trace_files(doc_info);
    
    // Initialise trace files
    problem.initialise_trace_files();
    
    // Output initial conditions
    problem.doc_solution(doc_info,false,false);
    
    // Increment counter for solutions 
    doc_info.number()++;
    
    // Set up periodic base flow and record how many period it takes
    const unsigned n_periods_to_set_up_periodic_base_flow
     = problem.set_up_periodic_base_flow(dt,nstep,doc_info,
                                         base_flow_norm_tolerance);
    
    // Set up storage for base state dofs over one period
    const unsigned n_entries_required = (27*(nstep/10))-3;
    Vector<DoubleVector> base_state_dofs_over_period(n_entries_required);
    
    // Perform power method and store number of iterations
    const unsigned n_power_method_iterations =
     problem.perform_power_method(dt,nstep,doc_info,power_method_tolerance,
                                  max_iter,eigenvalue,eigenvector,
                                  base_state_dofs_over_period);
    
    cout << "\nDominating eigenvalue is " << eigenvalue << endl;
    
    // The corresponding dominating eigenvector is "input"
    
    // Document in the global trace file
    global_trace << GlobalPhysicalVariables::MMC_Re_current << " "
                 << eigenvalue << " "
                 << n_periods_to_set_up_periodic_base_flow << " "
                 << n_power_method_iterations << std::endl;
    
    // The initial guess for the next power method will be the eigenvector
    // calculated by this power method
    
    // Clear and close trace files
    problem.close_trace_files();
    
   } // End loop over Reynolds number
 }
 
} // End of main

References GlobalPhysicalVariables::a, GlobalPhysicalVariables::AngularFrequency, oomph::CommandLineArgs::Argc, oomph::DocInfo::directory(), e(), GlobalPhysicalVariables::Epsilon, TestSoln::eta, GlobalPhysicalVariables::Eta, MergeRestartFiles::filename, GlobalPhysicalVariables::G, helpers::lower(), GlobalPhysicalVariables::MMC_Re_current, GlobalPhysicalVariables::MMC_Re_lower_limit, GlobalPhysicalVariables::MMC_Re_upper_limit, GlobalPhysicalVariables::Nparam_steps, oomph::DocInfo::number(), BiharmonicTestFunctions2::Pi, problem, oomph::DocInfo::set_directory(), and Flag_definition::setup().