#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 Tuckerman counter-rotating lids problem (no power method) More...
Variables
double	GlobalPhysicalVariables::St = 1.0
	Strouhal number. More...

double	GlobalPhysicalVariables::InvFr = 0.0
	Inverse of Froude number. More...

double	GlobalPhysicalVariables::Re_current

double	GlobalPhysicalVariables::ReSt_current

double	GlobalPhysicalVariables::ReInvFr_current

double	GlobalPhysicalVariables::Re_lower_limit = 300
	Loop information. More...

double	GlobalPhysicalVariables::Re_upper_limit = 302

unsigned	GlobalPhysicalVariables::Nparam_steps = 5

Vector< double >	GlobalPhysicalVariables::G (3)
	Direction of gravity. More...
Function Documentation

◆ main()

int main	(	int argc	,
		char *	argv[]
	)
Driver for Tuckerman counter-rotating lids problem (no power method)
////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////
 {
  // Store command line arguments
  CommandLineArgs::setup(argc,argv);
  
  // Set duration of timestep
  const double dt = 0.005;
  
  // Set number of elements in radial (r) direction
  const unsigned base_n_r = 25;
  const unsigned perturbed_n_r = 25;
  
  // Set number of elements in axial (z) direction
  const unsigned base_n_z = 25;
  const unsigned perturbed_n_z = 25;
  
  // Determine height of computational domain (this is just the aspect
  // ratio since we take the cylinder radius to always be one)
  const double domain_height = GlobalPhysicalVariables::Gamma;
  
  // Set direction of gravity (vertically downwards)
  GlobalPhysicalVariables::G[0] = 0.0;
  GlobalPhysicalVariables::G[1] = -1.0;
  GlobalPhysicalVariables::G[2] = 0.0;
  
  // Set maximum number of iterations of the power method
  const unsigned max_iter = 10000;
  
  // Set tolerance for power method convergence (this is usually machine
  // precision according to Bai, Demmel, Dongarra, Ruhe & Van der Vorst)
  double power_method_tolerance = 1e-8;
  
  // Cache upper and lower limits and number of steps for loop over Re
  double lower = GlobalPhysicalVariables::Re_lower_limit;
  double upper = GlobalPhysicalVariables::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 and use a lower tolerance
  if(CommandLineArgs::Argc>1)
   {
    n_param_steps = 1;
    power_method_tolerance = 1e-5;
   }
  
  // -----------------------------------------
  // 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,domain_height);
   
   // 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_k%i.dat",
           doc_info.directory().c_str(),
           GlobalPhysicalVariables::k);
   global_trace.open(filename);
   global_trace << "Re, dominating eigenvalue, "
                << "n_base_flow_continuation_steps, "
                << "n_power_method_iterations" << std::endl;
   
   // Begin loop over Reynolds number
   for(unsigned param=0;param<n_param_steps;param++)
    {
     // Set Re_current to the correct value
     if(lower==upper || n_param_steps<=1)
      {
       GlobalPhysicalVariables::Re_current = lower;
      }
     else
      {
       GlobalPhysicalVariables::Re_current =
        lower + (((upper-lower)/(n_param_steps-1))*param);
      }
     
     // Set ReSt and ReInvFr to the correct values
     GlobalPhysicalVariables::ReSt_current = 
      GlobalPhysicalVariables::Re_current*GlobalPhysicalVariables::St;
     GlobalPhysicalVariables::ReInvFr_current = 
      GlobalPhysicalVariables::Re_current*GlobalPhysicalVariables::InvFr;
     
     cout << "\n====================================================" << endl;
     cout << "Beginning parameter run " << param+1 << " of "
          << n_param_steps << ": Re = "
          << GlobalPhysicalVariables::Re_current << endl;
     cout << "====================================================" << endl;
     
     // Pass updated Reynolds number (and ReSt and ReInvFr) to elements
     problem.pass_updated_nondim_parameters_to_elements();
     
     // 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()++;
     
     // Find steady base flow
     // ---------------------
  
     // Determine number of steps to take when ramping up to actual
     // Reynolds number
     const unsigned n_base_flow_continuation_steps = 1;
     
     // If we're on the first parameter, need to get to steady state by
     // continuation
     if(param==0)
      {
       cout << "\nEntering steady base flow loop..." << endl;
       
       // Store the current parameter value
       const double Re_current_backup = GlobalPhysicalVariables::Re_current;
       
       // Loop over continuation steps
       for(unsigned i=0;i<=n_base_flow_continuation_steps;i++)
        {
         // Update current parameters
         GlobalPhysicalVariables::Re_current =
          i*(Re_current_backup/n_base_flow_continuation_steps);
         GlobalPhysicalVariables::ReSt_current = 
          GlobalPhysicalVariables::Re_current*GlobalPhysicalVariables::St;
         GlobalPhysicalVariables::ReInvFr_current = 
          GlobalPhysicalVariables::Re_current*GlobalPhysicalVariables::InvFr;
         
         cout << "   - Setting GlobalPhysicalVariables::Re_current = "
              << GlobalPhysicalVariables::Re_current << endl;
         
         // Pass updated Reynolds number (and ReSt and ReInvFr) to elements
         problem.pass_updated_nondim_parameters_to_elements();
         
         // Perform a steady base flow solve
         problem.base_flow_steady_newton_solve();
        }
      }
     // Otherwise, just require one steady base flow solve
     else { problem.base_flow_steady_newton_solve(); }
     
     // Enable Jacobian reuse in the perturbed state problem. Note that this
     // can only be done since the problem is linear and the base state is
     // steady, and so the Jacobian will be the same at each timestep.
     // Calling this function also resets the
     // Problem::Jacobian_has_been_computed flag to "false", so that during
     // the first solve it will be recomputed. Note that this is the reason
     // that we call this function here rather than in the problem constructor.
     // Were we to call it in the constructor then the Jacobian would not be
     // recomputed when we changed the problem parameters (Reynolds number etc.)
     problem.perturbed_state_problem_pt()->enable_jacobian_reuse();
     
     // Perform power method and store number of iterations
     const unsigned n_power_method_iterations =
      problem.perform_power_method(dt,doc_info,power_method_tolerance,
                                   max_iter,eigenvalue,eigenvector);
     
     cout << "\nDominating eigenvalue is " << eigenvalue << endl;
     
     // The corresponding dominating eigenvector is "input"
     
     // Document in the global trace file
     global_trace << GlobalPhysicalVariables::Re_current << " ";
     ios_base::fmtflags flags =  cout.flags(); // Save old flags
     global_trace.precision(14); // ten decimal places
     global_trace << eigenvalue << " ";
     global_trace.flags(flags);  // Set the flags to the way they were
     global_trace << n_base_flow_continuation_steps << " "
                  << 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,domain_height);
   
   // 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_k%i.dat",
           doc_info.directory().c_str(),
           GlobalPhysicalVariables::k);
   global_trace.open(filename);
   global_trace << "Re, dominating eigenvalue, "
                << "n_base_flow_continuation_steps, "
                << "n_power_method_iterations" << std::endl;
   
   // Begin loop over Reynolds number
   for(unsigned param=0;param<n_param_steps;param++)
    {
     // Set Re_current to the correct value
     if(lower==upper || n_param_steps<=1)
      {
       GlobalPhysicalVariables::Re_current = lower;
      }
     else
      {
       GlobalPhysicalVariables::Re_current =
        lower + (((upper-lower)/(n_param_steps-1))*param);
      }
     
     // Set ReSt and ReInvFr to the correct values
     GlobalPhysicalVariables::ReSt_current = 
      GlobalPhysicalVariables::Re_current*GlobalPhysicalVariables::St;
     GlobalPhysicalVariables::ReInvFr_current = 
      GlobalPhysicalVariables::Re_current*GlobalPhysicalVariables::InvFr;
     
     cout << "\n====================================================" << endl;
     cout << "Beginning parameter run " << param+1 << " of "
          << n_param_steps << ": Re = "
          << GlobalPhysicalVariables::Re_current << endl;
     cout << "====================================================" << endl;
     
     // Pass updated Reynolds number (and ReSt and ReInvFr) to elements
     problem.pass_updated_nondim_parameters_to_elements();
     
     // 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()++;
     
     // Find steady base flow
     // ---------------------
  
     // Determine number of steps to take when ramping up to actual
     // Reynolds number
     const unsigned n_base_flow_continuation_steps = 1;
     
     // If we're on the first parameter, need to get to steady state by
     // continuation
     if(param==0)
      {
       cout << "\nEntering steady base flow loop..." << endl;
       
       // Store the current parameter value
       const double Re_current_backup = GlobalPhysicalVariables::Re_current;
       
       // Loop over continuation steps
       for(unsigned i=0;i<=n_base_flow_continuation_steps;i++)
        {
         // Update current parameters
         GlobalPhysicalVariables::Re_current =
          i*(Re_current_backup/n_base_flow_continuation_steps);
         GlobalPhysicalVariables::ReSt_current = 
          GlobalPhysicalVariables::Re_current*GlobalPhysicalVariables::St;
         GlobalPhysicalVariables::ReInvFr_current = 
          GlobalPhysicalVariables::Re_current*GlobalPhysicalVariables::InvFr;
         
         cout << "   - Setting GlobalPhysicalVariables::Re_current = "
              << GlobalPhysicalVariables::Re_current << endl;
         
         // Pass updated Reynolds number (and ReSt and ReInvFr) to elements
         problem.pass_updated_nondim_parameters_to_elements();
         
         // Perform a steady base flow solve
         problem.base_flow_steady_newton_solve();
        }
      }
     // Otherwise, just require one steady base flow solve
     else { problem.base_flow_steady_newton_solve(); }
     
     // Enable Jacobian reuse in the perturbed state problem. Note that this
     // can only be done since the problem is linear and the base state is
     // steady, and so the Jacobian will be the same at each timestep.
     // Calling this function also resets the
     // Problem::Jacobian_has_been_computed flag to "false", so that during
     // the first solve it will be recomputed. Note that this is the reason
     // that we call this function here rather than in the problem constructor.
     // Were we to call it in the constructor then the Jacobian would not be
     // recomputed when we changed the problem parameters (Reynolds number etc.)
     problem.perturbed_state_problem_pt()->enable_jacobian_reuse();
     
     // Perform power method and store number of iterations
     const unsigned n_power_method_iterations =
      problem.perform_power_method(dt,doc_info,power_method_tolerance,
                                   max_iter,eigenvalue,eigenvector);
     
     cout << "\nDominating eigenvalue is " << eigenvalue << endl;
     
     // The corresponding dominating eigenvector is "input"
     
     // Document in the global trace file
     global_trace << GlobalPhysicalVariables::Re_current << " ";
     ios_base::fmtflags flags =  cout.flags(); // Save old flags
     global_trace.precision(14); // ten decimal places
     global_trace << eigenvalue << " ";
     global_trace.flags(flags);  // Set the flags to the way they were
     global_trace << n_base_flow_continuation_steps << " "
                  << 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
Classes

Namespaces

Functions

Variables

Function Documentation

◆ main()
