axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks_ref.cc File Reference

#include "generic.h"
#include "navier_stokes.h"
#include "refineable_linearised_axisym_navier_stokes_elements.h"
#include "refineable_linearised_axisym_navier_stokes_elements.cc"
#include "multi_domain_linearised_axisym_navier_stokes_elements.h"
#include "linearised_axisym_navier_stokes_elements.h"
#include "linearised_axisym_navier_stokes_elements.cc"
#include "meshes/rectangular_quadmesh.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. More...

Variables
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.

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

{
 // Number of elements in radial (r) direction
 const unsigned base_n_r = 4;
 const unsigned perturbed_n_r = 4;
 
 // Number of elements in axial (z) direction
 const unsigned base_n_z = 4;
 const unsigned perturbed_n_z = 4;
 
 // 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 up labels for output
 DocInfo doc_info;
 
 // Determine length of timestep
 const double dt = 0.01;
 
 // Set number of timesteps to perform during a single iteration of the
 // power method
 const unsigned nstep = 1;
 
 // Cache tolerance for power method convergence
 // Note that we are using a very large tolerance for the purposes of
 // this demo driver so that it converges in a very small number of steps
 const double power_method_tolerance = 1e-4;
 
 // Set maximum number of iterations of the power method
 const unsigned max_iter = 1000;
 
 // For the purposes of code validation we will solve this problem twice.
 // The first time we shall use Taylor-Hood elements and the second time
 // Crouzeix-Raviart elements.
 
 // --------------------------------------------------------------------
 // Taylor-Hood elements
 // --------------------------------------------------------------------
 {
  // Build stability problem (this creates base and perturbed state problems)
  StabilityProblem<RefineableAxisymmetricQTaylorHoodElement,
   RefineableLinearisedAxisymmetricQTaylorHoodMultiDomainElement>
   problem_TH(base_n_r,base_n_z,perturbed_n_r,perturbed_n_z,domain_height);
  
  // Output directory
  doc_info.set_directory("RESLT_TH");
 
  // Create and initialise trace file
  ofstream trace;
  char filename[256];
  sprintf(filename,"%s/trace.dat",
          doc_info.directory().c_str());
  trace.open(filename);
  trace << "Re, dominating eigenvalue, n_power_method_iterations" << std::endl;
  
  // Refine both base and perturbed state problems uniformly
  problem_TH.base_state_problem_pt()->refine_uniformly();
  problem_TH.perturbed_state_problem_pt()->refine_uniformly();
  
  // Set up vector of elements to be refined
  const unsigned temp_array[22]
   = {0,1,4,5,8,9,12,13,15,29,31,45,47,50,51,54,55,58,59,61,62,63};
  Vector<unsigned> els_to_refine;
  for(unsigned i=0;i<22;i++) { els_to_refine.push_back(temp_array[i]); }
  
  // Refine selected elements
  problem_TH.base_state_problem_pt()->refine_selected_elements(els_to_refine);
  problem_TH.perturbed_state_problem_pt()
   ->refine_selected_elements(els_to_refine);
  
  // Initialise timestep for perturbed state problem
  // (also sets weights for all timesteppers)
  problem_TH.perturbed_state_problem_pt()->initialise_dt(dt);
  
  // Set initial conditions
  problem_TH.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_TH.perturbed_state_problem_pt()->get_dofs(eigenvector);
  
  // Initialise counter for solutions
  doc_info.number()=0;
  
  // Set base state problem's boundary conditions
  problem_TH.base_state_problem_pt()->set_boundary_conditions();
  
  // Perform a steady base flow solve
  problem_TH.base_state_problem_pt()->steady_newton_solve();
 
  // Doc the base flow solution
  problem_TH.base_state_problem_pt()->doc_solution(doc_info);
  
  // 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.
  problem_TH.perturbed_state_problem_pt()->enable_jacobian_reuse();
  
  // Perform power method and store number of iterations
  const unsigned n_power_method_iterations =
   problem_TH.perform_power_method(dt,nstep,doc_info,power_method_tolerance,
                                   max_iter,eigenvalue,eigenvector);
  
  cout << "\nDominating eigenvalue is " << eigenvalue << endl;
  
  // Document in trace file
  trace << GlobalPhysicalVariables::Re << " "
        << eigenvalue << " "
        << n_power_method_iterations << std::endl;
 }
 
 
 // --------------------------------------------------------------------
 // Crouzeix-Raviart elements
 // --------------------------------------------------------------------
 {
  // Build stability problem (this creates base and perturbed state problems)
  StabilityProblem<RefineableAxisymmetricQCrouzeixRaviartElement,
   RefineableLinearisedAxisymmetricQCrouzeixRaviartMultiDomainElement>
   problem_CR(base_n_r,base_n_z,perturbed_n_r,perturbed_n_z,domain_height);
  
  // Output directory
  doc_info.set_directory("RESLT_CR");
 
  // Create and initialise trace file
  ofstream trace;
  char filename[256];
  sprintf(filename,"%s/trace.dat",
          doc_info.directory().c_str());
  trace.open(filename);
  trace << "Re, dominating eigenvalue, n_power_method_iterations" << std::endl;
  
  // Refine both base and perturbed state problems uniformly
  problem_CR.base_state_problem_pt()->refine_uniformly();
  problem_CR.perturbed_state_problem_pt()->refine_uniformly();
  
  // Set up vector of elements to be refined
  const unsigned temp_array[22]
   = {0,1,4,5,8,9,12,13,15,29,31,45,47,50,51,54,55,58,59,61,62,63};
  Vector<unsigned> els_to_refine;
  for(unsigned i=0;i<22;i++) { els_to_refine.push_back(temp_array[i]); }
  
  // Refine selected elements
  problem_CR.base_state_problem_pt()->refine_selected_elements(els_to_refine);
  problem_CR.perturbed_state_problem_pt()
   ->refine_selected_elements(els_to_refine);
  
  // Initialise timestep for perturbed state problem
  // (also sets weights for all timesteppers)
  problem_CR.perturbed_state_problem_pt()->initialise_dt(dt);
  
  // Set initial conditions
  problem_CR.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_CR.perturbed_state_problem_pt()->get_dofs(eigenvector);
  
  // Initialise counter for solutions
  doc_info.number()=0;
  
  // Set base state problem's boundary conditions
  problem_CR.base_state_problem_pt()->set_boundary_conditions();
  
  // Perform a steady base flow solve
  problem_CR.base_state_problem_pt()->steady_newton_solve();
 
  // Doc the base flow solution
  problem_CR.base_state_problem_pt()->doc_solution(doc_info);
  
  // 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.
  problem_CR.perturbed_state_problem_pt()->enable_jacobian_reuse();
  
  // Perform power method and store number of iterations
  const unsigned n_power_method_iterations =
   problem_CR.perform_power_method(dt,nstep,doc_info,power_method_tolerance,
                                   max_iter,eigenvalue,eigenvector);
  
  cout << "\nDominating eigenvalue is " << eigenvalue << endl;
  
  // Document in trace file
  trace << GlobalPhysicalVariables::Re << " "
        << eigenvalue << " "
        << n_power_method_iterations << std::endl;
 }
 
} // End of main

References oomph::DocInfo::directory(), e(), MergeRestartFiles::filename, GlobalPhysicalVariables::G, GlobalPhysicalVariables::Gamma, i, oomph::DocInfo::number(), GlobalPhysicalVariables::Re, and oomph::DocInfo::set_directory().