#include <continuation_system.h>

Inheritance diagram for libMesh::ContinuationSystem:

[legend]

Public Types
enum	Predictor { Euler, AB2, Invalid_Predictor }

typedef ContinuationSystem	sys_type

typedef FEMSystem	Parent

typedef bool(TimeSolver::*	TimeSolverResPtr) (bool, DiffContext &)

typedef std::map< std::string, SparseMatrix< Number > * >::iterator	matrices_iterator

typedef std::map< std::string, SparseMatrix< Number > * >::const_iterator	const_matrices_iterator

typedef Number(*	ValueFunctionPointer) (const Point &p, const Parameters &Parameters, const std::string &sys_name, const std::string &unknown_name)

typedef Gradient(*	GradientFunctionPointer) (const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name)

typedef std::map< std::string, NumericVector< Number > * >::iterator	vectors_iterator

typedef std::map< std::string, NumericVector< Number > * >::const_iterator	const_vectors_iterator

Public Member Functions
	ContinuationSystem (EquationSystems &es, const std::string &name, const unsigned int number)

virtual	~ContinuationSystem ()

virtual void	clear () override

virtual void	solve () override

void	continuation_solve ()

void	advance_arcstep ()

void	set_max_arclength_stepsize (Real maxds)

void	save_current_solution ()

virtual void	assembly (bool get_residual, bool get_jacobian, bool apply_heterogeneous_constraints=false, bool apply_no_constraints=false) override

void	mesh_position_get ()

void	mesh_position_set ()

virtual std::unique_ptr< DiffContext >	build_context () override

virtual void	init_context (DiffContext &) override

virtual void	postprocess () override

virtual void	assemble_qoi (const QoISet &indices=QoISet()) override

virtual void	assemble_qoi_derivative (const QoISet &qoi_indices=QoISet(), bool include_liftfunc=true, bool apply_constraints=true) override

Real	numerical_jacobian_h_for_var (unsigned int var_num) const

void	set_numerical_jacobian_h_for_var (unsigned int var_num, Real new_h)

void	numerical_jacobian (TimeSolverResPtr res, FEMContext &context) const

void	numerical_elem_jacobian (FEMContext &context) const

void	numerical_side_jacobian (FEMContext &context) const

void	numerical_nonlocal_jacobian (FEMContext &context) const

virtual void	reinit () override

virtual void	assemble () override

virtual LinearSolver< Number > *	get_linear_solver () const override

virtual std::pair< unsigned int, Real >	get_linear_solve_parameters () const override

virtual void	release_linear_solver (LinearSolver< Number > *) const override

virtual std::pair< unsigned int, Real >	adjoint_solve (const QoISet &qoi_indices=QoISet()) override

virtual std::unique_ptr< DifferentiablePhysics >	clone_physics () override

virtual std::unique_ptr< DifferentiableQoI >	clone () override

const DifferentiablePhysics *	get_physics () const

DifferentiablePhysics *	get_physics ()

void	attach_physics (DifferentiablePhysics *physics_in)

void	swap_physics (DifferentiablePhysics *&swap_physics)

const DifferentiableQoI *	get_qoi () const

DifferentiableQoI *	get_qoi ()

void	attach_qoi (DifferentiableQoI *qoi_in)

void	set_time_solver (std::unique_ptr< TimeSolver > _time_solver)

TimeSolver &	get_time_solver ()

const TimeSolver &	get_time_solver () const

virtual void	element_postprocess (DiffContext &)

virtual void	side_postprocess (DiffContext &)

unsigned int	get_second_order_dot_var (unsigned int var) const

bool	have_first_order_scalar_vars () const

bool	have_second_order_scalar_vars () const

sys_type &	system ()

virtual void	disable_cache () override

virtual std::string	system_type () const override

virtual void	assemble_residual_derivatives (const ParameterVector &parameters) override

virtual std::pair< unsigned int, Real >	sensitivity_solve (const ParameterVector &parameters) override

virtual std::pair< unsigned int, Real >	weighted_sensitivity_solve (const ParameterVector &parameters, const ParameterVector &weights) override

virtual std::pair< unsigned int, Real >	weighted_sensitivity_adjoint_solve (const ParameterVector &parameters, const ParameterVector &weights, const QoISet &qoi_indices=QoISet()) override

virtual void	adjoint_qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities) override

virtual void	forward_qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities) override

virtual void	qoi_parameter_hessian (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &hessian) override

virtual void	qoi_parameter_hessian_vector_product (const QoISet &qoi_indices, const ParameterVector &parameters, const ParameterVector &vector, SensitivityData &product) override

SparseMatrix< Number > &	add_matrix (const std::string &mat_name)

void	remove_matrix (const std::string &mat_name)

bool	have_matrix (const std::string &mat_name) const

const SparseMatrix< Number > *	request_matrix (const std::string &mat_name) const

SparseMatrix< Number > *	request_matrix (const std::string &mat_name)

const SparseMatrix< Number > &	get_matrix (const std::string &mat_name) const

SparseMatrix< Number > &	get_matrix (const std::string &mat_name)

virtual unsigned int	n_matrices () const override

void	init ()

virtual void	reinit_constraints ()

bool	is_initialized ()

virtual void	update ()

virtual void	restrict_solve_to (const SystemSubset *subset, const SubsetSolveMode subset_solve_mode=SUBSET_ZERO)

bool	is_adjoint_already_solved () const

void	set_adjoint_already_solved (bool setting)

virtual void	qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)

virtual bool	compare (const System &other_system, const Real threshold, const bool verbose) const

const std::string &	name () const

void	project_solution (FunctionBase< Number > f, FunctionBase< Gradient > g=nullptr) const

void	project_solution (FEMFunctionBase< Number > f, FEMFunctionBase< Gradient > g=nullptr) const

void	project_solution (ValueFunctionPointer fptr, GradientFunctionPointer gptr, const Parameters &parameters) const

void	project_vector (NumericVector< Number > &new_vector, FunctionBase< Number > f, FunctionBase< Gradient > g=nullptr, int is_adjoint=-1) const

void	project_vector (NumericVector< Number > &new_vector, FEMFunctionBase< Number > f, FEMFunctionBase< Gradient > g=nullptr, int is_adjoint=-1) const

void	project_vector (ValueFunctionPointer fptr, GradientFunctionPointer gptr, const Parameters &parameters, NumericVector< Number > &new_vector, int is_adjoint=-1) const

void	boundary_project_solution (const std::set< boundary_id_type > &b, const std::vector< unsigned int > &variables, FunctionBase< Number > f, FunctionBase< Gradient > g=nullptr)

void	boundary_project_solution (const std::set< boundary_id_type > &b, const std::vector< unsigned int > &variables, ValueFunctionPointer fptr, GradientFunctionPointer gptr, const Parameters &parameters)

void	boundary_project_vector (const std::set< boundary_id_type > &b, const std::vector< unsigned int > &variables, NumericVector< Number > &new_vector, FunctionBase< Number > f, FunctionBase< Gradient > g=nullptr, int is_adjoint=-1) const

void	boundary_project_vector (const std::set< boundary_id_type > &b, const std::vector< unsigned int > &variables, ValueFunctionPointer fptr, GradientFunctionPointer gptr, const Parameters &parameters, NumericVector< Number > &new_vector, int is_adjoint=-1) const

unsigned int	number () const

void	update_global_solution (std::vector< Number > &global_soln) const

void	update_global_solution (std::vector< Number > &global_soln, const processor_id_type dest_proc) const

const MeshBase &	get_mesh () const

MeshBase &	get_mesh ()

const DofMap &	get_dof_map () const

DofMap &	get_dof_map ()

const EquationSystems &	get_equation_systems () const

EquationSystems &	get_equation_systems ()

bool	active () const

void	activate ()

void	deactivate ()

void	set_basic_system_only ()

vectors_iterator	vectors_begin ()

const_vectors_iterator	vectors_begin () const

vectors_iterator	vectors_end ()

const_vectors_iterator	vectors_end () const

NumericVector< Number > &	add_vector (const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)

void	remove_vector (const std::string &vec_name)

bool &	project_solution_on_reinit (void)

bool	have_vector (const std::string &vec_name) const

const NumericVector< Number > *	request_vector (const std::string &vec_name) const

NumericVector< Number > *	request_vector (const std::string &vec_name)

const NumericVector< Number > *	request_vector (const unsigned int vec_num) const

NumericVector< Number > *	request_vector (const unsigned int vec_num)

const NumericVector< Number > &	get_vector (const std::string &vec_name) const

NumericVector< Number > &	get_vector (const std::string &vec_name)

const NumericVector< Number > &	get_vector (const unsigned int vec_num) const

NumericVector< Number > &	get_vector (const unsigned int vec_num)

const std::string &	vector_name (const unsigned int vec_num) const

const std::string &	vector_name (const NumericVector< Number > &vec_reference) const

void	set_vector_as_adjoint (const std::string &vec_name, int qoi_num)

int	vector_is_adjoint (const std::string &vec_name) const

void	set_vector_preservation (const std::string &vec_name, bool preserve)

bool	vector_preservation (const std::string &vec_name) const

NumericVector< Number > &	add_adjoint_solution (unsigned int i=0)

NumericVector< Number > &	get_adjoint_solution (unsigned int i=0)

const NumericVector< Number > &	get_adjoint_solution (unsigned int i=0) const

NumericVector< Number > &	add_sensitivity_solution (unsigned int i=0)

NumericVector< Number > &	get_sensitivity_solution (unsigned int i=0)

const NumericVector< Number > &	get_sensitivity_solution (unsigned int i=0) const

NumericVector< Number > &	add_weighted_sensitivity_adjoint_solution (unsigned int i=0)

NumericVector< Number > &	get_weighted_sensitivity_adjoint_solution (unsigned int i=0)

const NumericVector< Number > &	get_weighted_sensitivity_adjoint_solution (unsigned int i=0) const

NumericVector< Number > &	add_weighted_sensitivity_solution ()

NumericVector< Number > &	get_weighted_sensitivity_solution ()

const NumericVector< Number > &	get_weighted_sensitivity_solution () const

NumericVector< Number > &	add_adjoint_rhs (unsigned int i=0)

NumericVector< Number > &	get_adjoint_rhs (unsigned int i=0)

const NumericVector< Number > &	get_adjoint_rhs (unsigned int i=0) const

NumericVector< Number > &	add_sensitivity_rhs (unsigned int i=0)

NumericVector< Number > &	get_sensitivity_rhs (unsigned int i=0)

const NumericVector< Number > &	get_sensitivity_rhs (unsigned int i=0) const

unsigned int	n_vectors () const

unsigned int	n_vars () const

unsigned int	n_variable_groups () const

unsigned int	n_components () const

dof_id_type	n_dofs () const

dof_id_type	n_active_dofs () const

dof_id_type	n_constrained_dofs () const

dof_id_type	n_local_constrained_dofs () const

dof_id_type	n_local_dofs () const

unsigned int	add_variable (const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=nullptr)

unsigned int	add_variable (const std::string &var, const Order order=FIRST, const FEFamily=LAGRANGE, const std::set< subdomain_id_type > *const active_subdomains=nullptr)

unsigned int	add_variables (const std::vector< std::string > &vars, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=nullptr)

unsigned int	add_variables (const std::vector< std::string > &vars, const Order order=FIRST, const FEFamily=LAGRANGE, const std::set< subdomain_id_type > *const active_subdomains=nullptr)

const Variable &	variable (unsigned int var) const

const VariableGroup &	variable_group (unsigned int vg) const

bool	has_variable (const std::string &var) const

const std::string &	variable_name (const unsigned int i) const

unsigned short int	variable_number (const std::string &var) const

void	get_all_variable_numbers (std::vector< unsigned int > &all_variable_numbers) const

unsigned int	variable_scalar_number (const std::string &var, unsigned int component) const

unsigned int	variable_scalar_number (unsigned int var_num, unsigned int component) const

const FEType &	variable_type (const unsigned int i) const

const FEType &	variable_type (const std::string &var) const

bool	identify_variable_groups () const

void	identify_variable_groups (const bool)

Real	calculate_norm (const NumericVector< Number > &v, unsigned int var, FEMNormType norm_type, std::set< unsigned int > *skip_dimensions=nullptr) const

Real	calculate_norm (const NumericVector< Number > &v, const SystemNorm &norm, std::set< unsigned int > *skip_dimensions=nullptr) const

void	read_header (Xdr &io, const std::string &version, const bool read_header=true, const bool read_additional_data=true, const bool read_legacy_format=false)

void	read_legacy_data (Xdr &io, const bool read_additional_data=true)

template<typename ValType >
void	read_serialized_data (Xdr &io, const bool read_additional_data=true)

void	read_serialized_data (Xdr &io, const bool read_additional_data=true)

template<typename InValType >
std::size_t	read_serialized_vectors (Xdr &io, const std::vector< NumericVector< Number > *> &vectors) const

std::size_t	read_serialized_vectors (Xdr &io, const std::vector< NumericVector< Number > *> &vectors) const

template<typename InValType >
void	read_parallel_data (Xdr &io, const bool read_additional_data)

void	read_parallel_data (Xdr &io, const bool read_additional_data)

void	write_header (Xdr &io, const std::string &version, const bool write_additional_data) const

void	write_serialized_data (Xdr &io, const bool write_additional_data=true) const

std::size_t	write_serialized_vectors (Xdr &io, const std::vector< const NumericVector< Number > *> &vectors) const

void	write_parallel_data (Xdr &io, const bool write_additional_data) const

std::string	get_info () const

void	attach_init_function (void fptr(EquationSystems &es, const std::string &name))

void	attach_init_object (Initialization &init)

void	attach_assemble_function (void fptr(EquationSystems &es, const std::string &name))

void	attach_assemble_object (Assembly &assemble)

void	attach_constraint_function (void fptr(EquationSystems &es, const std::string &name))

void	attach_constraint_object (Constraint &constrain)

void	attach_QOI_function (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))

void	attach_QOI_object (QOI &qoi)

void	attach_QOI_derivative (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices, bool include_liftfunc, bool apply_constraints))

void	attach_QOI_derivative_object (QOIDerivative &qoi_derivative)

virtual void	user_initialization ()

virtual void	user_assembly ()

virtual void	user_constrain ()

virtual void	user_QOI (const QoISet &qoi_indices)

virtual void	user_QOI_derivative (const QoISet &qoi_indices=QoISet(), bool include_liftfunc=true, bool apply_constraints=true)

virtual void	re_update ()

virtual void	restrict_vectors ()

virtual void	prolong_vectors ()

Number	current_solution (const dof_id_type global_dof_number) const

unsigned int	n_qois () const

Number	point_value (unsigned int var, const Point &p, const bool insist_on_success=true) const

Number	point_value (unsigned int var, const Point &p, const Elem &e) const

Number	point_value (unsigned int var, const Point &p, const Elem *e) const

Gradient	point_gradient (unsigned int var, const Point &p, const bool insist_on_success=true) const

Gradient	point_gradient (unsigned int var, const Point &p, const Elem &e) const

Gradient	point_gradient (unsigned int var, const Point &p, const Elem *e) const

Tensor	point_hessian (unsigned int var, const Point &p, const bool insist_on_success=true) const

Tensor	point_hessian (unsigned int var, const Point &p, const Elem &e) const

Tensor	point_hessian (unsigned int var, const Point &p, const Elem *e) const

void	local_dof_indices (const unsigned int var, std::set< dof_id_type > &var_indices) const

void	zero_variable (NumericVector< Number > &v, unsigned int var_num) const

bool &	hide_output ()

void	projection_matrix (SparseMatrix< Number > &proj_mat) const

const Parallel::Communicator &	comm () const

processor_id_type	n_processors () const

processor_id_type	processor_id () const

virtual void	clear_physics ()

virtual void	init_physics (const System &sys)

virtual bool	element_time_derivative (bool request_jacobian, DiffContext &)

virtual bool	element_constraint (bool request_jacobian, DiffContext &)

virtual bool	side_time_derivative (bool request_jacobian, DiffContext &)

virtual bool	side_constraint (bool request_jacobian, DiffContext &)

virtual bool	nonlocal_time_derivative (bool request_jacobian, DiffContext &)

virtual bool	nonlocal_constraint (bool request_jacobian, DiffContext &)

virtual void	time_evolving (unsigned int var)

virtual void	time_evolving (unsigned int var, unsigned int order)

bool	is_time_evolving (unsigned int var) const

virtual bool	eulerian_residual (bool request_jacobian, DiffContext &)

virtual bool	eulerian_residual (bool request_jacobian, DiffContext &context) override

virtual bool	mass_residual (bool request_jacobian, DiffContext &)

virtual bool	mass_residual (bool request_jacobian, DiffContext &) override

virtual bool	side_mass_residual (bool request_jacobian, DiffContext &)

virtual bool	nonlocal_mass_residual (bool request_jacobian, DiffContext &c)

virtual bool	damping_residual (bool request_jacobian, DiffContext &)

virtual bool	side_damping_residual (bool request_jacobian, DiffContext &)

virtual bool	nonlocal_damping_residual (bool request_jacobian, DiffContext &)

virtual void	set_mesh_system (System *sys)

const System *	get_mesh_system () const

System *	get_mesh_system ()

virtual void	set_mesh_x_var (unsigned int var)

unsigned int	get_mesh_x_var () const

virtual void	set_mesh_y_var (unsigned int var)

unsigned int	get_mesh_y_var () const

virtual void	set_mesh_z_var (unsigned int var)

unsigned int	get_mesh_z_var () const

bool	_eulerian_time_deriv (bool request_jacobian, DiffContext &)

bool	have_first_order_vars () const

const std::set< unsigned int > &	get_first_order_vars () const

bool	is_first_order_var (unsigned int var) const

bool	have_second_order_vars () const

const std::set< unsigned int > &	get_second_order_vars () const

bool	is_second_order_var (unsigned int var) const

virtual void	init_qoi (std::vector< Number > &)

virtual void	clear_qoi ()

virtual void	element_qoi (DiffContext &, const QoISet &)

virtual void	element_qoi_derivative (DiffContext &, const QoISet &)

virtual void	side_qoi (DiffContext &, const QoISet &)

virtual void	side_qoi_derivative (DiffContext &, const QoISet &)

virtual void	thread_join (std::vector< Number > &qoi, const std::vector< Number > &other_qoi, const QoISet &qoi_indices)

virtual void	parallel_op (const Parallel::Communicator &communicator, std::vector< Number > &sys_qoi, std::vector< Number > &local_qoi, const QoISet &qoi_indices)

virtual void	finalize_derivative (NumericVector< Number > &derivatives, std::size_t qoi_index)

Static Public Member Functions
static std::string	get_info ()

static void	print_info (std::ostream &out=libMesh::out)

static unsigned int	n_objects ()

static void	enable_print_counter_info ()

static void	disable_print_counter_info ()

Public Attributes
Real *	continuation_parameter

bool	quiet

Real	continuation_parameter_tolerance

Real	solution_tolerance

Real	initial_newton_tolerance

Real	old_continuation_parameter

Real	min_continuation_parameter

Real	max_continuation_parameter

Real	Theta

Real	Theta_LOCA

unsigned int	n_backtrack_steps

unsigned int	n_arclength_reductions

Real	ds_min

Predictor	predictor

Real	newton_stepgrowth_aggressiveness

bool	newton_progress_check

bool	fe_reinit_during_postprocess

Real	numerical_jacobian_h

Real	verify_analytic_jacobians

std::unique_ptr< TimeSolver >	time_solver

Real	deltat

bool	postprocess_sides

bool	print_solution_norms

bool	print_solutions

bool	print_residual_norms

bool	print_residuals

bool	print_jacobian_norms

bool	print_jacobians

bool	print_element_solutions

bool	print_element_residuals

