libMesh::ContinuationSystem Class Reference

#include <continuation_system.h>

Inheritance diagram for libMesh::ContinuationSystem:

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 () libmesh_override
 
virtual void solve () libmesh_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) libmesh_override
 
void mesh_position_get ()
 
void mesh_position_set ()
 
virtual std::unique_ptr< DiffContextbuild_context () libmesh_override
 
virtual void init_context (DiffContext &) libmesh_override
 
virtual void postprocess () libmesh_override
 
virtual void assemble_qoi (const QoISet &indices=QoISet()) libmesh_override
 
virtual void assemble_qoi_derivative (const QoISet &qoi_indices=QoISet(), bool include_liftfunc=true, bool apply_constraints=true) libmesh_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 () libmesh_override
 
virtual void assemble () libmesh_override
 
virtual LinearSolver< Number > * get_linear_solver () const libmesh_override
 
virtual std::pair< unsigned int, Realget_linear_solve_parameters () const libmesh_override
 
virtual void release_linear_solver (LinearSolver< Number > *) const libmesh_override
 
virtual std::pair< unsigned int, Realadjoint_solve (const QoISet &qoi_indices=QoISet()) libmesh_override
 
virtual std::unique_ptr< DifferentiablePhysicsclone_physics () libmesh_override
 
virtual std::unique_ptr< DifferentiableQoIclone () libmesh_override
 
const DifferentiablePhysicsget_physics () const
 
DifferentiablePhysicsget_physics ()
 
void attach_physics (DifferentiablePhysics *physics_in)
 
void swap_physics (DifferentiablePhysics *&swap_physics)
 
const DifferentiableQoIget_qoi () const
 
DifferentiableQoIget_qoi ()
 
void attach_qoi (DifferentiableQoI *qoi_in)
 
void set_time_solver (std::unique_ptr< TimeSolver > _time_solver)
 
TimeSolverget_time_solver ()
 
const TimeSolverget_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_typesystem ()
 
virtual void disable_cache () libmesh_override
 
virtual std::string system_type () const libmesh_override
 
virtual void assemble_residual_derivatives (const ParameterVector &parameters) libmesh_override
 
virtual std::pair< unsigned int, Realsensitivity_solve (const ParameterVector &parameters) libmesh_override
 
virtual std::pair< unsigned int, Realweighted_sensitivity_solve (const ParameterVector &parameters, const ParameterVector &weights) libmesh_override
 
virtual std::pair< unsigned int, Realweighted_sensitivity_adjoint_solve (const ParameterVector &parameters, const ParameterVector &weights, const QoISet &qoi_indices=QoISet()) libmesh_override
 
virtual void adjoint_qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities) libmesh_override
 
virtual void forward_qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities) libmesh_override
 
virtual void qoi_parameter_hessian (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &hessian) libmesh_override
 
virtual void qoi_parameter_hessian_vector_product (const QoISet &qoi_indices, const ParameterVector &parameters, const ParameterVector &vector, SensitivityData &product) libmesh_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 libmesh_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=libmesh_nullptr) const
 
void project_solution (FEMFunctionBase< Number > *f, FEMFunctionBase< Gradient > *g=libmesh_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=libmesh_nullptr, int is_adjoint=-1) const
 
void project_vector (NumericVector< Number > &new_vector, FEMFunctionBase< Number > *f, FEMFunctionBase< Gradient > *g=libmesh_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=libmesh_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=libmesh_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 MeshBaseget_mesh () const
 
MeshBaseget_mesh ()
 
const DofMapget_dof_map () const
 
DofMapget_dof_map ()
 
const EquationSystemsget_equation_systems () const
 
EquationSystemsget_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=libmesh_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=libmesh_nullptr)
 
unsigned int add_variables (const std::vector< std::string > &vars, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_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=libmesh_nullptr)
 
const Variablevariable (unsigned int var) const
 
const VariableGroupvariable_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 FETypevariable_type (const unsigned int i) const
 
const FETypevariable_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=libmesh_nullptr) const
 
Real calculate_norm (const NumericVector< Number > &v, const SystemNorm &norm, std::set< unsigned int > *skip_dimensions=libmesh_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
 
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::Communicatorcomm () 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) libmesh_override
 
virtual bool mass_residual (bool request_jacobian, DiffContext &)
 
virtual bool mass_residual (bool request_jacobian, DiffContext &) libmesh_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 Systemget_mesh_system () const
 
Systemget_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

Realcontinuation_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< TimeSolvertime_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< Numberqoi
 
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 () libmesh_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
 
DifferentiableQoIdiff_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
 
NewtonSolvernewton_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

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

Definition at line 290 of file implicit_system.h.

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

Definition at line 735 of file system.h.

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 119 of file reference_counter.h.

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 510 of file system.h.

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

Matrix iterator typedefs.

Definition at line 289 of file implicit_system.h.

The type of the parent.

Definition at line 79 of file continuation_system.h.

The type of system.

Definition at line 74 of file continuation_system.h.

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.

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 506 of file system.h.

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

Vector iterator typedefs.

Definition at line 734 of file system.h.

Member Enumeration Documentation

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.

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.

Constructor & Destructor Documentation

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().

31  :
32  Parent(es, name_in, number_in),
34  quiet(true),
36  solution_tolerance(1.e-6),
41  Theta(1.),
42  Theta_LOCA(1.),
45  ds_min(1.e-8),
51  tangent_initialized(false),
53  dlambda_ds(0.707),
54  ds(0.1),
55  ds_current(0.1),
57  previous_ds(0.),
58  newton_step(0)
59 {
60  // Warn about using untested code
61  libmesh_experimental();
62 
63  if (libMesh::on_command_line("--solver_system_names"))
64  linear_solver->init((this->name()+"_").c_str());
65  else
66  linear_solver->init();
67 }
static std::unique_ptr< LinearSolver< Number > > build(const libMesh::Parallel::Communicator &comm_in, const SolverPackage solver_package=libMesh::default_solver_package())
std::unique_ptr< LinearSolver< Number > > linear_solver
const class libmesh_nullptr_t libmesh_nullptr
const std::string & name() const
Definition: system.h:1998
bool on_command_line(const std::string &arg)
Definition: libmesh.C:920
libMesh::ContinuationSystem::~ContinuationSystem ( )
virtual

Destructor.

Definition at line 72 of file continuation_system.C.

References clear().

73 {
74  this->clear();
75 }
virtual void clear() libmesh_override

Member Function Documentation

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(), libMesh::NewmarkSolver::element_residual(), and libMesh::DifferentiablePhysics::init_context().

void libMesh::System::activate ( )
inlineinherited

Activates the system. Only active systems are solved.

Definition at line 2054 of file system.h.

References libMesh::System::_active.

Referenced by libMesh::System::get_equation_systems().

2055 {
2056  _active = true;
2057 }
bool libMesh::System::active ( ) const
inlineinherited
Returns
true if the system is active, false otherwise. An active system will be solved.

Definition at line 2046 of file system.h.

References libMesh::System::_active.

Referenced by libMesh::System::get_equation_systems().

2047 {
2048  return _active;
2049 }
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 1022 of file system.C.

References libMesh::System::add_vector().

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

1023 {
1024  std::ostringstream adjoint_rhs_name;
1025  adjoint_rhs_name << "adjoint_rhs" << i;
1026 
1027  return this->add_vector(adjoint_rhs_name.str(), false);
1028 }
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Definition: system.C:662
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 958 of file system.C.

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

Referenced by libMesh::ImplicitSystem::adjoint_solve(), and libMesh::System::project_solution_on_reinit().

959 {
960  std::ostringstream adjoint_name;
961  adjoint_name << "adjoint_solution" << i;
962 
963  NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
964  this->set_vector_as_adjoint(adjoint_name.str(), i);
965  return returnval;
966 }
void set_vector_as_adjoint(const std::string &vec_name, int qoi_num)
Definition: system.C:886
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Definition: system.C:662
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 228 of file diff_system.C.

References libMesh::DofMap::add_dirichlet_boundary(), libMesh::DirichletBoundary::b, libMesh::DirichletBoundary::f, libMesh::DirichletBoundary::f_fem, libMesh::DofMap::get_dirichlet_boundaries(), libMesh::System::get_dof_map(), and libMesh::DirichletBoundary::variables.

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

230 {
231  // We're assuming that there could be a lot more variables than
232  // boundary conditions, so we search each of the boundary conditions
233  // for this variable rather than looping over boundary conditions
234  // in a separate loop and searching through all the variables.
235  const DirichletBoundaries * all_dbcs =
237 
238  if (all_dbcs)
239  {
240  // We need to cache the DBCs to be added so that we add them
241  // after looping over the existing DBCs. Otherwise, we're polluting
242  // the thing we're looping over.
243  std::vector<DirichletBoundary*> new_dbcs;
244 
245  DirichletBoundaries::const_iterator dbc_it = all_dbcs->begin();
246  for ( ; dbc_it != all_dbcs->end(); ++dbc_it )
247  {
248  libmesh_assert(*dbc_it);
249  DirichletBoundary & dbc = *(*dbc_it);
250 
251  // Look for second order variable in the current
252  // DirichletBoundary object
253  std::vector<unsigned int>::const_iterator dbc_var_it =
254  std::find( dbc.variables.begin(), dbc.variables.end(), var_idx );
255 
256  // If we found it, then we also need to add it's corresponding
257  // "dot" variable to a DirichletBoundary
258  std::vector<unsigned int> vars_to_add;
259  if (dbc_var_it != dbc.variables.end())
260  vars_to_add.push_back(dot_var_idx);
261 
262  if (!vars_to_add.empty())
263  {
264  // We need to check if the boundary condition is time-dependent.
265  // Currently, we cannot automatically differentiate w.r.t. time
266  // so if the user supplies a time-dependent Dirichlet BC, then
267  // we can't automatically support the Dirichlet BC for the
268  // "velocity" boundary condition, so we error. Otherwise,
269  // the "velocity boundary condition will just be zero.
270  bool is_time_evolving_bc = false;
271  if (dbc.f)
272  is_time_evolving_bc = dbc.f->is_time_dependent();
273  else if (dbc.f_fem)
274  // We it's a FEMFunctionBase object, it will be implicitly
275  // time-dependent since it is assumed to depend on the solution.
276  is_time_evolving_bc = true;
277  else
278  libmesh_error_msg("Could not find valid boundary function!");
279 
280  if (is_time_evolving_bc)
281  libmesh_error_msg("Cannot currently support time-dependent Dirichlet BC for dot variables!");
282 
283 
284  DirichletBoundary * new_dbc;
285 
286  if (dbc.f)
287  {
288  ZeroFunction<Number> zf;
289 
290  new_dbc = new DirichletBoundary(dbc.b, vars_to_add, zf);
291  }
292  else
293  libmesh_error();
294 
295  new_dbcs.push_back(new_dbc);
296  }
297  }
298 
299  // Now add the new DBCs for the "dot" vars to the DofMap
300  std::vector<DirichletBoundary*>::iterator new_dbc_it =
301  new_dbcs.begin();
302 
303  for ( ; new_dbc_it != new_dbcs.end(); ++new_dbc_it )
304  {
305  const DirichletBoundary & dbc = *(*new_dbc_it);
306  this->get_dof_map().add_dirichlet_boundary(dbc);
307  delete *new_dbc_it;
308  }
309 
310  } // if (all_dbcs)
311 }
const DofMap & get_dof_map() const
Definition: system.h:2030
const DirichletBoundaries * get_dirichlet_boundaries() const
Definition: dof_map.h:1161
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
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 210 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().

211 {
212  // only add matrices before initializing...
213  if (!_can_add_matrices)
214  libmesh_error_msg("ERROR: Too late. Cannot add matrices to the system after initialization"
215  << "\n any more. You should have done this earlier.");
216 
217  // Return the matrix if it is already there.
218  if (this->have_matrix(mat_name))
219  return *(_matrices[mat_name]);
220 
221  // Otherwise build the matrix and return it.
222  SparseMatrix<Number> * buf = SparseMatrix<Number>::build(this->comm()).release();
223  _matrices.insert (std::make_pair (mat_name, buf));
224 
225  return *buf;
226 }
static std::unique_ptr< SparseMatrix< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
std::map< std::string, SparseMatrix< Number > * > _matrices
bool have_matrix(const std::string &mat_name) const
const Parallel::Communicator & comm() const
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().

199 {
200  const std::set<unsigned int> & second_order_vars = this->get_second_order_vars();
201  if (!second_order_vars.empty())
202  {
203  for (std::set<unsigned int>::const_iterator var_it = second_order_vars.begin();
204  var_it != second_order_vars.end(); ++var_it)
205  {
206  const Variable & var = this->variable(*var_it);
207  std::string new_var_name = std::string("dot_")+var.name();
208 
209  unsigned int v_var_idx;
210 
211  if (var.active_subdomains().empty())
212  v_var_idx = this->add_variable( new_var_name, var.type() );
213  else
214  v_var_idx = this->add_variable( new_var_name, var.type(), &var.active_subdomains() );
215 
216  _second_order_dot_vars.insert( std::pair<unsigned int,unsigned int>(*var_it,v_var_idx) );
217 
218  // The new velocities are time evolving variables of first order
219  this->time_evolving( v_var_idx, 1 );
220 
221  // And if there are any boundary conditions set on the second order
222  // variable, we also need to set it on its velocity variable.
223  this->add_dot_var_dirichlet_bcs( *var_it, v_var_idx );
224  }
225  }
226 }
const std::set< unsigned int > & get_second_order_vars() const
Definition: diff_physics.h:530
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Definition: system.C:1082
void add_dot_var_dirichlet_bcs(unsigned int var_idx, unsigned int dot_var_idx)
Definition: diff_system.C:228
std::map< unsigned int, unsigned int > _second_order_dot_vars
Definition: diff_physics.h:571
const Variable & variable(unsigned int var) const
Definition: system.h:2114
virtual void time_evolving(unsigned int var)
Definition: diff_physics.h:249
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 1052 of file system.C.

References libMesh::System::add_vector().

Referenced by libMesh::ImplicitSystem::assemble_residual_derivatives(), and libMesh::System::project_solution_on_reinit().

1053 {
1054  std::ostringstream sensitivity_rhs_name;
1055  sensitivity_rhs_name << "sensitivity_rhs" << i;
1056 
1057  return this->add_vector(sensitivity_rhs_name.str(), false);
1058 }
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Definition: system.C:662
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 907 of file system.C.

References libMesh::System::add_vector().

Referenced by libMesh::System::project_solution_on_reinit(), and libMesh::ImplicitSystem::sensitivity_solve().

