Public Member Functions | List of all members
TwoDDGProblem< ELEMENT > Class Template Reference
+ Inheritance diagram for TwoDDGProblem< ELEMENT >:

Public Member Functions

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

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

Constructor & Destructor Documentation

◆ TwoDDGProblem() [1/5]

template<class ELEMENT >
TwoDDGProblem< ELEMENT >::TwoDDGProblem ( const unsigned &  Nx,
const unsigned &  Ny 
)
inline
244  {
245  this->set_explicit_time_stepper_pt(new RungeKutta<4>);
246  this->disable_info_in_newton_solve();
247 
248  this->enable_discontinuous_formulation();
249 
250  Problem::mesh_pt() = new TwoDDGMesh<ELEMENT>(Nx,Ny);
251 
252  unsigned n_element = Problem::mesh_pt()->nelement();
253  for(unsigned e=0;e<n_element;e++)
254  {
255  ELEMENT* cast_element_pt =
256  dynamic_cast<ELEMENT*>(Problem::mesh_pt()->element_pt(e));
257 
258  cast_element_pt->wind_fct_pt() = &Global::two_d_wind;
259  }
260 
261 
262  std::cout << "How many " << assign_eqn_numbers() << std::endl;
263  }
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
TwoDDGMesh
Definition: two_d_advection.cc:66
oomph::Problem::set_explicit_time_stepper_pt
void set_explicit_time_stepper_pt(ExplicitTimeStepper *const &explicit_time_stepper_pt)
Definition: problem.cc:1586
oomph::Problem::assign_eqn_numbers
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Definition: problem.cc:1989
oomph::Problem::disable_info_in_newton_solve
void disable_info_in_newton_solve()
Disable the output of information when in the newton solver.
Definition: problem.h:2941
oomph::Problem::enable_discontinuous_formulation
void enable_discontinuous_formulation()
Definition: problem.h:1727
oomph::RungeKutta
Definition: explicit_timesteppers.h:187
GlobalParameters::Nx
unsigned Nx
Number of elements in each direction (used by SimpleCubicMesh)
Definition: structured_cubic_point_source.cc:114
GlobalParameters::Ny
unsigned Ny
Definition: structured_cubic_point_source.cc:115
Global::two_d_wind
void two_d_wind(const oomph::Vector< double > &x, oomph::Vector< double > &wind)
Definition: two_d_advection.cc:47

References e(), oomph::Mesh::element_pt(), oomph::Problem::mesh_pt(), oomph::Mesh::nelement(), GlobalParameters::Nx, GlobalParameters::Ny, and Global::two_d_wind().

◆ ~TwoDDGProblem() [1/4]

template<class ELEMENT >
TwoDDGProblem< ELEMENT >::~TwoDDGProblem ( )
inline
381  {
382  delete Problem::mesh_pt();
383  delete this->explicit_time_stepper_pt();
384  }
oomph::Problem::explicit_time_stepper_pt
ExplicitTimeStepper *& explicit_time_stepper_pt()
Return a pointer to the explicit timestepper.
Definition: problem.h:1555

References oomph::Problem::mesh_pt().

◆ TwoDDGProblem() [2/5]