bool	print_element_jacobians

SparseMatrix< Number > *	matrix

bool	zero_out_matrix_and_rhs

NumericVector< Number > *	rhs

bool	assemble_before_solve

bool	use_fixed_solution

int	extra_quadrature_order

std::unique_ptr< NumericVector< Number > >	solution

std::unique_ptr< NumericVector< Number > >	current_local_solution

Real	time

std::vector< Number >	qoi

bool	compute_internal_sides

bool	assemble_qoi_sides

bool	assemble_qoi_internal_sides

bool	assemble_qoi_elements

Protected Types
enum	RHS_Mode { Residual, G_Lambda }

typedef std::map< std::string, std::pair< unsigned int, unsigned int > >	Counts

Protected Member Functions
virtual void	init_data () override

void	add_second_order_dot_vars ()

void	add_dot_var_dirichlet_bcs (unsigned int var_idx, unsigned int dot_var_idx)

virtual void	init_matrices ()

void	project_vector (NumericVector< Number > &, int is_adjoint=-1) const

void	project_vector (const NumericVector< Number > &, NumericVector< Number > &, int is_adjoint=-1) const

void	increment_constructor_count (const std::string &name)

void	increment_destructor_count (const std::string &name)

Protected Attributes
RHS_Mode	rhs_mode

DifferentiablePhysics *	_diff_physics

DifferentiableQoI *	diff_qoi

const Parallel::Communicator &	_communicator

System *	_mesh_sys

unsigned int	_mesh_x_var

unsigned int	_mesh_y_var

unsigned int	_mesh_z_var

std::vector< unsigned int >	_time_evolving

std::set< unsigned int >	_first_order_vars

std::set< unsigned int >	_second_order_vars

std::map< unsigned int, unsigned int >	_second_order_dot_vars

Static Protected Attributes
static Counts	_counts

static Threads::atomic< unsigned int >	_n_objects

static Threads::spin_mutex	_mutex

static bool	_enable_print_counter = true

Private Member Functions
void	initialize_tangent ()

void	solve_tangent ()

void	update_solution ()

void	set_Theta ()

void	set_Theta_LOCA ()

void	apply_predictor ()

Private Attributes
NumericVector< Number > *	du_ds

NumericVector< Number > *	previous_du_ds

NumericVector< Number > *	previous_u

NumericVector< Number > *	y

NumericVector< Number > *	y_old

NumericVector< Number > *	z

NumericVector< Number > *	delta_u

std::unique_ptr< LinearSolver< Number > >	linear_solver

bool	tangent_initialized

NewtonSolver *	newton_solver

Real	dlambda_ds

Real	ds

Real	ds_current

Real	previous_dlambda_ds

Real	previous_ds

unsigned int	newton_step

Detailed Description

This class inherits from the FEMSystem. It can be used to do arclength continuation. Most of the ideas and the notation here come from HB Keller's 1977 paper:

* @InProceedings{Kell-1977,
*   author    = {H.~B.~Keller},
*   title     = {{Numerical solution of bifurcation and nonlinear eigenvalue problems}},
*   booktitle = {Applications of Bifurcation Theory, P.~H.~Rabinowitz (ed.)},
*   year      = 1977,
*   publisher = {Academic Press},
*   pages     = {359--389},
*   notes     = {QA 3 U45 No.\ 38 (PMA)}
* }
*

Author: John W. Peterson

Date: 2007

Definition at line 55 of file continuation_system.h.

Member Typedef Documentation

◆ const_matrices_iterator

typedef std::map<std::string, SparseMatrix<Number> *>::const_iterator libMesh::ImplicitSystem::const_matrices_iterator

inherited

Definition at line 306 of file implicit_system.h.

◆ const_vectors_iterator

typedef std::map<std::string, NumericVector<Number> *>::const_iterator libMesh::System::const_vectors_iterator

inherited

Definition at line 749 of file system.h.

◆ Counts

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > libMesh::ReferenceCounter::Counts

protectedinherited

Data structure to log the information. The log is identified by the class name.

Definition at line 117 of file reference_counter.h.

◆ GradientFunctionPointer

typedef Gradient(* libMesh::System::GradientFunctionPointer) (const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name)

inherited

Definition at line 524 of file system.h.

◆ matrices_iterator

typedef std::map<std::string, SparseMatrix<Number> *>::iterator libMesh::ImplicitSystem::matrices_iterator

inherited

Matrix iterator typedefs.

Definition at line 305 of file implicit_system.h.

◆ Parent

typedef FEMSystem libMesh::ContinuationSystem::Parent

The type of the parent.

Definition at line 79 of file continuation_system.h.

◆ sys_type

typedef ContinuationSystem libMesh::ContinuationSystem::sys_type

The type of system.

Definition at line 74 of file continuation_system.h.

◆ TimeSolverResPtr

typedef bool(TimeSolver::* libMesh::FEMSystem::TimeSolverResPtr) (bool, DiffContext &)

inherited

Syntax sugar to make numerical_jacobian() declaration easier.

Definition at line 214 of file fem_system.h.

◆ ValueFunctionPointer

typedef Number(* libMesh::System::ValueFunctionPointer) (const Point &p, const Parameters &Parameters, const std::string &sys_name, const std::string &unknown_name)

inherited

Projects arbitrary functions onto the current solution. The function value fptr and its gradient gptr are represented by function pointers. A gradient gptr is only required/used for projecting onto finite element spaces with continuous derivatives.

Definition at line 520 of file system.h.

◆ vectors_iterator

typedef std::map<std::string, NumericVector<Number> *>::iterator libMesh::System::vectors_iterator

inherited

Vector iterator typedefs.

Definition at line 748 of file system.h.

Member Enumeration Documentation

◆ Predictor

enum libMesh::ContinuationSystem::Predictor

The code provides the ability to select from different predictor schemes for getting the initial guess for the solution at the next point on the solution arc.

Enumerator
Euler	First-order Euler predictor
AB2	Second-order explicit Adams-Bashforth predictor
Invalid_Predictor	Invalid predictor

Definition at line 221 of file continuation_system.h.

                  {
     Euler,
 
     AB2,
 
     Invalid_Predictor
   };

◆ RHS_Mode

enum libMesh::ContinuationSystem::RHS_Mode

protected

There are (so far) two different vectors which may be assembled using the assembly routine: 1.) The Residual = the normal PDE weighted residual 2.) G_Lambda = the derivative wrt the parameter lambda of the PDE weighted residual