908 {
909  std::ostringstream sensitivity_name;
910  sensitivity_name << "sensitivity_solution" << i;
911 
912  return this->add_vector(sensitivity_name.str());
913 }
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Definition: system.C:662
unsigned int libMesh::System::add_variable ( const std::string &  var,
const FEType type,
const std::set< subdomain_id_type > *const  active_subdomains = libmesh_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 1082 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_nullptr, 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::project_solution_on_reinit().

1085 {
1086  libmesh_assert(!this->is_initialized());
1087 
1088  // Make sure the variable isn't there already
1089  // or if it is, that it's the type we want
1090  for (unsigned int v=0; v<this->n_vars(); v++)
1091  if (this->variable_name(v) == var)
1092  {
1093  if (this->variable_type(v) == type)
1094  return _variables[v].number();
1095 
1096  libmesh_error_msg("ERROR: incompatible variable " << var << " has already been added for this system!");
1097  }
1098 
1099  // Optimize for VariableGroups here - if the user is adding multiple
1100  // variables of the same FEType and subdomain restriction, catch
1101  // that here and add them as members of the same VariableGroup.
1102  //
1103  // start by setting this flag to whatever the user has requested
1104  // and then consider the conditions which should negate it.
1105  bool should_be_in_vg = this->identify_variable_groups();
1106 
1107  // No variable groups, nothing to add to
1108  if (!this->n_variable_groups())
1109  should_be_in_vg = false;
1110 
1111  else
1112  {
1113  VariableGroup & vg(_variable_groups.back());
1114 
1115  // get a pointer to their subdomain restriction, if any.
1116  const std::set<subdomain_id_type> * const
1117  their_active_subdomains (vg.implicitly_active() ?
1118  libmesh_nullptr : &vg.active_subdomains());
1119 
1120  // Different types?
1121  if (vg.type() != type)
1122  should_be_in_vg = false;
1123 
1124  // they are restricted, we aren't?
1125  if (their_active_subdomains && !active_subdomains)
1126  should_be_in_vg = false;
1127 
1128  // they aren't restricted, we are?
1129  if (!their_active_subdomains && active_subdomains)
1130  should_be_in_vg = false;
1131 
1132  if (their_active_subdomains && active_subdomains)
1133  // restricted to different sets?
1134  if (*their_active_subdomains != *active_subdomains)
1135  should_be_in_vg = false;
1136 
1137  // OK, after all that, append the variable to the vg if none of the conditions
1138  // were violated
1139  if (should_be_in_vg)
1140  {
1141  const unsigned short curr_n_vars = cast_int<unsigned short>
1142  (this->n_vars());
1143 
1144  vg.append (var);
1145 
1146  _variables.push_back(vg(vg.n_variables()-1));
1147  _variable_numbers[var] = curr_n_vars;
1148  return curr_n_vars;
1149  }
1150  }
1151 
1152  // otherwise, fall back to adding a single variable group
1153  return this->add_variables (std::vector<std::string>(1, var),
1154  type,
1155  active_subdomains);
1156 }
std::map< std::string, unsigned short int > _variable_numbers
Definition: system.h:1903
const std::string & variable_name(const unsigned int i) const
Definition: system.h:2134
const class libmesh_nullptr_t libmesh_nullptr
std::vector< Variable > _variables
Definition: system.h:1892
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
unsigned int n_variable_groups() const
Definition: system.h:2094
std::vector< VariableGroup > _variable_groups
Definition: system.h:1897
bool is_initialized()
Definition: system.h:2070
bool identify_variable_groups() const
Definition: system.h:2182
unsigned int number() const
Definition: system.h:2006
unsigned int add_variables(const std::vector< std::string > &vars, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Definition: system.C:1172
unsigned int n_vars() const
Definition: system.h:2086
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 = libmesh_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 1160 of file system.C.

References libMesh::System::add_variable().

1164 {
1165  return this->add_variable(var,
1166  FEType(order, family),
1167  active_subdomains);
1168 }
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Definition: system.C:1082
unsigned int libMesh::System::add_variables ( const std::vector< std::string > &  vars,
const FEType type,
const std::set< subdomain_id_type > *const  active_subdomains = libmesh_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 1172 of file system.C.

References libMesh::System::_variable_groups, libMesh::System::_variable_numbers, libMesh::System::_variables, libMesh::System::is_initialized(), libmesh_nullptr, 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(), libMesh::System::add_variables(), and libMesh::System::project_solution_on_reinit().

1175 {
1176  libmesh_assert(!this->is_initialized());
1177 
1178  // Make sure the variable isn't there already
1179  // or if it is, that it's the type we want
1180  for (std::size_t ov=0; ov<vars.size(); ov++)
1181  for (unsigned int v=0; v<this->n_vars(); v++)
1182  if (this->variable_name(v) == vars[ov])
1183  {
1184  if (this->variable_type(v) == type)
1185  return _variables[v].number();
1186 
1187  libmesh_error_msg("ERROR: incompatible variable " << vars[ov] << " has already been added for this system!");
1188  }
1189 
1190  const unsigned short curr_n_vars = cast_int<unsigned short>
1191  (this->n_vars());
1192 
1193  const unsigned int next_first_component = this->n_components();
1194 
1195  // Add the variable group to the list
1196  _variable_groups.push_back((active_subdomains == libmesh_nullptr) ?
1197  VariableGroup(this, vars, curr_n_vars,
1198  next_first_component, type) :
1199  VariableGroup(this, vars, curr_n_vars,
1200  next_first_component, type, *active_subdomains));
1201 
1202  const VariableGroup & vg (_variable_groups.back());
1203 
1204  // Add each component of the group individually
1205  for (std::size_t v=0; v<vars.size(); v++)
1206  {
1207  _variables.push_back (vg(v));
1208  _variable_numbers[vars[v]] = cast_int<unsigned short>
1209  (curr_n_vars+v);
1210  }
1211 
1212  libmesh_assert_equal_to ((curr_n_vars+vars.size()), this->n_vars());
1213 
1214  // BSK - Defer this now to System::init_data() so we can detect
1215  // VariableGroups 12/28/2012
1216  // // Add the variable group to the _dof_map
1217  // _dof_map->add_variable_group (vg);
1218 
1219  // Return the number of the new variable
1220  return cast_int<unsigned int>(curr_n_vars+vars.size()-1);
1221 }
std::map< std::string, unsigned short int > _variable_numbers
Definition: system.h:1903
const std::string & variable_name(const unsigned int i) const
Definition: system.h:2134
const class libmesh_nullptr_t libmesh_nullptr
std::vector< Variable > _variables
Definition: system.h:1892
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
unsigned int n_components() const
Definition: system.h:2102
std::vector< VariableGroup > _variable_groups
Definition: system.h:1897
bool is_initialized()
Definition: system.h:2070
unsigned int number() const
Definition: system.h:2006
unsigned int n_vars() const
Definition: system.h:2086
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 = libmesh_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 1225 of file system.C.

References libMesh::System::add_variables().

1229 {
1230  return this->add_variables(vars,
1231  FEType(order, family),
1232  active_subdomains);
1233 }
unsigned int add_variables(const std::vector< std::string > &vars, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Definition: system.C:1172
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 662 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::FrequencySystem::set_frequencies(), libMesh::FrequencySystem::set_frequencies_by_range(), and libMesh::FrequencySystem::set_frequencies_by_steps().

665 {
666  // Return the vector if it is already there.
667  if (this->have_vector(vec_name))
668  return *(_vectors[vec_name]);
669 
670  // Otherwise build the vector
671  NumericVector<Number> * buf = NumericVector<Number>::build(this->comm()).release();
672  _vectors.insert (std::make_pair (vec_name, buf));
673  _vector_projections.insert (std::make_pair (vec_name, projections));
674 
675  _vector_types.insert (std::make_pair (vec_name, type));
676 
677  // Vectors are primal by default
678  _vector_is_adjoint.insert (std::make_pair (vec_name, -1));
679 
680  // Initialize it if necessary
681  if (_is_initialized)
682  {
683  if (type == GHOSTED)
684  {
685 #ifdef LIBMESH_ENABLE_GHOSTED
686  buf->init (this->n_dofs(), this->n_local_dofs(),
687  _dof_map->get_send_list(), false,
688  GHOSTED);
689 #else
690  libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
691 #endif
692  }
693  else
694  buf->init (this->n_dofs(), this->n_local_dofs(), false, type);
695  }
696 
697  return *buf;
698 }
std::map< std::string, ParallelType > _vector_types
Definition: system.h:1933
bool _is_initialized
Definition: system.h:1952
std::unique_ptr< DofMap > _dof_map
Definition: system.h:1865
bool have_vector(const std::string &vec_name) const
Definition: system.h:2206
std::map< std::string, int > _vector_is_adjoint
Definition: system.h:1928
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
dof_id_type n_local_dofs() const
Definition: system.C:185
std::map< std::string, NumericVector< Number > * > _vectors
Definition: system.h:1916
const Parallel::Communicator & comm() const
std::map< std::string, bool > _vector_projections
Definition: system.h:1922
dof_id_type n_dofs() const
Definition: system.C:148
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 990 of file system.C.

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

Referenced by libMesh::System::project_solution_on_reinit(), and libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve().

991 {
992  std::ostringstream adjoint_name;
993  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
994 
995  NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
996  this->set_vector_as_adjoint(adjoint_name.str(), i);
997  return returnval;
998 }
void set_vector_as_adjoint(const std::string &vec_name, int qoi_num)
Definition: system.C:886
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Definition: system.C:662
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 937 of file system.C.

References libMesh::System::add_vector().

Referenced by libMesh::System::project_solution_on_reinit(), and libMesh::ImplicitSystem::weighted_sensitivity_solve().

938 {
939  return this->add_vector("weighted_sensitivity_solution");
940 }
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Definition: system.C:662
void libMesh::ImplicitSystem::adjoint_qoi_parameter_sensitivity ( const QoISet qoi_indices,
const ParameterVector parameters,
SensitivityData sensitivities 
)
virtualinherited

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 704 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::qoi, libMesh::Real, libMesh::ParameterVector::size(), and libMesh::TOLERANCE.

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

707 {
708  ParameterVector & parameters =
709  const_cast<ParameterVector &>(parameters_in);
710 
711  const unsigned int Np = cast_int<unsigned int>
712  (parameters.size());
713  const unsigned int Nq = cast_int<unsigned int>
714  (qoi.size());
715 
716  // An introduction to the problem:
717  //
718  // Residual R(u(p),p) = 0
719  // partial R / partial u = J = system matrix
720  //
721  // This implies that:
722  // d/dp(R) = 0
723  // (partial R / partial p) +
724  // (partial R / partial u) * (partial u / partial p) = 0
725 
726  // We first do an adjoint solve:
727  // J^T * z = (partial q / partial u)
728  // if we havent already or dont have an initial condition for the adjoint
729  if (!this->is_adjoint_already_solved())
730  {
731  this->adjoint_solve(qoi_indices);
732  }
733 
734  this->assemble_residual_derivatives(parameters_in);
735 
736  // Get ready to fill in sensitivities:
737  sensitivities.allocate_data(qoi_indices, *this, parameters);
738 
739  // We use the identities:
740  // dq/dp = (partial q / partial p) + (partial q / partial u) *
741  // (partial u / partial p)
742  // dq/dp = (partial q / partial p) + (J^T * z) *
743  // (partial u / partial p)
744  // dq/dp = (partial q / partial p) + z * J *
745  // (partial u / partial p)
746 
747  // Leading to our final formula:
748  // dq/dp = (partial q / partial p) - z * (partial R / partial p)
749 
750  // In the case of adjoints with heterogenous Dirichlet boundary
751  // function phi, where
752  // q := S(u) - R(u,phi)
753  // the final formula works out to:
754  // dq/dp = (partial S / partial p) - z * (partial R / partial p)
755  // Because we currently have no direct access to
756  // (partial S / partial p), we use the identity
757  // (partial S / partial p) = (partial q / partial p) +
758  // phi * (partial R / partial p)
759  // to derive an equivalent equation:
760  // dq/dp = (partial q / partial p) - (z-phi) * (partial R / partial p)
761 
762  // Since z-phi degrees of freedom are zero for constrained indices,
763  // we can use the same constrained -(partial R / partial p) that we
764  // use for forward sensitivity solves, taking into account the
765  // differing sign convention.
766  //
767  // Since that vector is constrained, its constrained indices are
768  // zero, so its product with phi is zero, so we can neglect the
769  // evaluation of phi terms.
770 
771  for (unsigned int j=0; j != Np; ++j)
772  {
773  // We currently get partial derivatives via central differencing
774 
775  // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
776  // (partial R / partial p) ~= (rhs(p+dp) - rhs(p-dp))/(2*dp)
777 
778  Number old_parameter = *parameters[j];
779 
780  const Real delta_p =
781  TOLERANCE * std::max(std::abs(old_parameter), 1e-3);
782 
783  *parameters[j] = old_parameter - delta_p;
784  this->assemble_qoi(qoi_indices);
785  std::vector<Number> qoi_minus = this->qoi;
786 
787  NumericVector<Number> & neg_partialR_partialp = this->get_sensitivity_rhs(j);
788 
789  *parameters[j] = old_parameter + delta_p;
790  this->assemble_qoi(qoi_indices);
791  std::vector<Number> & qoi_plus = this->qoi;
792 
793  std::vector<Number> partialq_partialp(Nq, 0);
794  for (unsigned int i=0; i != Nq; ++i)
795  if (qoi_indices.has_index(i))
796  partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
797 
798  // Don't leave the parameter changed
799  *parameters[j] = old_parameter;
800 
801  for (unsigned int i=0; i != Nq; ++i)
802  if (qoi_indices.has_index(i))
803  sensitivities[i][j] = partialq_partialp[i] +
804  neg_partialR_partialp.dot(this->get_adjoint_solution(i));
805  }
806 
807  // All parameters have been reset.
808  // Reset the original qoi.
809 
810  this->assemble_qoi(qoi_indices);
811 }
double abs(double a)
virtual std::pair< unsigned int, Real > adjoint_solve(const QoISet &qoi_indices=QoISet()) libmesh_override
NumericVector< Number > & get_sensitivity_rhs(unsigned int i=0)
Definition: system.C:1062
static const Real TOLERANCE
long double max(long double a, double b)
std::vector< Number > qoi
Definition: system.h:1539
bool is_adjoint_already_solved() const
Definition: system.h:374
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void assemble_qoi(const QoISet &qoi_indices=QoISet()) libmesh_override
virtual void assemble_residual_derivatives(const ParameterVector &parameters) libmesh_override
std::pair< unsigned int, Real > libMesh::DifferentiableSystem::adjoint_solve ( const QoISet qoi_indices = QoISet())
virtualinherited

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().

165 {
166  // Get the time solver object associated with the system, and tell it that
167  // we are solving the adjoint problem
168  this->get_time_solver().set_is_adjoint(true);
169 
170  return this->ImplicitSystem::adjoint_solve(qoi_indices);
171 }
virtual std::pair< unsigned int, Real > adjoint_solve(const QoISet &qoi_indices=QoISet()) libmesh_override
void set_is_adjoint(bool _is_adjoint_value)
Definition: time_solver.h:240
TimeSolver & get_time_solver()
Definition: diff_system.h:415
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().

937 {
938  // Solve for the updated tangent du1/ds, d(lambda1)/ds
939  solve_tangent();
940 
941  // Advance the solution and the parameter to the next value.
942  update_solution();
943 }
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().

1376 {
1377  if (predictor == Euler)
1378  {
1379  // 1.) Euler Predictor
1380  // Predict next the solution
1381  solution->add(ds_current, *du_ds);
1382  solution->close();
1383 
1384  // Predict next parameter value
1386  }
1387 
1388 
1389  else if (predictor == AB2)
1390  {
1391  // 2.) 2nd-order explicit AB predictor
1392  libmesh_assert_not_equal_to (previous_ds, 0.0);
1393  const Real stepratio = ds_current/previous_ds;
1394 
1395  // Build up next solution value.
1396  solution->add( 0.5*ds_current*(2.+stepratio), *du_ds);
1397  solution->add(-0.5*ds_current*stepratio , *previous_du_ds);
1398  solution->close();
1399 
1400  // Next parameter value
1402  0.5*ds_current*((2.+stepratio)*dlambda_ds -
1403  stepratio*previous_dlambda_ds);
1404  }
1405 
1406  else
1407  libmesh_error_msg("Unknown predictor!");
1408 }
NumericVector< Number > * previous_du_ds
NumericVector< Number > * du_ds
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1509
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void libMesh::DifferentiableSystem::assemble ( )
virtualinherited

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().

146 {
147  this->assembly(true, true);
148 }
virtual void assembly(bool get_residual, bool get_jacobian, bool apply_heterogeneous_constraints=false, bool apply_no_constraints=false) libmesh_override=0
void libMesh::FEMSystem::assemble_qoi ( const QoISet indices = QoISet())
virtualinherited

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 1119 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::DifferentiableQoI::parallel_op(), libMesh::Threads::parallel_reduce(), libMesh::System::qoi, libMesh::StoredRange< iterator_type, object_type >::reset(), and libMesh::System::update().

1120 {
1121  LOG_SCOPE("assemble_qoi()", "FEMSystem");
1122 
1123  const MeshBase & mesh = this->get_mesh();
1124 
1125  this->update();
1126 
1127  const unsigned int Nq = cast_int<unsigned int>(qoi.size());
1128 
1129  // the quantity of interest is assumed to be a sum of element and
1130  // side terms
1131  for (unsigned int i=0; i != Nq; ++i)
1132  if (qoi_indices.has_index(i))
1133  qoi[i] = 0;
1134 
1135  // Create a non-temporary qoi_contributions object, so we can query
1136  // its results after the reduction
1137  QoIContributions qoi_contributions(*this, *(this->diff_qoi), qoi_indices);
1138 
1139  // Loop over every active mesh element on this processor
1140  Threads::parallel_reduce(elem_range.reset(mesh.active_local_elements_begin(),
1141  mesh.active_local_elements_end()),
1142  qoi_contributions);
1143 
1144  this->diff_qoi->parallel_op( this->comm(), this->qoi, qoi_contributions.qoi, qoi_indices );
1145 }
MeshBase & mesh
Base class for Mesh.
Definition: mesh_base.h:69
DifferentiableQoI * diff_qoi
Definition: diff_system.h:378
virtual void parallel_op(const Parallel::Communicator &communicator, std::vector< Number > &sys_qoi, std::vector< Number > &local_qoi, const QoISet &qoi_indices)
const MeshBase & get_mesh() const
Definition: system.h:2014
std::vector< Number > qoi
Definition: system.h:1539
virtual element_iterator active_local_elements_begin()=0
virtual void update()
Definition: system.C:406
const Parallel::Communicator & comm() const
void parallel_reduce(const Range &range, Body &body)
Definition: threads_none.h:101
virtual element_iterator active_local_elements_end()=0
void libMesh::FEMSystem::assemble_qoi_derivative ( const QoISet qoi_indices = QoISet(),
bool  include_liftfunc = true,
bool  apply_constraints = true 
)
virtualinherited

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 1149 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::Threads::parallel_for(), libMesh::System::qoi, libMesh::StoredRange< iterator_type, object_type >::reset(), libMesh::System::update(), and libMesh::NumericVector< T >::zero().

1152 {
1153  LOG_SCOPE("assemble_qoi_derivative()", "FEMSystem");
1154 
1155  const MeshBase & mesh = this->get_mesh();
1156 
1157  this->update();
1158 
1159  // The quantity of interest derivative assembly accumulates on
1160  // initially zero vectors
1161  for (std::size_t i=0; i != qoi.size(); ++i)
1162  if (qoi_indices.has_index(i))
1163  this->add_adjoint_rhs(i).zero();
1164 
1165  // Loop over every active mesh element on this processor
1166  Threads::parallel_for (elem_range.reset(mesh.active_local_elements_begin(),
1167  mesh.active_local_elements_end()),
1168  QoIDerivativeContributions(*this, qoi_indices,
1169  *(this->diff_qoi),
1170  include_liftfunc,
1171  apply_constraints));
1172 
1173  for (std::size_t i=0; i != qoi.size(); ++i)
1174  if (qoi_indices.has_index(i))
1175  this->diff_qoi->finalize_derivative(this->get_adjoint_rhs(i),i);
1176 }
virtual void finalize_derivative(NumericVector< Number > &derivatives, std::size_t qoi_index)
void parallel_for(const Range &range, const Body &body)
Definition: threads_none.h:73
MeshBase & mesh
Base class for Mesh.
Definition: mesh_base.h:69
virtual void zero()=0
DifferentiableQoI * diff_qoi
Definition: diff_system.h:378
const MeshBase & get_mesh() const
Definition: system.h:2014
std::vector< Number > qoi
Definition: system.h:1539
virtual element_iterator active_local_elements_begin()=0
NumericVector< Number > & add_adjoint_rhs(unsigned int i=0)
Definition: system.C:1022
virtual void update()
Definition: system.C:406
virtual element_iterator active_local_elements_end()=0
bool has_index(unsigned int) const
Definition: qoi_set.h:221
void libMesh::ImplicitSystem::assemble_residual_derivatives ( const ParameterVector parameters)
virtualinherited

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 661 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(), libMesh::ImplicitSystem::assembly(), and libMesh::ImplicitSystem::sensitivity_solve().

662 {
663  ParameterVector & parameters =
664  const_cast<ParameterVector &>(parameters_in);
665 
666  const unsigned int Np = cast_int<unsigned int>
667  (parameters.size());
668 
669  for (unsigned int p=0; p != Np; ++p)
670  {
671  NumericVector<Number> & sensitivity_rhs = this->add_sensitivity_rhs(p);
672 
673  // Approximate -(partial R / partial p) by
674  // (R(p-dp) - R(p+dp)) / (2*dp)
675 
676  Number old_parameter = *parameters[p];
677 
678  const Real delta_p =
679  TOLERANCE * std::max(std::abs(old_parameter), 1e-3);
680 
681  *parameters[p] -= delta_p;
682 
683  // this->assembly(true, false, true);
684  this->assembly(true, false, false);
685  this->rhs->close();
686  sensitivity_rhs = *this->rhs;
687 
688  *parameters[p] = old_parameter + delta_p;
689 
690  // this->assembly(true, false, true);
691  this->assembly(true, false, false);
692  this->rhs->close();
693 
694  sensitivity_rhs -= *this->rhs;
695  sensitivity_rhs /= (2*delta_p);
696  sensitivity_rhs.close();
697 
698  *parameters[p] = old_parameter;
699  }
700 }
double abs(double a)
NumericVector< Number > * rhs
NumericVector< Number > & add_sensitivity_rhs(unsigned int i=0)
Definition: system.C:1052
static const Real TOLERANCE
long double max(long double a, double b)
virtual void assembly(bool, bool, bool=false, bool=false)
virtual void close()=0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void libMesh::FEMSystem::assembly ( bool  get_residual,
bool  get_jacobian,
bool  apply_heterogeneous_constraints = false,
bool  apply_no_constraints = false 
)
virtualinherited

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 852 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::DenseMatrix< T >::l1_norm(), libMesh::SparseMatrix< T >::l1_norm(), libmesh_nullptr, 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().

855 {
856  libmesh_assert(get_residual || get_jacobian);
857 
858  // Log residual and jacobian and combined performance separately
859 #ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
860  const char * log_name;
861  if (get_residual && get_jacobian)
862  log_name = "assembly()";
863  else if (get_residual)
864  log_name = "assembly(get_residual)";
865  else
866  log_name = "assembly(get_jacobian)";
867 
868  LOG_SCOPE(log_name, "FEMSystem");
869 #endif
870 
871  const MeshBase & mesh = this->get_mesh();
872 
874  {
875  this->solution->close();
876 
877  std::streamsize old_precision = libMesh::out.precision();
879  libMesh::out << "|U| = "
880  << this->solution->l1_norm()
881  << std::endl;
882  libMesh::out.precision(old_precision);
883  }
884  if (print_solutions)
885  {
886  std::streamsize old_precision = libMesh::out.precision();
888  libMesh::out << "U = [" << *(this->solution)
889  << "];" << std::endl;
890  libMesh::out.precision(old_precision);
891  }
892 
893  // Is this definitely necessary? [RHS]
894  // Yes. [RHS 2012]
895  if (get_jacobian)
896  matrix->zero();
897  if (get_residual)
898  rhs->zero();
899 
900  // Stupid C++ lets you set *Real* verify_analytic_jacobians = true!
901  if (verify_analytic_jacobians > 0.5)
902  {
903  libMesh::err << "WARNING! verify_analytic_jacobians was set "
904  << "to absurdly large value of "
905  << verify_analytic_jacobians << std::endl;
906  libMesh::err << "Resetting to 1e-6!" << std::endl;
908  }
909 
910  // In time-dependent problems, the nonlinear function we're trying
911  // to solve at each timestep may depend on the particular solver
912  // we're using
913  libmesh_assert(time_solver.get());
914 
915  // Build the residual and jacobian contributions on every active
916  // mesh element on this processor
918  (elem_range.reset(mesh.active_local_elements_begin(),
920  AssemblyContributions(*this, get_residual, get_jacobian,
921  apply_heterogeneous_constraints,
922  apply_no_constraints));
923 
924  // Check and see if we have SCALAR variables
925  bool have_scalar = false;
926  for (unsigned int i=0; i != this->n_variable_groups(); ++i)
927  {
928  if (this->variable_group(i).type().family == SCALAR)
929  {
930  have_scalar = true;
931  break;
932  }
933  }
934 
935  // SCALAR dofs are stored on the last processor, so we'll evaluate
936  // their equation terms there and only if we have a SCALAR variable
937  if (this->processor_id() == (this->n_processors()-1) && have_scalar)
938  {
939  std::unique_ptr<DiffContext> con = this->build_context();
940  FEMContext & _femcontext = cast_ref<FEMContext &>(*con);
941  this->init_context(_femcontext);
942  _femcontext.pre_fe_reinit(*this, libmesh_nullptr);
943 
944  bool jacobian_computed =
945  this->time_solver->nonlocal_residual(get_jacobian, _femcontext);
946 
947  // Nonlocal residuals are likely to be length 0, in which case we
948  // don't need to do any more. And we shouldn't try to do any
949  // more; lots of DenseVector/DenseMatrix code assumes rank>0.
950  if (_femcontext.get_elem_residual().size())
951  {
952  // Compute a numeric jacobian if we have to
953  if (get_jacobian && !jacobian_computed)
954  {
955  // Make sure we didn't compute a jacobian and lie about it
956  libmesh_assert_equal_to (_femcontext.get_elem_jacobian().l1_norm(), 0.0);
957  // Logging of numerical jacobians is done separately
958  this->numerical_nonlocal_jacobian(_femcontext);
959  }
960 
961  // Compute a numeric jacobian if we're asked to verify the
962  // analytic jacobian we got
963  if (get_jacobian && jacobian_computed &&
964  this->verify_analytic_jacobians != 0.0)
965  {
966  DenseMatrix<Number> analytic_jacobian(_femcontext.get_elem_jacobian());
967 
968  _femcontext.get_elem_jacobian().zero();
969  // Logging of numerical jacobians is done separately
970  this->numerical_nonlocal_jacobian(_femcontext);
971 
972  Real analytic_norm = analytic_jacobian.l1_norm();
973  Real numerical_norm = _femcontext.get_elem_jacobian().l1_norm();
974 
975  // If we can continue, we'll probably prefer the analytic jacobian
976  analytic_jacobian.swap(_femcontext.get_elem_jacobian());
977 
978  // The matrix "analytic_jacobian" will now hold the error matrix
979  analytic_jacobian.add(-1.0, _femcontext.get_elem_jacobian());
980  Real error_norm = analytic_jacobian.l1_norm();
981 
982  Real relative_error = error_norm /
983  std::max(analytic_norm, numerical_norm);
984 
985  if (relative_error > this->verify_analytic_jacobians)
986  {
987  libMesh::err << "Relative error " << relative_error
988  << " detected in analytic jacobian on nonlocal dofs!"
989  << std::endl;
990 
991  std::streamsize old_precision = libMesh::out.precision();
993  libMesh::out << "J_analytic nonlocal = "
994  << _femcontext.get_elem_jacobian() << std::endl;
995  analytic_jacobian.add(1.0, _femcontext.get_elem_jacobian());
996  libMesh::out << "J_numeric nonlocal = "
997  << analytic_jacobian << std::endl;
998 
999  libMesh::out.precision(old_precision);
1000 
1001  libmesh_error_msg("Relative error too large, exiting!");
1002  }
1003  }
1004 
1005  add_element_system
1006  (*this, get_residual, get_jacobian,
1007  apply_heterogeneous_constraints, apply_no_constraints, _femcontext);
1008  }
1009  }
1010 
1011  if (get_residual && (print_residual_norms || print_residuals))
1012  this->rhs->close();
1013  if (get_residual && print_residual_norms)
1014  {
1015  std::streamsize old_precision = libMesh::out.precision();
1016  libMesh::out.precision(16);
1017  libMesh::out << "|F| = " << this->rhs->l1_norm() << std::endl;
1018  libMesh::out.precision(old_precision);
1019  }
1020  if (get_residual && print_residuals)
1021  {
1022  std::streamsize old_precision = libMesh::out.precision();
1023  libMesh::out.precision(16);
1024  libMesh::out << "F = [" << *(this->rhs) << "];" << std::endl;
1025  libMesh::out.precision(old_precision);
1026  }
1027 
1028  if (get_jacobian && (print_jacobian_norms || print_jacobians))
1029  this->matrix->close();
1030  if (get_jacobian && print_jacobian_norms)
1031  {
1032  std::streamsize old_precision = libMesh::out.precision();
1033  libMesh::out.precision(16);
1034  libMesh::out << "|J| = " << this->matrix->l1_norm() << std::endl;
1035  libMesh::out.precision(old_precision);
1036  }
1037  if (get_jacobian && print_jacobians)
1038  {
1039  std::streamsize old_precision = libMesh::out.precision();
1040  libMesh::out.precision(16);
1041  libMesh::out << "J = [" << *(this->matrix) << "];" << std::endl;
1042  libMesh::out.precision(old_precision);
1043  }
1044 }
FEFamily family
Definition: fe_type.h:204
const FEType & type() const
Definition: variable.h:119
void numerical_nonlocal_jacobian(FEMContext &context) const
Definition: fem_system.C:1305
Real verify_analytic_jacobians
Definition: fem_system.h:209
virtual void pre_fe_reinit(const System &, const Elem *e)
Definition: fem_context.C:1553
virtual void zero() libmesh_override
Definition: dense_matrix.h:792
void parallel_for(const Range &range, const Body &body)
Definition: threads_none.h:73
processor_id_type n_processors() const
std::unique_ptr< TimeSolver > time_solver
Definition: diff_system.h:234
MeshBase & mesh
Real l1_norm() const
Definition: dense_matrix.h:990
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
const VariableGroup & variable_group(unsigned int vg) const
Definition: system.h:2124
const DenseVector< Number > & get_elem_residual() const
Definition: diff_context.h:249
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
virtual Real l1_norm() const =0
long double max(long double a, double b)
Base class for Mesh.
Definition: mesh_base.h:69
virtual void zero()=0
const DenseMatrix< Number > & get_elem_jacobian() const
Definition: diff_context.h:283
unsigned int n_variable_groups() const
Definition: system.h:2094
const MeshBase & get_mesh() const
Definition: system.h:2014
virtual void zero()=0
virtual void init_context(DiffContext &) libmesh_override
Definition: fem_system.C:1337
virtual element_iterator active_local_elements_begin()=0
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseMatrix< T3 > &mat)
Definition: dense_matrix.h:883
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1509
OStreamProxy err(std::cerr)
virtual std::unique_ptr< DiffContext > build_context() libmesh_override
Definition: fem_system.C:1313
virtual void close()=0
virtual void close()=0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
SparseMatrix< Number > * matrix
OStreamProxy out(std::cout)
virtual element_iterator active_local_elements_end()=0
virtual Real l1_norm() const =0
std::streamsize precision() const
processor_id_type processor_id() const
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 1764 of file system.C.

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

Referenced by libMesh::System::read_parallel_data().

1766 {
1767  libmesh_assert(fptr);
1768 
1770  {
1771  libmesh_here();
1772  libMesh::out << "WARNING: Cannot specify both assembly function and object!"
1773  << std::endl;
1774 
1776  }
1777 
1779 }
Assembly * _assemble_system_object
Definition: system.h:1822
const class libmesh_nullptr_t libmesh_nullptr
void(* _assemble_system_function)(EquationSystems &es, const std::string &name)
Definition: system.h:1816
OStreamProxy out(std::cout)
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 1783 of file system.C.

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

Referenced by libMesh::System::read_parallel_data().

1784 {
1786  {
1787  libmesh_here();
1788  libMesh::out << "WARNING: Cannot specify both assembly object and function!"
1789  << std::endl;
1790 
1792  }
1793 
1794  _assemble_system_object = &assemble_in;
1795 }
Assembly * _assemble_system_object
Definition: system.h:1822
const class libmesh_nullptr_t libmesh_nullptr
void(* _assemble_system_function)(EquationSystems &es, const std::string &name)
Definition: system.h:1816
OStreamProxy out(std::cout)
void libMesh::System::attach_constraint_function ( void   fptrEquationSystems &es, const std::string &name)
inherited

Register a user function for imposing constraints.

Definition at line 1799 of file system.C.

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

Referenced by libMesh::System::read_parallel_data().

1801 {
1802  libmesh_assert(fptr);
1803 
1805  {
1806  libmesh_here();
1807  libMesh::out << "WARNING: Cannot specify both constraint function and object!"
1808  << std::endl;
1809 
1811  }
1812 
1814 }
void(* _constrain_system_function)(EquationSystems &es, const std::string &name)
Definition: system.h:1827
Constraint * _constrain_system_object
Definition: system.h:1833
const class libmesh_nullptr_t libmesh_nullptr
OStreamProxy out(std::cout)
void libMesh::System::attach_constraint_object ( System::Constraint constrain)
inherited

Register a user object for imposing constraints.

Definition at line 1818 of file system.C.

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

Referenced by libMesh::System::read_parallel_data().

1819 {
1821  {
1822  libmesh_here();
1823  libMesh::out << "WARNING: Cannot specify both constraint object and function!"
1824  << std::endl;
1825 
1827  }
1828 
1829  _constrain_system_object = &constrain;
1830 }
void(* _constrain_system_function)(EquationSystems &es, const std::string &name)
Definition: system.h:1827
Constraint * _constrain_system_object
Definition: system.h:1833
const class libmesh_nullptr_t libmesh_nullptr
OStreamProxy out(std::cout)
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 1729 of file system.C.

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

Referenced by libMesh::System::read_parallel_data().

1731 {
1732  libmesh_assert(fptr);
1733 
1735  {
1736  libmesh_here();
1737  libMesh::out << "WARNING: Cannot specify both initialization function and object!"
1738  << std::endl;
1739 
1741  }
1742 
1743  _init_system_function = fptr;
1744 }
const class libmesh_nullptr_t libmesh_nullptr
Initialization * _init_system_object
Definition: system.h:1811
void(* _init_system_function)(EquationSystems &es, const std::string &name)
Definition: system.h:1805
OStreamProxy out(std::cout)
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 1748 of file system.C.

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

Referenced by libMesh::System::read_parallel_data().

1749 {
1751  {
1752  libmesh_here();
1753  libMesh::out << "WARNING: Cannot specify both initialization object and function!"
1754  << std::endl;
1755 
1757  }
1758 
1759  _init_system_object = &init_in;
1760 }
const class libmesh_nullptr_t libmesh_nullptr
Initialization * _init_system_object
Definition: system.h:1811
void(* _init_system_function)(EquationSystems &es, const std::string &name)
Definition: system.h:1805
OStreamProxy out(std::cout)
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(), libMesh::DifferentiablePhysics::init_physics(), and libMesh::DifferentiableSystem::swap_physics().

198  { this->_diff_physics = (physics_in->clone_physics()).release();
199  this->_diff_physics->init_physics(*this);}
DifferentiablePhysics * _diff_physics
Definition: diff_system.h:371
virtual void init_physics(const System &sys)
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.

226  { this->diff_qoi = (qoi_in->clone()).release();
227  // User needs to resize qoi system qoi accordingly
228  this->diff_qoi->init_qoi( this->qoi );}
DifferentiableQoI * diff_qoi
Definition: diff_system.h:378
std::vector< Number > qoi
Definition: system.h:1539
virtual void init_qoi(std::vector< Number > &)
Definition: diff_qoi.h:69
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 1870 of file system.C.

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

Referenced by libMesh::System::read_parallel_data().

1872 {
1873  libmesh_assert(fptr);
1874 
1876  {
1877  libmesh_here();
1878  libMesh::out << "WARNING: Cannot specify both QOI derivative function and object!"
1879  << std::endl;
1880 
1882  }
1883 
1885 }
const class libmesh_nullptr_t libmesh_nullptr
QOIDerivative * _qoi_evaluate_derivative_object
Definition: system.h:1859
OStreamProxy out(std::cout)
void(* _qoi_evaluate_derivative_function)(EquationSystems &es, const std::string &name, const QoISet &qoi_indices, bool include_liftfunc, bool apply_constraints)
Definition: system.h:1850
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 1889 of file system.C.

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

Referenced by libMesh::System::read_parallel_data().

1890 {
1892  {
1893  libmesh_here();
1894  libMesh::out << "WARNING: Cannot specify both QOI derivative object and function!"
1895  << std::endl;
1896 
1898  }
1899 
1900  _qoi_evaluate_derivative_object = &qoi_derivative;
1901 }
const class libmesh_nullptr_t libmesh_nullptr
QOIDerivative * _qoi_evaluate_derivative_object
Definition: system.h:1859
OStreamProxy out(std::cout)
void(* _qoi_evaluate_derivative_function)(EquationSystems &es, const std::string &name, const QoISet &qoi_indices, bool include_liftfunc, bool apply_constraints)
Definition: system.h:1850
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 1834 of file system.C.

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

Referenced by libMesh::System::read_parallel_data().

1837 {
1838  libmesh_assert(fptr);
1839 
1841  {
1842  libmesh_here();
1843  libMesh::out << "WARNING: Cannot specify both QOI function and object!"
1844  << std::endl;
1845 
1847  }
1848 
1849  _qoi_evaluate_function = fptr;
1850 }
const class libmesh_nullptr_t libmesh_nullptr
void(* _qoi_evaluate_function)(EquationSystems &es, const std::string &name, const QoISet &qoi_indices)
Definition: system.h:1838
QOI * _qoi_evaluate_object
Definition: system.h:1845
OStreamProxy out(std::cout)
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 1854 of file system.C.

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

Referenced by libMesh::System::read_parallel_data().

1855 {
1857  {
1858  libmesh_here();
1859  libMesh::out << "WARNING: Cannot specify both QOI object and function!"
1860  << std::endl;
1861 
1863  }
1864 
1865  _qoi_evaluate_object = &qoi_in;
1866 }
const class libmesh_nullptr_t libmesh_nullptr
void(* _qoi_evaluate_function)(EquationSystems &es, const std::string &name, const QoISet &qoi_indices)
Definition: system.h:1838
QOI * _qoi_evaluate_object
Definition: system.h:1845
OStreamProxy out(std::cout)
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 = libmesh_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 986 of file system_projection.C.

990 {
991  this->boundary_project_vector(b, variables, *solution, f, g);
992 
993  solution->localize(*current_local_solution);
994 }
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=libmesh_nullptr, int is_adjoint=-1) const
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1509
std::unique_ptr< NumericVector< Number > > current_local_solution
Definition: system.h:1521
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 969 of file system_projection.C.

974 {
975  WrappedFunction<Number> f(*this, fptr, &parameters);
976  WrappedFunction<Gradient> g(*this, gptr, &parameters);
977  this->boundary_project_solution(b, variables, &f, &g);
978 }
void boundary_project_solution(const std::set< boundary_id_type > &b, const std::vector< unsigned int > &variables, FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr)
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 = libmesh_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 1022 of file system_projection.C.

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

1028 {
1029  LOG_SCOPE ("boundary_project_vector()", "System");
1030 
1032  (ConstElemRange (this->get_mesh().active_local_elements_begin(),
1033  this->get_mesh().active_local_elements_end() ),
1034  BoundaryProjectSolution(b, variables, *this, f, g,
1035  this->get_equation_systems().parameters,
1036  new_vector)
1037  );
1038 
1039  // We don't do SCALAR dofs when just projecting the boundary, so
1040  // we're done here.
1041 
1042  new_vector.close();
1043 
1044 #ifdef LIBMESH_ENABLE_CONSTRAINTS
1045  if (is_adjoint == -1)
1046  this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
1047  else if (is_adjoint >= 0)
1049  is_adjoint);
1050 #endif
1051 }
void parallel_for(const Range &range, const Body &body)
Definition: threads_none.h:73
StoredRange< MeshBase::const_element_iterator, const Elem * > ConstElemRange
Definition: elem_range.h:34
const MeshBase & get_mesh() const
Definition: system.h:2014
const DofMap & get_dof_map() const
Definition: system.h:2030
void enforce_constraints_exactly(const System &system, NumericVector< Number > *v=libmesh_nullptr, bool homogeneous=false) const
virtual void close()=0
const EquationSystems & get_equation_systems() const
Definition: system.h:698
void enforce_adjoint_constraints_exactly(NumericVector< Number > &v, unsigned int q) const
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 1004 of file system_projection.C.

1011 {
1012  WrappedFunction<Number> f(*this, fptr, &parameters);
1013  WrappedFunction<Gradient> g(*this, gptr, &parameters);
1014  this->boundary_project_vector(b, variables, new_vector, &f, &g,
1015  is_adjoint);
1016 }
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=libmesh_nullptr, int is_adjoint=-1) const
std::unique_ptr< DiffContext > libMesh::FEMSystem::build_context ( )
virtualinherited

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 1313 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().

1314 {
1315  FEMContext * fc = new FEMContext(*this);
1316 
1317  DifferentiablePhysics * phys = this->get_physics();
1318 
1319  libmesh_assert (phys);
1320 
1321  // If we are solving a moving mesh problem, tell that to the Context
1322  fc->set_mesh_system(phys->get_mesh_system());
1323  fc->set_mesh_x_var(phys->get_mesh_x_var());
1324  fc->set_mesh_y_var(phys->get_mesh_y_var());
1325  fc->set_mesh_z_var(phys->get_mesh_z_var());
1326 
1327  fc->set_deltat_pointer( &deltat );
1328 
1329  // If we are solving the adjoint problem, tell that to the Context
1330  fc->is_adjoint() = this->get_time_solver().is_adjoint();
1331 
1332  return std::unique_ptr<DiffContext>(fc);
1333 }
bool is_adjoint() const
Definition: time_solver.h:233
void set_mesh_z_var(unsigned int z_var)
Definition: fem_context.h:859
void set_mesh_x_var(unsigned int x_var)
Definition: fem_context.h:831
bool is_adjoint() const
Definition: diff_context.h:454
const System * get_mesh_system() const
Definition: diff_physics.h:624
unsigned int get_mesh_x_var() const
Definition: diff_physics.h:636
virtual void set_mesh_system(System *sys)
Definition: fem_context.h:805
const DifferentiablePhysics * get_physics() const
Definition: diff_system.h:182
unsigned int get_mesh_z_var() const
Definition: diff_physics.h:648
unsigned int get_mesh_y_var() const
Definition: diff_physics.h:642
void set_mesh_y_var(unsigned int y_var)
Definition: fem_context.h:845
void set_deltat_pointer(Real *dt)
Definition: diff_context.C:103
TimeSolver & get_time_solver()
Definition: diff_system.h:415
Real libMesh::System::calculate_norm ( const NumericVector< Number > &  v,
unsigned int  var,
FEMNormType  norm_type,
std::set< unsigned int > *  skip_dimensions = libmesh_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 1364 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(), libMesh::UnsteadySolver::du(), and libMesh::System::project_solution_on_reinit().

1368 {
1369  //short circuit to save time
1370  if (norm_type == DISCRETE_L1 ||
1371  norm_type == DISCRETE_L2 ||
1372  norm_type == DISCRETE_L_INF)
1373  return discrete_var_norm(v,var,norm_type);
1374 
1375  // Not a discrete norm
1376  std::vector<FEMNormType> norms(this->n_vars(), L2);
1377  std::vector<Real> weights(this->n_vars(), 0.0);
1378  norms[var] = norm_type;
1379  weights[var] = 1.0;
1380  Real val = this->calculate_norm(v, SystemNorm(norms, weights), skip_dimensions);
1381  return val;
1382 }
Real calculate_norm(const NumericVector< Number > &v, unsigned int var, FEMNormType norm_type, std::set< unsigned int > *skip_dimensions=libmesh_nullptr) const
Definition: system.C:1364
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned int n_vars() const
Definition: system.h:2086
Real discrete_var_norm(const NumericVector< Number > &v, unsigned int var, FEMNormType norm_type) const
Definition: system.C:1345
Real libMesh::System::calculate_norm ( const NumericVector< Number > &  v,
const SystemNorm norm,
std::set< unsigned int > *  skip_dimensions = libmesh_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 1386 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_nullptr, 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().

1389 {
1390  // This function must be run on all processors at once
1391  parallel_object_only();
1392 
1393  LOG_SCOPE ("calculate_norm()", "System");
1394 
1395  // Zero the norm before summation
1396  Real v_norm = 0.;
1397 
1398  if (norm.is_discrete())
1399  {
1400  //Check to see if all weights are 1.0 and all types are equal
1401  FEMNormType norm_type0 = norm.type(0);
1402  unsigned int check_var = 0;
1403  for (; check_var != this->n_vars(); ++check_var)
1404  if ((norm.weight(check_var) != 1.0) || (norm.type(check_var) != norm_type0))
1405  break;
1406 
1407  //All weights were 1.0 so just do the full vector discrete norm
1408  if (check_var == this->n_vars())
1409  {
1410  if (norm_type0 == DISCRETE_L1)
1411  return v.l1_norm();
1412  if (norm_type0 == DISCRETE_L2)
1413  return v.l2_norm();
1414  if (norm_type0 == DISCRETE_L_INF)
1415  return v.linfty_norm();
1416  else
1417  libmesh_error_msg("Invalid norm_type0 = " << norm_type0);
1418  }
1419 
1420  for (unsigned int var=0; var != this->n_vars(); ++var)
1421  {
1422  // Skip any variables we don't need to integrate
1423  if (norm.weight(var) == 0.0)
1424  continue;
1425 
1426  v_norm += norm.weight(var) * discrete_var_norm(v, var, norm.type(var));
1427  }
1428 
1429  return v_norm;
1430  }
1431 
1432  // Localize the potentially parallel vector
1433  std::unique_ptr<NumericVector<Number>> local_v = NumericVector<Number>::build(this->comm());
1434  local_v->init(v.size(), true, SERIAL);
1435  v.localize (*local_v, _dof_map->get_send_list());
1436 
1437  // I'm not sure how best to mix Hilbert norms on some variables (for
1438  // which we'll want to square then sum then square root) with norms
1439  // like L_inf (for which we'll just want to take an absolute value
1440  // and then sum).
1441  bool using_hilbert_norm = true,
1442  using_nonhilbert_norm = true;
1443 
1444  // Loop over all variables
1445  for (unsigned int var=0; var != this->n_vars(); ++var)
1446  {
1447  // Skip any variables we don't need to integrate
1448  Real norm_weight_sq = norm.weight_sq(var);
1449  if (norm_weight_sq == 0.0)
1450  continue;
1451  Real norm_weight = norm.weight(var);
1452 
1453  // Check for unimplemented norms (rather than just returning 0).
1454  FEMNormType norm_type = norm.type(var);
1455  if ((norm_type==H1) ||
1456  (norm_type==H2) ||
1457  (norm_type==L2) ||
1458  (norm_type==H1_SEMINORM) ||
1459  (norm_type==H2_SEMINORM))
1460  {
1461  if (!using_hilbert_norm)
1462  libmesh_not_implemented();
1463  using_nonhilbert_norm = false;
1464  }
1465  else if ((norm_type==L1) ||
1466  (norm_type==L_INF) ||
1467  (norm_type==W1_INF_SEMINORM) ||
1468  (norm_type==W2_INF_SEMINORM))
1469  {
1470  if (!using_nonhilbert_norm)
1471  libmesh_not_implemented();
1472  using_hilbert_norm = false;
1473  }
1474  else
1475  libmesh_not_implemented();
1476 
1477  const FEType & fe_type = this->get_dof_map().variable_type(var);
1478 
1479  // Allow space for dims 0-3, even if we don't use them all
1480  std::vector<std::unique_ptr<FEBase>> fe_ptrs(4);
1481  std::vector<std::unique_ptr<QBase>> q_rules(4);
1482 
1483  const std::set<unsigned char> & elem_dims = _mesh.elem_dimensions();
1484 
1485  // Prepare finite elements for each dimension present in the mesh
1486  for (std::set<unsigned char>::const_iterator d_it = elem_dims.begin();
1487  d_it != elem_dims.end(); ++d_it)
1488  {
1489  if (skip_dimensions && skip_dimensions->find(*d_it) != skip_dimensions->end())
1490  continue;
1491 
1492  // Construct quadrature and finite element objects
1493  q_rules[*d_it] = fe_type.default_quadrature_rule (*d_it);
1494  fe_ptrs[*d_it] = FEBase::build(*d_it, fe_type);
1495 
1496  // Attach quadrature rule to FE object
1497  fe_ptrs[*d_it]->attach_quadrature_rule (q_rules[*d_it].get());
1498  }
1499 
1500  std::vector<dof_id_type> dof_indices;
1501 
1502  // Begin the loop over the elements
1503  for (const auto & elem : this->get_mesh().active_local_element_ptr_range())
1504  {
1505  const unsigned int dim = elem->dim();
1506 
1507 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1508 
1509  // One way for implementing this would be to exchange the fe with the FEInterface- class.
1510  // However, it needs to be discussed whether integral-norms make sense for infinite elements.
1511  // or in which sense they could make sense.
1512  if (elem->infinite() )
1513  libmesh_not_implemented();
1514 
1515 #endif
1516 
1517  if (skip_dimensions && skip_dimensions->find(dim) != skip_dimensions->end())
1518  continue;
1519 
1520  FEBase * fe = fe_ptrs[dim].get();
1521  QBase * qrule = q_rules[dim].get();
1522  libmesh_assert(fe);
1523  libmesh_assert(qrule);
1524 
1525  const std::vector<Real> & JxW = fe->get_JxW();
1526  const std::vector<std::vector<Real>> * phi = libmesh_nullptr;
1527  if (norm_type == H1 ||
1528  norm_type == H2 ||
1529  norm_type == L2 ||
1530  norm_type == L1 ||
1531  norm_type == L_INF)
1532  phi = &(fe->get_phi());
1533 
1534  const std::vector<std::vector<RealGradient>> * dphi = libmesh_nullptr;
1535  if (norm_type == H1 ||
1536  norm_type == H2 ||
1537  norm_type == H1_SEMINORM ||
1538  norm_type == W1_INF_SEMINORM)
1539  dphi = &(fe->get_dphi());
1540 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1541  const std::vector<std::vector<RealTensor>> * d2phi = libmesh_nullptr;
1542  if (norm_type == H2 ||
1543  norm_type == H2_SEMINORM ||
1544  norm_type == W2_INF_SEMINORM)
1545  d2phi = &(fe->get_d2phi());
1546 #endif
1547 
1548  fe->reinit (elem);
1549 
1550  this->get_dof_map().dof_indices (elem, dof_indices, var);
1551 
1552  const unsigned int n_qp = qrule->n_points();
1553 
1554  const unsigned int n_sf = cast_int<unsigned int>
1555  (dof_indices.size());
1556 
1557  // Begin the loop over the Quadrature points.
1558  for (unsigned int qp=0; qp<n_qp; qp++)
1559  {
1560  if (norm_type == L1)
1561  {
1562  Number u_h = 0.;
1563  for (unsigned int i=0; i != n_sf; ++i)
1564  u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1565  v_norm += norm_weight *
1566  JxW[qp] * std::abs(u_h);
1567  }
1568 
1569  if (norm_type == L_INF)
1570  {
1571  Number u_h = 0.;
1572  for (unsigned int i=0; i != n_sf; ++i)
1573  u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1574  v_norm = std::max(v_norm, norm_weight * std::abs(u_h));
1575  }
1576 
1577  if (norm_type == H1 ||
1578  norm_type == H2 ||
1579  norm_type == L2)
1580  {
1581  Number u_h = 0.;
1582  for (unsigned int i=0; i != n_sf; ++i)
1583  u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1584  v_norm += norm_weight_sq *
1585  JxW[qp] * TensorTools::norm_sq(u_h);
1586  }
1587 
1588  if (norm_type == H1 ||
1589  norm_type == H2 ||
1590  norm_type == H1_SEMINORM)
1591  {
1592  Gradient grad_u_h;
1593  for (unsigned int i=0; i != n_sf; ++i)
1594  grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
1595  v_norm += norm_weight_sq *
1596  JxW[qp] * grad_u_h.norm_sq();
1597  }
1598 
1599  if (norm_type == W1_INF_SEMINORM)
1600  {
1601  Gradient grad_u_h;
1602  for (unsigned int i=0; i != n_sf; ++i)
1603  grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
1604  v_norm = std::max(v_norm, norm_weight * grad_u_h.norm());
1605  }
1606 
1607 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1608  if (norm_type == H2 ||
1609  norm_type == H2_SEMINORM)
1610  {
1611  Tensor hess_u_h;
1612  for (unsigned int i=0; i != n_sf; ++i)
1613  hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
1614  v_norm += norm_weight_sq *
1615  JxW[qp] * hess_u_h.norm_sq();
1616  }
1617 
1618  if (norm_type == W2_INF_SEMINORM)
1619  {
1620  Tensor hess_u_h;
1621  for (unsigned int i=0; i != n_sf; ++i)
1622  hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
1623  v_norm = std::max(v_norm, norm_weight * hess_u_h.norm());
1624  }
1625 #endif
1626  }
1627  }
1628  }
1629 
1630  if (using_hilbert_norm)
1631  {
1632  this->comm().sum(v_norm);
1633  v_norm = std::sqrt(v_norm);
1634  }
1635  else
1636  {
1637  this->comm().max(v_norm);
1638  }
1639 
1640  return v_norm;
1641 }
double abs(double a)
virtual numeric_index_type size() const =0
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1719
void add_scaled(const TypeTensor< T2 > &, const T)
Definition: type_tensor.h:777
void add_scaled(const TypeVector< T2 > &, const T)
Definition: type_vector.h:624
virtual Real l2_norm() const =0
const class libmesh_nullptr_t libmesh_nullptr
std::unique_ptr< DofMap > _dof_map
Definition: system.h:1865
virtual Real l1_norm() const =0
long double max(long double a, double b)
const MeshBase & get_mesh() const
Definition: system.h:2014
const DofMap & get_dof_map() const
Definition: system.h:2030
NumberVectorValue Gradient
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
const std::set< unsigned char > & elem_dimensions() const
Definition: mesh_base.h:207
FEGenericBase< Real > FEBase
virtual Real linfty_norm() const =0
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
NumberTensorValue Tensor
const Parallel::Communicator & comm() const
ParallelType type() const
unsigned int n_vars() const
Definition: system.h:2086
MeshBase & _mesh
Definition: system.h:1877
virtual void localize(std::vector< T > &v_local) const =0
Real discrete_var_norm(const NumericVector< Number > &v, unsigned int var, FEMNormType norm_type) const
Definition: system.C:1345
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:1924
void libMesh::ContinuationSystem::clear ( )
virtual

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().

81 {
82  // FIXME: Do anything here, e.g. zero vectors, etc?
83 
84  // Call the Parent's clear function
85  Parent::clear();
86 }
virtual void clear() libmesh_override
Definition: diff_system.C:69
virtual void libMesh::DifferentiablePhysics::clear_physics ( )
virtualinherited

Clear any data structures associated with the physics.

Referenced by libMesh::DifferentiableSystem::clear(), and libMesh::DifferentiablePhysics::DifferentiablePhysics().

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 {}
virtual std::unique_ptr<DifferentiableQoI> libMesh::DifferentiableSystem::clone ( )
inlinevirtualinherited

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().

170  {
171  libmesh_not_implemented();
172  // dummy
173  return std::unique_ptr<DifferentiableQoI>(this);
174  }
virtual std::unique_ptr<DifferentiablePhysics> libMesh::DifferentiableSystem::clone_physics ( )
inlinevirtualinherited

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

Implements libMesh::DifferentiablePhysics.

Definition at line 159 of file diff_system.h.

160  {
161  libmesh_not_implemented();
162  // dummy
163  return std::unique_ptr<DifferentiablePhysics>(this);
164  }
const Parallel::Communicator& libMesh::ParallelObject::comm ( ) const
inlineinherited
Returns
A reference to the Parallel::Communicator object used by this mesh.

Definition at line 87 of file parallel_object.h.

References libMesh::ParallelObject::_communicator.

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), libMesh::__libmesh_petsc_diff_solver_monitor(), libMesh::__libmesh_petsc_diff_solver_residual(), libMesh::__libmesh_petsc_snes_fd_residual(), libMesh::__libmesh_petsc_snes_jacobian(), libMesh::__libmesh_petsc_snes_mffd_residual(), libMesh::__libmesh_petsc_snes_postcheck(), libMesh::__libmesh_petsc_snes_residual(), libMesh::__libmesh_tao_equality_constraints(), libMesh::__libmesh_tao_equality_constraints_jacobian(), libMesh::__libmesh_tao_gradient(), libMesh::__libmesh_tao_hessian(), libMesh::__libmesh_tao_inequality_constraints(), libMesh::__libmesh_tao_inequality_constraints_jacobian(), libMesh::__libmesh_tao_objective(), libMesh::MeshRefinement::_coarsen_elements(), libMesh::ExactSolution::_compute_error(), libMesh::ParmetisPartitioner::_do_repartition(), libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::BoundaryInfo::_find_id_maps(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_mult(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult_add(), libMesh::EquationSystems::_read_impl(), libMesh::MeshRefinement::_refine_elements(), libMesh::MeshRefinement::_smooth_flags(), libMesh::ImplicitSystem::add_matrix(), libMesh::System::add_vector(), libMesh::EigenSparseLinearSolver< T >::adjoint_solve(), libMesh::UnstructuredMesh::all_second_order(), libMesh::MeshTools::Modification::all_tri(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::FEMSystem::assemble_qoi(), libMesh::MeshCommunication::assign_global_indices(), libMesh::ParmetisPartitioner::assign_partitioning(), libMesh::DofMap::attach_matrix(), libMesh::Parallel::BinSorter< KeyType, IdxType >::binsort(), libMesh::Parallel::Sort< KeyType, IdxType >::binsort(), libMesh::MeshCommunication::broadcast(), libMesh::SparseMatrix< T >::build(), libMesh::MeshTools::Generation::build_extrusion(), libMesh::Parallel::Histogram< KeyType, IdxType >::build_histogram(), libMesh::PetscNonlinearSolver< T >::build_mat_null_space(), libMesh::BoundaryInfo::build_node_list_from_side_list(), libMesh::EquationSystems::build_parallel_solution_vector(), libMesh::MeshBase::cache_elem_dims(), libMesh::System::calculate_norm(), libMesh::DofMap::check_dirichlet_bcid_consistency(), libMesh::DistributedVector< T >::clone(), libMesh::EigenSparseVector< T >::clone(), libMesh::LaspackVector< T >::clone(), libMesh::EpetraVector< T >::clone(), libMesh::PetscVector< T >::clone(), libMesh::EpetraVector< T >::close(), libMesh::Parallel::Sort< KeyType, IdxType >::communicate_bins(), libMesh::Nemesis_IO_Helper::compute_num_global_elem_blocks(), libMesh::Nemesis_IO_Helper::compute_num_global_nodesets(), libMesh::Nemesis_IO_Helper::compute_num_global_sidesets(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), libMesh::ExodusII_IO::copy_elemental_solution(), libMesh::MeshTools::correct_node_proc_ids(), libMesh::MeshTools::create_bounding_box(), libMesh::MeshTools::create_nodal_bounding_box(), libMesh::MeshRefinement::create_parent_error_vector(), libMesh::MeshTools::create_processor_bounding_box(), libMesh::MeshTools::create_subdomain_bounding_box(), libMesh::MeshCommunication::delete_remote_elements(), libMesh::DofMap::distribute_dofs(), DMlibMeshFunction(), DMlibMeshJacobian(), DMlibMeshSetSystem_libMesh(), DMVariableBounds_libMesh(), libMesh::MeshRefinement::eliminate_unrefined_patches(), libMesh::EpetraVector< T >::EpetraVector(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::MeshRefinement::flag_elements_by_elem_fraction(), libMesh::MeshRefinement::flag_elements_by_error_fraction(), libMesh::MeshRefinement::flag_elements_by_nelem_target(), libMesh::MeshCommunication::gather(), libMesh::MeshCommunication::gather_neighboring_elements(), libMesh::CondensedEigenSystem::get_eigenpair(), libMesh::DofMap::get_info(), libMesh::ImplicitSystem::get_linear_solver(), libMesh::EquationSystems::get_solution(), libMesh::LocationMap< T >::init(), libMesh::TopologyMap::init(), libMesh::TimeSolver::init(), libMesh::TaoOptimizationSolver< T >::init(), libMesh::PetscNonlinearSolver< T >::init(), libMesh::DistributedVector< T >::init(), libMesh::EpetraVector< T >::init(), libMesh::PetscVector< T >::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::EigenSystem::init_data(), libMesh::EigenSystem::init_matrices(), libMesh::ParmetisPartitioner::initialize(), libMesh::OptimizationSystem::initialize_equality_constraints_storage(), libMesh::OptimizationSystem::initialize_inequality_constraints_storage(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Elem >(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_topology_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_boundary_ids(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), libMesh::MeshTools::libmesh_assert_valid_neighbors(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_flags(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_object_ids(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_p_levels(), libMesh::MeshTools::libmesh_assert_valid_refinement_flags(), libMesh::MeshTools::libmesh_assert_valid_unique_ids(), libMesh::libmesh_petsc_snes_residual_helper(), libMesh::MeshRefinement::limit_level_mismatch_at_edge(), libMesh::MeshRefinement::limit_level_mismatch_at_node(), libMesh::MeshRefinement::limit_overrefined_boundary(), libMesh::MeshRefinement::limit_underrefined_boundary(), libMesh::MeshRefinement::make_coarsening_compatible(), libMesh::MeshCommunication::make_elems_parallel_consistent(), libMesh::MeshRefinement::make_flags_parallel_consistent(), libMesh::MeshCommunication::make_new_node_proc_ids_parallel_consistent(), libMesh::MeshCommunication::make_new_nodes_parallel_consistent(), libMesh::MeshCommunication::make_node_ids_parallel_consistent(), libMesh::MeshCommunication::make_node_proc_ids_parallel_consistent(), libMesh::MeshCommunication::make_node_unique_ids_parallel_consistent(), libMesh::MeshCommunication::make_nodes_parallel_consistent(), libMesh::MeshCommunication::make_p_levels_parallel_consistent(), libMesh::MeshRefinement::make_refinement_compatible(), libMesh::DistributedVector< T >::max(), libMesh::FEMSystem::mesh_position_set(), libMesh::MeshSerializer::MeshSerializer(), libMesh::DistributedVector< T >::min(), libMesh::DistributedMesh::n_active_elem(), libMesh::MeshTools::n_active_levels(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::BoundaryInfo::n_edge_conds(), libMesh::CondensedEigenSystem::n_global_non_condensed_dofs(), libMesh::MeshTools::n_levels(), libMesh::BoundaryInfo::n_nodeset_conds(), libMesh::MeshTools::n_p_levels(), libMesh::BoundaryInfo::n_shellface_conds(), libMesh::DistributedMesh::parallel_max_elem_id(), libMesh::DistributedMesh::parallel_max_node_id(), libMesh::ReplicatedMesh::parallel_max_unique_id(), libMesh::DistributedMesh::parallel_max_unique_id(), libMesh::DistributedMesh::parallel_n_elem(), libMesh::DistributedMesh::parallel_n_nodes(), libMesh::SparsityPattern::Build::parallel_sync(), libMesh::MeshTools::paranoid_n_levels(), libMesh::Partitioner::partition(), libMesh::LinearPartitioner::partition_range(), libMesh::MetisPartitioner::partition_range(), libMesh::Partitioner::partition_unpartitioned_elements(), libMesh::petsc_auto_fieldsplit(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::MeshBase::prepare_for_use(), libMesh::SparseMatrix< T >::print(), libMesh::Nemesis_IO::read(), libMesh::XdrIO::read(), libMesh::CheckpointIO::read_header(), libMesh::XdrIO::read_header(), libMesh::XdrIO::read_serialized_bc_names(), libMesh::XdrIO::read_serialized_bcs_helper(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::MeshBase::recalculate_n_partitions(), libMesh::MeshCommunication::redistribute(), libMesh::MeshRefinement::refine_and_coarsen_elements(), libMesh::DistributedMesh::renumber_dof_objects(), libMesh::CheckpointIO::select_split_config(), libMesh::MeshCommunication::send_coarse_ghosts(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::Partitioner::set_parent_processor_ids(), libMesh::PetscDiffSolver::setup_petsc_data(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::Parallel::Sort< KeyType, IdxType >::sort(), libMesh::split_mesh(), libMesh::MeshBase::subdomain_ids(), libMesh::BoundaryInfo::sync(), libMesh::Parallel::sync_element_data_by_parent_id(), libMesh::Parallel::sync_node_data_by_element_id(), libMesh::MeshRefinement::test_level_one(), libMesh::MeshRefinement::test_unflagged(), libMesh::MeshTools::total_weight(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::NameBasedIO::write(), libMesh::XdrIO::write(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), libMesh::DistributedVector< T >::zero_clone(), libMesh::EigenSparseVector< T >::zero_clone(), libMesh::LaspackVector< T >::zero_clone(), libMesh::EpetraVector< T >::zero_clone(), and libMesh::PetscVector< T >::zero_clone().

88  { return _communicator; }
const Parallel::Communicator & _communicator
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 512 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(), and libMesh::System::set_adjoint_already_solved().

515 {
516  // we do not care for matrices, but for vectors
517  libmesh_assert (_is_initialized);
518  libmesh_assert (other_system._is_initialized);
519 
520  if (verbose)
521  {
522  libMesh::out << " Systems \"" << _sys_name << "\"" << std::endl;
523  libMesh::out << " comparing matrices not supported." << std::endl;
524  libMesh::out << " comparing names...";
525  }
526 
527  // compare the name: 0 means identical
528  const int name_result = _sys_name.compare(other_system.name());
529  if (verbose)
530  {
531  if (name_result == 0)
532  libMesh::out << " identical." << std::endl;
533  else
534  libMesh::out << " names not identical." << std::endl;
535  libMesh::out << " comparing solution vector...";
536  }
537 
538 
539  // compare the solution: -1 means identical
540  const int solu_result = solution->compare (*other_system.solution.get(),
541  threshold);
542 
543  if (verbose)
544  {
545  if (solu_result == -1)
546  libMesh::out << " identical up to threshold." << std::endl;
547  else
548  libMesh::out << " first difference occurred at index = "
549  << solu_result << "." << std::endl;
550  }
551 
552 
553  // safety check, whether we handle at least the same number
554  // of vectors
555  std::vector<int> ov_result;
556 
557  if (this->n_vectors() != other_system.n_vectors())
558  {
559  if (verbose)
560  {
561  libMesh::out << " Fatal difference. This system handles "
562  << this->n_vectors() << " add'l vectors," << std::endl
563  << " while the other system handles "
564  << other_system.n_vectors()
565  << " add'l vectors." << std::endl
566  << " Aborting comparison." << std::endl;
567  }
568  return false;
569  }
570  else if (this->n_vectors() == 0)
571  {
572  // there are no additional vectors...
573  ov_result.clear ();
574  }
575  else
576  {
577  // compare other vectors
578  for (const_vectors_iterator pos = _vectors.begin();
579  pos != _vectors.end(); ++pos)
580  {
581  if (verbose)
582  libMesh::out << " comparing vector \""
583  << pos->first << "\" ...";
584 
585  // assume they have the same name
586  const NumericVector<Number> & other_system_vector =
587  other_system.get_vector(pos->first);
588 
589  ov_result.push_back(pos->second->compare (other_system_vector,
590  threshold));
591 
592  if (verbose)
593  {
594  if (ov_result[ov_result.size()-1] == -1)
595  libMesh::out << " identical up to threshold." << std::endl;
596  else
597  libMesh::out << " first difference occurred at" << std::endl
598  << " index = " << ov_result[ov_result.size()-1] << "." << std::endl;
599  }
600 
601  }
602 
603  } // finished comparing additional vectors
604 
605 
606  bool overall_result;
607 
608  // sum up the results
609  if ((name_result==0) && (solu_result==-1))
610  {
611  if (ov_result.size()==0)
612  overall_result = true;
613  else
614  {
615  bool ov_identical;
616  unsigned int n = 0;
617  do
618  {
619  ov_identical = (ov_result[n]==-1);
620  n++;
621  }
622  while (ov_identical && n<ov_result.size());
623  overall_result = ov_identical;
624  }
625  }
626  else
627  overall_result = false;
628 
629  if (verbose)
630  {
631  libMesh::out << " finished comparisons, ";
632  if (overall_result)
633  libMesh::out << "found no differences." << std::endl << std::endl;
634  else
635  libMesh::out << "found differences." << std::endl << std::endl;
636  }
637 
638  return overall_result;
639 }
bool _is_initialized
Definition: system.h:1952
std::map< std::string, NumericVector< Number > * >::const_iterator const_vectors_iterator
Definition: system.h:735
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1509
std::map< std::string, NumericVector< Number > * > _vectors
Definition: system.h:1916
const std::string _sys_name
Definition: system.h:1882
OStreamProxy out(std::cout)
unsigned int n_vectors() const
Definition: system.h:2214
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().

359 {
360  // Be sure the user has set the continuation parameter pointer
362  libmesh_error_msg("You must set the continuation_parameter pointer " \
363  << "to a member variable of the derived class, preferably in the " \
364  << "Derived class's init_data function. This is how the ContinuationSystem " \
365  << "updates the continuation parameter.");
366 
367  // Use extra precision for all the numbers printed in this function.
368  std::streamsize old_precision = libMesh::out.precision();
370  libMesh::out.setf(std::ios_base::scientific);
371 
372  // We can't start solving the augmented PDE system unless the tangent
373  // vectors have been initialized. This only needs to occur once.
374  if (!tangent_initialized)
376 
377  // Save the old value of -du/dlambda. This will be used after the Newton iterations
378  // to compute the angle between previous tangent vectors. This cosine of this angle is
379  //
380  // tau := abs( (du/d(lambda)_i , du/d(lambda)_{i-1}) / (||du/d(lambda)_i|| * ||du/d(lambda)_{i-1}||) )
381  //
382  // The scaling factor tau (which should vary between 0 and 1) is used to shrink the step-size ds
383  // when we are approaching a turning point. Note that it can only shrink the step size.
384  *y_old = *y;
385 
386  // Set pointer to underlying Newton solver
387  if (!newton_solver)
388  newton_solver = cast_ptr<NewtonSolver *> (this->time_solver->diff_solver().get());
389 
390  // A pair for catching return values from linear system solves.
391  std::pair<unsigned int, Real> rval;
392 
393  // Convergence flag for the entire arcstep
394  bool arcstep_converged = false;
395 
396  // Begin loop over arcstep reductions.
397  for (unsigned int ns=0; ns<n_arclength_reductions; ++ns)
398  {
399  if (!quiet)
400  {
401  libMesh::out << "Current arclength stepsize, ds_current=" << ds_current << std::endl;
402  libMesh::out << "Current parameter value, lambda=" << *continuation_parameter << std::endl;
403  }
404 
405  // Upon exit from the nonlinear loop, the newton_converged flag
406  // will tell us the convergence status of Newton's method.
407  bool newton_converged = false;
408 
409  // The nonlinear residual before *any* nonlinear steps have been taken.
410  Real nonlinear_residual_firststep = 0.;
411 
412  // The nonlinear residual from the current "k" Newton step, before the Newton step
413  Real nonlinear_residual_beforestep = 0.;
414 
415  // The nonlinear residual from the current "k" Newton step, after the Newton step
416  Real nonlinear_residual_afterstep = 0.;
417 
418  // The linear solver tolerance, can be updated dynamically at each Newton step.
419  Real current_linear_tolerance = 0.;
420 
421  // The nonlinear loop
423  {
424  libMesh::out << "\n === Starting Newton step " << newton_step << " ===" << std::endl;
425 
426  // Set the linear system solver tolerance
427  // // 1.) Set the current linear tolerance based as a multiple of the current residual of the system.
428  // const Real residual_multiple = 1.e-4;
429  // Real current_linear_tolerance = residual_multiple*nonlinear_residual_beforestep;
430 
431  // // But if the current residual isn't small, don't let the solver exit with zero iterations!
432  // if (current_linear_tolerance > 1.)
433  // current_linear_tolerance = residual_multiple;
434 
435  // 2.) Set the current linear tolerance based on the method based on technique of Eisenstat & Walker.
436  if (newton_step==0)
437  {
438  // At first step, only try reducing the residual by a small amount
439  current_linear_tolerance = initial_newton_tolerance;//0.01;
440  }
441 
442  else
443  {
444  // The new tolerance is based on the ratio of the most recent tolerances
445  const Real alp=0.5*(1.+std::sqrt(5.));
446  const Real gam=0.9;
447 
448  libmesh_assert_not_equal_to (nonlinear_residual_beforestep, 0.0);
449  libmesh_assert_not_equal_to (nonlinear_residual_afterstep, 0.0);
450 
451  current_linear_tolerance = std::min(gam*std::pow(nonlinear_residual_afterstep/nonlinear_residual_beforestep, alp),
452  current_linear_tolerance*current_linear_tolerance
453  );
454 
455  // Don't let it get ridiculously small!!
456  if (current_linear_tolerance < 1.e-12)
457  current_linear_tolerance = 1.e-12;
458  }
459 
460  if (!quiet)
461  libMesh::out << "Using current_linear_tolerance=" << current_linear_tolerance << std::endl;
462 
463 
464  // Assemble the residual (and Jacobian).
465  rhs_mode = Residual;
466  assembly(true, // Residual
467  true); // Jacobian
468  rhs->close();
469 
470  // Save the current nonlinear residual. We don't need to recompute the residual unless
471  // this is the first step, since it was already computed as part of the convergence check
472  // at the end of the last loop iteration.
473  if (newton_step==0)
474  {
475  nonlinear_residual_beforestep = rhs->l2_norm();
476 
477  // Store the residual before any steps have been taken. This will *not*
478  // be updated at each step, and can be used to see if any progress has
479  // been made from the initial residual at later steps.
480  nonlinear_residual_firststep = nonlinear_residual_beforestep;
481 
482  const Real old_norm_u = solution->l2_norm();
483  libMesh::out << " (before step) ||R||_{L2} = " << nonlinear_residual_beforestep << std::endl;
484  libMesh::out << " (before step) ||R||_{L2}/||u|| = " << nonlinear_residual_beforestep / old_norm_u << std::endl;
485 
486  // In rare cases (very small arcsteps), it's possible that the residual is
487  // already below our absolute linear tolerance.
488  if (nonlinear_residual_beforestep < solution_tolerance)
489  {
490  if (!quiet)
491  libMesh::out << "Initial guess satisfied linear tolerance, exiting with zero Newton iterations!" << std::endl;
492 
493  // Since we go straight from here to the solve of the next tangent, we
494  // have to close the matrix before it can be assembled again.
495  matrix->close();
496  newton_converged=true;
497  break; // out of Newton iterations, with newton_converged=true
498  }
499  }
500 
501  else
502  {
503  nonlinear_residual_beforestep = nonlinear_residual_afterstep;
504  }
505 
506 
507  // Solve the linear system G_u*z = G
508  // Initial guess?
509  z->zero(); // It seems to be extremely important to zero z here, otherwise the solver quits early.
510  z->close();
511 
512  // It's possible that we have selected the current_linear_tolerance so large that
513  // a guess of z=zero yields a linear system residual |Az + R| small enough that the
514  // linear solver exits in zero iterations. If this happens, we will reduce the
515  // current_linear_tolerance until the linear solver does at least 1 iteration.
516  do
517  {
518  rval =
519  linear_solver->solve(*matrix,
520  *z,
521  *rhs,
522  //1.e-12,
523  current_linear_tolerance,
524  newton_solver->max_linear_iterations); // max linear iterations
525 
526  if (rval.first==0)
527  {
528  if (newton_step==0)
529  {
530  libMesh::out << "Repeating initial solve with smaller linear tolerance!" << std::endl;
531  current_linear_tolerance *= initial_newton_tolerance; // reduce the linear tolerance to force the solver to do some work
532  }
533  else
534  // We shouldn't get here ... it means the linear solver did no work on a Newton
535  // step other than the first one. If this happens, we need to think more about our
536  // tolerance selection.
537  libmesh_error_msg("Linear solver did no work!");
538  }
539 
540  } while (rval.first==0);
541 
542 
543  if (!quiet)
544  libMesh::out << " G_u*z = G solver converged at step "
545  << rval.first
546  << " linear tolerance = "
547  << rval.second
548  << "."
549  << std::endl;
550 
551  // Sometimes (I am not sure why) the linear solver exits after zero iterations.
552  // Perhaps it is hitting PETSc's divergence tolerance dtol??? If this occurs,
553  // we should break out of the Newton iteration loop because nothing further is
554  // going to happen... Of course if the tolerance is already small enough after
555  // zero iterations (how can this happen?!) we should not quit.
556  if ((rval.first == 0) && (rval.second > current_linear_tolerance*nonlinear_residual_beforestep))
557  {
558  if (!quiet)
559  libMesh::out << "Linear solver exited in zero iterations!" << std::endl;
560 
561  // Try to find out the reason for convergence/divergence
562  linear_solver->print_converged_reason();
563 
564  break; // out of Newton iterations
565  }
566 
567  // Note: need to scale z by -1 since our code always solves Jx=R
568  // instead of Jx=-R.
569  z->scale(-1.);
570  z->close();
571 
572 
573 
574 
575 
576 
577  // Assemble the G_Lambda vector, skip residual.
578  rhs_mode = G_Lambda;
579 
580  // Assemble both rhs and Jacobian
581  assembly(true, // Residual
582  false); // Jacobian
583 
584  // Not sure if this is really necessary
585  rhs->close();
586  const Real yrhsnorm=rhs->l2_norm();
587  if (yrhsnorm == 0.0)
588  libmesh_error_msg("||G_Lambda|| = 0");
589 
590  // We select a tolerance for the y-system which is based on the inexact Newton
591  // tolerance but scaled by an extra term proportional to the RHS (which is not -> 0 in this case)
592  const Real ysystemtol=current_linear_tolerance*(nonlinear_residual_beforestep/yrhsnorm);
593  if (!quiet)
594  libMesh::out << "ysystemtol=" << ysystemtol << std::endl;
595 
596  // Solve G_u*y = G_{\lambda}
597  // FIXME: Initial guess? This is really a solve for -du/dlambda so we could try
598  // initializing it with the latest approximation to that... du/dlambda ~ du/ds * ds/dlambda
599  //*y = *solution;
600  //y->add(-1., *previous_u);
601  //y->scale(-1. / (*continuation_parameter - old_continuation_parameter)); // Be careful of divide by zero...
602  //y->close();
603 
604  // const unsigned int max_attempts=1;
605  // unsigned int attempt=0;
606  // do
607  // {
608  // if (!quiet)
609  // libMesh::out << "Trying to solve tangent system, attempt " << attempt << std::endl;
610 
611  rval =
612  linear_solver->solve(*matrix,
613  *y,
614  *rhs,
615  //1.e-12,
616  ysystemtol,
617  newton_solver->max_linear_iterations); // max linear iterations
618 
619  if (!quiet)
620  libMesh::out << " G_u*y = G_{lambda} solver converged at step "
621  << rval.first
622  << ", linear tolerance = "
623  << rval.second
624  << "."
625  << std::endl;
626 
627  // Sometimes (I am not sure why) the linear solver exits after zero iterations.
628  // Perhaps it is hitting PETSc's divergence tolerance dtol??? If this occurs,
629  // we should break out of the Newton iteration loop because nothing further is
630  // going to happen...
631  if ((rval.first == 0) && (rval.second > ysystemtol))
632  {
633  if (!quiet)
634  libMesh::out << "Linear solver exited in zero iterations!" << std::endl;
635 
636  break; // out of Newton iterations
637  }
638 
639  // ++attempt;
640  // } while ((attempt<max_attempts) && (rval.first==newton_solver->max_linear_iterations));
641 
642 
643 
644 
645 
646  // Compute N, the residual of the arclength constraint eqn.
647  // Note 1: N(u,lambda,s) := (u-u_{old}, du_ds) + (lambda-lambda_{old}, dlambda_ds) - _ds
648  // We temporarily use the delta_u vector as a temporary vector for this calculation.
649  *delta_u = *solution;
650  delta_u->add(-1., *previous_u);
651 
652  // First part of the arclength constraint
654  const Number N2 = ((*continuation_parameter) - old_continuation_parameter)*dlambda_ds;
655  const Number N3 = ds_current;
656 
657  if (!quiet)
658  {
659  libMesh::out << " N1=" << N1 << std::endl;
660  libMesh::out << " N2=" << N2 << std::endl;
661  libMesh::out << " N3=" << N3 << std::endl;
662  }
663 
664  // The arclength constraint value
665  const Number N = N1+N2-N3;
666 
667  if (!quiet)
668  libMesh::out << " N=" << N << std::endl;
669 
670  const Number duds_dot_z = du_ds->dot(*z);
671  const Number duds_dot_y = du_ds->dot(*y);
672 
673  //libMesh::out << "duds_dot_z=" << duds_dot_z << std::endl;
674  //libMesh::out << "duds_dot_y=" << duds_dot_y << std::endl;
675  //libMesh::out << "dlambda_ds=" << dlambda_ds << std::endl;
676 
677  const Number delta_lambda_numerator = -(N + Theta_LOCA*Theta_LOCA*Theta*duds_dot_z);
678  const Number delta_lambda_denominator = (dlambda_ds - Theta_LOCA*Theta_LOCA*Theta*duds_dot_y);
679 
680  libmesh_assert_not_equal_to (delta_lambda_denominator, 0.0);
681 
682  // Now, we are ready to compute the step delta_lambda
683  const Number delta_lambda_comp = delta_lambda_numerator /
684  delta_lambda_denominator;
685  // Lambda is real-valued
686  const Real delta_lambda = libmesh_real(delta_lambda_comp);
687 
688  // Knowing delta_lambda, we are ready to update delta_u
689  // delta_u = z - delta_lambda*y
690  delta_u->zero();
691  delta_u->add(1., *z);
692  delta_u->add(-delta_lambda, *y);
693  delta_u->close();
694 
695  // Update the system solution and the continuation parameter.
696  solution->add(1., *delta_u);
697  solution->close();
698  *continuation_parameter += delta_lambda;
699 
700  // Did the Newton step actually reduce the residual?
701  rhs_mode = Residual;
702  assembly(true, // Residual
703  false); // Jacobian
704  rhs->close();
705  nonlinear_residual_afterstep = rhs->l2_norm();
706 
707 
708  // In a "normal" Newton step, ||du||/||R|| > 1 since the most recent
709  // step is where you "just were" and the current residual is where
710  // you are now. It can occur that ||du||/||R|| < 1, but these are
711  // likely not good cases to attempt backtracking (?).
712  const Real norm_du_norm_R = delta_u->l2_norm() / nonlinear_residual_afterstep;
713  if (!quiet)
714  libMesh::out << " norm_du_norm_R=" << norm_du_norm_R << std::endl;
715 
716 
717  // Factor to decrease the stepsize by for backtracking
718  Real newton_stepfactor = 1.;
719 
720  const bool attempt_backtracking =
721  (nonlinear_residual_afterstep > solution_tolerance)
722  && (nonlinear_residual_afterstep > nonlinear_residual_beforestep)
723  && (n_backtrack_steps>0)
724  && (norm_du_norm_R > 1.)
725  ;
726 
727  // If residual is not reduced, do Newton back tracking.
728  if (attempt_backtracking)
729  {
730  if (!quiet)
731  libMesh::out << "Newton step did not reduce residual." << std::endl;
732 
733  // back off the previous step.
734  solution->add(-1., *delta_u);
735  solution->close();
736  *continuation_parameter -= delta_lambda;
737 
738  // Backtracking: start cutting the Newton stepsize by halves until
739  // the new residual is actually smaller...
740  for (unsigned int backtrack_step=0; backtrack_step<n_backtrack_steps; ++backtrack_step)
741  {
742  newton_stepfactor *= 0.5;
743 
744  if (!quiet)
745  libMesh::out << "Shrinking step size by " << newton_stepfactor << std::endl;
746 
747  // Take fractional step
748  solution->add(newton_stepfactor, *delta_u);
749  solution->close();
750  *continuation_parameter += newton_stepfactor*delta_lambda;
751 
752  rhs_mode = Residual;
753  assembly(true, // Residual
754  false); // Jacobian
755  rhs->close();
756  nonlinear_residual_afterstep = rhs->l2_norm();
757 
758  if (!quiet)
759  libMesh::out << "At shrink step "
760  << backtrack_step
761  << ", nonlinear_residual_afterstep="
762  << nonlinear_residual_afterstep
763  << std::endl;
764 
765  if (nonlinear_residual_afterstep < nonlinear_residual_beforestep)
766  {
767  if (!quiet)
768  libMesh::out << "Backtracking succeeded!" << std::endl;
769 
770  break; // out of backtracking loop
771  }
772 
773  else
774  {
775  // Back off that step
776  solution->add(-newton_stepfactor, *delta_u);
777  solution->close();
778  *continuation_parameter -= newton_stepfactor*delta_lambda;
779  }
780 
781  // Save a copy of the solution from before the Newton step.
782  //std::unique_ptr<NumericVector<Number>> prior_iterate = solution->clone();
783  }
784  } // end if (attempte_backtracking)
785 
786 
787  // If we tried backtracking but the residual is still not reduced, print message.
788  if ((attempt_backtracking) && (nonlinear_residual_afterstep > nonlinear_residual_beforestep))
789  {
790  //libMesh::err << "Backtracking failed." << std::endl;
791  libMesh::out << "Backtracking failed." << std::endl;
792 
793  // 1.) Quit, exit program.
794  //libmesh_error_msg("Backtracking failed!");
795 
796  // 2.) Continue with last newton_stepfactor
797  if (newton_step<3)
798  {
799  solution->add(newton_stepfactor, *delta_u);
800  solution->close();
801  *continuation_parameter += newton_stepfactor*delta_lambda;
802  if (!quiet)
803  libMesh::out << "Backtracking could not reduce residual ... continuing anyway!" << std::endl;
804  }
805 
806  // 3.) Break out of Newton iteration loop with newton_converged = false,
807  // reduce the arclength stepsize, and try again.
808  else
809  {
810  break; // out of Newton iteration loop, with newton_converged=false
811  }
812  }
813 
814  // Another type of convergence check: suppose the residual has not been reduced
815  // from its initial value after half of the allowed Newton steps have occurred.
816  // In our experience, this typically means that it isn't going to converge and
817  // we could probably save time by dropping out of the Newton iteration loop and
818  // trying a smaller arcstep.
819  if (this->newton_progress_check)
820  {
821  if ((nonlinear_residual_afterstep > nonlinear_residual_firststep) &&
822  (newton_step+1 > static_cast<unsigned int>(0.5*newton_solver->max_nonlinear_iterations)))
823  {
824  libMesh::out << "Progress check failed: the current residual: "
825  << nonlinear_residual_afterstep
826  << ", is\n"
827  << "larger than the initial residual, and half of the allowed\n"
828  << "number of Newton iterations have elapsed.\n"
829  << "Exiting Newton iterations with converged==false." << std::endl;
830 
831  break; // out of Newton iteration loop, newton_converged = false
832  }
833  }
834 
835  // Safety check: Check the current continuation parameter against user-provided min-allowable parameter value
837  {
838  libMesh::out << "Continuation parameter fell below min-allowable value." << std::endl;
839  break; // out of Newton iteration loop, newton_converged = false
840  }
841 
842  // Safety check: Check the current continuation parameter against user-provided max-allowable parameter value
843  if ( (max_continuation_parameter != 0.0) &&
845  {
846  libMesh::out << "Current continuation parameter value: "
848  << " exceeded max-allowable value."
849  << std::endl;
850  break; // out of Newton iteration loop, newton_converged = false
851  }
852 
853 
854  // Check the convergence of the parameter and the solution. If they are small
855  // enough, we can break out of the Newton iteration loop.
856  const Real norm_delta_u = delta_u->l2_norm();
857  const Real norm_u = solution->l2_norm();
858  libMesh::out << " delta_lambda = " << delta_lambda << std::endl;
859  libMesh::out << " newton_stepfactor*delta_lambda = " << newton_stepfactor*delta_lambda << std::endl;
860  libMesh::out << " lambda_current = " << *continuation_parameter << std::endl;
861  libMesh::out << " ||delta_u|| = " << norm_delta_u << std::endl;
862  libMesh::out << " ||delta_u||/||u|| = " << norm_delta_u / norm_u << std::endl;
863 
864 
865  // Evaluate the residual at the current Newton iterate. We don't want to detect
866  // convergence due to a small Newton step when the residual is still not small.
867  rhs_mode = Residual;
868  assembly(true, // Residual
869  false); // Jacobian
870  rhs->close();
871  const Real norm_residual = rhs->l2_norm();
872  libMesh::out << " ||R||_{L2} = " << norm_residual << std::endl;
873  libMesh::out << " ||R||_{L2}/||u|| = " << norm_residual / norm_u << std::endl;
874 
875 
876  // FIXME: The norm_delta_u tolerance (at least) should be relative.
877  // It doesn't make sense to converge a solution whose size is ~ 10^5 to
878  // a tolerance of 1.e-6. Oh, and we should also probably check the
879  // (relative) size of the residual as well, instead of just the step.
880  if ((std::abs(delta_lambda) < continuation_parameter_tolerance) &&
881  //(norm_delta_u < solution_tolerance) && // This is a *very* strict criterion we can probably skip
882  (norm_residual < solution_tolerance))
883  {
884  if (!quiet)
885  libMesh::out << "Newton iterations converged!" << std::endl;
886 
887  newton_converged = true;
888  break; // out of Newton iterations
889  }
890  } // end nonlinear loop
891 
892  if (!newton_converged)
893  {
894  libMesh::out << "Newton iterations of augmented system did not converge!" << std::endl;
895 
896  // Reduce ds_current, recompute the solution and parameter, and continue to next
897  // arcstep, if there is one.
898  ds_current *= 0.5;
899 
900  // Go back to previous solution and parameter value.
901  *solution = *previous_u;
903 
904  // Compute new predictor with smaller ds
905  apply_predictor();
906  }
907  else
908  {
909  // Set step convergence and break out
910  arcstep_converged=true;
911  break; // out of arclength reduction loop
912  }
913 
914  } // end loop over arclength reductions
915 
916  // Check for convergence of the whole arcstep. If not converged at this
917  // point, we have no choice but to quit.
918  if (!arcstep_converged)
919  libmesh_error_msg("Arcstep failed to converge after max number of reductions! Exiting...");
920 
921  // Print converged solution control parameter and max value.
922  libMesh::out << "lambda_current=" << *continuation_parameter << std::endl;
923  //libMesh::out << "u_max=" << solution->max() << std::endl;
924 
925  // Reset old stream precision and flags.
926  libMesh::out.precision(old_precision);
927  libMesh::out.unsetf(std::ios_base::scientific);
928 
929  // Note: we don't want to go on to the next guess yet, since the user may
930  // want to post-process this data. It's up to the user to call advance_arcstep()
931  // when they are ready to go on.
932 }
T libmesh_real(T a)
double abs(double a)
NumericVector< Number > * previous_u
std::unique_ptr< LinearSolver< Number > > linear_solver
virtual Real l2_norm() const =0
unsigned int max_nonlinear_iterations
Definition: diff_solver.h:157
std::unique_ptr< TimeSolver > time_solver
Definition: diff_system.h:234
NumericVector< Number > * delta_u
NumericVector< Number > * rhs
virtual void zero()=0
virtual void scale(const T factor)=0
NumericVector< Number > * du_ds
virtual void assembly(bool get_residual, bool get_jacobian, bool apply_heterogeneous_constraints=false, bool apply_no_constraints=false) libmesh_override
Definition: fem_system.C:852
NumericVector< Number > * z
unsigned int max_linear_iterations
Definition: diff_solver.h:149
NumericVector< Number > * y_old
NumericVector< Number > * y
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1509
void unsetf(std::ios_base::fmtflags mask)
double pow(double a, int b)
virtual void close()=0
virtual void close()=0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::ios_base::fmtflags setf(std::ios_base::fmtflags fmtfl)
SparseMatrix< Number > * matrix
virtual T dot(const NumericVector< T > &v) const =0
virtual void add(const numeric_index_type i, const T value)=0
OStreamProxy out(std::cout)
long double min(long double a, double b)
std::streamsize precision() const
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 192 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::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().

193 {
194  // Check the sizes
195  libmesh_assert_less (global_dof_number, _dof_map->n_dofs());
196  libmesh_assert_less (global_dof_number, current_local_solution->size());
197 
198  return (*current_local_solution)(global_dof_number);
199 }
std::unique_ptr< DofMap > _dof_map
Definition: system.h:1865
std::unique_ptr< NumericVector< Number > > current_local_solution
Definition: system.h:1521
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().

375  {
376  return request_jacobian;
377  }
void libMesh::System::deactivate ( )
inlineinherited

Deactivates the system. Only active systems are solved.

Definition at line 2062 of file system.h.

References libMesh::System::_active.

Referenced by libMesh::System::get_equation_systems().

2063 {
2064  _active = false;
2065 }
void libMesh::ImplicitSystem::disable_cache ( )
virtualinherited

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 313 of file implicit_system.C.

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

313  {
314  this->assemble_before_solve = true;
315  this->get_linear_solver()->reuse_preconditioner(false);
316 }
virtual LinearSolver< Number > * get_linear_solver() const
virtual void reuse_preconditioner(bool)
bool assemble_before_solve
Definition: system.h:1463
void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited
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().

143  {
144  return request_jacobian;
145  }
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 {}
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  {}
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  {}
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().

125  {
126  return request_jacobian;
127  }
void libMesh::ReferenceCounter::enable_print_counter_info ( )
staticinherited

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

Definition at line 101 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

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

102 {
103  _enable_print_counter = true;
104  return;
105 }
virtual bool libMesh::FEMPhysics::eulerian_residual ( bool  request_jacobian,
DiffContext context 
)
virtualinherited

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.

Referenced by libMesh::FEMPhysics::~FEMPhysics().

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.

295  {
296  return request_jacobian;
297  }
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(), and libMesh::DifferentiableQoI::init_context().

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

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 815 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::qoi, libMesh::Real, libMesh::ExplicitSystem::rhs, libMesh::ImplicitSystem::sensitivity_solve(), libMesh::ParameterVector::size(), and libMesh::TOLERANCE.

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

818 {
819  ParameterVector & parameters =
820  const_cast<ParameterVector &>(parameters_in);
821 
822  const unsigned int Np = cast_int<unsigned int>
823  (parameters.size());
824  const unsigned int Nq = cast_int<unsigned int>
825  (qoi.size());
826 
827  // An introduction to the problem:
828  //
829  // Residual R(u(p),p) = 0
830  // partial R / partial u = J = system matrix
831  //
832  // This implies that:
833  // d/dp(R) = 0
834  // (partial R / partial p) +
835  // (partial R / partial u) * (partial u / partial p) = 0
836 
837  // We first solve for (partial u / partial p) for each parameter:
838  // J * (partial u / partial p) = - (partial R / partial p)
839 
840  this->sensitivity_solve(parameters);
841 
842  // Get ready to fill in sensitivities:
843  sensitivities.allocate_data(qoi_indices, *this, parameters);
844 
845  // We use the identity:
846  // dq/dp = (partial q / partial p) + (partial q / partial u) *
847  // (partial u / partial p)
848 
849  // We get (partial q / partial u) from the user
850  this->assemble_qoi_derivative(qoi_indices,
851  /* include_liftfunc = */ true,
852  /* apply_constraints = */ false);
853 
854  // We don't need these to be closed() in this function, but libMesh
855  // standard practice is to have them closed() by the time the
856  // function exits
857  for (std::size_t i=0; i != this->qoi.size(); ++i)
858  if (qoi_indices.has_index(i))
859  this->get_adjoint_rhs(i).close();
860 
861  for (unsigned int j=0; j != Np; ++j)
862  {
863  // We currently get partial derivatives via central differencing
864 
865  // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
866 
867  Number old_parameter = *parameters[j];
868 
869  const Real delta_p =
870  TOLERANCE * std::max(std::abs(old_parameter), 1e-3);
871 
872  *parameters[j] = old_parameter - delta_p;
873  this->assemble_qoi(qoi_indices);
874  std::vector<Number> qoi_minus = this->qoi;
875 
876  *parameters[j] = old_parameter + delta_p;
877  this->assemble_qoi(qoi_indices);
878  std::vector<Number> & qoi_plus = this->qoi;
879 
880  std::vector<Number> partialq_partialp(Nq, 0);
881  for (unsigned int i=0; i != Nq; ++i)
882  if (qoi_indices.has_index(i))
883  partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
884 
885  // Don't leave the parameter changed
886  *parameters[j] = old_parameter;
887 
888  for (unsigned int i=0; i != Nq; ++i)
889  if (qoi_indices.has_index(i))
890  sensitivities[i][j] = partialq_partialp[i] +
891  this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(j));
892  }
893 
894  // All parameters have been reset.
895  // We didn't cache the original rhs or matrix for memory reasons,
896  // but we can restore them to a state consistent solution -
897  // principle of least surprise.
898  this->assembly(true, true);
899  this->rhs->close();
900  this->matrix->close();
901  this->assemble_qoi(qoi_indices);
902 }
double abs(double a)
NumericVector< Number > & get_sensitivity_solution(unsigned int i=0)
Definition: system.C:917
virtual std::pair< unsigned int, Real > sensitivity_solve(const ParameterVector &parameters) libmesh_override
NumericVector< Number > * rhs
static const Real TOLERANCE
long double max(long double a, double b)
virtual void assemble_qoi_derivative(const QoISet &qoi_indices=QoISet(), bool include_liftfunc=true, bool apply_constraints=true) libmesh_override
virtual void assembly(bool, bool, bool=false, bool=false)
std::vector< Number > qoi
Definition: system.h:1539
virtual void close()=0
virtual void close()=0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
SparseMatrix< Number > * matrix
virtual T dot(const NumericVector< T > &v) const =0
virtual void assemble_qoi(const QoISet &qoi_indices=QoISet()) libmesh_override
NumericVector< Number > & get_adjoint_rhs(unsigned int i=0)
Definition: system.C:1032
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 1032 of file system.C.

References libMesh::System::get_vector().

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

1033 {
1034  std::ostringstream adjoint_rhs_name;
1035  adjoint_rhs_name << "adjoint_rhs" << i;
1036 
1037  return this->get_vector(adjoint_rhs_name.str());
1038 }
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:775
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 1042 of file system.C.

References libMesh::System::get_vector().

1043 {
1044  std::ostringstream adjoint_rhs_name;
1045  adjoint_rhs_name << "adjoint_rhs" << i;
1046 
1047  return this->get_vector(adjoint_rhs_name.str());
1048 }
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:775
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 970 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::System::project_solution_on_reinit(), libMesh::ImplicitSystem::qoi_parameter_hessian(), libMesh::ImplicitSystem::qoi_parameter_hessian_vector_product(), and libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve().

971 {
972  std::ostringstream adjoint_name;
973  adjoint_name << "adjoint_solution" << i;
974 
975  return this->get_vector(adjoint_name.str());
976 }
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:775
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 980 of file system.C.

References libMesh::System::get_vector().

981 {
982  std::ostringstream adjoint_name;
983  adjoint_name << "adjoint_solution" << i;
984 
985  return this->get_vector(adjoint_name.str());
986 }
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:775
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 1259 of file system.C.

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

Referenced by libMesh::System::project_solution_on_reinit().

1260 {
1261  all_variable_numbers.resize(n_vars());
1262 
1263  // Make sure the variable exists
1264  std::map<std::string, unsigned short int>::const_iterator
1265  it = _variable_numbers.begin();
1266  std::map<std::string, unsigned short int>::const_iterator
1267  it_end = _variable_numbers.end();
1268 
1269  unsigned int count = 0;
1270  for ( ; it != it_end; ++it)
1271  {
1272  all_variable_numbers[count] = it->second;
1273  count++;
1274  }
1275 }
std::map< std::string, unsigned short int > _variable_numbers
Definition: system.h:1903
unsigned int n_vars() const
Definition: system.h:2086
const DofMap & libMesh::System::get_dof_map ( ) const
inlineinherited
Returns
A constant reference to this system's _dof_map.

Definition at line 2030 of file system.h.

References libMesh::System::_dof_map.

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), libMesh::__libmesh_petsc_diff_solver_residual(), libMesh::__libmesh_petsc_snes_postcheck(), libMesh::ExactSolution::_compute_error(), libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::DifferentiableSystem::add_dot_var_dirichlet_bcs(), libMesh::HPCoarsenTest::add_projection(), libMesh::UnsteadySolver::adjoint_advance_timestep(), libMesh::ImplicitSystem::adjoint_solve(), libMesh::NewmarkSolver::advance_timestep(), libMesh::UnsteadySolver::advance_timestep(), libMesh::EquationSystems::allgather(), libMesh::EquationSystems::build_discontinuous_solution_vector(), libMesh::EquationSystems::build_parallel_solution_vector(), libMesh::System::calculate_norm(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), DMCreateDomainDecomposition_libMesh(), DMCreateFieldDecomposition_libMesh(), DMlibMeshFunction(), DMlibMeshJacobian(), DMlibMeshSetSystem_libMesh(), libMesh::DofMap::enforce_constraints_exactly(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::System::get_info(), libMesh::EquationSystems::get_solution(), libMesh::SystemSubsetBySubdomain::init(), libMesh::SecondOrderUnsteadySolver::init_data(), libMesh::UnsteadySolver::init_data(), libMesh::EigenSystem::init_matrices(), libMesh::ImplicitSystem::init_matrices(), libMesh::CondensedEigenSystem::initialize_condensed_dofs(), libMesh::OptimizationSystem::initialize_equality_constraints_storage(), libMesh::OptimizationSystem::initialize_inequality_constraints_storage(), libMesh::System::local_dof_indices(), libMesh::DofMap::max_constraint_error(), libMesh::DGFEMContext::neighbor_side_fe_reinit(), libMesh::UnsteadySolver::old_nonlinear_solution(), libMesh::SecondOrderUnsteadySolver::old_solution_accel(), libMesh::SecondOrderUnsteadySolver::old_solution_rate(), libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::operator()(), libMesh::petsc_auto_fieldsplit(), libMesh::ErrorVector::plot_error(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::FEMContext::pre_fe_reinit(), libMesh::System::re_update(), libMesh::SecondOrderUnsteadySolver::reinit(), libMesh::UnsteadySolver::reinit(), libMesh::ImplicitSystem::reinit(), libMesh::EigenSystem::reinit(), libMesh::System::reinit_constraints(), libMesh::EquationSystems::reinit_solutions(), libMesh::UnsteadySolver::retrieve_timestep(), libMesh::HPCoarsenTest::select_refinement(), libMesh::ImplicitSystem::sensitivity_solve(), libMesh::NewtonSolver::solve(), libMesh::PetscDiffSolver::solve(), libMesh::PetscNonlinearSolver< T >::solve(), libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve(), libMesh::ImplicitSystem::weighted_sensitivity_solve(), libMesh::EnsightIO::write_scalar_ascii(), and libMesh::EnsightIO::write_vector_ascii().

2031 {
2032  return *_dof_map;
2033 }
std::unique_ptr< DofMap > _dof_map
Definition: system.h:1865
DofMap & libMesh::System::get_dof_map ( )
inlineinherited
Returns
A writable reference to this system's _dof_map.

Definition at line 2038 of file system.h.

References libMesh::System::_dof_map.

2039 {
2040  return *_dof_map;
2041 }
std::unique_ptr< DofMap > _dof_map
Definition: system.h:1865
EquationSystems& libMesh::System::get_equation_systems ( )
inlineinherited
Returns
A reference to this system's parent EquationSystems object.

Definition at line 703 of file system.h.

References libMesh::System::_equation_systems, libMesh::System::activate(), libMesh::System::active(), libMesh::System::deactivate(), and libMesh::System::set_basic_system_only().

703 { return _equation_systems; }
EquationSystems & _equation_systems
Definition: system.h:1871
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; }
std::set< unsigned int > _first_order_vars
Definition: diff_physics.h:559
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().

48 {
49 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
50 
51  std::ostringstream oss;
52 
53  oss << '\n'
54  << " ---------------------------------------------------------------------------- \n"
55  << "| Reference count information |\n"
56  << " ---------------------------------------------------------------------------- \n";
57 
58  for (Counts::iterator it = _counts.begin();
59  it != _counts.end(); ++it)
60  {
61  const std::string name(it->first);
62  const unsigned int creations = it->second.first;
63  const unsigned int destructions = it->second.second;
64 
65  oss << "| " << name << " reference count information:\n"
66  << "| Creations: " << creations << '\n'
67  << "| Destructions: " << destructions << '\n';
68  }
69 
70  oss << " ---------------------------------------------------------------------------- \n";
71 
72  return oss.str();
73 
74 #else
75 
76  return "";
77 
78 #endif
79 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
std::string libMesh::System::get_info ( ) const
inherited
Returns
A string containing information about the system.

Definition at line 1645 of file system.C.

References libMesh::FEType::family, libMesh::System::get_dof_map(), libMesh::DofMap::get_info(), libMesh::FEType::inf_map, libMesh::System::n_constrained_dofs(), libMesh::System::n_dofs(), libMesh::System::n_local_constrained_dofs(), libMesh::System::n_local_dofs(), libMesh::System::n_matrices(), libMesh::System::n_variable_groups(), libMesh::VariableGroup::n_variables(), libMesh::System::n_vectors(), libMesh::VariableGroup::name(), libMesh::System::name(), libMesh::System::number(), libMesh::FEType::order, libMesh::FEType::radial_family, libMesh::FEType::radial_order, libMesh::System::system_type(), libMesh::Variable::type(), libMesh::DofMap::variable_group(), and libMesh::System::variable_group().

Referenced by libMesh::System::read_parallel_data().

1646 {
1647  std::ostringstream oss;
1648 
1649 
1650  const std::string & sys_name = this->name();
1651 
1652  oss << " System #" << this->number() << ", \"" << sys_name << "\"\n"
1653  << " Type \"" << this->system_type() << "\"\n"
1654  << " Variables=";
1655 
1656  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1657  {
1658  const VariableGroup & vg_description (this->variable_group(vg));
1659 
1660  if (vg_description.n_variables() > 1) oss << "{ ";
1661  for (unsigned int vn=0; vn<vg_description.n_variables(); vn++)
1662  oss << "\"" << vg_description.name(vn) << "\" ";
1663  if (vg_description.n_variables() > 1) oss << "} ";
1664  }
1665 
1666  oss << '\n';
1667 
1668  oss << " Finite Element Types=";
1669 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
1670  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1671  oss << "\""
1672  << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
1673  << "\" ";
1674 #else
1675  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1676  {
1677  oss << "\""
1678  << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
1679  << "\", \""
1680  << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().radial_family)
1681  << "\" ";
1682  }
1683 
1684  oss << '\n' << " Infinite Element Mapping=";
1685  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1686  oss << "\""
1687  << Utility::enum_to_string<InfMapType>(this->get_dof_map().variable_group(vg).type().inf_map)
1688  << "\" ";
1689 #endif
1690 
1691  oss << '\n';
1692 
1693  oss << " Approximation Orders=";
1694  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1695  {
1696 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
1697  oss << "\""
1698  << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
1699  << "\" ";
1700 #else
1701  oss << "\""
1702  << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
1703  << "\", \""
1704  << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().radial_order)
1705  << "\" ";
1706 #endif
1707  }
1708 
1709  oss << '\n';
1710 
1711  oss << " n_dofs()=" << this->n_dofs() << '\n';
1712  oss << " n_local_dofs()=" << this->n_local_dofs() << '\n';
1713 #ifdef LIBMESH_ENABLE_CONSTRAINTS
1714  oss << " n_constrained_dofs()=" << this->n_constrained_dofs() << '\n';
1715  oss << " n_local_constrained_dofs()=" << this->n_local_constrained_dofs() << '\n';
1716 #endif
1717 
1718  oss << " " << "n_vectors()=" << this->n_vectors() << '\n';
1719  oss << " " << "n_matrices()=" << this->n_matrices() << '\n';
1720  // oss << " " << "n_additional_matrices()=" << this->n_additional_matrices() << '\n';
1721 
1722  oss << this->get_dof_map().get_info();
1723 
1724  return oss.str();
1725 }
FEFamily family
Definition: fe_type.h:204
dof_id_type n_constrained_dofs() const
Definition: system.C:155
const FEType & type() const
Definition: variable.h:119
std::string get_info() const
Definition: dof_map.C:2642
OrderWrapper radial_order
Definition: fe_type.h:237
OrderWrapper order
Definition: fe_type.h:198
const VariableGroup & variable_group(unsigned int vg) const
Definition: system.h:2124
const VariableGroup & variable_group(const unsigned int c) const
Definition: dof_map.h:1679
const std::string & name() const
Definition: system.h:1998
unsigned int n_variable_groups() const
Definition: system.h:2094
const DofMap & get_dof_map() const
Definition: system.h:2030
InfMapType inf_map
Definition: fe_type.h:258
virtual unsigned int n_matrices() const
Definition: system.h:2220
FEFamily radial_family
Definition: fe_type.h:250
dof_id_type n_local_dofs() const
Definition: system.C:185
virtual std::string system_type() const
Definition: system.h:473
dof_id_type n_local_constrained_dofs() const
Definition: system.C:170
unsigned int number() const
Definition: system.h:2006
dof_id_type n_dofs() const
Definition: system.C:148
unsigned int n_vectors() const
Definition: system.h:2214
std::pair< unsigned int, Real > libMesh::DifferentiableSystem::get_linear_solve_parameters ( ) const
virtualinherited
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.

185 {
186  libmesh_assert(time_solver.get());
187  libmesh_assert_equal_to (&(time_solver->system()), this);
188  return std::make_pair(this->time_solver->diff_solver()->max_linear_iterations,
189  this->time_solver->diff_solver()->relative_residual_tolerance);
190 }
std::unique_ptr< TimeSolver > time_solver
Definition: diff_system.h:234
LinearSolver< Number > * libMesh::DifferentiableSystem::get_linear_solver ( ) const
virtualinherited
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.

176 {
177  libmesh_assert(time_solver.get());
178  libmesh_assert_equal_to (&(time_solver->system()), this);
179  return this->time_solver->linear_solver().get();
180 }
std::unique_ptr< TimeSolver > time_solver
Definition: diff_system.h:234
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 270 of file implicit_system.C.

References libMesh::ImplicitSystem::_matrices.

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

271 {
272  // Make sure the matrix exists
273  const_matrices_iterator pos = _matrices.find (mat_name);
274 
275  if (pos == _matrices.end())
276  libmesh_error_msg("ERROR: matrix " << mat_name << " does not exist in this system!");
277 
278  return *(pos->second);
279 }
std::map< std::string, SparseMatrix< Number > * > _matrices
std::map< std::string, SparseMatrix< Number > * >::const_iterator const_matrices_iterator
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 283 of file implicit_system.C.

References libMesh::ImplicitSystem::_matrices.

284 {
285  // Make sure the matrix exists
286  matrices_iterator pos = _matrices.find (mat_name);
287 
288  if (pos == _matrices.end())
289  libmesh_error_msg("ERROR: matrix " << mat_name << " does not exist in this system!");
290 
291  return *(pos->second);
292 }
std::map< std::string, SparseMatrix< Number > * > _matrices
std::map< std::string, SparseMatrix< Number > * >::iterator matrices_iterator
const MeshBase & libMesh::System::get_mesh ( ) const
inlineinherited
Returns
A constant reference to this systems's _mesh.

Definition at line 2014 of file system.h.

References libMesh::System::_mesh.

Referenced by libMesh::ExactSolution::_compute_error(), libMesh::HPCoarsenTest::add_projection(), libMesh::FEMSystem::assemble_qoi(), libMesh::FEMSystem::assemble_qoi_derivative(), libMesh::FEMSystem::assembly(), libMesh::System::calculate_norm(), DMCreateDomainDecomposition_libMesh(), DMCreateFieldDecomposition_libMesh(), DMlibMeshSetSystem_libMesh(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointResidualErrorEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::SystemSubsetBySubdomain::init(), libMesh::System::init_data(), libMesh::EigenSystem::init_matrices(), libMesh::ImplicitSystem::init_matrices(), libMesh::System::local_dof_indices(), libMesh::DofMap::max_constraint_error(), libMesh::FEMSystem::mesh_position_get(), libMesh::FEMSystem::mesh_position_set(), libMesh::petsc_auto_fieldsplit(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::FEMSystem::postprocess(), libMesh::ImplicitSystem::reinit(), libMesh::EigenSystem::reinit(), libMesh::HPSingularity::select_refinement(), libMesh::HPCoarsenTest::select_refinement(), and libMesh::System::zero_variable().

2015 {
2016  return _mesh;
2017 }
MeshBase & _mesh
Definition: system.h:1877
MeshBase & libMesh::System::get_mesh ( )
inlineinherited
Returns
A reference to this systems's _mesh.

Definition at line 2022 of file system.h.

References libMesh::System::_mesh.

2023 {
2024  return _mesh;
2025 }
MeshBase & _mesh
Definition: system.h:1877
const System * libMesh::DifferentiablePhysics::get_mesh_system ( ) const
inlineinherited
Returns
A const reference to the system with variables corresponding to mesh nodal coordinates, or NULL 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(), and libMesh::DifferentiablePhysics::init_context().

625 {
626  return _mesh_sys;
627 }
System * libMesh::DifferentiablePhysics::get_mesh_system ( )
inlineinherited
Returns
A reference to the system with variables corresponding to mesh nodal coordinates, or NULL if the mesh is fixed.

Definition at line 630 of file diff_physics.h.

References libMesh::DifferentiablePhysics::_mesh_sys.

631 {
632  return _mesh_sys;
633 }
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(), and libMesh::DifferentiablePhysics::init_context().

637 {
638  return _mesh_x_var;
639 }
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(), and libMesh::DifferentiablePhysics::init_context().

643 {
644  return _mesh_y_var;
645 }
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(), and libMesh::DifferentiablePhysics::init_context().

649 {
650  return _mesh_z_var;
651 }
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; }
DifferentiablePhysics * _diff_physics
Definition: diff_system.h:371
const DifferentiableQoI* libMesh::DifferentiableSystem::get_qoi ( ) const