#include <fenv.h>
#include "generic.h"
#include "generalised_newtonian_axisym_navier_stokes.h"
#include "solid.h"
#include "constitutive.h"
#include "fluid_interface.h"
#include "meshes/triangle_mesh.h"
#include "overloaded_element_body.h"

Classes
class	oomph::MyTaylorHoodElement< DIM >
	Overload TaylorHood element to modify output. More...

class	oomph::FaceGeometry< MyTaylorHoodElement >

class	oomph::FaceGeometry< FaceGeometry< MyTaylorHoodElement > >

class	AxisymmetricVibratingShellProblem< ELEMENT >
	Problem class to simulate the settling of a viscous drop. More...

Namespaces
	oomph
	DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).

	oomph::Problem_Parameter
	Namespace for Problem Parameter.

Functions
void	oomph::Problem_Parameter::oscillation (const double &time, const Vector< double > &x, Vector< double > &force)
	Function describing the oscillation. More...

double	oomph::Problem_Parameter::free_surface_profile (const double r)

int	main (int argc, char **argv)
	Driver code for vibrating drop problem. More...

Variables
unsigned	oomph::Problem_Parameter::BDF_type =2
	Which bdf timestepper do we use? More...

unsigned	oomph::Problem_Parameter::Nprev_for_extrapolation_of_strain_rate =2

double	oomph::Problem_Parameter::Uniform_element_area =0.004
	Uniform initial element area. More...

std::string	oomph::Problem_Parameter::Directory = "RESLT"
	Output directory. More...

std::string	oomph::Problem_Parameter::Restart_file =""
	Name of restart file. More...

DocInfo	oomph::Problem_Parameter::Doc_info_trace
	Doc info trace file. More...

DocInfo	oomph::Problem_Parameter::Doc_info_soln
	Doc info solutions. More...

double	oomph::Problem_Parameter::External_pressure =0.0
	Set external pressure. More...

double	oomph::Problem_Parameter::Re =1.967e-1
	Reynolds number. More...

double	oomph::Problem_Parameter::ReInvFr =1.0

double	oomph::Problem_Parameter::Ca = 42.19
	Capillary number. More...

double	oomph::Problem_Parameter::St = 1.440e0
	Strouhal number. More...

double	oomph::Problem_Parameter::Re_St = Re*St
	Reynolds multiplied Strouhal. More...

double	oomph::Problem_Parameter::F_max = 40.24
	Maximum force of the applied oscillation (non-dimensionalised) More...

double	oomph::Problem_Parameter::Oscillating_body_force_multiplier =1.0
	Body force multiplier. More...

double	oomph::Problem_Parameter::Initial_settling_time = 0.25
	Time after which oscillation starts. More...

double	oomph::Problem_Parameter::Dt_min =1.0e-6
	Minimum allowable timestep. More...

double	oomph::Problem_Parameter::Nu =0.3
	Pseudo-solid Poisson ratio. More...

double	oomph::Problem_Parameter::C1 =1.0
	Pseudo-solid Mooney-Rivlin parameter. More...

double	oomph::Problem_Parameter::E =2.2
	Pseudo-solid Young's modulus. More...

double	oomph::Problem_Parameter::Shell_radius =1.0
	Radius of the shell. More...

ConstitutiveLaw *	oomph::Problem_Parameter::Constitutive_law_pt =0
	Constitutive law used to determine the mesh deformation. More...

double	oomph::Problem_Parameter::Yield_stress = 2.861e-1
	have a generic yield stress, required for the output function More...

double	oomph::Problem_Parameter::C_yield_stress =2.526e-1
	Yield stress. More...

double	oomph::Problem_Parameter::C_Re =7.827e-1
	Reynolds number. More...

double	oomph::Problem_Parameter::C_St = 2.966e-1
	Strouhal number. More...

double	oomph::Problem_Parameter::C_Re_St = C_Re*C_St
	Reynolds multiplied Strouhal. More...

double	oomph::Problem_Parameter::HB_yield_stress =2.861e-1
	Yield stress. More...

