#include <algorithm>
#include "generic.h"
#include "polar_navier_stokes.h"
#include "meshes/rectangular_quadmesh.h"
#include "./refineable_r_mesh.h"
#include "./streamfunction_include.h"
#include "./jh_includes.h"

Classes
class	PolarNSProblem< ELEMENT >
	Driven cavity problem in rectangular domain. More...

Namespaces
	Global_Physical_Variables
	Global variables.

Functions
void	Global_Physical_Variables::traction_function (const double &time, const Vector< double > &x, Vector< double > &traction)
	Unused (but assigned) function to specify tractions. More...

int	main ()

Variables
double	Global_Physical_Variables::R_l =0.01

double	Global_Physical_Variables::R_r =1.

int	Global_Physical_Variables::xmesh =30

int	Global_Physical_Variables::ymesh =15

double	Global_Physical_Variables::Rstep_prestart =30.0

double	Global_Physical_Variables::Rmax_prestart =94.

double	Global_Physical_Variables::Rstep =0.1

double	Global_Physical_Variables::Rmax =100.

double	Global_Physical_Variables::epsilon =0.1

bool	Global_Physical_Variables::inlet_traction =false

double	Global_Physical_Variables::eta_inlet =1.0

bool	Global_Physical_Variables::outlet_traction =true

double	Global_Physical_Variables::eta_outlet =0.0

bool	Global_Physical_Variables::pinv =true

bool	Global_Physical_Variables::stokes =false

bool	Global_Physical_Variables::log_mesh =true

bool	Global_Physical_Variables::new_outlet_region =true

int	Global_Physical_Variables::uniform = 0

int	Global_Physical_Variables::adaptive = 0

Function Documentation

◆ main()

int main ( )

