Public Member Functions | Private Types | Private Attributes | List of all members
NavierStokesProblem< ELEMENT > Class Template Reference

Entry flow problem in quarter tube domain. More...

+ Inheritance diagram for NavierStokesProblem< ELEMENT >:

Public Member Functions

 NavierStokesProblem (DocInfo &doc_info, const string &node_file_name, const string &element_file_name, const string &face_file_name)
 Constructor: Pass DocInfo object and file names. More...
 
 ~NavierStokesProblem ()
 Destructor (empty) More...
 
void actions_after_newton_solve ()
 Doc the solution after solve. More...
 
void actions_before_newton_solve ()
 Update the problem specs before solve. More...
 
void doc_solution ()
 Doc the solution. More...
 
TetgenMesh< ELEMENT > * mesh_pt ()
 
 NavierStokesProblem ()
 Constructor. More...
 
 ~NavierStokesProblem ()
 Destructor. More...
 
void actions_before_newton_solve ()
 Update the problem specs before solve (empty) More...
 
void actions_after_newton_solve ()
 Update the problem specs before solve (empty) More...
 
void actions_after_parameter_increase (double *const &parameter_pt)
 Update the problem specs after an increase in a parameter. More...
 
void create_spacetime_mesh ()
 
void apply_boundary_conditions ()
 
void pin_redundant_temporal_nodes ()
 
void assign_time_slice_id ()
 Assign the time slice number to each element. More...
 
void set_up_spacetime_solver ()
 Set up the solver for this Problem. More...
 
void complete_problem_setup ()
 Complete problem setup: pass pointers to physical variables. More...
 
void run_natural_continuation (double &parameter, const double &parameter_target, const unsigned &max_n_parameter_step)
 
void doc_solution (const bool &doc_spacetime_soln=false)
 Doc the solution. More...
 
- Public Member Functions inherited from oomph::Problem
virtual void debug_hook_fct (const unsigned &i)
 
void set_analytic_dparameter (double *const &parameter_pt)
 
void unset_analytic_dparameter (double *const &parameter_pt)
 
bool is_dparameter_calculated_analytically (double *const &parameter_pt)
 
void set_analytic_hessian_products ()
 
void unset_analytic_hessian_products ()
 
bool are_hessian_products_calculated_analytically ()
 
void set_pinned_values_to_zero ()
 
bool distributed () const
 
virtual void actions_before_adapt ()
 
virtual void actions_after_adapt ()
 Actions that are to be performed after a mesh adaptation. More...
 
OomphCommunicator * communicator_pt ()
 access function to the oomph-lib communicator More...
 
const OomphCommunicator * communicator_pt () const
 access function to the oomph-lib communicator, const version More...
 
 Problem ()
 
 Problem (const Problem &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const Problem &)=delete
 Broken assignment operator. More...
 
virtual ~Problem ()
 Virtual destructor to clean up memory. More...
 
Mesh *& mesh_pt ()
 Return a pointer to the global mesh. More...
 
Mesh *const & mesh_pt () const
 Return a pointer to the global mesh (const version) More...
 
Mesh *& mesh_pt (const unsigned &imesh)
 
Mesh *const & mesh_pt (const unsigned &imesh) const
 Return a pointer to the i-th submesh (const version) More...
 
unsigned nsub_mesh () const
 Return number of submeshes. More...
 
unsigned add_sub_mesh (Mesh *const &mesh_pt)
 
void flush_sub_meshes ()
 
void build_global_mesh ()
 
void rebuild_global_mesh ()
 
LinearSolver *& linear_solver_pt ()
 Return a pointer to the linear solver object. More...
 
LinearSolver *const & linear_solver_pt () const
 Return a pointer to the linear solver object (const version) More...
 
LinearSolver *& mass_matrix_solver_for_explicit_timestepper_pt ()
 
LinearSolver * mass_matrix_solver_for_explicit_timestepper_pt () const
 
EigenSolver *& eigen_solver_pt ()
 Return a pointer to the eigen solver object. More...
 
EigenSolver *const & eigen_solver_pt () const
 Return a pointer to the eigen solver object (const version) More...
 
Time *& time_pt ()
 Return a pointer to the global time object. More...
 
Time * time_pt () const
 Return a pointer to the global time object (const version). More...
 
double & time ()
 Return the current value of continuous time. More...
 
double time () const
 Return the current value of continuous time (const version) More...
 
TimeStepper *& time_stepper_pt ()
 
const TimeStepper * time_stepper_pt () const
 
TimeStepper *& time_stepper_pt (const unsigned &i)
 Return a pointer to the i-th timestepper. More...
 
ExplicitTimeStepper *& explicit_time_stepper_pt ()
 Return a pointer to the explicit timestepper. More...
 
unsigned long set_timestepper_for_all_data (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data=false)
 
virtual void shift_time_values ()
 Shift all values along to prepare for next timestep. More...
 
AssemblyHandler *& assembly_handler_pt ()
 Return a pointer to the assembly handler object. More...
 
AssemblyHandler *const & assembly_handler_pt () const
 Return a pointer to the assembly handler object (const version) More...
 
double & minimum_dt ()
 Access function to min timestep in adaptive timestepping. More...
 
double & maximum_dt ()
 Access function to max timestep in adaptive timestepping. More...
 
unsigned & max_newton_iterations ()
 Access function to max Newton iterations before giving up. More...
 
void problem_is_nonlinear (const bool &prob_lin)
 Access function to Problem_is_nonlinear. More...
 
double & max_residuals ()
 
bool & time_adaptive_newton_crash_on_solve_fail ()
 Access function for Time_adaptive_newton_crash_on_solve_fail. More...
 
double & newton_solver_tolerance ()
 
void add_time_stepper_pt (TimeStepper *const &time_stepper_pt)
 
void set_explicit_time_stepper_pt (ExplicitTimeStepper *const &explicit_time_stepper_pt)
 
void initialise_dt (const double &dt)
 
void initialise_dt (const Vector< double > &dt)
 
Data *& global_data_pt (const unsigned &i)
 Return a pointer to the the i-th global data object. More...
 
void add_global_data (Data *const &global_data_pt)
 
void flush_global_data ()
 
LinearAlgebraDistribution *const & dof_distribution_pt () const
 Return the pointer to the dof distribution (read-only) More...
 
unsigned long ndof () const
 Return the number of dofs. More...
 
unsigned ntime_stepper () const
 Return the number of time steppers. More...
 
unsigned nglobal_data () const
 Return the number of global data values. More...
 
unsigned self_test ()
 Self-test: Check meshes and global data. Return 0 for OK. More...
 
void enable_store_local_dof_pt_in_elements ()
 
void disable_store_local_dof_pt_in_elements ()
 
unsigned long assign_eqn_numbers (const bool &assign_local_eqn_numbers=true)
 
void describe_dofs (std::ostream &out= *(oomph_info.stream_pt())) const
 
void enable_discontinuous_formulation ()
 
void disable_discontinuous_formulation ()
 
void get_dofs (DoubleVector &dofs) const
 
void get_dofs (const unsigned &t, DoubleVector &dofs) const
 Return vector of the t'th history value of all dofs. More...
 
void set_dofs (const DoubleVector &dofs)
 Set the values of the dofs. More...
 
void set_dofs (const unsigned &t, DoubleVector &dofs)
 Set the history values of the dofs. More...
 
void set_dofs (const unsigned &t, Vector< double * > &dof_pt)
 
void add_to_dofs (const double &lambda, const DoubleVector &increment_dofs)
 Add lambda x incremenet_dofs[l] to the l-th dof. More...
 
double * global_dof_pt (const unsigned &i)
 
double & dof (const unsigned &i)
 i-th dof in the problem More...
 
double dof (const unsigned &i) const
 i-th dof in the problem (const version) More...
 
double *& dof_pt (const unsigned &i)
 Pointer to i-th dof in the problem. More...
 
double * dof_pt (const unsigned &i) const
 Pointer to i-th dof in the problem (const version) More...
 
virtual void get_inverse_mass_matrix_times_residuals (DoubleVector &Mres)
 
virtual void get_dvaluesdt (DoubleVector &f)
 
virtual void get_residuals (DoubleVector &residuals)
 Get the total residuals Vector for the problem. More...
 
virtual void get_jacobian (DoubleVector &residuals, DenseDoubleMatrix &jacobian)
 
virtual void get_jacobian (DoubleVector &residuals, CRDoubleMatrix &jacobian)
 
virtual void get_jacobian (DoubleVector &residuals, CCDoubleMatrix &jacobian)
 
virtual void get_jacobian (DoubleVector &residuals, SumOfMatrices &jacobian)
 
void get_fd_jacobian (DoubleVector &residuals, DenseMatrix< double > &jacobian)
 Get the full Jacobian by finite differencing. More...
 
void get_derivative_wrt_global_parameter (double *const &parameter_pt, DoubleVector &result)
 
void get_hessian_vector_products (DoubleVectorWithHaloEntries const &Y, Vector< DoubleVectorWithHaloEntries > const &C, Vector< DoubleVectorWithHaloEntries > &product)
 
void solve_eigenproblem (const unsigned &n_eval, Vector< std::complex< double >> &eigenvalue, Vector< DoubleVector > &eigenvector, const bool &steady=true)
 Solve the eigenproblem. More...
 
void solve_eigenproblem (const unsigned &n_eval, Vector< std::complex< double >> &eigenvalue, const bool &steady=true)
 
virtual void get_eigenproblem_matrices (CRDoubleMatrix &mass_matrix, CRDoubleMatrix &main_matrix, const double &shift=0.0)
 
void assign_eigenvector_to_dofs (DoubleVector &eigenvector)
 Assign the eigenvector passed to the function to the dofs. More...
 
void add_eigenvector_to_dofs (const double &epsilon, const DoubleVector &eigenvector)
 
void store_current_dof_values ()
 Store the current values of the degrees of freedom. More...
 
void restore_dof_values ()
 Restore the stored values of the degrees of freedom. More...
 
void enable_jacobian_reuse ()
 
void disable_jacobian_reuse ()
 Disable recycling of Jacobian in Newton iteration. More...
 
bool jacobian_reuse_is_enabled ()
 Is recycling of Jacobian in Newton iteration enabled? More...
 
bool & use_predictor_values_as_initial_guess ()
 
void newton_solve ()
 Use Newton method to solve the problem. More...
 
void enable_globally_convergent_newton_method ()
 enable globally convergent Newton method More...
 
void disable_globally_convergent_newton_method ()
 disable globally convergent Newton method More...
 
void newton_solve (unsigned const &max_adapt)
 
void steady_newton_solve (unsigned const &max_adapt=0)
 
void copy (Problem *orig_problem_pt)
 
virtual Problem * make_copy ()
 
virtual void read (std::ifstream &restart_file, bool &unsteady_restart)
 
virtual void read (std::ifstream &restart_file)
 
virtual void dump (std::ofstream &dump_file) const
 