double	oomph::Problem_Parameter::HB_flow_index = 0.828
	Power law exponent. More...

double	oomph::Problem_Parameter::HB_Re =2.430e-1
	Reynolds number. More...

double	oomph::Problem_Parameter::HB_St =1.295e0
	Strouhal number. More...

double	oomph::Problem_Parameter::HB_Re_St = HB_Re*HB_St
	Reynolds multiplied Strouhal. More...

double	oomph::Problem_Parameter::S_alpha =3.181e-1
	Pre-factor in Sisko model K / ( mu_inf T^(n-1) ) More...

double	oomph::Problem_Parameter::S_flow_index =4.89e-6
	Flow index in Sisko model. More...

double	oomph::Problem_Parameter::S_Re =1.967e-1
	Reynolds number. More...

double	oomph::Problem_Parameter::S_St = 1.440e0
	Strouhal number. More...

double	oomph::Problem_Parameter::S_Re_St = S_Re*S_St
	Reynolds multiplied Strouhal. More...

double	oomph::Problem_Parameter::Critical_strain_rate =1.4e-14
	Critical strain rate. More...

GeneralisedNewtonianConstitutiveEquation< 3 > *	oomph::Problem_Parameter::Const_eqn_pt =0
	Fluid constitutive equation. More...

ofstream	oomph::Problem_Parameter::Trace_file
	Trace file. More...

ofstream	oomph::Problem_Parameter::Norm_file
	File to document the norm of the solution (for validation purposes) More...

Function Documentation

◆ main()

int main	(	int argc	,
		char **	argv
	)