template<class ELEMENT >
TwoDDGProblem< ELEMENT >::TwoDDGProblem ( const unsigned &  Nx,
const unsigned &  Ny 
)
inline
388  {
389  this->set_explicit_time_stepper_pt(new RungeKutta<4>);
390  this->disable_info_in_newton_solve();
391 
392  this->enable_discontinuous_formulation();
393 
394  Problem::mesh_pt() = new TwoDDGMesh<ELEMENT>(Nx,Ny);
395 
396  unsigned n_element = Problem::mesh_pt()->nelement();
397  for(unsigned e=0;e<n_element;e++)
398  {
399  ELEMENT* cast_element_pt =
400  dynamic_cast<ELEMENT*>(Problem::mesh_pt()->element_pt(e));
401 
402  cast_element_pt->gamma_pt() = &Global::Gamma;
403 
404  }
405 
406  //Set the boundary conditions
407  Vector<double> u_initial(4);
408  Vector<double> x(2);
409  //Pin all boundaries
410  for(unsigned b=0;b<2;b++)
411  {
412  const unsigned n_node = mesh_pt()->nboundary_node(b);
413  for(unsigned n=0;n<n_node;n++)
414  {
415  Node* nod_pt = mesh_pt()->boundary_node_pt(b,n);
416  const unsigned n_value = nod_pt->nvalue();
417  for(unsigned i=0;i<n_value;i++) {nod_pt->pin(i);}
418  //Now set the values
419  for(unsigned i=0;i<2;i++) {x[i] = nod_pt->x(i);}
420  Global::exact_solution(0.0,x,u_initial);
421  for(unsigned i=0;i<n_value;i++) {nod_pt->set_value(i,u_initial[i]);}
422  }
423  }
424 
425 
426  std::cout << "How many " << assign_eqn_numbers() << std::endl;
427  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
b
Scalar * b
Definition: benchVecAdd.cpp:17
TwoDDGProblem::mesh_pt
TwoDDGMesh< ELEMENT > * mesh_pt()
Definition: two_d_advection.cc:239
oomph::Data::pin
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
oomph::Data::set_value
void set_value(const unsigned &i, const double &value_)
Definition: nodes.h:271
oomph::Data::nvalue
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
oomph::Node
Definition: nodes.h:906
oomph::Node::x
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
Global::Gamma
double Gamma
Definition: couette.cc:44
Global::exact_solution
void exact_solution(const double &t, const Vector< double > &x, Vector< double > &u)
Function that determines the initial conditions.
Definition: couette.cc:47
plotDoE.x
list x
Definition: plotDoE.py:28

References b, e(), oomph::Mesh::element_pt(), Global::exact_solution(), Global::Gamma, i, oomph::Problem::mesh_pt(), n, oomph::Mesh::nelement(), oomph::Data::nvalue(), GlobalParameters::Nx, GlobalParameters::Ny, oomph::Data::pin(), oomph::Data::set_value(), plotDoE::x, and oomph::Node::x().

◆ ~TwoDDGProblem() [2/4]

template<class ELEMENT >
TwoDDGProblem< ELEMENT >::~TwoDDGProblem ( )
inline
849  {
850  delete Problem::mesh_pt();
851  delete this->explicit_time_stepper_pt();
852  }

References oomph::Problem::mesh_pt().

◆ TwoDDGProblem() [3/5]

template<class ELEMENT >
TwoDDGProblem< ELEMENT >::TwoDDGProblem ( const unsigned &  Nx,
const unsigned &  Ny 
)
inline
856  {
857  this->set_explicit_time_stepper_pt(new SSP_RungeKutta<2>);
858  this->disable_info_in_newton_solve();
859 
860  this->enable_discontinuous_formulation();
861 
862  Problem::mesh_pt() = new TwoDDGMesh<ELEMENT>(Nx,Ny);
863 
864  unsigned n_element = Problem::mesh_pt()->nelement();
865  for(unsigned e=0;e<n_element;e++)
866  {
867  ELEMENT* cast_element_pt =
868  dynamic_cast<ELEMENT*>(Problem::mesh_pt()->element_pt(e));
869 
870  cast_element_pt->gamma_pt() = &Global::Gamma;
871 
872  }
873 
874  //Pin the inlet boundary
875  {
876  unsigned b=3;
877  {
878  const unsigned n_node = mesh_pt()->nboundary_node(b);
879  for(unsigned n=0;n<n_node;n++)
880  {
881  Node* nod_pt = mesh_pt()->boundary_node_pt(b,n);
882  const unsigned n_value = nod_pt->nvalue();
883  for(unsigned i=0;i<n_value;i++) {nod_pt->pin(i);}
884  }
885  }
886  }
887 
888 
889  std::cout << "How many " << assign_eqn_numbers() << std::endl;
890  }
SSP_RungeKutta
Definition: ff_step.cc:45

References b, e(), oomph::Mesh::element_pt(), Global::Gamma, i, oomph::Problem::mesh_pt(), n, oomph::Mesh::nelement(), oomph::Data::nvalue(), GlobalParameters::Nx, GlobalParameters::Ny, and oomph::Data::pin().

◆ ~TwoDDGProblem() [3/4]

template<class ELEMENT >
TwoDDGProblem< ELEMENT >::~TwoDDGProblem ( )
inline
486  {
487  delete Problem::mesh_pt();
488  delete this->explicit_time_stepper_pt();
489  }

References oomph::Problem::mesh_pt().

◆ TwoDDGProblem() [4/5]

template<class ELEMENT >
TwoDDGProblem< ELEMENT >::TwoDDGProblem ( const unsigned &  Nx,
const unsigned &  Ny 
)
inline
493  {
494  this->set_explicit_time_stepper_pt(new RungeKutta<4>);
495  this->disable_info_in_newton_solve();
496 
497  this->enable_discontinuous_formulation();
498 
499  Problem::mesh_pt() = new TwoDDGMesh<ELEMENT>(Nx,Ny);
500 
501  unsigned n_element = Problem::mesh_pt()->nelement();
502  for(unsigned e=0;e<n_element;e++)
503  {
504  ELEMENT* cast_element_pt =
505  dynamic_cast<ELEMENT*>(Problem::mesh_pt()->element_pt(e));
506 
507  cast_element_pt->gamma_pt() = &Global::Gamma;
508 
509  }
510 
511  //Pin the inlet boundary
512  {
513  unsigned b=3;
514  {
515  const unsigned n_node = mesh_pt()->nboundary_node(b);
516  for(unsigned n=0;n<n_node;n++)
517  {
518  Node* nod_pt = mesh_pt()->boundary_node_pt(b,n);
519  const unsigned n_value = nod_pt->nvalue();
520  for(unsigned i=0;i<n_value;i++) {nod_pt->pin(i);}
521  }
522  }
523  }
524 
525 
526  std::cout << "How many " << assign_eqn_numbers() << std::endl;
527  }

References b, e(), oomph::Mesh::element_pt(), Global::Gamma, i, oomph::Problem::mesh_pt(), n, oomph::Mesh::nelement(), oomph::Data::nvalue(), GlobalParameters::Nx, GlobalParameters::Ny, and oomph::Data::pin().

◆ ~TwoDDGProblem() [4/4]

template<class ELEMENT >
TwoDDGProblem< ELEMENT >::~TwoDDGProblem ( )
inline
410  {
411  delete Problem::mesh_pt();
412  delete this->explicit_time_stepper_pt();
413  }

References oomph::Problem::mesh_pt().

◆ TwoDDGProblem() [5/5]

template<class ELEMENT >
TwoDDGProblem< ELEMENT >::TwoDDGProblem ( const unsigned &  Nx,
const unsigned &  Ny 
)
inline
417  {
418  this->set_explicit_time_stepper_pt(new RungeKutta<4>);
419  this->disable_info_in_newton_solve();
420 
421  this->enable_discontinuous_formulation();
422 
423  Problem::mesh_pt() = new TwoDDGMesh<ELEMENT>(Nx,Ny);
424 
425  unsigned n_element = Problem::mesh_pt()->nelement();
426  for(unsigned e=0;e<n_element;e++)
427  {
428  ELEMENT* cast_element_pt =
429  dynamic_cast<ELEMENT*>(Problem::mesh_pt()->element_pt(e));
430 
431  cast_element_pt->gamma_pt() = &Global::Gamma;
432 
433  }
434 
435  //Set the boundary conditions
436  Vector<double> u_initial(4);
437  Vector<double> x(2);
438  //Pin all boundaries
439  for(unsigned b=0;b<4;b++)
440  {
441  const unsigned n_node = mesh_pt()->nboundary_node(b);
442  for(unsigned n=0;n<n_node;n++)
443  {
444  Node* nod_pt = mesh_pt()->boundary_node_pt(b,n);
445  const unsigned n_value = nod_pt->nvalue();
446  for(unsigned i=0;i<n_value;i++) {nod_pt->pin(i);}
447  //Now set the values
448  for(unsigned i=0;i<2;i++) {x[i] = nod_pt->x(i);}
449  Global::exact_solution(0.0,x,u_initial);
450  for(unsigned i=0;i<n_value;i++) {nod_pt->set_value(i,u_initial[i]);}
451  }
452  }
453 
454 
455  std::cout << "How many " << assign_eqn_numbers() << std::endl;
456  }

References b, e(), oomph::Mesh::element_pt(), Global::exact_solution(), Global::Gamma, i, oomph::Problem::mesh_pt(), n, oomph::Mesh::nelement(), oomph::Data::nvalue(), GlobalParameters::Nx, GlobalParameters::Ny, oomph::Data::pin(), oomph::Data::set_value(), plotDoE::x, and oomph::Node::x().

Member Function Documentation

◆ actions_after_explicit_stage()

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::actions_after_explicit_stage ( )
inlinevirtual

Empty virtual function that should be overloaded to update any dependent data or boundary conditions that should be advanced after each stage of an explicit time step (Runge-Kutta steps contain multiple stages per time step, most others only contain one).

Reimplemented from oomph::ExplicitTimeSteppableObject.

505  {
506  apply_boundary_conditions(this->time());
507  }
TwoDDGProblem::apply_boundary_conditions
void apply_boundary_conditions(const double &t)
Definition: two_d_euler.cc:483
oomph::Problem::time
double & time()
Return the current value of continuous time.
Definition: problem.cc:11531

◆ apply_boundary_conditions()

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::apply_boundary_conditions ( const double &  t)
inline
484  {
485  //Set the boundary conditions
486  Vector<double> u_initial(4);
487  Vector<double> x(2);
488  //Pin all boundaries
489  for(unsigned b=0;b<4;b++)
490  {
491  const unsigned n_node = mesh_pt()->nboundary_node(b);
492  for(unsigned n=0;n<n_node;n++)
493  {
494  Node* nod_pt = mesh_pt()->boundary_node_pt(b,n);
495  const unsigned n_value = nod_pt->nvalue();
496  //Now set the values
497  for(unsigned i=0;i<2;i++) {x[i] = nod_pt->x(i);}
498  Global::exact_solution(t,x,u_initial);
499  for(unsigned i=0;i<n_value;i++) {nod_pt->set_value(i,u_initial[i]);}
500  }
501  }
502  }
plotPSD.t
t
Definition: plotPSD.py:36

References b, Global::exact_solution(), i, n, oomph::Data::nvalue(), oomph::Data::set_value(), plotPSD::t, plotDoE::x, and oomph::Node::x().

◆ apply_initial_conditions() [1/5]

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::apply_initial_conditions ( )
inline
288  {
289  //Storage for the coordinates
290  Vector<double> x(2);
291  //Storage for the initial condition
292  Vector<double> initial_u(1);
293  //Loop over all the nodes in the mesh
294  unsigned n_node = mesh_pt()->nnode();
295  for(unsigned n=0;n<n_node;n++)
296  {
297  x[0] = mesh_pt()->node_pt(n)->x(0);
298  x[1] = mesh_pt()->node_pt(n)->x(1);
299  //Get the initial conditions
300  Global::initial_condition(0.0,x,initial_u);
301 
302  mesh_pt()->node_pt(n)->set_value(0,initial_u[0]);
303  }
304  }
Global::initial_condition
void initial_condition(const double &time, const Vector< double > &x, Vector< double > &u)
Function that determines the initial conditions.
Definition: one_d_advection.cc:51

References Global::initial_condition(), n, and plotDoE::x.

◆ apply_initial_conditions() [2/5]

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::apply_initial_conditions ( )
inline
454  {
455  //Storage for the coordinates
456  Vector<double> x(2);
457  //Storage for the initial condition
458  Vector<double> initial_u(4);
459  //Loop over all the nodes in the mesh
460  unsigned n_node = mesh_pt()->nnode();
461  for(unsigned n=0;n<n_node;n++)
462  {
463  Node* nod_pt = mesh_pt()->node_pt(n);
464  x[0] = nod_pt->x(0);
465  x[1] = nod_pt->x(1);
466  //Get the initial conditions
467  Global::exact_solution(0.0,x,initial_u);
468  //Set them
469  for(unsigned i=0;i<4;i++)
470  {
471  nod_pt->set_value(i,initial_u[i]);
472  }
473  }
474  }

References Global::exact_solution(), i, n, oomph::Data::set_value(), plotDoE::x, and oomph::Node::x().

◆ apply_initial_conditions() [3/5]

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::apply_initial_conditions ( )
inline
893  {
894  const double gamma = Global::Gamma;
895  const double E = 1.0/(gamma-1.0) + 4.5*gamma;
896  //Loop over all the nodes in the mesh
897  unsigned n_node = mesh_pt()->nnode();
898  for(unsigned n=0;n<n_node;n++)
899  {
900  Node* nod_pt = mesh_pt()->node_pt(n);
901  //Set the values
902  nod_pt->set_value(0,gamma);
903  nod_pt->set_value(1,E);
904  nod_pt->set_value(2,3*gamma);
905  nod_pt->set_value(3,0.0);
906  }
907  }
Global_Physical_Variables::E
double E
Elastic modulus.
Definition: TwenteMeshGluing.cpp:68
mathsFunc::gamma
Mdouble gamma(Mdouble gamma_in)
This is the gamma function returns the true value for the half integer value.
Definition: ExtendedMath.cc:116

References Global_Physical_Variables::E, mathsFunc::gamma(), Global::Gamma, n, and oomph::Data::set_value().

◆ apply_initial_conditions() [4/5]

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::apply_initial_conditions ( )
inline
530  {
531  double M = 2.0;
532  const double gamma = Global::Gamma;
533  const double E = 1.0/(gamma-1.0) + 0.5*M*M*gamma;
534  //Loop over all the nodes in the mesh
535  unsigned n_node = mesh_pt()->nnode();
536  for(unsigned n=0;n<n_node;n++)
537  {
538  Node* nod_pt = mesh_pt()->node_pt(n);
539  //Set the values
540  nod_pt->set_value(0,gamma);
541  nod_pt->set_value(1,E);
542  nod_pt->set_value(2,M*gamma);
543  nod_pt->set_value(3,0.0);
544  }
545  }
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: Eigen/Eigen/src/Core/Matrix.h:186

References Global_Physical_Variables::E, mathsFunc::gamma(), Global::Gamma, n, and oomph::Data::set_value().

◆ apply_initial_conditions() [5/5]

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::apply_initial_conditions ( )
inline
511  {
512  //Storage for the coordinates
513  Vector<double> x(2);
514  //Storage for the initial condition
515  Vector<double> initial_u(4);
516  //Loop over all the nodes in the mesh
517  unsigned n_node = mesh_pt()->nnode();
518  for(unsigned n=0;n<n_node;n++)
519  {
520  Node* nod_pt = mesh_pt()->node_pt(n);
521  x[0] = nod_pt->x(0);
522  x[1] = nod_pt->x(1);
523  //Get the initial conditions
524  Global::exact_solution(0.0,x,initial_u);
525  //Set them
526  for(unsigned i=0;i<4;i++)
527  {
528  nod_pt->set_value(i,initial_u[i]);
529  }
530  }
531  }

References Global::exact_solution(), i, n, oomph::Data::set_value(), plotDoE::x, and oomph::Node::x().

◆ compute_error() [1/3]

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::compute_error ( const double &  t,
Vector< double > &  error 
)
inline

Compute the complete errors in the problem.

267  {
268  error[0] = 0.0;
269 
270  Vector<double> local_error(1);
271  Vector<double> local_norm(1);
272 
273  const unsigned n_element = Problem::mesh_pt()->nelement();
274  //Do the timestep
275  for(unsigned e=0;e<n_element;e++)
276  {
277  dynamic_cast<ScalarAdvectionEquations<2>*>(
278  Problem::mesh_pt()->element_pt(e))
279  ->compute_error(std::cout,
280  Global::initial_condition,t,local_error,
281  local_norm);
282  error[0] += local_error[0];
283  }
284  }
TwoDDGProblem::compute_error
void compute_error(const double &t, Vector< double > &error)
Compute the complete errors in the problem.
Definition: two_d_advection.cc:266
oomph::ScalarAdvectionEquations
Base class for advection equations.
Definition: scalar_advection_elements.h:48
calibrate.error
int error
Definition: calibrate.py:297

References oomph::compute_error(), e(), oomph::Mesh::element_pt(), calibrate::error, Global::initial_condition(), oomph::Problem::mesh_pt(), oomph::Mesh::nelement(), and plotPSD::t.

◆ compute_error() [2/3]

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::compute_error ( const double &  t,
Vector< double > &  error 
)
inline

Compute the complete errors in the problem.

431  {
432  error.initialise(0.0);
433 
434  Vector<double> local_error(4);
435  Vector<double> local_norm(4);
436 
437  const unsigned n_element = Problem::mesh_pt()->nelement();
438  //Do the timestep
439  for(unsigned e=0;e<n_element;e++)
440  {
441  dynamic_cast<ELEMENT*>(
442  Problem::mesh_pt()->element_pt(e))
443  ->compute_error(std::cout,
444  Global::exact_solution,t,local_error,
445  local_norm);
446  for(unsigned i=0;i<4;i++)
447  {
448  error[i] += local_error[i];
449  }
450  }
451  }

References oomph::compute_error(), e(), oomph::Mesh::element_pt(), calibrate::error, Global::exact_solution(), i, oomph::Problem::mesh_pt(), oomph::Mesh::nelement(), and plotPSD::t.

◆ compute_error() [3/3]

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::compute_error ( const double &  t,
Vector< double > &  error 
)
inline

Compute the complete errors in the problem.

460  {
461  error.initialise(0.0);
462 
463  Vector<double> local_error(4);
464  Vector<double> local_norm(4);
465 
466  const unsigned n_element = Problem::mesh_pt()->nelement();
467  //Do the timestep
468  for(unsigned e=0;e<n_element;e++)
469  {
470  dynamic_cast<ELEMENT*>(
471  Problem::mesh_pt()->element_pt(e))
472  ->compute_error(std::cout,
473  Global::exact_solution,t,local_error,
474  local_norm);
475  for(unsigned i=0;i<4;i++)
476  {
477  error[i] += local_error[i];
478  }
479  }
480  }

References oomph::compute_error(), e(), oomph::Mesh::element_pt(), calibrate::error, Global::exact_solution(), i, oomph::Problem::mesh_pt(), oomph::Mesh::nelement(), and plotPSD::t.

◆ limit() [1/2]

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::limit ( )
inline
910  {
911  const double eps = 1.0e-14;
912  //First stage is to calculate the averages of all elements
913  const unsigned n_element = this->mesh_pt()->nelement();
914  for(unsigned e=0;e<n_element;e++)
915  {
916  //Calculate the averages
917  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))
918  ->allocate_memory_for_averages();
919  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))->calculate_averages();
920  }
921 
922  Vector<double> s(2,0.0), x(2);
923 
924  const unsigned n_flux = 4;
925 
926  const unsigned n_dim = 2;
927 
928 
929  //Store pointers to neighbours
930  DenseMatrix<ELEMENT*> bulk_neighbour(n_element,4);
931  //Storage for the unknowns
932  Vector<double> interpolated_u(n_flux);
933 
934  //Now we need to actually do the limiting, but let's check things first
935  for(unsigned e=0;e<n_element;e++)
936  {
937  //Store the element
938  ELEMENT* cast_element_pt =
939  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
940  //Find the centre of the element
941  cast_element_pt->interpolated_x(s,x);
942 
943 
944  //Storage for the approximation to the central gradient from all four
945  //faces
946  RankThreeTensor<double> grad(4,n_flux,2);
947  //Storage for the neighbours size
948  Vector<double> size(4);
949 
950  //Get the local coordinate of the nodes
951  Vector<double> s_face(1), s_bulk(2);
952  //Get the central value
953  DenseMatrix<double> prim_centre(2,n_flux);
954  //My Centre
955  for(unsigned i=0;i<n_flux;i++)
956  {interpolated_u[i] = cast_element_pt->average_value(i);}
957  prim_centre(0,0) = interpolated_u[0];
958  prim_centre(0,1) = cast_element_pt->pressure(interpolated_u);
959  prim_centre(0,2) = interpolated_u[2]/interpolated_u[0];
960  prim_centre(0,3) = interpolated_u[3]/interpolated_u[0];
961 
962  //Now store the average primitivee values
963  for(unsigned i=0;i<n_flux;i++)
964  {cast_element_pt->average_prim_value(i) = prim_centre(0,i);}
965 
966  //Loop over all the faces
967  for(unsigned f=0;f<4;f++)
968  {
969  DGFaceElement* face_element_pt =
970  dynamic_cast<DGFaceElement*>(cast_element_pt->face_element_pt(f));
971 
972  //Get the average of the nodal points on the face
973  //(***NOT GENERAL**)
974  const unsigned n_face_node = face_element_pt->nnode();
975 
976  //Storage for the average nodal values
977  DenseMatrix<double> face_average(n_face_node,n_flux);
978 
979  //Storage for the neighbouring cell's centre average
980  DenseMatrix<double> cell_average(n_face_node,n_flux);
981 
982  //Storage for the sizes of the neighbouring bulk elements
983  Vector<double> neighbour_size(n_face_node);
984 
985  //Now the storage for the nodal values and centre values
986  DenseMatrix<double> nodal_x(n_face_node,n_dim);
987  //Now the storage for the centre values
988  DenseMatrix<double> centre_x(n_face_node,n_dim);
989 
990  //Neighbours face
991  FaceElement* neighbour_face_pt=0;
992 
993  //Now calculate those average values
994  for(unsigned n=0;n<n_face_node;n++)
995  {
996  //Get the local coordinates of the node
997  face_element_pt->local_coordinate_of_node(n,s_face);
998  //Get the global coordinate of the node
999  for(unsigned i=0;i<n_dim;i++)
1000  {
1001  nodal_x(n,i) = face_element_pt->interpolated_x(s_face,i);
1002  }
1003 
1004  //Get the flux at the node in this face
1005  face_element_pt->interpolated_u(s_face,interpolated_u);
1006 
1007  /*std::cout << "Values at face :";
1008  for(unsigned i=0;i<n_flux;i++)
1009  {
1010  std::cout << interpolated_u[i] << " ";
1011  }
1012  std::cout << std::endl;*/
1013 
1014  for(unsigned i=0;i<n_flux;i++)
1015  {
1016  face_average(n,i) = interpolated_u[i];
1017  }
1018 
1019  //Now get the bulk coordinate
1020  face_element_pt->get_local_coordinate_in_bulk(s_face,s_bulk);
1021  //Now get the neighbour
1022  cast_element_pt->
1023  get_neighbouring_face_and_local_coordinate(
1024  face_element_pt->face_index(),s_bulk,neighbour_face_pt,s_face);
1025 
1026  //Get the fluxes of the neighbour
1027  dynamic_cast<DGFaceElement*>(neighbour_face_pt)
1028  ->interpolated_u(s_face,interpolated_u);
1029  //Add the flux from the neighbour
1030  ELEMENT* neighbour_bulk_element_pt = dynamic_cast<ELEMENT*>(
1031  neighbour_face_pt->bulk_element_pt());
1032 
1033  /*std::cout << "Values at face neighbour :";
1034  for(unsigned i=0;i<n_flux;i++)
1035  {
1036  std::cout << interpolated_u[i] << " ";
1037  }
1038  std::cout << std::endl;
1039 
1040  std::cout << "SANITY CHECK " << face_element_pt << " "
1041  << neighbour_face_pt << " "
1042  << cast_element_pt << " "
1043  << neighbour_bulk_element_pt << "\n";*/
1044 
1045 
1046  //Get the centre of the neighbouring bulk element
1047  for(unsigned i=0;i<n_dim;i++) {s_bulk[i] = 0.0;}
1048  //Do it
1049  for(unsigned i=0;i<n_dim;i++)
1050  {
1051  centre_x(n,i) = neighbour_bulk_element_pt
1052  ->interpolated_x(s_bulk,i);
1053  }
1054 
1055  for(unsigned i=0;i<n_flux;i++)
1056  {
1057  face_average(n,i) += interpolated_u[i];
1058  face_average(n,i) *= 0.5;
1059  //Get the neighbour
1060  cell_average(n,i) = neighbour_bulk_element_pt->average_value(i);
1061  }
1062 
1063  //Get the size of the neighbour
1064  neighbour_size[n] = neighbour_bulk_element_pt->average_value(n_flux);
1065 
1066  //Only bother to store the first one
1067  if(n==0)
1068  {
1069  bulk_neighbour(e,f) = neighbour_bulk_element_pt;
1070  }
1071  }
1072 
1073  //Now we have the face averages and the cell averages
1074  //Check that we have a conforming mesh
1075  for(unsigned n=1;n<n_face_node;n++)
1076  {
1077  if((std::abs(centre_x(0,0) - centre_x(n,0)) > eps) ||
1078  (std::abs(centre_x(0,1) - centre_x(n,1)) > eps))
1079  {
1080  throw OomphLibError("Mesh is not conforming in limiter",
1081  OOMPH_CURRENT_FUNCTION,
1082  OOMPH_EXCEPTION_LOCATION);
1083  }
1084  }
1085 
1086 
1087  //Set centre_x(0,0) to be the current element
1088  for(unsigned i=0;i<n_dim;i++)
1089  {
1090  centre_x(0,i) = x[i];
1091  }
1092 
1093  //Set the size of the neighbour
1094  size[f] = neighbour_size[0];
1095 
1096  //Convert our values to primitive variables
1097  DenseMatrix<double> prim_nodal(2,n_flux);
1098 
1099  //Face averages (corner nodes)
1100  //std::cout << "Nodal values \n";
1101  for(unsigned n=0;n<2;n++)
1102  {
1103  //Load the fluxes
1104  for(unsigned i=0;i<n_flux;i++)
1105  {interpolated_u[i] = face_average(n,i);}
1106 
1107  //Output the fluxes
1108  /*std::cout << n << " : ";
1109  for(unsigned i=0;i<n_flux;i++)
1110  {std::cout << interpolated_u[i] << " ";}
1111  std::cout << std::endl;*/
1112 
1113  prim_nodal(n,0) = interpolated_u[0];
1114  prim_nodal(n,1) = cast_element_pt->pressure(interpolated_u);
1115  prim_nodal(n,2) = interpolated_u[2]/interpolated_u[0];
1116  prim_nodal(n,3) = interpolated_u[3]/interpolated_u[0];
1117  }
1118 
1119  //Neighbours centre (use value from first node
1120  for(unsigned i=0;i<n_flux;i++)
1121  {interpolated_u[i] = cell_average(0,i);}
1122 
1123  /*std::cout << "Neighbour centre : ";
1124  for(unsigned i=0;i<n_flux;i++)
1125  {std::cout << interpolated_u[i] << " ";}
1126  std::cout << std::endl;*/
1127 
1128  prim_centre(1,0) = interpolated_u[0];
1129  //Should really use neighbour here
1130  prim_centre(1,1) = cast_element_pt->pressure(interpolated_u);
1131  prim_centre(1,2) = interpolated_u[2]/interpolated_u[0];
1132  prim_centre(1,3) = interpolated_u[3]/interpolated_u[0];
1133 
1134 
1135  //We can now estimate the gradients
1136 
1137 
1138  //Output first
1139  /*std::cout << "Centre " << x[0] << " " << x[1] << " : ";
1140  for(unsigned i=0;i<n_flux;i++)
1141  {std::cout << prim_centre(0,i) << " ";}
1142 
1143  std::cout << "\n Neighbour " << centre_x(0,0) << " "
1144  << centre_x(0,1) << " : ";
1145  for(unsigned i=0;i<n_flux;i++)
1146  {std::cout << prim_centre(1,i) << " ";}
1147 
1148  for(unsigned n=0;n<2;n++)
1149  {
1150  std::cout << "\n Node " << nodal_x(n,0) << " "
1151  << nodal_x(n,1) << " : ";
1152  for(unsigned i=0;i<n_flux;i++)
1153  {std::cout << prim_nodal(n,i) << " ";}
1154  }
1155  std::cout << std::endl;*/
1156 
1157 
1158  double x_len = centre_x(1,0) - centre_x(0,0);
1159  double y_len = centre_x(1,1) - centre_x(0,1);
1160 
1161  double X_len = nodal_x(1,0) - nodal_x(0,0);
1162  double Y_len = nodal_x(1,1) - nodal_x(0,1);
1163 
1164  //If we have ghost cell's handle it!
1165  if((std::abs(x_len) < eps) && (std::abs(y_len) < eps))
1166  {
1167  //The distance is twice the distance from the centre to
1168  //the midpoint of the edges
1169  double X_mid = nodal_x(0,0) + 0.5*X_len;
1170  double Y_mid = nodal_x(0,1) + 0.5*Y_len;
1171 
1172  x_len = 2.0*(X_mid - centre_x(0,0));
1173  y_len = 2.0*(Y_mid - centre_x(0,1));
1174  }
1175 
1176  double A = X_len*y_len + x_len*Y_len;
1177 
1178  for(unsigned i=0;i<n_flux;i++)
1179  {
1180  grad(f,i,0) =
1181  ((prim_centre(1,i) - prim_centre(0,i))*Y_len
1182  - (prim_nodal(1,i) - prim_nodal(0,i))*y_len)/A;
1183  grad(f,i,1) =
1184  ((prim_centre(1,i) - prim_centre(0,i))*X_len
1185  - (prim_nodal(1,i) - prim_nodal(0,i))*x_len)/A;
1186  }
1187  }
1188 
1189  //Let's output the face gradients
1190  /*for(unsigned f=0;f<4;f++)
1191  {
1192  std::cout << "Neighbour " << size[f] << " : ";
1193  for(unsigned i=0;i<n_flux;i++)
1194  {
1195  std::cout << grad(f,i,0) << " " << grad(f,i,1) << " ";
1196  }
1197  std::cout << std::endl;
1198  }*/
1199 
1200  //Zero out the gradients
1201  for(unsigned i=0;i<n_flux;i++)
1202  {
1203  for(unsigned j=0;j<2;j++)
1204  {
1205  cast_element_pt->average_gradient(i,j) = 0.0;
1206  }
1207  }
1208 
1209  double sum = 0.0;
1210  for(unsigned f=0;f<4;f++)
1211  {
1212  const double A = size[f];
1213  sum += A;
1214  for(unsigned i=0;i<n_flux;i++)
1215  {
1216  for(unsigned j=0;j<2;j++)
1217  {
1218  cast_element_pt->average_gradient(i,j) += A*grad(f,i,j);
1219  }
1220  }
1221  }
1222 
1223  //Divide by the combined sum
1224  for(unsigned i=0;i<n_flux;i++)
1225  {
1226  for(unsigned j=0;j<2;j++)
1227  {
1228  cast_element_pt->average_gradient(i,j) /= sum;
1229  }
1230  }
1231  }
1232 
1233  //Let's have a look at it
1234  /*for(unsigned e=0;e<n_element;e++)
1235  {
1236  std::cout << "Approximate gradient in element " << e << " : ";
1237  for(unsigned i=0;i<n_flux;i++)
1238  {
1239  for(unsigned j=0;j<2;j++)
1240  {
1241  std::cout << dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))->
1242  average_gradient(i,j) << " " ;
1243  }
1244  }
1245  std::cout << "\n";
1246  }*/
1247 
1248  //OK now we use the weighted gradients to construct our limited gradient
1249  const double eps2 = 1.0e-10;
1250 
1251  for(unsigned e=0;e<n_element;e++)
1252  {
1253  ELEMENT* cast_element_pt =
1254  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
1255 
1256  /*std::cout << "Base element \n";
1257  cast_element_pt->output(std::cout);*/
1258 
1259  DenseMatrix<double> g(n_flux,4,0.0);
1260 
1261  //Get the neighbours
1262  for(unsigned f=0;f<4;f++)
1263  {
1264  ELEMENT* neighbour_bulk_element_pt = bulk_neighbour(e,f);
1265 
1266  /*std::cout << "Neighbour face" << f << "\n";
1267  neighbour_bulk_element_pt->output(std::cout);*/
1268 
1269 
1270  //Now sort it out
1271  for(unsigned i=0;i<n_flux;i++)
1272  {
1273  for(unsigned j=0;j<2;j++)
1274  {
1275  g(i,f) +=
1276  std::pow(neighbour_bulk_element_pt->average_gradient(i,j),2.0);
1277  }
1278  }
1279  }
1280 
1281  DenseMatrix<double> limited_gradient(n_flux,2,0.0);
1282 
1283  //Let's now see it
1284  for(unsigned i=0;i<n_flux;i++)
1285  {
1286 
1287  //Calcaulte the denominators
1288  double denom = 4.0*eps2;
1289  for(unsigned f=0;f<4;f++) {denom += g(i,f)*g(i,f)*g(i,f);}
1290 
1291  //Calculate the weights
1292  Vector<double> w(4,0.0);
1293  w[0] = (g(i,1)*g(i,2)*g(i,3) + eps2)/denom;
1294  w[1] = (g(i,0)*g(i,2)*g(i,3) + eps2)/denom;
1295  w[2] = (g(i,0)*g(i,1)*g(i,3) + eps2)/denom;
1296  w[3] = (g(i,0)*g(i,1)*g(i,2) + eps2)/denom;
1297 
1298  //Now reconstruct the limited gradient
1299  for(unsigned f=0;f<4;f++)
1300  {
1301  for(unsigned j=0;j<2;j++)
1302  {
1303  limited_gradient(i,j) +=
1304  w[f]*bulk_neighbour(e,f)->average_gradient(i,j);
1305  }
1306  }
1307  }
1308 
1309  //std::cout << "In element " << e << " :\n ";
1310  //We should now have the limited gradients for all fluxe
1311  /*for(unsigned i=0;i<n_flux;i++)
1312  {
1313  std::cout << cast_element_pt->average_gradient(i,0) << " "
1314  << cast_element_pt->average_gradient(i,1) << " : ";
1315  std::cout << limited_gradient(i,0) << " "
1316  << limited_gradient(i,1) << "\n";
1317  }
1318  std::cout << "\n";*/
1319 
1320  //Now we have limited the primitive variables we have to evaluate the
1321  //values at each node
1322  const unsigned n_node = cast_element_pt->nnode();
1323  //Find the centre of the element
1324  cast_element_pt->interpolated_x(s,x);
1325 
1326 
1327  //Work out distances to the nodes
1328  DenseMatrix<double> dx(n_node,2);
1329  for(unsigned n=0;n<n_node;n++)
1330  {
1331  Node* nod_pt = cast_element_pt->node_pt(n);
1332  for(unsigned i=0;i<n_dim;i++) {dx(n,i) = nod_pt->x(i) - x[i];}
1333  }
1334 
1335  //New gradients of primitive values
1336  DenseMatrix<double> dprim(n_node,n_flux);
1337 
1338 
1339  //Reconstruct limited primitive value gradients
1340  for(unsigned n=0;n<n_node;n++)
1341  {
1342  for(unsigned i=0;i<n_flux;i++)
1343  {
1344  dprim(n,i) = dx(n,0)*limited_gradient(i,0)
1345  + dx(n,1)*limited_gradient(i,1);
1346  }
1347  }
1348 
1349  //Limited values
1350  DenseMatrix<double> limited_value(n_node,n_flux);
1351  //Density
1352  for(unsigned n=0;n<n_node;n++)
1353  {
1354  limited_value(n,0) =
1355  cast_element_pt->average_prim_value(0) + dprim(n,0);
1356  }
1357 
1358  //If any are too small correct
1359  double min = 1.0e-2;
1360  bool loop_flag=false;
1361  do
1362  {
1363  bool limit_more=false;
1364  for(unsigned n=0;n<n_node;n++)
1365  {
1366  //If we're too small, reduce the gradient
1367  if(limited_value(n,0) < min)
1368  {
1369  limit_more=true;
1370  break;
1371  }
1372  }
1373 
1374  //Now do we limit more
1375  if(limit_more)
1376  {
1377  std::cout << "Limiting density gradient further \n";
1378  limited_gradient(0,0) *= 0.5;
1379  limited_gradient(0,1) *= 0.5;
1380  //Calculated the densities again
1381  for(unsigned n=0;n<n_node;n++)
1382  {
1383  limited_value(n,0) = cast_element_pt->average_prim_value(0)
1384  + dx(n,0)*limited_gradient(0,0)
1385  + dx(n,1)*limited_gradient(0,1);
1386  }
1387  //Loop again
1388  loop_flag = true;
1389  }
1390  } while(loop_flag);
1391 
1392  //Reconstruct the momentum
1393  for(unsigned n=0;n<n_node;n++)
1394  {
1395  for(unsigned i=0;i<2;i++)
1396  {
1397  limited_value(n,2+i) =
1398  cast_element_pt->average_value(2+i)
1399  + cast_element_pt->average_prim_value(0)*dprim(n,2+i)
1400  + cast_element_pt->average_prim_value(2+i)*dprim(n,0);
1401  }
1402  }
1403 
1404  //Now reconstruct the energy
1405  const double gamma = cast_element_pt->gamma();
1406 
1407  for(unsigned n=0;n<n_node;n++)
1408  {
1409  limited_value(n,1) = cast_element_pt->average_value(1)
1410  + (1.0/(gamma-1.0))*dprim(n,1)
1411  + 0.5*dprim(n,0)*(cast_element_pt->average_prim_value(2)*
1412  cast_element_pt->average_prim_value(2) +
1413  cast_element_pt->average_prim_value(3)*
1414  cast_element_pt->average_prim_value(3))
1415  + cast_element_pt->average_prim_value(0)*(
1416  cast_element_pt->average_prim_value(2)*dprim(2,n) +
1417  cast_element_pt->average_prim_value(3)*dprim(3,n));
1418  }
1419 
1420  //Check for negative pressures
1421  bool negative_pressure = false;
1422  for(unsigned n=0;n<n_node;n++)
1423  {
1424  for(unsigned i=0;i<n_flux;i++)
1425  {interpolated_u[i] = limited_value(n,i);}
1426 
1427  //Get the pressure
1428  double p = cast_element_pt->pressure(interpolated_u);
1429  if(p < eps2) {negative_pressure = true; break;}
1430  }
1431 
1432  //Correct negative pressure by limiting energy gradient to zero
1433  if(negative_pressure)
1434  {
1435  std::cout << "Correcting negative pressure\n";
1436  //If the average value is zero, then set the average value from
1437  //the average velocities
1438  double average_E = cast_element_pt->average_value(1);
1439  double rho_u = cast_element_pt->average_value(2);
1440  double rho_v = cast_element_pt->average_value(3);
1441  double rho = cast_element_pt->average_value(0);
1442 
1443  //Find the average kinetic energy
1444  double kin = 0.5*(rho_u*rho_u + rho_v*rho_v)/rho;
1445 
1446  //If the average energy is less than the average kinetic energy
1447  //we will have negative pressures (we don't want this!)
1448  if(average_E < kin)
1449  {
1450  if(average_E < 0)
1451  {
1452  throw OomphLibError("Real problem energy negative\n",
1453  OOMPH_CURRENT_FUNCTION,
1454  OOMPH_EXCEPTION_LOCATION);
1455  }
1456 
1457  //Keep the ratio of the momentum transport the same
1458  double ratio = rho_u/rho_v;
1459 
1460  //Find the sign of rho_v
1461  int sign = 1;
1462  if(rho_v < 0.0) {sign *= -1;}
1463 
1464  double sum_square_momenta = 2.0*average_E*rho;
1465 
1466  //Reset the values preserving the sign of v
1467  rho_v = sign*sqrt(sum_square_momenta/(1.0 + ratio*ratio));
1468  rho_u = ratio*rho_v;
1469  }
1470 
1471  for(unsigned n=0;n<n_node;n++)
1472  {
1473  //Set consistent averages
1474  limited_value(n,0) = rho;
1475  limited_value(n,1) = average_E;
1476  limited_value(n,2) = rho_u;
1477  limited_value(n,3) = rho_v;
1478  }
1479  }
1480 
1481  //Will have to think about boundary conditions
1482 
1483  //Finally set the limited values
1484  for(unsigned n=0;n<n_node;n++)
1485  {
1486  //Cache the node
1487  Node* nod_pt = cast_element_pt->node_pt(n);
1488  for(unsigned i=0;i<n_flux;i++)
1489  {
1490  //Only limit if it's not pinned!
1491  if(!nod_pt->is_pinned(i))
1492  {
1493  /*std::cout << "Limiting " <<
1494  nod_pt->value(i) << " to " <<
1495  limited_value(n,i) << "\n";*/
1496  //Do it
1497  nod_pt->set_value(i,limited_value(n,i));
1498  }
1499  }
1500  }
1501  }
1502 
1503  }
abs
AnnoyingScalar abs(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:135
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
p
float * p
Definition: Tutorial_Map_using.cpp:9
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
A
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:47
oomph::DGFaceElement
Definition: dg_elements.h:50
oomph::DGFaceElement::interpolated_u
virtual void interpolated_u(const Vector< double > &s, Vector< double > &f)
Return the interpolated values of the unknown fluxes.
Definition: dg_elements.cc:196
oomph::Data::is_pinned
bool is_pinned(const unsigned &i) const
Test whether the i-th variable is pinned (1: true; 0: false).
Definition: nodes.h:417
oomph::DenseMatrix
Definition: matrices.h:386
oomph::FaceElement
Definition: elements.h:4338
oomph::FaceElement::face_index
int & face_index()
Definition: elements.h:4626
oomph::FaceElement::bulk_element_pt
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4735
oomph::FaceElement::interpolated_x
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Definition: elements.h:4528
oomph::FaceElement::get_local_coordinate_in_bulk
void get_local_coordinate_in_bulk(const Vector< double > &s, Vector< double > &s_bulk) const
Definition: elements.cc:6384
oomph::FiniteElement::local_coordinate_of_node
virtual void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Definition: elements.h:1842
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
oomph::OomphLibError
Definition: oomph_definitions.h:222
oomph::RankThreeTensor
A Rank 3 Tensor class.
Definition: matrices.h:1370
f
static int f(const TensorMap< Tensor< int, 3 > > &tensor)
Definition: cxx11_tensor_map.cpp:237
min
#define min(a, b)
Definition: datatypes.h:22
s
RealScalar s
Definition: level1_cplx_impl.h:130
CRBond_Bessel::eps
double eps
Definition: crbond_bessel.cc:24
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
SYCL::sign
T sign(T x)
Definition: cxx11_tensor_builtins_sycl.cpp:172
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References abs(), oomph::FaceElement::bulk_element_pt(), e(), CRBond_Bessel::eps, f(), oomph::FaceElement::face_index(), mathsFunc::gamma(), oomph::FaceElement::get_local_coordinate_in_bulk(), i, oomph::DGFaceElement::interpolated_u(), oomph::FaceElement::interpolated_x(), oomph::Data::is_pinned(), j, oomph::FiniteElement::local_coordinate_of_node(), min, n, oomph::FiniteElement::nnode(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, p, Eigen::bfloat16_impl::pow(), s, oomph::Data::set_value(), SYCL::sign(), size, sqrt(), w, plotDoE::x, and oomph::Node::x().

◆ limit() [2/2]

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::limit ( )
inline
549  {
550  const double eps = 1.0e-14;
551  //First stage is to calculate the averages of all elements
552  const unsigned n_element = this->mesh_pt()->nelement();
553  for(unsigned e=0;e<n_element;e++)
554  {
555  //Calculate the averages
556  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))
557  ->allocate_memory_for_averages();
558  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))->calculate_averages();
559  }
560 
561  Vector<double> s(2,0.0), x(2);
562 
563  const unsigned n_flux = 4;
564 
565  const unsigned n_dim = 2;
566 
567 
568  //Store pointers to neighbours
569  DenseMatrix<ELEMENT*> bulk_neighbour(n_element,4);
570  //Storage for the unknowns
571  Vector<double> interpolated_u(n_flux);
572 
573  //Now we need to actually do the limiting, but let's check things first
574  for(unsigned e=0;e<n_element;e++)
575  {
576  //Store the element
577  ELEMENT* cast_element_pt =
578  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
579  //Find the centre of the element
580  cast_element_pt->interpolated_x(s,x);
581 
582 
583  //Storage for the approximation to the central gradient from all four
584  //faces
585  RankThreeTensor<double> grad(4,n_flux,2);
586  //Storage for the neighbours size
587  Vector<double> size(4);
588 
589  //Get the local coordinate of the nodes
590  Vector<double> s_face(1), s_bulk(2);
591  //Get the central value
592  DenseMatrix<double> prim_centre(2,n_flux);
593  //My Centre
594  for(unsigned i=0;i<n_flux;i++)
595  {interpolated_u[i] = cast_element_pt->average_value(i);}
596  prim_centre(0,0) = interpolated_u[0];
597  prim_centre(0,1) = cast_element_pt->pressure(interpolated_u);
598  prim_centre(0,2) = interpolated_u[2]/interpolated_u[0];
599  prim_centre(0,3) = interpolated_u[3]/interpolated_u[0];
600 
601  //Now store the average primitivee values
602  for(unsigned i=0;i<n_flux;i++)
603  {cast_element_pt->average_prim_value(i) = prim_centre(0,i);}
604 
605  //Loop over all the faces
606  for(unsigned f=0;f<4;f++)
607  {
608  DGFaceElement* face_element_pt =
609  dynamic_cast<DGFaceElement*>(cast_element_pt->face_element_pt(f));
610 
611  //Get the average of the nodal points on the face
612  //(***NOT GENERAL**)
613  const unsigned n_face_node = face_element_pt->nnode();
614 
615  //Storage for the average nodal values
616  DenseMatrix<double> face_average(n_face_node,n_flux);
617 
618  //Storage for the neighbouring cell's centre average
619  DenseMatrix<double> cell_average(n_face_node,n_flux);
620 
621  //Storage for the sizes of the neighbouring bulk elements
622  Vector<double> neighbour_size(n_face_node);
623 
624  //Now the storage for the nodal values and centre values
625  DenseMatrix<double> nodal_x(n_face_node,n_dim);
626  //Now the storage for the centre values
627  DenseMatrix<double> centre_x(n_face_node,n_dim);
628 
629  //Neighbours face
630  FaceElement* neighbour_face_pt=0;
631 
632  //Now calculate those average values
633  for(unsigned n=0;n<n_face_node;n++)
634  {
635  //Get the local coordinates of the node
636  face_element_pt->local_coordinate_of_node(n,s_face);
637  //Get the global coordinate of the node
638  for(unsigned i=0;i<n_dim;i++)
639  {
640  nodal_x(n,i) = face_element_pt->interpolated_x(s_face,i);
641  }
642 
643  //Get the flux at the node in this face
644  face_element_pt->interpolated_u(s_face,interpolated_u);
645 
646  /*std::cout << "Values at face :";
647  for(unsigned i=0;i<n_flux;i++)
648  {
649  std::cout << interpolated_u[i] << " ";
650  }
651  std::cout << std::endl;*/
652 
653  for(unsigned i=0;i<n_flux;i++)
654  {
655  face_average(n,i) = interpolated_u[i];
656  }
657 
658  //Now get the bulk coordinate
659  face_element_pt->get_local_coordinate_in_bulk(s_face,s_bulk);
660  //Now get the neighbour
661  cast_element_pt->
662  get_neighbouring_face_and_local_coordinate(
663  face_element_pt->face_index(),s_bulk,neighbour_face_pt,s_face);
664 
665  //Get the fluxes of the neighbour
666  dynamic_cast<DGFaceElement*>(neighbour_face_pt)
667  ->interpolated_u(s_face,interpolated_u);
668  //Add the flux from the neighbour
669  ELEMENT* neighbour_bulk_element_pt = dynamic_cast<ELEMENT*>(
670  neighbour_face_pt->bulk_element_pt());
671 
672  /*std::cout << "Values at face neighbour :";
673  for(unsigned i=0;i<n_flux;i++)
674  {
675  std::cout << interpolated_u[i] << " ";
676  }
677  std::cout << std::endl;
678 
679  std::cout << "SANITY CHECK " << face_element_pt << " "
680  << neighbour_face_pt << " "
681  << cast_element_pt << " "
682  << neighbour_bulk_element_pt << "\n";*/
683 
684 
685  //Get the centre of the neighbouring bulk element
686  for(unsigned i=0;i<n_dim;i++) {s_bulk[i] = 0.0;}
687  //Do it
688  for(unsigned i=0;i<n_dim;i++)
689  {
690  centre_x(n,i) = neighbour_bulk_element_pt
691  ->interpolated_x(s_bulk,i);
692  }
693 
694  for(unsigned i=0;i<n_flux;i++)
695  {
696  face_average(n,i) += interpolated_u[i];
697  face_average(n,i) *= 0.5;
698  //Get the neighbour
699  cell_average(n,i) = neighbour_bulk_element_pt->average_value(i);
700  }
701 
702  //Get the size of the neighbour
703  neighbour_size[n] = neighbour_bulk_element_pt->average_value(n_flux);
704 
705  //Only bother to store the first one
706  if(n==0)
707  {
708  bulk_neighbour(e,f) = neighbour_bulk_element_pt;
709  }
710  }
711 
712  //Now we have the face averages and the cell averages
713  //Check that we have a conforming mesh
714  for(unsigned n=1;n<n_face_node;n++)
715  {
716  if((std::abs(centre_x(0,0) - centre_x(n,0)) > eps) ||
717  (std::abs(centre_x(0,1) - centre_x(n,1)) > eps))
718  {
719  throw OomphLibError("Mesh is not conforming in limiter",
720  OOMPH_CURRENT_FUNCTION,
721  OOMPH_EXCEPTION_LOCATION);
722  }
723  }
724 
725 
726  //Set centre_x(0,0) to be the current element
727  for(unsigned i=0;i<n_dim;i++)
728  {
729  centre_x(0,i) = x[i];
730  }
731 
732  //Set the size of the neighbour
733  size[f] = neighbour_size[0];
734 
735  //Convert our values to primitive variables
736  DenseMatrix<double> prim_nodal(2,n_flux);
737 
738  //Face averages (corner nodes)
739  //std::cout << "Nodal values \n";
740  for(unsigned n=0;n<2;n++)
741  {
742  //Load the fluxes
743  for(unsigned i=0;i<n_flux;i++)
744  {interpolated_u[i] = face_average(n,i);}
745 
746  //Output the fluxes
747  /*std::cout << n << " : ";
748  for(unsigned i=0;i<n_flux;i++)
749  {std::cout << interpolated_u[i] << " ";}
750  std::cout << std::endl;*/
751 
752  prim_nodal(n,0) = interpolated_u[0];
753  prim_nodal(n,1) = cast_element_pt->pressure(interpolated_u);
754  prim_nodal(n,2) = interpolated_u[2]/interpolated_u[0];
755  prim_nodal(n,3) = interpolated_u[3]/interpolated_u[0];
756  }
757 
758  //Neighbours centre (use value from first node
759  for(unsigned i=0;i<n_flux;i++)
760  {interpolated_u[i] = cell_average(0,i);}
761 
762  /*std::cout << "Neighbour centre : ";
763  for(unsigned i=0;i<n_flux;i++)
764  {std::cout << interpolated_u[i] << " ";}
765  std::cout << std::endl;*/
766 
767  prim_centre(1,0) = interpolated_u[0];
768  //Should really use neighbour here
769  prim_centre(1,1) = cast_element_pt->pressure(interpolated_u);
770  prim_centre(1,2) = interpolated_u[2]/interpolated_u[0];
771  prim_centre(1,3) = interpolated_u[3]/interpolated_u[0];
772 
773 
774  //We can now estimate the gradients
775 
776 
777  //Output first
778  /*std::cout << "Centre " << x[0] << " " << x[1] << " : ";
779  for(unsigned i=0;i<n_flux;i++)
780  {std::cout << prim_centre(0,i) << " ";}
781 
782  std::cout << "\n Neighbour " << centre_x(0,0) << " "
783  << centre_x(0,1) << " : ";
784  for(unsigned i=0;i<n_flux;i++)
785  {std::cout << prim_centre(1,i) << " ";}
786 
787  for(unsigned n=0;n<2;n++)
788  {
789  std::cout << "\n Node " << nodal_x(n,0) << " "
790  << nodal_x(n,1) << " : ";
791  for(unsigned i=0;i<n_flux;i++)
792  {std::cout << prim_nodal(n,i) << " ";}
793  }
794  std::cout << std::endl;*/
795 
796 
797  double x_len = centre_x(1,0) - centre_x(0,0);
798  double y_len = centre_x(1,1) - centre_x(0,1);
799 
800  double X_len = nodal_x(1,0) - nodal_x(0,0);
801  double Y_len = nodal_x(1,1) - nodal_x(0,1);
802 
803  //If we have ghost cell's handle it!
804  if((std::abs(x_len) < eps) && (std::abs(y_len) < eps))
805  {
806  //The distance is twice the distance from the centre to
807  //the midpoint of the edges
808  double X_mid = nodal_x(0,0) + 0.5*X_len;
809  double Y_mid = nodal_x(0,1) + 0.5*Y_len;
810 
811  x_len = 2.0*(X_mid - centre_x(0,0));
812  y_len = 2.0*(Y_mid - centre_x(0,1));
813  }
814 
815  double A = X_len*y_len + x_len*Y_len;
816 
817  for(unsigned i=0;i<n_flux;i++)
818  {
819  grad(f,i,0) =
820  ((prim_centre(1,i) - prim_centre(0,i))*Y_len
821  - (prim_nodal(1,i) - prim_nodal(0,i))*y_len)/A;
822  grad(f,i,1) =
823  ((prim_centre(1,i) - prim_centre(0,i))*X_len
824  - (prim_nodal(1,i) - prim_nodal(0,i))*x_len)/A;
825  }
826  }
827 
828  //Let's output the face gradients
829  /*for(unsigned f=0;f<4;f++)
830  {
831  std::cout << "Neighbour " << size[f] << " : ";
832  for(unsigned i=0;i<n_flux;i++)
833  {
834  std::cout << grad(f,i,0) << " " << grad(f,i,1) << " ";
835  }
836  std::cout << std::endl;
837  }*/
838 
839  //Zero out the gradients
840  for(unsigned i=0;i<n_flux;i++)
841  {
842  for(unsigned j=0;j<2;j++)
843  {
844  cast_element_pt->average_gradient(i,j) = 0.0;
845  }
846  }
847 
848  double sum = 0.0;
849  for(unsigned f=0;f<4;f++)
850  {
851  const double A = size[f];
852  sum += A;
853  for(unsigned i=0;i<n_flux;i++)
854  {
855  for(unsigned j=0;j<2;j++)
856  {
857  cast_element_pt->average_gradient(i,j) += A*grad(f,i,j);
858  }
859  }
860  }
861 
862  //Divide by the combined sum
863  for(unsigned i=0;i<n_flux;i++)
864  {
865  for(unsigned j=0;j<2;j++)
866  {
867  cast_element_pt->average_gradient(i,j) /= sum;
868  }
869  }
870  }
871 
872  //Let's have a look at it
873  /*for(unsigned e=0;e<n_element;e++)
874  {
875  std::cout << "Approximate gradient in element " << e << " : ";
876  for(unsigned i=0;i<n_flux;i++)
877  {
878  for(unsigned j=0;j<2;j++)
879  {
880  std::cout << dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e))->
881  average_gradient(i,j) << " " ;
882  }
883  }
884  std::cout << "\n";
885  }*/
886 
887  //OK now we use the weighted gradients to construct our limited gradient
888  const double eps2 = 1.0e-10;
889 
890  for(unsigned e=0;e<n_element;e++)
891  {
892  ELEMENT* cast_element_pt =
893  dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
894 
895  /*std::cout << "Base element \n";
896  cast_element_pt->output(std::cout);*/
897 
898  DenseMatrix<double> g(n_flux,4,0.0);
899 
900  //Get the neighbours
901  for(unsigned f=0;f<4;f++)
902  {
903  ELEMENT* neighbour_bulk_element_pt = bulk_neighbour(e,f);
904 
905  /*std::cout << "Neighbour face" << f << "\n";
906  neighbour_bulk_element_pt->output(std::cout);*/
907 
908 
909  //Now sort it out
910  for(unsigned i=0;i<n_flux;i++)
911  {
912  for(unsigned j=0;j<2;j++)
913  {
914  g(i,f) +=
915  std::pow(neighbour_bulk_element_pt->average_gradient(i,j),2.0);
916  }
917  }
918  }
919 
920  DenseMatrix<double> limited_gradient(n_flux,2,0.0);
921 
922  //Let's now see it
923  for(unsigned i=0;i<n_flux;i++)
924  {
925 
926  //Calcaulte the denominators
927  double denom = 4.0*eps2;
928  for(unsigned f=0;f<4;f++) {denom += g(i,f)*g(i,f)*g(i,f);}
929 
930  //Calculate the weights
931  Vector<double> w(4,0.0);
932  w[0] = (g(i,1)*g(i,2)*g(i,3) + eps2)/denom;
933  w[1] = (g(i,0)*g(i,2)*g(i,3) + eps2)/denom;
934  w[2] = (g(i,0)*g(i,1)*g(i,3) + eps2)/denom;
935  w[3] = (g(i,0)*g(i,1)*g(i,2) + eps2)/denom;
936 
937  //Now reconstruct the limited gradient
938  for(unsigned f=0;f<4;f++)
939  {
940  for(unsigned j=0;j<2;j++)
941  {
942  limited_gradient(i,j) +=
943  w[f]*bulk_neighbour(e,f)->average_gradient(i,j);
944  }
945  }
946  }
947 
948  //std::cout << "In element " << e << " :\n ";
949  //We should now have the limited gradients for all fluxe
950  /*for(unsigned i=0;i<n_flux;i++)
951  {
952  std::cout << cast_element_pt->average_gradient(i,0) << " "
953  << cast_element_pt->average_gradient(i,1) << " : ";
954  std::cout << limited_gradient(i,0) << " "
955  << limited_gradient(i,1) << "\n";
956  }
957  std::cout << "\n";*/
958 
959  //Now we have limited the primitive variables we have to evaluate the
960  //values at each node
961  const unsigned n_node = cast_element_pt->nnode();
962  //Find the centre of the element
963  cast_element_pt->interpolated_x(s,x);
964 
965 
966  //Work out distances to the nodes
967  DenseMatrix<double> dx(n_node,2);
968  for(unsigned n=0;n<n_node;n++)
969  {
970  Node* nod_pt = cast_element_pt->node_pt(n);
971  for(unsigned i=0;i<n_dim;i++) {dx(n,i) = nod_pt->x(i) - x[i];}
972  }
973 
974  //New gradients of primitive values
975  DenseMatrix<double> dprim(n_node,n_flux);
976 
977 
978  //Reconstruct limited primitive value gradients
979  for(unsigned n=0;n<n_node;n++)
980  {
981  for(unsigned i=0;i<n_flux;i++)
982  {
983  dprim(n,i) = dx(n,0)*limited_gradient(i,0)
984  + dx(n,1)*limited_gradient(i,1);
985  }
986  }
987 
988  //Limited values
989  DenseMatrix<double> limited_value(n_node,n_flux);
990  //Density
991  for(unsigned n=0;n<n_node;n++)
992  {
993  limited_value(n,0) =
994  cast_element_pt->average_prim_value(0) + dprim(n,0);
995  }
996 
997  //If any are too small correct
998  double min = 1.0e-2;
999  bool loop_flag=false;
1000  do
1001  {
1002  bool limit_more=false;
1003  for(unsigned n=0;n<n_node;n++)
1004  {
1005  //If we're too small, reduce the gradient
1006  if(limited_value(n,0) < min)
1007  {
1008  limit_more=true;
1009  break;
1010  }
1011  }
1012 
1013  //Now do we limit more
1014  if(limit_more)
1015  {
1016  std::cout << "Limiting density gradient further \n";
1017  limited_gradient(0,0) *= 0.5;
1018  limited_gradient(0,1) *= 0.5;
1019  //Calculated the densities again
1020  for(unsigned n=0;n<n_node;n++)
1021  {
1022  limited_value(n,0) = cast_element_pt->average_prim_value(0)
1023  + dx(n,0)*limited_gradient(0,0)
1024  + dx(n,1)*limited_gradient(0,1);
1025  }
1026  //Loop again
1027  loop_flag = true;
1028  }
1029  } while(loop_flag);
1030 
1031  //Reconstruct the momentum
1032  for(unsigned n=0;n<n_node;n++)
1033  {
1034  for(unsigned i=0;i<2;i++)
1035  {
1036  limited_value(n,2+i) =
1037  cast_element_pt->average_value(2+i)
1038  + cast_element_pt->average_prim_value(0)*dprim(n,2+i)
1039  + cast_element_pt->average_prim_value(2+i)*dprim(n,0);
1040  }
1041  }
1042 
1043  //Now reconstruct the energy
1044  const double gamma = cast_element_pt->gamma();
1045 
1046  for(unsigned n=0;n<n_node;n++)
1047  {
1048  limited_value(n,1) = cast_element_pt->average_value(1)
1049  + (1.0/(gamma-1.0))*dprim(n,1)
1050  + 0.5*dprim(n,0)*(cast_element_pt->average_prim_value(2)*
1051  cast_element_pt->average_prim_value(2) +
1052  cast_element_pt->average_prim_value(3)*
1053  cast_element_pt->average_prim_value(3))
1054  + cast_element_pt->average_prim_value(0)*(
1055  cast_element_pt->average_prim_value(2)*dprim(2,n) +
1056  cast_element_pt->average_prim_value(3)*dprim(3,n));
1057  }
1058 
1059  //Check for negative pressures
1060  bool negative_pressure = false;
1061  for(unsigned n=0;n<n_node;n++)
1062  {
1063  for(unsigned i=0;i<n_flux;i++)
1064  {interpolated_u[i] = limited_value(n,i);}
1065 
1066  //Get the pressure
1067  double p = cast_element_pt->pressure(interpolated_u);
1068  if(p < eps2) {negative_pressure = true; break;}
1069  }
1070 
1071  //Correct negative pressure by limiting energy gradient to zero
1072  if(negative_pressure)
1073  {
1074  std::cout << "Correcting negative pressure\n";
1075  //If the average value is zero, then set the average value from
1076  //the average velocities
1077  double average_E = cast_element_pt->average_value(1);
1078  double rho_u = cast_element_pt->average_value(2);
1079  double rho_v = cast_element_pt->average_value(3);
1080  double rho = cast_element_pt->average_value(0);
1081 
1082  //Find the average kinetic energy
1083  double kin = 0.5*(rho_u*rho_u + rho_v*rho_v)/rho;
1084 
1085  //If the average energy is less than the average kinetic energy
1086  //we will have negative pressures (we don't want this!)
1087  if(average_E < kin)
1088  {
1089  if(average_E < 0)
1090  {
1091  throw OomphLibError("Real problem energy negative\n",
1092  OOMPH_CURRENT_FUNCTION,
1093  OOMPH_EXCEPTION_LOCATION);
1094  }
1095 
1096  //Keep the ratio of the momentum transport the same
1097  double ratio = rho_u/rho_v;
1098 
1099  //Find the sign of rho_v
1100  int sign = 1;
1101  if(rho_v < 0.0) {sign *= -1;}
1102 
1103  double sum_square_momenta = 2.0*average_E*rho;
1104 
1105  //Reset the values preserving the sign of v
1106  rho_v = sign*sqrt(sum_square_momenta/(1.0 + ratio*ratio));
1107  rho_u = ratio*rho_v;
1108  }
1109 
1110  for(unsigned n=0;n<n_node;n++)
1111  {
1112  //Set consistent averages
1113  limited_value(n,0) = rho;
1114  limited_value(n,1) = average_E;
1115  limited_value(n,2) = rho_u;
1116  limited_value(n,3) = rho_v;
1117  }
1118  }
1119 
1120  //Will have to think about boundary conditions
1121 
1122  //Finally set the limited values
1123  for(unsigned n=0;n<n_node;n++)
1124  {
1125  //Cache the node
1126  Node* nod_pt = cast_element_pt->node_pt(n);
1127  for(unsigned i=0;i<n_flux;i++)
1128  {
1129  //Only limit if it's not pinned!
1130  if(!nod_pt->is_pinned(i))
1131  {
1132  /*std::cout << "Limiting " <<
1133  nod_pt->value(i) << " to " <<
1134  limited_value(n,i) << "\n";*/
1135  //Do it
1136  nod_pt->set_value(i,limited_value(n,i));
1137  }
1138  }
1139  }
1140  }
1141 
1142  }