Driver for PolarNS test problem – test drive with two different types of element.

 { 
  // Set up doc info
  // ---------------
  DocInfo doc_info;
  doc_info.set_directory("RESLT");
  doc_info.number()=0; 
  // And again for eigenmode
  DocInfo eigenmode;
  eigenmode.set_directory("RESLT");
  eigenmode.number() = 100;
  // ---------------
  
  using namespace Global_Physical_Variables;
  Vector<double> shear_stress(2,0.0);
  
  // Build the problem
  PolarNSProblem<RefineablePolarCrouzeixRaviartElement > problem;
  //PolarNSProblem<RefineablePolarTaylorHoodElement > problem;
   
  cout << endl << "# xmesh: " << xmesh << " ymesh: " << ymesh << endl;
  
  //Increace the Reynolds number in large increments 
  for(;Re<Rmax_prestart;Re+=Rstep_prestart)
    {
      cout << endl << "Prestart solving for Re: " << Re << endl << endl; 
  
      problem.newton_solve(); 
      problem.doc_solution(doc_info);
      doc_info.number()+=1;
      //problem.output_streamfunction(doc_info,false);
    }
  
  exit(1);
  
  //Always start from same Reynolds number
  Re=Rmax_prestart;
  
  problem.newton_solve();
  problem.doc_solution(doc_info);
  doc_info.number()+=1;
  
  bool sym=true;
  
  //Storage for eigenvalues and eigenvectors 
  Vector<complex<double> > eigenvalues;
  Vector<DoubleVector> eigenvectors;
  //Desired number eigenvalues
  unsigned n_eval=6;
  //Set shift
  static_cast<ARPACK*>( problem.eigen_solver_pt() ) -> set_shift( 0.4 );
  
  //Solve the eigenproblem
  cout << "facking" << endl;
  problem.solve_eigenproblem(n_eval,eigenvalues,eigenvectors);
  for(unsigned k=0;k<n_eval;k++)  cout << "eigenvalue: " << eigenvalues[k] << endl;
  cout << "arpack" << endl;
  //Activate pitchfork tracking
  DoubleVector eigenvector(eigenvectors[0]);
  
  //Compute symmetry vector
  DoubleVector symmetry;
  double norm = 1.e-2;
  problem.get_symmetry(symmetry,norm);
  
  ofstream out("scrap_up.dat");
  out << "# Output { Re, Alpha, Int du_dphi(-1), Int du_dphi(1), Delta_Shear_Stress, Jacobian sign }" << endl;
  out << "# Radius ratio: " << R_l << " xmesh: " << xmesh << " ymesh " << ymesh << endl;
  out << "# Alpha = " << Alpha << endl; 
  out << "# New outlet region: " << new_outlet_region << endl; 
  out << "# Inlet traction: " << inlet_traction << " Outlet traction: " << outlet_traction << endl;
  out << "# eta_inlet (Homotopy parameter): " << eta_inlet << " eta_outlet: " << eta_outlet << endl;
  out << "# Pin v at inlet: " << pinv << endl;
  out << "# Using symmetry for pitchfork solve: " << sym << endl << endl;
  
  //Activate pitchfork tracking
  if(sym)
   {
    cout << endl << "Using symmetry vector from get_symmetry routine.  Re = " << Re << endl << endl;
    problem.activate_pitchfork_tracking(&Global_Physical_Variables::Re,symmetry);
   }
  else
   {
    cout << endl << "Using eigenvector for pitchfork detection.  Re = " << Re << endl << endl;
    problem.activate_pitchfork_tracking(&Global_Physical_Variables::Re,eigenvector);
   }
  
  //Find an intitial pitchfork
  problem.newton_solve();
  
  //Check on slack parameter
  unsigned long n_dof = problem.ndof();cout << "ndof: " << n_dof << endl;
  DoubleVector soln;
  problem.get_dofs(soln);
  unsigned long half=(n_dof/2)-1;
  cout << "According to dofs Re is: " << problem.dof(half) 
       << " and slack parameter is: " << problem.dof(n_dof-1) << endl;
  
  double slack=problem.dof(n_dof-1);
  
  problem.deactivate_bifurcation_tracking();
  
  doc_info.number()=10;
  problem.doc_solution(doc_info);
  problem.get_shear_stress(shear_stress);
  cout << "Shear stress is: Lower = " << shear_stress[0] << " Upper = " << shear_stress[1] 
       << " Delta = " << (shear_stress[0]+shear_stress[1]) << endl;
  
  cout << endl << "Pitchfork detected at Re = " << Re << " and Alpha = " << Alpha << endl;
  cout << "========================================================" << endl;
  
  out << Re << " " << Alpha << " " << shear_stress[0] << " " << shear_stress[1] << " " 
      << (shear_stress[0]+shear_stress[1]) << " " << 0 << " " << slack << endl;
  
  
  cout << "========================================================" << endl;
  cout << endl << "Computing eigenvalue at pitchfork.  Re = " << Re << endl << endl;
  //Solve the eigenproblem
  problem.solve_eigenproblem(n_eval,eigenvalues,eigenvectors);
  for(unsigned k=0;k<n_eval;k++)  cout << "eigenvalue: " << eigenvalues[k] << endl;
  
  // OUTPUT THE EIGENMODE
  eigenmode.number()=100;
  problem.store_current_dof_values();
  problem.pin_boundaries_to_zero();
  problem.assign_eigenvector_to_dofs( eigenvectors[0] );
  problem.doc_solution( eigenmode ); 
  problem.restore_dof_values(); 
  
  //Store solution at bifurcation
  n_dof = problem.ndof(); cout << "ndof: " << n_dof << endl;
  DoubleVector soln_at_bifurcation;
  problem.get_dofs(soln_at_bifurcation);
  DoubleVector critical_eigenvector(eigenvectors[0]);
  
  // part one ---------------------------------------------------------------------------------------
  
  double multiplier=-5.0;
  int sign=0;
  
  cout << endl << "Adding eigenvector[0] with multiplier " << multiplier << " to degrees of freedom" << endl; 
  for(unsigned long d=0;d<n_dof;d++)
   {
    problem.dof(d)=soln_at_bifurcation[d]+multiplier*critical_eigenvector[d];
   }
  
  //Now we increace/decreace the Reynolds number slightly and compute a solution
  Re+=epsilon;
  cout << "Adjusting Reynolds number by " << epsilon << " to Re: " << Re
       << " and newton solving" << endl << endl;
  problem.newton_solve();
  problem.get_shear_stress(shear_stress);
  sign=problem.get_Jacobian_sign();
  cout << "Shear stress is: Lower = " << shear_stress[0] << " Upper = " << shear_stress[1] 
       << " Delta = " << (shear_stress[0]+shear_stress[1]) << endl;
  
  out << Re << " " << Alpha << " " << shear_stress[0] << " " << shear_stress[1] << " " << (shear_stress[0]+shear_stress[1]) << " " << sign << endl;
  
  //Now we step in Re and try to continue the bifurcated state
  double ds=Rstep;
  Re-=ds*0.9;
  
  ds=problem.arc_length_step_solve( &Re, ds );
  cout << endl << "Re: " << Re << endl;
  problem.get_shear_stress(shear_stress);
  sign=problem.get_Jacobian_sign();
  cout << "Shear stress is: Lower = " << shear_stress[0] << " Upper = " << shear_stress[1] 
       << " Delta = " << (shear_stress[0]+shear_stress[1]) << endl;
  if(inlet_traction && outlet_traction) cout << " pext: " << problem.get_pext() << endl;
  out << Re << " " << Alpha << " " << shear_stress[0] << " " << shear_stress[1] << " " << (shear_stress[0]+shear_stress[1]) << " " << sign << endl;
  
  for(;Re<Rmax;)
   {
    problem.arc_length_step_solve( &Re, ds );
    cout << endl << "Re: " << Re << endl;
    problem.get_shear_stress(shear_stress);
  sign=problem.get_Jacobian_sign();
    cout << "Shear stress is: Lower = " << shear_stress[0] << " Upper = " << shear_stress[1] 
     << " Delta = " << (shear_stress[0]+shear_stress[1]) << endl;
    if(inlet_traction && outlet_traction) cout << " pext: " << problem.get_pext() << endl;
    out << Re << " " << Alpha << " " << shear_stress[0] << " " << shear_stress[1] << " " << (shear_stress[0]+shear_stress[1]) << " " << sign << endl;
   }
  
  //Compute the picture at Re=Rmax/Re=Rmin
  Re=Rmax;
  cout << endl << "Computing picture at Re: " << Re << endl << endl;
  problem.newton_solve(2);
  doc_info.number()=20;
  problem.doc_solution(doc_info);
  problem.output_streamfunction(doc_info,false);
  doc_info.number()+=1;
  problem.get_shear_stress(shear_stress);
  sign=problem.get_Jacobian_sign();
  out << Re << " " << Alpha << " " << shear_stress[0] << " " << shear_stress[1] << " " << (shear_stress[0]+shear_stress[1]) << " " << sign << endl;
  
  out.close();
  
 } // end_of_main

Classes

Namespaces

Functions

Variables

Function Documentation

◆ main()