two_layer_interface_nonaxisym_perturbations.cc File Reference

#include "generic.h"
#include "linearised_axisym_navier_stokes_elements.h"
#include "linearised_axisym_navier_stokes_elements.cc"
#include "multi_domain_linearised_axisym_navier_stokes_elements.h"
#include "fluid_interface.h"
#include "linearised_axisymmetric_fluid_interface_elements.h"
#include "linearised_axisymmetric_fluid_interface_elements.cc"
#include "perturbed_spines.h"
#include "perturbed_spines.cc"
#include "oomph_crbond_bessel.h"
#include "meshes/two_layer_spine_mesh.h"
#include "two_layer_perturbed_spine_mesh.template.h"
#include "two_layer_perturbed_spine_mesh.template.cc"
#include <fenv.h>

Classes
class	ElementCmp
	Function-type-object to perform comparison of elements. More...
class	BaseStateProblem< BASE_ELEMENT >
	Base state problem class for Tuckerman counter-rotating lids problem. More...
class	PerturbedStateProblem< BASE_ELEMENT, PERTURBED_ELEMENT >

Namespaces
	GlobalPhysicalVariables
	Namespace for physical parameters.

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

Variables
double	GlobalPhysicalVariables::Viscosity_Ratio = 0.1
double	GlobalPhysicalVariables::Density_Ratio = 1.0
double	GlobalPhysicalVariables::Ca = 1.0
	Capillary number. More...
Vector< int >	GlobalPhysicalVariables::Vector_of_azimuthal_mode_numbers
	Vector of azimuthal mode numbers to investigate. More...
Vector< double >	GlobalPhysicalVariables::G (3)
	Direction of gravity. More...

Function Documentation

main()

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

////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// Driver for two-layer rotating cylinder problem, where we directly simulate the instability by kicking off the base and perturbed state problems together and watching to see whether the perturbation grows or decays.

