Public Member Functions | Private Member Functions | Private Attributes | List of all members
StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT > Class Template Reference

Public Member Functions

 StabilityProblem (const unsigned &base_n_r, const unsigned &base_n_z, const unsigned &perturbed_n_r, const unsigned &perturbed_n_z, const double &domain_height)
 Constructor: Build base and perturbed state problems. More...
 
 ~StabilityProblem ()
 Destructor (empty) More...
 
void set_initial_condition ()
 Set initial conditions. More...
 
void unsteady_run (const double &dt, const unsigned &n_timesteps, DocInfo &doc_info)
 Integrate forwards in time with timestep dt for n_timesteps. More...
 
unsigned perform_power_method (const double &dt, const unsigned &n_timesteps, DocInfo &doc_info, const double &tolerance, const unsigned &max_iter, double &calc_eigenvalue, DoubleVector &input)
 
BaseStateProblem< BASE_ELEMENT > * base_state_problem_pt () const
 Access function for base state problem. More...
 
PerturbedStateProblem< BASE_ELEMENT, PERTURBED_ELEMENT > * perturbed_state_problem_pt () const
 Access function for perturbed state problem. More...
 
 StabilityProblem (const unsigned &base_n_r, const unsigned &base_n_z, const unsigned &perturbed_n_r, const unsigned &perturbed_n_z, const double &domain_height)
 Constructor: Build base and perturbed state problems. More...
 
 ~StabilityProblem ()
 Destructor (empty) More...
 
void set_initial_condition ()
 Set initial conditions. More...
 
void unsteady_run (const double &dt, const unsigned &n_timesteps, DocInfo &doc_info)
 Integrate forwards in time with timestep dt for n_timesteps. More...
 
unsigned perform_power_method (const double &dt, const unsigned &n_timesteps, DocInfo &doc_info, const double &tolerance, const unsigned &max_iter, double &calc_eigenvalue, DoubleVector &input)
 
BaseStateProblem< BASE_ELEMENT > * base_state_problem_pt () const
 Access function for base state problem. More...
 
PerturbedStateProblem< BASE_ELEMENT, PERTURBED_ELEMENT > * perturbed_state_problem_pt () const
 Access function for perturbed state problem. More...
 
 StabilityProblem (const unsigned &base_n_r, const unsigned &base_n_z, const unsigned &perturbed_n_r, const unsigned &perturbed_n_z, const double &domain_height)
 Constructor: Build base and perturbed state problems. More...
 
 ~StabilityProblem ()
 Destructor (empty) More...
 
void create_trace_files (DocInfo &doc_info)
 Create trace files. More...
 
void initialise_trace_files ()
 Initialise trace files (creates column headings) More...
 
void close_trace_files ()
 Clear and close trace files. More...
 
void initialise_dt (const double &dt)
 Initialise timestep (single dt version) More...
 
void initialise_dt (const Vector< double > &dt)
 Initialise timestep (vector of dts version) More...
 
void set_initial_condition ()
 Set initial conditions. More...
 
void doc_solution (DocInfo &doc_info, const bool &output_base_soln, const bool &output_perturbed_soln)
 Document base and perturbed state solutions. More...
 
void pass_updated_nondim_parameters_to_elements ()
 
unsigned perform_power_method (const double &dt, DocInfo &doc_info, const double &tolerance, const unsigned &max_iter, double &calc_eigenvalue, DoubleVector &input)
 
void base_flow_steady_newton_solve ()
 
void get_perturbed_state_problem_dofs (DoubleVector &d)
 Get perturbed state dofs. More...
 
void reset_global_time_to_zero ()
 (Re-)set the global time in both problems to zero More...
 
BaseStateProblem< BASE_ELEMENT > * base_state_problem_pt () const
 Access function for base state problem. More...
 
PerturbedStateProblem< BASE_ELEMENT, PERTURBED_ELEMENT > * perturbed_state_problem_pt () const
 Access function for perturbed state problem. More...
 
 StabilityProblem (const unsigned &base_n_r, const unsigned &base_n_z, const unsigned &perturbed_n_r, const unsigned &perturbed_n_z, const double &domain_height)
 Constructor: Build base and perturbed state problems. More...
 
 ~StabilityProblem ()
 Destructor (empty) More...
 
void create_trace_files (DocInfo &doc_info)
 Create trace files. More...
 
void initialise_trace_files ()
 Initialise trace files (creates column headings) More...
 
void close_trace_files ()
 Clear and close trace files. More...
 
void initialise_dt (const double &dt)
 Initialise timestep (single dt version) More...
 
void initialise_dt (const Vector< double > &dt)
 Initialise timestep (vector of dts version) More...
 
void set_initial_condition ()
 Set initial conditions. More...
 
void doc_solution (DocInfo &doc_info, const bool &output_base_soln, const bool &output_perturbed_soln)
 Document base and perturbed state solutions. More...
 
void pass_updated_nondim_parameters_to_elements ()
 
unsigned perform_power_method (const double &dt, DocInfo &doc_info, const double &tolerance, const unsigned &max_iter, double &calc_eigenvalue, DoubleVector &input)
 
void base_flow_steady_newton_solve ()
 
void get_perturbed_state_problem_dofs (DoubleVector &d)
 Get perturbed state dofs. More...
 
void reset_global_time_to_zero ()
 (Re-)set the global time in both problems to zero More...
 
BaseStateProblem< BASE_ELEMENT > * base_state_problem_pt () const
 Access function for base state problem. More...
 
PerturbedStateProblem< BASE_ELEMENT, PERTURBED_ELEMENT > * perturbed_state_problem_pt () const
 Access function for perturbed state problem. More...
 
 StabilityProblem (const unsigned &base_n_r, const unsigned &base_n_z, const unsigned &perturbed_n_r, const unsigned &perturbed_n_z, const double &radius_inner_cylinder, const double &radius_outer_cylinder, const double &l_z)
 Constructor: Build base and perturbed state problems. More...
 
 ~StabilityProblem ()
 Destructor (empty) More...
 
void create_trace_files (DocInfo &doc_info)
 Create trace files. More...
 
void initialise_trace_files ()
 Initialise trace files (creates column headings) More...
 
void close_trace_files ()
 Clear and close trace files. More...
 
void initialise_dt (const double &dt)
 Initialise timestep (single dt version) More...
 
void initialise_dt (const Vector< double > &dt)
 Initialise timestep (vector of dts version) More...
 
void set_initial_condition ()
 Set initial conditions. More...
 
void doc_solution (DocInfo &doc_info, const bool &output_base_soln, const bool &output_perturbed_soln)
 Document base and perturbed state solutions. More...
 
void unsteady_run (const double &dt, const unsigned &n_timesteps, DocInfo &doc_info, const bool &output_base_soln, Vector< DoubleVector > &perturbed_state_dofs_over_period, Vector< DoubleVector > &base_state_dofs_over_period)
 Integrate forwards in time with timestep dt for n_timesteps. More...
 
void unsteady_run_for_power_method (const double &dt, const unsigned &n_timesteps, DocInfo &doc_info, const DoubleVector &input_perturbed_state_dofs, DoubleVector &output_perturbed_state_dofs, const bool &output_base_soln, Vector< DoubleVector > &perturbed_state_dofs_over_period, Vector< DoubleVector > &base_state_dofs_over_period)
 
unsigned set_up_periodic_base_flow (const double &dt, const unsigned &n_timesteps, DocInfo &doc_info, const double &tolerance)
 
unsigned perform_power_method (const double &dt, const unsigned &n_timesteps, DocInfo &doc_info, const double &tolerance, const unsigned &max_iter, double &calc_eigenvalue, DoubleVector &input, Vector< DoubleVector > &base_state_dofs_over_period)
 
void get_perturbed_state_problem_dofs (DoubleVector &d)
 Get perturbed state dofs. More...
 
void reset_global_time_to_zero ()
 (Re-)set the global time in both problems to zero More...
 
const BaseStateProblem< BASE_ELEMENT > * base_state_problem_pt ()
 Access function for base state problem. More...
 
const PerturbedStateProblem< BASE_ELEMENT, PERTURBED_ELEMENT > * perturbed_state_problem_pt ()
 Access function for perturbed state problem. More...
 
void set_first_power_method_iteration_flag_to_true ()
 Switch on flag. More...
 
void set_first_power_method_iteration_flag_to_false ()
 Switch off flag. More...
 

Private Member Functions

void set_perturbed_state_time_equal_to_base_state_time ()
 
void set_perturbed_state_time_equal_to_base_state_time ()
 
void set_perturbed_state_time_equal_to_base_state_time ()
 

Private Attributes

BaseStateProblem< BASE_ELEMENT > * Base_state_problem_pt
 Pointer to base state problem class. More...
 
PerturbedStateProblem< BASE_ELEMENT, PERTURBED_ELEMENT > * Perturbed_state_problem_pt
 Pointer to perturbed state problem class. More...
 
ofstream Trace_file_periodic_base_flow
 
