adaptive_hopf_with_separate_meshes.cc File Reference

#include <iostream>
#include <fstream>
#include <cstdio>
#include "generic.h"
#include "navier_stokes.h"
#include "linearised_navier_stokes.h"
#include "assert.h"

Classes
class	GeneralEllipse
	My own Ellipse class. More...
class	RectangleWithHoleDomain
	Rectangular domain with circular whole. More...
class	RectangleWithHoleMesh< ELEMENT >
	Domain-based mesh for rectangular mesh with circular hole. More...
class	RefineableRectangleWithHoleMesh< ELEMENT >
class	MyNavierStokesElement
class	FlowAroundCylinderProblem< ELEMENT >
	Flow around a cylinder in rectangular domain. More...

Namespaces
	Global_Parameters
	Namespace for global parameters.

Functions
std::complex< double >	Global_Parameters::Eigenfunction_normalisation (1.0, 0.0)
int	main ()
	Driver. More...

Function Documentation

main()

int main

(

)

Driver.

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

{
 // Length and height of domain
 double length = 14.0;
 double height = 1.0/Global_Parameters::B;
 
 //Create a new ellipse object with equal semi-minor and semi-major axes
 //to be the central cylinder
 GeneralEllipse* cylinder_pt = 
  new GeneralEllipse(4.5,0.0,0.5,0.5);
 
 // Create Problem
 FlowAroundCylinderProblem 
  <MyNavierStokesElement/*RefineableQCrouzeixRaviartElement<2>*/,
 RefineableLinearisedQCrouzeixRaviartMultiDomainElement>
  problem(cylinder_pt,length,height);
 
 //Refine the problem uniformly a couple of times
 problem.refine_uniformly(); problem.refine_uniformly();
 
 // Solve adaptively with up to two rounds of refinement
 problem.newton_solve(2);
 
 //Now increase the Reynolds number to 90 in steps of 5
 //without any further adaptation
 for(unsigned i=0;i<3;i++)
  {
   Global_Parameters::Re += 5.0;
   problem.newton_solve();
  }
 
 //Assign memory for the eigenvalues and eigenvectors
 Vector<DoubleVector> eigenvectors;
 double frequency = 0.0;
 
 //If we are reading in from the disk
 if(Global_Parameters::Read_in_eigenfunction_from_disk)
  {
   //Read in the  eigenvalue from the data file
   std::ifstream input("eigen.dat");
   input >> frequency;
   
   //Read in the eigenvector from the data file
   const unsigned n_dof = problem.ndof();
   eigenvectors.resize(2);
   LinearAlgebraDistribution dist(problem.communicator_pt(),n_dof,false);
   //Rebuild the vector
   eigenvectors[0].build(&dist,0.0);
   eigenvectors[1].build(&dist,0.0);
 
   for(unsigned n=0;n<n_dof;n++)
    {
     input >> eigenvectors[0][n];
     input >> eigenvectors[1][n];
    }
   input.close();
  }
 //Otherwise solve the eigenproblem
 else
  {
   Vector<std::complex<double> > eigenvalues;
   //Now solve the eigenproblem
   problem.solve_eigenproblem(6,eigenvalues,eigenvectors);
   frequency = eigenvalues[0].imag();
  }
 
 
 //Create the perturbation problem and solve the eigenproblem on it
 //FlowAroundCylinderLinPertProblem 
 // <RefineableQCrouzeixRaviartElement<2>,RefineableLinearisedQCrouzeixRaviartMu//ltiDomainElement
 // > pert_problem(problem.mesh_pt());
 
 problem.add_eigenproblem();
 
 //Transfer eigenfunction across to perturbation mesh
 problem.transfer_eigenfunction_as_initial_condition(eigenvectors[0],
                                                     eigenvectors[1]);
 
 problem.set_eigenvalue(0.0,frequency);
 
 problem.newton_solve(2); //Slow!
 
 problem.doc_solution();
 
 
 exit(1);
 //Try to find the Hopf exactly by using the initial
 //guesses for the eigenvalue and eigenvector from
 //the data file
 problem.activate_hopf_tracking(&Global_Parameters::Re,
                                frequency,
                                eigenvectors[0],      
                                eigenvectors[1]);      
 //Solve the problem 
 problem.newton_solve();
 //Report the value of the bifurcation
 std::cout << "Hopf bifurcation found at  " << Global_Parameters::Re << " " <<
  problem.dof(problem.ndof()-1) << std::endl;
 
 //Solve it again with three rounds of adaptation
 problem.newton_solve(2);
 
 std::cout << "Hopf bifurcation found at " << Global_Parameters::Re << " " <<
  problem.dof(problem.ndof()-1) << std::endl;
 problem.mesh_pt()->output("bif_soln.dat");
 
 //Output the value of the critical Reynolds number into a data file
 std::ofstream trace("trace.dat");
 trace << Global_Parameters::Re << " " <<
  problem.dof(problem.ndof()-1) << "\n";
 trace.close();
}

References Global_Parameters::B, Global_Physical_Variables::height(), i, n, problem, Global_Parameters::Re, and Global_Parameters::Read_in_eigenfunction_from_disk.