void dump (const std::string &dump_file_name) const
 
void delete_all_external_storage ()
 
virtual void symmetrise_eigenfunction_for_adaptive_pitchfork_tracking ()
 
double * bifurcation_parameter_pt () const
 
void get_bifurcation_eigenfunction (Vector< DoubleVector > &eigenfunction)
 
void activate_fold_tracking (double *const &parameter_pt, const bool &block_solve=true)
 
void activate_bifurcation_tracking (double *const &parameter_pt, const DoubleVector &eigenvector, const bool &block_solve=true)
 
void activate_bifurcation_tracking (double *const &parameter_pt, const DoubleVector &eigenvector, const DoubleVector &normalisation, const bool &block_solve=true)
 
void activate_pitchfork_tracking (double *const &parameter_pt, const DoubleVector &symmetry_vector, const bool &block_solve=true)
 
void activate_hopf_tracking (double *const &parameter_pt, const bool &block_solve=true)
 
void activate_hopf_tracking (double *const &parameter_pt, const double &omega, const DoubleVector &null_real, const DoubleVector &null_imag, const bool &block_solve=true)
 
void deactivate_bifurcation_tracking ()
 
void reset_assembly_handler_to_default ()
 Reset the system to the standard non-augemented state. More...
 
double arc_length_step_solve (double *const &parameter_pt, const double &ds, const unsigned &max_adapt=0)
 
double arc_length_step_solve (Data *const &data_pt, const unsigned &data_index, const double &ds, const unsigned &max_adapt=0)
 
void reset_arc_length_parameters ()
 
int & sign_of_jacobian ()
 
void explicit_timestep (const double &dt, const bool &shift_values=true)
 Take an explicit timestep of size dt. More...
 
void unsteady_newton_solve (const double &dt)
 
void unsteady_newton_solve (const double &dt, const bool &shift_values)
 
void unsteady_newton_solve (const double &dt, const unsigned &max_adapt, const bool &first, const bool &shift=true)
 
double doubly_adaptive_unsteady_newton_solve (const double &dt, const double &epsilon, const unsigned &max_adapt, const bool &first, const bool &shift=true)
 
double doubly_adaptive_unsteady_newton_solve (const double &dt, const double &epsilon, const unsigned &max_adapt, const unsigned &suppress_resolve_after_spatial_adapt_flag, const bool &first, const bool &shift=true)
 
double adaptive_unsteady_newton_solve (const double &dt_desired, const double &epsilon)
 
double adaptive_unsteady_newton_solve (const double &dt_desired, const double &epsilon, const bool &shift_values)
 
void assign_initial_values_impulsive ()
 
void assign_initial_values_impulsive (const double &dt)
 
void calculate_predictions ()
 Calculate predictions. More...
 
void enable_mass_matrix_reuse ()
 
void disable_mass_matrix_reuse ()
 
bool mass_matrix_reuse_is_enabled ()
 Return whether the mass matrix is being reused. More...
 
void refine_uniformly (const Vector< unsigned > &nrefine_for_mesh)
 
