AirwayReopeningProblem< ELEMENT > Class Template Reference

Bretherton problem. More...

+ Inheritance diagram for AirwayReopeningProblem< ELEMENT >:

Public Member Functions

 AirwayReopeningProblem ()
 Constructor: More...
 
 ~AirwayReopeningProblem ()
 Destructor, clean up all allocated memory. More...
 
void actions_after_change_in_global_parameter (double *const &parameter_pt)
 
void actions_before_newton_convergence_check ()
 
void actions_before_newton_solve ()
 Update before solve: empty. More...
 
void actions_after_newton_solve ()
 
void fix_pressure (const unsigned &e, const unsigned &l, const double &pvalue)
 Fix pressure value l in element e to value p_value. More...
 
void parameter_study (const unsigned &nsteps, const bool &restart)
 Run a parameter study; perform specified number of steps. More...
 
void doc_solution (DocInfo &doc_info)
 Doc the solution. More...
 
void construct_poisson_outlet_elements ()
 Set up the actual outlet elements. More...
 
void delete_outlet_elements ()
 Delete the outlet elements. More...
 
void construct_pressure_gradient_outlet_elements ()
 Set up the inital outlet elements. More...
 
double get_outlet_flux ()
 
void change_boundary_conditions ()
 
void dump (ofstream &dump_file) const
 Dump the entire problem. More...
 
void read (ifstream &restart_file)
 Read the entire problem. More...
 
void connect_walls (Mesh *const &wall_mesh_pt, GeomObject *const &geom_object_pt, map< double, pair< GeomObject *, Vector< double > > > &globalmap)
 Setup the coupling between the upper and lower walls. 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...
 
OomphCommunicatorcommunicator_pt ()
 access function to the oomph-lib communicator More...
 
const OomphCommunicatorcommunicator_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 ()
 
LinearSolvermass_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...
 
Timetime_pt () const
 Return a pointer to the global time object (const version). More...
 
doubletime ()
 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 TimeSteppertime_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...
 
doubleminimum_dt ()
 Access function to min timestep in adaptive timestepping. More...
 
doublemaximum_dt ()
 Access function to max timestep in adaptive timestepping. More...
 
unsignedmax_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...
 
doublemax_residuals ()
 
booltime_adaptive_newton_crash_on_solve_fail ()
 Access function for Time_adaptive_newton_crash_on_solve_fail. More...
 
doublenewton_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...
 
doubleglobal_dof_pt (const unsigned &i)
 