It is up to the derived class to handle writing separate assembly code for the different cases. Usually something like: switch (rhs_mode) { case Residual: { Fu(i) += ... // normal PDE residual break; }

case G_Lambda: { Fu(i) += ... // derivative wrt control parameter break; }

Enumerator
Residual
G_Lambda

Definition at line 286 of file continuation_system.h.

286 {Residual,

287 G_Lambda};

libMesh::ContinuationSystem::Residual

Definition: continuation_system.h:286

libMesh::ContinuationSystem::G_Lambda

Definition: continuation_system.h:287

Constructor & Destructor Documentation

◆ ContinuationSystem()

libMesh::ContinuationSystem::ContinuationSystem	(	EquationSystems &	es,
		const std::string &	name,
		const unsigned int	number
	)

Constructor. Optionally initializes required data structures.

Definition at line 29 of file continuation_system.C.

References linear_solver, libMesh::System::name(), and libMesh::on_command_line().

                                                                       :
   Parent(es, name_in, number_in),
   continuation_parameter(nullptr),
   quiet(true),
   continuation_parameter_tolerance(1.e-6),
   solution_tolerance(1.e-6),
   initial_newton_tolerance(0.01),
   old_continuation_parameter(0.),
   min_continuation_parameter(0.),
   max_continuation_parameter(0.),
   Theta(1.),
   Theta_LOCA(1.),
   n_backtrack_steps(5),
   n_arclength_reductions(5),
   ds_min(1.e-8),
   predictor(Euler),
   newton_stepgrowth_aggressiveness(1.),
   newton_progress_check(true),
   rhs_mode(Residual),
   linear_solver(LinearSolver<Number>::build(es.comm())),
   tangent_initialized(false),
   newton_solver(nullptr),
   dlambda_ds(0.707),
   ds(0.1),
   ds_current(0.1),
   previous_dlambda_ds(0.),
   previous_ds(0.),
   newton_step(0)
 {
   // Warn about using untested code
   libmesh_experimental();
 
   if (libMesh::on_command_line("--solver-system-names"))
     linear_solver->init((this->name()+"_").c_str());
   else
     linear_solver->init();
 }

◆ ~ContinuationSystem()

libMesh::ContinuationSystem::~ContinuationSystem ( )

virtual

Destructor.

Definition at line 72 of file continuation_system.C.

References clear().

 {
   this->clear();
 }

Member Function Documentation

◆ _eulerian_time_deriv()

bool libMesh::DifferentiablePhysics::_eulerian_time_deriv	(	bool	request_jacobian,
		DiffContext &
	)

inherited

This method simply combines element_time_derivative() and eulerian_residual(), which makes its address useful as a pointer-to-member-function when refactoring.

Referenced by libMesh::EulerSolver::element_residual(), libMesh::Euler2Solver::element_residual(), and libMesh::NewmarkSolver::element_residual().

◆ activate()

void libMesh::System::activate ( )

inlineinherited

Activates the system. Only active systems are solved.

Definition at line 2073 of file system.h.

References libMesh::System::_active.

 {
   _active = true;
 }

◆ active()

bool libMesh::System::active ( ) const

inlineinherited

Returns: true if the system is active, false otherwise. An active system will be solved.

Definition at line 2065 of file system.h.

References libMesh::System::_active.

 {
   return _active;
 }

◆ add_adjoint_rhs()

NumericVector< Number > & libMesh::System::add_adjoint_rhs ( unsigned int i = 0 )

inherited

Returns: A reference to one of the system's adjoint rhs vectors, by default the one corresponding to the first qoi. Creates the vector if it doesn't already exist.

Definition at line 1021 of file system.C.

References libMesh::System::add_vector().

Referenced by libMesh::ExplicitSystem::assemble_qoi_derivative(), and libMesh::FEMSystem::assemble_qoi_derivative().

 {
   std::ostringstream adjoint_rhs_name;
   adjoint_rhs_name << "adjoint_rhs" << i;
 
   return this->add_vector(adjoint_rhs_name.str(), false);
 }

◆ add_adjoint_solution()

NumericVector< Number > & libMesh::System::add_adjoint_solution ( unsigned int i = 0 )

inherited

Returns: A reference to one of the system's adjoint solution vectors, by default the one corresponding to the first qoi. Creates the vector if it doesn't already exist.

Definition at line 957 of file system.C.

References libMesh::System::add_vector(), and libMesh::System::set_vector_as_adjoint().

Referenced by libMesh::ImplicitSystem::adjoint_solve().

 {
   std::ostringstream adjoint_name;
   adjoint_name << "adjoint_solution" << i;
 
   NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
   this->set_vector_as_adjoint(adjoint_name.str(), i);
   return returnval;
 }

◆ add_dot_var_dirichlet_bcs()

void libMesh::DifferentiableSystem::add_dot_var_dirichlet_bcs	(	unsigned int	var_idx,
		unsigned int	dot_var_idx
	)

protectedinherited

Helper function to and Dirichlet boundary conditions to "dot" variable cousins of second order variables in the system. The function takes the second order variable index, it's corresponding "dot" variable index and then searches for DirichletBoundary objects for var_idx and then adds a DirichletBoundary object for dot_var_idx using the same boundary ids and functors for the var_idx DirichletBoundary.

Definition at line 227 of file diff_system.C.

References libMesh::DofMap::add_dirichlet_boundary(), libMesh::DofMap::get_dirichlet_boundaries(), and libMesh::System::get_dof_map().

Referenced by libMesh::DifferentiableSystem::add_second_order_dot_vars().

 {
   // We're assuming that there could be a lot more variables than
   // boundary conditions, so we search each of the boundary conditions
   // for this variable rather than looping over boundary conditions
   // in a separate loop and searching through all the variables.
   const DirichletBoundaries * all_dbcs =
     this->get_dof_map().get_dirichlet_boundaries();
 
   if (all_dbcs)
     {
       // We need to cache the DBCs to be added so that we add them
       // after looping over the existing DBCs. Otherwise, we're polluting
       // the thing we're looping over.
       std::vector<DirichletBoundary *> new_dbcs;
 
       for (const auto & dbc : *all_dbcs)
         {
           libmesh_assert(dbc);
 
           // Look for second order variable in the current
           // DirichletBoundary object
           std::vector<unsigned int>::const_iterator dbc_var_it =
             std::find( dbc->variables.begin(), dbc->variables.end(), var_idx );
 
           // If we found it, then we also need to add it's corresponding
           // "dot" variable to a DirichletBoundary
           std::vector<unsigned int> vars_to_add;
           if (dbc_var_it != dbc->variables.end())
             vars_to_add.push_back(dot_var_idx);
 
           if (!vars_to_add.empty())
             {
               // We need to check if the boundary condition is time-dependent.
               // Currently, we cannot automatically differentiate w.r.t. time
               // so if the user supplies a time-dependent Dirichlet BC, then
               // we can't automatically support the Dirichlet BC for the
               // "velocity" boundary condition, so we error. Otherwise,
               // the "velocity boundary condition will just be zero.
               bool is_time_evolving_bc = false;
               if (dbc->f)
                 is_time_evolving_bc = dbc->f->is_time_dependent();
               else if (dbc->f_fem)
                 // We it's a FEMFunctionBase object, it will be implicitly
                 // time-dependent since it is assumed to depend on the solution.
                 is_time_evolving_bc = true;
               else
                 libmesh_error_msg("Could not find valid boundary function!");
 
               if (is_time_evolving_bc)
                 libmesh_error_msg("Cannot currently support time-dependent Dirichlet BC for dot variables!");
 
 
               DirichletBoundary * new_dbc;
 
               if (dbc->f)
                 {
                   ZeroFunction<Number> zf;
 
                   new_dbc = new DirichletBoundary(dbc->b, vars_to_add, zf);
                 }
               else
                 libmesh_error();
 
               new_dbcs.push_back(new_dbc);
             }
         }
 
       // Let the DofMap make its own deep copy of the DirichletBC objects and delete our copy.
       for (const auto & dbc : new_dbcs)
         {
           this->get_dof_map().add_dirichlet_boundary(*dbc);
           delete dbc;
         }
 
     } // if (all_dbcs)
 }

◆ add_matrix()

SparseMatrix< Number > & libMesh::ImplicitSystem::add_matrix ( const std::string & mat_name )

inherited

Adds the additional matrix mat_name to this system. Only allowed prior to assemble(). All additional matrices have the same sparsity pattern as the matrix used during solution. When not System but the user wants to initialize the mayor matrix, then all the additional matrices, if existent, have to be initialized by the user, too.

Definition at line 203 of file implicit_system.C.

References libMesh::ImplicitSystem::_can_add_matrices, libMesh::ImplicitSystem::_matrices, libMesh::SparseMatrix< T >::build(), libMesh::ParallelObject::comm(), and libMesh::ImplicitSystem::have_matrix().

Referenced by libMesh::ImplicitSystem::add_system_matrix(), libMesh::EigenTimeSolver::init(), and libMesh::NewmarkSystem::NewmarkSystem().

 {
   // only add matrices before initializing...
   if (!_can_add_matrices)
     libmesh_error_msg("ERROR: Too late.  Cannot add matrices to the system after initialization"
                       << "\n any more.  You should have done this earlier.");
 
   // Return the matrix if it is already there.
   if (this->have_matrix(mat_name))
     return *(_matrices[mat_name]);
 
   // Otherwise build the matrix and return it.
   SparseMatrix<Number> * buf = SparseMatrix<Number>::build(this->comm()).release();
   _matrices.insert (std::make_pair (mat_name, buf));
 
   return *buf;
 }

◆ add_second_order_dot_vars()

void libMesh::DifferentiableSystem::add_second_order_dot_vars ( )

protectedinherited

Helper function to add "velocity" variables that are cousins to second order-in-time variables in the DifferentiableSystem. This function is only called if the TimeSolver is a FirstOrderUnsteadySolver.

Definition at line 198 of file diff_system.C.

References libMesh::DifferentiablePhysics::_second_order_dot_vars, libMesh::Variable::active_subdomains(), libMesh::DifferentiableSystem::add_dot_var_dirichlet_bcs(), libMesh::System::add_variable(), libMesh::DifferentiablePhysics::get_second_order_vars(), libMesh::Variable::name(), libMesh::DifferentiablePhysics::time_evolving(), libMesh::Variable::type(), and libMesh::System::variable().

Referenced by libMesh::DifferentiableSystem::init_data().

 {
   const std::set<unsigned int> & second_order_vars = this->get_second_order_vars();
   if (!second_order_vars.empty())
     {
       for (const auto & var_id : second_order_vars)
         {
           const Variable & var = this->variable(var_id);
           std::string new_var_name = std::string("dot_")+var.name();
 
           unsigned int v_var_idx;
 
           if (var.active_subdomains().empty())
             v_var_idx = this->add_variable( new_var_name, var.type() );
           else
             v_var_idx = this->add_variable( new_var_name, var.type(), &var.active_subdomains() );
 
           _second_order_dot_vars.insert(std::pair<unsigned int, unsigned int>(var_id, v_var_idx));
 
           // The new velocities are time evolving variables of first order
           this->time_evolving( v_var_idx, 1 );
 
           // And if there are any boundary conditions set on the second order
           // variable, we also need to set it on its velocity variable.
           this->add_dot_var_dirichlet_bcs(var_id, v_var_idx);
         }
     }
 }

◆ add_sensitivity_rhs()

NumericVector< Number > & libMesh::System::add_sensitivity_rhs ( unsigned int i = 0 )

inherited

Returns: A reference to one of the system's sensitivity rhs vectors, by default the one corresponding to the first parameter. Creates the vector if it doesn't already exist.

Definition at line 1051 of file system.C.

References libMesh::System::add_vector().

Referenced by libMesh::ImplicitSystem::assemble_residual_derivatives().

 {
   std::ostringstream sensitivity_rhs_name;
   sensitivity_rhs_name << "sensitivity_rhs" << i;
 
   return this->add_vector(sensitivity_rhs_name.str(), false);
 }

◆ add_sensitivity_solution()

NumericVector< Number > & libMesh::System::add_sensitivity_solution ( unsigned int i = 0 )

inherited

Returns: A reference to one of the system's solution sensitivity vectors, by default the one corresponding to the first parameter. Creates the vector if it doesn't already exist.

Definition at line 906 of file system.C.

References libMesh::System::add_vector().

Referenced by libMesh::ImplicitSystem::sensitivity_solve().

 {
   std::ostringstream sensitivity_name;
   sensitivity_name << "sensitivity_solution" << i;
 
   return this->add_vector(sensitivity_name.str());
 }

◆ add_variable() [1/2]

unsigned int libMesh::System::add_variable	(	const std::string &	var,
		const FEType &	type,
		const std::set< subdomain_id_type > *const	active_subdomains = `nullptr`
	)

inherited

Adds the variable var to the list of variables for this system.

Returns: The index number for the new variable.

Definition at line 1081 of file system.C.

References libMesh::System::_variable_groups, libMesh::System::_variable_numbers, libMesh::System::_variables, libMesh::System::add_variables(), libMesh::VariableGroup::append(), libMesh::System::identify_variable_groups(), libMesh::System::is_initialized(), libMesh::System::n_variable_groups(), libMesh::System::n_vars(), libMesh::System::number(), libMesh::System::variable_name(), and libMesh::System::variable_type().

Referenced by libMesh::DifferentiableSystem::add_second_order_dot_vars(), libMesh::System::add_variable(), libMesh::ErrorVector::plot_error(), and libMesh::System::read_header().

 {
   libmesh_assert(!this->is_initialized());
 
   // Make sure the variable isn't there already
   // or if it is, that it's the type we want
   for (unsigned int v=0; v<this->n_vars(); v++)
     if (this->variable_name(v) == var)
       {
         if (this->variable_type(v) == type)
           return _variables[v].number();
 
         libmesh_error_msg("ERROR: incompatible variable " << var << " has already been added for this system!");
       }
 
   // Optimize for VariableGroups here - if the user is adding multiple
   // variables of the same FEType and subdomain restriction, catch
   // that here and add them as members of the same VariableGroup.
   //
   // start by setting this flag to whatever the user has requested
   // and then consider the conditions which should negate it.
   bool should_be_in_vg = this->identify_variable_groups();
 
   // No variable groups, nothing to add to
   if (!this->n_variable_groups())
     should_be_in_vg = false;
 
   else
     {
       VariableGroup & vg(_variable_groups.back());
 
       // get a pointer to their subdomain restriction, if any.
       const std::set<subdomain_id_type> * const
         their_active_subdomains (vg.implicitly_active() ?
                                  nullptr : &vg.active_subdomains());
 
       // Different types?
       if (vg.type() != type)
         should_be_in_vg = false;
 
       // they are restricted, we aren't?
       if (their_active_subdomains && !active_subdomains)
         should_be_in_vg = false;
 
       // they aren't restricted, we are?
       if (!their_active_subdomains && active_subdomains)
         should_be_in_vg = false;
 
       if (their_active_subdomains && active_subdomains)
         // restricted to different sets?
         if (*their_active_subdomains != *active_subdomains)
           should_be_in_vg = false;
 
       // OK, after all that, append the variable to the vg if none of the conditions
       // were violated
       if (should_be_in_vg)
         {
           const unsigned short curr_n_vars = cast_int<unsigned short>
             (this->n_vars());
 
           vg.append (var);
 
           _variables.push_back(vg(vg.n_variables()-1));
           _variable_numbers[var] = curr_n_vars;
           return curr_n_vars;
         }
     }
 
   // otherwise, fall back to adding a single variable group
   return this->add_variables (std::vector<std::string>(1, var),
                               type,
                               active_subdomains);
 }

◆ add_variable() [2/2]

unsigned int libMesh::System::add_variable	(	const std::string &	var,
		const Order	order = `FIRST`,
		const FEFamily	family = `LAGRANGE`,
		const std::set< subdomain_id_type > *const	active_subdomains = `nullptr`
	)

inherited

Adds the variable var to the list of variables for this system. Same as before, but assumes LAGRANGE as default value for FEType.family.

Definition at line 1159 of file system.C.

References libMesh::System::add_variable().

 {
   return this->add_variable(var,
                             FEType(order, family),
                             active_subdomains);
 }

◆ add_variables() [1/2]

unsigned int libMesh::System::add_variables	(	const std::vector< std::string > &	vars,
		const FEType &	type,
		const std::set< subdomain_id_type > *const	active_subdomains = `nullptr`
	)

inherited

Adds the variable var to the list of variables for this system.

Returns: The index number for the new variable.

Definition at line 1171 of file system.C.

References libMesh::System::_variable_groups, libMesh::System::_variable_numbers, libMesh::System::_variables, libMesh::System::is_initialized(), libMesh::System::n_components(), libMesh::System::n_vars(), libMesh::System::number(), libMesh::System::variable_name(), and libMesh::System::variable_type().

Referenced by libMesh::System::add_variable(), and libMesh::System::add_variables().

 {
   libmesh_assert(!this->is_initialized());
 
   // Make sure the variable isn't there already
   // or if it is, that it's the type we want
   for (std::size_t ov=0; ov<vars.size(); ov++)
     for (unsigned int v=0; v<this->n_vars(); v++)
       if (this->variable_name(v) == vars[ov])
         {
           if (this->variable_type(v) == type)
             return _variables[v].number();
 
           libmesh_error_msg("ERROR: incompatible variable " << vars[ov] << " has already been added for this system!");
         }
 
   const unsigned short curr_n_vars = cast_int<unsigned short>
     (this->n_vars());
 
   const unsigned int next_first_component = this->n_components();
 
   // Add the variable group to the list
   _variable_groups.push_back((active_subdomains == nullptr) ?
                              VariableGroup(this, vars, curr_n_vars,
                                            next_first_component, type) :
                              VariableGroup(this, vars, curr_n_vars,
                                            next_first_component, type, *active_subdomains));
 
   const VariableGroup & vg (_variable_groups.back());
 
   // Add each component of the group individually
   for (auto v : IntRange<unsigned int>(0, vars.size()))
     {
       _variables.push_back (vg(v));
       _variable_numbers[vars[v]] = cast_int<unsigned short>
         (curr_n_vars+v);
     }
 
   libmesh_assert_equal_to ((curr_n_vars+vars.size()), this->n_vars());
 
   // BSK - Defer this now to System::init_data() so we can detect
   // VariableGroups 12/28/2012
   // // Add the variable group to the _dof_map
   // _dof_map->add_variable_group (vg);
 
   // Return the number of the new variable
   return cast_int<unsigned int>(curr_n_vars+vars.size()-1);
 }

◆ add_variables() [2/2]

unsigned int libMesh::System::add_variables	(	const std::vector< std::string > &	vars,
		const Order	order = `FIRST`,
		const FEFamily	family = `LAGRANGE`,
		const std::set< subdomain_id_type > *const	active_subdomains = `nullptr`
	)

inherited

Adds the variable var to the list of variables for this system. Same as before, but assumes LAGRANGE as default value for FEType.family.

Definition at line 1224 of file system.C.

References libMesh::System::add_variables().

 {
   return this->add_variables(vars,
                              FEType(order, family),
                              active_subdomains);
 }

◆ add_vector()

NumericVector< Number > & libMesh::System::add_vector	(	const std::string &	vec_name,
		const bool	projections = `true`,
		const ParallelType	type = `PARALLEL`
	)

inherited

Adds the additional vector vec_name to this system. All the additional vectors are similarly distributed, like the solution, and initialized to zero.

By default vectors added by add_vector are projected to changed grids by reinit(). To zero them instead (more efficient), pass "false" as the second argument

Definition at line 661 of file system.C.

References libMesh::System::_dof_map, libMesh::System::_is_initialized, libMesh::System::_vector_is_adjoint, libMesh::System::_vector_projections, libMesh::System::_vector_types, libMesh::System::_vectors, libMesh::NumericVector< T >::build(), libMesh::ParallelObject::comm(), libMesh::GHOSTED, libMesh::System::have_vector(), libMesh::NumericVector< T >::init(), libMesh::System::n_dofs(), and libMesh::System::n_local_dofs().

Referenced by libMesh::System::add_adjoint_rhs(), libMesh::System::add_adjoint_solution(), libMesh::System::add_sensitivity_rhs(), libMesh::System::add_sensitivity_solution(), libMesh::ExplicitSystem::add_system_rhs(), libMesh::System::add_weighted_sensitivity_adjoint_solution(), libMesh::System::add_weighted_sensitivity_solution(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::SecondOrderUnsteadySolver::init(), libMesh::UnsteadySolver::init(), libMesh::OptimizationSystem::init_data(), init_data(), libMesh::NewmarkSystem::NewmarkSystem(), libMesh::System::read_header(), libMesh::FrequencySystem::set_frequencies(), libMesh::FrequencySystem::set_frequencies_by_range(), and libMesh::FrequencySystem::set_frequencies_by_steps().

 {
   // Return the vector if it is already there.
   if (this->have_vector(vec_name))
     return *(_vectors[vec_name]);
 
   // Otherwise build the vector
   NumericVector<Number> * buf = NumericVector<Number>::build(this->comm()).release();
   _vectors.insert (std::make_pair (vec_name, buf));
   _vector_projections.insert (std::make_pair (vec_name, projections));
 
   _vector_types.insert (std::make_pair (vec_name, type));
 
   // Vectors are primal by default
   _vector_is_adjoint.insert (std::make_pair (vec_name, -1));
 
   // Initialize it if necessary
   if (_is_initialized)
     {
       if (type == GHOSTED)
         {
 #ifdef LIBMESH_ENABLE_GHOSTED
           buf->init (this->n_dofs(), this->n_local_dofs(),
                      _dof_map->get_send_list(), false,
                      GHOSTED);
 #else
           libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
 #endif
         }
       else
         buf->init (this->n_dofs(), this->n_local_dofs(), false, type);
     }
 
   return *buf;
 }

◆ add_weighted_sensitivity_adjoint_solution()

NumericVector< Number > & libMesh::System::add_weighted_sensitivity_adjoint_solution ( unsigned int i = 0 )

inherited

Returns: A reference to one of the system's weighted sensitivity adjoint solution vectors, by default the one corresponding to the first qoi. Creates the vector if it doesn't already exist.

Definition at line 989 of file system.C.

References libMesh::System::add_vector(), and libMesh::System::set_vector_as_adjoint().

Referenced by libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve().

 {
   std::ostringstream adjoint_name;
   adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
 
   NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
   this->set_vector_as_adjoint(adjoint_name.str(), i);
   return returnval;
 }

◆ add_weighted_sensitivity_solution()

NumericVector< Number > & libMesh::System::add_weighted_sensitivity_solution ( )

inherited

Returns: A reference to the solution of the last weighted sensitivity solve Creates the vector if it doesn't already exist.

Definition at line 936 of file system.C.

References libMesh::System::add_vector().

Referenced by libMesh::ImplicitSystem::weighted_sensitivity_solve().

 {
   return this->add_vector("weighted_sensitivity_solution");
 }

◆ adjoint_qoi_parameter_sensitivity()

void libMesh::ImplicitSystem::adjoint_qoi_parameter_sensitivity	(	const QoISet &	qoi_indices,
		const ParameterVector &	parameters,
		SensitivityData &	sensitivities
	)

overridevirtualinherited

Solves for the derivative of each of the system's quantities of interest q in qoi[qoi_indices] with respect to each parameter in parameters, placing the result for qoi i and parameter j into sensitivities[i][j].

Uses adjoint_solve() and the adjoint sensitivity method.

Currently uses finite differenced derivatives (partial q / partial p) and (partial R / partial p).

Reimplemented from libMesh::System.

Definition at line 697 of file implicit_system.C.

References std::abs(), libMesh::ImplicitSystem::adjoint_solve(), libMesh::SensitivityData::allocate_data(), libMesh::ExplicitSystem::assemble_qoi(), libMesh::ImplicitSystem::assemble_residual_derivatives(), libMesh::NumericVector< T >::dot(), libMesh::System::get_sensitivity_rhs(), libMesh::QoISet::has_index(), libMesh::System::is_adjoint_already_solved(), std::max(), libMesh::System::n_qois(), libMesh::System::qoi, libMesh::Real, libMesh::ParameterVector::size(), and libMesh::TOLERANCE.

 {
   ParameterVector & parameters =
     const_cast<ParameterVector &>(parameters_in);
 
   const unsigned int Np = cast_int<unsigned int>
     (parameters.size());
   const unsigned int Nq = this->n_qois();
 
   // An introduction to the problem:
   //
   // Residual R(u(p),p) = 0
   // partial R / partial u = J = system matrix
   //
   // This implies that:
   // d/dp(R) = 0
   // (partial R / partial p) +
   // (partial R / partial u) * (partial u / partial p) = 0
 
   // We first do an adjoint solve:
   // J^T * z = (partial q / partial u)
   // if we havent already or dont have an initial condition for the adjoint
   if (!this->is_adjoint_already_solved())
     {
       this->adjoint_solve(qoi_indices);
     }
 
   this->assemble_residual_derivatives(parameters_in);
 
   // Get ready to fill in sensitivities:
   sensitivities.allocate_data(qoi_indices, *this, parameters);
 
   // We use the identities:
   // dq/dp = (partial q / partial p) + (partial q / partial u) *
   //         (partial u / partial p)
   // dq/dp = (partial q / partial p) + (J^T * z) *
   //         (partial u / partial p)
   // dq/dp = (partial q / partial p) + z * J *
   //         (partial u / partial p)
 
   // Leading to our final formula:
   // dq/dp = (partial q / partial p) - z * (partial R / partial p)
 
   // In the case of adjoints with heterogenous Dirichlet boundary
   // function phi, where
   // q :=  S(u) - R(u,phi)
   // the final formula works out to:
   // dq/dp = (partial S / partial p) - z * (partial R / partial p)
   // Because we currently have no direct access to
   // (partial S / partial p), we use the identity
   // (partial S / partial p) = (partial q / partial p) +
   //                           phi * (partial R / partial p)
   // to derive an equivalent equation:
   // dq/dp = (partial q / partial p) - (z-phi) * (partial R / partial p)
 
   // Since z-phi degrees of freedom are zero for constrained indices,
   // we can use the same constrained -(partial R / partial p) that we
   // use for forward sensitivity solves, taking into account the
   // differing sign convention.
   //
   // Since that vector is constrained, its constrained indices are
   // zero, so its product with phi is zero, so we can neglect the
   // evaluation of phi terms.
 
   for (unsigned int j=0; j != Np; ++j)
     {
       // We currently get partial derivatives via central differencing
 
       // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
       // (partial R / partial p) ~= (rhs(p+dp) - rhs(p-dp))/(2*dp)
 
       Number old_parameter = *parameters[j];
 
       const Real delta_p =
         TOLERANCE * std::max(std::abs(old_parameter), 1e-3);
 
       *parameters[j] = old_parameter - delta_p;
       this->assemble_qoi(qoi_indices);
       std::vector<Number> qoi_minus = this->qoi;
 
       NumericVector<Number> & neg_partialR_partialp = this->get_sensitivity_rhs(j);
 
       *parameters[j] = old_parameter + delta_p;
       this->assemble_qoi(qoi_indices);
       std::vector<Number> & qoi_plus = this->qoi;
 
       std::vector<Number> partialq_partialp(Nq, 0);
       for (unsigned int i=0; i != Nq; ++i)
         if (qoi_indices.has_index(i))
           partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
 
       // Don't leave the parameter changed
       *parameters[j] = old_parameter;
 
       for (unsigned int i=0; i != Nq; ++i)
         if (qoi_indices.has_index(i))
           sensitivities[i][j] = partialq_partialp[i] +
             neg_partialR_partialp.dot(this->get_adjoint_solution(i));
     }
 
   // All parameters have been reset.
   // Reset the original qoi.
 
   this->assemble_qoi(qoi_indices);
 }

◆ adjoint_solve()

std::pair< unsigned int, Real > libMesh::DifferentiableSystem::adjoint_solve ( const QoISet & qoi_indices = QoISet() )

overridevirtualinherited

This function sets the _is_adjoint boolean member of TimeSolver to true and then calls the adjoint_solve in implicit system

Reimplemented from libMesh::ImplicitSystem.

Definition at line 164 of file diff_system.C.

References libMesh::ImplicitSystem::adjoint_solve(), libMesh::DifferentiableSystem::get_time_solver(), and libMesh::TimeSolver::set_is_adjoint().

 {
   // Get the time solver object associated with the system, and tell it that
   // we are solving the adjoint problem
   this->get_time_solver().set_is_adjoint(true);
 
   return this->ImplicitSystem::adjoint_solve(qoi_indices);
 }

◆ advance_arcstep()

void libMesh::ContinuationSystem::advance_arcstep ( )

Call this function after a continuation solve to compute the tangent and get the next initial guess.

Definition at line 936 of file continuation_system.C.

References solve_tangent(), and update_solution().

 {
   // Solve for the updated tangent du1/ds, d(lambda1)/ds
   solve_tangent();
 
   // Advance the solution and the parameter to the next value.
   update_solution();
 }

◆ apply_predictor()

void libMesh::ContinuationSystem::apply_predictor ( )

private

Applies the predictor to the current solution to get a guess for the next solution.

Definition at line 1375 of file continuation_system.C.

References AB2, continuation_parameter, dlambda_ds, ds_current, du_ds, Euler, predictor, previous_dlambda_ds, previous_ds, previous_du_ds, libMesh::Real, and libMesh::System::solution.

Referenced by continuation_solve(), and update_solution().

 {
   if (predictor == Euler)
     {
       // 1.) Euler Predictor
       // Predict next the solution
       solution->add(ds_current, *du_ds);
       solution->close();
 
       // Predict next parameter value
       *continuation_parameter += ds_current*dlambda_ds;
     }
 
 
   else if (predictor == AB2)
     {
       // 2.) 2nd-order explicit AB predictor
       libmesh_assert_not_equal_to (previous_ds, 0.0);
       const Real stepratio = ds_current/previous_ds;
 
       // Build up next solution value.
       solution->add( 0.5*ds_current*(2.+stepratio), *du_ds);
       solution->add(-0.5*ds_current*stepratio     , *previous_du_ds);
       solution->close();
 
       // Next parameter value
       *continuation_parameter +=
         0.5*ds_current*((2.+stepratio)*dlambda_ds -
                         stepratio*previous_dlambda_ds);
     }
 
   else
     libmesh_error_msg("Unknown predictor!");
 }

◆ assemble()

void libMesh::DifferentiableSystem::assemble ( )

overridevirtualinherited

Prepares matrix and rhs for matrix assembly. Users should not reimplement this. Note that in some cases only current_local_solution is used during assembly, and, therefore, if solution has been altered without update() being called, then the user must call update() before calling this function.

Reimplemented from libMesh::ImplicitSystem.

Definition at line 145 of file diff_system.C.

References libMesh::DifferentiableSystem::assembly().

 {
   this->assembly(true, true);
 }

◆ assemble_qoi()

void libMesh::FEMSystem::assemble_qoi ( const QoISet & indices = QoISet() )

overridevirtualinherited

Runs a qoi assembly loop over all elements, and if assemble_qoi_sides is true over all sides.

Users may have to override this function if they have any quantities of interest that are not expressible as a sum of element qois.

Reimplemented from libMesh::ExplicitSystem.

Definition at line 1121 of file fem_system.C.

References libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::ParallelObject::comm(), libMesh::DifferentiableSystem::diff_qoi, libMesh::System::get_mesh(), libMesh::QoISet::has_index(), mesh, libMesh::System::n_qois(), libMesh::DifferentiableQoI::parallel_op(), libMesh::Threads::parallel_reduce(), libMesh::System::qoi, libMesh::StoredRange< iterator_type, object_type >::reset(), and libMesh::System::update().

 {
   LOG_SCOPE("assemble_qoi()", "FEMSystem");
 
   const MeshBase & mesh = this->get_mesh();
 
   this->update();
 
   const unsigned int Nq = this->n_qois();
 
   // the quantity of interest is assumed to be a sum of element and
   // side terms
   for (unsigned int i=0; i != Nq; ++i)
     if (qoi_indices.has_index(i))
       qoi[i] = 0;
 
   // Create a non-temporary qoi_contributions object, so we can query
   // its results after the reduction
   QoIContributions qoi_contributions(*this, *(this->diff_qoi), qoi_indices);
 
   // Loop over every active mesh element on this processor
   Threads::parallel_reduce(elem_range.reset(mesh.active_local_elements_begin(),
                                             mesh.active_local_elements_end()),
                            qoi_contributions);
 
   this->diff_qoi->parallel_op( this->comm(), this->qoi, qoi_contributions.qoi, qoi_indices );
 }

◆ assemble_qoi_derivative()

void libMesh::FEMSystem::assemble_qoi_derivative	(	const QoISet &	qoi_indices = `QoISet()`,
		bool	include_liftfunc = `true`,
		bool	apply_constraints = `true`
	)

overridevirtualinherited

Runs a qoi derivative assembly loop over all elements, and if assemble_qoi_sides is true over all sides.

Users may have to override this function for quantities of interest that are not expressible as a sum of element qois.

Reimplemented from libMesh::ExplicitSystem.

Definition at line 1151 of file fem_system.C.

References libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::System::add_adjoint_rhs(), libMesh::DifferentiableSystem::diff_qoi, libMesh::DifferentiableQoI::finalize_derivative(), libMesh::System::get_mesh(), libMesh::QoISet::has_index(), mesh, libMesh::System::n_qois(), libMesh::Threads::parallel_for(), libMesh::StoredRange< iterator_type, object_type >::reset(), libMesh::System::update(), and libMesh::NumericVector< T >::zero().

 {
   LOG_SCOPE("assemble_qoi_derivative()", "FEMSystem");
 
   const MeshBase & mesh = this->get_mesh();
 
   this->update();
 
   // The quantity of interest derivative assembly accumulates on
   // initially zero vectors
   for (unsigned int i=0; i != this->n_qois(); ++i)
     if (qoi_indices.has_index(i))
       this->add_adjoint_rhs(i).zero();
 
   // Loop over every active mesh element on this processor
   Threads::parallel_for (elem_range.reset(mesh.active_local_elements_begin(),
                                           mesh.active_local_elements_end()),
                          QoIDerivativeContributions(*this, qoi_indices,
                                                     *(this->diff_qoi),
                                                     include_liftfunc,
                                                     apply_constraints));
 
   for (unsigned int i=0; i != this->n_qois(); ++i)
     if (qoi_indices.has_index(i))
       this->diff_qoi->finalize_derivative(this->get_adjoint_rhs(i),i);
 }

◆ assemble_residual_derivatives()

void libMesh::ImplicitSystem::assemble_residual_derivatives ( const ParameterVector & parameters )

overridevirtualinherited

Residual parameter derivative function.

Uses finite differences by default.

This will assemble the sensitivity rhs vectors to hold -(partial R / partial p_i), making them ready to solve the forward sensitivity equation.

Can be overridden in derived classes.

Reimplemented from libMesh::System.

Definition at line 654 of file implicit_system.C.

References std::abs(), libMesh::System::add_sensitivity_rhs(), libMesh::ImplicitSystem::assembly(), libMesh::NumericVector< T >::close(), std::max(), libMesh::Real, libMesh::ExplicitSystem::rhs, libMesh::ParameterVector::size(), and libMesh::TOLERANCE.

Referenced by libMesh::ImplicitSystem::adjoint_qoi_parameter_sensitivity(), and libMesh::ImplicitSystem::sensitivity_solve().

 {
   ParameterVector & parameters =
     const_cast<ParameterVector &>(parameters_in);
 
   const unsigned int Np = cast_int<unsigned int>
     (parameters.size());
 
   for (unsigned int p=0; p != Np; ++p)
     {
       NumericVector<Number> & sensitivity_rhs = this->add_sensitivity_rhs(p);
 
       // Approximate -(partial R / partial p) by
       // (R(p-dp) - R(p+dp)) / (2*dp)
 
       Number old_parameter = *parameters[p];
 
       const Real delta_p =
         TOLERANCE * std::max(std::abs(old_parameter), 1e-3);
 
       *parameters[p] -= delta_p;
 
       //      this->assembly(true, false, true);
       this->assembly(true, false, false);
       this->rhs->close();
       sensitivity_rhs = *this->rhs;
 
       *parameters[p] = old_parameter + delta_p;
 
       //      this->assembly(true, false, true);
       this->assembly(true, false, false);
       this->rhs->close();
 
       sensitivity_rhs -= *this->rhs;
       sensitivity_rhs /= (2*delta_p);
       sensitivity_rhs.close();
 
       *parameters[p] = old_parameter;
     }
 }

◆ assembly()

void libMesh::FEMSystem::assembly	(	bool	get_residual,
		bool	get_jacobian,
		bool	apply_heterogeneous_constraints = `false`,
		bool	apply_no_constraints = `false`
	)

overridevirtualinherited

Prepares matrix or rhs for matrix assembly. Users may reimplement this to add pre- or post-assembly code before or after calling FEMSystem::assembly(). Note that in some cases only current_local_solution is used during assembly, and, therefore, if solution has been altered without update() being called, then the user must call update() before calling this function.

Implements libMesh::DifferentiableSystem.

Definition at line 854 of file fem_system.C.

References libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::DenseMatrix< T >::add(), libMesh::FEMSystem::build_context(), libMesh::NumericVector< T >::close(), libMesh::SparseMatrix< T >::close(), libMesh::err, libMesh::FEType::family, libMesh::DiffContext::get_elem_jacobian(), libMesh::DiffContext::get_elem_residual(), libMesh::System::get_mesh(), libMesh::FEMSystem::init_context(), libMesh::NumericVector< T >::l1_norm(), libMesh::SparseMatrix< T >::l1_norm(), libMesh::DenseMatrix< T >::l1_norm(), libMesh::ImplicitSystem::matrix, std::max(), mesh, libMesh::ParallelObject::n_processors(), libMesh::System::n_variable_groups(), libMesh::FEMSystem::numerical_nonlocal_jacobian(), libMesh::out, libMesh::Threads::parallel_for(), libMesh::FEMContext::pre_fe_reinit(), libMesh::BasicOStreamProxy< charT, traits >::precision(), libMesh::DifferentiableSystem::print_jacobian_norms, libMesh::DifferentiableSystem::print_jacobians, libMesh::DifferentiableSystem::print_residual_norms, libMesh::DifferentiableSystem::print_residuals, libMesh::DifferentiableSystem::print_solution_norms, libMesh::DifferentiableSystem::print_solutions, libMesh::ParallelObject::processor_id(), libMesh::Real, libMesh::StoredRange< iterator_type, object_type >::reset(), libMesh::ExplicitSystem::rhs, libMesh::SCALAR, libMesh::DenseVector< T >::size(), libMesh::System::solution, libMesh::DifferentiableSystem::time_solver, libMesh::Variable::type(), libMesh::System::variable_group(), libMesh::FEMSystem::verify_analytic_jacobians, libMesh::DenseMatrix< T >::zero(), libMesh::SparseMatrix< T >::zero(), and libMesh::NumericVector< T >::zero().

Referenced by continuation_solve(), and solve_tangent().

 {
   libmesh_assert(get_residual || get_jacobian);
 
   // Log residual and jacobian and combined performance separately
 #ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
   const char * log_name;
   if (get_residual && get_jacobian)
     log_name = "assembly()";
   else if (get_residual)
     log_name = "assembly(get_residual)";
   else
     log_name = "assembly(get_jacobian)";
 
   LOG_SCOPE(log_name, "FEMSystem");
 #endif
 
   const MeshBase & mesh = this->get_mesh();
 
   if (print_solution_norms)
     {
       this->solution->close();
 
       std::streamsize old_precision = libMesh::out.precision();
       libMesh::out.precision(16);
       libMesh::out << "|U| = "
                    << this->solution->l1_norm()
                    << std::endl;
       libMesh::out.precision(old_precision);
     }
   if (print_solutions)
     {
       std::streamsize old_precision = libMesh::out.precision();
       libMesh::out.precision(16);
       libMesh::out << "U = [" << *(this->solution)
                    << "];" << std::endl;
       libMesh::out.precision(old_precision);
     }
 
   // Is this definitely necessary? [RHS]
   // Yes. [RHS 2012]
   if (get_jacobian)
     matrix->zero();
   if (get_residual)
     rhs->zero();
 
   // Stupid C++ lets you set *Real* verify_analytic_jacobians = true!
   if (verify_analytic_jacobians > 0.5)
     {
       libMesh::err << "WARNING!  verify_analytic_jacobians was set "
                    << "to absurdly large value of "
                    << verify_analytic_jacobians << std::endl;
       libMesh::err << "Resetting to 1e-6!" << std::endl;
       verify_analytic_jacobians = 1e-6;
     }
 
   // In time-dependent problems, the nonlinear function we're trying
   // to solve at each timestep may depend on the particular solver
   // we're using
   libmesh_assert(time_solver.get());
 
   // Build the residual and jacobian contributions on every active
   // mesh element on this processor
   Threads::parallel_for
     (elem_range.reset(mesh.active_local_elements_begin(),
                       mesh.active_local_elements_end()),
      AssemblyContributions(*this, get_residual, get_jacobian,
                            apply_heterogeneous_constraints,
                            apply_no_constraints));
 
   // Check and see if we have SCALAR variables
   bool have_scalar = false;
   for (unsigned int i=0; i != this->n_variable_groups(); ++i)
     {
       if (this->variable_group(i).type().family == SCALAR)
         {
           have_scalar = true;
           break;
         }
     }
 
   // SCALAR dofs are stored on the last processor, so we'll evaluate
   // their equation terms there and only if we have a SCALAR variable
   if (this->processor_id() == (this->n_processors()-1) && have_scalar)
     {
       std::unique_ptr<DiffContext> con = this->build_context();
       FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
       this->init_context(_femcontext);
       _femcontext.pre_fe_reinit(*this, nullptr);
 
       bool jacobian_computed =
         this->time_solver->nonlocal_residual(get_jacobian, _femcontext);
 
       // Nonlocal residuals are likely to be length 0, in which case we
       // don't need to do any more.  And we shouldn't try to do any
       // more; lots of DenseVector/DenseMatrix code assumes rank>0.
       if (_femcontext.get_elem_residual().size())
         {
           // Compute a numeric jacobian if we have to
           if (get_jacobian && !jacobian_computed)
             {
               // Make sure we didn't compute a jacobian and lie about it
               libmesh_assert_equal_to (_femcontext.get_elem_jacobian().l1_norm(), 0.0);
               // Logging of numerical jacobians is done separately
               this->numerical_nonlocal_jacobian(_femcontext);
             }
 
           // Compute a numeric jacobian if we're asked to verify the
           // analytic jacobian we got
           if (get_jacobian && jacobian_computed &&
               this->verify_analytic_jacobians != 0.0)
             {
               DenseMatrix<Number> analytic_jacobian(_femcontext.get_elem_jacobian());
 
               _femcontext.get_elem_jacobian().zero();
               // Logging of numerical jacobians is done separately
               this->numerical_nonlocal_jacobian(_femcontext);
 
               Real analytic_norm = analytic_jacobian.l1_norm();
               Real numerical_norm = _femcontext.get_elem_jacobian().l1_norm();
 
               // If we can continue, we'll probably prefer the analytic jacobian
               analytic_jacobian.swap(_femcontext.get_elem_jacobian());
 
               // The matrix "analytic_jacobian" will now hold the error matrix
               analytic_jacobian.add(-1.0, _femcontext.get_elem_jacobian());
               Real error_norm = analytic_jacobian.l1_norm();
 
               Real relative_error = error_norm /
                 std::max(analytic_norm, numerical_norm);
 
               if (relative_error > this->verify_analytic_jacobians)
                 {
                   libMesh::err << "Relative error " << relative_error
                                << " detected in analytic jacobian on nonlocal dofs!"
                                << std::endl;
 
                   std::streamsize old_precision = libMesh::out.precision();
                   libMesh::out.precision(16);
                   libMesh::out << "J_analytic nonlocal = "
                                << _femcontext.get_elem_jacobian() << std::endl;
                   analytic_jacobian.add(1.0, _femcontext.get_elem_jacobian());
                   libMesh::out << "J_numeric nonlocal = "
                                << analytic_jacobian << std::endl;
 
                   libMesh::out.precision(old_precision);
 
                   libmesh_error_msg("Relative error too large, exiting!");
                 }
             }
 
           add_element_system
             (*this, get_residual, get_jacobian,
              apply_heterogeneous_constraints, apply_no_constraints, _femcontext);
         }
     }
 
   if (get_residual && (print_residual_norms || print_residuals))
     this->rhs->close();
   if (get_residual && print_residual_norms)
     {
       std::streamsize old_precision = libMesh::out.precision();
       libMesh::out.precision(16);
       libMesh::out << "|F| = " << this->rhs->l1_norm() << std::endl;
       libMesh::out.precision(old_precision);
     }
   if (get_residual && print_residuals)
     {
       std::streamsize old_precision = libMesh::out.precision();
       libMesh::out.precision(16);
       libMesh::out << "F = [" << *(this->rhs) << "];" << std::endl;
       libMesh::out.precision(old_precision);
     }
 
   if (get_jacobian && (print_jacobian_norms || print_jacobians))
     this->matrix->close();
   if (get_jacobian && print_jacobian_norms)
     {
       std::streamsize old_precision = libMesh::out.precision();
       libMesh::out.precision(16);
       libMesh::out << "|J| = " << this->matrix->l1_norm() << std::endl;
       libMesh::out.precision(old_precision);
     }
   if (get_jacobian && print_jacobians)
     {
       std::streamsize old_precision = libMesh::out.precision();
       libMesh::out.precision(16);
       libMesh::out << "J = [" << *(this->matrix) << "];" << std::endl;
       libMesh::out.precision(old_precision);
     }
 }

◆ attach_assemble_function()

void libMesh::System::attach_assemble_function ( void fptrEquationSystems &es, const std::string &name )

inherited

Register a user function to use in assembling the system matrix and RHS.

Definition at line 1777 of file system.C.

References libMesh::System::_assemble_system_function, libMesh::System::_assemble_system_object, and libMesh::out.

 {
   libmesh_assert(fptr);
 
   if (_assemble_system_object != nullptr)
     {
       libmesh_here();
       libMesh::out << "WARNING:  Cannot specify both assembly function and object!"
                    << std::endl;
 
       _assemble_system_object = nullptr;
     }
 
   _assemble_system_function = fptr;
 }

◆ attach_assemble_object()

void libMesh::System::attach_assemble_object ( System::Assembly & assemble_in )

inherited

Register a user object to use in assembling the system matrix and RHS.

Definition at line 1796 of file system.C.

References libMesh::System::_assemble_system_function, libMesh::System::_assemble_system_object, and libMesh::out.

 {
   if (_assemble_system_function != nullptr)
     {
       libmesh_here();
       libMesh::out << "WARNING:  Cannot specify both assembly object and function!"
                    << std::endl;
 
       _assemble_system_function = nullptr;
     }
 
   _assemble_system_object = &assemble_in;
 }

◆ attach_constraint_function()

void libMesh::System::attach_constraint_function ( void fptrEquationSystems &es, const std::string &name )

inherited

Register a user function for imposing constraints.

Definition at line 1812 of file system.C.

References libMesh::System::_constrain_system_function, libMesh::System::_constrain_system_object, and libMesh::out.

 {
   libmesh_assert(fptr);
 
   if (_constrain_system_object != nullptr)
     {
       libmesh_here();
       libMesh::out << "WARNING:  Cannot specify both constraint function and object!"
                    << std::endl;
 
       _constrain_system_object = nullptr;
     }
 
   _constrain_system_function = fptr;
 }

◆ attach_constraint_object()

void libMesh::System::attach_constraint_object ( System::Constraint & constrain )

inherited

Register a user object for imposing constraints.

Definition at line 1831 of file system.C.

References libMesh::System::_constrain_system_function, libMesh::System::_constrain_system_object, and libMesh::out.

 {
   if (_constrain_system_function != nullptr)
     {
       libmesh_here();
       libMesh::out << "WARNING:  Cannot specify both constraint object and function!"
                    << std::endl;
 
       _constrain_system_function = nullptr;
     }
 
   _constrain_system_object = &constrain;
 }

◆ attach_init_function()

void libMesh::System::attach_init_function ( void fptrEquationSystems &es, const std::string &name )

inherited

Register a user function to use in initializing the system.

Definition at line 1742 of file system.C.

References libMesh::System::_init_system_function, libMesh::System::_init_system_object, and libMesh::out.

 {
   libmesh_assert(fptr);
 
   if (_init_system_object != nullptr)
     {
       libmesh_here();
       libMesh::out << "WARNING:  Cannot specify both initialization function and object!"
                    << std::endl;
 
       _init_system_object = nullptr;
     }
 
   _init_system_function = fptr;
 }

◆ attach_init_object()

void libMesh::System::attach_init_object ( System::Initialization & init_in )

inherited

Register a user class to use to initialize the system.

Note: This is exclusive with the attach_init_function.

Definition at line 1761 of file system.C.

References libMesh::System::_init_system_function, libMesh::System::_init_system_object, and libMesh::out.

 {
   if (_init_system_function != nullptr)
     {
       libmesh_here();
       libMesh::out << "WARNING:  Cannot specify both initialization object and function!"
                    << std::endl;
 
       _init_system_function = nullptr;
     }
 
   _init_system_object = &init_in;
 }

◆ attach_physics()

void libMesh::DifferentiableSystem::attach_physics ( DifferentiablePhysics * physics_in )

inlineinherited

Attach external Physics object.

Definition at line 197 of file diff_system.h.

References libMesh::DifferentiableSystem::_diff_physics, libMesh::DifferentiablePhysics::clone_physics(), and libMesh::DifferentiablePhysics::init_physics().

198 { this->_diff_physics = (physics_in->clone_physics()).release();

199 this->_diff_physics->init_physics(*this);}

libMesh::DifferentiableSystem::_diff_physics

DifferentiablePhysics * _diff_physics

Definition: diff_system.h:371

libMesh::DifferentiablePhysics::init_physics

virtual void init_physics(const System &sys)

◆ attach_qoi()

void libMesh::DifferentiableSystem::attach_qoi ( DifferentiableQoI * qoi_in )

inlineinherited

Attach external QoI object.

Definition at line 225 of file diff_system.h.

References libMesh::DifferentiableQoI::clone(), libMesh::DifferentiableSystem::diff_qoi, libMesh::DifferentiableQoI::init_qoi(), and libMesh::System::qoi.

   { this->diff_qoi = (qoi_in->clone()).release();
     // User needs to resize qoi system qoi accordingly
     this->diff_qoi->init_qoi( this->qoi );}

◆ attach_QOI_derivative()

void libMesh::System::attach_QOI_derivative ( void fptrEquationSystems &es, const std::string &name, const QoISet &qoi_indices, bool include_liftfunc, bool apply_constraints )

inherited

Register a user function for evaluating derivatives of a quantity of interest with respect to test functions, whose values should be placed in System::rhs

Definition at line 1883 of file system.C.

References libMesh::System::_qoi_evaluate_derivative_function, libMesh::System::_qoi_evaluate_derivative_object, and libMesh::out.

 {
   libmesh_assert(fptr);
 
   if (_qoi_evaluate_derivative_object != nullptr)
     {
       libmesh_here();
       libMesh::out << "WARNING:  Cannot specify both QOI derivative function and object!"
                    << std::endl;
 
       _qoi_evaluate_derivative_object = nullptr;
     }
 
   _qoi_evaluate_derivative_function = fptr;
 }

◆ attach_QOI_derivative_object()

void libMesh::System::attach_QOI_derivative_object ( QOIDerivative & qoi_derivative )

inherited

Register a user object for evaluating derivatives of a quantity of interest with respect to test functions, whose values should be placed in System::rhs

Definition at line 1902 of file system.C.

References libMesh::System::_qoi_evaluate_derivative_function, libMesh::System::_qoi_evaluate_derivative_object, and libMesh::out.

 {
   if (_qoi_evaluate_derivative_function != nullptr)
     {
       libmesh_here();
       libMesh::out << "WARNING:  Cannot specify both QOI derivative object and function!"
                    << std::endl;
 
       _qoi_evaluate_derivative_function = nullptr;
     }
 
   _qoi_evaluate_derivative_object = &qoi_derivative;
 }

◆ attach_QOI_function()

void libMesh::System::attach_QOI_function ( void fptrEquationSystems &es, const std::string &name, const QoISet &qoi_indices )

inherited

Register a user function for evaluating the quantities of interest, whose values should be placed in System::qoi

Definition at line 1847 of file system.C.

References libMesh::System::_qoi_evaluate_function, libMesh::System::_qoi_evaluate_object, and libMesh::out.

 {
   libmesh_assert(fptr);
 
   if (_qoi_evaluate_object != nullptr)
     {
       libmesh_here();
       libMesh::out << "WARNING:  Cannot specify both QOI function and object!"
                    << std::endl;
 
       _qoi_evaluate_object = nullptr;
     }
 
   _qoi_evaluate_function = fptr;
 }

◆ attach_QOI_object()

void libMesh::System::attach_QOI_object ( QOI & qoi )

inherited

Register a user object for evaluating the quantities of interest, whose values should be placed in System::qoi

Definition at line 1867 of file system.C.

References libMesh::System::_qoi_evaluate_function, libMesh::System::_qoi_evaluate_object, and libMesh::out.

 {
   if (_qoi_evaluate_function != nullptr)
     {
       libmesh_here();
       libMesh::out << "WARNING:  Cannot specify both QOI object and function!"
                    << std::endl;
 
       _qoi_evaluate_function = nullptr;
     }
 
   _qoi_evaluate_object = &qoi_in;
 }

◆ boundary_project_solution() [1/2]

void libMesh::System::boundary_project_solution	(	const std::set< boundary_id_type > &	b,
		const std::vector< unsigned int > &	variables,
		FunctionBase< Number > *	f,
		FunctionBase< Gradient > *	g = `nullptr`
	)

inherited

Projects arbitrary boundary functions onto a vector of degree of freedom values for the current system. Only degrees of freedom which affect the function's trace on a boundary in the set b are affected. Only degrees of freedom associated with the variables listed in the vector variables are projected. The function value f and its gradient g are user-provided cloneable functors. A gradient g is only required/used for projecting onto finite element spaces with continuous derivatives. If non-default Parameters are to be used, they can be provided in the parameters argument.

This method projects an arbitrary boundary function onto the solution via L2 projections and nodal interpolations on each element.

Definition at line 985 of file system_projection.C.

 {
   this->boundary_project_vector(b, variables, *solution, f, g);
 
   solution->localize(*current_local_solution);
 }

◆ boundary_project_solution() [2/2]

void libMesh::System::boundary_project_solution	(	const std::set< boundary_id_type > &	b,
		const std::vector< unsigned int > &	variables,
		ValueFunctionPointer	fptr,
		GradientFunctionPointer	gptr,
		const Parameters &	parameters
	)

inherited

Projects arbitrary boundary functions onto a vector of degree of freedom values for the current system. Only degrees of freedom which affect the function's trace on a boundary in the set b are affected. Only degrees of freedom associated with the variables listed in the vector variables are projected. The function value fptr and its gradient gptr are represented by function pointers. A gradient gptr is only required/used for projecting onto finite element spaces with continuous derivatives.

This method projects components of an arbitrary boundary function onto the solution via L2 projections and nodal interpolations on each element.

Definition at line 968 of file system_projection.C.

 {
   WrappedFunction<Number> f(*this, fptr, &parameters);
   WrappedFunction<Gradient> g(*this, gptr, &parameters);
   this->boundary_project_solution(b, variables, &f, &g);
 }

◆ boundary_project_vector() [1/2]

void libMesh::System::boundary_project_vector	(	const std::set< boundary_id_type > &	b,
		const std::vector< unsigned int > &	variables,
		NumericVector< Number > &	new_vector,
		FunctionBase< Number > *	f,
		FunctionBase< Gradient > *	g = `nullptr`,
		int	is_adjoint = `-1`
	)		const

inherited

Projects arbitrary boundary functions onto a vector of degree of freedom values for the current system. Only degrees of freedom which affect the function's trace on a boundary in the set b are affected. Only degrees of freedom associated with the variables listed in the vector variables are projected. The function value f and its gradient g are user-provided cloneable functors. A gradient g is only required/used for projecting onto finite element spaces with continuous derivatives. If non-default Parameters are to be used, they can be provided in the parameters argument.

Constrain the new vector using the requested adjoint rather than primal constraints if is_adjoint is non-negative.

This method projects an arbitrary function via L2 projections and nodal interpolations on each element.

Definition at line 1021 of file system_projection.C.

References libMesh::NumericVector< T >::close(), and libMesh::Threads::parallel_for().

 {
   LOG_SCOPE ("boundary_project_vector()", "System");
 
   Threads::parallel_for
     (ConstElemRange (this->get_mesh().active_local_elements_begin(),
                      this->get_mesh().active_local_elements_end() ),
      BoundaryProjectSolution(b, variables, *this, f, g,
                              this->get_equation_systems().parameters,
                              new_vector)
      );
 
   // We don't do SCALAR dofs when just projecting the boundary, so
   // we're done here.
 
   new_vector.close();
 
 #ifdef LIBMESH_ENABLE_CONSTRAINTS
   if (is_adjoint == -1)
     this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
   else if (is_adjoint >= 0)
     this->get_dof_map().enforce_adjoint_constraints_exactly(new_vector,
                                                             is_adjoint);
 #endif
 }

◆ boundary_project_vector() [2/2]

void libMesh::System::boundary_project_vector	(	const std::set< boundary_id_type > &	b,
		const std::vector< unsigned int > &	variables,
		ValueFunctionPointer	fptr,
		GradientFunctionPointer	gptr,
		const Parameters &	parameters,
		NumericVector< Number > &	new_vector,
		int	is_adjoint = `-1`
	)		const

inherited

Projects arbitrary boundary functions onto a vector of degree of freedom values for the current system. Only degrees of freedom which affect the function's trace on a boundary in the set b are affected. Only degrees of freedom associated with the variables listed in the vector variables are projected. The function value fptr and its gradient gptr are represented by function pointers. A gradient gptr is only required/used for projecting onto finite element spaces with continuous derivatives.

Constrain the new vector using the requested adjoint rather than primal constraints if is_adjoint is non-negative.

This method projects an arbitrary boundary function via L2 projections and nodal interpolations on each element.

Definition at line 1003 of file system_projection.C.

 {
   WrappedFunction<Number> f(*this, fptr, &parameters);
   WrappedFunction<Gradient> g(*this, gptr, &parameters);
   this->boundary_project_vector(b, variables, new_vector, &f, &g,
                                 is_adjoint);
 }

◆ build_context()

std::unique_ptr< DiffContext > libMesh::FEMSystem::build_context ( )

overridevirtualinherited

Builds a FEMContext object with enough information to do evaluations on each element.

For most problems, the default FEMSystem implementation is correct; users who subclass FEMContext will need to also reimplement this method to build it.

Reimplemented from libMesh::DifferentiableSystem.

Definition at line 1318 of file fem_system.C.

References libMesh::DifferentiableSystem::deltat, libMesh::DifferentiablePhysics::get_mesh_system(), libMesh::DifferentiablePhysics::get_mesh_x_var(), libMesh::DifferentiablePhysics::get_mesh_y_var(), libMesh::DifferentiablePhysics::get_mesh_z_var(), libMesh::DifferentiableSystem::get_physics(), libMesh::DifferentiableSystem::get_time_solver(), libMesh::TimeSolver::is_adjoint(), libMesh::DiffContext::is_adjoint(), libMesh::DiffContext::set_deltat_pointer(), libMesh::FEMContext::set_mesh_system(), libMesh::FEMContext::set_mesh_x_var(), libMesh::FEMContext::set_mesh_y_var(), and libMesh::FEMContext::set_mesh_z_var().

Referenced by libMesh::FEMSystem::assembly(), libMesh::FEMSystem::mesh_position_get(), and libMesh::FEMSystem::mesh_position_set().

 {
   FEMContext * fc = new FEMContext(*this);
 
   DifferentiablePhysics * phys = this->get_physics();
 
   libmesh_assert (phys);
 
   // If we are solving a moving mesh problem, tell that to the Context
   fc->set_mesh_system(phys->get_mesh_system());
   fc->set_mesh_x_var(phys->get_mesh_x_var());
   fc->set_mesh_y_var(phys->get_mesh_y_var());
   fc->set_mesh_z_var(phys->get_mesh_z_var());
 
   fc->set_deltat_pointer( &deltat );
 
   // If we are solving the adjoint problem, tell that to the Context
   fc->is_adjoint() = this->get_time_solver().is_adjoint();
 
   return std::unique_ptr<DiffContext>(fc);
 }

◆ calculate_norm() [1/2]

Real libMesh::System::calculate_norm	(	const NumericVector< Number > &	v,
		unsigned int	var,
		FEMNormType	norm_type,
		std::set< unsigned int > *	skip_dimensions = `nullptr`
	)		const

inherited

Returns: A norm of variable var in the vector v, in the specified norm (e.g. L2, L_INF, H1)

Definition at line 1378 of file system.C.

References libMesh::DISCRETE_L1, libMesh::DISCRETE_L2, libMesh::DISCRETE_L_INF, libMesh::System::discrete_var_norm(), libMesh::L2, libMesh::System::n_vars(), and libMesh::Real.

Referenced by libMesh::AdaptiveTimeSolver::calculate_norm(), and libMesh::UnsteadySolver::du().

 {
   //short circuit to save time
   if (norm_type == DISCRETE_L1 ||
       norm_type == DISCRETE_L2 ||
       norm_type == DISCRETE_L_INF)
     return discrete_var_norm(v,var,norm_type);
 
   // Not a discrete norm
   std::vector<FEMNormType> norms(this->n_vars(), L2);
   std::vector<Real> weights(this->n_vars(), 0.0);
   norms[var] = norm_type;
   weights[var] = 1.0;
   Real val = this->calculate_norm(v, SystemNorm(norms, weights), skip_dimensions);
   return val;
 }

◆ calculate_norm() [2/2]

Real libMesh::System::calculate_norm	(	const NumericVector< Number > &	v,
		const SystemNorm &	norm,
		std::set< unsigned int > *	skip_dimensions = `nullptr`
	)		const

inherited

Returns: A norm of the vector v, using component_norm and component_scale to choose and weight the norms of each variable.

Definition at line 1400 of file system.C.

References libMesh::System::_dof_map, libMesh::System::_mesh, std::abs(), libMesh::TypeVector< T >::add_scaled(), libMesh::TypeTensor< T >::add_scaled(), libMesh::FEGenericBase< OutputType >::build(), libMesh::NumericVector< T >::build(), libMesh::ParallelObject::comm(), libMesh::FEType::default_quadrature_rule(), libMesh::DISCRETE_L1, libMesh::DISCRETE_L2, libMesh::DISCRETE_L_INF, libMesh::System::discrete_var_norm(), libMesh::DofMap::dof_indices(), libMesh::MeshBase::elem_dimensions(), libMesh::FEGenericBase< OutputType >::get_d2phi(), libMesh::System::get_dof_map(), libMesh::FEGenericBase< OutputType >::get_dphi(), libMesh::FEAbstract::get_JxW(), libMesh::System::get_mesh(), libMesh::FEGenericBase< OutputType >::get_phi(), libMesh::H1, libMesh::H1_SEMINORM, libMesh::H2, libMesh::H2_SEMINORM, libMesh::SystemNorm::is_discrete(), libMesh::L1, libMesh::NumericVector< T >::l1_norm(), libMesh::L2, libMesh::NumericVector< T >::l2_norm(), libMesh::L_INF, libMesh::NumericVector< T >::linfty_norm(), libMesh::NumericVector< T >::localize(), std::max(), libMesh::Parallel::Communicator::max(), libMesh::QBase::n_points(), libMesh::System::n_vars(), libMesh::TypeVector< T >::norm(), libMesh::TypeTensor< T >::norm(), libMesh::TensorTools::norm_sq(), libMesh::TypeVector< T >::norm_sq(), libMesh::TypeTensor< T >::norm_sq(), libMesh::Real, libMesh::FEAbstract::reinit(), libMesh::SERIAL, libMesh::NumericVector< T >::size(), libMesh::Parallel::Communicator::sum(), libMesh::SystemNorm::type(), libMesh::DofMap::variable_type(), libMesh::W1_INF_SEMINORM, libMesh::W2_INF_SEMINORM, libMesh::SystemNorm::weight(), and libMesh::SystemNorm::weight_sq().

 {
   // This function must be run on all processors at once
   parallel_object_only();
 
   LOG_SCOPE ("calculate_norm()", "System");
 
   // Zero the norm before summation
   Real v_norm = 0.;
 
   if (norm.is_discrete())
     {
       //Check to see if all weights are 1.0 and all types are equal
       FEMNormType norm_type0 = norm.type(0);
       unsigned int check_var = 0;
       for (; check_var != this->n_vars(); ++check_var)
         if ((norm.weight(check_var) != 1.0) || (norm.type(check_var) != norm_type0))
           break;
 
       //All weights were 1.0 so just do the full vector discrete norm
       if (check_var == this->n_vars())
         {
           if (norm_type0 == DISCRETE_L1)
             return v.l1_norm();
           if (norm_type0 == DISCRETE_L2)
             return v.l2_norm();
           if (norm_type0 == DISCRETE_L_INF)
             return v.linfty_norm();
           else
             libmesh_error_msg("Invalid norm_type0 = " << norm_type0);
         }
 
       for (unsigned int var=0; var != this->n_vars(); ++var)
         {
           // Skip any variables we don't need to integrate
           if (norm.weight(var) == 0.0)
             continue;
 
           v_norm += norm.weight(var) * discrete_var_norm(v, var, norm.type(var));
         }
 
       return v_norm;
     }
 
   // Localize the potentially parallel vector
   std::unique_ptr<NumericVector<Number>> local_v = NumericVector<Number>::build(this->comm());
   local_v->init(v.size(), true, SERIAL);
   v.localize (*local_v, _dof_map->get_send_list());
 
   // I'm not sure how best to mix Hilbert norms on some variables (for
   // which we'll want to square then sum then square root) with norms
   // like L_inf (for which we'll just want to take an absolute value
   // and then sum).
   bool using_hilbert_norm = true,
     using_nonhilbert_norm = true;
 
   // Loop over all variables
   for (unsigned int var=0; var != this->n_vars(); ++var)
     {
       // Skip any variables we don't need to integrate
       Real norm_weight_sq = norm.weight_sq(var);
       if (norm_weight_sq == 0.0)
         continue;
       Real norm_weight = norm.weight(var);
 
       // Check for unimplemented norms (rather than just returning 0).
       FEMNormType norm_type = norm.type(var);
       if ((norm_type==H1) ||
           (norm_type==H2) ||
           (norm_type==L2) ||
           (norm_type==H1_SEMINORM) ||
           (norm_type==H2_SEMINORM))
         {
           if (!using_hilbert_norm)
             libmesh_not_implemented();
           using_nonhilbert_norm = false;
         }
       else if ((norm_type==L1) ||
                (norm_type==L_INF) ||
                (norm_type==W1_INF_SEMINORM) ||
                (norm_type==W2_INF_SEMINORM))
         {
           if (!using_nonhilbert_norm)
             libmesh_not_implemented();
           using_hilbert_norm = false;
         }
       else
         libmesh_not_implemented();
 
       const FEType & fe_type = this->get_dof_map().variable_type(var);
 
       // Allow space for dims 0-3, even if we don't use them all
       std::vector<std::unique_ptr<FEBase>> fe_ptrs(4);
       std::vector<std::unique_ptr<QBase>> q_rules(4);
 
       const std::set<unsigned char> & elem_dims = _mesh.elem_dimensions();
 
       // Prepare finite elements for each dimension present in the mesh
       for (const auto & dim : elem_dims)
         {
           if (skip_dimensions && skip_dimensions->find(dim) != skip_dimensions->end())
             continue;
 
           // Construct quadrature and finite element objects
           q_rules[dim] = fe_type.default_quadrature_rule (dim);
           fe_ptrs[dim] = FEBase::build(dim, fe_type);
 
           // Attach quadrature rule to FE object
           fe_ptrs[dim]->attach_quadrature_rule (q_rules[dim].get());
         }
 
       std::vector<dof_id_type> dof_indices;
 
       // Begin the loop over the elements
       for (const auto & elem : this->get_mesh().active_local_element_ptr_range())
         {
           const unsigned int dim = elem->dim();
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
           // One way for implementing this would be to exchange the fe with the FEInterface- class.
           // However, it needs to be discussed whether integral-norms make sense for infinite elements.
           // or in which sense they could make sense.
           if (elem->infinite() )
             libmesh_not_implemented();
 
 #endif
 
           if (skip_dimensions && skip_dimensions->find(dim) != skip_dimensions->end())
             continue;
 
           FEBase * fe = fe_ptrs[dim].get();
           QBase * qrule = q_rules[dim].get();
           libmesh_assert(fe);
           libmesh_assert(qrule);
 
           const std::vector<Real> &               JxW = fe->get_JxW();
           const std::vector<std::vector<Real>> * phi = nullptr;
           if (norm_type == H1 ||
               norm_type == H2 ||
               norm_type == L2 ||
               norm_type == L1 ||
               norm_type == L_INF)
             phi = &(fe->get_phi());
 
           const std::vector<std::vector<RealGradient>> * dphi = nullptr;
           if (norm_type == H1 ||
               norm_type == H2 ||
               norm_type == H1_SEMINORM ||
               norm_type == W1_INF_SEMINORM)
             dphi = &(fe->get_dphi());
 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
           const std::vector<std::vector<RealTensor>> *   d2phi = nullptr;
           if (norm_type == H2 ||
               norm_type == H2_SEMINORM ||
               norm_type == W2_INF_SEMINORM)
             d2phi = &(fe->get_d2phi());
 #endif
 
           fe->reinit (elem);
 
           this->get_dof_map().dof_indices (elem, dof_indices, var);
 
           const unsigned int n_qp = qrule->n_points();
 
           const unsigned int n_sf = cast_int<unsigned int>
             (dof_indices.size());
 
           // Begin the loop over the Quadrature points.
           for (unsigned int qp=0; qp<n_qp; qp++)
             {
               if (norm_type == L1)
                 {
                   Number u_h = 0.;
                   for (unsigned int i=0; i != n_sf; ++i)
                     u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
                   v_norm += norm_weight *
                     JxW[qp] * std::abs(u_h);
                 }
 
               if (norm_type == L_INF)
                 {
                   Number u_h = 0.;
                   for (unsigned int i=0; i != n_sf; ++i)
                     u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
                   v_norm = std::max(v_norm, norm_weight * std::abs(u_h));
                 }
 
               if (norm_type == H1 ||
                   norm_type == H2 ||
                   norm_type == L2)
                 {
                   Number u_h = 0.;
                   for (unsigned int i=0; i != n_sf; ++i)
                     u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
                   v_norm += norm_weight_sq *
                     JxW[qp] * TensorTools::norm_sq(u_h);
                 }
 
               if (norm_type == H1 ||
                   norm_type == H2 ||
                   norm_type == H1_SEMINORM)
                 {
                   Gradient grad_u_h;
                   for (unsigned int i=0; i != n_sf; ++i)
                     grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
                   v_norm += norm_weight_sq *
                     JxW[qp] * grad_u_h.norm_sq();
                 }
 
               if (norm_type == W1_INF_SEMINORM)
                 {
                   Gradient grad_u_h;
                   for (unsigned int i=0; i != n_sf; ++i)
                     grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
                   v_norm = std::max(v_norm, norm_weight * grad_u_h.norm());
                 }
 
 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
               if (norm_type == H2 ||
                   norm_type == H2_SEMINORM)
                 {
                   Tensor hess_u_h;
                   for (unsigned int i=0; i != n_sf; ++i)
                     hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
                   v_norm += norm_weight_sq *
                     JxW[qp] * hess_u_h.norm_sq();
                 }
 
               if (norm_type == W2_INF_SEMINORM)
                 {
                   Tensor hess_u_h;
                   for (unsigned int i=0; i != n_sf; ++i)
                     hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
                   v_norm = std::max(v_norm, norm_weight * hess_u_h.norm());
                 }
 #endif
             }
         }
     }
 
   if (using_hilbert_norm)
     {
       this->comm().sum(v_norm);
       v_norm = std::sqrt(v_norm);
     }
   else
     {
       this->comm().max(v_norm);
     }
 
   return v_norm;
 }

◆ clear()

void libMesh::ContinuationSystem::clear ( )

overridevirtual

Clear all the data structures associated with the system.

Reimplemented from libMesh::DifferentiableSystem.

Definition at line 80 of file continuation_system.C.

References libMesh::DifferentiableSystem::clear().

Referenced by ~ContinuationSystem().

 {
   // FIXME: Do anything here, e.g. zero vectors, etc?
 
   // Call the Parent's clear function
   Parent::clear();
 }

◆ clear_physics()

virtual void libMesh::DifferentiablePhysics::clear_physics ( )

virtualinherited

Clear any data structures associated with the physics.

Referenced by libMesh::DifferentiableSystem::clear().

◆ clear_qoi()

virtual void libMesh::DifferentiableQoI::clear_qoi ( )

inlinevirtualinherited

Clear all the data structures associated with the QoI.

Definition at line 75 of file diff_qoi.h.

Referenced by libMesh::DifferentiableSystem::clear().

75 {}

◆ clone()

virtual std::unique_ptr<DifferentiableQoI> libMesh::DifferentiableSystem::clone ( )

inlineoverridevirtualinherited

We don't allow systems to be attached to each other

Implements libMesh::DifferentiableQoI.

Definition at line 169 of file diff_system.h.

Referenced by libMesh::AdjointRefinementEstimator::estimate_error().

   {
     libmesh_not_implemented();
     // dummy
     return std::unique_ptr<DifferentiableQoI>(this);
   }

◆ clone_physics()

virtual std::unique_ptr<DifferentiablePhysics> libMesh::DifferentiableSystem::clone_physics ( )

inlineoverridevirtualinherited

We don't allow systems to be attached to each other

Implements libMesh::DifferentiablePhysics.

Definition at line 159 of file diff_system.h.

   {
     libmesh_not_implemented();
     // dummy
     return std::unique_ptr<DifferentiablePhysics>(this);
   }

◆ comm()

const Parallel::Communicator& libMesh::ParallelObject::comm ( ) const

inlineinherited

Returns: A reference to the Parallel::Communicator object used by this mesh.

Definition at line 89 of file parallel_object.h.

References libMesh::ParallelObject::_communicator.

90 { return _communicator; }

libMesh::ParallelObject::_communicator

const Parallel::Communicator & _communicator

Definition: parallel_object.h:107

◆ compare()

bool libMesh::System::compare	(	const System &	other_system,
		const Real	threshold,
		const bool	verbose
	)		const

virtualinherited

Returns: true when the other system contains identical data, up to the given threshold. Outputs some diagnostic info when verbose is set.

Definition at line 514 of file system.C.

References libMesh::System::_is_initialized, libMesh::System::_sys_name, libMesh::System::_vectors, libMesh::System::get_vector(), libMesh::System::n_vectors(), libMesh::System::name(), libMesh::out, and libMesh::System::solution.

Referenced by libMesh::EquationSystems::compare().

 {
   // we do not care for matrices, but for vectors
   libmesh_assert (_is_initialized);
   libmesh_assert (other_system._is_initialized);
 
   if (verbose)
     {
       libMesh::out << "  Systems \"" << _sys_name << "\"" << std::endl;
       libMesh::out << "   comparing matrices not supported." << std::endl;
       libMesh::out << "   comparing names...";
     }
 
   // compare the name: 0 means identical
   const int name_result = _sys_name.compare(other_system.name());
   if (verbose)
     {
       if (name_result == 0)
         libMesh::out << " identical." << std::endl;
       else
         libMesh::out << "  names not identical." << std::endl;
       libMesh::out << "   comparing solution vector...";
     }
 
 
   // compare the solution: -1 means identical
   const int solu_result = solution->compare (*other_system.solution.get(),
                                              threshold);
 
   if (verbose)
     {
       if (solu_result == -1)
         libMesh::out << " identical up to threshold." << std::endl;
       else
         libMesh::out << "  first difference occurred at index = "
                      << solu_result << "." << std::endl;
     }
 
 
   // safety check, whether we handle at least the same number
   // of vectors
   std::vector<int> ov_result;
 
   if (this->n_vectors() != other_system.n_vectors())
     {
       if (verbose)
         {
           libMesh::out << "   Fatal difference. This system handles "
                        << this->n_vectors() << " add'l vectors," << std::endl
                        << "   while the other system handles "
                        << other_system.n_vectors()
                        << " add'l vectors." << std::endl
                        << "   Aborting comparison." << std::endl;
         }
       return false;
     }
   else if (this->n_vectors() == 0)
     {
       // there are no additional vectors...
       ov_result.clear ();
     }
   else
     {
       // compare other vectors
       for (auto & pr : _vectors)
         {
           if (verbose)
             libMesh::out << "   comparing vector \""
                          << pr.first << "\" ...";
 
           // assume they have the same name
           const NumericVector<Number> & other_system_vector =
             other_system.get_vector(pr.first);
 
           ov_result.push_back(pr.second->compare (other_system_vector,
                                                   threshold));
 
           if (verbose)
             {
               if (ov_result[ov_result.size()-1] == -1)
                 libMesh::out << " identical up to threshold." << std::endl;
               else
                 libMesh::out << " first difference occurred at" << std::endl
                              << "   index = " << ov_result[ov_result.size()-1] << "." << std::endl;
             }
         }
     } // finished comparing additional vectors
 
 
   bool overall_result;
 
   // sum up the results
   if ((name_result==0) && (solu_result==-1))
     {
       if (ov_result.size()==0)
         overall_result = true;
       else
         {
           bool ov_identical;
           unsigned int n    = 0;
           do
             {
               ov_identical = (ov_result[n]==-1);
               n++;
             }
           while (ov_identical && n<ov_result.size());
           overall_result = ov_identical;
         }
     }
   else
     overall_result = false;
 
   if (verbose)
     {
       libMesh::out << "   finished comparisons, ";
       if (overall_result)
         libMesh::out << "found no differences." << std::endl << std::endl;
       else
         libMesh::out << "found differences." << std::endl << std::endl;
     }
 
   return overall_result;
 }

◆ continuation_solve()

void libMesh::ContinuationSystem::continuation_solve ( )

Perform a continuation solve of the system. In general, you can only begin the continuation solves after either reading in or solving for two previous values of the control parameter. The prior two solutions are required for starting up the continuation method.

Definition at line 358 of file continuation_system.C.

References std::abs(), libMesh::NumericVector< T >::add(), apply_predictor(), libMesh::FEMSystem::assembly(), libMesh::NumericVector< T >::close(), libMesh::SparseMatrix< T >::close(), continuation_parameter, continuation_parameter_tolerance, delta_u, dlambda_ds, libMesh::NumericVector< T >::dot(), ds_current, du_ds, G_Lambda, initial_newton_tolerance, initialize_tangent(), libMesh::NumericVector< T >::l2_norm(), libMesh::libmesh_real(), linear_solver, libMesh::ImplicitSystem::matrix, max_continuation_parameter, libMesh::DiffSolver::max_linear_iterations, libMesh::DiffSolver::max_nonlinear_iterations, std::min(), min_continuation_parameter, n_arclength_reductions, n_backtrack_steps, newton_progress_check, newton_solver, newton_step, old_continuation_parameter, libMesh::out, std::pow(), libMesh::BasicOStreamProxy< charT, traits >::precision(), previous_u, quiet, libMesh::Real, Residual, libMesh::ExplicitSystem::rhs, rhs_mode, libMesh::NumericVector< T >::scale(), libMesh::BasicOStreamProxy< charT, traits >::setf(), libMesh::System::solution, solution_tolerance, tangent_initialized, Theta, Theta_LOCA, libMesh::DifferentiableSystem::time_solver, libMesh::BasicOStreamProxy< charT, traits >::unsetf(), y, y_old, z, and libMesh::NumericVector< T >::zero().

 {
   // Be sure the user has set the continuation parameter pointer
   if (!continuation_parameter)
     libmesh_error_msg("You must set the continuation_parameter pointer " \
                       << "to a member variable of the derived class, preferably in the " \
                       << "Derived class's init_data function.  This is how the ContinuationSystem " \
                       << "updates the continuation parameter.");
 
   // Use extra precision for all the numbers printed in this function.
   std::streamsize old_precision = libMesh::out.precision();
   libMesh::out.precision(16);
   libMesh::out.setf(std::ios_base::scientific);
 
   // We can't start solving the augmented PDE system unless the tangent
   // vectors have been initialized.  This only needs to occur once.
   if (!tangent_initialized)
     initialize_tangent();
 
   // Save the old value of -du/dlambda.  This will be used after the Newton iterations
   // to compute the angle between previous tangent vectors.  This cosine of this angle is
   //
   // tau := abs( (du/d(lambda)_i , du/d(lambda)_{i-1}) / (||du/d(lambda)_i|| * ||du/d(lambda)_{i-1}||) )
   //
   // The scaling factor tau (which should vary between 0 and 1) is used to shrink the step-size ds
   // when we are approaching a turning point.  Note that it can only shrink the step size.
   *y_old = *y;
 
   // Set pointer to underlying Newton solver
   if (!newton_solver)
     newton_solver = cast_ptr<NewtonSolver *> (this->time_solver->diff_solver().get());
 
   // A pair for catching return values from linear system solves.
   std::pair<unsigned int, Real> rval;
 
   // Convergence flag for the entire arcstep
   bool arcstep_converged = false;
 
   // Begin loop over arcstep reductions.
   for (unsigned int ns=0; ns<n_arclength_reductions; ++ns)
     {
       if (!quiet)
         {
           libMesh::out << "Current arclength stepsize, ds_current=" << ds_current << std::endl;
           libMesh::out << "Current parameter value, lambda=" << *continuation_parameter << std::endl;
         }
 
       // Upon exit from the nonlinear loop, the newton_converged flag
       // will tell us the convergence status of Newton's method.
       bool newton_converged = false;
 
       // The nonlinear residual before *any* nonlinear steps have been taken.
       Real nonlinear_residual_firststep = 0.;
 
       // The nonlinear residual from the current "k" Newton step, before the Newton step
       Real nonlinear_residual_beforestep = 0.;
 
       // The nonlinear residual from the current "k" Newton step, after the Newton step
       Real nonlinear_residual_afterstep = 0.;
 
       // The linear solver tolerance, can be updated dynamically at each Newton step.
       Real current_linear_tolerance = 0.;
 
       // The nonlinear loop
       for (newton_step=0; newton_step<newton_solver->max_nonlinear_iterations; ++newton_step)
         {
           libMesh::out << "\n === Starting Newton step " << newton_step << " ===" << std::endl;
 
           // Set the linear system solver tolerance
           //   // 1.) Set the current linear tolerance based as a multiple of the current residual of the system.
           //   const Real residual_multiple = 1.e-4;
           //   Real current_linear_tolerance = residual_multiple*nonlinear_residual_beforestep;
 
           //   // But if the current residual isn't small, don't let the solver exit with zero iterations!
           //   if (current_linear_tolerance > 1.)
           //     current_linear_tolerance = residual_multiple;
 
           // 2.) Set the current linear tolerance based on the method based on technique of Eisenstat & Walker.
           if (newton_step==0)
             {
               // At first step, only try reducing the residual by a small amount
               current_linear_tolerance = initial_newton_tolerance;//0.01;
             }
 
           else
             {
               // The new tolerance is based on the ratio of the most recent tolerances
               const Real alp=0.5*(1.+std::sqrt(5.));
               const Real gam=0.9;
 
               libmesh_assert_not_equal_to (nonlinear_residual_beforestep, 0.0);
               libmesh_assert_not_equal_to (nonlinear_residual_afterstep, 0.0);
 
               current_linear_tolerance = std::min(gam*std::pow(nonlinear_residual_afterstep/nonlinear_residual_beforestep, alp),
                                                   current_linear_tolerance*current_linear_tolerance
                                                   );
 
               // Don't let it get ridiculously small!!
               if (current_linear_tolerance < 1.e-12)
                 current_linear_tolerance = 1.e-12;
             }
 
           if (!quiet)
             libMesh::out << "Using current_linear_tolerance=" << current_linear_tolerance << std::endl;
 
 
           // Assemble the residual (and Jacobian).
           rhs_mode = Residual;
           assembly(true,   // Residual
                    true); // Jacobian
           rhs->close();
 
           // Save the current nonlinear residual.  We don't need to recompute the residual unless
           // this is the first step, since it was already computed as part of the convergence check
           // at the end of the last loop iteration.
           if (newton_step==0)
             {
               nonlinear_residual_beforestep = rhs->l2_norm();
 
               // Store the residual before any steps have been taken.  This will *not*
               // be updated at each step, and can be used to see if any progress has
               // been made from the initial residual at later steps.
               nonlinear_residual_firststep = nonlinear_residual_beforestep;
 
               const Real old_norm_u = solution->l2_norm();
               libMesh::out << "  (before step) ||R||_{L2} = " << nonlinear_residual_beforestep << std::endl;
               libMesh::out << "  (before step) ||R||_{L2}/||u|| = " << nonlinear_residual_beforestep / old_norm_u << std::endl;
 
               // In rare cases (very small arcsteps), it's possible that the residual is
               // already below our absolute linear tolerance.
               if (nonlinear_residual_beforestep  < solution_tolerance)
                 {
                   if (!quiet)
                     libMesh::out << "Initial guess satisfied linear tolerance, exiting with zero Newton iterations!" << std::endl;
 
                   // Since we go straight from here to the solve of the next tangent, we
                   // have to close the matrix before it can be assembled again.
                   matrix->close();
                   newton_converged=true;
                   break; // out of Newton iterations, with newton_converged=true
                 }
             }
 
           else
             {
               nonlinear_residual_beforestep = nonlinear_residual_afterstep;
             }
 
 
           // Solve the linear system G_u*z = G
           // Initial guess?
           z->zero(); // It seems to be extremely important to zero z here, otherwise the solver quits early.
           z->close();
 
           // It's possible that we have selected the current_linear_tolerance so large that
           // a guess of z=zero yields a linear system residual |Az + R| small enough that the
           // linear solver exits in zero iterations.  If this happens, we will reduce the
           // current_linear_tolerance until the linear solver does at least 1 iteration.
           do
             {
               rval =
                 linear_solver->solve(*matrix,
                                      *z,
                                      *rhs,
                                      //1.e-12,
                                      current_linear_tolerance,
                                      newton_solver->max_linear_iterations);   // max linear iterations
 
               if (rval.first==0)
                 {
                   if (newton_step==0)
                     {
                       libMesh::out << "Repeating initial solve with smaller linear tolerance!" << std::endl;
                       current_linear_tolerance *= initial_newton_tolerance; // reduce the linear tolerance to force the solver to do some work
                     }
                   else
                     // We shouldn't get here ... it means the linear solver did no work on a Newton
                     // step other than the first one.  If this happens, we need to think more about our
                     // tolerance selection.
                     libmesh_error_msg("Linear solver did no work!");
                 }
 
             } while (rval.first==0);
 
 
           if (!quiet)
             libMesh::out << "  G_u*z = G solver converged at step "
                          << rval.first
                          << " linear tolerance = "
                          << rval.second
                          << "."
                          << std::endl;
 
           // Sometimes (I am not sure why) the linear solver exits after zero iterations.
           // Perhaps it is hitting PETSc's divergence tolerance dtol???  If this occurs,
           // we should break out of the Newton iteration loop because nothing further is
           // going to happen...  Of course if the tolerance is already small enough after
           // zero iterations (how can this happen?!) we should not quit.
           if ((rval.first == 0) && (rval.second > current_linear_tolerance*nonlinear_residual_beforestep))
             {
               if (!quiet)
                 libMesh::out << "Linear solver exited in zero iterations!" << std::endl;
 
               // Try to find out the reason for convergence/divergence
               linear_solver->print_converged_reason();
 
               break; // out of Newton iterations
             }
 
           // Note: need to scale z by -1 since our code always solves Jx=R
           // instead of Jx=-R.
           z->scale(-1.);
           z->close();
 
 
 
 
 
 
           // Assemble the G_Lambda vector, skip residual.
           rhs_mode = G_Lambda;
 
           // Assemble both rhs and Jacobian
           assembly(true,  // Residual
                    false); // Jacobian
 
           // Not sure if this is really necessary
           rhs->close();
           const Real yrhsnorm=rhs->l2_norm();
           if (yrhsnorm == 0.0)
             libmesh_error_msg("||G_Lambda|| = 0");
 
           // We select a tolerance for the y-system which is based on the inexact Newton
           // tolerance but scaled by an extra term proportional to the RHS (which is not -> 0 in this case)
           const Real ysystemtol=current_linear_tolerance*(nonlinear_residual_beforestep/yrhsnorm);
           if (!quiet)
             libMesh::out << "ysystemtol=" << ysystemtol << std::endl;
 
           // Solve G_u*y = G_{\lambda}
           // FIXME: Initial guess?  This is really a solve for -du/dlambda so we could try
           // initializing it with the latest approximation to that... du/dlambda ~ du/ds * ds/dlambda
           //*y = *solution;
           //y->add(-1., *previous_u);
           //y->scale(-1. / (*continuation_parameter - old_continuation_parameter)); // Be careful of divide by zero...
           //y->close();
 
           //  const unsigned int max_attempts=1;
           // unsigned int attempt=0;
           //   do
           //     {
           //       if (!quiet)
           // libMesh::out << "Trying to solve tangent system, attempt " << attempt << std::endl;
 
           rval =
             linear_solver->solve(*matrix,
                                  *y,
                                  *rhs,
                                  //1.e-12,
                                  ysystemtol,
                                  newton_solver->max_linear_iterations);   // max linear iterations
 
           if (!quiet)
             libMesh::out << "  G_u*y = G_{lambda} solver converged at step "
                          << rval.first
                          << ", linear tolerance = "
                          << rval.second
                          << "."
                          << std::endl;
 
           // Sometimes (I am not sure why) the linear solver exits after zero iterations.
           // Perhaps it is hitting PETSc's divergence tolerance dtol???  If this occurs,
           // we should break out of the Newton iteration loop because nothing further is
           // going to happen...
           if ((rval.first == 0) && (rval.second > ysystemtol))
             {
               if (!quiet)
                 libMesh::out << "Linear solver exited in zero iterations!" << std::endl;
 
               break; // out of Newton iterations
             }
 
           //       ++attempt;
           //     } while ((attempt<max_attempts) && (rval.first==newton_solver->max_linear_iterations));
 
 
 
 
 
           // Compute N, the residual of the arclength constraint eqn.
           // Note 1: N(u,lambda,s) := (u-u_{old}, du_ds) + (lambda-lambda_{old}, dlambda_ds) - _ds
           // We temporarily use the delta_u vector as a temporary vector for this calculation.
           *delta_u = *solution;
           delta_u->add(-1., *previous_u);
 
           // First part of the arclength constraint
           const Number N1 = Theta_LOCA*Theta_LOCA*Theta*delta_u->dot(*du_ds);
           const Number N2 = ((*continuation_parameter) - old_continuation_parameter)*dlambda_ds;
           const Number N3 = ds_current;
 
           if (!quiet)
             {
               libMesh::out << "  N1=" << N1 << std::endl;
               libMesh::out << "  N2=" << N2 << std::endl;
               libMesh::out << "  N3=" << N3 << std::endl;
             }
 
           // The arclength constraint value
           const Number N = N1+N2-N3;
 
           if (!quiet)
             libMesh::out << "  N=" << N << std::endl;
 
           const Number duds_dot_z = du_ds->dot(*z);
           const Number duds_dot_y = du_ds->dot(*y);
 
           //libMesh::out << "duds_dot_z=" << duds_dot_z << std::endl;
           //libMesh::out << "duds_dot_y=" << duds_dot_y << std::endl;
           //libMesh::out << "dlambda_ds=" << dlambda_ds << std::endl;
 
           const Number delta_lambda_numerator   = -(N          + Theta_LOCA*Theta_LOCA*Theta*duds_dot_z);
           const Number delta_lambda_denominator =  (dlambda_ds - Theta_LOCA*Theta_LOCA*Theta*duds_dot_y);
 
           libmesh_assert_not_equal_to (delta_lambda_denominator, 0.0);
 
           // Now, we are ready to compute the step delta_lambda
           const Number delta_lambda_comp = delta_lambda_numerator /
             delta_lambda_denominator;
           // Lambda is real-valued
           const Real delta_lambda = libmesh_real(delta_lambda_comp);
 
           // Knowing delta_lambda, we are ready to update delta_u
           // delta_u = z - delta_lambda*y
           delta_u->zero();
           delta_u->add(1., *z);
           delta_u->add(-delta_lambda, *y);
           delta_u->close();
 
           // Update the system solution and the continuation parameter.
           solution->add(1., *delta_u);
           solution->close();
           *continuation_parameter += delta_lambda;
 
           // Did the Newton step actually reduce the residual?
           rhs_mode = Residual;
           assembly(true,   // Residual
                    false); // Jacobian
           rhs->close();
           nonlinear_residual_afterstep = rhs->l2_norm();
 
 
           // In a "normal" Newton step, ||du||/||R|| > 1 since the most recent
           // step is where you "just were" and the current residual is where
           // you are now.  It can occur that ||du||/||R|| < 1, but these are
           // likely not good cases to attempt backtracking (?).
           const Real norm_du_norm_R = delta_u->l2_norm() / nonlinear_residual_afterstep;
           if (!quiet)
             libMesh::out << "  norm_du_norm_R=" << norm_du_norm_R << std::endl;
 
 
           // Factor to decrease the stepsize by for backtracking
           Real newton_stepfactor = 1.;
 
           const bool attempt_backtracking =
             (nonlinear_residual_afterstep > solution_tolerance)
             && (nonlinear_residual_afterstep > nonlinear_residual_beforestep)
             && (n_backtrack_steps>0)
             && (norm_du_norm_R > 1.)
             ;
 
           // If residual is not reduced, do Newton back tracking.
           if (attempt_backtracking)
             {
               if (!quiet)
                 libMesh::out << "Newton step did not reduce residual." << std::endl;
 
               // back off the previous step.
               solution->add(-1., *delta_u);
               solution->close();
               *continuation_parameter -= delta_lambda;
 
               // Backtracking: start cutting the Newton stepsize by halves until
               // the new residual is actually smaller...
               for (unsigned int backtrack_step=0; backtrack_step<n_backtrack_steps; ++backtrack_step)
                 {
                   newton_stepfactor *= 0.5;
 
                   if (!quiet)
                     libMesh::out << "Shrinking step size by " << newton_stepfactor << std::endl;
 
                   // Take fractional step
                   solution->add(newton_stepfactor, *delta_u);
                   solution->close();
                   *continuation_parameter += newton_stepfactor*delta_lambda;
 
                   rhs_mode = Residual;
                   assembly(true,   // Residual
                            false); // Jacobian
                   rhs->close();
                   nonlinear_residual_afterstep = rhs->l2_norm();
 
                   if (!quiet)
                     libMesh::out << "At shrink step "
                                  << backtrack_step
                                  << ", nonlinear_residual_afterstep="
                                  << nonlinear_residual_afterstep
                                  << std::endl;
 
                   if (nonlinear_residual_afterstep < nonlinear_residual_beforestep)
                     {
                       if (!quiet)
                         libMesh::out << "Backtracking succeeded!" << std::endl;
 
                       break; // out of backtracking loop
                     }
 
                   else
                     {
                       // Back off that step
                       solution->add(-newton_stepfactor, *delta_u);
                       solution->close();
                       *continuation_parameter -= newton_stepfactor*delta_lambda;
                     }
 
                   // Save a copy of the solution from before the Newton step.
                   //std::unique_ptr<NumericVector<Number>> prior_iterate = solution->clone();
                 }
             } // end if (attempte_backtracking)
 
 
           // If we tried backtracking but the residual is still not reduced, print message.
           if ((attempt_backtracking) && (nonlinear_residual_afterstep > nonlinear_residual_beforestep))
             {
               //libMesh::err << "Backtracking failed." << std::endl;
               libMesh::out << "Backtracking failed." << std::endl;
 
               // 1.) Quit, exit program.
               //libmesh_error_msg("Backtracking failed!");
 
               // 2.) Continue with last newton_stepfactor
               if (newton_step<3)
                 {
                   solution->add(newton_stepfactor, *delta_u);
                   solution->close();
                   *continuation_parameter += newton_stepfactor*delta_lambda;
                   if (!quiet)
                     libMesh::out << "Backtracking could not reduce residual ... continuing anyway!" << std::endl;
                 }
 
               // 3.) Break out of Newton iteration loop with newton_converged = false,
               //     reduce the arclength stepsize, and try again.
               else
                 {
                   break; // out of Newton iteration loop, with newton_converged=false
                 }
             }
 
           // Another type of convergence check: suppose the residual has not been reduced
           // from its initial value after half of the allowed Newton steps have occurred.
           // In our experience, this typically means that it isn't going to converge and
           // we could probably save time by dropping out of the Newton iteration loop and
           // trying a smaller arcstep.
           if (this->newton_progress_check)
             {
               if ((nonlinear_residual_afterstep > nonlinear_residual_firststep) &&
                   (newton_step+1 > static_cast<unsigned int>(0.5*newton_solver->max_nonlinear_iterations)))
                 {
                   libMesh::out << "Progress check failed: the current residual: "
                                << nonlinear_residual_afterstep
                                << ", is\n"
                                << "larger than the initial residual, and half of the allowed\n"
                                << "number of Newton iterations have elapsed.\n"
                                << "Exiting Newton iterations with converged==false." << std::endl;
 
                   break; // out of Newton iteration loop, newton_converged = false
                 }
             }
 
           // Safety check: Check the current continuation parameter against user-provided min-allowable parameter value
           if (*continuation_parameter < min_continuation_parameter)
             {
               libMesh::out << "Continuation parameter fell below min-allowable value." << std::endl;
               break; // out of Newton iteration loop, newton_converged = false
             }
 
           // Safety check: Check the current continuation parameter against user-provided max-allowable parameter value
           if ( (max_continuation_parameter != 0.0) &&
                (*continuation_parameter > max_continuation_parameter) )
             {
               libMesh::out << "Current continuation parameter value: "
                            << *continuation_parameter
                            << " exceeded max-allowable value."
                            << std::endl;
               break; // out of Newton iteration loop, newton_converged = false
             }
 
 
           // Check the convergence of the parameter and the solution.  If they are small
           // enough, we can break out of the Newton iteration loop.
           const Real norm_delta_u = delta_u->l2_norm();
           const Real norm_u = solution->l2_norm();
           libMesh::out << "  delta_lambda                   = " << delta_lambda << std::endl;
           libMesh::out << "  newton_stepfactor*delta_lambda = " << newton_stepfactor*delta_lambda << std::endl;
           libMesh::out << "  lambda_current                 = " << *continuation_parameter << std::endl;
           libMesh::out << "  ||delta_u||                    = " << norm_delta_u << std::endl;
           libMesh::out << "  ||delta_u||/||u||              = " << norm_delta_u / norm_u << std::endl;
 
 
           // Evaluate the residual at the current Newton iterate.  We don't want to detect
           // convergence due to a small Newton step when the residual is still not small.
           rhs_mode = Residual;
           assembly(true,   // Residual
                    false); // Jacobian
           rhs->close();
           const Real norm_residual = rhs->l2_norm();
           libMesh::out << "  ||R||_{L2} = " << norm_residual << std::endl;
           libMesh::out << "  ||R||_{L2}/||u|| = " << norm_residual / norm_u << std::endl;
 
 
           // FIXME: The norm_delta_u tolerance (at least) should be relative.
           // It doesn't make sense to converge a solution whose size is ~ 10^5 to
           // a tolerance of 1.e-6.  Oh, and we should also probably check the
           // (relative) size of the residual as well, instead of just the step.
           if ((std::abs(delta_lambda) < continuation_parameter_tolerance) &&
               //(norm_delta_u       < solution_tolerance)               && // This is a *very* strict criterion we can probably skip
               (norm_residual      < solution_tolerance))
             {
               if (!quiet)
                 libMesh::out << "Newton iterations converged!" << std::endl;
 
               newton_converged = true;
               break; // out of Newton iterations
             }
         } // end nonlinear loop
 
       if (!newton_converged)
         {
           libMesh::out << "Newton iterations of augmented system did not converge!" << std::endl;
 
           // Reduce ds_current, recompute the solution and parameter, and continue to next
           // arcstep, if there is one.
           ds_current *= 0.5;
 
           // Go back to previous solution and parameter value.
           *solution = *previous_u;
           *continuation_parameter = old_continuation_parameter;
 
           // Compute new predictor with smaller ds
           apply_predictor();
         }
       else
         {
           // Set step convergence and break out
           arcstep_converged=true;
           break; // out of arclength reduction loop
         }
 
     } // end loop over arclength reductions
 
   // Check for convergence of the whole arcstep.  If not converged at this
   // point, we have no choice but to quit.
   if (!arcstep_converged)
     libmesh_error_msg("Arcstep failed to converge after max number of reductions! Exiting...");
 
   // Print converged solution control parameter and max value.
   libMesh::out << "lambda_current=" << *continuation_parameter << std::endl;
   //libMesh::out << "u_max=" << solution->max() << std::endl;
 
   // Reset old stream precision and flags.
   libMesh::out.precision(old_precision);
   libMesh::out.unsetf(std::ios_base::scientific);
 
   // Note: we don't want to go on to the next guess yet, since the user may
   // want to post-process this data.  It's up to the user to call advance_arcstep()
   // when they are ready to go on.
 }

◆ current_solution()

Number libMesh::System::current_solution ( const dof_id_type global_dof_number ) const

inherited

Returns: The current solution for the specified global DOF.

Definition at line 194 of file system.C.

References libMesh::System::_dof_map, and libMesh::System::current_local_solution.

Referenced by libMesh::ExactSolution::_compute_error(), libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::HPCoarsenTest::add_projection(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::HPCoarsenTest::select_refinement(), libMesh::EnsightIO::write_scalar_ascii(), and libMesh::EnsightIO::write_vector_ascii().

 {
   // Check the sizes
   libmesh_assert_less (global_dof_number, _dof_map->n_dofs());
   libmesh_assert_less (global_dof_number, current_local_solution->size());
 
   return (*current_local_solution)(global_dof_number);
 }

◆ damping_residual()

virtual bool libMesh::DifferentiablePhysics::damping_residual	(	bool	request_jacobian,
		DiffContext &
	)

inlinevirtualinherited

Subtracts a damping vector contribution on elem from elem_residual. This method is not used in first-order-in-time problems. For second-order-in-time problems, this is the $C(u,\ddot{u})\ddot{u}$ term. This method is only called for UnsteadySolver-based TimeSolvers.

If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.

If the problem has no damping, the default "do-nothing" is correct. Otherwise, this must be reimplemented.

Definition at line 374 of file diff_physics.h.

Referenced by libMesh::EulerSolver::element_residual(), libMesh::Euler2Solver::element_residual(), and libMesh::NewmarkSolver::element_residual().

                                                 {
     return request_jacobian;
   }

◆ deactivate()

void libMesh::System::deactivate ( )

inlineinherited

Deactivates the system. Only active systems are solved.

Definition at line 2081 of file system.h.

References libMesh::System::_active.

 {
   _active = false;
 }

◆ disable_cache()

void libMesh::ImplicitSystem::disable_cache ( )

overridevirtualinherited

Avoids use of any cached data that might affect any solve result. Should be overridden in derived systems.

Reimplemented from libMesh::System.

Definition at line 306 of file implicit_system.C.

References libMesh::System::assemble_before_solve, libMesh::ImplicitSystem::get_linear_solver(), and libMesh::LinearSolver< T >::reuse_preconditioner().

                                     {
   this->assemble_before_solve = true;
   this->get_linear_solver()->reuse_preconditioner(false);
 }

◆ disable_print_counter_info()

void libMesh::ReferenceCounter::disable_print_counter_info ( )

staticinherited

Definition at line 106 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

Referenced by libMesh::LibMeshInit::LibMeshInit().

 {
   _enable_print_counter = false;
   return;
 }

◆ element_constraint()

virtual bool libMesh::DifferentiablePhysics::element_constraint	(	bool	request_jacobian,
		DiffContext &
	)

inlinevirtualinherited

Adds the constraint contribution on elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.

Users may need to reimplement this for their particular PDE.

To implement the constraint 0 = G(u), the user should examine u = elem_solution and add (G(u), phi_i) to elem_residual in elem_constraint().

Definition at line 142 of file diff_physics.h.

Referenced by libMesh::EulerSolver::element_residual(), libMesh::Euler2Solver::element_residual(), libMesh::SteadySolver::element_residual(), libMesh::EigenTimeSolver::element_residual(), and libMesh::NewmarkSolver::element_residual().

                                                   {
     return request_jacobian;
   }

◆ element_postprocess()

virtual void libMesh::DifferentiableSystem::element_postprocess ( DiffContext & )

inlinevirtualinherited

Does any work that needs to be done on elem in a postprocessing loop.

Definition at line 282 of file diff_system.h.

282 {}

◆ element_qoi()

virtual void libMesh::DifferentiableQoI::element_qoi	(	DiffContext &	,
		const QoISet &
	)

inlinevirtualinherited

Does any work that needs to be done on elem in a quantity of interest assembly loop, outputting to elem_qoi.

Only qois included in the supplied QoISet need to be assembled.

Definition at line 108 of file diff_qoi.h.

110 {}

◆ element_qoi_derivative()

virtual void libMesh::DifferentiableQoI::element_qoi_derivative	(	DiffContext &	,
		const QoISet &
	)

inlinevirtualinherited

Does any work that needs to be done on elem in a quantity of interest derivative assembly loop, outputting to elem_qoi_derivative

Only qois included in the supplied QoISet need their derivatives assembled.

Definition at line 120 of file diff_qoi.h.

122 {}

◆ element_time_derivative()

virtual bool libMesh::DifferentiablePhysics::element_time_derivative	(	bool	request_jacobian,
		DiffContext &
	)

inlinevirtualinherited

Adds the time derivative contribution on elem to elem_residual. If this method receives request_jacobian = true, then it should compute elem_jacobian and return true if possible. If elem_jacobian has not been computed then the method should return false.

Users need to reimplement this for their particular PDE.

To implement the physics model du/dt = F(u), the user should examine u = elem_solution and add (F(u), phi_i) to elem_residual in elem_time_derivative().

Definition at line 124 of file diff_physics.h.

Referenced by libMesh::SteadySolver::element_residual(), and libMesh::EigenTimeSolver::element_residual().

                                                        {
     return request_jacobian;
   }

◆ enable_print_counter_info()

void libMesh::ReferenceCounter::enable_print_counter_info ( )

staticinherited

Methods to enable/disable the reference counter output from print_info()

Definition at line 100 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

 {
   _enable_print_counter = true;
   return;
 }

◆ eulerian_residual() [1/2]

virtual bool libMesh::FEMPhysics::eulerian_residual	(	bool	request_jacobian,
		DiffContext &	context
	)

overridevirtualinherited

Adds a pseudo-convection contribution on elem to elem_residual, if the nodes of elem are being translated by a moving mesh.

This function assumes that the user's time derivative equations (except for any equations involving unknown mesh xyz coordinates themselves) are expressed in an Eulerian frame of reference, and that the user is satisfied with an unstabilized convection term. Lagrangian equations will probably require overriding eulerian_residual() with a blank function; ALE or stabilized formulations will require reimplementing eulerian_residual() entirely.

Reimplemented from libMesh::DifferentiablePhysics.

◆ eulerian_residual() [2/2]

virtual bool libMesh::DifferentiablePhysics::eulerian_residual	(	bool	request_jacobian,
		DiffContext &
	)

inlinevirtualinherited

Adds a pseudo-convection contribution on elem to elem_residual, if the nodes of elem are being translated by a moving mesh.

The library provides a basic implementation in FEMPhysics::eulerian_residual()

Reimplemented in libMesh::FEMPhysics.

Definition at line 294 of file diff_physics.h.

                                                  {
     return request_jacobian;
   }

◆ finalize_derivative()

virtual void libMesh::DifferentiableQoI::finalize_derivative	(	NumericVector< Number > &	derivatives,
		std::size_t	qoi_index
	)

virtualinherited

Method to finalize qoi derivatives which require more than just a simple sum of element contributions.

Referenced by libMesh::FEMSystem::assemble_qoi_derivative().

◆ forward_qoi_parameter_sensitivity()

void libMesh::ImplicitSystem::forward_qoi_parameter_sensitivity	(	const QoISet &	qoi_indices,
		const ParameterVector &	parameters,
		SensitivityData &	sensitivities
	)

overridevirtualinherited

Solves for the derivative of each of the system's quantities of interest q in qoi[qoi_indices] with respect to each parameter in parameters, placing the result for qoi i and parameter j into sensitivities[i][j].

Uses the forward sensitivity method.

Currently uses finite differenced derivatives (partial q / partial p) and (partial R / partial p).

Reimplemented from libMesh::System.

Definition at line 807 of file implicit_system.C.

References std::abs(), libMesh::SensitivityData::allocate_data(), libMesh::ExplicitSystem::assemble_qoi(), libMesh::ExplicitSystem::assemble_qoi_derivative(), libMesh::ImplicitSystem::assembly(), libMesh::NumericVector< T >::close(), libMesh::SparseMatrix< T >::close(), libMesh::NumericVector< T >::dot(), libMesh::System::get_adjoint_rhs(), libMesh::System::get_sensitivity_solution(), libMesh::QoISet::has_index(), libMesh::ImplicitSystem::matrix, std::max(), libMesh::System::n_qois(), libMesh::System::qoi, libMesh::Real, libMesh::ExplicitSystem::rhs, libMesh::ImplicitSystem::sensitivity_solve(), libMesh::ParameterVector::size(), and libMesh::TOLERANCE.

 {
   ParameterVector & parameters =
     const_cast<ParameterVector &>(parameters_in);
 
   const unsigned int Np = cast_int<unsigned int>
     (parameters.size());
   const unsigned int Nq = this->n_qois();
 
   // An introduction to the problem:
   //
   // Residual R(u(p),p) = 0
   // partial R / partial u = J = system matrix
   //
   // This implies that:
   // d/dp(R) = 0
   // (partial R / partial p) +
   // (partial R / partial u) * (partial u / partial p) = 0
 
   // We first solve for (partial u / partial p) for each parameter:
   // J * (partial u / partial p) = - (partial R / partial p)
 
   this->sensitivity_solve(parameters);
 
   // Get ready to fill in sensitivities:
   sensitivities.allocate_data(qoi_indices, *this, parameters);
 
   // We use the identity:
   // dq/dp = (partial q / partial p) + (partial q / partial u) *
   //         (partial u / partial p)
 
   // We get (partial q / partial u) from the user
   this->assemble_qoi_derivative(qoi_indices,
                                 /* include_liftfunc = */ true,
                                 /* apply_constraints = */ false);
 
   // We don't need these to be closed() in this function, but libMesh
   // standard practice is to have them closed() by the time the
   // function exits
   for (unsigned int i=0; i != this->n_qois(); ++i)
     if (qoi_indices.has_index(i))
       this->get_adjoint_rhs(i).close();
 
   for (unsigned int j=0; j != Np; ++j)
     {
       // We currently get partial derivatives via central differencing
 
       // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
 
       Number old_parameter = *parameters[j];
 
       const Real delta_p =
         TOLERANCE * std::max(std::abs(old_parameter), 1e-3);
 
       *parameters[j] = old_parameter - delta_p;
       this->assemble_qoi(qoi_indices);
       std::vector<Number> qoi_minus = this->qoi;
 
       *parameters[j] = old_parameter + delta_p;
       this->assemble_qoi(qoi_indices);
       std::vector<Number> & qoi_plus = this->qoi;
 
       std::vector<Number> partialq_partialp(Nq, 0);
       for (unsigned int i=0; i != Nq; ++i)
         if (qoi_indices.has_index(i))
           partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
 
       // Don't leave the parameter changed
       *parameters[j] = old_parameter;
 
       for (unsigned int i=0; i != Nq; ++i)
         if (qoi_indices.has_index(i))
           sensitivities[i][j] = partialq_partialp[i] +
             this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(j));
     }
 
   // All parameters have been reset.
   // We didn't cache the original rhs or matrix for memory reasons,
   // but we can restore them to a state consistent solution -
   // principle of least surprise.
   this->assembly(true, true);
   this->rhs->close();
   this->matrix->close();
   this->assemble_qoi(qoi_indices);
 }

◆ get_adjoint_rhs() [1/2]

NumericVector< Number > & libMesh::System::get_adjoint_rhs ( unsigned int i = 0 )

inherited

Returns: A reference to one of the system's adjoint rhs vectors, by default the one corresponding to the first qoi. This what the user's QoI derivative code should assemble when setting up an adjoint problem

Definition at line 1031 of file system.C.

References libMesh::System::get_vector().

Referenced by libMesh::ImplicitSystem::adjoint_solve(), libMesh::ImplicitSystem::forward_qoi_parameter_sensitivity(), libMesh::ImplicitSystem::qoi_parameter_hessian(), libMesh::ImplicitSystem::qoi_parameter_hessian_vector_product(), and libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve().

 {
   std::ostringstream adjoint_rhs_name;
   adjoint_rhs_name << "adjoint_rhs" << i;
 
   return this->get_vector(adjoint_rhs_name.str());
 }

◆ get_adjoint_rhs() [2/2]

const NumericVector< Number > & libMesh::System::get_adjoint_rhs ( unsigned int i = 0 ) const

inherited

Returns: A reference to one of the system's adjoint rhs vectors, by default the one corresponding to the first qoi.

Definition at line 1041 of file system.C.

References libMesh::System::get_vector().

 {
   std::ostringstream adjoint_rhs_name;
   adjoint_rhs_name << "adjoint_rhs" << i;
 
   return this->get_vector(adjoint_rhs_name.str());
 }

◆ get_adjoint_solution() [1/2]

NumericVector< Number > & libMesh::System::get_adjoint_solution ( unsigned int i = 0 )

inherited

Returns: A reference to one of the system's adjoint solution vectors, by default the one corresponding to the first qoi.

Definition at line 969 of file system.C.

References libMesh::System::get_vector().

Referenced by libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::ImplicitSystem::adjoint_solve(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::AdjointResidualErrorEstimator::estimate_error(), libMesh::ImplicitSystem::qoi_parameter_hessian(), libMesh::ImplicitSystem::qoi_parameter_hessian_vector_product(), and libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve().

 {
   std::ostringstream adjoint_name;
   adjoint_name << "adjoint_solution" << i;
 
   return this->get_vector(adjoint_name.str());
 }

◆ get_adjoint_solution() [2/2]

const NumericVector< Number > & libMesh::System::get_adjoint_solution ( unsigned int i = 0 ) const

inherited

Returns: A reference to one of the system's adjoint solution vectors, by default the one corresponding to the first qoi.

Definition at line 979 of file system.C.

References libMesh::System::get_vector().

 {
   std::ostringstream adjoint_name;
   adjoint_name << "adjoint_solution" << i;
 
   return this->get_vector(adjoint_name.str());
 }

◆ get_all_variable_numbers()

void libMesh::System::get_all_variable_numbers ( std::vector< unsigned int > & all_variable_numbers ) const

inherited

Fills all_variable_numbers with all the variable numbers for the variables that have been added to this system.

Definition at line 1258 of file system.C.

References libMesh::System::_variable_numbers, and libMesh::System::n_vars().

 {
   all_variable_numbers.resize(n_vars());
 
   // Make sure the variable exists
   std::map<std::string, unsigned short int>::const_iterator
     it = _variable_numbers.begin();
   std::map<std::string, unsigned short int>::const_iterator
     it_end = _variable_numbers.end();
 
   unsigned int count = 0;
   for ( ; it != it_end; ++it)
     {
       all_variable_numbers[count] = it->second;
       count++;
     }
 }

◆ get_dof_map() [1/2]

const DofMap & libMesh::System::get_dof_map ( ) const

inlineinherited

Returns: A constant reference to this system's _dof_map.

Definition at line 2049 of file system.h.

References libMesh::System::_dof_map.

 {
   return *_dof_map;
 }

◆ get_dof_map() [2/2]

DofMap & libMesh::System::get_dof_map ( )

inlineinherited

Returns: A writable reference to this system's _dof_map.

Definition at line 2057 of file system.h.

References libMesh::System::_dof_map.

 {
   return *_dof_map;
 }

◆ get_equation_systems() [1/2]

const EquationSystems& libMesh::System::get_equation_systems ( ) const

inlineinherited

Returns: A constant reference to this system's parent EquationSystems object.

Definition at line 712 of file system.h.

References libMesh::System::_equation_systems.

712 { return _equation_systems; }

libMesh::System::_equation_systems

EquationSystems & _equation_systems

Definition: system.h:1890

◆ get_equation_systems() [2/2]

EquationSystems& libMesh::System::get_equation_systems ( )

inlineinherited

Returns: A reference to this system's parent EquationSystems object.

Definition at line 717 of file system.h.

References libMesh::System::_equation_systems.

717 { return _equation_systems; }

libMesh::System::_equation_systems

EquationSystems & _equation_systems

Definition: system.h:1890

◆ get_first_order_vars()

const std::set<unsigned int>& libMesh::DifferentiablePhysics::get_first_order_vars ( ) const

inlineinherited

Returns: The set of first order in time variable indices. May be empty.

Definition at line 517 of file diff_physics.h.

References libMesh::DifferentiablePhysics::_first_order_vars.

Referenced by libMesh::DifferentiableSystem::have_first_order_scalar_vars().

518 { return _first_order_vars; }

libMesh::DifferentiablePhysics::_first_order_vars

std::set< unsigned int > _first_order_vars

Definition: diff_physics.h:559

◆ get_info() [1/2]

std::string libMesh::ReferenceCounter::get_info ( )

staticinherited

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

References libMesh::ReferenceCounter::_counts, and libMesh::Quality::name().

Referenced by libMesh::ReferenceCounter::print_info().

 {
 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
 
   std::ostringstream oss;
 
   oss << '\n'
       << " ---------------------------------------------------------------------------- \n"
       << "| Reference count information                                                |\n"
       << " ---------------------------------------------------------------------------- \n";
 
   for (const auto & pr : _counts)
     {
       const std::string name(pr.first);
       const unsigned int creations    = pr.second.first;
       const unsigned int destructions = pr.second.second;
 
       oss << "| " << name << " reference count information:\n"
           << "|  Creations:    " << creations    << '\n'
           << "|  Destructions: " << destructions << '\n';
     }
 
   oss << " ---------------------------------------------------------------------------- \n";
 
   return oss.str();
 
 #else
 
   return "";
 
 #endif
 }

◆ get_info() [2/2]

std::string libMesh::System::get_info ( ) const

inherited

Returns: A string containing information about the system.

Definition at line 1658 of file system.C.

 {
   std::ostringstream oss;
 
 
   const std::string & sys_name = this->name();
 
   oss << "   System #"  << this->number() << ", \"" << sys_name << "\"\n"
       << "    Type \""  << this->system_type() << "\"\n"
       << "    Variables=";
 
   for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
     {
       const VariableGroup & vg_description (this->variable_group(vg));
 
       if (vg_description.n_variables() > 1) oss << "{ ";
       for (unsigned int vn=0; vn<vg_description.n_variables(); vn++)
         oss << "\"" << vg_description.name(vn) << "\" ";
       if (vg_description.n_variables() > 1) oss << "} ";
     }
 
   oss << '\n';
 
   oss << "    Finite Element Types=";
 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
   for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
     oss << "\""
         << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
         << "\" ";
 #else
   for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
     {
       oss << "\""
           << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
           << "\", \""
           << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().radial_family)
           << "\" ";
     }
 
   oss << '\n' << "    Infinite Element Mapping=";
   for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
     oss << "\""
         << Utility::enum_to_string<InfMapType>(this->get_dof_map().variable_group(vg).type().inf_map)
         << "\" ";
 #endif
 
   oss << '\n';
 
   oss << "    Approximation Orders=";
   for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
     {
 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
       oss << "\""
           << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
           << "\" ";
 #else
       oss << "\""
           << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
           << "\", \""
           << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().radial_order)
           << "\" ";
 #endif
     }
 
   oss << '\n';
 
   oss << "    n_dofs()="             << this->n_dofs()             << '\n';
   oss << "    n_local_dofs()="       << this->n_local_dofs()       << '\n';
 #ifdef LIBMESH_ENABLE_CONSTRAINTS
   oss << "    n_constrained_dofs()=" << this->n_constrained_dofs() << '\n';
   oss << "    n_local_constrained_dofs()=" << this->n_local_constrained_dofs() << '\n';
 #endif
 
   oss << "    " << "n_vectors()="  << this->n_vectors()  << '\n';
   oss << "    " << "n_matrices()="  << this->n_matrices()  << '\n';
   //   oss << "    " << "n_additional_matrices()=" << this->n_additional_matrices() << '\n';
 
   oss << this->get_dof_map().get_info();
 
   return oss.str();
 }

◆ get_linear_solve_parameters()

std::pair< unsigned int, Real > libMesh::DifferentiableSystem::get_linear_solve_parameters ( ) const

overridevirtualinherited

Returns: An integer corresponding to the upper iteration count limit and a Real corresponding to the convergence tolerance to be used in linear adjoint and/or sensitivity solves

Reimplemented from libMesh::ImplicitSystem.

Definition at line 184 of file diff_system.C.

References libMesh::DifferentiableSystem::time_solver.

 {
   libmesh_assert(time_solver.get());
   libmesh_assert_equal_to (&(time_solver->system()), this);
   return std::make_pair(this->time_solver->diff_solver()->max_linear_iterations,
                         this->time_solver->diff_solver()->relative_residual_tolerance);
 }

◆ get_linear_solver()

LinearSolver< Number > * libMesh::DifferentiableSystem::get_linear_solver ( ) const

overridevirtualinherited

Returns: A pointer to a linear solver appropriate for use in adjoint and/or sensitivity solves

Reimplemented from libMesh::ImplicitSystem.

Definition at line 175 of file diff_system.C.

References libMesh::DifferentiableSystem::time_solver.

 {
   libmesh_assert(time_solver.get());
   libmesh_assert_equal_to (&(time_solver->system()), this);
   return this->time_solver->linear_solver().get();
 }

◆ get_matrix() [1/2]

const SparseMatrix< Number > & libMesh::ImplicitSystem::get_matrix ( const std::string & mat_name ) const

inherited

Returns: A const reference to this system's additional matrix named mat_name.

None of these matrices is involved in the solution process. Access is only granted when the matrix is already properly initialized.

Definition at line 263 of file implicit_system.C.

References libMesh::ImplicitSystem::_matrices.

Referenced by libMesh::NewmarkSystem::compute_matrix(), libMesh::EigenTimeSolver::solve(), and libMesh::NewmarkSystem::update_rhs().

 {
   // Make sure the matrix exists
   const_matrices_iterator pos = _matrices.find (mat_name);
 
   if (pos == _matrices.end())
     libmesh_error_msg("ERROR: matrix " << mat_name << " does not exist in this system!");
 
   return *(pos->second);
 }

◆ get_matrix() [2/2]

SparseMatrix< Number > & libMesh::ImplicitSystem::get_matrix ( const std::string & mat_name )

inherited

Returns: A writable reference to this system's additional matrix named mat_name.

None of these matrices is involved in the solution process. Access is only granted when the matrix is already properly initialized.

Definition at line 276 of file implicit_system.C.

References libMesh::ImplicitSystem::_matrices.

 {
   // Make sure the matrix exists
   matrices_iterator pos = _matrices.find (mat_name);
 
   if (pos == _matrices.end())
     libmesh_error_msg("ERROR: matrix " << mat_name << " does not exist in this system!");
 
   return *(pos->second);
 }

◆ get_mesh() [1/2]

const MeshBase & libMesh::System::get_mesh ( ) const

inlineinherited

Returns: A constant reference to this systems's _mesh.

Definition at line 2033 of file system.h.

References libMesh::System::_mesh.

 {
   return _mesh;
 }

◆ get_mesh() [2/2]

MeshBase & libMesh::System::get_mesh ( )

inlineinherited

Returns: A reference to this systems's _mesh.

Definition at line 2041 of file system.h.

References libMesh::System::_mesh.

 {
   return _mesh;
 }

◆ get_mesh_system() [1/2]

const System * libMesh::DifferentiablePhysics::get_mesh_system ( ) const

inlineinherited

Returns: A const reference to the system with variables corresponding to mesh nodal coordinates, or nullptr if the mesh is fixed. Useful for ALE calculations.

Definition at line 624 of file diff_physics.h.

References libMesh::DifferentiablePhysics::_mesh_sys.

Referenced by libMesh::FEMSystem::build_context().

 {
   return _mesh_sys;
 }

◆ get_mesh_system() [2/2]

System * libMesh::DifferentiablePhysics::get_mesh_system ( )

inlineinherited

Returns: A reference to the system with variables corresponding to mesh nodal coordinates, or nullptr if the mesh is fixed.

Definition at line 630 of file diff_physics.h.

References libMesh::DifferentiablePhysics::_mesh_sys.

 {
   return _mesh_sys;
 }

◆ get_mesh_x_var()

unsigned int libMesh::DifferentiablePhysics::get_mesh_x_var ( ) const

inlineinherited

Returns: The variable number corresponding to the mesh x coordinate. Useful for ALE calculations.

Definition at line 636 of file diff_physics.h.

References libMesh::DifferentiablePhysics::_mesh_x_var.

Referenced by libMesh::FEMSystem::build_context().

 {
   return _mesh_x_var;
 }

◆ get_mesh_y_var()

unsigned int libMesh::DifferentiablePhysics::get_mesh_y_var ( ) const

inlineinherited

Returns: The variable number corresponding to the mesh y coordinate. Useful for ALE calculations.

Definition at line 642 of file diff_physics.h.

References libMesh::DifferentiablePhysics::_mesh_y_var.

Referenced by libMesh::FEMSystem::build_context().

 {
   return _mesh_y_var;
 }

◆ get_mesh_z_var()

unsigned int libMesh::DifferentiablePhysics::get_mesh_z_var ( ) const

inlineinherited

Returns: The variable number corresponding to the mesh z coordinate. Useful for ALE calculations.

Definition at line 648 of file diff_physics.h.

References libMesh::DifferentiablePhysics::_mesh_z_var.

Referenced by libMesh::FEMSystem::build_context().

 {
   return _mesh_z_var;
 }

◆ get_physics() [1/2]

const DifferentiablePhysics* libMesh::DifferentiableSystem::get_physics ( ) const

inlineinherited

Returns: A const reference to a DifferentiablePhysics object.

Note: If no external Physics object is attached, the default is this.

Definition at line 182 of file diff_system.h.

References libMesh::DifferentiableSystem::_diff_physics.

Referenced by libMesh::EulerSolver::_general_residual(), libMesh::Euler2Solver::_general_residual(), libMesh::SteadySolver::_general_residual(), libMesh::NewmarkSolver::_general_residual(), libMesh::FEMSystem::build_context(), libMesh::EulerSolver::element_residual(), libMesh::Euler2Solver::element_residual(), libMesh::EigenTimeSolver::element_residual(), libMesh::FEMSystem::init_context(), libMesh::EigenTimeSolver::nonlocal_residual(), and libMesh::EigenTimeSolver::side_residual().

183 { return this->_diff_physics; }

libMesh::DifferentiableSystem::_diff_physics

DifferentiablePhysics * _diff_physics

Definition: diff_system.h:371

◆ get_physics() [2/2]

DifferentiablePhysics* libMesh::DifferentiableSystem::get_physics ( )

inlineinherited

Returns: A reference to a DifferentiablePhysics object.

Note: If no external Physics object is attached, the default is this.

Definition at line 191 of file diff_system.h.

References libMesh::DifferentiableSystem::_diff_physics.

192 { return this->_diff_physics; }

libMesh::DifferentiableSystem::_diff_physics

DifferentiablePhysics * _diff_physics

Definition: diff_system.h:371

Public Types

Public Member Functions

Static Public Member Functions

Public Attributes

Protected Types

Protected Member Functions

Protected Attributes

Static Protected Attributes

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

◆ const_matrices_iterator

◆ const_vectors_iterator

◆ Counts

◆ GradientFunctionPointer

◆ matrices_iterator

◆ Parent

◆ sys_type

◆ TimeSolverResPtr

◆ ValueFunctionPointer

◆ vectors_iterator

Member Enumeration Documentation

◆ Predictor

◆ RHS_Mode

Constructor & Destructor Documentation

◆ ContinuationSystem()

◆ ~ContinuationSystem()

Member Function Documentation

◆ _eulerian_time_deriv()

◆ activate()

◆ active()

◆ add_adjoint_rhs()

◆ add_adjoint_solution()

◆ add_dot_var_dirichlet_bcs()

◆ add_matrix()

◆ add_second_order_dot_vars()

◆ add_sensitivity_rhs()

◆ add_sensitivity_solution()

◆ add_variable() [1/2]

◆ add_variable() [2/2]

◆ add_variables() [1/2]

◆ add_variables() [2/2]

◆ add_vector()

◆ add_weighted_sensitivity_adjoint_solution()

◆ add_weighted_sensitivity_solution()

◆ adjoint_qoi_parameter_sensitivity()

◆ adjoint_solve()

◆ advance_arcstep()

◆ apply_predictor()

◆ assemble()

◆ assemble_qoi()

◆ assemble_qoi_derivative()

◆ assemble_residual_derivatives()

◆ assembly()

◆ attach_assemble_function()

◆ attach_assemble_object()

◆ attach_constraint_function()

◆ attach_constraint_object()

◆ attach_init_function()

◆ attach_init_object()

◆ attach_physics()

◆ attach_qoi()

◆ attach_QOI_derivative()

◆ attach_QOI_derivative_object()

◆ attach_QOI_function()

◆ attach_QOI_object()

◆ boundary_project_solution() [1/2]

◆ boundary_project_solution() [2/2]

◆ boundary_project_vector() [1/2]

◆ boundary_project_vector() [2/2]

◆ build_context()

◆ calculate_norm() [1/2]

◆ calculate_norm() [2/2]

◆ clear()

◆ clear_physics()

◆ clear_qoi()

◆ clone()

◆ clone_physics()

◆ comm()