three_d_fp.cc File Reference

#include <fenv.h>
#include "generic.h"
#include "navier_stokes.h"
#include "meshes/simple_cubic_mesh.h"
#include "meshes/simple_cubic_tet_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::Through_flow }
	Enumeration for the problem ids. More...

Functions
void	Global_Variables::prescribed_traction (const double &t, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
	Traction at the outflow boundary. More...
int	main (int argc, char **argv)
	Driver for Fp preconditioner. More...

Variables
Vector< double >	Global_Variables::Linear_solver_time
	Storage for linear solver times during Newton steps. 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;
   
 //Loop over problems: Driven cavity and step
 for (unsigned problem_id=0;problem_id<2;problem_id++) 
  {
   
   if (problem_id==0)
    {
     out_file.open("three_d_iter_driven_cavity.dat");
    }
   else 
    {      
     out_file.open("three_d_iter_through_flow.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
   bool use_lsc=true;
   //bool use_robin=true;
   for (unsigned iprec=0;iprec<2;iprec++)
    {
     
     // Loop over three cases (tets, non-refineable/refineable bricks)
     for (unsigned icase=0;icase<=2;icase++) 
      {
       bool use_tets=false;
       bool use_adaptivity=false;
       
       switch(icase)
        {
         
        case 0:
         
         header1=" Tets";
         use_tets=true;
         use_adaptivity=false;
         break;
         
        case 1:
         
         header1=" Bricks";
         use_tets=false;
         use_adaptivity=false;
         break;
         
        case 2:
         
         header1=" Refineable Bricks";
         use_tets=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";
        }   
       
        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=4; 
        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_tets,use_adaptivity,use_lsc,
           use_hypre_for_pressure,use_block_diagonal_for_momentum,
           use_hypre_for_momentum_diagonals,problem_id);
           
          
          
          // Refine a few times
          if (use_adaptivity)
           {             
            // Note: This manually chosen refinement pattern introduces
            // doubly hanging nodes on boundary and in the interior
            Vector<unsigned> elements_to_be_refined;
            unsigned nel=problem.bulk_mesh_pt()->nelement();
            unsigned e=0;
            while (e<nel)
             {
              elements_to_be_refined.push_back(e);
              e+=7;
             }
            problem.refine_selected_elements(0,elements_to_be_refined); 
            elements_to_be_refined.clear();
            elements_to_be_refined.push_back(1);
            problem.refine_selected_elements(0,elements_to_be_refined); 
            elements_to_be_refined.clear();
            elements_to_be_refined.push_back(4);
            problem.refine_selected_elements(0,elements_to_be_refined); 
           }
          // Attach traction elements now with the problem in its
          // most refined state
          if (use_tets)
           {
            problem.create_traction_elements<TTaylorHoodElement<3> >();
           }
          else 
           {
            if (use_adaptivity)
             {
              problem.create_refineable_traction_elements
               <RefineableQTaylorHoodElement<3> >();
             }
            else
             {
              problem.create_traction_elements<QTaylorHoodElement<3> >();
             }
           }
          
          
          
          // Loop over Reynolds numbers (limited during validation)
          double start_re = 0.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 e(), Global_Variables::Iterations, j, Global_Variables::Linear_solver_time, oomph::DocInfo::number(), oomph::oomph_info, problem, Global_Variables::Re, oomph::DocInfo::set_directory(), and oomph::Global_string_for_annotation::string().