References abs(), oomph::FaceElement::bulk_element_pt(), e(), CRBond_Bessel::eps, f(), oomph::FaceElement::face_index(), mathsFunc::gamma(), oomph::FaceElement::get_local_coordinate_in_bulk(), i, oomph::DGFaceElement::interpolated_u(), oomph::FaceElement::interpolated_x(), oomph::Data::is_pinned(), j, oomph::FiniteElement::local_coordinate_of_node(), min, n, oomph::FiniteElement::nnode(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, p, Eigen::bfloat16_impl::pow(), s, oomph::Data::set_value(), SYCL::sign(), size, sqrt(), w, plotDoE::x, and oomph::Node::x().

◆ mesh_pt() [1/5]

template<class ELEMENT >
TwoDDGMesh<ELEMENT>* TwoDDGProblem< ELEMENT >::mesh_pt ( )
inline
240  {return dynamic_cast<TwoDDGMesh<ELEMENT>*>(Problem::mesh_pt());}

References oomph::Problem::mesh_pt().

◆ mesh_pt() [2/5]

template<class ELEMENT >
TwoDDGMesh<ELEMENT>* TwoDDGProblem< ELEMENT >::mesh_pt ( )
inline
377  {return dynamic_cast<TwoDDGMesh<ELEMENT>*>(Problem::mesh_pt());}

References oomph::Problem::mesh_pt().

◆ mesh_pt() [3/5]

template<class ELEMENT >
TwoDDGMesh<ELEMENT>* TwoDDGProblem< ELEMENT >::mesh_pt ( )
inline
845  {return dynamic_cast<TwoDDGMesh<ELEMENT>*>(Problem::mesh_pt());}

References oomph::Problem::mesh_pt().

◆ mesh_pt() [4/5]

template<class ELEMENT >
TwoDDGMesh<ELEMENT>* TwoDDGProblem< ELEMENT >::mesh_pt ( )
inline
482  {return dynamic_cast<TwoDDGMesh<ELEMENT>*>(Problem::mesh_pt());}

References oomph::Problem::mesh_pt().

◆ mesh_pt() [5/5]

template<class ELEMENT >
TwoDDGMesh<ELEMENT>* TwoDDGProblem< ELEMENT >::mesh_pt ( )
inline
406  {return dynamic_cast<TwoDDGMesh<ELEMENT>*>(Problem::mesh_pt());}

References oomph::Problem::mesh_pt().

◆ parameter_study() [1/5]

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::parameter_study ( std::ostream &  trace,
const bool &  disc 
)
inline
307  {
308  this->enable_mass_matrix_reuse();
309  double dt = 0.001;
310  apply_initial_conditions();
311 
312  Vector<double> error(1,0.0);
313  char filename[100];
314 
315  unsigned count=1;
316  for(unsigned t=0;t<500;t++)
317  {
318  explicit_timestep(dt);
319  if(count==250)
320  {
321  if(disc)
322  {
323  sprintf(filename,"disc_%li_time%g.dat",mesh_pt()->nelement(),time());
324  }
325  else
326  {
327  sprintf(filename,"cont_%li_time%g.dat",mesh_pt()->nelement(),time());
328  }
329  std::ofstream outfile(filename);
330  mesh_pt()->output(outfile);
331  count=0;
332  }
333 
334  compute_error(this->time(),error);
335  trace << mesh_pt()->nelement() << " "
336  << this->time() << " " << sqrt(error[0]) << std::endl;
337  ++count;
338  }
339  }
TwoDDGProblem::apply_initial_conditions
void apply_initial_conditions()
Definition: two_d_advection.cc:287
oomph::Problem::explicit_timestep
void explicit_timestep(const double &dt, const bool &shift_values=true)
Take an explicit timestep of size dt.
Definition: problem.cc:10918
oomph::Problem::enable_mass_matrix_reuse
void enable_mass_matrix_reuse()
Definition: problem.cc:11807
MergeRestartFiles.filename
string filename
Definition: MergeRestartFiles.py:39

References oomph::compute_error(), calibrate::error, MergeRestartFiles::filename, sqrt(), and plotPSD::t.

◆ parameter_study() [2/5]

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::parameter_study ( std::ostream &  trace,
const bool &  disc 
)
inline
477  {
478  this->enable_mass_matrix_reuse();
479  double dt = 0.001;
480  apply_initial_conditions();
481 
482  Vector<double> error(4,0.0);
483  char filename[100];
484 
485  unsigned count=1;
486  for(unsigned t=0;t<200;t++)
487  {
488  explicit_timestep(dt);
489  if(count==50)
490  {
491  if(disc)
492  {
493  sprintf(filename,"disc_%li_time%g.dat",mesh_pt()->nelement(),time());
494  }
495  else
496  {
497  sprintf(filename,"cont_%li_time%g.dat",mesh_pt()->nelement(),time());
498  }
499  std::ofstream outfile(filename);
500  mesh_pt()->output(outfile);
501  count=0;
502  }
503 
504  compute_error(this->time(),error);
505  trace << mesh_pt()->nelement() << " "
506  << this->time() << " " << sqrt(error[0]) <<
507  " " << sqrt(error[1]) << " " << sqrt(error[2])
508  << " " << sqrt(error[3]) << std::endl;
509  ++count;
510  }
511  }

References oomph::compute_error(), calibrate::error, MergeRestartFiles::filename, sqrt(), and plotPSD::t.

◆ parameter_study() [3/5]

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::parameter_study ( std::ostream &  trace,
const bool &  disc 
)
inline
1508  {
1509  this->enable_mass_matrix_reuse();
1510  double dt = 0.001;
1511  apply_initial_conditions();
1512 
1513  Vector<double> error(4,0.0);
1514  char filename[100];
1515 
1516  unsigned count=1;
1517  for(unsigned t=0;t<5000;t++)
1518  {
1519  std::cout << "Taking timestep " << t << "\n";
1520  explicit_timestep(dt);
1521  //Now limit
1522  limit();
1523  if(count==50)
1524  {
1525  if(disc)
1526  {
1527  sprintf(filename,"disc_%li_time%g.dat",mesh_pt()->nelement(),time());
1528  }
1529  else
1530  {
1531  sprintf(filename,"cont_%li_time%g.dat",mesh_pt()->nelement(),time());
1532  }
1533  std::ofstream outfile(filename);
1534  mesh_pt()->output(outfile);
1535  count=0;
1536  }
1537  trace << mesh_pt()->nelement() << " "
1538  << this->time() << std::endl;
1539  ++count;
1540  }
1541  }
TwoDDGProblem::limit
void limit()
Definition: ff_step.cc:909

References calibrate::error, MergeRestartFiles::filename, and plotPSD::t.

◆ parameter_study() [4/5]

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::parameter_study ( std::ostream &  trace,
const bool &  disc 
)
inline
1147  {
1148  this->enable_mass_matrix_reuse();
1149  double dt = 0.0001;
1150  apply_initial_conditions();
1151 
1152  Vector<double> error(4,0.0);
1153  char filename[100];
1154 
1155  unsigned count=1;
1156  for(unsigned t=0;t<5000;t++)
1157  {
1158  explicit_timestep(dt);
1159  //Now limit
1160  //limit();
1161  if(count==50)
1162  {
1163  if(disc)
1164  {
1165  sprintf(filename,"disc_%li_time%g.dat",mesh_pt()->nelement(),time());
1166  }
1167  else
1168  {
1169  sprintf(filename,"cont_%li_time%g.dat",mesh_pt()->nelement(),time());
1170  }
1171  std::ofstream outfile(filename);
1172  mesh_pt()->output(outfile);
1173  count=0;
1174  }
1175 
1176  trace << mesh_pt()->nelement() << " "
1177  << this->time() << std::endl;
1178  ++count;
1179  }
1180  }

References calibrate::error, MergeRestartFiles::filename, and plotPSD::t.

◆ parameter_study() [5/5]

template<class ELEMENT >
void TwoDDGProblem< ELEMENT >::parameter_study ( std::ostream &  trace,
const bool &  disc 
)
inline
534  {
535  this->enable_mass_matrix_reuse();
536  double dt = 0.001;
537  apply_initial_conditions();
538 
539  Vector<double> error(4,0.0);
540  char filename[100];
541 
542  unsigned count=1;
543  for(unsigned t=0;t<50;t++)
544  {
545  explicit_timestep(dt);
546  if(count==50)
547  {
548  if(disc)
549  {
550  sprintf(filename,"disc_%li_time%g.dat",mesh_pt()->nelement(),time());
551  }
552  else
553  {
554  sprintf(filename,"cont_%li_time%g.dat",mesh_pt()->nelement(),time());
555  }
556  std::ofstream outfile(filename);
557  mesh_pt()->output(outfile);
558  count=0;
559  }
560 
561  compute_error(this->time(),error);
562  trace << mesh_pt()->nelement() << " "
563  << this->time() << " " << sqrt(error[0]) <<
564  " " << sqrt(error[1]) << " " << sqrt(error[2])
565  << " " << sqrt(error[3]) << std::endl;
566  ++count;
567  }
568  }

References oomph::compute_error(), calibrate::error, MergeRestartFiles::filename, sqrt(), and plotPSD::t.

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