doubledof (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...
 
doubledof_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...
 
booluse_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 Problemmake_copy ()
 
virtual void read (std::ifstream &restart_file, bool &unsteady_restart)
 
void dump (const std::string &dump_file_name) const
 
void delete_all_external_storage ()
 
virtual void symmetrise_eigenfunction_for_adaptive_pitchfork_tracking ()
 
doublebifurcation_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 ()
 
intsign_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 Attributes

ELEMENT * Control_element_pt
 Pointer to control element. More...
 
ofstream Trace_file
 Trace file. More...
 
BrethertonSpineMesh< ELEMENT, SpineLineFluidInterfaceElement< ELEMENT > > * Bulk_mesh_pt
 Pointer to bulk mesh. More...
 
OneDLagrangianMesh< FSIHermiteBeamElement > * Upper_wall_mesh_pt
 Pointer to Wall mesh. More...
 
OneDLagrangianMesh< FSIHermiteBeamElement > * Lower_wall_mesh_pt
 Pointer to lower wall mesh. More...
 
MeshConstraint_mesh_pt
 Pointer to the constraint mesh. More...
 
GeomObjectUndeformed_upper_wall_geom_pt
 Pointer to the Undeformed wall. More...
 
GeomObjectUndeformed_lower_wall_geom_pt
 Pointer to the Undeformed wall. More...
 
double Fixed_spine_height
 Fixed height value. More...
 
DataBubble_pressure_data_pt
 Data value that represents the bubble pressure. More...
 
DataMesh_fraction_at_transition_pt
 
Vector< FiniteElement * > Gravity_traction_element_pt
 Vector of pointers to the traction elements. More...
 
Vector< FiniteElement * > Outlet_traction_element_pt
 Vector of pointers to the outlet traciton elements. More...
 
FluxConstraintFlux_constraint_pt
 Flux constraint element. More...
 
bool Frontal_solver
 Boolean flag to specify the use of the frontal solver. More...
 

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_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_bifurcation_parameter ()
 
virtual void actions_after_parameter_increase (double *const &parameter_pt)
 
doubledof_derivative (const unsigned &i)
 
doubledof_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 *, boolCalculate_dparameter_analytic
 
bool Calculate_hessian_products_analytic
 
LinearAlgebraDistributionDof_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< doubleMax_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< doubleDof_derivative
 Storage for the derivative of the problem variables wrt arc-length. More...
 
Vector< doubleDof_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
 
OomphCommunicatorCommunicator_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 AirwayReopeningProblem< ELEMENT >

Bretherton problem.

Constructor & Destructor Documentation

◆ AirwayReopeningProblem()

template<class ELEMENT >
AirwayReopeningProblem< ELEMENT >::AirwayReopeningProblem

Constructor:

Problem constructor.

1606 {
1607  Frontal_solver = false;
1608 
1609  // Number of elements in the deposited film region
1610  unsigned nx1=100;//300;
1611 
1612  // Number of elements in the bottom part of the transition region
1613  unsigned nx2=6;//12;
1614 
1615  // Number of elements in channel region
1616  unsigned nx3=100;//300;
1617 
1618  // Number of elements in the vertical part of the transition region
1619  // (=half the number of elements in the liquid filled region ahead
1620  // of the finger tip)
1621  unsigned nhalf=2;//4;
1622 
1623  // Number of elements through thickness of deposited film
1624  unsigned nh=2;//3;
1625 
1626  // Number of elements in each wall
1627  unsigned nwall = 100;// 450;
1628 
1629  // Thickness of deposited film
1630  double h = 0.128;
1631 
1632  // Start coordinate on wall
1633  double xi0= -300.0;
1634 
1635  double start_transition = -1.0;
1636 
1637  // End of transition region on wall
1638  double xi1 = 1.5;
1639 
1640  // End of liquid filled region (inflow) on wall
1641  double xi2 = 150.0;//300.0
1642 
1643  //Set the frontal solver, if desired
1644  if(Frontal_solver)
1645  {
1646  linear_solver_pt() = new HSL_MA42;
1647  }
1648 
1649  //Create a bubble presure data
1650  Bubble_pressure_data_pt = new Data(1);
1651  //This will be global data
1653 
1654  //Create a mesh fraction data
1657  //Pin the value
1660 
1661  // The underformed wall is a straight line at y = 1.0
1663 
1664  // Create the Lagrangian Mesh
1666  new OneDLagrangianMesh<FSIHermiteBeamElement>(nwall,xi0,xi2,
1668 
1669  //Loop over the wall elements
1670  unsigned n_upper_wall_element = Upper_wall_mesh_pt->nelement();
1671  for(unsigned e=0;e<n_upper_wall_element;e++)
1672  {
1673  FSIHermiteBeamElement *element_pt = dynamic_cast<FSIHermiteBeamElement*>(
1675 
1676  //Assign the undeformed beam shape
1678 
1679  //Set the load vector
1681 
1682  //Set the wall thickness
1683  element_pt->h_pt() = &Global_Physical_Variables::H;
1684 
1685  //Set the axial pre-stress
1686  element_pt->sigma0_pt() = &Global_Physical_Variables::T;
1687 
1688  //Function that specifies the load ratios
1689  element_pt->q_pt() = &Global_Physical_Variables::Q;
1690 
1691  //In this case the normal is directed out of the fluid
1692  element_pt->set_normal_pointing_out_of_fluid();
1693  }
1694 
1695  //Create a geometric object that represents the wall geometry
1696  MeshAsGeomObject* upper_wall_element_pt =
1698 
1699  // The underformed lower wall is a straight line at y = -1.0
1701 
1702  // Create the Lagrangian Mesh
1704  new OneDLagrangianMesh<FSIHermiteBeamElement>(nwall,xi0,xi2,
1706 
1707  //Loop over the wall elements
1708  unsigned n_lower_wall_element = Lower_wall_mesh_pt->nelement();
1709  for(unsigned e=0;e<n_lower_wall_element;e++)
1710  {
1711  FSIHermiteBeamElement *element_pt = dynamic_cast<FSIHermiteBeamElement*>(
1713 
1714  //Assign the undeformed beam shape
1716 
1717  //Set the load vector
1718  element_pt->load_vector_fct_pt() =
1720 
1721  //Set the wall thickness
1722  element_pt->h_pt() = &Global_Physical_Variables::H;
1723 
1724  //Set the axial pre-stress
1725  element_pt->sigma0_pt() = &Global_Physical_Variables::T;
1726 
1727  //Function that specifies the load ratios
1728  element_pt->q_pt() = &Global_Physical_Variables::Q;
1729 
1730  //In this case the normal is OK
1731  //element_pt->set_normal_pointing_into_fluid();
1732  }
1733 
1734  //Create a geometric object that represents the wall geometry
1735  MeshAsGeomObject* lower_wall_element_pt =
1737 
1738  // Create wall geom objects
1739  //GeomObject* lower_wall_pt = new StraightLine(-1.0);
1740  GeomObject* lower_wall_pt = lower_wall_element_pt;
1742 
1743 
1744  //GeomObject* upper_wall_pt=new StraightLine( 1.0);
1745  GeomObject* upper_wall_pt = upper_wall_element_pt;
1747 
1748  //Now create the mesh
1749  Bulk_mesh_pt = new BrethertonSpineMesh<ELEMENT,
1751  (nx1,nx2,nx3,nh,nhalf,h,lower_wall_pt,upper_wall_pt,xi0,start_transition,
1752  xi1,xi2);
1753 
1754  //Set the vertical fraction
1757 
1758  // Store the control element
1760 
1761  //Create a fixed element using the central spine
1762  FixSpineHeightElement* fix_spine_element_pt =
1764  static_cast<SpineNode*>(Control_element_pt->node_pt(8)));
1765 
1766  //Set the fixed spine height
1768  static_cast<SpineNode*>(Control_element_pt
1769  ->node_pt(8))->spine_pt()->height();
1770 
1771  //Add the fixed height to it
1772  fix_spine_element_pt->height_pt() = &Fixed_spine_height;
1773 
1774  //Set the pressure data
1775  fix_spine_element_pt->set_traded_pressure_data(Bubble_pressure_data_pt);
1776 
1777  //Add the Fixed element to the mesh
1778  Bulk_mesh_pt->add_element_pt(fix_spine_element_pt);
1779 
1780  //Loop over the elements next to the inflow boundarys
1781  for(unsigned b=3;b<=5;b+=2)
1782  {
1783  unsigned nbound_element = Bulk_mesh_pt->nboundary_element(b);
1784  for(unsigned e=0;e<nbound_element;e++)
1785  {
1790 
1791  //Add the elements to the local storage
1792  Gravity_traction_element_pt.push_back(inlet_element);
1793 
1794  //Add the elements to the mesh
1795  Bulk_mesh_pt->add_element_pt(inlet_element);
1796  }
1797  }
1798 
1799  //Set a prescribed flux element
1801 
1802  Constraint_mesh_pt = new Mesh;
1803 
1805 
1807 
1809 
1811 
1812  {
1813  //Set the boundary coordinates on the upper wall
1814  unsigned nbound = Bulk_mesh_pt->nboundary_node(2);
1815  Vector<double> zeta(1);
1816  for(unsigned n=0;n<nbound;n++)
1817  {
1818  Node* node_pt = Bulk_mesh_pt->boundary_node_pt(2,n);
1819  zeta[0] = node_pt->x(0);
1820  node_pt->set_coordinates_on_boundary(2,zeta);
1821  }
1822  }
1823 
1824  //Sort out the interaction the upper wall
1825  FSI_functions::setup_fluid_load_info_for_solid_elements<ELEMENT,2>
1827 
1828  {
1829  //Set the boundary coordinates on the lower wall
1830  unsigned nbound = Bulk_mesh_pt->nboundary_node(0);
1831  Vector<double> zeta(1);
1832  for(unsigned n=0;n<nbound;n++)
1833  {
1834  Node* node_pt = Bulk_mesh_pt->boundary_node_pt(0,n);
1835  zeta[0] = node_pt->x(0);
1836  node_pt->set_coordinates_on_boundary(0,zeta);
1837  }
1838  }
1839 
1840  //Sort out the interaction the lower wall
1841  FSI_functions::setup_fluid_load_info_for_solid_elements<ELEMENT,2>
1843 
1845 
1847 
1849 
1851 
1853 
1854 
1855  // Connect the upper mesh
1856  connect_walls(Lower_wall_mesh_pt,upper_wall_pt,
1858 
1859  connect_walls(Upper_wall_mesh_pt,lower_wall_pt,
1861 
1862  // Set the boundary conditions for this problem: All nodes are
1863  // free by default -- just pin the ones that have Dirichlet conditions
1864  // here
1865 
1866  //Pin one end of the mesh (prevent transvere slide)
1869 
1870  //Pin both ends of the lower mesh
1873 
1874  //Pin the lower mesh vertically
1876 
1877  // No slip on boundaries 0 1 (inflow) and 2
1878  for(unsigned long ibound=0;ibound<=2;ibound++)
1879  {
1880  unsigned num_nod=Bulk_mesh_pt->nboundary_node(ibound);
1881  for (unsigned inod=0;inod<num_nod;inod++)
1882  {
1883  //Leave the inflow traction free in the both directions
1884  if(ibound!=1)
1885  {
1886  Bulk_mesh_pt->boundary_node_pt(ibound,inod)->pin(0);
1887  Bulk_mesh_pt->boundary_node_pt(ibound,inod)->pin(1);
1888  }
1889  }
1890  }
1891 
1892  // Set the values of the velocities on the wall boundaries
1893  for (unsigned ibound=0;ibound<=2;ibound+=2)
1894  {
1895  unsigned num_nod=Bulk_mesh_pt->nboundary_node(ibound);
1896  for (unsigned inod=0;inod<num_nod;inod++)
1897  {
1898  // Parallel flow
1899  Bulk_mesh_pt->boundary_node_pt(ibound,inod)->set_value(0,-1.0);
1900  Bulk_mesh_pt->boundary_node_pt(ibound,inod)->set_value(1, 0.0);
1901  }
1902  }
1903 
1904  //Loop over the upper wall
1905  {
1906  unsigned b=2;
1907  unsigned num_nod=Bulk_mesh_pt->nboundary_node(b);
1908  for (unsigned inod=0;inod<num_nod;inod++)
1909  {
1910  //cast to a spine node
1911  SpineNode* spine_node_pt = static_cast<SpineNode*>
1912  (Bulk_mesh_pt->boundary_node_pt(b,inod));
1913 
1914  //Find the id of the node
1915  unsigned id = spine_node_pt->node_update_fct_id();
1916  switch(id)
1917  {
1918 
1919  case 4:
1920  case 5:
1921  //Add the auxilliary function, so that the no-slip condition is
1922  //applied properly
1923  spine_node_pt->set_auxiliary_node_update_fct_pt(
1925  break;
1926 
1927  case 6:
1928  spine_node_pt->set_auxiliary_node_update_fct_pt(
1930  break;
1931 
1932  default:
1933  cout << "Should not get here" << std::endl;
1934  break;
1935  }
1936  }
1937  }
1938 
1939  //Loop over the lower wall
1940  {
1941  unsigned b=0;
1942  unsigned num_nod=Bulk_mesh_pt->nboundary_node(b);
1943  for (unsigned inod=0;inod<num_nod;inod++)
1944  {
1945  //cast to a spine node
1946  SpineNode* spine_node_pt = static_cast<SpineNode*>
1947  (Bulk_mesh_pt->boundary_node_pt(b,inod));
1948 
1949  //Find the id of the node
1950  unsigned id = spine_node_pt->node_update_fct_id();
1951  switch(id)
1952  {
1953 
1954  case 0:
1955  case 1:
1956  case 6:
1957  //Add the auxilliary function, so that the no-slip condition is
1958  //applied properly
1959  spine_node_pt->set_auxiliary_node_update_fct_pt(
1961  break;
1962 
1963  default:
1964  cout << "Should not get here" << std::endl;
1965  break;
1966  }
1967  }
1968  }
1969 
1970  //Parallel, uniform outflow on boundaries 3 and 5
1971  for (unsigned ibound=3;ibound<=5;ibound+=2)
1972  {
1973  unsigned num_nod=Bulk_mesh_pt->nboundary_node(ibound);
1974  for (unsigned inod=0;inod<num_nod;inod++)
1975  {
1976  // Parallel inflow
1977  Bulk_mesh_pt->boundary_node_pt(ibound,inod)->set_value(0, -1.0);
1978  Bulk_mesh_pt->boundary_node_pt(ibound,inod)->set_value(1, 0.0);
1979  }
1980  }
1981 
1982  //Loop over the elements in the layer
1983  unsigned long n_bulk=Bulk_mesh_pt->nbulk();
1984  for(unsigned long i=0;i<n_bulk;i++)
1985  {
1986  //Cast to a fluid element
1987  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->
1988  bulk_element_pt(i));
1989  //Set the Reynolds number, etc
1990  el_pt->re_pt() = &Global_Physical_Variables::Re;
1991  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
1992  el_pt->re_invfr_pt() = &Global_Physical_Variables::ReInvFr;
1993  el_pt->g_pt() = &Global_Physical_Variables::G;
1994  }
1995 
1996  //Loop over 1D interface elements and set capillary number
1997  unsigned interface_element_pt_range = Bulk_mesh_pt->ninterface_element();
1998  for(unsigned i=0;i<interface_element_pt_range;i++)
1999  {
2000  //Cast to a interface element
2004 
2005  //Set the Capillary number
2007  //Set the external pressure data
2009 
2010  //We need to make sure that we hijack the nodes on the boundaries
2011  if(i==0) {el_pt->hijack_nodal_value(0,0);}
2012  if(i==interface_element_pt_range-1)
2013  {el_pt->hijack_nodal_value(el_pt->nnode()-1,0);}
2014  }
2015 
2016  //Loop over the traction elements
2017  unsigned n_traction = Gravity_traction_element_pt.size();
2018  for(unsigned e=0;e<n_traction;e++)
2019  {
2020  //Cast to the traction element
2024 
2025  //Set the Reynolds number divided by the Froude number
2027  el_pt->g_pt() = &Global_Physical_Variables::G;
2028  }
2029 
2030  //Do equation numbering
2031  cout << "Number of unknowns: " << assign_eqn_numbers() << std::endl;
2032 
2033  //Sort the elements if using the frontal solver
2034  //If this is not done, it's very slow indeed
2035  //If it is done, then it's very fast compared to SuperLU
2036  if(Frontal_solver)
2037  {
2038  std::sort(mesh_pt()->element_pt().begin(),
2039  mesh_pt()->element_pt().end(),
2040  ElementCmp());
2041  }
2042 }
int i
Definition: BiCGSTAB_step_by_step.cpp:9
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
Array< double, 1, 3 > e(1./3., 0.5, 2.)
Scalar * b
Definition: benchVecAdd.cpp:17
double Fixed_spine_height
Fixed height value.
Definition: elastic_bretherton.cc:1577
Mesh * Constraint_mesh_pt
Pointer to the constraint mesh.
Definition: elastic_bretherton.cc:1568
void connect_walls(Mesh *const &wall_mesh_pt, GeomObject *const &geom_object_pt, map< double, pair< GeomObject *, Vector< double > > > &globalmap)
Setup the coupling between the upper and lower walls.
Definition: elastic_bretherton.cc:1486
Data * Bubble_pressure_data_pt
Data value that represents the bubble pressure.
Definition: elastic_bretherton.cc:1580
ELEMENT * Control_element_pt
Pointer to control element.
Definition: elastic_bretherton.cc:1551
OneDLagrangianMesh< FSIHermiteBeamElement > * Lower_wall_mesh_pt
Pointer to lower wall mesh.
Definition: elastic_bretherton.cc:1565
BrethertonSpineMesh< ELEMENT, SpineLineFluidInterfaceElement< ELEMENT > > * Bulk_mesh_pt
Pointer to bulk mesh.
Definition: elastic_bretherton.cc:1559
Vector< FiniteElement * > Gravity_traction_element_pt
Vector of pointers to the traction elements.
Definition: elastic_bretherton.cc:1588
Data * Mesh_fraction_at_transition_pt
Definition: elastic_bretherton.cc:1585
OneDLagrangianMesh< FSIHermiteBeamElement > * Upper_wall_mesh_pt
Pointer to Wall mesh.
Definition: elastic_bretherton.cc:1562
bool Frontal_solver
Boolean flag to specify the use of the frontal solver.
Definition: elastic_bretherton.cc:1597
FluxConstraint * Flux_constraint_pt
Flux constraint element.
Definition: elastic_bretherton.cc:1594
GeomObject * Undeformed_lower_wall_geom_pt
Pointer to the Undeformed wall.
Definition: elastic_bretherton.cc:1574
void construct_pressure_gradient_outlet_elements()
Set up the inital outlet elements.
Definition: elastic_bretherton.cc:1345
GeomObject * Undeformed_upper_wall_geom_pt
Pointer to the Undeformed wall.
Definition: elastic_bretherton.cc:1571
Function-type-object to perform comparison of elements.
Definition: elastic_bretherton.cc:317
Point element that is used to set the bubble pressure.
Definition: elastic_bretherton.cc:994
double *& height_pt()
Access function to the prescribed spine height.
Definition: elastic_bretherton.cc:1040
void set_traded_pressure_data(Data *traded_pressure_data_pt)
Definition: elastic_bretherton.cc:1065
Definition: elastic_bretherton.cc:344
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
Definition: elastic_bretherton.cc:522
double *& re_invfr_pt()
Pointer to Reynolds number divided by Froude number.
Definition: elastic_bretherton.cc:513
Definition: bretherton_spine_mesh.template.h:46
void set_spine_centre_fraction_pt(double *const &fraction_pt)
Set the pointer to the spine centre's vertial fraction.
Definition: bretherton_spine_mesh.template.h:200
unsigned long nbulk() const
Number of elements in bulk.
Definition: bretherton_spine_mesh.template.h:99
FiniteElement *& interface_element_pt(const unsigned long &i)
Access functions for pointers to interface elements.
Definition: bretherton_spine_mesh.template.h:81
ELEMENT * control_element_pt()
Definition: bretherton_spine_mesh.template.h:187
unsigned long ninterface_element() const
Number of elements on interface.
Definition: bretherton_spine_mesh.template.h:87
Definition: nodes.h:86
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
double * value_pt(const unsigned &i) const
Definition: nodes.h:324
void set_value(const unsigned &i, const double &value_)
Definition: nodes.h:271
Definition: beam_elements.h:352
void set_normal_pointing_out_of_fluid()
Definition: beam_elements.h:384
double *& q_pt()
Definition: fsi.h:248
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
double *& ca_pt()
Pointer to the Capillary number.
Definition: interface_elements.h:492
void set_external_pressure_data(Data *external_pressure_data_pt)
Definition: interface_elements.h:539
Generalised element used to specify the mass flux (=1)
Definition: jh_includes.h:36
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
Definition: geom_objects.h:101
Definition: frontal_solver.h:56
Data * hijack_nodal_value(const unsigned &n, const unsigned &i, const bool &return_data=true)
Definition: hijacked_elements.h:214
void(*&)(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load) load_vector_fct_pt()
Reference to the load vector function pointer.
Definition: beam_elements.h:169
double *& sigma0_pt()
Return a pointer to axial prestress.
Definition: beam_elements.h:233
double *& h_pt()
Return a pointer to non-dim. wall thickness.
Definition: beam_elements.h:240
GeomObject *& undeformed_beam_pt()
Definition: beam_elements.h:253
Definition: mesh_as_geometric_object.h:93
Definition: mesh.h:67
unsigned long nboundary_node(const unsigned &ibound) const
Return number of nodes on a particular boundary.
Definition: mesh.h:833
int face_index_at_boundary(const unsigned &b, const unsigned &e) const
Definition: mesh.h:896
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Definition: mesh.h:878
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
Definition: mesh.h:840
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
void add_element_pt(GeneralisedElement *const &element_pt)
Add a (pointer to) an element to the mesh.
Definition: mesh.h:617
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Definition: mesh.h:493
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Definition: nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
void set_auxiliary_node_update_fct_pt(AuxNodeUpdateFctPt aux_node_update_fct_pt)
Definition: nodes.h:1596
virtual void set_coordinates_on_boundary(const unsigned &b, const unsigned &k, const Vector< double > &boundary_zeta)
Definition: nodes.cc:2394
unsigned add_sub_mesh(Mesh *const &mesh_pt)
Definition: problem.h:1330
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
Definition: problem.h:1466
void add_global_data(Data *const &global_data_pt)
Definition: problem.h:1654
void build_global_mesh()
Definition: problem.cc:1493
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Definition: problem.cc:1989
SolidNode * boundary_node_pt(const unsigned &b, const unsigned &n)
Return n-th SolidNodes on b-th boundary.
Definition: mesh.h:2612
void pin_position(const unsigned &i)
Pin the nodal position.
Definition: nodes.h:1816
Definition: specific_node_update_interface_elements.h:556
void node_update(const bool &update_all_solid_nodes=false)
Definition: spines.cc:84
Definition: spines.h:328
unsigned & node_update_fct_id()
Access function to ID of node update function (within specific mesh)
Definition: spines.h:384
Definition: geom_objects.h:452
static constexpr lastp1_t end
Definition: IndexedViewHelper.h:79
EIGEN_STRONG_INLINE const Eigen::CwiseBinaryOp< Eigen::internal::scalar_zeta_op< typename DerivedX::Scalar >, const DerivedX, const DerivedQ > zeta(const Eigen::ArrayBase< DerivedX > &x, const Eigen::ArrayBase< DerivedQ > &q)
Definition: SpecialFunctionsArrayAPI.h:152
double ReSt
Womersley number.
Definition: rayleigh_instability.cc:56
void spring_load_lower(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Load function for the wall elements.
Definition: elastic_bretherton.cc:216
double T
Period of oscillations.
Definition: elastic_bretherton.cc:113
GeomObject * Lower_wall_pt
Pointer to the lower wall.
Definition: elastic_bretherton.cc:143
void spring_load(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Load function for the wall elements.
Definition: elastic_bretherton.cc:160
GeomObject * Upper_wall_pt
Pointer to the upper wall.
Definition: elastic_bretherton.cc:140
double Q
Ratio of scales.
Definition: elastic_bretherton.cc:131
double Ca
Capillary number.
Definition: fibre.cc:61
map< double, pair< GeomObject *, Vector< double > > > lower_map
Lower wall map.
Definition: elastic_bretherton.cc:149
double ReInvFr
Product of Reynolds number and inverse of Froude number.
Definition: fibre.cc:58
double Re
Reynolds number.
Definition: fibre.cc:55
map< double, pair< GeomObject *, Vector< double > > > upper_map
Upper wall map.
Definition: elastic_bretherton.cc:146
Vector< double > G(3)
Direction of gravity.
Definition: spherical_shell_convection.cc:62
double H
Nondim thickness.
Definition: steady_ring.cc:50
void no_slip_condition_first(Node *node_pt)
Function that is used to set and update the no-slip boundary condition.
Definition: elastic_bretherton.cc:277
void no_slip_condition_second(Node *node_pt)
Definition: elastic_bretherton.cc:294

References oomph::Mesh::add_element_pt(), b, Global_Physical_Variables::Ca, oomph::FluidInterfaceElement::ca_pt(), e(), Eigen::placeholders::end, Global_Physical_Variables::G, SpineGravityTractionElement< ELEMENT >::g_pt(), Global_Physical_Variables::H, oomph::KirchhoffLoveBeamEquations::h_pt(), FixSpineHeightElement::height_pt(), oomph::Hijacked< ELEMENT >::hijack_nodal_value(), i, oomph::KirchhoffLoveBeamEquations::load_vector_fct_pt, Global_Physical_Variables::lower_map, Global_Physical_Variables::Lower_wall_pt, n, oomph::FiniteElement::nnode(), No_Slip::no_slip_condition_first(), No_Slip::no_slip_condition_second(), oomph::SpineNode::node_update_fct_id(), Global_Physical_Variables::Q, oomph::FSIWallElement::q_pt(), Global_Physical_Variables::Re, SpineGravityTractionElement< ELEMENT >::re_invfr_pt(), Global_Physical_Variables::ReInvFr, Global_Physical_Variables::ReSt, oomph::Node::set_auxiliary_node_update_fct_pt(), oomph::Node::set_coordinates_on_boundary(), oomph::FluidInterfaceElement::set_external_pressure_data(), oomph::FSIHermiteBeamElement::set_normal_pointing_out_of_fluid(), FixSpineHeightElement::set_traded_pressure_data(), oomph::KirchhoffLoveBeamEquations::sigma0_pt(), Global_Physical_Variables::spring_load(), Global_Physical_Variables::spring_load_lower(), Global_Physical_Variables::T, oomph::KirchhoffLoveBeamEquations::undeformed_beam_pt(), Global_Physical_Variables::upper_map, Global_Physical_Variables::Upper_wall_pt, oomph::Node::x(), and Eigen::zeta().

◆ ~AirwayReopeningProblem()

template<class ELEMENT >
AirwayReopeningProblem< ELEMENT >::~AirwayReopeningProblem ( )
inline

Destructor, clean up all allocated memory.

1224  {
1225  //Delete objects created in constructor in reverse order
1226  delete Constraint_mesh_pt;
1227  delete Bulk_mesh_pt;
1228  //Delete objects associated with lower wall mesh
1230  delete Lower_wall_mesh_pt;
1232  //Delete objects associated with upper wall mesh
1234  delete Upper_wall_mesh_pt;
1236 
1237  //Delete global data
1239  delete Bubble_pressure_data_pt;
1240  //Delete the linear solver, if allocated
1241  if(Frontal_solver) {delete linear_solver_pt();}
1242  }

References Global_Physical_Variables::Lower_wall_pt, and Global_Physical_Variables::Upper_wall_pt.

Member Function Documentation

◆ actions_after_change_in_global_parameter()

template<class ELEMENT >
void AirwayReopeningProblem< ELEMENT >::actions_after_change_in_global_parameter ( double *const &  parameter_pt)
inlinevirtual

Overload the continuation actions because we're continuing in Ca which does not affect the mesh

Reimplemented from oomph::Problem.

1248  {
1249  //Check that the multipliers have worked
1250  using namespace Global_Physical_Variables;
1251  ReInvFr = Bo/Ca;
1252  Re = ReCa*Ca;
1253  Q = Ca*Gamma;
1254  T = 100.0*Gamma/H;
1255  }
MatrixXf H
Definition: HessenbergDecomposition_matrixH.cpp:4
double ReInvFr
Product of Rynolds number and inverse of Froude number.
Definition: axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:67
double Re
Reynolds number.
Definition: axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:61
Global variables.
Definition: TwenteMeshGluing.cpp:60
double Bo
Bond number.
Definition: elastic_bretherton.cc:97
double Gamma
Definition: elastic_bretherton.cc:110
double ReCa
Reynolds divided by Capillary number.
Definition: elastic_bretherton.cc:91

References Global_Physical_Variables::Bo, Global_Physical_Variables::Ca, Global_Physical_Variables::Gamma, H, Global_Physical_Variables::Q, GlobalPhysicalVariables::Re, Global_Physical_Variables::ReCa, and GlobalPhysicalVariables::ReInvFr.

◆ actions_after_newton_solve()

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

Update after solve can remain empty, because the update is performed automatically after every Newton step.

Reimplemented from oomph::Problem.

1288 {}

◆ actions_before_newton_convergence_check()

template<class ELEMENT >
void AirwayReopeningProblem< ELEMENT >::actions_before_newton_convergence_check ( )
inlinevirtual

Spine heights/lengths are unknowns in the problem so their values get corrected during each Newton step. However, changing their value does not automatically change the nodal positions, so we need to update all of them

Reimplemented from oomph::Problem.

1262  {
1264 
1265  // Mesh Update
1267 
1268  // This driver code cannot be allowed to use the analytical form of
1269  // get_dresidual_dnodal_coordinates(...) that is implemented in the
1270  // NavierStokesEquations class, since the elemental residuals have
1271  // contributions from external data which is not taken into account
1272  // by that routine. We therefore force the bulk elements to use the
1273  // fully-finite differenced version.
1274  const unsigned n_bulk_element = Bulk_mesh_pt->nelement();
1275  for(unsigned e=0;e<n_bulk_element;e++)
1276  {
1277  ElementWithMovingNodes* el_pt =
1278  dynamic_cast<ElementWithMovingNodes*>(Bulk_mesh_pt->element_pt(e));
1280  }
1281  }
void actions_after_change_in_global_parameter(double *const &parameter_pt)
Definition: elastic_bretherton.cc:1246
Definition: element_with_moving_nodes.h:65
void evaluate_shape_derivs_by_direct_fd()
Evaluate shape derivatives by direct finite differencing.
Definition: element_with_moving_nodes.h:210

References Global_Physical_Variables::Ca, e(), and oomph::ElementWithMovingNodes::evaluate_shape_derivs_by_direct_fd().

◆ actions_before_newton_solve()

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

Update before solve: empty.

Reimplemented from oomph::Problem.

1284 {}

◆ change_boundary_conditions()

template<class ELEMENT >
void AirwayReopeningProblem< ELEMENT >::change_boundary_conditions ( )
inline
1387  {
1388  //Delete the traction elements and set what we really want
1390 
1391  //Construct the new elements
1393 
1394  //Rebuild the global mesh
1396 
1397  //Now set a consistent flux
1398  Flux_constraint_pt->set_flux(get_outlet_flux());
1399  //Unpin the flux constraint
1401 
1402  //Hijack the nodes again
1403  for(unsigned i=0;i<Outlet_traction_element_pt.size();i++)
1404  {
1405  dynamic_cast<SpineGravityTractionElement<ELEMENT>*>(
1406  Outlet_traction_element_pt[i])->hijack_all_nodes();
1407  }
1408 
1409  //Hijack the nodes again
1410  for(unsigned i=0;i<Gravity_traction_element_pt.size();i++)
1411  {
1412  dynamic_cast<SpineGravityTractionElement<ELEMENT>*>(
1413  Gravity_traction_element_pt[i])->hijack_all_nodes();
1414  }
1415 
1416 
1417  //Do equation numbering
1418  cout << "Number of unknowns: " << assign_eqn_numbers() << std::endl;
1419 
1420 
1421  //Sort the elements if using the frontal solver
1422  if(Frontal_solver)
1423  {
1424  std::sort(mesh_pt()->element_pt().begin(),
1425  mesh_pt()->element_pt().end(),
1426  ElementCmp());
1427  }
1428 
1429  }
double get_outlet_flux()
Definition: elastic_bretherton.cc:1369
void delete_outlet_elements()
Delete the outlet elements.
Definition: elastic_bretherton.cc:1331
Vector< FiniteElement * > Outlet_traction_element_pt
Vector of pointers to the outlet traciton elements.
Definition: elastic_bretherton.cc:1591
void construct_poisson_outlet_elements()
Set up the actual outlet elements.
Definition: elastic_bretherton.cc:1307
void unpin(const unsigned &i)
Unpin the i-th stored variable.
Definition: nodes.h:391
void rebuild_global_mesh()
Definition: problem.cc:1533

References Eigen::placeholders::end, and i.

◆ connect_walls()

template<class ELEMENT >
void AirwayReopeningProblem< ELEMENT >::connect_walls ( Mesh *const &  wall_mesh_pt,
GeomObject *const &  geom_object_pt,
map< double, pair< GeomObject *, Vector< double > > > &  globalmap 
)
inline

Setup the coupling between the upper and lower walls.

1489  {
1490  //Loop over the elements in the wall mesh
1491  unsigned n_wall_element = wall_mesh_pt->nelement();
1492  for(unsigned e=0;e<n_wall_element;e++)
1493  {
1494  //A set of geometric data to be added to the element
1495  set<Data*> extra_data;
1496 
1497  //Cast each element to an FSIWallElement
1498  FSIHermiteBeamElement *solid_element_pt =
1499  dynamic_cast<FSIHermiteBeamElement*>(wall_mesh_pt->element_pt(e));
1500 
1501  //Find the number of Gauss points of the element
1502  unsigned nintpt = solid_element_pt->integral_pt()->nweight();
1503  //Find the dimension of the element
1504  unsigned el_dim = solid_element_pt->dim();
1505  //Set storage for the local coordinates in the solid and bulk elements
1506  Vector<double> s(el_dim), s_bulk(el_dim);
1507  //Set storage for the zeta (intrinsic coordinate) at the Gauss point
1509  //Define storage for the geometric object that contains the Gauss point
1510  GeomObject *sub_obj_pt;
1511 
1512  //Loop over the integration points
1513  for(unsigned ipt=0;ipt<nintpt;ipt++)
1514  {
1515  //Loop over the dimension of the solid element and find the local
1516  //coordinates of the Gauss points
1517  for(unsigned i=0;i<el_dim;i++)
1518  {s[i] = solid_element_pt->integral_pt()->knot(ipt,i);}
1519  //Get the value of zeta at the integration point
1520  //Note that zeta coincides with xi (the Lagrangian coordinate)
1521  //in solid elements
1522  solid_element_pt->interpolated_xi(s,zeta);
1523 
1524  //Find the geometric object and local coordinate within that
1525  //object for the given value of the intrinsic coordinate, zeta.
1526  geom_object_pt->locate_zeta(zeta, sub_obj_pt, s_bulk);
1527 
1528  //Stick the result in the map
1529  globalmap[zeta[0]] = make_pair(sub_obj_pt,s_bulk);
1530 
1531  //Now sort add the external data to the set
1532  unsigned n_sub_geom_data = sub_obj_pt->ngeom_data();
1533  for(unsigned i=0;i<n_sub_geom_data;i++)
1534  {
1535  extra_data.insert(sub_obj_pt->geom_data_pt(i));
1536  }
1537  }
1538 
1539  //Add the data to the external data of the object
1540  for(set<Data*>::iterator it=extra_data.begin();
1541  it!=extra_data.end();++it)
1542  {
1543  solid_element_pt->add_external_data(*it);
1544  }
1545  }
1546  }
unsigned dim() const
Definition: elements.h:2611
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Definition: elements.cc:307
virtual unsigned ngeom_data() const
Definition: geom_objects.h:209
virtual void locate_zeta(const Vector< double > &zeta, GeomObject *&sub_geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
Definition: geom_objects.h:381
virtual Data * geom_data_pt(const unsigned &j)
Definition: geom_objects.h:233
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double interpolated_xi(const Vector< double > &s, const unsigned &i) const
Definition: elements.cc:7104
RealScalar s
Definition: level1_cplx_impl.h:130
unsigned el_dim
dimension
Definition: overloaded_cartesian_element_body.h:30

References oomph::GeneralisedElement::add_external_data(), oomph::FiniteElement::dim(), e(), el_dim, oomph::Mesh::element_pt(), oomph::GeomObject::geom_data_pt(), i, oomph::FiniteElement::integral_pt(), oomph::SolidFiniteElement::interpolated_xi(), oomph::Integral::knot(), oomph::GeomObject::locate_zeta(), oomph::Mesh::nelement(), oomph::GeomObject::ngeom_data(), oomph::Integral::nweight(), s, and Eigen::zeta().

◆ construct_poisson_outlet_elements()

template<class ELEMENT >
void AirwayReopeningProblem< ELEMENT >::construct_poisson_outlet_elements ( )
inline

Set up the actual outlet elements.

1308  {
1309  //Loop over the elements on the outlet boundary
1310  unsigned b=1;
1311  unsigned nbound_element = Bulk_mesh_pt->nboundary_element(b);
1312  for(unsigned e=0;e<nbound_element;e++)
1313  {
1314  SpineGravityTractionElement<ELEMENT>* outlet_element
1318 
1319  //Add the flux constraint
1320  outlet_element->set_delta_p_pt(Flux_constraint_pt->internal_data_pt(0));
1321 
1322  //Add the elements to the local storage
1323  Outlet_traction_element_pt.push_back(outlet_element);
1324 
1325  //Add the elements to the mesh
1326  Bulk_mesh_pt->add_element_pt(outlet_element);
1327  }
1328  }
void set_delta_p_pt(Data *const &delta_p_pt)
Access function to the pointer to pressure gradient data.
Definition: elastic_bretherton.cc:481

References b, e(), and SpineGravityTractionElement< ELEMENT >::set_delta_p_pt().

◆ construct_pressure_gradient_outlet_elements()

template<class ELEMENT >
void AirwayReopeningProblem< ELEMENT >::construct_pressure_gradient_outlet_elements ( )
inline

Set up the inital outlet elements.

1346  {
1347  //Loop over the elements next to the inflow boundarys
1348  unsigned b=1;
1349  unsigned nbound_element = Bulk_mesh_pt->nboundary_element(b);
1350  for(unsigned e=0;e<nbound_element;e++)
1351  {
1356 
1357  //Add the elements to the local storage
1358  Outlet_traction_element_pt.push_back(outlet_element);
1359 
1360  //Add the traction function
1361  outlet_element->traction_fct_pt() =
1363 
1364  //Add the elements to the mesh
1365  Bulk_mesh_pt->add_element_pt(outlet_element);
1366  }
1367  }
Definition: spines.h:477
void hydrostatic_pressure(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Function that prescribes the hydrostatic pressure field at the outlet.
Definition: elastic_bretherton.cc:152

References b, e(), and Global_Physical_Variables::hydrostatic_pressure().

◆ delete_outlet_elements()

template<class ELEMENT >
void AirwayReopeningProblem< ELEMENT >::delete_outlet_elements ( )
inline

Delete the outlet elements.

1332  {
1333  unsigned n_outlet = Outlet_traction_element_pt.size();
1334  for(unsigned e=n_outlet;e>0;e--)
1335  {
1336  //Remove and delete the outlet elements
1337  Bulk_mesh_pt->element_pt().pop_back();
1338  delete Outlet_traction_element_pt[e-1];
1340  }
1342  }

References e().

◆ doc_solution()

template<class ELEMENT >
void AirwayReopeningProblem< ELEMENT >::doc_solution ( DocInfo doc_info)

Doc the solution.

2049 {
2050 
2051  ofstream some_file;
2052  char filename[100];
2053 
2054  // Number of plot points
2055  unsigned npts=5;
2056 
2057  // Number of interface elements
2058 // unsigned ninterface=Bulk_mesh_pt->ninterface_element();
2059 
2060  // Number of spines
2061  //unsigned nspine=Bulk_mesh_pt->nspine();
2062 
2063  // Control coordinate: At bubble tip
2064  Vector<double> s(2);
2065  s[0]=1.0;
2066  s[1]=1.0;
2067 
2068  // Last proper spine
2069  unsigned last_spine=Bulk_mesh_pt->nfree_surface_spines()-1;
2070 
2071  // Doc
2076  Trace_file << " " << get_outlet_flux();
2077  Trace_file << " " <<
2079  Trace_file << " " <<
2081  Trace_file << " " <<
2083  Trace_file << " " <<
2086  Trace_file << " " << Bulk_mesh_pt->spine_pt(0)->height();
2087  Trace_file << " " << Bulk_mesh_pt->spine_pt(last_spine)->height();
2088  Trace_file << " " << 1.3375*pow(Global_Physical_Variables::Ca,2.0/3.0);
2089  Trace_file << " " <<
2090  (Bubble_pressure_data_pt->value(0)-Control_element_pt->interpolated_p_nst(s))*
2092  Trace_file << " " << 1.0+3.8*pow(Global_Physical_Variables::Ca,2.0/3.0);
2093  Trace_file << " " << Control_element_pt->interpolated_u_nst(s,0);
2094  Trace_file << " " << Control_element_pt->interpolated_u_nst(s,1);
2095  Trace_file << " "
2096 // << abs(dynamic_cast<SpineLineFluidInterfaceElement<ELEMENT>*>(
2097 // Bulk_mesh_pt->interface_element_pt(0))->
2098 // actual_contact_angle_left())*
2099 // 180.0/MathematicalConstants::Pi << " " ;
2100 // Trace_file << " "
2101 // << abs(dynamic_cast<SpineLineFluidInterfaceElement<ELEMENT>*>(
2102 // Bulk_mesh_pt->interface_element_pt(ninterface-1))->
2103 // actual_contact_angle_right())*
2104 // 180.0/MathematicalConstants::Pi
2105  << " ";
2106  //Trace_file << get_outlet_flux() << " ";
2107  Trace_file << std::endl;
2108 
2109 
2110  // Output solution
2111  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
2112  doc_info.number());
2113  some_file.open(filename);
2114  unsigned n_bulk=Bulk_mesh_pt->nbulk();
2115  for(unsigned i=0;i<n_bulk;i++)
2116  {
2117  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->
2118  bulk_element_pt(i));
2119  el_pt->output(some_file,npts);
2120  }
2121  some_file.close();
2122 
2123 }
ofstream Trace_file
Trace file.
Definition: elastic_bretherton.cc:1554
unsigned nfree_surface_spines()
Definition: bretherton_spine_mesh.template.h:107
double value(const unsigned &i) const
Definition: nodes.h:293
std::string directory() const
Output directory.
Definition: oomph_utilities.h:524
unsigned & number()
Number used (e.g.) for labeling output files.
Definition: oomph_utilities.h:554
Spine *& spine_pt(const unsigned long &i)
Return the i-th spine in the mesh.
Definition: spines.h:623
double & height()
Access function to spine height.
Definition: spines.h:150
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
double Tube_width
Natural width of the open tube.
Definition: elastic_bretherton.cc:125
string filename
Definition: MergeRestartFiles.py:39

References Global_Physical_Variables::Bo, Global_Physical_Variables::Ca, oomph::DocInfo::directory(), MergeRestartFiles::filename, i, oomph::DocInfo::number(), Eigen::bfloat16_impl::pow(), Global_Physical_Variables::Re, s, oomph::Problem_Parameter::Trace_file, and Global_Physical_Variables::Tube_width.

◆ dump()

template<class ELEMENT >
void AirwayReopeningProblem< ELEMENT >::dump ( ofstream &  dump_file) const
inlinevirtual

Dump the entire problem.

Reimplemented from oomph::Problem.

1433  {
1434  using namespace Global_Physical_Variables;
1435  dump_file << Re << std::endl;
1436  dump_file << ReSt << std::endl;
1437  dump_file << ReInvFr << std::endl;
1438  dump_file << Ca << std::endl;
1439  dump_file << Bo << std::endl;
1440  dump_file << Pext << std::endl;
1441  dump_file << G[0] << " " << G[1] << std::endl;
1442  dump_file << Nu << std::endl;
1443  dump_file << H << std::endl;
1444  dump_file << Gamma << std::endl;
1445  dump_file << T << std::endl;
1446  dump_file << Kstiff << std::endl;
1447  dump_file << Q << std::endl;
1448  dump_file << Kseparation << std::endl;
1449  dump_file << Tube_width << std::endl;
1450  dump_file << Ktable << std::endl;
1451 
1452  dump_file << Fixed_spine_height << std::endl;
1453  //Standard problem dumper
1454  Problem::dump(dump_file);
1455  }
JacobiRotation< float > G
Definition: Jacobi_makeGivens.cpp:2
double Nu
Definition: ConstraintElementsUnitTest.cpp:40
double ReSt
Product of Reynolds and Strouhal numbers.
Definition: axisym_navier_stokes/counter_rotating_disks/counter_rotating_disks.cc:64
double Pext
External Pressure.
Definition: elastic_bretherton.cc:99
double Kstiff
Spring stiffness.
Definition: elastic_bretherton.cc:116
double Kseparation
Huge stiffness for hitting.
Definition: elastic_bretherton.cc:119
double Ktable
Huge stiffness for the table.
Definition: elastic_bretherton.cc:122

References Global_Physical_Variables::Bo, Global_Physical_Variables::Ca, G, Global_Physical_Variables::Gamma, H, Global_Physical_Variables::Kseparation, Global_Physical_Variables::Kstiff, Global_Physical_Variables::Ktable, Constitutive_Parameters::Nu, Global_Physical_Variables::Pext, Global_Physical_Variables::Q, GlobalPhysicalVariables::Re, GlobalPhysicalVariables::ReInvFr, GlobalPhysicalVariables::ReSt, and Global_Physical_Variables::Tube_width.

◆ fix_pressure()

template<class ELEMENT >
void AirwayReopeningProblem< ELEMENT >::fix_pressure ( const unsigned e,
const unsigned l,
const double pvalue 
)
inline

Fix pressure value l in element e to value p_value.

1293  {
1294  //Fix the pressure at that element
1295  dynamic_cast<ELEMENT *>(Bulk_mesh_pt->element_pt(e))->
1296  fix_pressure(l,pvalue);
1297  }
void fix_pressure(const unsigned &e, const unsigned &l, const double &pvalue)
Fix pressure value l in element e to value p_value.
Definition: elastic_bretherton.cc:1291

References e().

◆ get_outlet_flux()

template<class ELEMENT >
double AirwayReopeningProblem< ELEMENT >::get_outlet_flux ( )
inline
1370  {
1371  if(dynamic_cast<SpineGravityTractionElement<ELEMENT>*>
1373  {return 0.0;}
1374 
1375  double flux = 0.0;
1376  //Loop over the outlet and calculate the flux
1377  for(unsigned e=0;e<Outlet_traction_element_pt.size();e++)
1378  {
1379  flux += dynamic_cast<SpineGravityTractionElement<ELEMENT>*>
1381  }
1382  return flux;
1383  }
void flux(const double &time, const Vector< double > &x, double &flux)
Get flux applied along boundary x=0.
Definition: pretend_melt.cc:59
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
Definition: gen_axisym_advection_diffusion_elements.h:424

References e(), ProblemParameters::flux(), and oomph::get_flux().

◆ parameter_study()

template<class ELEMENT >
void AirwayReopeningProblem< ELEMENT >::parameter_study ( const unsigned nsteps,
const bool restart 
)

Run a parameter study; perform specified number of steps.

Parameter study.

2134 {
2135  double flux = 0.0;
2136  // Increase maximum residual
2137  Problem::Max_residuals=100.0;
2139  //Problem::Newton_solver_tolerance=1.0e-6;
2140 
2141  // Set output directory
2142  DocInfo doc_info;
2143  doc_info.set_directory("RESLT");
2144  doc_info.number()=0;
2145 
2146  /*double zeta_trans_start_lo = -1.0;
2147  double zeta_trans_end_lo = 1.5;
2148  double fraction = 0.45;*/
2149 
2150  // Open trace file
2151  char filename[100], dumpfile[100];
2152  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
2153  Trace_file.open(filename);
2154 
2155  Trace_file << "VARIABLES=\"Ca\",";
2156  Trace_file << "\"h<sub>bottom</sub>\",\"h<sub>too</sub>\",";
2157  Trace_file << "\"h<sub>Bretherton</sub>\",\"p<sub>tip</sub>\",";
2158  Trace_file << "\"p<sub>tip (Bretherton)</sub>\",\"u<sub>stag</sub>\",";
2159  Trace_file << "\"v<sub>stag</sub>\"";
2160  Trace_file << "\"<greek>a</greek><sub>bottom</sub>\",";
2161  Trace_file << "\"<greek>a</greek><sub>top</sub>\"";
2162  Trace_file << std::endl;
2163 
2164  if(restart)
2165  {
2166  ifstream dumpfile("dump.dat");
2167  read(dumpfile);
2168  dumpfile.close();
2169  }
2170 
2171  // Initial scaling factor for Ca reduction
2172  //double factor=1.5;
2173 
2174  // Update so the initial domain shape is plotted properly
2176 
2177  //Doc initial solution
2178  doc_solution(doc_info);
2179 
2180  //Loop over the steps
2181  for(unsigned step=1;step<=nsteps;step++)
2182  {
2183  cout << std::endl << "STEP " << step << std::endl;
2184 
2185  //Solve
2186  cout << "Solving for capillary number: "
2188  << " Bo = "
2189  << Global_Physical_Variables::Bo << std::endl;
2190 
2191  if((!restart) || (step > 1))
2192  {
2193  using namespace Global_Physical_Variables;
2194  newton_solve();
2195  }
2196  else
2197  {
2198  cout << "Restart skipping" << std::endl;
2199  }
2200  cout << "Bubble pressure is " << Bubble_pressure_data_pt->value(0) << std::endl;
2201 
2202  // Doc solution
2203  doc_info.number()++;
2204  doc_solution(doc_info);
2205 
2206  //Dump if we're solving a real problem
2207  if(step > 1)
2208  {
2209  using namespace Global_Physical_Variables;
2210  sprintf(dumpfile,
2211  "%s/dump.Ca_%g.Bo_%g.Re_%g.H_%g.TW_%g_Frac_%g.Kstiff_%g.dat",
2212  doc_info.directory().c_str(),
2213  Ca,Bo,Re,H,Tube_width,
2215  ofstream dumpstream(dumpfile);
2216  dump(dumpstream);
2217  dumpstream.close();
2218  cout << "Actual flux is " << get_outlet_flux() << std::endl;
2219  }
2220 
2221  //First time round, unpin things
2222  if(step==1)
2223  {
2225  flux = Flux_constraint_pt->read_flux();
2226  }
2227  //Otherwise do our parameter study
2228  else
2229  {
2230  using namespace Global_Physical_Variables;
2231  //Now start to increase the flux
2232  //The computed flux is -0.389 or something so
2233  //if we start from -0.35 and subtract 0.05 we get 0.4
2234  //and then continue from the first time
2235  if(step==2) {flux = -0.35;}
2236 
2237  //Then worry about the flux
2238  if(flux > -2.0) {flux -= 0.05; Flux_constraint_pt->set_flux(flux);}
2239  else
2240  {
2241  if(flux < -2.0) {flux = -2.0; Flux_constraint_pt->set_flux(flux);}
2242 
2243  //Now increase the capillary number
2244  if(Ca < 0.1) {Ca += 0.01;}
2245  else
2246  {
2247  Ca = 0.1;
2248  //Finally reduce the spring stiffness
2249  Kstiff -= 0.5*1.0e-7;
2250  }
2251  }
2252  }
2253  }
2254 }
void doc_solution(DocInfo &doc_info)
Doc the solution.
Definition: elastic_bretherton.cc:2048
void read(ifstream &restart_file)
Read the entire problem.
Definition: elastic_bretherton.cc:1458
void change_boundary_conditions()
Definition: elastic_bretherton.cc:1386
void dump(ofstream &dump_file) const
Dump the entire problem.
Definition: elastic_bretherton.cc:1432
Definition: oomph_utilities.h:499
void set_directory(const std::string &directory)
Definition: oomph_utilities.cc:298
void newton_solve()
Use Newton method to solve the problem.
Definition: problem.cc:8783
Definition: restart2.cpp:8
unsigned Max_newton_iterations
Maximum number of newton iterations.
Definition: elements.cc:1654

References Global_Physical_Variables::Bo, Global_Physical_Variables::Ca, oomph::DocInfo::directory(), MergeRestartFiles::filename, ProblemParameters::flux(), H, Global_Physical_Variables::Kstiff, oomph::Locate_zeta_helpers::Max_newton_iterations, oomph::DocInfo::number(), GlobalPhysicalVariables::Re, Eigen::TensorSycl::internal::read(), oomph::DocInfo::set_directory(), oomph::Problem_Parameter::Trace_file, and Global_Physical_Variables::Tube_width.

◆ read()

template<class ELEMENT >
void AirwayReopeningProblem< ELEMENT >::read ( ifstream &  restart_file)
inlinevirtual

Read the entire problem.

Reimplemented from oomph::Problem.

1459  {
1460  using namespace Global_Physical_Variables;
1461  restart_file >> Re;
1462  restart_file >> ReSt;
1463  restart_file >> ReInvFr;
1464  restart_file >> Ca;
1465  restart_file >> Bo;
1466  restart_file >> Pext;
1467  restart_file >> G[0];
1468  restart_file >> G[1];
1469  restart_file >> Nu;
1470  restart_file >> H;
1471  restart_file >> Gamma;
1472  restart_file >> T;
1473  restart_file >> Kstiff;
1474  restart_file >> Q;
1475  restart_file >> Kseparation;
1476  restart_file >> Tube_width;
1477  restart_file >> Ktable;
1478 
1479  restart_file >> Fixed_spine_height;
1480 
1481  //Standard problem reader
1482  Problem::read(restart_file);
1483  }
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::enable_if_t< PacketLoad, PacketType > read(const TensorMapper &tensorMapper, const StorageIndex &NCIndex, const StorageIndex &CIndex, const StorageIndex &ld)
read, a template function used for loading the data from global memory. This function is used to guar...
Definition: TensorContractionSycl.h:162

References Global_Physical_Variables::Bo, Global_Physical_Variables::Ca, G, Global_Physical_Variables::Gamma, H, Global_Physical_Variables::Kseparation, Global_Physical_Variables::Kstiff, Global_Physical_Variables::Ktable, Constitutive_Parameters::Nu, Global_Physical_Variables::Pext, Global_Physical_Variables::Q, GlobalPhysicalVariables::Re, Eigen::TensorSycl::internal::read(), GlobalPhysicalVariables::ReInvFr, GlobalPhysicalVariables::ReSt, Global_Physical_Variables::T, and Global_Physical_Variables::Tube_width.

Member Data Documentation

◆ Bubble_pressure_data_pt

template<class ELEMENT >
Data* AirwayReopeningProblem< ELEMENT >::Bubble_pressure_data_pt
private

Data value that represents the bubble pressure.

◆ Bulk_mesh_pt

template<class ELEMENT >
BrethertonSpineMesh<ELEMENT, SpineLineFluidInterfaceElement<ELEMENT> >* AirwayReopeningProblem< ELEMENT >::Bulk_mesh_pt
private

Pointer to bulk mesh.

◆ Constraint_mesh_pt

template<class ELEMENT >
Mesh* AirwayReopeningProblem< ELEMENT >::Constraint_mesh_pt
private

Pointer to the constraint mesh.

◆ Control_element_pt

template<class ELEMENT >
ELEMENT* AirwayReopeningProblem< ELEMENT >::Control_element_pt
private

Pointer to control element.

◆ Fixed_spine_height

template<class ELEMENT >
double AirwayReopeningProblem< ELEMENT >::Fixed_spine_height
private

Fixed height value.

◆ Flux_constraint_pt

template<class ELEMENT >
FluxConstraint* AirwayReopeningProblem< ELEMENT >::Flux_constraint_pt
private

Flux constraint element.

◆ Frontal_solver

template<class ELEMENT >
bool AirwayReopeningProblem< ELEMENT >::Frontal_solver
private

Boolean flag to specify the use of the frontal solver.

◆ Gravity_traction_element_pt

template<class ELEMENT >
Vector<FiniteElement*> AirwayReopeningProblem< ELEMENT >::Gravity_traction_element_pt
private

Vector of pointers to the traction elements.

◆ Lower_wall_mesh_pt

template<class ELEMENT >
OneDLagrangianMesh<FSIHermiteBeamElement>* AirwayReopeningProblem< ELEMENT >::Lower_wall_mesh_pt
private

Pointer to lower wall mesh.

◆ Mesh_fraction_at_transition_pt

template<class ELEMENT >
Data* AirwayReopeningProblem< ELEMENT >::Mesh_fraction_at_transition_pt
private

Data value that represents the vertical fraction between the upper and lower walls of the centre of the spines at the point of transition of the spine mesh

◆ Outlet_traction_element_pt

template<class ELEMENT >
Vector<FiniteElement*> AirwayReopeningProblem< ELEMENT >::Outlet_traction_element_pt
private

Vector of pointers to the outlet traciton elements.

◆ Trace_file

template<class ELEMENT >
ofstream AirwayReopeningProblem< ELEMENT >::Trace_file
private

Trace file.

◆ Undeformed_lower_wall_geom_pt

template<class ELEMENT >
GeomObject* AirwayReopeningProblem< ELEMENT >::Undeformed_lower_wall_geom_pt
private

Pointer to the Undeformed wall.

◆ Undeformed_upper_wall_geom_pt

template<class ELEMENT >
GeomObject* AirwayReopeningProblem< ELEMENT >::Undeformed_upper_wall_geom_pt
private

Pointer to the Undeformed wall.

◆ Upper_wall_mesh_pt

template<class ELEMENT >
OneDLagrangianMesh<FSIHermiteBeamElement>* AirwayReopeningProblem< ELEMENT >::Upper_wall_mesh_pt
private

Pointer to Wall mesh.


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