void refine_uniformly (const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
 
void refine_uniformly_and_prune (const Vector< unsigned > &nrefine_for_mesh)
 
void refine_uniformly_and_prune (const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
 
void refine_uniformly (DocInfo &doc_info)
 
void refine_uniformly_and_prune (DocInfo &doc_info)
 
void refine_uniformly ()
 
void refine_uniformly (const unsigned &i_mesh, DocInfo &doc_info)
 Do uniform refinement for submesh i_mesh with documentation. More...
 
void refine_uniformly (const unsigned &i_mesh)
 Do uniform refinement for submesh i_mesh without documentation. More...
 
void p_refine_uniformly (const Vector< unsigned > &nrefine_for_mesh)
 
void p_refine_uniformly (const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
 
void p_refine_uniformly_and_prune (const Vector< unsigned > &nrefine_for_mesh)
 
void p_refine_uniformly_and_prune (const Vector< unsigned > &nrefine_for_mesh, DocInfo &doc_info)
 
void p_refine_uniformly (DocInfo &doc_info)
 
void p_refine_uniformly_and_prune (DocInfo &doc_info)
 
void p_refine_uniformly ()
 
void p_refine_uniformly (const unsigned &i_mesh, DocInfo &doc_info)
 Do uniform p-refinement for submesh i_mesh with documentation. More...
 
void p_refine_uniformly (const unsigned &i_mesh)
 Do uniform p-refinement for submesh i_mesh without documentation. More...
 
void refine_selected_elements (const Vector< unsigned > &elements_to_be_refined)
 
void refine_selected_elements (const Vector< RefineableElement * > &elements_to_be_refined_pt)
 
void refine_selected_elements (const unsigned &i_mesh, const Vector< unsigned > &elements_to_be_refined)
 
void refine_selected_elements (const unsigned &i_mesh, const Vector< RefineableElement * > &elements_to_be_refined_pt)
 
void refine_selected_elements (const Vector< Vector< unsigned >> &elements_to_be_refined)
 
void refine_selected_elements (const Vector< Vector< RefineableElement * >> &elements_to_be_refined_pt)
 
void p_refine_selected_elements (const Vector< unsigned > &elements_to_be_refined)
 
void p_refine_selected_elements (const Vector< PRefineableElement * > &elements_to_be_refined_pt)
 
void p_refine_selected_elements (const unsigned &i_mesh, const Vector< unsigned > &elements_to_be_refined)
 
void p_refine_selected_elements (const unsigned &i_mesh, const Vector< PRefineableElement * > &elements_to_be_refined_pt)
 
void p_refine_selected_elements (const Vector< Vector< unsigned >> &elements_to_be_refined)
 
void p_refine_selected_elements (const Vector< Vector< PRefineableElement * >> &elements_to_be_refined_pt)
 
unsigned unrefine_uniformly ()
 
unsigned unrefine_uniformly (const unsigned &i_mesh)
 
void p_unrefine_uniformly (DocInfo &doc_info)
 
void p_unrefine_uniformly (const unsigned &i_mesh, DocInfo &doc_info)
 Do uniform p-unrefinement for submesh i_mesh without documentation. More...
 
void adapt (unsigned &n_refined, unsigned &n_unrefined)
 
void adapt ()
 
void p_adapt (unsigned &n_refined, unsigned &n_unrefined)
 
void p_adapt ()
 
void adapt_based_on_error_estimates (unsigned &n_refined, unsigned &n_unrefined, Vector< Vector< double >> &elemental_error)
 
void adapt_based_on_error_estimates (Vector< Vector< double >> &elemental_error)
 
void get_all_error_estimates (Vector< Vector< double >> &elemental_error)
 
void doc_errors (DocInfo &doc_info)
 Get max and min error for all elements in submeshes. More...
 
void doc_errors ()
 Get max and min error for all elements in submeshes. More...
 
void enable_info_in_newton_solve ()
 
void disable_info_in_newton_solve ()
 Disable the output of information when in the newton solver. More...
 
- Public Member Functions inherited from oomph::ExplicitTimeSteppableObject
 ExplicitTimeSteppableObject ()
 Empty constructor. More...
 
 ExplicitTimeSteppableObject (const ExplicitTimeSteppableObject &)=delete
 Broken copy constructor. More...
 
void operator= (const ExplicitTimeSteppableObject &)=delete
 Broken assignment operator. More...
 
virtual ~ExplicitTimeSteppableObject ()
 Empty destructor. More...
 
virtual void actions_before_explicit_stage ()
 
virtual void actions_after_explicit_stage ()
 

Private Types

enum  { Diagonal_preconditioner =0 , Lower_triangular_preconditioner =1 }
 
enum  {
  Lower_wall_boundary_id =0 , Outflow_boundary_id =1 , Upper_wall_boundary_id =2 , Inflow_boundary_id =3 ,
  Cylinder_surface_boundary_id =4 , Initial_time_boundary_id =5 , Final_time_boundary_id =6
}
 

Private Attributes

DocInfo Doc_info
 Doc info object. More...
 
IterativeLinearSolver * Solver_pt
 Oomph-lib iterative linear solver. More...
 
Preconditioner * Prec_pt
 Preconditioner. More...
 
ExtrudedCubeMeshFromQuadMesh< ELEMENT > * Bulk_mesh_pt
 Pointer to the space-time mesh. More...
 
RefineableQuadMeshWithMovingCylinder< MyRefineableQTaylorHoodElement > * Spatial_mesh_pt
 
bool Periodicity_has_been_enforced
 

Additional Inherited Members

- Public Types inherited from oomph::Problem
typedef void(* SpatialErrorEstimatorFctPt) (Mesh *&mesh_pt, Vector< double > &elemental_error)
 Function pointer for spatial error estimator. More...
 
typedef void(* SpatialErrorEstimatorWithDocFctPt) (Mesh *&mesh_pt, Vector< double > &elemental_error, DocInfo &doc_info)
 Function pointer for spatial error estimator with doc. More...
 
- Public Attributes inherited from oomph::Problem
bool Shut_up_in_newton_solve
 
- Static Public Attributes inherited from oomph::Problem
static bool Suppress_warning_about_actions_before_read_unstructured_meshes
 
- Protected Types inherited from oomph::Problem
enum  Assembly_method {
  Perform_assembly_using_vectors_of_pairs , Perform_assembly_using_two_vectors , Perform_assembly_using_maps , Perform_assembly_using_lists ,
  Perform_assembly_using_two_arrays
}
 Enumerated flags to determine which sparse assembly method is used. More...
 
- Protected Member Functions inherited from oomph::Problem
unsigned setup_element_count_per_dof ()
 
virtual void sparse_assemble_row_or_column_compressed (Vector< int * > &column_or_row_index, Vector< int * > &row_or_column_start, Vector< double * > &value, Vector< unsigned > &nnz, Vector< double * > &residual, bool compressed_row_flag)
 
virtual void actions_before_newton_convergence_check ()
 
virtual void actions_before_newton_step ()
 
virtual void actions_after_newton_step ()
 
virtual void actions_before_implicit_timestep ()
 
virtual void actions_after_implicit_timestep ()
 
virtual void actions_after_implicit_timestep_and_error_estimation ()
 
virtual void actions_before_explicit_timestep ()
 Actions that should be performed before each explicit time step. More...
 
virtual void actions_after_explicit_timestep ()
 Actions that should be performed after each explicit time step. More...
 
virtual void actions_before_read_unstructured_meshes ()
 
virtual void actions_after_read_unstructured_meshes ()
 
virtual void actions_after_change_in_global_parameter (double *const &parameter_pt)
 
virtual void actions_after_change_in_bifurcation_parameter ()
 
double & dof_derivative (const unsigned &i)
 
double & dof_current (const unsigned &i)
 
virtual void set_initial_condition ()
 
virtual double global_temporal_error_norm ()
 
unsigned newton_solve_continuation (double *const &parameter_pt)
 
unsigned newton_solve_continuation (double *const &parameter_pt, DoubleVector &z)
 
void calculate_continuation_derivatives (double *const &parameter_pt)
 
void calculate_continuation_derivatives (const DoubleVector &z)
 
void calculate_continuation_derivatives_fd (double *const &parameter_pt)
 
bool does_pointer_correspond_to_problem_data (double *const &parameter_pt)
 
void set_consistent_pinned_values_for_continuation ()
 
- Protected Attributes inherited from oomph::Problem
Vector< Problem * > Copy_of_problem_pt
 
std::map< double *, bool > Calculate_dparameter_analytic
 
bool Calculate_hessian_products_analytic
 
LinearAlgebraDistribution * Dof_distribution_pt
 
Vector< double * > Dof_pt
 Vector of pointers to dofs. More...
 
DoubleVectorWithHaloEntries Element_count_per_dof
 
double Relaxation_factor
 
double Newton_solver_tolerance
 
unsigned Max_newton_iterations
 Maximum number of Newton iterations. More...
 
unsigned Nnewton_iter_taken
 
Vector< double > Max_res
 Maximum residuals at start and after each newton iteration. More...
 
double Max_residuals
 
bool Time_adaptive_newton_crash_on_solve_fail
 
bool Jacobian_reuse_is_enabled
 Is re-use of Jacobian in Newton iteration enabled? Default: false. More...
 
bool Jacobian_has_been_computed
 
bool Problem_is_nonlinear
 
bool Pause_at_end_of_sparse_assembly
 
bool Doc_time_in_distribute
 
unsigned Sparse_assembly_method
 
unsigned Sparse_assemble_with_arrays_initial_allocation
 
unsigned Sparse_assemble_with_arrays_allocation_increment
 
Vector< Vector< unsigned > > Sparse_assemble_with_arrays_previous_allocation
 
double Numerical_zero_for_sparse_assembly
 
double FD_step_used_in_get_hessian_vector_products
 
bool Mass_matrix_reuse_is_enabled
 
bool Mass_matrix_has_been_computed
 
bool Discontinuous_element_formulation
 
double Minimum_dt
 Minimum desired dt: if dt falls below this value, exit. More...
 
double Maximum_dt
 Maximum desired dt. More...
 
double DTSF_max_increase
 
double DTSF_min_decrease
 
double Minimum_dt_but_still_proceed
 
bool Scale_arc_length
 Boolean to control whether arc-length should be scaled. More...
 
double Desired_proportion_of_arc_length
 Proportion of the arc-length to taken by the parameter. More...
 
double Theta_squared
 
int Sign_of_jacobian
 Storage for the sign of the global Jacobian. More...
 
double Continuation_direction
 
double Parameter_derivative
 Storage for the derivative of the global parameter wrt arc-length. More...
 
double Parameter_current
 Storage for the present value of the global parameter. More...
 
bool Use_continuation_timestepper
 Boolean to control original or new storage of dof stuff. More...
 
unsigned Dof_derivative_offset
 
unsigned Dof_current_offset
 
Vector< double > Dof_derivative
 Storage for the derivative of the problem variables wrt arc-length. More...
 
Vector< double > Dof_current
 Storage for the present values of the variables. More...
 
double Ds_current
 Storage for the current step value. More...
 
unsigned Desired_newton_iterations_ds
 
double Minimum_ds
 Minimum desired value of arc-length. More...
 
bool Bifurcation_detection
 Boolean to control bifurcation detection via determinant of Jacobian. More...
 
bool Bisect_to_find_bifurcation
 Boolean to control wheter bisection is used to located bifurcation. More...
 
bool First_jacobian_sign_change
 Boolean to indicate whether a sign change has occured in the Jacobian. More...
 
bool Arc_length_step_taken
 Boolean to indicate whether an arc-length step has been taken. More...
 
bool Use_finite_differences_for_continuation_derivatives
 
OomphCommunicator * Communicator_pt
 The communicator for this problem. More...
 
bool Always_take_one_newton_step
 
double Timestep_reduction_factor_after_nonconvergence
 
bool Keep_temporal_error_below_tolerance
 
- Static Protected Attributes inherited from oomph::Problem
static ContinuationStorageScheme Continuation_time_stepper
 Storage for the single static continuation timestorage object. More...
 

Detailed Description

template<class ELEMENT>
class NavierStokesProblem< ELEMENT >

Entry flow problem in quarter tube domain.

NavierStokes problem.

Member Enumeration Documentation

◆ anonymous enum

template<class ELEMENT >
anonymous enum
private
Enumerator
Diagonal_preconditioner 
Lower_triangular_preconditioner 
1076  {
1077  Diagonal_preconditioner=0,
1078  Lower_triangular_preconditioner=1
1079  };
NavierStokesProblem::Lower_triangular_preconditioner
@ Lower_triangular_preconditioner
Definition: space_time_oscillating_cylinder.cc:1078
NavierStokesProblem::Diagonal_preconditioner
@ Diagonal_preconditioner
Definition: space_time_oscillating_cylinder.cc:1077

◆ anonymous enum

template<class ELEMENT >
anonymous enum
private
Enumerator
Lower_wall_boundary_id 
Outflow_boundary_id 
Upper_wall_boundary_id 
Inflow_boundary_id 
Cylinder_surface_boundary_id 
Initial_time_boundary_id 
Final_time_boundary_id 
1096  {
1097  Lower_wall_boundary_id=0,
1098  Outflow_boundary_id=1,
1099  Upper_wall_boundary_id=2,
1100  Inflow_boundary_id=3,
1101  Cylinder_surface_boundary_id=4,
1102  Initial_time_boundary_id=5,
1103  Final_time_boundary_id=6
1104  };
NavierStokesProblem::Initial_time_boundary_id
@ Initial_time_boundary_id
Definition: space_time_oscillating_cylinder.cc:1102
NavierStokesProblem::Upper_wall_boundary_id
@ Upper_wall_boundary_id
Definition: space_time_oscillating_cylinder.cc:1099
NavierStokesProblem::Lower_wall_boundary_id
@ Lower_wall_boundary_id
Definition: space_time_oscillating_cylinder.cc:1097
NavierStokesProblem::Inflow_boundary_id
@ Inflow_boundary_id
Definition: space_time_oscillating_cylinder.cc:1100
NavierStokesProblem::Cylinder_surface_boundary_id
@ Cylinder_surface_boundary_id
Definition: space_time_oscillating_cylinder.cc:1101
NavierStokesProblem::Final_time_boundary_id
@ Final_time_boundary_id
Definition: space_time_oscillating_cylinder.cc:1103
NavierStokesProblem::Outflow_boundary_id
@ Outflow_boundary_id
Definition: space_time_oscillating_cylinder.cc:1098

Constructor & Destructor Documentation

◆ NavierStokesProblem() [1/2]

template<class ELEMENT >
NavierStokesProblem< ELEMENT >::NavierStokesProblem ( DocInfo &  doc_info,
const string &  node_file_name,
const string &  element_file_name,
const string &  face_file_name 
)

Constructor: Pass DocInfo object and file names.

Constructor: Pass DocInfo object mesh files.

111  : Doc_info(doc_info)
112 {
113  //Create mesh
114  Problem::mesh_pt() = new TetgenMesh<ELEMENT>(node_file_name,
115  element_file_name,
116  face_file_name);
117 
118  //Doc the boundaries
119  ofstream some_file;
120  char filename[100];
121  sprintf(filename,"boundaries.dat");
122  some_file.open(filename);
123  mesh_pt()->output_boundaries(some_file);
124  some_file.close();
125 
126 
127  // Set the boundary conditions for this problem: All nodal values are
128  // free by default -- just pin the ones that have Dirichlet conditions
129 
130  // Pin transverse velocities on outer walls -- the outer boundary
131  // behaves like a channel (open along the x-axis) with slippery walls
132  // on the transverse boundaries
133  {
134  unsigned ibound=0;
135  {
136  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
137  for (unsigned inod=0;inod<num_nod;inod++)
138  {
139  mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
140  mesh_pt()->boundary_node_pt(ibound,inod)->pin(2);
141  }
142  }
143  }
144 
145 
146  // Pin all velocity components on internal block -- behaves like
147  // a rigid body
148  {
149  unsigned ibound=1;
150  {
151  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
152  for (unsigned inod=0;inod<num_nod;inod++)
153  {
154  mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
155  mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
156  mesh_pt()->boundary_node_pt(ibound,inod)->pin(2);
157  }
158  }
159  }
160 
161 
162  // Loop over the elements to set up element-specific
163  // things that cannot be handled by constructor
164  unsigned n_element = mesh_pt()->nelement();
165  for(unsigned i=0;i<n_element;i++)
166  {
167  // Upcast from GeneralisedElement to the present element
168  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
169 
170  //Set the Reynolds number, etc
171  el_pt->re_pt() = &Global_Physical_Variables::Re;
172  }
173 
174 
175  //Do equation numbering
176  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
177 
178 } // end_of_constructor
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
NavierStokesProblem::mesh_pt
TetgenMesh< ELEMENT > * mesh_pt()
Definition: mesh_from_tetgen_navier_stokes.cc:86
NavierStokesProblem::Doc_info
DocInfo Doc_info
Doc info object.
Definition: mesh_from_tetgen_navier_stokes.cc:96
oomph::Problem::assign_eqn_numbers
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Definition: problem.cc:1989
oomph::TetgenMesh
Definition: tetgen_mesh.template.h:52
Global_Physical_Variables::Re
double Re
Reynolds number.
Definition: fibre.cc:55
MergeRestartFiles.filename
string filename
Definition: MergeRestartFiles.py:39

References oomph::Problem::assign_eqn_numbers(), MergeRestartFiles::filename, i, NavierStokesProblem< ELEMENT >::mesh_pt(), and Global_Physical_Variables::Re.

◆ ~NavierStokesProblem() [1/2]

template<class ELEMENT >
NavierStokesProblem< ELEMENT >::~NavierStokesProblem
inline

Destructor (empty)

Destructor for NavierStokes problem in cubic domain.

67 {}

◆ NavierStokesProblem() [2/2]

template<class ELEMENT >
NavierStokesProblem< ELEMENT >::NavierStokesProblem

Constructor.

Constructor for NavierStokes problem in cubic domain.

1112  :
1113  Solver_pt(0),
1114  Prec_pt(0),
1115  Bulk_mesh_pt(0),
1116  Spatial_mesh_pt(0),
1117  Periodicity_has_been_enforced(false)
1118 {
1119  // Update the mesh-related parameters
1120  GlobalParameters::update_mesh_parameters();
1121 
1122  // Update the physical parameters
1123  GlobalParameters::update_physical_parameters();
1124 
1125  // Generate the space-time mesh (stored as Bulk_mesh_pt)
1126  create_spacetime_mesh();
1127 
1128  // Assign the appropriate boundary conditions
1129  apply_boundary_conditions();
1130 
1131  // Assign the time-slice ID to each element
1132  // NOTE: This must be done after the boundary conditions have been assigned
1133  // otherwise it won't know about the periodicity.)
1134  assign_time_slice_id();
1135 
1136  // Complete the problem setup to make the elements fully functional
1137  complete_problem_setup();
1138 
1139  // Make it pseudo-traction-free
1140  SpaceTimeNavierStokesMixedOrderEquations<2>::Gamma[0]=0.0;
1141  SpaceTimeNavierStokesMixedOrderEquations<2>::Gamma[1]=0.0;
1142 
1143  // Use the block lower triangular preconditioner
1144  GlobalParameters::Preconditioner=Lower_triangular_preconditioner;
1145 
1146  // Call the auxiliary solver setup function
1147  set_up_spacetime_solver();
1148 
1149  // Assign equation numbers
1150  oomph_info << "\nNumber of equations: " << assign_eqn_numbers() << std::endl;
1151 } // End of NavierStokesProblem
NavierStokesProblem::Periodicity_has_been_enforced
bool Periodicity_has_been_enforced
Definition: space_time_oscillating_cylinder.cc:1092
NavierStokesProblem::Bulk_mesh_pt
ExtrudedCubeMeshFromQuadMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the space-time mesh.
Definition: space_time_oscillating_cylinder.cc:1082
NavierStokesProblem::set_up_spacetime_solver
void set_up_spacetime_solver()
Set up the solver for this Problem.
Definition: space_time_oscillating_cylinder.cc:1395
NavierStokesProblem::Spatial_mesh_pt
RefineableQuadMeshWithMovingCylinder< MyRefineableQTaylorHoodElement > * Spatial_mesh_pt
Definition: space_time_oscillating_cylinder.cc:1089
NavierStokesProblem::complete_problem_setup
void complete_problem_setup()
Complete problem setup: pass pointers to physical variables.
Definition: space_time_oscillating_cylinder.cc:1491
NavierStokesProblem::assign_time_slice_id
void assign_time_slice_id()
Assign the time slice number to each element.
Definition: space_time_oscillating_cylinder.cc:1517
NavierStokesProblem::create_spacetime_mesh
void create_spacetime_mesh()
Definition: space_time_oscillating_cylinder.cc:1248
NavierStokesProblem::Prec_pt
Preconditioner * Prec_pt
Preconditioner.
Definition: space_time_oscillating_cylinder.cc:1072
NavierStokesProblem::apply_boundary_conditions
void apply_boundary_conditions()
Definition: space_time_oscillating_cylinder.cc:1585
NavierStokesProblem::Solver_pt
IterativeLinearSolver * Solver_pt
Oomph-lib iterative linear solver.
Definition: space_time_oscillating_cylinder.cc:1069
oomph::SpaceTimeNavierStokesMixedOrderEquations
Definition: mixed_order_petrov_galerkin_space_time_navier_stokes_elements.h:470
GlobalParameters::update_mesh_parameters
void update_mesh_parameters()
Definition: space_time_oscillating_cylinder.cc:701
GlobalParameters::Preconditioner
unsigned Preconditioner
----------------------—Domain Properties------------------------—
Definition: space_time_oscillating_cylinder.cc:725
GlobalParameters::update_physical_parameters
void update_physical_parameters()
Definition: space_time_oscillating_cylinder.cc:655
oomph::oomph_info
OomphInfo oomph_info
Definition: oomph_definitions.cc:319

References NavierStokesProblem< ELEMENT >::apply_boundary_conditions(), oomph::Problem::assign_eqn_numbers(), NavierStokesProblem< ELEMENT >::assign_time_slice_id(), NavierStokesProblem< ELEMENT >::complete_problem_setup(), NavierStokesProblem< ELEMENT >::create_spacetime_mesh(), NavierStokesProblem< ELEMENT >::Lower_triangular_preconditioner, oomph::oomph_info, GlobalParameters::Preconditioner, NavierStokesProblem< ELEMENT >::set_up_spacetime_solver(), GlobalParameters::update_mesh_parameters(), and GlobalParameters::update_physical_parameters().

◆ ~NavierStokesProblem() [2/2]

template<class ELEMENT >
NavierStokesProblem< ELEMENT >::~NavierStokesProblem ( )

Destructor.

Member Function Documentation

◆ actions_after_newton_solve() [1/2]

template<class ELEMENT >
void NavierStokesProblem< ELEMENT >::actions_after_newton_solve ( )
inlinevirtual

Doc the solution after solve.

Reimplemented from oomph::Problem.

71  {
72  // Doc solution after solving
73  doc_solution();
74 
75  // Increment label for output files
76  Doc_info.number()++;
77  }
NavierStokesProblem::doc_solution
void doc_solution()
Doc the solution.
Definition: mesh_from_tetgen_navier_stokes.cc:205
oomph::DocInfo::number
unsigned & number()
Number used (e.g.) for labeling output files.
Definition: oomph_utilities.h:554

References GlobalParameters::Doc_info, and oomph::DocInfo::number().

◆ actions_after_newton_solve() [2/2]

template<class ELEMENT >
void NavierStokesProblem< ELEMENT >::actions_after_newton_solve ( )
inlinevirtual

Update the problem specs before solve (empty)

Reimplemented from oomph::Problem.

1019 {}

◆ actions_after_parameter_increase()

template<class ELEMENT >
void NavierStokesProblem< ELEMENT >::actions_after_parameter_increase ( double *const &  parameter_pt)
virtual

Update the problem specs after an increase in a parameter.

Reimplemented from oomph::Problem.

1210 {
1211  // If we have been passed a pointer to the Reynolds number or period ratio
1212  if ((parameter_pt==&GlobalParameters::Re)||
1213  (parameter_pt==&GlobalParameters::Period_ratio))
1214  {
1215  // Call the physical parameters update function
1216  GlobalParameters::update_physical_parameters();
1217 
1218  // Update the boundary conditions (for the update of the velocities on
1219  // the cylinder boundary)
1220  apply_boundary_conditions();
1221  }
1222  else if (parameter_pt==&GlobalParameters::Amplitude)
1223  {
1224  // Update the nodal positions
1225  Bulk_mesh_pt->node_update();
1226 
1227  // Update the boundary conditions (for the update of the velocities on
1228  // the cylinder boundary)
1229  apply_boundary_conditions();
1230  }
1231  else
1232  {
1233  // Tell the user
1234  oomph_info << "\nCalled actions_after_parameter_increase(...) but I"
1235  << "\ndon't know what to do for this parameter. I'm going to "
1236  << "\nassume I'm not meant to do anything here. I hope you know"
1237  << "\nwhat you're doing...\n" << std::endl;
1238  }
1239 } // End of actions_after_parameter_increase
GlobalParameters::Amplitude
double Amplitude
Amplitude of the cylinder motion.
Definition: extrude_with_macro_element_representation.cc:254
GlobalParameters::Period_ratio
double Period_ratio
The ratio T_e/T_s.
Definition: flow_past_oscillating_cylinder.cc:212
GlobalParameters::Re
double Re
Reynolds number.
Definition: extrude_triangle_generated_mesh.cc:48

References GlobalParameters::Amplitude, oomph::oomph_info, GlobalParameters::Period_ratio, GlobalParameters::Re, and GlobalParameters::update_physical_parameters().

◆ actions_before_newton_solve() [1/2]

template<class ELEMENT >
void NavierStokesProblem< ELEMENT >::actions_before_newton_solve
virtual

Update the problem specs before solve.

Set the inflow boundary conditions.

Reimplemented from oomph::Problem.

186 {
187 
188  // Apply conditions on obstacle
189  unsigned ibound=1;
190  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
191  for (unsigned inod=0;inod<num_nod;inod++)
192  {
193  // Block moves in z-direction
194  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(2,1.0);
195 
196  }
197 
198 } // end_of_actions_before_newton_solve

◆ actions_before_newton_solve() [2/2]

template<class ELEMENT >
void NavierStokesProblem< ELEMENT >::actions_before_newton_solve ( )
inlinevirtual

Update the problem specs before solve (empty)

Reimplemented from oomph::Problem.

1016 {}

◆ apply_boundary_conditions()

template<class ELEMENT >
void NavierStokesProblem< ELEMENT >::apply_boundary_conditions

Assign the appropriate boundary conditions and enforce periodicity in the time direction

Assign the appropriate boundary conditions, i.e. periodicity in the t-direction. In the x and y-direction apply Dirichlet boundary conditions. In summary: Boundary 0 (t=0) – Periodic in time (with boundary 5) Boundary 1 (y=0) – Dirichlet Boundary 2 (x=1) – Dirichlet Boundary 3 (y=1) – Dirichlet Boundary 4 (x=0) – Dirichlet Boundary 5 (t=1) – Periodic in time (with boundary 0)

1586 {
1587  // Get the start time
1588  double start_t=TimingHelpers::timer();
1589 
1590  // Number of SPATIAL dimensions
1591  unsigned n_dim=3;
1592 
1593  // Storage for the spatial coordinates
1594  Vector<double> spatial_coordinates(n_dim-1,0.0);
1595 
1596  // Storage for the solution (also passes the pressure field)
1597  Vector<double> u_exact(n_dim,0.0);
1598 
1599  // Number of nodes on t=0 boundary
1600  unsigned n_boundary0_node=Bulk_mesh_pt->
1601  nboundary_node(Initial_time_boundary_id);
1602 
1603  // Number of nodes on t=1 boundary
1604  unsigned n_boundary1_node=Bulk_mesh_pt->
1605  nboundary_node(Final_time_boundary_id);
1606 
1607  // Get the number of boundaries in the mesh
1608  unsigned n_boundary=Bulk_mesh_pt->nboundary();
1609 
1610  // Loop over the boundaries
1611  for (unsigned b=0; b<n_boundary; b++)
1612  {
1613  // Get the number of nodes on the b-th boundary
1614  unsigned n_node=Bulk_mesh_pt->nboundary_node(b);
1615 
1616  // If we're on either of the time boundaries
1617  if ((b==Initial_time_boundary_id)||(b==Final_time_boundary_id))
1618  {
1619  // If we need to set up periodicity
1620  if (!Periodicity_has_been_enforced)
1621  {
1622  //----------------------------
1623  // Establish nodal periodicity
1624  //----------------------------
1625  // Make sure there are as many nodes on boundary 0 as there are on
1626  // boundary (n_boundary-1)
1627  if (n_boundary0_node!=n_boundary1_node)
1628  {
1629  // Create an output stream
1630  std::ofstream output_file;
1631 
1632  // Open a file
1633  output_file.open("RESLT/nodes_b0.csv");
1634 
1635  // Loop over the nodes on the t=0 boundary
1636  for (unsigned i=0; i<n_boundary0_node; i++)
1637  {
1638  // Output the coordinates of the i-th node on the t=0 boundary
1639  Bulk_mesh_pt->boundary_node_pt(Initial_time_boundary_id,i)->
1640  output(output_file);
1641  }
1642 
1643  // Close the file
1644  output_file.close();
1645 
1646  // Open a file
1647  output_file.open("RESLT/nodes_b1.csv");
1648 
1649  // Loop over the nodes on the t=0 boundary
1650  for (unsigned i=0; i<n_boundary1_node; i++)
1651  {
1652  // Output the coordinates of the i-th node on the t=1 boundary
1653  Bulk_mesh_pt->boundary_node_pt(Final_time_boundary_id,i)->
1654  output(output_file);
1655  }
1656 
1657  // Close the file
1658  output_file.close();
1659 
1660  // Create an output stream
1661  std::ostringstream error_message_stream;
1662 
1663  // Create an error message
1664  error_message_stream << "Different number of nodes on t=0 and t=1 "
1665  << "boundary!\nThere are " << n_boundary0_node
1666  << " nodes on boundary " << Initial_time_boundary_id
1667  << " and " << n_boundary1_node
1668  << " nodes on boundary " << Final_time_boundary_id
1669  << "!" << std::endl;
1670 
1671  // Throw an error
1672  throw OomphLibError(error_message_stream.str(),
1673  OOMPH_CURRENT_FUNCTION,
1674  OOMPH_EXCEPTION_LOCATION);
1675  }
1676 
1677  // Loop over the nodes on the t=0 boundary
1678  for (unsigned i=0; i<n_boundary0_node; i++)
1679  {
1680  // Get the pointer to the associated node
1681  Node* node0_pt=Bulk_mesh_pt->boundary_node_pt(Initial_time_boundary_id,i);
1682 
1683  // Boolean to indicate whether or not the neighbour has been found
1684  bool has_neighbour_node_been_found=false;
1685 
1686  // Loop over the nodes on the t=1 boundary
1687  for (unsigned j=0; i<n_boundary1_node; j++)
1688  {
1689  // Get the pointer to the associated node
1690  Node* node1_pt=Bulk_mesh_pt->boundary_node_pt(Final_time_boundary_id,j);
1691 
1692  // Distance value
1693  double distance=0.0;
1694 
1695  // Loop over the entries of x
1696  for (unsigned k=0; k<n_dim-1; k++)
1697  {
1698  // Update the distance (2 norm)
1699  distance+=pow(((node0_pt->x(k))-(node1_pt->x(k))),2.0);
1700  }
1701 
1702  // Square root it
1703  distance=std::sqrt(distance);
1704 
1705  // Check if it matches to within a reasonable tolerance
1706  if (std::fabs(distance)<Tree::max_neighbour_finding_tolerance())
1707  {
1708  // Make the nodes periodic; the node on the t=1 boundary now points
1709  // to the node on the t=0 boundary.
1710  node1_pt->make_periodic(node0_pt);
1711 
1712  // We've found the neighbouring node
1713  has_neighbour_node_been_found=true;
1714 
1715  // We're done; break out!
1716  break;
1717  }
1718  } // for (unsigned i=0;i<n_boundary0_node;i++)
1719 
1720  // If we get here and we haven't found the neighbouring node, something's
1721  // wrong so throw an error
1722  if (!has_neighbour_node_been_found)
1723  {
1724  // Throw an error
1725  throw OomphLibError("Couldn't find neighbouring node!",
1726  OOMPH_CURRENT_FUNCTION,
1727  OOMPH_EXCEPTION_LOCATION);
1728  }
1729  } // for (unsigned i=0;i<n_boundary0_node;i++)
1730 
1731  // If we've got here then the periodicity has been set up
1732  Periodicity_has_been_enforced=true;
1733  } // if (!Periodicity_has_been_enforced)
1734  }
1735  // If we're not on the time boundaries
1736  else
1737  {
1738  // If we're on the lower/upper wall
1739  if ((b==Lower_wall_boundary_id)||(b==Upper_wall_boundary_id))
1740  {
1741  // Loop over the nodes on the b-th boundary
1742  for (unsigned n=0; n<n_node; n++)
1743  {
1744  // Get a pointer to the n-th node on the b-th boundary
1745  Node* node_pt=Bulk_mesh_pt->boundary_node_pt(b,n);
1746 
1747  // Loop over the velocity components
1748  for (unsigned i=0; i<n_dim-1; i++)
1749  {
1750  // Pin the i-th velocity component
1751  node_pt->pin(i);
1752  }
1753 
1754  // Apply tow-tank boundary conditions (nonzero horizontal velocity)
1755  node_pt->set_value(0,1.0);
1756 
1757  // Apply tow-tank boundary conditions (no vertical velocity)
1758  node_pt->set_value(1,0.0);
1759  } // for (unsigned n=0;n<n_node;n++)
1760  }
1761  // If we're at the outflow boundary
1762  else if (b==Outflow_boundary_id)
1763  {
1764  // Don't actually do anything; it's better than forcing parallel flow
1765  // at the outlet (no numerical oscillations/boundary layer near exit)
1766  }
1767  // If we're at the inflow boundary
1768  else if (b==Inflow_boundary_id)
1769  {
1770  // Loop over the nodes
1771  for (unsigned n=0; n<n_node; n++)
1772  {
1773  // Pointer to the boundary node
1774  Node* boundary_node_pt=Bulk_mesh_pt->boundary_node_pt(b,n);
1775 
1776  // Pin the horizontal velocity
1777  boundary_node_pt->pin(0);
1778 
1779  // Pin the vertical velocity
1780  boundary_node_pt->pin(1);
1781 
1782  // Uniform flow
1783  boundary_node_pt->set_value(0,1.0);
1784 
1785  // Parallel flow means no vertical velocity component
1786  boundary_node_pt->set_value(1,0.0);
1787  }
1788  }
1789  // If we're dealing with the cylinder surface (apply no-slip conditions)
1790  else if (b==Cylinder_surface_boundary_id)
1791  {
1792  // Loop over the nodes
1793  for (unsigned i_nod=0; i_nod<n_node; i_nod++)
1794  {
1795  // Pointer to the boundary node
1796  Node* boundary_node_pt=Bulk_mesh_pt->boundary_node_pt(b,i_nod);
1797 
1798  // Pin the horizontal velocity at every node
1799  boundary_node_pt->pin(0);
1800 
1801  // Pin the vertical velocity at every node
1802  boundary_node_pt->pin(1);
1803 
1804  // Get the time value associated with this node
1805  double time=boundary_node_pt->x(n_dim-1);
1806 
1807  // Storage for the cylinder velocity
1808  Vector<double> u(2,0.0);
1809 
1810  // Get the velocity of the cylinder at this point in time
1811  GlobalParameters::Cylinder_pt->velocity(time,u);
1812 
1813  //--------------------------------------------------------------------
1814  // Apply no-slip condition for NS on a moving wall node noting:
1815  // Uu = (D/T_e) dR/dt => u = St dR/dt = St dR_cyl/dt,
1816  // where the Strouhal number is defined as St = D/(UT_e).
1817  //--------------------------------------------------------------------
1818  // Set the horizontal velocity value
1819  boundary_node_pt->set_value(0,GlobalParameters::St*u[0]);
1820 
1821  // Set the vertical velocity value
1822  boundary_node_pt->set_value(1,GlobalParameters::St*u[1]);
1823  }
1824  } // if ((b==Lower_wall_boundary_id)||(b==Upper_wall_boundary_id))
1825  } // if ((b==Initial_time_boundary_id)||(b==Final_time_boundary_id))
1826  } // for (unsigned b=0;b<n_bound;b++)
1827 
1828  // Record the end time
1829  double end_t=TimingHelpers::timer();
1830 
1831  // Compute the time taken
1832  double boundary_condition_application_time=end_t-start_t;
1833 
1834  // Output the setup time to the screen
1835  oomph_info << "Time taken for application of boundary conditions [sec]: "
1836  << boundary_condition_application_time << std::endl;
1837 } // End of apply_boundary_conditions
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
b
Scalar * b
Definition: benchVecAdd.cpp:17
OscillatingCylinder::velocity
virtual void velocity(const double &t, Vector< double > &u) const
Velocity at any given point on the rigid cylinder at time, t.
Definition: extrude_with_macro_element_representation.cc:151
oomph::Data::pin
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
oomph::Data::set_value
void set_value(const unsigned &i, const double &value_)
Definition: nodes.h:271
oomph::Node
Definition: nodes.h:906
oomph::Node::x
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
oomph::Node::make_periodic
virtual void make_periodic(Node *const &node_pt)
Definition: nodes.cc:2257
oomph::OomphLibError
Definition: oomph_definitions.h:222
oomph::Problem::time
double & time()
Return the current value of continuous time.
Definition: problem.cc:11531
k
char char char int int * k
Definition: level2_impl.h:374
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
GlobalParameters::Cylinder_pt
OscillatingCylinder * Cylinder_pt
---------------------------------—TIME-INTEGRATION PARAMETERS---—
Definition: extrude_with_macro_element_representation.cc:248
GlobalParameters::St
double St
The default Strouhal number (overloaded with input value if given)
Definition: flow_past_oscillating_cylinder.cc:223
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86
timer
double timer
Definition: oomph_metis_from_parmetis_3.1.1/struct.h:210
output
void output(std::ostream &outfile, const unsigned &nplot)
Overload output function.
Definition: overloaded_element_body.h:490
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References b, GlobalParameters::Cylinder_pt, boost::multiprecision::fabs(), i, j, k, oomph::Node::make_periodic(), n, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, output(), oomph::Data::pin(), Eigen::ArrayBase< Derived >::pow(), oomph::Data::set_value(), sqrt(), GlobalParameters::St, OscillatingCylinder::velocity(), and oomph::Node::x().

Referenced by NavierStokesProblem< ELEMENT >::NavierStokesProblem().

◆ assign_time_slice_id()

template<class ELEMENT >
void NavierStokesProblem< ELEMENT >::assign_time_slice_id

Assign the time slice number to each element.

Set up the map which maps a given element e to it's (i,j,k) coordinates in the mesh where i,j and k are integers (such that 0<=i<N_x, 0<=j<N_y, 0<=k<N_t).

1518 {
1519  // Number of dimensions
1520  unsigned n_dim=3;
1521 
1522  // Space for the local coordinates (at the center of the element)
1523  Vector<double> s(n_dim,0.0);
1524 
1525  // The ID of the coordinate we want, i.e. the ID of the time-direction
1526  unsigned time_index=n_dim-1;
1527 
1528  // Storage for the time coordinate
1529  double t=0.0;
1530 
1531  // Get the number of elements in the mesh
1532  unsigned n_element=Bulk_mesh_pt->nelement();
1533 
1534  // Loop over the elements
1535  for (unsigned i=0; i<n_element; i++)
1536  {
1537  // Get the Eulerian coordinates at the center of the element
1538  t=Bulk_mesh_pt->finite_element_pt(i)->interpolated_x(s,time_index);
1539 
1540  // The length of an element in the time direction
1541  double el_length=GlobalParameters::L_t/double(GlobalParameters::N_t);
1542 
1543  // Subtract the length of half an element off (since we're at the centre
1544  // of the element using s=(0,0,0))
1545  t-=el_length/2.0;
1546 
1547  // Divide by the length of an element
1548  t/=el_length;
1549 
1550  // Store it (have to round first otherwise 0.999 would become 0)
1551  unsigned id=unsigned(GlobalParameters::round(t));
1552 
1553  // Upcast the element
1554  ELEMENT* el_pt=dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(i));
1555 
1556  // Assign the time slice ID
1557  el_pt->set_time_slab_id(id);
1558 
1559  // There are 5 unique dofs per element 3 in common with neighbouring
1560  // time slices (i.e. u, v and p) so 3 unique ones per element. We can
1561  // either use an impulsive start (i.e. pin dofs in the first time slice)
1562  // or periodic boundary conditions (the dofs in the final time slice
1563  // are the dofs in the first time slice) so the actual number of dof
1564  // types is always the same
1565  GlobalParameters::N_dof_type=3*GlobalParameters::N_t;
1566 
1567  // Finally, tell it how many dof types there are in the mesh
1568  el_pt->set_ndof_types(GlobalParameters::N_dof_type);
1569  } // for (unsigned i=0;i<n_element;i++)
1570 } // End of assign_time_slice_id
s
RealScalar s
Definition: level1_cplx_impl.h:130
GlobalParameters::L_t
double L_t
The length of the mesh in the time direction.
Definition: space_time_oscillating_cylinder.cc:691
GlobalParameters::N_t
unsigned N_t
The number of elements in the time direction.
Definition: space_time_oscillating_cylinder.cc:694
GlobalParameters::round
double round(const double &d)
---------------------------------------—DOCUMENTATION HELPERS---—
Definition: flow_past_oscillating_cylinder.cc:410
GlobalParameters::N_dof_type
unsigned N_dof_type
Definition: space_time_oscillating_cylinder.cc:729
plotPSD.t
t
Definition: plotPSD.py:36

References i, GlobalParameters::L_t, GlobalParameters::N_dof_type, GlobalParameters::N_t, GlobalParameters::round(), s, and plotPSD::t.

Referenced by NavierStokesProblem< ELEMENT >::NavierStokesProblem().

◆ complete_problem_setup()

template<class ELEMENT >
void NavierStokesProblem< ELEMENT >::complete_problem_setup

Complete problem setup: pass pointers to physical variables.

Complete problem setup: do anything else that's needed to make the elements fully functional (e.g. pass pointers to problem parameters)

1492 {
1493  // Get the number of elements in the bulk mesh
1494  unsigned n_bulk_element=Bulk_mesh_pt->nelement();
1495 
1496  // Loop over the bulk elements
1497  for (unsigned e=0; e<n_bulk_element; e++)
1498  {
1499  // Upcast to a fluid element
1500  ELEMENT* el_pt=dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
1501 
1502  // Set the Reynolds number (set equal to the Womersley number)
1503  el_pt->re_pt()=&GlobalParameters::Re;
1504 
1505  // Set the Strouhal number
1506  el_pt->re_st_pt()=&GlobalParameters::ReSt;
1507  }
1508 } // End of complete_problem_setup
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
GlobalParameters::ReSt
double ReSt
The Womersley number.
Definition: flow_past_oscillating_cylinder.cc:226

References e(), GlobalParameters::Re, and GlobalParameters::ReSt.

Referenced by NavierStokesProblem< ELEMENT >::NavierStokesProblem().

◆ create_spacetime_mesh()

template<class ELEMENT >
void NavierStokesProblem< ELEMENT >::create_spacetime_mesh

Create the space-time mesh with the chosen number of elements in the time direction and the chosen spatial resolution

Helper function to create the space-time mesh (to be assigned to Bulk_mesh_pt) with the chosen number of elements in the time direction and an appropriate spatial resolution (to capture the time-periodic solution properly).

1249 {
1250  // Storage for the start time
1251  double start_t=0.0;
1252 
1253  //--------------------------
1254  // Generate 2D spatial mesh:
1255  //--------------------------
1256  // Record the start time
1257  start_t=TimingHelpers::timer();
1258 
1259  // Use BDF2
1260  add_time_stepper_pt(new Steady<0>);
1261 
1262  // Make a new cylinder
1263  GlobalParameters::Cylinder_pt=
1264  new OscillatingCylinder(&GlobalParameters::Radius,
1265  &GlobalParameters::Amplitude,
1266  time_pt());
1267 
1268  // Make a new mesh and assign its pointer
1269  Spatial_mesh_pt=
1270  new RefineableQuadMeshWithMovingCylinder<MyRefineableQTaylorHoodElement>(
1271  GlobalParameters::Cylinder_pt,
1272  GlobalParameters::Annular_region_radius,
1273  GlobalParameters::Length_of_central_box,
1274  GlobalParameters::X_left,
1275  GlobalParameters::X_right,
1276  GlobalParameters::Height);
1277 
1278  // Loop over the refinements
1279  for (unsigned i=0; i<GlobalParameters::N_uniform_refinement_before_solve; i++)
1280  {
1281  // Refine the mesh
1282  Spatial_mesh_pt->refine_uniformly();
1283  }
1284 
1285  // Quick statistics about the mesh and it's setup
1286  oomph_info << "\nNumber of nodes in spatial mesh: "
1287  << Spatial_mesh_pt->nnode()
1288  << "\nNumber of elements in spatial mesh: "
1289  << Spatial_mesh_pt->nelement()
1290  << "\nTime taken to generate refined spatial mesh [sec]: "
1291  << TimingHelpers::timer()-start_t << std::endl;
1292 
1293  //---------------------------------
1294  // Generate the 3D space-time mesh:
1295  //---------------------------------
1296  // Indicate that we want the mesh extrusion time to be doc-ed
1297  MeshExtrusionHelpers::Mesh_extrusion_helper.enable_doc_mesh_setup_time();
1298 
1299  // Create the extruded mesh
1300  Bulk_mesh_pt=new ExtrudedCubeMeshFromQuadMesh<ELEMENT>
1301  (Spatial_mesh_pt,GlobalParameters::N_t,GlobalParameters::L_t);
1302 
1303  // The created space-time mesh is the only mesh so assign it
1304  Problem::mesh_pt()=Bulk_mesh_pt;
1305 
1306  // Record the start time
1307  start_t=TimingHelpers::timer();
1308 
1309  // Pin the redundant temporal nodes
1310  pin_redundant_temporal_nodes();
1311 
1312  // Reorder the nodes (don't use the regular setup ordering)
1313  NodeReordering::reorder_nodes(Bulk_mesh_pt,false);
1314 
1315  // Document the setup time
1316  oomph_info << "\nTime taken for redundant node pinning/node reordering [sec]: "
1317  << TimingHelpers::timer()-start_t << std::endl;
1318 
1319  // Check the maximum deformation of the mesh due to the cylinder motion
1320  GlobalParameters::doc_maximum_central_box_deformation();
1321 
1322  // Tell the user we've finished
1323  oomph_info << "\nCompleted mesh generation and documentation!" << std::endl;
1324 } // End of create_spacetime_mesh
NavierStokesProblem::pin_redundant_temporal_nodes
void pin_redundant_temporal_nodes()
Definition: space_time_oscillating_cylinder.cc:1335
OscillatingCylinder
Oscillating cylinder class.
Definition: extrude_with_macro_element_representation.cc:46
oomph::ExtrudedCubeMeshFromQuadMesh
Definition: extruded_cube_mesh_from_quad_mesh_with_macro_elements.template.h:99
oomph::MeshExtrusionHelpers::ExtrusionHelper::enable_doc_mesh_setup_time
void enable_doc_mesh_setup_time()
Enable doc-ing of the mesh setup.
Definition: extruded_cube_mesh_from_quad_mesh_with_macro_elements.template.h:73
oomph::Mesh::nnode
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
oomph::Mesh::nelement
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
oomph::Problem::add_time_stepper_pt
void add_time_stepper_pt(TimeStepper *const &time_stepper_pt)
Definition: problem.cc:1545
oomph::Problem::time_pt
Time *& time_pt()
Return a pointer to the global time object.
Definition: problem.h:1504
oomph::TreeBasedRefineableMeshBase::refine_uniformly
void refine_uniformly(DocInfo &doc_info)
Refine mesh uniformly and doc process.
Definition: refineable_mesh.cc:1772
GlobalParameters::Height
double Height
Height of domain.
Definition: flow_past_oscillating_cylinder.cc:316
GlobalParameters::Annular_region_radius
double Annular_region_radius
The radius of the annular region surrounding the cylinder.
Definition: extrude_with_macro_element_representation.cc:271
GlobalParameters::doc_maximum_central_box_deformation
void doc_maximum_central_box_deformation()
Document the maximum deformation inside the central box.
Definition: flow_past_oscillating_cylinder.cc:348
GlobalParameters::Radius
double Radius
Radius of the cylinder.
Definition: extrude_with_macro_element_representation.cc:251
GlobalParameters::X_left
double X_left
X-coordinate of upstream end of domain.
Definition: flow_past_oscillating_cylinder.cc:319
GlobalParameters::N_uniform_refinement_before_solve
unsigned N_uniform_refinement_before_solve
Number of uniform refinements before any solve.
Definition: extrude_with_macro_element_representation.cc:268
GlobalParameters::Length_of_central_box
double Length_of_central_box
Side-length of the square box in the mesh surrounding the cylinder.
Definition: flow_past_oscillating_cylinder.cc:325
GlobalParameters::X_right
double X_right
X-coordinate of downstream end of domain.
Definition: flow_past_oscillating_cylinder.cc:322
NodeReordering::reorder_nodes
void reorder_nodes(Mesh *mesh_pt, const bool &use_old_ordering)
Definition: space_time_oscillating_cylinder.cc:449
oomph::MeshExtrusionHelpers::Mesh_extrusion_helper
class oomph::MeshExtrusionHelpers::ExtrusionHelper Mesh_extrusion_helper

References GlobalParameters::Amplitude, GlobalParameters::Annular_region_radius, GlobalParameters::Cylinder_pt, GlobalParameters::doc_maximum_central_box_deformation(), oomph::MeshExtrusionHelpers::ExtrusionHelper::enable_doc_mesh_setup_time(), GlobalParameters::Height, i, GlobalParameters::L_t, GlobalParameters::Length_of_central_box, oomph::MeshExtrusionHelpers::Mesh_extrusion_helper, GlobalParameters::N_t, GlobalParameters::N_uniform_refinement_before_solve, oomph::oomph_info, GlobalParameters::Radius, NodeReordering::reorder_nodes(), GlobalParameters::X_left, and GlobalParameters::X_right.

Referenced by NavierStokesProblem< ELEMENT >::NavierStokesProblem().

◆ doc_solution() [1/2]

template<class ELEMENT >
void NavierStokesProblem< ELEMENT >::doc_solution

Doc the solution.

206 {
207 
208  ofstream some_file;
209  char filename[100];
210 
211  // Number of plot points
212  unsigned npts;
213  npts=5;
214 
215  // Output solution
216  sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
217  Doc_info.number());
218  some_file.open(filename);
219  mesh_pt()->output(some_file,npts);
220  some_file.close();
221 
222 } // end_of_doc_solution
oomph::DocInfo::directory
std::string directory() const
Output directory.
Definition: oomph_utilities.h:524

References oomph::DocInfo::directory(), GlobalParameters::Doc_info, MergeRestartFiles::filename, and oomph::DocInfo::number().

◆ doc_solution() [2/2]

template<class ELEMENT >
void NavierStokesProblem< ELEMENT >::doc_solution ( const bool &  doc_spacetime_soln = false)

Doc the solution.

Document the solution.

1993 {
1994  // Start the clock
1995  double timer_s=TimingHelpers::timer();
1996 
1997  // Make an ofstream object to output the solution
1998  std::ofstream some_file;
1999 
2000  // Storage for the filename
2001  char filename[100];
2002 
2003  // Number of plot points to use for the big space-time solution
2004  unsigned n_plot_point=3;
2005 
2006  //-----------------
2007  // Output solution:
2008  //-----------------
2009  // Create the filename suffix
2010  sprintf(filename,"%s/soln%i",
2011  GlobalParameters::Doc_info.directory().c_str(),
2012  GlobalParameters::Doc_info.number());
2013 
2014  // Convert it to a string to allow ease of replacing dots
2015  std::string filename_as_string(filename);
2016 
2017  // Replace all dots with the string "pt"
2018  std::replace(filename_as_string.begin(),filename_as_string.end(),'.','p');
2019 
2020  // Now append the filename extension
2021  filename_as_string+=".dat";
2022 
2023  // Open a file with the constructed filename
2024  some_file.open(filename_as_string.c_str());
2025 
2026  // Set the precision of the outputted data
2027  some_file.precision(20);
2028 
2029  // Output the (numerically) approximated solution
2030  Bulk_mesh_pt->output(some_file,n_plot_point);
2031 
2032  // We're done; close the file
2033  some_file.close();
2034 
2035  // Finally, output the time taken
2036  oomph_info << "Total time for documentation [sec]: "
2037  << TimingHelpers::timer()-timer_s << std::endl;
2038 
2039  // Increment counter for solutions
2040  GlobalParameters::Doc_info.number()++;
2041 } // End of doc_solution
GlobalParameters::Doc_info
DocInfo Doc_info
Helper for documenting.
Definition: extrude_triangle_generated_mesh.cc:57
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286

References GlobalParameters::Doc_info, MergeRestartFiles::filename, oomph::DocInfo::number(), oomph::oomph_info, and oomph::Global_string_for_annotation::string().

◆ mesh_pt()

template<class ELEMENT >
TetgenMesh<ELEMENT>* NavierStokesProblem< ELEMENT >::mesh_pt ( )
inline
87  {
88  return dynamic_cast<TetgenMesh<ELEMENT>*>(Problem::mesh_pt());
89  }

Referenced by NavierStokesProblem< ELEMENT >::NavierStokesProblem().

◆ pin_redundant_temporal_nodes()

template<class ELEMENT >
void NavierStokesProblem< ELEMENT >::pin_redundant_temporal_nodes

The mixed order elements use linear interpolation in time so the only nodes which contribute to the unknowns in the system are those that lie on the temporal boundaries of the elements. Thus, all nodes that do not lie on these boundaries need to be pinned (otherwise we'd get zero rows in the system matrix making it singular...).

1336 {
1337  // Number of nodes in each direction
1338  unsigned n_node_1d=Bulk_mesh_pt->finite_element_pt(0)->nnode_1d();
1339 
1340  // Number of nodes in a space-time element
1341  unsigned n_el_node=Bulk_mesh_pt->finite_element_pt(0)->nnode();
1342 
1343  // Sanity check: only works for 3D space-time elements (2D space + 1D time)
1344  if (n_el_node!=std::pow(n_node_1d,3))
1345  {
1346  // Throw an error
1347  throw OomphLibError("Can currently only deal with 3D space-time elements!",
1348  OOMPH_CURRENT_FUNCTION,
1349  OOMPH_EXCEPTION_LOCATION);
1350  }
1351 
1352  // Get the number of elements in the mesh
1353  unsigned n_element=Bulk_mesh_pt->nelement();
1354 
1355  // Loop over the elements
1356  for (unsigned i=0; i<n_element; i++)
1357  {
1358  // Loop over the nodes
1359  for (unsigned j=0; j<n_el_node; j++)
1360  {
1361  // Storage for the local time slice ID (0<=i_temporal<=NNODE_1D-1)
1362  unsigned j_temporal=0;
1363 
1364  // The spatial node number
1365  unsigned j_spatial=j%(n_node_1d*n_node_1d);
1366 
1367  // Which local time slice are we in?
1368  j_temporal=(j-j_spatial)/(n_node_1d*n_node_1d);
1369 
1370  // If we're not on first/final elemental time slice
1371  if ((j_temporal!=0)&&(j_temporal!=n_node_1d-1))
1372  {
1373  // Get a pointer to the j-th node in the i-th element
1374  Node* node_pt=Bulk_mesh_pt->finite_element_pt(i)->node_pt(j);
1375 
1376  // Get the number of unknowns at this node
1377  unsigned n_value=node_pt->nvalue();
1378 
1379  // Loop over the unknowns
1380  for (unsigned k=0; k<n_value; k++)
1381  {
1382  // Pin the k-th unknown at this node
1383  node_pt->pin(k);
1384  }
1385  } // if ((j_temporal!=0)&&(j_temporal!=n_node_1d-1))
1386  } // for (unsigned j=0;j<n_el_node;j++)
1387  } // for (unsigned i=0;i<n_element;i++)
1388 } // End of pin_redundant_temporal_nodes
oomph::Data::nvalue
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483

References i, j, k, oomph::Data::nvalue(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::Data::pin(), and Eigen::bfloat16_impl::pow().

◆ run_natural_continuation()

template<class ELEMENT >
void NavierStokesProblem< ELEMENT >::run_natural_continuation ( double &  parameter,
const double &  parameter_target,
const unsigned &  max_n_parameter_step 
)

Use continuation in a particular parameter. This should make things sufficiently generic so natural continuation can be used interchangeably with ease for any given input parameter. NOTE: It's important that the function: 

       actions_after_parameter_increase(...)

be overloaded so that anything that needs to be updated after a change in the input parameter is indeed updated.

Use continuation in a particular parameter. This should make things sufficiently generic so natural/arc-length continuation can be used interchangeably with ease for any given input parameter. NOTE: It's important that the function: 

       actions_after_parameter_increase(...)

be overloaded so that anything that needs to be updated after a change in the input parameter is indeed updated.

Inputs:

(1) The parameter that we're running continuation in; (2) The target parameter value; (4) The max. number of steps used to reach parameter_target; (5) (Optional) Whether or not to document the solution at each step.

1864 {
1865  // If we can actually do anything
1866  if ((std::fabs(parameter_target-parameter)>=1.0e-14)&&
1867  (max_n_parameter_step!=0))
1868  {
1869  // Store the string used to denote the input parameter
1870  std::string parameter_string=
1871  GlobalParameters::parameter_to_string(&parameter);
1872 
1873  // Copy the initial parameter value
1874  double initial_parameter=parameter;
1875 
1876  // The size of the amplitude value increment
1877  double parameter_increment=((parameter_target-parameter)/
1878  double(max_n_parameter_step));
1879 
1880  // Tell the user what we're doing
1881  oomph_info << ANSIEscapeCode::Red
1882  << "\nStarting natural continuation process for "
1883  << parameter_string << "!" << ANSIEscapeCode::Reset
1884  << "\n\nInitial " << parameter_string << " value: " << parameter
1885  << "\nTarget " << parameter_string << " value: "
1886  << parameter_target << "\nMax. number of parameter steps: "
1887  << max_n_parameter_step << "\n\n"
1888  << ANSIEscapeCode::Red
1889  << "Parameter value cases (" << parameter_string << "): "
1890  << ANSIEscapeCode::Reset << std::endl;
1891 
1892  // Loop over the increments
1893  for (unsigned i=0; i<max_n_parameter_step+1; i++)
1894  {
1895  // Tell the user
1896  oomph_info << " (" << i << ") "
1897  << initial_parameter+i*parameter_increment
1898  << std::endl;
1899  } // for (unsigned i=0;i<max_n_parameter_step+1;i++)
1900 
1901  // Vector to contain the Problem dofs
1902  DoubleVector dofs_backup;
1903 
1904  // The maximum number of times tono halve the parameter increment
1905  unsigned max_n_reattempt=30;
1906 
1907  // The number of times we've halved the parameter increment
1908  unsigned n_reattempt=0;
1909 
1910  // Loop over the increments
1911  while (std::fabs(parameter_target-parameter)>1.0e-14)
1912  {
1913  // Get the dofs
1914  get_dofs(dofs_backup);
1915 
1916  // Increment the parameter value
1917  parameter+=parameter_increment;
1918 
1919  // Tell the user
1920  oomph_info << "\n" << ANSIEscapeCode::Red
1921  << "Solving for " << parameter_string << " value: "
1922  << ANSIEscapeCode::Reset << parameter << std::endl;
1923 
1924  // Update anything that needs to be changed after a change in this parameter
1925  actions_after_parameter_increase(&parameter);
1926 
1927  // Try doing a solve
1928  try
1929  {
1930  // Solve for this parameter value
1931  newton_solve();
1932  }
1933  // If the simulation threw up an error
1934  catch (NewtonSolverError& error)
1935  {
1936  // Make sure we haven't had to try too many times
1937  if (n_reattempt<max_n_reattempt)
1938  {
1939  // Decrement the parameter value
1940  parameter-=parameter_increment;
1941 
1942  // Make all the necessary updates after a change in this parameter
1943  actions_after_parameter_increase(&parameter);
1944 
1945  // Reset the dofs
1946  set_dofs(dofs_backup);
1947 
1948  // Half the parameter increment
1949  parameter_increment*=0.5;
1950 
1951  // Tell the user
1952  oomph_info << "\n" << ANSIEscapeCode::Red
1953  << "Solve failed! Re-attempt " << n_reattempt+1
1954  << " -- halving parameter increment (for "
1955  << parameter_string << ") to: "
1956  << ANSIEscapeCode::Reset;
1957 
1958  // The parameter increment is actually just this increment
1959  oomph_info << parameter_increment << std::endl;
1960 
1961  // Indicate that we're trying again
1962  n_reattempt++;
1963 
1964  // Try again...
1965  continue;
1966  }
1967  // If we've had to try too many times
1968  else
1969  {
1970  // Create an output stream
1971  std::ostringstream error_message_stream;
1972 
1973  // Create the error message
1974  error_message_stream << "Re-attempted a solve too many times. "
1975  << "Exiting here!" << std::endl;
1976 
1977  // Throw an error
1978  throw OomphLibError(error_message_stream.str(),
1979  OOMPH_CURRENT_FUNCTION,
1980  OOMPH_EXCEPTION_LOCATION);
1981  } // if (n_reattempt<max_n_reattempt)
1982  } // try
1983  } // while (std::fabs(parameter_target-parameter)<1.0e-08)
1984  } // if ((std::fabs(parameter_target-parameter)>1.0e-10)&&(max...
1985 } // End of run_natural_continuation
NavierStokesProblem::actions_after_parameter_increase
void actions_after_parameter_increase(double *const &parameter_pt)
Update the problem specs after an increase in a parameter.
Definition: space_time_oscillating_cylinder.cc:1209
oomph::DoubleVector
Definition: double_vector.h:58
oomph::NewtonSolverError
A class to handle errors in the Newton solver.
Definition: problem.h:2952
oomph::Problem::set_dofs
void set_dofs(const DoubleVector &dofs)
Set the values of the dofs.
Definition: problem.cc:3411
oomph::Problem::get_dofs
void get_dofs(DoubleVector &dofs) const
Definition: problem.cc:2479
oomph::Problem::newton_solve
void newton_solve()
Use Newton method to solve the problem.
Definition: problem.cc:8783
GlobalParameters::parameter_to_string
std::string parameter_to_string(const double *const parameter_pt)
Definition: space_time_oscillating_cylinder.cc:866
calibrate.error
int error
Definition: calibrate.py:297
oomph::ANSIEscapeCode::Red
std::string Red
Definition: oomph_utilities.cc:71
oomph::ANSIEscapeCode::Reset
std::string Reset
Definition: oomph_utilities.cc:77

References calibrate::error, boost::multiprecision::fabs(), i, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, GlobalParameters::parameter_to_string(), oomph::ANSIEscapeCode::Red, oomph::ANSIEscapeCode::Reset, and oomph::Global_string_for_annotation::string().

◆ set_up_spacetime_solver()

template<class ELEMENT >
void NavierStokesProblem< ELEMENT >::set_up_spacetime_solver

Set up the solver for this Problem.

Assign the chosen solver to this Problem (and preconditioner if so desired)

1396 {
1397  // Create oomph-lib iterative linear solver
1398  Solver_pt=new GMRES<CRDoubleMatrix>;
1399 
1400  // Use RHS preconditioning
1401  //dynamic_cast<GMRES<CRDoubleMatrix>*>(Solver_pt)->set_preconditioner_RHS();
1402  dynamic_cast<GMRES<CRDoubleMatrix>*>(Solver_pt)->set_preconditioner_LHS();
1403 
1404  // Set the tolerance
1405  Solver_pt->tolerance()=1.0e-10;
1406 
1407  // Maximum number of iterations
1408  Solver_pt->max_iter()=200;
1409 
1410  // Set linear solver
1411  linear_solver_pt()=Solver_pt;
1412 
1413  //-----------------------------------------
1414  // Create the master-level preconditioners:
1415  //-----------------------------------------
1416  // Do we want to document the memory usage?
1417  bool document_memory_usage=true;
1418 
1419  // Solve the diagonal blocks associated with each time-slice separately
1420  if (GlobalParameters::Preconditioner==Diagonal_preconditioner)
1421  {
1422  // Create a new instance of the space-time preconditioner
1423  Prec_pt=new BlockDiagonalPreconditioner<CRDoubleMatrix>;
1424  }
1425  // Solve the block lower-triangular part of the system matrix
1426  else if (GlobalParameters::Preconditioner==Lower_triangular_preconditioner)
1427  {
1428  // Create a new instance of the space-time preconditioner
1429  Prec_pt=new BandedBlockTriangularPreconditioner<CRDoubleMatrix>;
1430 
1431  // Indicate that we're using a block lower triangular solve
1432  dynamic_cast<BandedBlockTriangularPreconditioner<CRDoubleMatrix>*>
1433  (Prec_pt)->lower_triangular();
1434 
1435  // The order of the interpolation in the time direction (linear)
1436  unsigned temporal_order=1;
1437 
1438  // Provide the bandwidth; only subdiagonal block entries
1439  dynamic_cast<BandedBlockTriangularPreconditioner<CRDoubleMatrix>*>
1440  (Prec_pt)->set_block_bandwidth(temporal_order);
1441 
1442  // If we want to document the memory usage
1443  if (document_memory_usage)
1444  {
1445  dynamic_cast<BandedBlockTriangularPreconditioner<CRDoubleMatrix>*>
1446  (Prec_pt)->enable_doc_memory_usage();
1447  }
1448  }
1449  // If the user provided an invalid input
1450  else
1451  {
1452  // Throw an error
1453  throw OomphLibError("Invalid choice of preconditioner.",
1454  OOMPH_CURRENT_FUNCTION,
1455  OOMPH_EXCEPTION_LOCATION);
1456  }
1457 
1458  // Allocate space for the DOF to block map; this tells the master block
1459  // preconditioner which DOF types to aggregate
1460  Vector<unsigned> dof_to_block_map;
1461 
1462  // Call the auxiliary function which sets up the mapping
1463  GlobalParameters::set_up_dof_to_block_mapping(dof_to_block_map);
1464 
1465  // Create an upcasted pointer to the master preconditioner
1466  GeneralPurposeBlockPreconditioner<CRDoubleMatrix>* upcasted_master_prec_pt=
1467  dynamic_cast<GeneralPurposeBlockPreconditioner<CRDoubleMatrix>*>(Prec_pt);
1468 
1469  // Build silently!
1470  upcasted_master_prec_pt->enable_silent_preconditioner_setup();
1471 
1472  // Pass the DOF-to-block map to the preconditioner
1473  upcasted_master_prec_pt->set_dof_to_block_map(dof_to_block_map);
1474 
1475  // Specify the subsidiary block preconditioner
1476  upcasted_master_prec_pt->set_subsidiary_preconditioner_function(
1477  SubsidiaryPreconditionerHelper::get_new_preconditioner);
1478 
1479  // Pass a pointer to the (space-time) mesh
1480  upcasted_master_prec_pt->add_mesh(Bulk_mesh_pt);
1481 
1482  // Now assign the preconditioner to the linear solver
1483  Solver_pt->preconditioner_pt()=Prec_pt;
1484 } // End of set_up_spacetime_solver
oomph::BandedBlockTriangularPreconditioner
Definition: general_purpose_space_time_block_preconditioner.h:199
oomph::BlockDiagonalPreconditioner
Definition: general_purpose_block_preconditioners.h:321
oomph::GMRES
The GMRES method.
Definition: iterative_linear_solver.h:1227
oomph::GeneralPurposeBlockPreconditioner
Definition: general_purpose_block_preconditioners.h:77
oomph::GeneralPurposeBlockPreconditioner::set_dof_to_block_map
void set_dof_to_block_map(Vector< unsigned > &dof_to_block_map)
Specify a DOF to block map.
Definition: general_purpose_block_preconditioners.h:182
oomph::GeneralPurposeBlockPreconditioner::add_mesh
void add_mesh(const Mesh *mesh_pt, const bool &allow_multiple_element_type_in_mesh=false)
Definition: general_purpose_block_preconditioners.h:191
oomph::GeneralPurposeBlockPreconditioner::set_subsidiary_preconditioner_function
void set_subsidiary_preconditioner_function(SubsidiaryPreconditionerFctPt sub_prec_fn)
access function to set the subsidiary preconditioner function.
Definition: general_purpose_block_preconditioners.h:136
oomph::IterativeLinearSolver::preconditioner_pt
Preconditioner *& preconditioner_pt()
Access function to preconditioner.
Definition: iterative_linear_solver.h:95
oomph::IterativeLinearSolver::tolerance
double & tolerance()
Access to convergence tolerance.
Definition: iterative_linear_solver.h:107
oomph::IterativeLinearSolver::max_iter
unsigned & max_iter()
Access to max. number of iterations.
Definition: iterative_linear_solver.h:113
oomph::Preconditioner::enable_silent_preconditioner_setup
void enable_silent_preconditioner_setup()
Set up the block preconditioner quietly!
Definition: preconditioner.h:127
oomph::Problem::linear_solver_pt
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
Definition: problem.h:1466
GlobalParameters::set_up_dof_to_block_mapping
void set_up_dof_to_block_mapping(Vector< unsigned > &dof_to_block_map)
Definition: space_time_oscillating_cylinder.cc:735
SubsidiaryPreconditionerHelper::get_new_preconditioner
Preconditioner * get_new_preconditioner()
Definition: space_time_oscillating_cylinder.cc:993

References oomph::GeneralPurposeBlockPreconditioner< MATRIX >::add_mesh(), oomph::Preconditioner::enable_silent_preconditioner_setup(), SubsidiaryPreconditionerHelper::get_new_preconditioner(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, GlobalParameters::Preconditioner, oomph::GeneralPurposeBlockPreconditioner< MATRIX >::set_dof_to_block_map(), oomph::GeneralPurposeBlockPreconditioner< MATRIX >::set_subsidiary_preconditioner_function(), and GlobalParameters::set_up_dof_to_block_mapping().

Referenced by NavierStokesProblem< ELEMENT >::NavierStokesProblem().

Member Data Documentation

◆ Bulk_mesh_pt

template<class ELEMENT >
ExtrudedCubeMeshFromQuadMesh<ELEMENT>* NavierStokesProblem< ELEMENT >::Bulk_mesh_pt
private

Pointer to the space-time mesh.

◆ Doc_info

template<class ELEMENT >
DocInfo NavierStokesProblem< ELEMENT >::Doc_info
private

Doc info object.

◆ Periodicity_has_been_enforced

template<class ELEMENT >
bool NavierStokesProblem< ELEMENT >::Periodicity_has_been_enforced
private

◆ Prec_pt

template<class ELEMENT >
Preconditioner* NavierStokesProblem< ELEMENT >::Prec_pt
private

Preconditioner.

◆ Solver_pt

template<class ELEMENT >
IterativeLinearSolver* NavierStokesProblem< ELEMENT >::Solver_pt
private

Oomph-lib iterative linear solver.

◆ Spatial_mesh_pt

template<class ELEMENT >
RefineableQuadMeshWithMovingCylinder< MyRefineableQTaylorHoodElement>* NavierStokesProblem< ELEMENT >::Spatial_mesh_pt
private

Pointer to the spatial mesh. NOTE: We keep this so that we retain a pointer to the Domain which will tell the ExtrudedDomain how to move the nodes when we want to do a parameter sweep in the cylinder oscillation amplitude...

The documentation for this class was generated from the following files: