two_d_fp.cc File Reference

#include "generic.h"
#include "navier_stokes.h"
#include "meshes/simple_rectangular_tri_mesh.h"
#include "meshes/simple_rectangular_quadmesh.h"
#include "meshes/rectangular_quadmesh.h"
#include "meshes/backward_step_mesh.h"

Classes
class	FpTestProblem
	Test problem for Fp/PCD preconditioner. More...

Namespaces
	Global_Variables
	Namespace for physical parameters.

Enumerations
enum	{ Global_Variables::Driven_cavity , Global_Variables::Step }
	Enumeration for the problem ids. More...

Functions
int	main (int argc, char **argv)
	Driver for Fp preconditioner. More...

Function Documentation

main()

int main	(	int argc	,
		char **	argv
	)

Driver for Fp preconditioner.

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

{
 
 
#ifdef OOMPH_HAS_MPI
 MPI_Helpers::init(argc,argv);
#endif
 
 
 //Label for output
 DocInfo doc_info;
 
 //Set output directory
 doc_info.set_directory("RESLT");
 
 //Doc number of gmres iterations
 ofstream out_file;
 
 // Set flags
 bool use_hypre_for_pressure=true;
 bool use_block_diagonal_for_momentum=false;
 bool use_hypre_for_momentum_diagonals=false;
 
 // Outermost loop over stress divergence or simple form
 for (unsigned do_stress_div=0;do_stress_div<2;do_stress_div++) 
  {
   if (do_stress_div)
    {     
     oomph_info << "Doing stress divergence form\n";
     NavierStokesEquations<2>::Gamma[0]=1.0;
     NavierStokesEquations<2>::Gamma[1]=1.0;
    }
   else
    {     
     oomph_info << "Doing simple form of viscous terms\n";
     NavierStokesEquations<2>::Gamma[0]=0.0;
     NavierStokesEquations<2>::Gamma[1]=0.0;
    }
   
   //Loop over problems: Driven cavity and step
   //unsigned problem_id=Global_Variables::Step;
   for (unsigned problem_id=0;problem_id<2;problem_id++) 
   {
    
    if (problem_id==Global_Variables::Driven_cavity)
     {
      if (do_stress_div==1)
       {
        out_file.open("two_d_iter_driven_cavity_stress_div.dat");
       }
      else
       {
        out_file.open("two_d_iter_driven_cavity_simple.dat");
       }      
     }
    else if (problem_id==Global_Variables::Step)
     {      
      if (do_stress_div==1)
       {
        out_file.open("two_d_iter_step_stress_div.dat");
       }
      else
       {
        out_file.open("two_d_iter_step_simple.dat");
       }
     }
    
    out_file
     << "VARIABLES=\"nel_1d\","
     << "\"ndof\"," 
     << "\"Re\"," 
     << "\" Newton iteration\","
     << "\"GMRES iterations\","
     << "\"Linear solver time\","
     << "\"doc number\""
     << std::endl;
    
    std::string header1;
    std::string header2;
    
     
    // Loop over preconditioners: iprec=0: LSC
    //                            iprec=1: Fp
    //                            iprec=2: Fp without Robin
    bool use_lsc=true;
    bool use_robin=true;
    for (unsigned iprec=0;iprec<2;iprec++)
     {
      
      // Loop over three cases (triangles, non-refineable/refineable quads)
      unsigned icase_lo=0;
      unsigned icase_hi=2;
      if (problem_id==Global_Variables::Step)
       {
        // No triangles for backward step
        icase_lo=1;
        icase_hi=2; 
       }
      for (unsigned icase=icase_lo;icase<=icase_hi;icase++) 
       {
        bool use_triangles=false;
        bool use_adaptivity=false;
        
        switch(icase)
         {
          
         case 0:
          
          header1=" Triangles";
          use_triangles=true;
          use_adaptivity=false;
          break;
          
         case 1:
           
          header1=" Quads";
          use_triangles=false;
          use_adaptivity=false;
          break;
           
         case 2:
           
          header1=" Refineable Quads";
          use_triangles=false;
          use_adaptivity=true;
          break;
           
         default:
          break;
         }
 
        oomph_info << "Doing it with " << header1 << std::endl;
         
        // Set preconditioner
        if (iprec==0)
         {
          use_lsc=true;
          header2=", LSC";
         }
        else if (iprec==1)
         {
          use_lsc=false;
          use_robin=true;
          header2=", Fp";
         }   
        else if (iprec==2)
         {
          use_lsc=false;
          use_robin=false;
          header2=", Fp without Robin";
         }   
         
        oomph_info << "Doing it with " << header2 << " preconditioner\n";
 
        // Write tecplot header
        string header="ZONE T=\""+header1+header2+"\"\n";
        out_file << header;
                 
        // Number of elements in x/y directions (reduced for validaton)
        unsigned max_nel_1d=32; 
        if (argc>1) max_nel_1d=2;
        for (unsigned nel_1d = 2; nel_1d <= max_nel_1d; nel_1d*=2) 
         {
           
          // Build the problem 
          FpTestProblem problem(
           nel_1d,use_triangles,use_adaptivity,use_lsc,use_robin,
           use_hypre_for_pressure,use_block_diagonal_for_momentum,
           use_hypre_for_momentum_diagonals,problem_id);
           
          // Refine a few times
          if (use_adaptivity)
           {
            Vector<unsigned> elements_to_be_refined;
            unsigned nel=problem.mesh_pt()->nelement();
            unsigned e=0;
            while (e<nel)
             {
              elements_to_be_refined.push_back(e);
              e+=7;
             }
            problem.refine_selected_elements(elements_to_be_refined); 
            elements_to_be_refined.clear();
            elements_to_be_refined.push_back(1);
            problem.refine_selected_elements(elements_to_be_refined); 
            elements_to_be_refined.clear();
            elements_to_be_refined.push_back(4);
            problem.refine_selected_elements(elements_to_be_refined); 
 
           }
 
//            DocInfo my_doc_info;           
//            switch(icase)
//             {
             
//             case 0:
             
//              problem.validate_fp<TTaylorHoodElement<2> >();
//              pause("done validation");
 
//              break;
             
//             case 1:
             
//              problem.validate_fp<QTaylorHoodElement<2> >();
//              pause("done validation");
 
//              break;
             
//             case 2:
             
 
//              problem.validate_fp<RefineableQTaylorHoodElement<2> >();
//              pause("done validation");
 
//              break;
             
//             default:
//              break;
//             }
           
 
          // Loop over Reynolds numbers (limited during validation)
          double start_re = 50.0; 
          if (argc>1) start_re=50;
          double end_re = 50.0; 
          for (double re = start_re; re <= end_re; re+=50.0)
           {
            
            // Set Reynolds
            Global_Variables::Re=re;
 
            // Solve the problem 
            problem.newton_solve();
           
            // Doc solution
            problem.doc_solution(doc_info);
            doc_info.number()++;
           
            // Doc iteration counts for each Newton iteration
            unsigned ndof = problem.ndof();
            unsigned iter = Global_Variables::Iterations.size();
 
            // Doc for all Newton iterations or only the last one?
            unsigned j_lo=0;
            //j_lo=iter-1;
            for (unsigned j = j_lo; j < iter; j++)
             {
              out_file
               << nel_1d << " "
               << ndof << " "
               << re << " " 
               << j << " "
               << Global_Variables::Iterations[j] << " "
               << Global_Variables::Linear_solver_time[j] << " "
               << doc_info.number()-1 << " "
               << std::endl;
             }
 
           }
         }     
       }
     }
   
    out_file.close();
   
   }
 
  }
 
 
#ifdef OOMPH_HAS_MPI
   MPI_Helpers::finalize();
#endif
   
   
  } // end_of_main

References Global_Variables::Driven_cavity, e(), Global_Variables::Iterations, j, Global_Variables::Linear_solver_time, oomph::DocInfo::number(), oomph::oomph_info, problem, Global_Variables::Re, oomph::DocInfo::set_directory(), Global_Variables::Step, and oomph::Global_string_for_annotation::string().