ofstream Trace_file_power_method
 Trace file which documents the convergence of the power method. More...
 
bool First_power_method_iteration
 

Detailed Description

template<class BASE_ELEMENT, class PERTURBED_ELEMENT>
class StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >

////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// Container class for the Tuckerman counter-rotating lids stability problem

////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////// Container class for the time-periodic Taylor–Couette stability problem

Constructor & Destructor Documentation

◆ StabilityProblem() [1/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::StabilityProblem ( const unsigned &  base_n_r,
const unsigned &  base_n_z,
const unsigned &  perturbed_n_r,
const unsigned &  perturbed_n_z,
const double &  domain_height 
)
inline

Constructor: Build base and perturbed state problems.

712  {
713  // Build base state problem
714  Base_state_problem_pt = new BaseStateProblem<BASE_ELEMENT>
715  (base_n_r,base_n_z,domain_height);
716 
717  // Build perturbed state problem
718  Perturbed_state_problem_pt = new PerturbedStateProblem
719  <BASE_ELEMENT,PERTURBED_ELEMENT>
720  (perturbed_n_r,perturbed_n_z,domain_height,
721  Base_state_problem_pt->mesh_pt());
722  }
BaseStateProblem
Base state problem class for Tuckerman counter-rotating lids problem.
Definition: axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:93
PerturbedStateProblem
Definition: axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:371
StabilityProblem::Perturbed_state_problem_pt
PerturbedStateProblem< BASE_ELEMENT, PERTURBED_ELEMENT > * Perturbed_state_problem_pt
Pointer to perturbed state problem class.
Definition: axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:770
StabilityProblem::Base_state_problem_pt
BaseStateProblem< BASE_ELEMENT > * Base_state_problem_pt
Pointer to base state problem class.
Definition: axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:766

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ ~StabilityProblem() [1/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::~StabilityProblem ( )
inline

Destructor (empty)

725 {}

◆ StabilityProblem() [2/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::StabilityProblem ( const unsigned &  base_n_r,
const unsigned &  base_n_z,
const unsigned &  perturbed_n_r,
const unsigned &  perturbed_n_z,
const double &  domain_height 
)
inline

Constructor: Build base and perturbed state problems.

773  {
774  // Build base state problem
775  Base_state_problem_pt = new BaseStateProblem<BASE_ELEMENT>
776  (base_n_r,base_n_z,domain_height);
777 
778  // Build perturbed state problem
779  Perturbed_state_problem_pt = new PerturbedStateProblem
780  <BASE_ELEMENT,PERTURBED_ELEMENT>
781  (perturbed_n_r,perturbed_n_z,domain_height,
782  Base_state_problem_pt->mesh_pt());
783  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ ~StabilityProblem() [2/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::~StabilityProblem ( )
inline

Destructor (empty)

786 {}

◆ StabilityProblem() [3/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::StabilityProblem ( const unsigned &  base_n_r,
const unsigned &  base_n_z,
const unsigned &  perturbed_n_r,
const unsigned &  perturbed_n_z,
const double &  domain_height 
)
inline

Constructor: Build base and perturbed state problems.

870  {
871  // Build base state problem
872  Base_state_problem_pt = new BaseStateProblem<BASE_ELEMENT>
873  (base_n_r,base_n_z,domain_height);
874 
875  // Build perturbed state problem
876  Perturbed_state_problem_pt = new PerturbedStateProblem
877  <BASE_ELEMENT,PERTURBED_ELEMENT>
878  (perturbed_n_r,perturbed_n_z,domain_height,
879  Base_state_problem_pt->mesh_pt());
880  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ ~StabilityProblem() [3/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::~StabilityProblem ( )
inline

Destructor (empty)

883 {}

◆ StabilityProblem() [4/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::StabilityProblem ( const unsigned &  base_n_r,
const unsigned &  base_n_z,
const unsigned &  perturbed_n_r,
const unsigned &  perturbed_n_z,
const double &  domain_height 
)
inline

Constructor: Build base and perturbed state problems.

928  {
929  // Build base state problem
930  Base_state_problem_pt = new BaseStateProblem<BASE_ELEMENT>
931  (base_n_r,base_n_z,domain_height);
932 
933  // Build perturbed state problem
934  Perturbed_state_problem_pt = new PerturbedStateProblem
935  <BASE_ELEMENT,PERTURBED_ELEMENT>
936  (perturbed_n_r,perturbed_n_z,domain_height,
937  Base_state_problem_pt->mesh_pt());
938  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ ~StabilityProblem() [4/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::~StabilityProblem ( )
inline

Destructor (empty)

941 {}

◆ StabilityProblem() [5/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::StabilityProblem ( const unsigned &  base_n_r,
const unsigned &  base_n_z,
const unsigned &  perturbed_n_r,
const unsigned &  perturbed_n_z,
const double &  radius_inner_cylinder,
const double &  radius_outer_cylinder,
const double &  l_z 
)
inline

Constructor: Build base and perturbed state problems.

875  {
876  // Build base state problem
877  Base_state_problem_pt = new BaseStateProblem<BASE_ELEMENT>
878  (base_n_r,base_n_z,radius_inner_cylinder,radius_outer_cylinder,l_z);
879 
880  // Build perturbed state problem
881  Perturbed_state_problem_pt =
882  new PerturbedStateProblem<BASE_ELEMENT,PERTURBED_ELEMENT>
883  (perturbed_n_r,perturbed_n_z,radius_inner_cylinder,
884  radius_outer_cylinder,l_z,Base_state_problem_pt->mesh_pt());
885  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ ~StabilityProblem() [5/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::~StabilityProblem ( )
inline

Destructor (empty)

888 {}

Member Function Documentation

◆ base_flow_steady_newton_solve() [1/2]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::base_flow_steady_newton_solve ( )
inline
981  {
982  Base_state_problem_pt->set_boundary_conditions();
983  Base_state_problem_pt->steady_newton_solve();
984  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt.

◆ base_flow_steady_newton_solve() [2/2]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::base_flow_steady_newton_solve ( )
inline
1039  {
1040  Base_state_problem_pt->set_boundary_conditions();
1041  Base_state_problem_pt->steady_newton_solve();
1042  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt.

◆ base_state_problem_pt() [1/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
const BaseStateProblem<BASE_ELEMENT>* StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::base_state_problem_pt ( )
inline

Access function for base state problem.

1045  { return Base_state_problem_pt; }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt.

◆ base_state_problem_pt() [2/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
BaseStateProblem<BASE_ELEMENT>* StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::base_state_problem_pt ( ) const
inline

Access function for base state problem.

757 { return Base_state_problem_pt; }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt.

◆ base_state_problem_pt() [3/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
BaseStateProblem<BASE_ELEMENT>* StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::base_state_problem_pt ( ) const
inline

Access function for base state problem.

818 { return Base_state_problem_pt; }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt.

◆ base_state_problem_pt() [4/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
BaseStateProblem<BASE_ELEMENT>* StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::base_state_problem_pt ( ) const
inline

Access function for base state problem.

1003 { return Base_state_problem_pt; }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt.

◆ base_state_problem_pt() [5/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
BaseStateProblem<BASE_ELEMENT>* StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::base_state_problem_pt ( ) const
inline

Access function for base state problem.

1061 { return Base_state_problem_pt; }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt.

◆ close_trace_files() [1/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::close_trace_files ( )
inline

Clear and close trace files.

915  {
916  // Initialise base and perturbed state trace files
917  Base_state_problem_pt->close_trace_file();
918  Perturbed_state_problem_pt->close_trace_file();
919 
920  // Initialise stability problem trace files
921  Trace_file_power_method.clear();
922  Trace_file_power_method.close();
923  }
StabilityProblem::Trace_file_power_method
ofstream Trace_file_power_method
Trace file which documents the convergence of the power method.
Definition: linearised_axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:1030

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Trace_file_power_method.

◆ close_trace_files() [2/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::close_trace_files ( )
inline

Clear and close trace files.

973  {
974  // Initialise base and perturbed state trace files
975  Base_state_problem_pt->close_trace_file();
976  Perturbed_state_problem_pt->close_trace_file();
977 
978  // Initialise stability problem trace files
979  Trace_file_power_method.clear();
980  Trace_file_power_method.close();
981  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Trace_file_power_method.

◆ close_trace_files() [3/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::close_trace_files ( )
inline

Clear and close trace files.

927  {
928  // Initialise base and perturbed state trace files
929  Base_state_problem_pt->close_trace_file();
930  Perturbed_state_problem_pt->close_trace_file();
931 
932  // Initialise stability problem trace files
933  Trace_file_periodic_base_flow.clear();
934  Trace_file_periodic_base_flow.close();
935  Trace_file_power_method.clear();
936  Trace_file_power_method.close();
937  }
StabilityProblem::Trace_file_periodic_base_flow
ofstream Trace_file_periodic_base_flow
Definition: linearised_axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:1027

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt, StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Trace_file_periodic_base_flow, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Trace_file_power_method.

◆ create_trace_files() [1/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::create_trace_files ( DocInfo &  doc_info)
inline

Create trace files.

887  {
888  // Set up base and perturbed state trace files
889  Base_state_problem_pt->create_trace_file(doc_info);
890  Perturbed_state_problem_pt->create_trace_file(doc_info);
891 
892  // Set up stability problem trace files
893  char filename[256];
894  sprintf(filename,"%s/power_method_trace_k%i_Re%4.2f.dat",
895  doc_info.directory().c_str(),
896  GlobalPhysicalVariables::k,
897  GlobalPhysicalVariables::Re_current);
898  Trace_file_power_method.open(filename);
899  }
oomph::DocInfo::directory
std::string directory() const
Output directory.
Definition: oomph_utilities.h:524
GlobalPhysicalVariables::k
int k
Azimuthal mode number k in e^ik(theta) decomposition.
Definition: axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:73
GlobalPhysicalVariables::Re_current
double Re_current
Definition: linearised_axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:70
MergeRestartFiles.filename
string filename
Definition: MergeRestartFiles.py:39

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, oomph::DocInfo::directory(), MergeRestartFiles::filename, GlobalPhysicalVariables::k, StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt, GlobalPhysicalVariables::Re_current, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Trace_file_power_method.

◆ create_trace_files() [2/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::create_trace_files ( DocInfo &  doc_info)
inline

Create trace files.

945  {
946  // Set up base and perturbed state trace files
947  Base_state_problem_pt->create_trace_file(doc_info);
948  Perturbed_state_problem_pt->create_trace_file(doc_info);
949 
950  // Set up stability problem trace files
951  char filename[256];
952  sprintf(filename,"%s/power_method_trace_k%i_Re%4.2f.dat",
953  doc_info.directory().c_str(),
954  GlobalPhysicalVariables::k,
955  GlobalPhysicalVariables::Re_current);
956  Trace_file_power_method.open(filename);
957  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, oomph::DocInfo::directory(), MergeRestartFiles::filename, GlobalPhysicalVariables::k, StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt, GlobalPhysicalVariables::Re_current, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Trace_file_power_method.

◆ create_trace_files() [3/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::create_trace_files ( DocInfo &  doc_info)
inline

Create trace files.

892  {
893  // Set up base and perturbed state trace files
894  Base_state_problem_pt->create_trace_file(doc_info);
895  Perturbed_state_problem_pt->create_trace_file(doc_info);
896 
897  // Set up stability problem trace files
898  char filename[256];
899  sprintf(filename,"%s/periodic_base_flow_trace_epsilon%2.1f_Re%4.2f.dat",
900  doc_info.directory().c_str(),
901  GlobalPhysicalVariables::Epsilon,
902  GlobalPhysicalVariables::MMC_Re_current);
903  Trace_file_periodic_base_flow.open(filename);
904 
905  sprintf(filename,"%s/power_method_trace_epsilon%2.1f_Re%4.2f.dat",
906  doc_info.directory().c_str(),
907  GlobalPhysicalVariables::Epsilon,
908  GlobalPhysicalVariables::MMC_Re_current);
909  Trace_file_power_method.open(filename);
910  }
GlobalPhysicalVariables::Epsilon
double Epsilon
Dimensionless modulation amplitude (epsilon)
Definition: time_periodic_taylor_couette.cc:84
GlobalPhysicalVariables::MMC_Re_current
double MMC_Re_current
Definition: time_periodic_taylor_couette.cc:73

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, oomph::DocInfo::directory(), GlobalPhysicalVariables::Epsilon, MergeRestartFiles::filename, GlobalPhysicalVariables::MMC_Re_current, StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt, StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Trace_file_periodic_base_flow, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Trace_file_power_method.

◆ doc_solution() [1/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::doc_solution ( DocInfo &  doc_info,
const bool &  output_base_soln,
const bool &  output_perturbed_soln 
)
inline

Document base and perturbed state solutions.

949  {
950  Base_state_problem_pt->doc_solution(doc_info,output_base_soln);
951  Perturbed_state_problem_pt->doc_solution(doc_info,output_perturbed_soln);
952  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

Referenced by StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::doc_solution().

◆ doc_solution() [2/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::doc_solution ( DocInfo &  doc_info,
const bool &  output_base_soln,
const bool &  output_perturbed_soln 
)
inline

Document base and perturbed state solutions.

1007  {
1008  Base_state_problem_pt->doc_solution(doc_info,output_base_soln);
1009  Perturbed_state_problem_pt->doc_solution(doc_info,output_perturbed_soln);
1010  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ doc_solution() [3/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::doc_solution ( DocInfo &  doc_info,
const bool &  output_base_soln,
const bool &  output_perturbed_soln 
)
inline

Document base and perturbed state solutions.

964  {
965  Base_state_problem_pt->
966  doc_solution(doc_info,output_base_soln);
967  Perturbed_state_problem_pt->doc_solution(doc_info,output_perturbed_soln);
968  }
StabilityProblem::doc_solution
void doc_solution(DocInfo &doc_info, const bool &output_base_soln, const bool &output_perturbed_soln)
Document base and perturbed state solutions.
Definition: linearised_axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:947

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::doc_solution(), and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ get_perturbed_state_problem_dofs() [1/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::get_perturbed_state_problem_dofs ( DoubleVector &  d)
inline

Get perturbed state dofs.

988  {
989  Perturbed_state_problem_pt->get_dofs(d);
990  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ get_perturbed_state_problem_dofs() [2/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::get_perturbed_state_problem_dofs ( DoubleVector &  d)
inline

Get perturbed state dofs.

1046  {
1047  Perturbed_state_problem_pt->get_dofs(d);
1048  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ get_perturbed_state_problem_dofs() [3/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::get_perturbed_state_problem_dofs ( DoubleVector &  d)
inline

Get perturbed state dofs.

1030  {
1031  Perturbed_state_problem_pt->get_dofs(d);
1032  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ initialise_dt() [1/6]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::initialise_dt ( const double &  dt)
inline

Initialise timestep (single dt version)

927  {
928  Base_state_problem_pt->initialise_dt(dt);
929  Perturbed_state_problem_pt->initialise_dt(dt);
930  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ initialise_dt() [2/6]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::initialise_dt ( const double &  dt)
inline

Initialise timestep (single dt version)

985  {
986  Base_state_problem_pt->initialise_dt(dt);
987  Perturbed_state_problem_pt->initialise_dt(dt);
988  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ initialise_dt() [3/6]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::initialise_dt ( const double &  dt)
inline

Initialise timestep (single dt version)

941  {
942  Base_state_problem_pt->initialise_dt(dt);
943  Perturbed_state_problem_pt->initialise_dt(dt);
944  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ initialise_dt() [4/6]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::initialise_dt ( const Vector< double > &  dt)
inline

Initialise timestep (vector of dts version)

934  {
935  Base_state_problem_pt->initialise_dt(dt);
936  Perturbed_state_problem_pt->initialise_dt(dt);
937  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ initialise_dt() [5/6]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::initialise_dt ( const Vector< double > &  dt)
inline

Initialise timestep (vector of dts version)

992  {
993  Base_state_problem_pt->initialise_dt(dt);
994  Perturbed_state_problem_pt->initialise_dt(dt);
995  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ initialise_dt() [6/6]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::initialise_dt ( const Vector< double > &  dt)
inline

Initialise timestep (vector of dts version)

948  {
949  Base_state_problem_pt->initialise_dt(dt);
950  Perturbed_state_problem_pt->initialise_dt(dt);
951  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ initialise_trace_files() [1/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::initialise_trace_files ( )
inline

Initialise trace files (creates column headings)

903  {
904  // Initialise base and perturbed state trace files
905  Base_state_problem_pt->initialise_trace_file();
906  Perturbed_state_problem_pt->initialise_trace_file();
907 
908  // Initialise stability problem trace files
909  Trace_file_power_method << "time, convergence_criterion, "
910  << "tolerance*abs(eigenvalue)" << std::endl;
911  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Trace_file_power_method.

◆ initialise_trace_files() [2/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::initialise_trace_files ( )
inline

Initialise trace files (creates column headings)

961  {
962  // Initialise base and perturbed state trace files
963  Base_state_problem_pt->initialise_trace_file();
964  Perturbed_state_problem_pt->initialise_trace_file();
965 
966  // Initialise stability problem trace files
967  Trace_file_power_method << "time, convergence_criterion, "
968  << "tolerance*abs(eigenvalue)" << std::endl;
969  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Trace_file_power_method.

◆ initialise_trace_files() [3/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::initialise_trace_files ( )
inline

Initialise trace files (creates column headings)

914  {
915  // Initialise base and perturbed state trace files
916  Base_state_problem_pt->initialise_trace_file();
917  Perturbed_state_problem_pt->initialise_trace_file();
918 
919  // Initialise stability problem trace files
920  Trace_file_periodic_base_flow << "time, norm, tolerance" << std::endl;
921  Trace_file_power_method << "time, convergence_criterion, "
922  << "tolerance*abs(eigenvalue)" << std::endl;
923  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt, StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Trace_file_periodic_base_flow, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Trace_file_power_method.

◆ pass_updated_nondim_parameters_to_elements() [1/2]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::pass_updated_nondim_parameters_to_elements ( )
inline

Pass updated nondimensional parameters to the elements of both problems

957  {
958  Base_state_problem_pt->pass_updated_nondim_parameters_to_elements();
959  Perturbed_state_problem_pt->pass_updated_nondim_parameters_to_elements();
960  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ pass_updated_nondim_parameters_to_elements() [2/2]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::pass_updated_nondim_parameters_to_elements ( )
inline

Pass updated nondimensional parameters to the elements of both problems

1015  {
1016  Base_state_problem_pt->pass_updated_nondim_parameters_to_elements();
1017  Perturbed_state_problem_pt->pass_updated_nondim_parameters_to_elements();
1018  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ perform_power_method() [1/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
unsigned StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::perform_power_method ( const double &  dt,
const unsigned &  n_timesteps,
DocInfo &  doc_info,
const double &  tolerance,
const unsigned &  max_iter,
double &  calc_eigenvalue,
DoubleVector &  input 
)

Perform a power method with a maximum of max_iter iterations. This function requires an initial guess for the dominant eigenvector "input"and a tolerance which controls at which point the dominant eigenvalue will be deemed converged. At this point the calculated value of this eigenvalue "calc_eigenvalue" and its corresponding eigenvector "input" will be returned. The return value of this function is the number of power method iterations which were performed.

Perform a power method with a maximum of max_iter iterations. This function requires an initial guess for the dominant eigenvector "input" and a tolerance which controls at which point the dominant eigenvalue will be deemed converged. At this point the calculated value of this eigenvalue "calc_eigenvalue" and its corresponding eigenvector "input" will be returned. The return value of this function is the number of power method iterations which were performed.

820 {
821  // Determine number of degrees of freedom in perturbed state problem
822  const unsigned n_dof_perturbed = Perturbed_state_problem_pt->ndof();
823 
824  // Begin power method loop
825  for(unsigned p=0;p<max_iter;p++)
826  {
827  cout << "\nThis is power method loop number " << p+1
828  << " of a maximum " << max_iter << "." << endl;
829 
830  // Normalise input vector
831  const double inv_norm = 1.0/input.norm();
832  for(unsigned i=0;i<n_dof_perturbed;i++) { input[i] *= inv_norm; }
833 
834  // Set up storage for output vector
835  DoubleVector output;
836 
837  // Re-set the starting dof vector for the perturbed state problem
838  Perturbed_state_problem_pt->set_dofs(input);
839 
840  // Perform unsteady run for n_timesteps
841  this->unsteady_run(dt,n_timesteps,doc_info);
842 
843  // Get dofs from perturbed state problem
844  Perturbed_state_problem_pt->get_dofs(output);
845 
846  // Calculate eigenvalue (scalar product of input and output vectors)
847  calc_eigenvalue = input.dot(output);
848 
849  // Calculate L2 norm of (output - eigenvalue*input)
850  // divided by the square root of the number of degrees of freedom
851  double convergence_criterion = 0.0;
852  for(unsigned i=0;i<n_dof_perturbed;i++)
853  {
854  convergence_criterion +=
855  (output[i] - calc_eigenvalue*input[i])
856  *(output[i] - calc_eigenvalue*input[i]);
857  }
858  convergence_criterion = sqrt(convergence_criterion/n_dof_perturbed);
859 
860  // Print convergence criterion and threshold
861  cout << "\nConvergence criterion = " << convergence_criterion
862  << "\ntolerance*abs(eigenvalue) = " << tolerance*abs(calc_eigenvalue)
863  << endl;
864 
865  // Check for convergence
866  if(convergence_criterion <= tolerance*abs(calc_eigenvalue))
867  {
868  cout << "\nPower method has converged to within a tolerance of "
869  << tolerance << "\nin " << p+1
870  << " iterations of the power method.\n" << endl;
871 
872  // Output the dominant eigenvector
873  Perturbed_state_problem_pt->doc_solution(doc_info);
874 
875  // Return the number of iterations of the power method
876  return p+1;
877  }
878  else
879  {
880  input = output;
881  output.clear();
882  }
883  }
884 
885  // If we reach here, the power method has failed to converge in max_iter
886  // iterations.
887  cout << "\nPower method has failed to converge to within a tolerance of "
888  << tolerance << "\nin " << max_iter
889  << " iterations of the power method" << endl;
890 
891  // Output the dominant eigenvector
892  Perturbed_state_problem_pt->doc_solution(doc_info);
893 
894  return max_iter;
895 
896 } // End of perform_power_method
abs
AnnoyingScalar abs(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:135
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
p
float * p
Definition: Tutorial_Map_using.cpp:9
StabilityProblem::unsteady_run
void unsteady_run(const double &dt, const unsigned &n_timesteps, DocInfo &doc_info)
Integrate forwards in time with timestep dt for n_timesteps.
Definition: axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:781
oomph::DoubleVector
Definition: double_vector.h:58
oomph::DoubleVector::norm
double norm() const
compute the 2 norm of this vector
Definition: double_vector.cc:867
oomph::DoubleVector::dot
double dot(const DoubleVector &vec) const
compute the dot product of this vector with the vector vec.
Definition: double_vector.cc:805
output
void output(std::ostream &outfile, const unsigned &nplot)
Overload output function.
Definition: overloaded_element_body.h:490

References abs(), oomph::DoubleVector::dot(), i, oomph::DoubleVector::norm(), output(), p, and sqrt().

◆ perform_power_method() [2/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
unsigned StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::perform_power_method ( const double &  dt,
const unsigned &  n_timesteps,
DocInfo &  doc_info,
const double &  tolerance,
const unsigned &  max_iter,
double &  calc_eigenvalue,
DoubleVector &  input 
)

Perform a power method with a maximum of max_iter iterations. This function requires an initial guess for the dominant eigenvector "input"and a tolerance which controls at which point the dominant eigenvalue will be deemed converged. At this point the calculated value of this eigenvalue "calc_eigenvalue" and its corresponding eigenvector "input" will be returned. The return value of this function is the number of power method iterations which were performed.

◆ perform_power_method() [3/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
unsigned StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::perform_power_method ( const double &  dt,
const unsigned &  n_timesteps,
DocInfo &  doc_info,
const double &  tolerance,
const unsigned &  max_iter,
double &  calc_eigenvalue,
DoubleVector &  input,
Vector< DoubleVector > &  base_state_dofs_over_period 
)

Perform a power method with a maximum of max_iter iterations. This function requires an initial guess for the dominant eigenvector "input"and a tolerance which controls at which point the dominant eigenvalue will be deemed converged. At this point the calculated value of this eigenvalue "calc_eigenvalue" and its corresponding eigenvector "input" will be returned. The return value of this function is the number of power method iterations which were performed.

1494 {
1495  // Initialise output_base_state_solution flag to true
1496  bool output_base_state_solution = true;
1497 
1498  // Reset the solution counter to zero
1499  doc_info.number()=0;
1500 
1501  // Determine number of degrees of freedom in perturbed state problem
1502  const unsigned n_dof_perturbed = Perturbed_state_problem_pt->ndof();
1503 
1504  // Set up external storage for the dofs over one period
1505  Vector<DoubleVector> dofs_over_period(n_timesteps);
1506 
1507  // Set flag to true
1508  this->set_first_power_method_iteration_flag_to_true();
1509 
1510  // Begin power method loop
1511  for(unsigned p=0;p<max_iter;p++)
1512  {
1513  cout << "\nThis is power method loop number " << p+1
1514  << " of a maximum " << max_iter << "." << endl;
1515 
1516  // Normalise input vector
1517  const double inv_norm = 1.0/input.norm();
1518  for(unsigned i=0;i<n_dof_perturbed;i++) { input[i] *= inv_norm; }
1519 
1520  // Set up storage for output vector from Arnoldi method
1521  DoubleVector output;
1522 
1523  // Perform unsteady run over this period
1524  unsteady_run_for_power_method(dt,n_timesteps,doc_info,input,output,
1525  output_base_state_solution,dofs_over_period,
1526  base_state_dofs_over_period);
1527 
1528  // Calculate eigenvalue (scalar product of input and output vectors)
1529  calc_eigenvalue = input.dot(output);
1530 
1531  // Calculate L2 norm of (output - eigenvalue*input)
1532  double convergence_criterion = 0.0;
1533  for(unsigned i=0;i<n_dof_perturbed;i++)
1534  {
1535  convergence_criterion +=
1536  (output[i] - calc_eigenvalue*input[i])
1537  *(output[i] - calc_eigenvalue*input[i]);
1538  }
1539  convergence_criterion = sqrt(convergence_criterion);
1540 
1541  cout << "\nConvergence criterion = " << convergence_criterion << endl;
1542  cout << "tolerance*abs(eigenvalue) = " << tolerance*abs(calc_eigenvalue)
1543  << endl;
1544 
1545  // Document convergence criterion and tolerance*eigenvalue in trace file
1546  Trace_file_power_method << Perturbed_state_problem_pt->time_pt()->time()
1547  << " " << convergence_criterion << " "
1548  << tolerance*abs(calc_eigenvalue) << std::endl;
1549 
1550  // We only wish to output the base state solution during the first
1551  // iteration of the power method (so that we can see what it looks like
1552  // when converged)
1553  output_base_state_solution = false;
1554 
1555  // Check for convergence
1556  if(convergence_criterion <= tolerance*abs(calc_eigenvalue))
1557  {
1558  cout << "\nPower method has converged to within a tolerance of "
1559  << tolerance << "\nin " << p+1
1560  << " iterations of the power method.\n" << endl;
1561 
1562  // Reset the solution counter to zero
1563  doc_info.number()=0;
1564 
1565  // Output the dominant eigenvector over an entire period
1566  for(unsigned i=0;i<n_timesteps;i++)
1567  {
1568  // Populate perturbed state problem with dofs which have been stored
1569  // over the last period
1570  Perturbed_state_problem_pt->
1571  assign_eigenvector_to_dofs(dofs_over_period[i]);
1572 
1573  // Output the perturbed state solution
1574  this->doc_solution(doc_info,false,true);
1575 
1576  // Increment counter for solutions
1577  doc_info.number()++;
1578  }
1579 
1580  // Return the number of iterations of the power method
1581  return p+1;
1582  }
1583  else
1584  {
1585  input = output;
1586  output.clear();
1587  }
1588 
1589  // Set flag to false
1590  this->set_first_power_method_iteration_flag_to_false();
1591  }
1592 
1593  // If we reach here, the power method has failed to converge in max_iter
1594  // iterations.
1595  cout << "\nPower method has failed to converge to within a tolerance of "
1596  << tolerance << "\nin " << max_iter
1597  << " iterations of the power method" << endl;
1598 
1599  // Reset the solution counter to zero
1600  doc_info.number()=0;
1601 
1602  // Output the dominant eigenvector over an entire period
1603  for(unsigned i=0;i<n_timesteps;i++)
1604  {
1605  // Populate perturbed state problem with dofs which have been stored
1606  // over the last period
1607  Perturbed_state_problem_pt->
1608  assign_eigenvector_to_dofs(dofs_over_period[i]);
1609 
1610  // Output the perturbed state solution
1611  this->doc_solution(doc_info,false,true);
1612 
1613  // Increment counter for solutions
1614  doc_info.number()++;
1615  }
1616 
1617  return max_iter;
1618 
1619 } // End of perform_power_method
StabilityProblem::unsteady_run_for_power_method
void unsteady_run_for_power_method(const double &dt, const unsigned &n_timesteps, DocInfo &doc_info, const DoubleVector &input_perturbed_state_dofs, DoubleVector &output_perturbed_state_dofs, const bool &output_base_soln, Vector< DoubleVector > &perturbed_state_dofs_over_period, Vector< DoubleVector > &base_state_dofs_over_period)
Definition: time_periodic_taylor_couette.cc:984
StabilityProblem::set_first_power_method_iteration_flag_to_false
void set_first_power_method_iteration_flag_to_false()
Switch off flag.
Definition: time_periodic_taylor_couette.cc:1056
StabilityProblem::set_first_power_method_iteration_flag_to_true
void set_first_power_method_iteration_flag_to_true()
Switch on flag.
Definition: time_periodic_taylor_couette.cc:1052
oomph::DocInfo::number
unsigned & number()
Number used (e.g.) for labeling output files.
Definition: oomph_utilities.h:554
oomph::Vector
Definition: oomph-lib/src/generic/Vector.h:58

References abs(), oomph::DoubleVector::dot(), i, oomph::DoubleVector::norm(), oomph::DocInfo::number(), output(), p, and sqrt().

◆ perform_power_method() [4/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
unsigned StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::perform_power_method ( const double &  dt,
DocInfo &  doc_info,
const double &  tolerance,
const unsigned &  max_iter,
double &  calc_eigenvalue,
DoubleVector &  input 
)

Perform a power method with a maximum of max_iter iterations. This function requires an initial guess for the dominant eigenvector "input"and a tolerance which controls at which point the dominant eigenvalue will be deemed converged. At this point the calculated value of this eigenvalue "calc_eigenvalue" and its corresponding eigenvector "input" will be returned. The return value of this function is the number of power method iterations which were performed.

1061 {
1062  // Initialise output_base_state_solution flag to true
1063  bool output_base_state_solution = true;
1064 
1065  // Reset the solution counter to zero
1066  doc_info.number()=0;
1067 
1068  // Determine number of degrees of freedom in perturbed state problem
1069  const unsigned n_dof_perturbed = Perturbed_state_problem_pt->ndof();
1070 
1071  // Begin power method loop
1072  for(unsigned p=0;p<max_iter;p++)
1073  {
1074  cout << "\nThis is power method loop number " << p+1
1075  << " of a maximum " << max_iter << "." << endl;
1076 
1077  // Normalise input vector
1078  const double inv_norm = 1.0/input.norm();
1079  for(unsigned i=0;i<n_dof_perturbed;i++) { input[i] *= inv_norm; }
1080 
1081  // Set up storage for output vector
1082  DoubleVector output;
1083 
1084  // Re-set the starting dof vector for the perturbed state problem
1085  Perturbed_state_problem_pt->set_dofs(input);
1086 
1087  // Set the perturbed state problem's timestepper to work in BDF1 mode
1088  Perturbed_state_problem_pt->time_stepper_pt()->turn_on_bdf1_mode();
1089 
1090  // Initialise timestep (also sets timestepper weights)
1091  this->initialise_dt(dt);
1092 
1093  // Take timestep
1094  Perturbed_state_problem_pt->unsteady_newton_solve(dt);
1095 
1096  // Doc info in trace file and output base state solution
1097  // (if flag is set appropriately)
1098  this->doc_solution(doc_info,output_base_state_solution,false);
1099 
1100  // Increment counter for solutions
1101  doc_info.number()++;
1102 
1103  // Get dofs from perturbed state problem
1104  Perturbed_state_problem_pt->get_dofs(output);
1105 
1106  // Calculate eigenvalue (scalar product of input and output vectors)
1107  calc_eigenvalue = input.dot(output);
1108 
1109  // Calculate L2 norm of (output - eigenvalue*input)
1110  // divided by the square root of the number of degrees of freedom
1111  double convergence_criterion = 0.0;
1112  for(unsigned i=0;i<n_dof_perturbed;i++)
1113  {
1114  convergence_criterion +=
1115  (output[i] - calc_eigenvalue*input[i])
1116  *(output[i] - calc_eigenvalue*input[i]);
1117  }
1118  convergence_criterion = sqrt(convergence_criterion/n_dof_perturbed);
1119 
1120  // Print convergence criterion and threshold
1121  cout << "\nConvergence criterion = " << convergence_criterion << endl;
1122  cout << "tolerance*abs(eigenvalue) = " << tolerance*abs(calc_eigenvalue)
1123  << endl;
1124 
1125  // Document convergence criterion and threshold in trace file
1126  Trace_file_power_method << Perturbed_state_problem_pt->time_pt()->time()
1127  << " " << convergence_criterion << " "
1128  << tolerance*abs(calc_eigenvalue) << std::endl;
1129 
1130  // We only wish to output the base state solution during the first
1131  // iteration of the power method (so that we can see what it looks like
1132  // when converged)
1133  output_base_state_solution = false;
1134 
1135  // Check for convergence
1136  if(convergence_criterion <= tolerance*abs(calc_eigenvalue))
1137  {
1138  cout << "\nPower method has converged to within a tolerance of "
1139  << tolerance << "\nin " << p+1
1140  << " iterations of the power method.\n" << endl;
1141 
1142  // Reset the solution counter to zero
1143  doc_info.number()=0;
1144 
1145  // Output the dominant eigenvector (perturbed state solution)
1146  this->doc_solution(doc_info,false,true);
1147 
1148  // Return the number of iterations of the power method
1149  return p+1;
1150  }
1151  else
1152  {
1153  input = output;
1154  output.clear();
1155  }
1156  }
1157 
1158  // If we reach here, the power method has failed to converge in max_iter
1159  // iterations.
1160  cout << "\nPower method has failed to converge to within a tolerance of "
1161  << tolerance << "\nin " << max_iter
1162  << " iterations of the power method" << endl;
1163 
1164  // Reset the solution counter to zero
1165  doc_info.number()=0;
1166 
1167  // Output the dominant eigenvector (perturbed state solution)
1168  this->doc_solution(doc_info,false,true);
1169 
1170  return max_iter;
1171 
1172 } // End of perform_power_method
StabilityProblem::initialise_dt
void initialise_dt(const double &dt)
Initialise timestep (single dt version)
Definition: linearised_axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:926

References abs(), oomph::DoubleVector::dot(), i, oomph::DoubleVector::norm(), oomph::DocInfo::number(), output(), p, and sqrt().

◆ perform_power_method() [5/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
unsigned StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::perform_power_method ( const double &  dt,
DocInfo &  doc_info,
const double &  tolerance,
const unsigned &  max_iter,
double &  calc_eigenvalue,
DoubleVector &  input 
)

Perform a power method with a maximum of max_iter iterations. This function requires an initial guess for the dominant eigenvector "input"and a tolerance which controls at which point the dominant eigenvalue will be deemed converged. At this point the calculated value of this eigenvalue "calc_eigenvalue" and its corresponding eigenvector "input" will be returned. The return value of this function is the number of power method iterations which were performed.

◆ perturbed_state_problem_pt() [1/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
const PerturbedStateProblem<BASE_ELEMENT,PERTURBED_ELEMENT>* StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::perturbed_state_problem_pt ( )
inline

Access function for perturbed state problem.

1049 { return Perturbed_state_problem_pt; }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ perturbed_state_problem_pt() [2/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
PerturbedStateProblem<BASE_ELEMENT,PERTURBED_ELEMENT>* StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::perturbed_state_problem_pt ( ) const
inline

Access function for perturbed state problem.

761 { return Perturbed_state_problem_pt; }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ perturbed_state_problem_pt() [3/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
PerturbedStateProblem<BASE_ELEMENT,PERTURBED_ELEMENT>* StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::perturbed_state_problem_pt ( ) const
inline

Access function for perturbed state problem.

822 { return Perturbed_state_problem_pt; }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ perturbed_state_problem_pt() [4/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
PerturbedStateProblem<BASE_ELEMENT,PERTURBED_ELEMENT>* StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::perturbed_state_problem_pt ( ) const
inline

Access function for perturbed state problem.

1007 { return Perturbed_state_problem_pt; }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ perturbed_state_problem_pt() [5/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
PerturbedStateProblem<BASE_ELEMENT,PERTURBED_ELEMENT>* StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::perturbed_state_problem_pt ( ) const
inline

Access function for perturbed state problem.

1065 { return Perturbed_state_problem_pt; }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ reset_global_time_to_zero() [1/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::reset_global_time_to_zero ( )
inline

(Re-)set the global time in both problems to zero

994  {
995  Base_state_problem_pt->time_pt()->time() = 0.0;
996  Perturbed_state_problem_pt->time_pt()->time() = 0.0;
997 
998  cout << "\nReset global time to zero in both problems.\n" << endl;
999  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ reset_global_time_to_zero() [2/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::reset_global_time_to_zero ( )
inline

(Re-)set the global time in both problems to zero

1052  {
1053  Base_state_problem_pt->time_pt()->time() = 0.0;
1054  Perturbed_state_problem_pt->time_pt()->time() = 0.0;
1055 
1056  cout << "\nReset global time to zero in both problems.\n" << endl;
1057  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ reset_global_time_to_zero() [3/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::reset_global_time_to_zero ( )
inline

(Re-)set the global time in both problems to zero

1036  {
1037  Base_state_problem_pt->time_pt()->time() = 0.0;
1038  Perturbed_state_problem_pt->time_pt()->time() = 0.0;
1039 
1040  cout << "\nReset global time to zero in both problems.\n" << endl;
1041  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ set_first_power_method_iteration_flag_to_false()

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_first_power_method_iteration_flag_to_false ( )
inline

Switch off flag.

1057  { First_power_method_iteration = false; }
StabilityProblem::First_power_method_iteration
bool First_power_method_iteration
Definition: time_periodic_taylor_couette.cc:1091

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::First_power_method_iteration.

◆ set_first_power_method_iteration_flag_to_true()

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_first_power_method_iteration_flag_to_true ( )
inline

Switch on flag.

1053  { First_power_method_iteration = true; }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::First_power_method_iteration.

◆ set_initial_condition() [1/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_initial_condition ( )
inline

Set initial conditions.

729  {
730  Base_state_problem_pt->set_initial_condition();
731  Perturbed_state_problem_pt->set_initial_condition();
732  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ set_initial_condition() [2/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_initial_condition ( )
inline

Set initial conditions.

790  {
791  Base_state_problem_pt->set_initial_condition();
792  Perturbed_state_problem_pt->set_initial_condition();
793  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ set_initial_condition() [3/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_initial_condition ( )
inline

Set initial conditions.

941  {
942  Base_state_problem_pt->set_initial_condition();
943  Perturbed_state_problem_pt->set_initial_condition();
944  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ set_initial_condition() [4/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_initial_condition ( )
inline

Set initial conditions.

999  {
1000  Base_state_problem_pt->set_initial_condition();
1001  Perturbed_state_problem_pt->set_initial_condition();
1002  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ set_initial_condition() [5/5]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_initial_condition ( )
inline

Set initial conditions.

955  {
956  Base_state_problem_pt->set_initial_condition();
957  Perturbed_state_problem_pt->set_initial_condition();
958  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ set_perturbed_state_time_equal_to_base_state_time() [1/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_perturbed_state_time_equal_to_base_state_time ( )
inlineprivate

(Re-)set the global time in the perturbed state problem to be equal to the global time in the base state problem

1014  {
1015  // Get current global time according to base state problem
1016  const double base_state_time = Base_state_problem_pt->time_pt()->time();
1017 
1018  // Set perturbed state problem's global time
1019  Perturbed_state_problem_pt->time_pt()->time() = base_state_time;
1020 
1021  cout << "\nReset perturbed state problem's global time to "
1022  << base_state_time << endl;
1023  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ set_perturbed_state_time_equal_to_base_state_time() [2/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_perturbed_state_time_equal_to_base_state_time ( )
inlineprivate

(Re-)set the global time in the perturbed state problem to be equal to the global time in the base state problem

1072  {
1073  // Get current global time according to base state problem
1074  const double base_state_time = Base_state_problem_pt->time_pt()->time();
1075 
1076  // Set perturbed state problem's global time
1077  Perturbed_state_problem_pt->time_pt()->time() = base_state_time;
1078 
1079  cout << "\nReset perturbed state problem's global time to "
1080  << base_state_time << endl;
1081  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ set_perturbed_state_time_equal_to_base_state_time() [3/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_perturbed_state_time_equal_to_base_state_time ( )
inlineprivate

(Re-)set the global time in the perturbed state problem to be equal to the global time in the base state problem

1064  {
1065  // Get current global time according to base state problem
1066  const double base_state_time = Base_state_problem_pt->time_pt()->time();
1067 
1068  // Set perturbed state problem's global time
1069  Perturbed_state_problem_pt->time_pt()->time() = base_state_time;
1070 
1071  cout << "\nReset perturbed state problem's global time to "
1072  << base_state_time << endl;
1073  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt.

◆ set_up_periodic_base_flow()

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
unsigned StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_up_periodic_base_flow ( const double &  dt,
const unsigned &  n_timesteps,
DocInfo &  doc_info,
const double &  tolerance 
)

Integrate just the base state problem forwards in time with timestep dt for n_timesteps. Once a periodic state has been reached, set perturbed state problem's global time to be equal to the base state problem's global time and begin solving the perturbed state problem also. This function returns the number of periods taken to satisfy the "periodic state" exit criteria.

1364 {
1365  // Check that base state timestepper is working in BDF2 mode
1366  if(Base_state_problem_pt->time_stepper_pt()->bdf1_mode())
1367  {
1368  std::cout << "Base state problem's timestepper is not in BDF2 mode."
1369  << "\n Exiting...\n" << std::endl;
1370  exit(1);
1371  }
1372 
1373  // Initialise timestep (also sets timestepper weights)
1374  this->initialise_dt(dt);
1375 
1376  // Set up storage for base state problem's dof vectors at the end of the
1377  // previous period and at the end of the current period
1378  Vector<double> previous_dofs;
1379  Vector<double> current_dofs;
1380 
1381  // Determine number of nodes in base state problem
1382  const unsigned n_node = Base_state_problem_pt->mesh_pt()->nnode();
1383 
1384  // Set up storage for norm of vector of differences between current
1385  // dofs and dofs from one period ago (initialise to 1.0)
1386  double norm_over_ndof = 1.0;
1387 
1388  // Set up storage for loop counter
1389  unsigned loop_counter = 0;
1390 
1391  // While this norm is greater than the base state convergence criterion
1392  while(norm_over_ndof > tolerance)
1393  {
1394  // Record the dofs corresponding to the azimuthal velocity at the end
1395  // of the previous period (equivalent to recording them at the start
1396  // of the current period, as we do here)
1397  for(unsigned n=0;n<n_node;n++)
1398  {
1399  previous_dofs.push_back
1400  (Base_state_problem_pt->mesh_pt()->node_pt(n)->value(2));
1401  }
1402 
1403  // Timestepping loop
1404  for(unsigned i=0;i<n_timesteps;i++)
1405  {
1406  // Output timestep and global time
1407  std::cout << " Timestep " << i+1 << " of " << n_timesteps
1408  << "; Time is " << Base_state_problem_pt->time_pt()->time()
1409  << " (base)" << std::endl;
1410 
1411  // Take timestep
1412  Base_state_problem_pt->unsteady_newton_solve(dt);
1413 
1414  // Output solution
1415  Base_state_problem_pt->doc_solution(doc_info,false);
1416 
1417  // Increment counter for solutions
1418  doc_info.number()++;
1419  }
1420 
1421  // Get dofs from base state problem at the end of this period
1422  for(unsigned n=0;n<n_node;n++)
1423  {
1424  current_dofs.push_back
1425  (Base_state_problem_pt->mesh_pt()->node_pt(n)->value(2));
1426  }
1427 
1428  // Find difference between each current dof and the corresponding
1429  // previous dof from one period ago
1430  for(unsigned i=0;i<n_node;i++) { current_dofs[i] -= previous_dofs[i]; }
1431 
1432  // Evaluate L2 norm of this vector over number of dofs
1433  norm_over_ndof = 0.0;
1434  for(unsigned i=0;i<n_node;i++)
1435  {
1436  norm_over_ndof += (current_dofs[i]*current_dofs[i]);
1437  }
1438  norm_over_ndof = sqrt(norm_over_ndof);
1439  norm_over_ndof = norm_over_ndof/n_node;
1440 
1441  ios_base::fmtflags flags = cout.flags(); // Save old flags
1442  cout << " ==> norm_over_ndof = " << scientific << norm_over_ndof
1443  << ", tolerance = " << tolerance << endl;
1444  cout.flags(flags); // Set the flags to the way they were
1445 
1446  // Document norm and tolerance in trace file
1447  Trace_file_periodic_base_flow << Base_state_problem_pt->time_pt()->time()
1448  << " " << norm_over_ndof
1449  << " " << tolerance << std::endl;
1450 
1451  // Clear the Vector<double>s
1452  previous_dofs.clear();
1453  current_dofs.clear();
1454 
1455  // Increment loop counter
1456  loop_counter++;
1457 
1458  } // End of while loop
1459 
1460  cout << "\nI have broken out of the base state loop after "
1461  << loop_counter << " periods.\n" << endl;
1462 
1463  // (Re-)set the perturbed state problem's global time to correspond
1464  // to the base state problem's global time
1465  this->set_perturbed_state_time_equal_to_base_state_time();
1466 
1467  // Return the number of periods required to reach a periodic flow
1468  return loop_counter;
1469 
1470 } // End of set_up_periodic_base_flow
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
StabilityProblem::set_perturbed_state_time_equal_to_base_state_time
void set_perturbed_state_time_equal_to_base_state_time()
Definition: linearised_axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:1013

References i, n, oomph::DocInfo::number(), and sqrt().

◆ unsteady_run() [1/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::unsteady_run ( const double &  dt,
const unsigned &  n_timesteps,
DocInfo &  doc_info 
)

Integrate forwards in time with timestep dt for n_timesteps.

784 {
785  // Timestepping loop
786  for(unsigned i=0;i<n_timesteps;i++)
787  {
788  // Output timestep and global time
789  std::cout << " Timestep " << i+1 << " of " << n_timesteps
790  << "; Time is " << Perturbed_state_problem_pt->time_pt()->time()
791  << std::endl;
792 
793  // Take timestep
794  Perturbed_state_problem_pt->unsteady_newton_solve(dt);
795 
796  } // End of loop over timesteps
797 } // End of unsteady_run

References i.

Referenced by StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::unsteady_run_for_power_method().

◆ unsteady_run() [2/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::unsteady_run ( const double &  dt,
const unsigned &  n_timesteps,
DocInfo &  doc_info 
)

Integrate forwards in time with timestep dt for n_timesteps.

◆ unsteady_run() [3/3]

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::unsteady_run ( const double &  dt,
const unsigned &  n_timesteps,
DocInfo &  doc_info,
const bool &  output_base_soln,
Vector< DoubleVector > &  perturbed_state_dofs_over_period,
Vector< DoubleVector > &  base_state_dofs_over_period 
)

Integrate forwards in time with timestep dt for n_timesteps.

Integrate forwards in time with timestep dt for n_timesteps Record an entire period's worth of the perturbed state problem solution in perturbed_state_dofs_over_period

1106 {
1107  // "Outer" timestepping loop
1108  for(unsigned i=0;i<n_timesteps;i++)
1109  {
1110  // Output timestep and global time
1111  std::cout << " Timestep " << i+1 << " of " << n_timesteps
1112  << "; Time is " << Base_state_problem_pt->time_pt()->time()
1113  << " (base), " << Perturbed_state_problem_pt->time_pt()->time()
1114  << " (perturbed)" << std::endl;
1115 
1116  // On the first "timestep", perform n_timesteps BDF1 timesteps
1117  if(i==0)
1118  {
1119  // Set the perturbed state problem's timestepper to work in BDF1 mode
1120  Base_state_problem_pt->time_stepper_pt()->turn_on_bdf1_mode();
1121  Perturbed_state_problem_pt->time_stepper_pt()->turn_on_bdf1_mode();
1122 
1123  // Determine the BDF1 timestep
1124  const double dt_bdf1 = dt/n_timesteps;
1125 
1126  // Initialise timestep (also sets timestepper weights)
1127  this->initialise_dt(dt_bdf1);
1128 
1129  // BDF1 timestepping loop
1130  for(unsigned j=0;j<n_timesteps;j++)
1131  {
1132 
1133  cout << " - BDF1 loop: Timestep " << j+1 << " of " << n_timesteps
1134  << "; Time (before advance) is "
1135  << Base_state_problem_pt->time_pt()->time()
1136  << " (base), "
1137  << Perturbed_state_problem_pt->time_pt()->time()
1138  << " (perturbed)" << std::endl;
1139 
1140  // Take timestep
1141  if(First_power_method_iteration)
1142  {
1143  Base_state_problem_pt->unsteady_newton_solve(dt_bdf1);
1144  Base_state_problem_pt->get_dofs(base_state_dofs_over_period[j]);
1145  }
1146  else
1147  {
1148  Base_state_problem_pt
1149  ->assign_eigenvector_to_dofs(base_state_dofs_over_period[j]);
1150  }
1151  Perturbed_state_problem_pt->unsteady_newton_solve(dt_bdf1);
1152 
1153  // Document in trace file on all steps other than the last
1154  // (when documentation is performed in outer loop)
1155  if(j!=(n_timesteps-1))
1156  {
1157  this->doc_solution(doc_info,false,false);
1158  }
1159 
1160  } // End of BDF1 timestepping loop
1161 
1162  // Now set the timestepper back to BDF2 mode
1163  Base_state_problem_pt->time_stepper_pt()->turn_off_bdf1_mode();
1164  Perturbed_state_problem_pt->time_stepper_pt()->turn_off_bdf1_mode();
1165  }
1166 
1167  // On the second "timestep", use a timestep that's twice as big as
1168  // used during the BDF1 timestepping loop
1169  else if(i==1)
1170  {
1171  // Determine the timestep and number of timesteps for this loop
1172  const double dt_thisloop = 2*dt/n_timesteps;
1173  const double n_timesteps_thisloop = 0.5*n_timesteps;
1174 
1175  // Determine the timestep for the previous loop
1176  const double dt_prevloop = dt/n_timesteps;
1177 
1178  Vector<double> dt_vector;
1179  dt_vector.push_back(dt_prevloop);
1180  dt_vector.push_back(dt_thisloop);
1181 
1182  // Initialise timestep for first time around this timestepping loop
1183  // (also sets timestepper weights)
1184  this->initialise_dt(dt_vector);
1185 
1186  // Timestepping loop
1187  for(unsigned j=0;j<n_timesteps_thisloop;j++)
1188  {
1189  std::cout << " - nstep/2 loop: Timestep " << j+1
1190  << " of " << n_timesteps_thisloop
1191  << "; Time (before advance) is "
1192  << Base_state_problem_pt->time_pt()->time()
1193  << " (base), "
1194  << Perturbed_state_problem_pt->time_pt()->time()
1195  << " (perturbed)" << std::endl;
1196 
1197  // Take timestep
1198  if(First_power_method_iteration)
1199  {
1200  Base_state_problem_pt->unsteady_newton_solve(dt_thisloop);
1201  Base_state_problem_pt
1202  ->get_dofs(base_state_dofs_over_period[(n_timesteps+j)]);
1203  }
1204  else
1205  {
1206  Base_state_problem_pt->assign_eigenvector_to_dofs
1207  (base_state_dofs_over_period[(n_timesteps+j)]);
1208  }
1209  Perturbed_state_problem_pt->unsteady_newton_solve(dt_thisloop);
1210 
1211  // Document in trace file on all steps other than the last
1212  // (when documentation is performed in outer loop)
1213  if(j!=(n_timesteps_thisloop-1))
1214  {
1215  this->doc_solution(doc_info,false,false);
1216  }
1217 
1218  // On the first time through this loop, initialise the timestep
1219  // again with equal timesteps
1220  if(j==0) { this->initialise_dt(dt_thisloop); }
1221  }
1222  }
1223 
1224  // On the third "timestep", use a timestep that's five times as big as
1225  // used during the BDF1 timestepping loop
1226  else if(i==2)
1227  {
1228  // Determine the timestep and number of timesteps for this loop
1229  const double dt_thisloop = 5*dt/n_timesteps;
1230  const double n_timesteps_thisloop = 0.2*n_timesteps;
1231 
1232  // Determine the timestep for the previous loop
1233  const double dt_prevloop = 2*dt/n_timesteps;
1234 
1235  Vector<double> dt_vector;
1236  dt_vector.push_back(dt_prevloop);
1237  dt_vector.push_back(dt_thisloop);
1238 
1239  // Initialise timestep for first time around this timestepping loop
1240  // (also sets timestepper weights)
1241  this->initialise_dt(dt_vector);
1242 
1243  // Timestepping loop
1244  for(unsigned j=0;j<n_timesteps_thisloop;j++)
1245  {
1246  std::cout << " - nstep/2 loop: Timestep " << j+1
1247  << " of " << n_timesteps_thisloop
1248  << "; Time (before advance) is "
1249  << Base_state_problem_pt->time_pt()->time()
1250  << " (base), "
1251  << Perturbed_state_problem_pt->time_pt()->time()
1252  << " (perturbed)" << std::endl;
1253 
1254  // Take timestep
1255  if(First_power_method_iteration)
1256  {
1257  Base_state_problem_pt->unsteady_newton_solve(dt_thisloop);
1258  Base_state_problem_pt
1259  ->get_dofs(base_state_dofs_over_period[((n_timesteps*1.5)+j)]);
1260  }
1261  else
1262  {
1263  Base_state_problem_pt->assign_eigenvector_to_dofs
1264  (base_state_dofs_over_period[((n_timesteps*1.5)+j)]);
1265  }
1266  Perturbed_state_problem_pt->unsteady_newton_solve(dt_thisloop);
1267 
1268  // Document in trace file on all steps other than the last
1269  // (when documentation is performed in outer loop)
1270  if(j!=(n_timesteps_thisloop-1))
1271  {
1272  this->doc_solution(doc_info,false,false);
1273  }
1274 
1275  // On the first time through this loop, initialise the timestep
1276  // again with equal timesteps
1277  if(j==0) { this->initialise_dt(dt_thisloop); }
1278  }
1279  }
1280 
1281  // On the fourth "timestep", use the actual timestep dt but need to
1282  // specify past dts as changing step size from last time
1283  else if(i==3)
1284  {
1285  // Determine the timestep for the previous loop
1286  const double dt_prevloop = 5*dt/n_timesteps;
1287 
1288  Vector<double> dt_vector;
1289  dt_vector.push_back(dt_prevloop);
1290  dt_vector.push_back(dt);
1291 
1292  // Initialise timestep for first time around this timestepping loop
1293  // (also sets timestepper weights)
1294  this->initialise_dt(dt_vector);
1295 
1296  // Take timestep
1297  if(First_power_method_iteration)
1298  {
1299  Base_state_problem_pt->unsteady_newton_solve(dt);
1300  Base_state_problem_pt
1301  ->get_dofs(base_state_dofs_over_period
1302  [((n_timesteps*1.5)+(n_timesteps*0.2))]);
1303  }
1304  else
1305  {
1306  Base_state_problem_pt->assign_eigenvector_to_dofs
1307  (base_state_dofs_over_period[((n_timesteps*1.5)+(n_timesteps*0.2))]);
1308  }
1309  Perturbed_state_problem_pt->unsteady_newton_solve(dt);
1310 
1311  // Initialise timestep back to equal timesteps of length dt
1312  // (also sets timestepper weights)
1313  this->initialise_dt(dt);
1314  }
1315 
1316  // On all subsequent timesteps perform normal BDF2 timestepping
1317  else
1318  {
1319  // Take timestep
1320  if(First_power_method_iteration)
1321  {
1322  Base_state_problem_pt->unsteady_newton_solve(dt);
1323  Base_state_problem_pt
1324  ->get_dofs(base_state_dofs_over_period
1325  [((n_timesteps*1.5)+(n_timesteps*0.2)+(i-3))]);
1326  }
1327  else
1328  {
1329  Base_state_problem_pt->assign_eigenvector_to_dofs
1330  (base_state_dofs_over_period
1331  [((n_timesteps*1.5)+(n_timesteps*0.2)+(i-3))]);
1332  }
1333  Perturbed_state_problem_pt->unsteady_newton_solve(dt);
1334  }
1335 
1336  // At end of each "outer" timestep, document info in trace file
1337  // and also output base state solution (if flag is set appropriately)
1338  this->doc_solution(doc_info,output_base_soln,false);
1339 
1340  // Store perturbed state problems dofs so that eigenvector over one
1341  // period can be outputted at the end of the power method
1342  Perturbed_state_problem_pt->get_dofs(perturbed_state_dofs_over_period[i]);
1343 
1344  // Increment counter for solutions
1345  doc_info.number()++;
1346 
1347  } // End of loop over "outer timesteps"
1348 } // End of unsteady run
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References i, j, and oomph::DocInfo::number().

◆ unsteady_run_for_power_method()

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
void StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::unsteady_run_for_power_method ( const double &  dt,
const unsigned &  n_timesteps,
DocInfo &  doc_info,
const DoubleVector &  input_perturbed_state_dofs,
DoubleVector &  output_perturbed_state_dofs,
const bool &  output_base_soln,
Vector< DoubleVector > &  perturbed_state_dofs_over_period,
Vector< DoubleVector > &  base_state_dofs_over_period 
)
inline

Integrate forwards in time with timestep dt for n_timesteps. Re-set the initial conditions for the perturbed state problem to be input_perturbed_state_dofs and, at the end of the timestepping loop, return the perturbed state problem's dofs

991  {
992  // Re-set the starting dof vector for the perturbed state problem
993  Perturbed_state_problem_pt->set_dofs(input_perturbed_state_dofs);
994 
995  // Perform standard unsteady run
996  this->unsteady_run(dt,n_timesteps,doc_info,
997  output_base_soln,perturbed_state_dofs_over_period,
998  base_state_dofs_over_period);
999 
1000  // Get dofs from perturbed state problem
1001  Perturbed_state_problem_pt->get_dofs(output_perturbed_state_dofs);
1002  }

References StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt, and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::unsteady_run().

Member Data Documentation

◆ Base_state_problem_pt

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
BaseStateProblem< BASE_ELEMENT > * StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Base_state_problem_pt
private

Pointer to base state problem class.

Referenced by StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::base_flow_steady_newton_solve(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::base_state_problem_pt(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::close_trace_files(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::create_trace_files(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::doc_solution(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::initialise_dt(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::initialise_trace_files(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::pass_updated_nondim_parameters_to_elements(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::reset_global_time_to_zero(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_initial_condition(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_perturbed_state_time_equal_to_base_state_time(), and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::StabilityProblem().

◆ First_power_method_iteration

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
bool StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::First_power_method_iteration
private

Flag to indicate whether or not we are doing the first iteration of the power method

Referenced by StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_first_power_method_iteration_flag_to_false(), and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_first_power_method_iteration_flag_to_true().

◆ Perturbed_state_problem_pt

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
PerturbedStateProblem< BASE_ELEMENT, PERTURBED_ELEMENT > * StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Perturbed_state_problem_pt
private

Pointer to perturbed state problem class.

Referenced by StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::close_trace_files(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::create_trace_files(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::doc_solution(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::get_perturbed_state_problem_dofs(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::initialise_dt(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::initialise_trace_files(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::pass_updated_nondim_parameters_to_elements(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::perturbed_state_problem_pt(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::reset_global_time_to_zero(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_initial_condition(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::set_perturbed_state_time_equal_to_base_state_time(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::StabilityProblem(), and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::unsteady_run_for_power_method().

◆ Trace_file_periodic_base_flow

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
ofstream StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Trace_file_periodic_base_flow
private

Trace file which documents the convergence of the base flow to a periodic state

Referenced by StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::close_trace_files(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::create_trace_files(), and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::initialise_trace_files().

◆ Trace_file_power_method

template<class BASE_ELEMENT , class PERTURBED_ELEMENT >
ofstream StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::Trace_file_power_method
private

Trace file which documents the convergence of the power method.

Referenced by StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::close_trace_files(), StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::create_trace_files(), and StabilityProblem< BASE_ELEMENT, PERTURBED_ELEMENT >::initialise_trace_files().

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