Driver code for vibrating drop problem.

 {
  
  ToleranceForVertexMismatchInPolygons::Tolerable_error=1e-3;
  
  // Store command line arguments
  CommandLineArgs::setup(argc,argv);
  
  // Define possible command line arguments and parse the ones that
  // were actually specified
  
  // Suppress writing of restart files (they're big!)
  CommandLineArgs::specify_command_line_flag("--suppress_restart_files");
  
  
  // Use extrapolated strain rate when determining viscosity?
  CommandLineArgs::specify_command_line_flag
   ("--use_current_strainrate_for_viscosity");
  
  // Name of restart file
  CommandLineArgs::specify_command_line_flag
   ("--restart_file",
    &Problem_Parameter::Restart_file);
   
  // Output directory
  CommandLineArgs::specify_command_line_flag(
    "--dir",
    &Problem_Parameter::Directory);
  
  // Multiplier for oscillatory body force (hierher MH hack)
  CommandLineArgs::specify_command_line_flag(
   "--osc_body_force_multiplier",
   &Problem_Parameter::Oscillating_body_force_multiplier);
  
  // BDF TYPE (1,2 or 4)
  CommandLineArgs::specify_command_line_flag(
    "--bdf_type",
    &Problem_Parameter::BDF_type);
  
  // Number of previous values to be used for extrapolation
  // of strain rate (<= Problem_Parameter::BDF_type)
  CommandLineArgs::specify_command_line_flag(
    "--nprev_for_extrapolation_of_strain_rate",
    &Problem_Parameter::Nprev_for_extrapolation_of_strain_rate);
  
  // Min number of steps period (as timestep limitation)
  unsigned min_ntstep_per_period=40;
  CommandLineArgs::specify_command_line_flag(
    "--min_ntstep_per_period",&min_ntstep_per_period);
  
  // Initial number of steps period 
  unsigned initial_ntstep_per_period=40;
  CommandLineArgs::specify_command_line_flag(
    "--initial_ntstep_per_period",&initial_ntstep_per_period);
  
  // Manual adaptation?
  unsigned nadapt_interval=1; 
  CommandLineArgs::specify_command_line_flag(
   "--manual_adapt",&nadapt_interval);
  
  // Uniform element area
  CommandLineArgs::specify_command_line_flag(
   "--uniform_element_area",
   &Problem_Parameter::Uniform_element_area);
  
  // Use Newtonian constitutive equation (and specify its 
  // viscosity ratio)
  double newtonian_viscosity_ratio=1.0;
  CommandLineArgs::specify_command_line_flag(
   "--use_newtonian_constitutive_equation",
   &newtonian_viscosity_ratio);
  
  // Use Casson const eqn
  CommandLineArgs::specify_command_line_flag(
   "--use_casson_constitutive_equation",
   &Problem_Parameter::Critical_strain_rate);
  
  // Use Herschel Bulkley const eqn
  CommandLineArgs::specify_command_line_flag(
   "--use_herschel_bulkley_constitutive_equation",
   &Problem_Parameter::Critical_strain_rate);
  
  // Use Sisko const eqn
  CommandLineArgs::specify_command_line_flag(
   "--use_sisko_constitutive_equation",
   &Problem_Parameter::Critical_strain_rate);
  
  // Temporal accuray
  double epsilon_t=2.0e-7;
  CommandLineArgs::specify_command_line_flag(
   "--epsilon_t",&epsilon_t);
  
  // Flag indicating whether we use fixed point iteration or not 
  // Number of fixed point iterations before we use Aitken extrapolation
  unsigned fixed_point_steps_before_aitken_extrapolation=0;
  CommandLineArgs::specify_command_line_flag(
   "--use_fixed_point_iteration",
   &fixed_point_steps_before_aitken_extrapolation);
  
  // Validation command line flag
  CommandLineArgs::specify_command_line_flag("--validation");
  
  // Parse command line
  CommandLineArgs::parse_and_assign(); 
  
  // Doc what has actually been specified on the command line
  CommandLineArgs::doc_specified_flags();
  
  // Choose constitutive equation
  if (CommandLineArgs::command_line_flag_has_been_set("--use_newtonian_constitutive_equation"))
   {
    oomph_info << "Using Newtonian constitutive equation with viscosity ratio: "
               << newtonian_viscosity_ratio << std::endl;
    Problem_Parameter::Const_eqn_pt = 
     new NewtonianConstitutiveEquation<3>(newtonian_viscosity_ratio);
   }
  else if (CommandLineArgs::command_line_flag_has_been_set
           ("--use_casson_constitutive_equation"))
   {
    oomph_info << "Using Casson constitutive equation"
               << " with cutoff strain rate: "
               << Problem_Parameter::Critical_strain_rate 
               << std::endl;
    Problem_Parameter::Const_eqn_pt = 
     new CassonTanMilRegWithBlendingConstitutiveEquation<3>
     (&Problem_Parameter::C_yield_stress,
      &Problem_Parameter::Critical_strain_rate);
  
    Problem_Parameter::Re=Problem_Parameter::C_Re;
    Problem_Parameter::St=Problem_Parameter::C_St;
    Problem_Parameter::Re_St=Problem_Parameter::C_Re_St;
    Problem_Parameter::Yield_stress=Problem_Parameter::C_yield_stress;
   }
  else if (CommandLineArgs::command_line_flag_has_been_set
           ("--use_herschel_bulkley_constitutive_equation"))
   {
    oomph_info << "Using Herschel Bulkley constitutive equation"
               << " with cutoff strain rate: "
               << Problem_Parameter::Critical_strain_rate 
               << std::endl;
    Problem_Parameter::Const_eqn_pt = 
     new HerschelBulkleyTanMilRegWithBlendingConstitutiveEquation<3>
     (&Problem_Parameter::HB_yield_stress,
      &Problem_Parameter::HB_flow_index,
      &Problem_Parameter::Critical_strain_rate);
  
    Problem_Parameter::Re=Problem_Parameter::HB_Re;
    Problem_Parameter::St=Problem_Parameter::HB_St;
    Problem_Parameter::Re_St=Problem_Parameter::HB_Re_St;
    Problem_Parameter::Yield_stress=Problem_Parameter::HB_yield_stress;
   }
  else if( (CommandLineArgs::command_line_flag_has_been_set
           ("--use_sisko_constitutive_equation")) ||
           (CommandLineArgs::command_line_flag_has_been_set
            ("--validation")) )
   {
    oomph_info << "Using Sisko constitutive equation"
               << " with cutoff strain rate: "
               << Problem_Parameter::Critical_strain_rate 
               << std::endl;
    Problem_Parameter::Const_eqn_pt = 
     new SiskoTanMilRegWithBlendingConstitutiveEquation<3>
     (&Problem_Parameter::S_alpha,
      &Problem_Parameter::S_flow_index,
      &Problem_Parameter::Critical_strain_rate);
  
    Problem_Parameter::Re=Problem_Parameter::S_Re;
    Problem_Parameter::St=Problem_Parameter::S_St;
    Problem_Parameter::Re_St=Problem_Parameter::S_Re_St;
    Problem_Parameter::Yield_stress=0.0;
   }
  else
   {
    oomph_info << "Please select constitutive equation!\n";
    abort();
   }
  
  
  if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
   {
    Problem_Parameter::Initial_settling_time=0.1;
    nadapt_interval=5;
    fixed_point_steps_before_aitken_extrapolation=500;
    Problem_Parameter::Directory="RESLT/";
    Problem_Parameter::Critical_strain_rate=1.0e-12;
   }
  
  
  // Create generalised Hookean constitutive equations
  Problem_Parameter::Constitutive_law_pt = 
   new GeneralisedHookean(&Problem_Parameter::Nu);
  
  // Open trace file
  Problem_Parameter::Trace_file.open((Problem_Parameter::Directory+"/trace.dat").c_str());
  
  // Write trace file
  Problem_Parameter::Trace_file  
   << "VARIABLES="
   << "\"time\", "  // 1
   << "\"drop height\"," // 2 
   << "\"body force\","  // 3
   << "\"volume\"," // 4
   << "\"number of elements\","// 5 
   << "\"max element area \"," // 6 
   << "\"min element area \"," // 7 
   << "\"max spatial error\"," // 8
   << "\"min spatial error\"," // 9
   << "\"norm of strain invariant\","  // 10
   << "\"norm of extrapolated strain invariant\"," // 11
   << "\"norm of error in extrapolated strain invariant\"," // 12
   << "\"min_invariant \"," // 13
   << "\"max_invariant \"," // 14
   << "\"min_viscosity \"," // 15
   << "\"max_viscosity \"," // 16
   << "\"norm of viscosity \"," // 17
   << "\"dt\"," // 18
   << "\"global temporal error norm\"," // 19
   << "\"doc number\"" // 20
   << std::endl; 
  
  
  // Open norm file
  Problem_Parameter::Norm_file.open((Problem_Parameter::Directory+"/norm.dat").c_str()); // hierher needed? Maybe keep 
  // alive for distant future when this becomes a demo driver code...
  
  // Output directory
  Problem_Parameter::Doc_info_trace.set_directory(Problem_Parameter::Directory);
  Problem_Parameter::Doc_info_soln.set_directory(Problem_Parameter::Directory);
  
  
  // Doc constitutive equation
  ofstream constitutive_file;
  constitutive_file.open((Problem_Parameter::Directory+"/constitutive.dat").c_str());
  double invariant_min=1.0e-15;
  double invariant_max=1.0e2;
  unsigned nstep_per_decade=10;
  double ratio=pow(0.1,double(1.0/nstep_per_decade));
  double invariant=invariant_max;
  while (invariant>invariant_min)
   {
    constitutive_file << invariant << " " 
                      << Problem_Parameter::Const_eqn_pt->viscosity(invariant) 
                      << std::endl;
    invariant*=ratio;
   }
  constitutive_file.close();
  
  // Create problem in initial configuration
  AxisymmetricVibratingShellProblem<
  GeneralisedNewtonianProjectableAxisymmetricTaylorHoodElement
   <MyTaylorHoodElement> > problem;
  
  // Timestepping parameters:
  //-------------------------
  
  oomph_info << "Tolerance for adaptive timestepping: epsilon_t = "
             << epsilon_t << std::endl;
  
  // Initial tiemstep
  double dt_new=1.0/double(initial_ntstep_per_period); 
  
  // max. permitted timestep
  double dt_max=1.0/double(min_ntstep_per_period);
  
  // Do timestepping until tmax
  double t_max=500.25; //2.0
  
  // Tolerance for fixed point iteration in %
  double fixed_point_tol=1.0; //1.0e-6;
  
  // Restart?
  if (CommandLineArgs::command_line_flag_has_been_set("--restart_file"))
   {
    problem.restart();
    oomph_info << "Done restart\n";
    
    // Have read in next timestep
    dt_new=problem.next_dt();
    
    // Doc the read-in initial conditions 
    problem.doc_solution();
   }
  else
   {
    oomph_info << "Not doing restart\n";
    
    // Initialise timestepper
    problem.next_dt()=dt_new;
    problem.initialise_dt(problem.next_dt());
    
    // Perform impulsive start from current state
    problem.assign_initial_values_impulsive();
  
    oomph_info << "zeroth order extrapol (one prev value)" << std::endl;
    problem.set_nprev_for_extrapolation_of_strain_rate_for_all_elements(1);
  
    // Take one timestep and doc 
    problem.unsteady_newton_solve(dt_new);
    problem.doc_solution();
    
    // Take a few timesteps with fixed timestep and without adaptation
    // to get some history values under our belt
    unsigned nstep=unsigned(t_max/problem.next_dt()); // hierher //5; 
    oomph_info << "Doing " << nstep 
               << " steps with constant timestep and mesh " << std::endl;
    //for (unsigned i=0;i<nstep;i++)
    unsigned count = 0;
    while (problem.time_pt()->time()<t_max)
     {
      // Adapt?
      if( (CommandLineArgs::command_line_flag_has_been_set("--manual_adapt") ||
           CommandLineArgs::command_line_flag_has_been_set("--validation") ) &&
          ((Problem_Parameter::Doc_info_trace.number())%nadapt_interval==0))
       {
        oomph_info << "hierher: time for another adaptation at doc number ="
                   << Problem_Parameter::Doc_info_trace.number() << std::endl;
        problem.adapt();
       }
  
      if (count==3)
       {
        oomph_info << "next order extrapol (two prev values)" << std::endl;
        problem.set_nprev_for_extrapolation_of_strain_rate_for_all_elements(2);
       }
  
      oomph_info << "HIERHER: SKIPPING FOURTH ORDER STUFF!\n";
      // if (i==7)
      //  {
      //   oomph_info << "next order extrapol (four prev values)" << std::endl;
      //   problem.set_nprev_for_extrapolation_of_strain_rate_for_all_elements(4);
      //  }
  
  
      // problem.unsteady_newton_solve(dt_new);     
  
      // Adaptive timestep with specified temporal tolerance epsilon_t
      dt_new=problem.adaptive_unsteady_newton_solve(problem.next_dt(),epsilon_t);
  
      // Next timestep 
      cout << "Suggested new dt: " << dt_new << std::endl;
      if(dt_new > dt_max)
       {
        dt_new=dt_max;
        cout<<"dt limited (from above) to: "<<dt_new<<endl;
       }
      if(dt_new < Problem_Parameter::Dt_min)
       {
        dt_new=Problem_Parameter::Dt_min;
        cout<<"dt limited (from below) to: "<<Problem_Parameter::Dt_min<<endl;
       }   
      problem.next_dt()=dt_new;
  
      //problem.doc_solution();
  
      //if(i<3) continue;
  
      // Fixed point iteration?
      if ( CommandLineArgs::command_line_flag_has_been_set("--use_fixed_point_iteration")||
           CommandLineArgs::command_line_flag_has_been_set("--validation") )
       {
        oomph_info << "DOING FIXED POINT ITERATION\n";
        
        problem.enable_fixed_point_iteration_for_strain_rate_for_all_elements();
  
        unsigned i_fp=0;
        double error_fixed_point_iteration = DBL_MAX;
  
        //for (unsigned i_fp=0;i_fp<nfixed_point;i_fp++)
        while(error_fixed_point_iteration > fixed_point_tol)
         {
          if(i_fp == fixed_point_steps_before_aitken_extrapolation)
           {
            oomph_info << "DOING AITKEN EXTRAPOLATION\n";
            problem.enable_aitken_extrapolation_for_all_elements();
           }
          oomph_info << "FIXED POINT RE-SOLVE " << i_fp << " FOR t = " 
                     << problem.time_pt()->time() << "\n";
          problem.newton_solve();
          //problem.doc_solution();
          error_fixed_point_iteration=
           problem.calculate_error_of_fixed_point_iteration();
          
          // if( error_fixed_point_iteration < fixed_point_tol)
          //  {   
          //   oomph_info << "Error in fixed point iteration is less than the "
          //              << "specified tolerance of "<<fixed_point_tol<<"%\n";
            
          //   oomph_info << "\nSTOPPING FIXED POINT ITERATION\n";
          //   problem.doc_solution();
          //   break;
          //  }
          
          // if(i_fp==nfixed_point-1)
          //  {
          //   oomph_info << "Maximum number of fixed point iterations ("
          //              <<nfixed_point<<") reached.\n";
  
          //   oomph_info << "\nSTOPPING FIXED POINT ITERATION\n";
          //   problem.doc_solution();
          //   break;
          //  }
  
          // Update latest guess
          problem.update_latest_fixed_point_iteration_guess_for_strain_rate_for_all_elements();
  
          i_fp++;
         }
  
        problem.doc_solution();
  
        // Switch off
        problem.disable_fixed_point_iteration_for_strain_rate_for_all_elements();
       }
      else
       {
        problem.doc_solution();
       }
      
      count++;     
      
      if( (count==20) &&
          CommandLineArgs::command_line_flag_has_been_set("--validation") )
       {
        break;
       }
      
  
     }
    oomph_info << "hierher bailing out after const timesteps\n";
    exit(0);
   }
  
  
  // Now do temporally adaptive solves, explicitly adaptiving every specified number
  //--------------------------------------------------------------------------------
  // of steps
  //---------
  if (CommandLineArgs::command_line_flag_has_been_set("--manual_adapt"))
   {
  while (problem.time_pt()->time()<t_max)
   {
      // Adaptive timestep with specified temporal tolerance epsilon_t
      dt_new=problem.adaptive_unsteady_newton_solve(problem.next_dt(),epsilon_t);
  
    // Next timestep 
    cout << "Suggested new dt: " << dt_new << std::endl;
    if(dt_new > dt_max)
     {
      dt_new=dt_max;
      cout<<"dt limited (from above) to: "<<dt_new<<endl;
     }
    if(dt_new < Problem_Parameter::Dt_min)
     {
      dt_new=Problem_Parameter::Dt_min;
      cout<<"dt limited from below to: "<<Problem_Parameter::Dt_min<<endl;
     }   
    problem.next_dt()=dt_new;
    
    // Doc it...
    problem.doc_solution();
  
    // Adapt?
    if ((Problem_Parameter::Doc_info_trace.number())%nadapt_interval==0)
     {
      oomph_info << "hierher: time for another adaptation at doc number ="
                 << Problem_Parameter::Doc_info_trace.number() << std::endl;
      problem.adapt();
     }
     }
   }
  // Doubly adaptive run
  //--------------------
  else
   {
    unsigned max_adapt=1;
    bool first=false;
    dt_new=problem.doubly_adaptive_unsteady_newton_solve(problem.next_dt(),epsilon_t,
                                                         max_adapt,first);
  
    // Next timestep 
    cout << "Suggested new dt: " << dt_new << std::endl;
    if(dt_new > dt_max)
     {
      dt_new=dt_max;
      cout<<"dt limited (from above) to: "<<dt_new<<endl;
     }
    if(dt_new < Problem_Parameter::Dt_min)
     {
      dt_new=Problem_Parameter::Dt_min;
      cout<<"dt limited from below to: "<<Problem_Parameter::Dt_min<<endl;
     }   
    problem.next_dt()=dt_new;
    
    // Doc it...
    problem.doc_solution();
   }
  
  
  
  
  
  
  
  
  exit(0);
  
 } //End of main

Classes

Namespaces

Functions

Variables

Function Documentation

◆ main()