{
 // Store command line arguments
 CommandLineArgs::setup(argc,argv);
 
 // Set duration of timestep
 const double dt = 0.05;
 
 // Set number of elements in radial (r) direction
 const unsigned base_n_r = 10;
 const unsigned perturbed_n_r = 10;
 
 // Set number of elements in axial (z) direction in lower fluid
 const unsigned base_n_z1 = 10;
 const unsigned perturbed_n_z1 = 10;
 
 // Set number of elements in axial (z) direction in upper fluid
 const unsigned base_n_z2 = 10;
 const unsigned perturbed_n_z2 = 10;
 
 // Height of lower fluid layer
 const double h1 = 1.0;
   
 // Height of upper fluid layer
 const double h2 = 1.0;
 
 // Compute the Womersley number
 GlobalPhysicalVariables::ReSt =
  GlobalPhysicalVariables::Re*GlobalPhysicalVariables::St;
 
 // Set direction of gravity (vertically downwards)
 GlobalPhysicalVariables::G[0] = 0.0;
 GlobalPhysicalVariables::G[1] = -1.0;
 GlobalPhysicalVariables::G[2] = 0.0;
 
 // Set azimuthal mode numbers to investigate
 for(int i=0;i<3;i++)
  {
   GlobalPhysicalVariables::Vector_of_azimuthal_mode_numbers.push_back(i);
  }
 
 // Set maximum time to run simulation for
 double t_max = 8.0;
 
 // If we are doing a validation run, set maximum time to be lower
 if(CommandLineArgs::Argc>1) { t_max = 0.1; }
 
 // Build base state problem
 BaseStateProblem<SpineElement<AxisymmetricQCrouzeixRaviartElement >,BDF<2> >
  base_problem(base_n_r,base_n_z1,base_n_z2,h1,h2);
 
 // Determine the size of the azimuthal mode number vector
 const unsigned n_perturbed_problems = 
  GlobalPhysicalVariables::Vector_of_azimuthal_mode_numbers.size();
 
 // Initialise vector of pointers to perturbation problems
 Vector<PerturbedStateProblem<SpineElement
  <AxisymmetricQCrouzeixRaviartElement >,
  PerturbedSpineElement
  <LinearisedAxisymmetricQCrouzeixRaviartMultiDomainElement >,
  BDF<2> >* > perturbed_problem_pt(n_perturbed_problems);
 
 // Build perturbed state problems
 for(unsigned i=0;i<n_perturbed_problems;i++)
  {
   perturbed_problem_pt[i] = new
    PerturbedStateProblem<SpineElement<AxisymmetricQCrouzeixRaviartElement >,
    PerturbedSpineElement
    <LinearisedAxisymmetricQCrouzeixRaviartMultiDomainElement >,
    BDF<2> > (perturbed_n_r,perturbed_n_z1,perturbed_n_z2,h1,h2,
              base_problem.bulk_mesh_pt(),
              GlobalPhysicalVariables::Vector_of_azimuthal_mode_numbers[i]);
  }
 
 // Create interface elements (this has to be done after
 // interaction is set up in perturbed state problem constructor)
 for(unsigned i=0;i<n_perturbed_problems;i++)
  {
   perturbed_problem_pt[i]->create_interface_elements();
  }
 
 // Rebuild global mesh from the sub-meshes
 base_problem.rebuild_global_mesh();
 for(unsigned i=0;i<n_perturbed_problems;i++)
  {
   perturbed_problem_pt[i]->rebuild_global_mesh();
  }
 
 // Assign equation numbers
 base_problem.assign_eqn_numbers();
 for(unsigned i=0;i<n_perturbed_problems;i++)
  {
   perturbed_problem_pt[i]->assign_eqn_numbers();
  }
 
 // Sort the elements in the base state problem (for frontal solver)
 std::sort(base_problem.mesh_pt()->element_pt().begin(),
           base_problem.mesh_pt()->element_pt().end(),
           ElementCmp());
 
 // Sort the elements in the perturbed state problem (for frontal solver)
 for(unsigned i=0;i<n_perturbed_problems;i++)
  {
   std::sort(perturbed_problem_pt[i]->mesh_pt()->element_pt().begin(),
             perturbed_problem_pt[i]->mesh_pt()->element_pt().end(),
             ElementCmp());
  }
 
 // Doc number of elements in both problems
 cout << "Number of elements:" << endl;
 cout << " - base, bulk         = "
      << base_problem.bulk_mesh_pt()->nelement() << endl;
 cout << " - base, surface      = "
      << base_problem.surface_mesh_pt()->nelement() << endl;
 for(unsigned i=0;i<n_perturbed_problems;i++)
  {
   cout << " - perturbed (" << i+1 << " of " 
        << n_perturbed_problems << "), bulk    = "
        << perturbed_problem_pt[i]->bulk_mesh_pt()->nelement() << endl;
   cout << " - perturbed (" << i+1 << " of " 
        << n_perturbed_problems << "), surface = "
        << perturbed_problem_pt[i]->surface_mesh_pt()->nelement() << endl;
  }
 
 // Initialise doc_info object
 DocInfo* doc_info_pt = new DocInfo();
 
 // Set output directory for solutions
 doc_info_pt->set_directory("RESLT");
 
 // Initialise counter for solutions
 doc_info_pt->number()=0;
 
 // Create and initialise trace file
 for(unsigned i=0;i<n_perturbed_problems;i++)
  {
   perturbed_problem_pt[i]->create_trace_file(doc_info_pt);
  }
 
 // Initialise timestep (also sets weights for all timesteppers)
 base_problem.initialise_dt(dt);
 for(unsigned i=0;i<n_perturbed_problems;i++)
  {
   perturbed_problem_pt[i]->initialise_dt(dt);
  }
 
 // Set initial conditions
 base_problem.set_initial_condition();
 for(unsigned i=0;i<n_perturbed_problems;i++)
  {
   perturbed_problem_pt[i]->set_initial_condition();
  }
 
 // Create info file
 ofstream info_file;
 char info_filename[100];
 sprintf(info_filename,"%s/info.dat",doc_info_pt->directory().c_str());
 info_file.open(info_filename);
 info_file.close();
 
 // Output base state solution
 base_problem.doc_solution(doc_info_pt);
 
 // Doc perturbed state initial conditions
 for(unsigned i=0;i<n_perturbed_problems;i++)
  {
   perturbed_problem_pt[i]->doc_solution(doc_info_pt,true);
  }
 
 // Increment counter for solutions
 doc_info_pt->number()++;
 
 // Determine number of timesteps
 const unsigned n_timestep = unsigned(t_max/dt);
 
 // Timestepping loop
 for(unsigned t=0;t<n_timestep;t++)
  {
   // Output timestep and global time to screen
   cout << "  Timestep " << t << " of " << n_timestep
        << "; Time is " << perturbed_problem_pt[0]->time_pt()->time()
        << endl;
   
   // Output timestep and global time to oomph output
   oomph_info << "================================"
              << "==============================\n"
              << "  Timestep " << t << " of " << n_timestep
              << "; Time is " << perturbed_problem_pt[0]->time_pt()->time()
              << "\n=============================="
              << "================================"
              << std::endl;
 
   // Solve perturbed problems
   for(unsigned i=0;i<n_perturbed_problems;i++)
    {
     oomph_info << "\nPerturbed problem (" << i+1 << " of " 
                << n_perturbed_problems 
                << "):\n-----------------------------\n" << endl;
     perturbed_problem_pt[i]->unsteady_newton_solve(dt);
    }
   
   // Doc all problems
   base_problem.doc_solution(doc_info_pt);
   for(unsigned i=0;i<n_perturbed_problems;i++)
    {
     perturbed_problem_pt[i]->doc_solution(doc_info_pt,true);
    }
   
   // Increment counter for solutions
   doc_info_pt->number()++;
   
  } // End of loop over timesteps
 
 // Clear and close trace files
 for(unsigned i=0;i<n_perturbed_problems;i++)
  {
   perturbed_problem_pt[i]->close_trace_file();
  }
 
} // End of main

References oomph::CommandLineArgs::Argc, oomph::Problem::assign_eqn_numbers(), BaseStateProblem< BASE_ELEMENT >::bulk_mesh_pt(), oomph::DocInfo::directory(), BaseStateProblem< BASE_ELEMENT >::doc_solution(), Eigen::placeholders::end, GlobalPhysicalVariables::G, i, oomph::Problem::initialise_dt(), BaseStateProblem< BASE_ELEMENT >::mesh_pt(), oomph::Mesh::nelement(), oomph::DocInfo::number(), oomph::oomph_info, GlobalPhysicalVariables::Re, oomph::Problem::rebuild_global_mesh(), GlobalPhysicalVariables::ReSt, oomph::DocInfo::set_directory(), BaseStateProblem< BASE_ELEMENT >::set_initial_condition(), Flag_definition::setup(), GlobalPhysicalVariables::St, BaseStateProblem< BASE_ELEMENT >::surface_mesh_pt(), plotPSD::t, and GlobalPhysicalVariables::Vector_of_azimuthal_mode_numbers.