mpi/distribution/bifurcation_tracking/pitchfork.cc File Reference

#include <iostream>
#include <fstream>
#include <cmath>
#include <typeinfo>
#include <algorithm>
#include <cstdio>
#include "generic.h"

Classes
class	Mesh1D< ELEMENT >
class	SSPorousChannelEquations
class	SSPorousChannelElement
class	UniformTranspiration< ELEMENT >

Namespaces
	Global_Physical_Variables
	Global variables.

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

Function Documentation

main()

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

{
#ifdef OOMPH_HAS_MPI
 MPI_Helpers::init(argc,argv);
#endif
 
 //Load the namespace
 using namespace Global_Physical_Variables;
 
 //Set the number of elements in each section of the mesh.
 unsigned Nx1=50, Nx2=50;
 //Set up the problem
 UniformTranspiration<SSPorousChannelElement> problem(Nx1,Nx2);
 
 //Open output file
 std::ostringstream trace_filename;
 trace_filename << "trace" << problem.communicator_pt()->my_rank()
                << ".dat";
 ofstream trace(trace_filename.str().c_str());
 
 //Track the node in the middle
 //The element Nx1 contains this node
 SSPorousChannelElement *Test_pt = 
  dynamic_cast<SSPorousChannelElement *>(problem.mesh_pt()->element_pt(Nx1));
  //Set the position of the node within the element
 Vector<double> s(1); s[0] = 1.0;
 
#ifdef OOMPH_HAS_MPI
 //Set up a dummy partition
 unsigned n_element = problem.mesh_pt()->nelement();
 Vector<unsigned> element_partition(n_element);
 for(unsigned e=0;e<n_element/2;e++) {element_partition[e]=0;}
 for(unsigned e=n_element/2;e<n_element;e++) {element_partition[e]=1;}
 
 //DocInfo mesh_doc_info;
 //bool report_stats=true;
 //mesh_doc_info.set_directory("RESLT_MESH");
 problem.distribute(element_partition);//,mesh_doc_info,report_stats);
 //problem.check_halo_schemes(mesh_doc_info);
 //problem.distribute();
#endif
 
 //Step up in Reynolds number
 for(unsigned i=0;i<6;i++)
  {
   //Increase the Reynolds number by 0.5 each time
   Re = i*0.5;
   problem.steady_newton_solve();
   //Output to the trace file
   trace << Re << " " 
        << Test_pt->interpolated_f(s,1) << " " 
        << Test_pt->interpolated_f(s,0) << std::endl;
  }
 
 
 //Temp to get the distribution
/* DoubleVector res;
 problem.get_residuals(res);
 
 Vector<unsigned> desired_global_eqns;
 for(unsigned i=0;i<400;i++)
  {
   desired_global_eqns.push_back(i);
  }
 
 DoubleVectorHaloScheme test(res.distribution_pt(),desired_global_eqns);
 
 DoubleVectorWithHaloEntries junk(res);
 
 junk.build_halo_scheme(&test);
 
 //Now let's check the synchronisation
 const unsigned n_row_local = junk.nrow_local();
 for(unsigned n=0;n<n_row_local;n++)
  {
   junk[n] = 100.0;
  }
 junk.global_value(200) = 100.0;
 junk.gather();
 
 for(unsigned n=0;n<problem.ndof();n++)
  {
   oomph_info << n << " " << junk.global_value(n) << "\n";
   }
 
#ifdef OOMPH_HAS_MPI
 MPI_Helpers::finalize();
#endif
exit(1);*/
 
 
 
 //Specify the symmetry the hard way
 //Need to consider the possibility of distribution, so
 //only work with local dofs
 unsigned n_dof_local = problem.dof_distribution_pt()->nrow_local();
 Vector<double> backup(n_dof_local);
  for(unsigned n=0;n<n_dof_local;n++)
   {
    backup[n] = problem.dof(n);
   }
  
  //Now sort out the problem
  unsigned n_node=problem.mesh_pt()->nnode();
  for(unsigned n=0;n<n_node;n++)
   {
    Node* nod_pt = problem.mesh_pt()->node_pt(n);
    double x = nod_pt->x(0);
    if(!nod_pt->is_pinned(1))
     {
      nod_pt->set_value(1,x*x);
     }
    if(!nod_pt->is_pinned(0))
     {
      nod_pt->set_value(0,0.0);
     }
   }
  
  //The symmetry vector needs to have the DOF distribution
 
  //LinearAlgebraDistribution dist(problem.communicator_pt(),n_dof,false);
  DoubleVector symm(problem.dof_distribution_pt());
  for(unsigned n=0;n<n_dof_local;n++)
   {
    symm[n] = problem.dof(n);
    problem.dof(n) = backup[n];
   }
 
 
 //Let's try to find the pitchfork it
 problem.activate_pitchfork_tracking(&Re,symm,false);
 
 problem.steady_newton_solve();
 
 std::cout << "Pitchfork detected at " << Re  << std::endl;
 if(problem.communicator_pt()->my_rank()==0)
  {
   std::cout << "The slack parameter is " 
             << problem.dof(2*n_dof_local+1) << std::endl;
  }
 
 //Output to the trace file
 trace << Re << " " 
       << Test_pt->interpolated_f(s,1) << " " 
       << Test_pt->interpolated_f(s,0) << std::endl;
 
 
 // problem.deactivate_bifurcation_tracking();
 // Re += 0.01;
 // problem.steady_newton_solve();
 
#ifdef OOMPH_HAS_MPI
 MPI_Helpers::finalize();
#endif
 }

References e(), i, SSPorousChannelEquations::interpolated_f(), oomph::Data::is_pinned(), n, problem, GlobalPhysicalVariables::Re, s, oomph::Data::set_value(), symm(), plotDoE::x, and oomph::Node::x().