libMesh
Public Types | Public Member Functions | Static Public Member Functions | Public Attributes | Protected Types | Protected Member Functions | Static Protected Member Functions | Protected Attributes | Static Protected Attributes | List of all members
ElasticityRBConstruction Class Reference

#include <rb_classes.h>

Inheritance diagram for ElasticityRBConstruction:
[legend]

Public Types

typedef ElasticityRBConstruction sys_type
 The type of system. More...
 
typedef RBConstruction Parent
 The type of the parent. More...
 
typedef std::map< std::string, SparseMatrix< Number > * >::iterator matrices_iterator
 Matrix iterator typedefs. More...
 
typedef std::map< std::string, SparseMatrix< Number > * >::const_iterator const_matrices_iterator
 
typedef std::map< std::string, NumericVector< Number > * >::iterator vectors_iterator
 Vector iterator typedefs. More...
 
typedef std::map< std::string, NumericVector< Number > * >::const_iterator const_vectors_iterator
 

Public Member Functions

 ElasticityRBConstruction (EquationSystems &es, const std::string &name_in, const unsigned int number_in)
 
virtual ~ElasticityRBConstruction ()
 Destructor. More...
 
virtual void init_data ()
 Initialize data structures. More...
 
virtual void init_context (FEMContext &c)
 Pre-request all relevant element data. More...
 
virtual void solve_for_matrix_and_rhs (LinearSolver< Number > &input_solver, SparseMatrix< Number > &input_matrix, NumericVector< Number > &input_rhs)
 Assembles & solves the linear system A*x=b for the specified matrix input_matrix and right-hand side rhs. More...
 
void set_rb_evaluation (RBEvaluation &rb_eval_in)
 Set the RBEvaluation object. More...
 
RBEvaluation & get_rb_evaluation ()
 Get a reference to the RBEvaluation object. More...
 
bool is_rb_eval_initialized () const
 
RBThetaExpansion & get_rb_theta_expansion ()
 Get a reference to the RBThetaExpansion object that that belongs to rb_eval. More...
 
void set_rb_assembly_expansion (RBAssemblyExpansion &rb_assembly_expansion_in)
 Set the rb_assembly_expansion object. More...
 
RBAssemblyExpansion & get_rb_assembly_expansion ()
 
sys_typesystem ()
 
virtual void clear () libmesh_override
 Clear all the data structures associated with the system. More...
 
virtual std::string system_type () const libmesh_override
 
virtual Real truth_solve (int plot_solution)
 Perform a "truth" solve, i.e. More...
 
virtual Real train_reduced_basis (const bool resize_rb_eval_data=true)
 Train the reduced basis. More...
 
virtual Real compute_max_error_bound ()
 (i) Compute the a posteriori error bound for each set of parameters in the training set, (ii) set current_parameters to the parameters that maximize the error bound, and (iii) return the maximum error bound. More...
 
const RBParameters & get_greedy_parameter (unsigned int i)
 Return the parameters chosen during the i^th step of the Greedy algorithm. More...
 
void set_rel_training_tolerance (Real new_training_tolerance)
 Get/set the relative tolerance for the basis training. More...
 
Real get_rel_training_tolerance ()
 
void set_abs_training_tolerance (Real new_training_tolerance)
 Get/set the absolute tolerance for the basis training. More...
 
Real get_abs_training_tolerance ()
 
void set_normalize_rb_bound_in_greedy (bool normalize_rb_bound_in_greedy)
 Get/set the boolean to indicate if we normalize the RB error in the greedy. More...
 
bool get_normalize_rb_bound_in_greedy ()
 
unsigned int get_Nmax () const
 Get/set Nmax, the maximum number of RB functions we are willing to compute. More...
 
virtual void set_Nmax (unsigned int Nmax)
 
void set_quiet_mode (bool quiet_mode_in)
 Set the quiet_mode flag. More...
 
bool is_quiet () const
 Is the system in quiet mode? More...
 
virtual void load_basis_function (unsigned int i)
 Load the i^th RB function into the RBConstruction solution vector. More...
 
virtual void load_rb_solution ()
 Load the RB solution from the most recent solve with rb_eval into this system's solution vector. More...
 
SparseMatrix< Number > * get_inner_product_matrix ()
 Get a pointer to inner_product_matrix. More...
 
SparseMatrix< Number > * get_Aq (unsigned int q)
 Get a pointer to Aq. More...
 
SparseMatrix< Number > * get_non_dirichlet_Aq (unsigned int q)
 Get a pointer to non_dirichlet_Aq. More...
 
virtual void initialize_rb_construction (bool skip_matrix_assembly=false, bool skip_vector_assembly=false)
 Allocate all the data structures necessary for the construction stage of the RB method. More...
 
NumericVector< Number > * get_Fq (unsigned int q)
 Get a pointer to Fq. More...
 
NumericVector< Number > * get_non_dirichlet_Fq (unsigned int q)
 Get a pointer to non-Dirichlet Fq. More...
 
NumericVector< Number > * get_output_vector (unsigned int n, unsigned int q_l)
 Get a pointer to the n^th output. More...
 
NumericVector< Number > * get_non_dirichlet_output_vector (unsigned int n, unsigned int q_l)
 Get a pointer to non-Dirichlet output vector. More...
 
virtual void get_all_matrices (std::map< std::string, SparseMatrix< Number > * > &all_matrices)
 Get a map that stores pointers to all of the matrices. More...
 
virtual void get_all_vectors (std::map< std::string, NumericVector< Number > * > &all_vectors)
 Get a map that stores pointers to all of the vectors. More...
 
virtual void get_output_vectors (std::map< std::string, NumericVector< Number > * > &all_vectors)
 Get a map that stores pointers to all of the vectors. More...
 
virtual void assemble_affine_expansion (bool skip_matrix_assembly, bool skip_vector_assembly)
 Assemble the matrices and vectors for this system. More...
 
void assemble_inner_product_matrix (SparseMatrix< Number > *input_matrix, bool apply_dof_constraints=true)
 Assemble the inner product matrix and store it in input_matrix. More...
 
void assemble_Aq_matrix (unsigned int q, SparseMatrix< Number > *input_matrix, bool apply_dof_constraints=true)
 Assemble the q^th affine matrix and store it in input_matrix. More...
 
void assemble_Fq_vector (unsigned int q, NumericVector< Number > *input_vector, bool apply_dof_constraints=true)
 Assemble the q^th affine vector and store it in input_matrix. More...
 
void add_scaled_Aq (Number scalar, unsigned int q_a, SparseMatrix< Number > *input_matrix, bool symmetrize)
 Add the scaled q^th affine matrix to input_matrix. More...
 
virtual void write_riesz_representors_to_files (const std::string &riesz_representors_dir, const bool write_binary_residual_representors)
 Write out all the Riesz representor data to files. More...
 
virtual void read_riesz_representors_from_files (const std::string &riesz_representors_dir, const bool write_binary_residual_representors)
 Read in all the Riesz representor data from files. More...
 
virtual void recompute_all_residual_terms (const bool compute_inner_products=true)
 This function computes all of the residual representors, can be useful when restarting a basis training computation. More...
 
virtual void process_parameters_file (const std::string &parameters_filename)
 Read in from the file specified by parameters_filename and set the this system's member variables accordingly. More...
 
void set_rb_construction_parameters (unsigned int n_training_samples_in, bool deterministic_training_in, unsigned int training_parameters_random_seed_in, bool quiet_mode_in, unsigned int Nmax_in, Real rel_training_tolerance_in, Real abs_training_tolerance_in, bool normalize_rb_error_bound_in_greedy_in, RBParameters mu_min_in, RBParameters mu_max_in, std::map< std::string, std::vector< Real >> discrete_parameter_values_in, std::map< std::string, bool > log_scaling)
 Set the state of this RBConstruction object based on the arguments to this function. More...
 
virtual void print_info ()
 Print out info that describes the current setup of this RBConstruction. More...
 
void print_basis_function_orthogonality ()
 Print out a matrix that shows the orthogonality of the RB basis functions. More...
 
unsigned int get_delta_N () const
 Get delta_N, the number of basis functions we add to the RB space per iteration of the greedy algorithm. More...
 
void set_inner_product_assembly (ElemAssembly &inner_product_assembly_in)
 Set the rb_assembly_expansion object. More...
 
ElemAssembly & get_inner_product_assembly ()
 
void zero_constrained_dofs_on_vector (NumericVector< Number > &vector)
 It is sometimes useful to be able to zero vector entries that correspond to constrained dofs. More...
 
void set_convergence_assertion_flag (bool flag)
 Setter for the flag determining if convergence should be checked after each solve. More...
 
numeric_index_type get_n_training_samples () const
 Get the total number of training samples. More...
 
numeric_index_type get_local_n_training_samples () const
 Get the total number of training samples local to this processor. More...
 
numeric_index_type get_first_local_training_index () const
 Get the first local index of the training parameters. More...
 
numeric_index_type get_last_local_training_index () const
 Get the last local index of the training parameters. More...
 
virtual void initialize_training_parameters (const RBParameters &mu_min, const RBParameters &mu_max, unsigned int n_training_parameters, std::map< std::string, bool > log_param_scale, bool deterministic=true)
 Initialize the parameter ranges and indicate whether deterministic or random training parameters should be used and whether or not we want the parameters to be scaled logarithmically. More...
 
virtual void load_training_set (std::map< std::string, std::vector< Number >> &new_training_set)
 Overwrite the training parameters with new_training_set. More...
 
void broadcast_parameters (unsigned int proc_id)
 Broadcasts parameters on processor proc_id to all processors. More...
 
void set_training_random_seed (unsigned int seed)
 Set the seed that is used to randomly generate training parameters. More...
 
void set_deterministic_training_parameter_name (const std::string &name)
 In some cases we only want to allow discrete parameter values, instead of parameters that may take any value in a specified interval. More...
 
const std::string & get_deterministic_training_parameter_name () const
 Get the name of the parameter that we will generate deterministic training parameters for. More...
 
void set_deterministic_training_parameter_repeats (unsigned int repeats)
 Set the number of times each sample of the deterministic training parameter is repeated. More...
 
unsigned int get_deterministic_training_parameter_repeats () const
 Get the number of times each sample of the deterministic training parameter is repeated. More...
 
virtual void reinit () libmesh_override
 Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be used. More...
 
virtual void assemble () libmesh_override
 Prepares matrix and _dof_map for matrix assembly. More...
 
virtual void restrict_solve_to (const SystemSubset *subset, const SubsetSolveMode subset_solve_mode=SUBSET_ZERO) libmesh_override
 After calling this method, any solve will be limited to the given subset. More...
 
virtual void solve () libmesh_override
 Assembles & solves the linear system A*x=b. More...
 
virtual LinearSolver< Number > * get_linear_solver () const libmesh_override
 
virtual void release_linear_solver (LinearSolver< Number > *) const libmesh_override
 Releases a pointer to a linear solver acquired by this->get_linear_solver() More...
 
virtual void assembly (bool get_residual, bool get_jacobian, bool apply_heterogeneous_constraints=false, bool apply_no_constraints=false) libmesh_override
 Assembles a residual in rhs and/or a jacobian in matrix, as requested. More...
 
unsigned int n_linear_iterations () const
 
Real final_linear_residual () const
 
void attach_shell_matrix (ShellMatrix< Number > *shell_matrix)
 This function enables the user to provide a shell matrix, i.e. More...
 
void detach_shell_matrix ()
 Detaches a shell matrix. More...
 
ShellMatrix< Number > * get_shell_matrix ()
 
virtual void disable_cache () libmesh_override
 Avoids use of any cached data that might affect any solve result. More...
 
virtual std::pair< unsigned int, Real > get_linear_solve_parameters () const
 
virtual void assemble_residual_derivatives (const ParameterVector &parameters) libmesh_override
 Residual parameter derivative function. More...
 
virtual std::pair< unsigned int, Real > sensitivity_solve (const ParameterVector &parameters) libmesh_override
 Assembles & solves the linear system(s) (dR/du)*u_p = -dR/dp, for those parameters contained within parameters. More...
 
virtual std::pair< unsigned int, Real > weighted_sensitivity_solve (const ParameterVector &parameters, const ParameterVector &weights) libmesh_override
 Assembles & solves the linear system(s) (dR/du)*u_w = sum(w_p*-dR/dp), for those parameters p contained within parameters weighted by the values w_p found within weights. More...
 
virtual std::pair< unsigned int, Real > adjoint_solve (const QoISet &qoi_indices=QoISet()) libmesh_override
 Assembles & solves the linear system (dR/du)^T*z = dq/du, for those quantities of interest q specified by qoi_indices. More...
 
virtual std::pair< unsigned int, Real > weighted_sensitivity_adjoint_solve (const ParameterVector &parameters, const ParameterVector &weights, const QoISet &qoi_indices=QoISet()) libmesh_override
 Assembles & solves the linear system(s) (dR/du)^T*z_w = sum(w_p*(d^2q/dudp - d^2R/dudp*z)), for those parameters p contained within parameters, weighted by the values w_p found within weights. More...
 
virtual void adjoint_qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities) libmesh_override
 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]. More...
 
virtual void forward_qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities) libmesh_override
 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]. More...
 
virtual void qoi_parameter_hessian (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &hessian) libmesh_override
 For each of the system's quantities of interest q in qoi[qoi_indices], and for a vector of parameters p, the parameter sensitivity Hessian H_ij is defined as H_ij = (d^2 q)/(d p_i d p_j) This Hessian is the output of this method, where for each q_i, H_jk is stored in hessian.second_derivative(i,j,k). More...
 
virtual void qoi_parameter_hessian_vector_product (const QoISet &qoi_indices, const ParameterVector &parameters, const ParameterVector &vector, SensitivityData &product) libmesh_override
 For each of the system's quantities of interest q in qoi[qoi_indices], and for a vector of parameters p, the parameter sensitivity Hessian H_ij is defined as H_ij = (d^2 q)/(d p_i d p_j) The Hessian-vector product, for a vector v_k in parameter space, is S_j = H_jk v_k This product is the output of this method, where for each q_i, S_j is stored in sensitivities[i][j]. More...
 
SparseMatrix< Number > & add_matrix (const std::string &mat_name)
 Adds the additional matrix mat_name to this system. More...
 
void remove_matrix (const std::string &mat_name)
 Removes the additional matrix mat_name from this system. More...
 
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
 
virtual void assemble_qoi (const QoISet &qoi_indices=QoISet()) libmesh_override
 Prepares qoi for quantity of interest assembly, then calls user qoi function. More...
 
virtual void assemble_qoi_derivative (const QoISet &qoi_indices=QoISet(), bool include_liftfunc=true, bool apply_constraints=true) libmesh_override
 Prepares adjoint_rhs for quantity of interest derivative assembly, then calls user qoi derivative function. More...
 
void init ()
 Initializes degrees of freedom on the current mesh. More...
 
virtual void reinit_constraints ()
 Reinitializes the constraints for this system. More...
 
bool is_initialized ()
 
virtual void update ()
 Update the local values to reflect the solution on neighboring processors. More...
 
bool is_adjoint_already_solved () const
 Accessor for the adjoint_already_solved boolean. More...
 
void set_adjoint_already_solved (bool setting)
 Setter for the adjoint_already_solved boolean. More...
 
virtual void qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
 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]. More...
 
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
 Projects arbitrary functions onto the current solution. More...
 
void project_solution (FEMFunctionBase< Number > *f, FEMFunctionBase< Gradient > *g=libmesh_nullptr) const
 Projects arbitrary functions onto the current solution. More...
 
void project_solution (Number fptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Gradient gptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), const Parameters &parameters) const
 Projects arbitrary functions onto the current solution. More...
 
void project_vector (NumericVector< Number > &new_vector, FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr, int is_adjoint=-1) const
 Projects arbitrary functions onto a vector of degree of freedom values for the current system. More...
 
void project_vector (NumericVector< Number > &new_vector, FEMFunctionBase< Number > *f, FEMFunctionBase< Gradient > *g=libmesh_nullptr, int is_adjoint=-1) const
 Projects arbitrary functions onto a vector of degree of freedom values for the current system. More...
 
void project_vector (Number fptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Gradient gptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), const Parameters &parameters, NumericVector< Number > &new_vector, int is_adjoint=-1) const
 Projects arbitrary functions onto a vector of degree of freedom values for the current system. More...
 
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)
 Projects arbitrary boundary functions onto a vector of degree of freedom values for the current system. More...
 
void boundary_project_solution (const std::set< boundary_id_type > &b, const std::vector< unsigned int > &variables, Number fptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Gradient gptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), const Parameters &parameters)
 Projects arbitrary boundary functions onto a vector of degree of freedom values for the current system. More...
 
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
 Projects arbitrary boundary functions onto a vector of degree of freedom values for the current system. More...
 
void boundary_project_vector (const std::set< boundary_id_type > &b, const std::vector< unsigned int > &variables, Number fptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), Gradient gptr(const Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name), const Parameters &parameters, NumericVector< Number > &new_vector, int is_adjoint=-1) const
 Projects arbitrary boundary functions onto a vector of degree of freedom values for the current system. More...
 
unsigned int number () const
 
void update_global_solution (std::vector< Number > &global_soln) const
 Fill the input vector global_soln so that it contains the global solution on all processors. More...
 
void update_global_solution (std::vector< Number > &global_soln, const processor_id_type dest_proc) const
 Fill the input vector global_soln so that it contains the global solution on processor dest_proc. More...
 
const MeshBase & get_mesh () const
 
MeshBase & get_mesh ()
 
const DofMap & get_dof_map () const
 
DofMap & get_dof_map ()
 
const EquationSystems & get_equation_systems () const
 
EquationSystems & get_equation_systems ()
 
bool active () const
 
void activate ()
 Activates the system. More...
 
void deactivate ()
 Deactivates the system. More...
 
void set_basic_system_only ()
 Sets the system to be "basic only": i.e. More...
 
vectors_iterator vectors_begin ()
 Beginning of vectors container. More...
 
const_vectors_iterator vectors_begin () const
 Beginning of vectors container. More...
 
vectors_iterator vectors_end ()
 End of vectors container. More...
 
const_vectors_iterator vectors_end () const
 End of vectors container. More...
 
NumericVector< Number > & add_vector (const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
 Adds the additional vector vec_name to this system. More...
 
void remove_vector (const std::string &vec_name)
 Removes the additional vector vec_name from this system. More...
 
bool & project_solution_on_reinit (void)
 Tells the System whether or not to project the solution vector onto new grids when the system is reinitialized. More...
 
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)
 Allows one to set the QoI index controlling whether the vector identified by vec_name represents a solution from the adjoint (qoi_num >= 0) or primal (qoi_num == -1) space. More...
 
int vector_is_adjoint (const std::string &vec_name) const
 
void set_vector_preservation (const std::string &vec_name, bool preserve)
 Allows one to set the boolean controlling whether the vector identified by vec_name should be "preserved": projected to new meshes, saved, etc. More...
 
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)
 Adds the variable var to the list of variables for this system. More...
 
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)
 Adds the variable var to the list of variables for this system. More...
 
unsigned int add_variables (const std::vector< std::string > &vars, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
 Adds the variable var to the list of variables for this system. More...
 
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)
 Adds the variable var to the list of variables for this system. More...
 
const Variable & variable (unsigned int var) const
 Return a constant reference to Variable var. More...
 
const VariableGroup & variable_group (unsigned int vg) const
 Return a constant reference to VariableGroup vg. More...
 
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
 Fills all_variable_numbers with all the variable numbers for the variables that have been added to this system. More...
 
unsigned int variable_scalar_number (const std::string &var, unsigned int component) const
 
unsigned int variable_scalar_number (unsigned int var_num, unsigned int component) const
 
const FEType & variable_type (const unsigned int i) const
 
const FEType & variable_type (const std::string &var) const
 
bool identify_variable_groups () const
 
void identify_variable_groups (const bool)
 Toggle automatic VariableGroup identification. More...
 
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)
 Reads the basic data header for this System. More...
 
void read_legacy_data (Xdr &io, const bool read_additional_data=true)
 Reads additional data, namely vectors, for this System. More...
 
template<typename ValType >
void read_serialized_data (Xdr &io, const bool read_additional_data=true)
 Reads additional data, namely vectors, for this System. More...
 
void read_serialized_data (Xdr &io, const bool read_additional_data=true)
 Non-templated version for backward compatibility. More...
 
template<typename InValType >
std::size_t read_serialized_vectors (Xdr &io, const std::vector< NumericVector< Number > * > &vectors) const
 Read a number of identically distributed vectors. More...
 
std::size_t read_serialized_vectors (Xdr &io, const std::vector< NumericVector< Number > * > &vectors) const
 Non-templated version for backward compatibility. More...
 
template<typename InValType >
void read_parallel_data (Xdr &io, const bool read_additional_data)
 Reads additional data, namely vectors, for this System. More...
 
void read_parallel_data (Xdr &io, const bool read_additional_data)
 Non-templated version for backward compatibility. More...
 
void write_header (Xdr &io, const std::string &version, const bool write_additional_data) const
 Writes the basic data header for this System. More...
 
void write_serialized_data (Xdr &io, const bool write_additional_data=true) const
 Writes additional data, namely vectors, for this System. More...
 
std::size_t write_serialized_vectors (Xdr &io, const std::vector< const NumericVector< Number > * > &vectors) const
 Serialize & write a number of identically distributed vectors. More...
 
void write_parallel_data (Xdr &io, const bool write_additional_data) const
 Writes additional data, namely vectors, for this System. More...
 
std::string get_info () const
 
void attach_init_function (void fptr(EquationSystems &es, const std::string &name))
 Register a user function to use in initializing the system. More...
 
void attach_init_object (Initialization &init)
 Register a user class to use to initialize the system. More...
 
void attach_assemble_function (void fptr(EquationSystems &es, const std::string &name))
 Register a user function to use in assembling the system matrix and RHS. More...
 
void attach_assemble_object (Assembly &assemble)
 Register a user object to use in assembling the system matrix and RHS. More...
 
void attach_constraint_function (void fptr(EquationSystems &es, const std::string &name))
 Register a user function for imposing constraints. More...
 
void attach_constraint_object (Constraint &constrain)
 Register a user object for imposing constraints. More...
 
void attach_QOI_function (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))
 Register a user function for evaluating the quantities of interest, whose values should be placed in System::qoi. More...
 
void attach_QOI_object (QOI &qoi)
 Register a user object for evaluating the quantities of interest, whose values should be placed in System::qoi. More...
 
void attach_QOI_derivative (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices, bool include_liftfunc, bool apply_constraints))
 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. More...
 
void attach_QOI_derivative_object (QOIDerivative &qoi_derivative)
 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. More...
 
virtual void user_initialization ()
 Calls user's attached initialization function, or is overridden by the user in derived classes. More...
 
virtual void user_assembly ()
 Calls user's attached assembly function, or is overridden by the user in derived classes. More...
 
virtual void user_constrain ()
 Calls user's attached constraint function, or is overridden by the user in derived classes. More...
 
virtual void user_QOI (const QoISet &qoi_indices)
 Calls user's attached quantity of interest function, or is overridden by the user in derived classes. More...
 
virtual void user_QOI_derivative (const QoISet &qoi_indices=QoISet(), bool include_liftfunc=true, bool apply_constraints=true)
 Calls user's attached quantity of interest derivative function, or is overridden by the user in derived classes. More...
 
virtual void re_update ()
 Re-update the local values when the mesh has changed. More...
 
virtual void restrict_vectors ()
 Restrict vectors after the mesh has coarsened. More...
 
virtual void prolong_vectors ()
 Prolong vectors after the mesh has refined. More...
 
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
 Calls the version of point_value() which takes a reference. More...
 
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
 Calls the version of point_gradient() which takes a reference. More...
 
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
 Calls the version of point_hessian() which takes a reference. More...
 
void local_dof_indices (const unsigned int var, std::set< dof_id_type > &var_indices) const
 Fills the std::set with the degrees of freedom on the local processor corresponding the the variable number passed in. More...
 
void zero_variable (NumericVector< Number > &v, unsigned int var_num) const
 Zeroes all dofs in v that correspond to variable number var_num. More...
 
bool & hide_output ()
 
const Parallel::Communicator & comm () const
 
processor_id_type n_processors () const
 
processor_id_type processor_id () const
 
void initialize_parameters (const RBParameters &mu_min_in, const RBParameters &mu_max_in, const std::map< std::string, std::vector< Real >> &discrete_parameter_values)
 Initialize the parameter ranges and set current_parameters. More...
 
void initialize_parameters (const RBParametrized &rb_parametrized)
 Initialize the parameter ranges and set current_parameters. More...
 
unsigned int get_n_params () const
 Get the number of parameters. More...
 
unsigned int get_n_continuous_params () const
 Get the number of continuous parameters. More...
 
unsigned int get_n_discrete_params () const
 Get the number of discrete parameters. More...
 
std::set< std::string > get_parameter_names () const
 Get a set that stores the parameter names. More...
 
const RBParameters & get_parameters () const
 Get the current parameters. More...
 
void set_parameters (const RBParameters &params)
 Set the current parameters to params. More...
 
const RBParameters & get_parameters_min () const
 Get an RBParameters object that specifies the minimum allowable value for each parameter. More...
 
const RBParameters & get_parameters_max () const
 Get an RBParameters object that specifies the maximum allowable value for each parameter. More...
 
Real get_parameter_min (const std::string &param_name) const
 Get minimum allowable value of parameter param_name. More...
 
Real get_parameter_max (const std::string &param_name) const
 Get maximum allowable value of parameter param_name. More...
 
void print_parameters () const
 Print the current parameters. More...
 
void write_parameter_data_to_files (const std::string &continuous_param_file_name, const std::string &discrete_param_file_name, const bool write_binary_data)
 Write out the parameter ranges to files. More...
 
void read_parameter_data_from_files (const std::string &continuous_param_file_name, const std::string &discrete_param_file_name, const bool read_binary_data)
 Read in the parameter ranges from files. More...
 
bool is_discrete_parameter (const std::string &mu_name) const
 Is parameter mu_name discrete? More...
 
const std::map< std::string, std::vector< Real > > & get_discrete_parameter_values () const
 Get a const reference to the discrete parameter values. More...
 
void print_discrete_parameter_values () const
 Print out all the discrete parameter values. More...
 

Static Public Member Functions

static void print_info (std::ostream &out=libMesh::out)
 Prints the reference information, by default to libMesh::out. More...
 
static void print_info (std::ostream &out=libMesh::out)
 Prints the reference information, by default to libMesh::out. More...
 
static UniquePtr< DirichletBoundary > build_zero_dirichlet_boundary_object ()
 It's helpful to be able to generate a DirichletBoundary that stores a ZeroFunction in order to impose Dirichlet boundary conditions. More...
 
static std::string get_info ()
 Gets a string containing the reference information. More...
 
static std::string get_info ()
 Gets a string containing the reference information. More...
 
static unsigned int n_objects ()
 Prints the number of outstanding (created, but not yet destroyed) objects. More...
 
static unsigned int n_objects ()
 Prints the number of outstanding (created, but not yet destroyed) objects. More...
 
static void enable_print_counter_info ()
 Methods to enable/disable the reference counter output from print_info() More...
 
static void enable_print_counter_info ()
 Methods to enable/disable the reference counter output from print_info() More...
 
static void disable_print_counter_info ()
 
static void disable_print_counter_info ()
 
static Real get_closest_value (Real value, const std::vector< Real > &list_of_values)
 

Public Attributes

unsigned int u_var
 Variable numbers. More...
 
unsigned int v_var
 
unsigned int w_var
 
ElasticityAssemblyExpansion elasticity_assembly_expansion
 The object that stores the "assembly" expansion of the parameter dependent PDE. More...
 
InnerProductAssembly ip_assembly
 Object to assemble the inner product matrix. More...
 
UniquePtr< DirichletBoundarydirichlet_bc
 The object that defines which degrees of freedom are on a Dirichlet boundary. More...
 
std::vector< Real > training_error_bounds
 Vector storing the values of the error bound for each parameter in the training set — the parameter giving the largest error bound is chosen for the next snapshot in the Greedy basis training. More...
 
UniquePtr< LinearSolver< Number > > inner_product_solver
 We store an extra linear solver object which we can optionally use for solving all systems in which the system matrix is set to inner_product_matrix. More...
 
LinearSolver< Number > * extra_linear_solver
 Also, we store a pointer to an extra linear solver. More...
 
UniquePtr< SparseMatrix< Number > > inner_product_matrix
 The inner product matrix. More...
 
std::vector< Number > truth_outputs
 Vector storing the truth output values from the most recent truth solve. More...
 
std::vector< std::vector< Number > > output_dual_innerprods
 The vector storing the dual norm inner product terms for each output. More...
 
std::vector< NumericVector< Number > * > Fq_representor
 Vector storing the residual representors associated with the right-hand side. More...
 
std::vector< Number > Fq_representor_innerprods
 Vectors storing the residual representor inner products to be used in computing the residuals online. More...
 
bool exit_on_repeated_greedy_parameters
 Boolean flag to indicate whether we exit the greedy if we select the same parameters twice in a row. More...
 
bool impose_internal_fluxes
 Boolean flag to indicate whether we impose "fluxes" (i.e. More...
 
bool compute_RB_inner_product
 Boolean flag to indicate whether we compute the RB_inner_product_matrix. More...
 
bool store_non_dirichlet_operators
 Boolean flag to indicate whether we store a second copy of each affine operator and vector which does not have Dirichlet bcs enforced. More...
 
bool use_empty_rb_solve_in_greedy
 A boolean flag to indicate whether or not we initialize the Greedy algorithm by performing rb_solves on the training set with an "empty" (i.e. More...
 
bool Fq_representor_innerprods_computed
 A boolean flag to indicate whether or not the Fq representor norms have already been computed — used to make sure that we don't recompute them unnecessarily. More...
 
UniquePtr< LinearSolver< Number > > linear_solver
 The LinearSolver defines the interface used to solve the linear_implicit system. More...
 
SparseMatrix< Number > * matrix
 The system matrix. More...
 
bool zero_out_matrix_and_rhs
 By default, the system will zero out the matrix and the right hand side. More...
 
NumericVector< Number > * rhs
 The system matrix. More...
 
bool assemble_before_solve
 Flag which tells the system to whether or not to call the user assembly function during each call to solve(). More...
 
bool use_fixed_solution
 A boolean to be set to true by systems using elem_fixed_solution, for optional use by e.g. More...
 
int extra_quadrature_order
 A member int that can be employed to indicate increased or reduced quadrature order. More...
 
UniquePtr< NumericVector< Number > > solution
 Data structure to hold solution values. More...
 
UniquePtr< NumericVector< Number > > current_local_solution
 All the values I need to compute my contribution to the simulation at hand. More...
 
Real time
 For time-dependent problems, this is the time t at the beginning of the current timestep. More...
 
std::vector< Number > qoi
 Values of the quantities of interest. More...
 
bool verbose_mode
 Public boolean to toggle verbose mode. More...
 

Protected Types

typedef std::map< std::string, std::pair< unsigned int, unsigned int > > Counts
 Data structure to log the information. More...
 
typedef std::map< std::string, std::pair< unsigned int, unsigned int > > Counts
 Data structure to log the information. More...
 

Protected Member Functions

virtual void allocate_data_structures ()
 Helper function that actually allocates all the data structures required by this class. More...
 
virtual void truth_assembly ()
 Assemble the truth matrix and right-hand side for current_parameters. More...
 
virtual UniquePtr< DGFEMContext > build_context ()
 Builds a DGFEMContext object with enough information to do evaluations on each element. More...
 
virtual SparseMatrix< Number > & get_matrix_for_output_dual_solves ()
 Return the matrix for the output residual dual norm solves. More...
 
virtual bool greedy_termination_test (Real abs_greedy_error, Real initial_greedy_error, int count)
 Function that indicates when to terminate the Greedy basis training. More...
 
void update_greedy_param_list ()
 Update the list of Greedily chosen parameters with current_parameters. More...
 
void add_scaled_matrix_and_vector (Number scalar, ElemAssembly *elem_assembly, SparseMatrix< Number > *input_matrix, NumericVector< Number > *input_vector, bool symmetrize=false, bool apply_dof_constraints=true)
 This function loops over the mesh and applies the specified interior and/or boundary assembly routines, then adds the scaled result to input_matrix and/or input_vector. More...
 
virtual void set_context_solution_vec (NumericVector< Number > &vec)
 Set current_local_solution = vec so that we can access vec from FEMContext during assembly. More...
 
virtual void assemble_misc_matrices ()
 Assemble and store all the inner-product matrix, the constraint matrix (for constrained problems) and the mass matrix (for time-dependent problems). More...
 
virtual void assemble_all_affine_operators ()
 Assemble and store all Q_a affine operators as well as the inner-product matrix. More...
 
virtual void assemble_all_affine_vectors ()
 Assemble and store the affine RHS vectors. More...
 
virtual void assemble_all_output_vectors ()
 Assemble and store the output vectors. More...
 
virtual void compute_output_dual_innerprods ()
 Compute and store the dual norm of each output functional. More...
 
virtual void compute_Fq_representor_innerprods (bool compute_inner_products=true)
 Compute the terms that are combined `online' to determine the dual norm of the residual. More...
 
virtual void enrich_RB_space ()
 Add a new basis function to the RB space. More...
 
virtual void update_system ()
 Update the system after enriching the RB space; this calls a series of functions to update the system properly. More...
 
virtual Real get_RB_error_bound ()
 
virtual void update_RB_system_matrices ()
 Compute the reduced basis matrices for the current basis. More...
 
virtual void update_residual_terms (bool compute_inner_products=true)
 Compute the terms that are combined `online' to determine the dual norm of the residual. More...
 
bool get_convergence_assertion_flag () const
 Getter for the flag determining if convergence should be checked after each solve. More...
 
void check_convergence (LinearSolver< Number > &input_solver)
 Check if the linear solver reports convergence. More...
 
RBParameters get_params_from_training_set (unsigned int index)
 Return the RBParameters in index index of training set. More...
 
void set_params_from_training_set (unsigned int index)
 Set parameters to the RBParameters stored in index index of the training set. More...
 
virtual void set_params_from_training_set_and_broadcast (unsigned int index)
 Load the specified training parameter and then broadcast to all processors. More...
 
virtual void init_matrices ()
 Initializes the matrices associated with this system. More...
 
void project_vector (NumericVector< Number > &, int is_adjoint=-1) const
 Projects the vector defined on the old mesh onto the new mesh. More...
 
void project_vector (const NumericVector< Number > &, NumericVector< Number > &, int is_adjoint=-1) const
 Projects the vector defined on the old mesh onto the new mesh. More...
 
void increment_constructor_count (const std::string &name)
 Increments the construction counter. More...
 
void increment_constructor_count (const std::string &name)
 Increments the construction counter. More...
 
void increment_destructor_count (const std::string &name)
 Increments the destruction counter. More...
 
void increment_destructor_count (const std::string &name)
 Increments the destruction counter. More...
 

Static Protected Member Functions

static void get_global_max_error_pair (const Parallel::Communicator &communicator, std::pair< numeric_index_type, Real > &error_pair)
 Static function to return the error pair (index,error) that is corresponds to the largest error on all processors. More...
 
static void generate_training_parameters_random (const Parallel::Communicator &communicator, std::map< std::string, bool > log_param_scale, std::map< std::string, NumericVector< Number > * > &training_parameters_in, unsigned int n_training_samples_in, const RBParameters &min_parameters, const RBParameters &max_parameters, int training_parameters_random_seed=-1, bool serial_training_set=false)
 Static helper function for generating a randomized set of parameters. More...
 
static void generate_training_parameters_deterministic (const Parallel::Communicator &communicator, std::map< std::string, bool > log_param_scale, std::map< std::string, NumericVector< Number > * > &training_parameters_in, unsigned int n_training_samples_in, const RBParameters &min_parameters, const RBParameters &max_parameters, bool serial_training_set=false)
 Static helper function for generating a deterministic set of parameters. More...
 

Protected Attributes

unsigned int Nmax
 Maximum number of reduced basis functions we are willing to use. More...
 
unsigned int delta_N
 The number of basis functions that we add at each greedy step. More...
 
bool quiet_mode
 Flag to indicate whether we print out extra information during the Offline stage. More...
 
bool output_dual_innerprods_computed
 A boolean flag to indicate whether or not the output dual norms have already been computed — used to make sure that we don't recompute them unnecessarily. More...
 
bool assert_convergence
 A boolean flag to indicate whether to check for proper convergence after each solve. More...
 
bool serial_training_set
 This boolean flag indicates whether or not the training set should be the same on all processors. More...
 
UniquePtr< NumericVector< Number > > inner_product_storage_vector
 We keep an extra temporary vector that is useful for performing inner products (avoids unnecessary memory allocation/deallocation). More...
 
unsigned int _n_linear_iterations
 The number of linear iterations required to solve the linear system Ax=b. More...
 
Real _final_linear_residual
 The final residual for the linear system Ax=b. More...
 
ShellMatrix< Number > * _shell_matrix
 User supplies shell matrix or NULL if no shell matrix is used. More...
 
const SystemSubset * _subset
 The current subset on which to solve (or NULL if none). More...
 
SubsetSolveMode _subset_solve_mode
 If restrict-solve-to-subset mode is active, this member decides what happens with the dofs outside the subset. More...
 
const Parallel::Communicator & _communicator
 

Static Protected Attributes

static Counts _counts
 Actually holds the data. More...
 
static Counts _counts
 Actually holds the data. More...
 
static Threads::atomic< unsigned int_n_objects
 The number of objects. More...
 
static Threads::atomic< unsigned int_n_objects
 The number of objects. More...
 
static Threads::spin_mutex _mutex
 Mutual exclusion object to enable thread-safe reference counting. More...
 
static Threads::spin_mutex _mutex
 Mutual exclusion object to enable thread-safe reference counting. More...
 
static bool _enable_print_counter = true
 Flag to control whether reference count information is printed when print_info is called. More...
 
static bool _enable_print_counter = true
 Flag to control whether reference count information is printed when print_info is called. More...
 

Detailed Description

Definition at line 45 of file rb_classes.h.

Member Typedef Documentation

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

Definition at line 283 of file implicit_system.h.

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

Definition at line 749 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 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 std::map<std::string, SparseMatrix<Number> *>::iterator libMesh::ImplicitSystem::matrices_iterator
inherited

Matrix iterator typedefs.

Definition at line 282 of file implicit_system.h.

The type of the parent.

Definition at line 70 of file rb_classes.h.

The type of system.

Definition at line 65 of file rb_classes.h.

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

Vector iterator typedefs.

Definition at line 748 of file system.h.

Constructor & Destructor Documentation

ElasticityRBConstruction::ElasticityRBConstruction ( EquationSystems es,
const std::string &  name_in,
const unsigned int  number_in 
)

Definition at line 49 of file rb_classes.h.

51  :
52  Parent(es, name_in, number_in),
54  ip_assembly(*this)
55  {}
InnerProductAssembly ip_assembly
Object to assemble the inner product matrix.
Definition: rb_classes.h:133
RBConstruction Parent
The type of the parent.
Definition: rb_classes.h:70
ElasticityAssemblyExpansion elasticity_assembly_expansion
The object that stores the "assembly" expansion of the parameter dependent PDE.
Definition: rb_classes.h:128
virtual ElasticityRBConstruction::~ElasticityRBConstruction ( )
virtual

Destructor.

Definition at line 60 of file rb_classes.h.

60 {}

Member Function Documentation

void libMesh::System::activate ( )
inherited

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 _active
Flag stating if the system is active or not.
Definition: system.h:1908
bool libMesh::System::active ( ) const
inherited
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 }
bool _active
Flag stating if the system is active or not.
Definition: system.h:1908
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 1041 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().

1042 {
1043  std::ostringstream adjoint_rhs_name;
1044  adjoint_rhs_name << "adjoint_rhs" << i;
1045 
1046  return this->add_vector(adjoint_rhs_name.str(), false);
1047 }
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:681
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 977 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().

978 {
979  std::ostringstream adjoint_name;
980  adjoint_name << "adjoint_solution" << i;
981 
982  NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
983  this->set_vector_as_adjoint(adjoint_name.str(), i);
984  return returnval;
985 }
void set_vector_as_adjoint(const std::string &vec_name, int qoi_num)
Allows one to set the QoI index controlling whether the vector identified by vec_name represents a so...
Definition: system.C:905
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:681
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(), main(), 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 UniquePtr< SparseMatrix< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a SparseMatrix<T> using the linear solver package specified by solver_package.
std::map< std::string, SparseMatrix< Number > * > _matrices
Some systems need an arbitrary number of matrices.
bool have_matrix(const std::string &mat_name) const
bool _can_add_matrices
true when additional matrices may still be added, false otherwise.
const Parallel::Communicator & comm() const
void libMesh::RBConstruction::add_scaled_Aq ( Number  scalar,
unsigned int  q_a,
SparseMatrix< Number > *  input_matrix,
bool  symmetrize 
)
inherited

Add the scaled q^th affine matrix to input_matrix.

If symmetrize==true, then we symmetrize Aq before adding it.

Definition at line 918 of file rb_construction.C.

References libMesh::SparseMatrix< T >::add(), libMesh::RBConstruction::add_scaled_matrix_and_vector(), libMesh::SparseMatrix< T >::close(), libMesh::RBAssemblyExpansion::get_A_assembly(), libMesh::RBConstruction::get_Aq(), libMesh::RBConstruction::get_rb_theta_expansion(), libmesh_nullptr, and libMesh::RBConstruction::rb_assembly_expansion.

Referenced by libMesh::RBSCMConstruction::add_scaled_symm_Aq(), and libMesh::RBConstruction::is_quiet().

922 {
923  LOG_SCOPE("add_scaled_Aq()", "RBConstruction");
924 
925  if (q_a >= get_rb_theta_expansion().get_n_A_terms())
926  libmesh_error_msg("Error: We must have q < Q_a in add_scaled_Aq.");
927 
928  if (!symmetrize)
929  {
930  input_matrix->add(scalar, *get_Aq(q_a));
931  input_matrix->close();
932  }
933  else
934  {
937  input_matrix,
939  symmetrize);
940  }
941 }
SparseMatrix< Number > * get_Aq(unsigned int q)
Get a pointer to Aq.
const class libmesh_nullptr_t libmesh_nullptr
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value)=0
Add value to the element (i,j).
RBThetaExpansion & get_rb_theta_expansion()
Get a reference to the RBThetaExpansion object that that belongs to rb_eval.
void add_scaled_matrix_and_vector(Number scalar, ElemAssembly *elem_assembly, SparseMatrix< Number > *input_matrix, NumericVector< Number > *input_vector, bool symmetrize=false, bool apply_dof_constraints=true)
This function loops over the mesh and applies the specified interior and/or boundary assembly routine...
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
ElemAssembly & get_A_assembly(unsigned int q)
Return a reference to the specified A_assembly object.
RBAssemblyExpansion * rb_assembly_expansion
This member holds the (parameter independent) assembly functors that define the "affine expansion" of...
void libMesh::RBConstruction::add_scaled_matrix_and_vector ( Number  scalar,
ElemAssembly elem_assembly,
SparseMatrix< Number > *  input_matrix,
NumericVector< Number > *  input_vector,
bool  symmetrize = false,
bool  apply_dof_constraints = true 
)
protectedinherited

This function loops over the mesh and applies the specified interior and/or boundary assembly routines, then adds the scaled result to input_matrix and/or input_vector.

If symmetrize==true then we assemble the symmetric part of the matrix, 0.5*(A + A^T)

Definition at line 604 of file rb_construction.C.

References libMesh::DofMap::_dof_coupling, libMesh::MeshBase::active_local_element_ptr_range(), libMesh::NumericVector< T >::add(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::FEAbstract::attach_quadrature_rule(), libMesh::ConstCouplingRow::begin(), libMesh::ElemAssembly::boundary_assembly(), libMesh::RBConstruction::build_context(), libMesh::BoundaryInfo::build_node_list(), libMesh::NumericVector< T >::close(), libMesh::SparseMatrix< T >::close(), libMesh::ParallelObject::comm(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::FEType::default_quadrature_rule(), libMesh::DGFEMContext::dg_terms_are_active(), libMesh::DenseMatrix< T >::el(), libMesh::FEMContext::elem_fe_reinit(), end, libMesh::ConstCouplingRow::end(), libMesh::err, libMesh::MeshBase::get_boundary_info(), libMesh::DiffContext::get_dof_indices(), libMesh::System::get_dof_map(), libMesh::FEMContext::get_elem(), libMesh::DGFEMContext::get_elem_elem_jacobian(), libMesh::DiffContext::get_elem_jacobian(), libMesh::DGFEMContext::get_elem_neighbor_jacobian(), libMesh::DiffContext::get_elem_residual(), libMesh::FEMContext::get_element_fe(), libMesh::FEAbstract::get_fe_type(), libMesh::System::get_mesh(), libMesh::DGFEMContext::get_neighbor_dof_indices(), libMesh::DGFEMContext::get_neighbor_elem_jacobian(), libMesh::DGFEMContext::get_neighbor_neighbor_jacobian(), libMesh::ElemAssembly::get_nodal_rhs_values(), libMesh::ElemAssembly::get_nodal_values(), libMesh::FEMContext::get_side(), libMesh::DenseMatrix< T >::get_transpose(), libMesh::DofObject::id(), libMesh::RBConstruction::impose_internal_fluxes, libMesh::RBConstruction::init_context(), libMesh::ElemAssembly::interior_assembly(), libMesh::ElemAssembly::is_nodal_rhs_values_overriden, libmesh_nullptr, libMesh::DenseMatrixBase< T >::m(), mesh, libMesh::DenseMatrixBase< T >::n(), libMesh::BoundaryInfo::n_nodeset_conds(), libMesh::Elem::n_sides(), libMesh::System::n_vars(), libMesh::Elem::neighbor_ptr(), libMesh::MeshBase::node_ref(), libMesh::FEMContext::pre_fe_reinit(), libMesh::DofObject::processor_id(), libMesh::Parallel::Communicator::rank(), libMesh::FEMContext::side, libMesh::DGFEMContext::side_fe_reinit(), libMesh::TRI3SUBDIVISION, and value.

Referenced by libMesh::RBConstruction::add_scaled_Aq(), libMesh::RBConstruction::assemble_all_output_vectors(), libMesh::RBConstruction::assemble_Aq_matrix(), libMesh::RBConstruction::assemble_Fq_vector(), libMesh::RBConstruction::assemble_inner_product_matrix(), libMesh::TransientRBConstruction::assemble_L2_matrix(), and libMesh::TransientRBConstruction::assemble_Mq_matrix().

610 {
611  LOG_SCOPE("add_scaled_matrix_and_vector()", "RBConstruction");
612 
613  bool assemble_matrix = (input_matrix != libmesh_nullptr);
614  bool assemble_vector = (input_vector != libmesh_nullptr);
615 
616  if (!assemble_matrix && !assemble_vector)
617  return;
618 
619  const MeshBase & mesh = this->get_mesh();
620 
621  // First add any node-based terms (e.g. point loads)
622  // We only enter this loop if we have at least one
623  // nodeset, since we use nodesets to indicate
624  // where to impose the node-based terms.
625  if (mesh.get_boundary_info().n_nodeset_conds() > 0)
626  {
627  std::vector<numeric_index_type> node_id_list;
628  std::vector<boundary_id_type> bc_id_list;
629 
630  // Get the list of nodes with boundary IDs
631  mesh.get_boundary_info().build_node_list(node_id_list, bc_id_list);
632 
633  for (std::size_t i=0; i<node_id_list.size(); i++)
634  {
635  const Node & node = mesh.node_ref(node_id_list[i]);
636 
637  // If node is on this processor, then all dofs on node are too
638  // so we can do the add below safely
639  if (node.processor_id() == this->comm().rank())
640  {
641  // Get the values to add to the rhs vector
642  std::vector<dof_id_type> nodal_dof_indices;
643  DenseMatrix<Number> nodal_matrix;
644  DenseVector<Number> nodal_rhs;
645  elem_assembly->get_nodal_values(
646  nodal_dof_indices, nodal_matrix, nodal_rhs, *this, node);
647 
648  if(!nodal_dof_indices.empty())
649  {
650  if(assemble_vector)
651  {
652  input_vector->add_vector(nodal_rhs, nodal_dof_indices);
653  }
654 
655  if(assemble_matrix)
656  {
657  input_matrix->add_matrix(nodal_matrix, nodal_dof_indices);
658  }
659  }
660 #ifdef LIBMESH_ENABLE_DEPRECATED
661  else if(elem_assembly->is_nodal_rhs_values_overriden)
662  {
663  // This is the "old" implementation, to be deprecated soon.
664  // Only enter this if ElemAssembly::get_nodal_rhs_values has
665  // a non-default implementation.
666  libmesh_deprecated();
667 
668  std::map<numeric_index_type, Number> rhs_values;
669  elem_assembly->get_nodal_rhs_values(rhs_values, *this, node);
670 
671  std::map<numeric_index_type, Number>::const_iterator it =
672  rhs_values.begin();
673  const std::map<numeric_index_type, Number>::const_iterator it_end =
674  rhs_values.end();
675  for ( ; it != it_end; ++it)
676  {
677  numeric_index_type dof_index = it->first;
678  Number value = it->second;
679  input_vector->add( dof_index, value);
680  }
681  }
682 #endif
683  }
684  }
685  }
686 
687  UniquePtr<DGFEMContext> c = this->build_context();
688  DGFEMContext & context = cast_ref<DGFEMContext &>(*c);
689 
690  this->init_context(context);
691 
692  for (const auto & elem : mesh.active_local_element_ptr_range())
693  {
694  // Subdivision elements need special care:
695  // - skip ghost elements
696  // - init special quadrature rule
697  UniquePtr<QBase> qrule;
698  if (elem->type() == TRI3SUBDIVISION)
699  {
700  const Tri3Subdivision * gh_elem = static_cast<const Tri3Subdivision *> (elem);
701  if (gh_elem->is_ghost())
702  continue ;
703  // A Gauss quadrature rule for numerical integration.
704  // For subdivision shell elements, a single Gauss point per
705  // element is sufficient, hence we use extraorder = 0.
706  const int extraorder = 0;
707  FEBase * elem_fe = libmesh_nullptr;
708  context.get_element_fe( 0, elem_fe );
709 
710  qrule = elem_fe->get_fe_type().default_quadrature_rule (2, extraorder);
711 
712  // Tell the finite element object to use our quadrature rule.
713  elem_fe->attach_quadrature_rule (qrule.get());
714  }
715 
716  context.pre_fe_reinit(*this, elem);
717  context.elem_fe_reinit();
718  elem_assembly->interior_assembly(context);
719 
720  for (context.side = 0;
721  context.side != context.get_elem().n_sides();
722  ++context.side )
723  {
724  // May not need to apply fluxes on non-boundary elements
725  if ((context.get_elem().neighbor_ptr(context.get_side()) != libmesh_nullptr) && !impose_internal_fluxes)
726  continue;
727 
728  bool reinit_succeeded = false;
729 
730  // Exceptions can be thrown in side_fe_reinit, e.g. due to a zero or negative
731  // Jacobian on an element side. This is sometimes intentional, e.g. some people
732  // use "degenerate HEX" elements when generating meshes. In order to support this
733  // case we just print a message and skip the side assembly on this element
734  // if an exception occurs during side_fe_reinit.
735  try
736  {
737  context.side_fe_reinit();
738  reinit_succeeded = true;
739  }
740  catch(std::exception& e)
741  {
742  libMesh::err << std::endl << "Error detected when computing element data on side "
743  << static_cast<int>(context.side) << " of element "
744  << context.get_elem().id() << std::endl;
745  libMesh::err << "Skipping assembly on this element side" << std::endl << std::endl;
746  }
747 
748  if(reinit_succeeded)
749  {
750  elem_assembly->boundary_assembly(context);
751  }
752 
753  if (context.dg_terms_are_active())
754  {
755  input_matrix->add_matrix (context.get_elem_elem_jacobian(),
756  context.get_dof_indices(),
757  context.get_dof_indices());
758 
759  input_matrix->add_matrix (context.get_elem_neighbor_jacobian(),
760  context.get_dof_indices(),
761  context.get_neighbor_dof_indices());
762 
763  input_matrix->add_matrix (context.get_neighbor_elem_jacobian(),
764  context.get_neighbor_dof_indices(),
765  context.get_dof_indices());
766 
767  input_matrix->add_matrix (context.get_neighbor_neighbor_jacobian(),
768  context.get_neighbor_dof_indices(),
769  context.get_neighbor_dof_indices());
770  }
771  }
772 
773  // Need to symmetrize before imposing
774  // periodic constraints
775  if (assemble_matrix && symmetrize)
776  {
777  DenseMatrix<Number> Ke_transpose;
778  context.get_elem_jacobian().get_transpose(Ke_transpose);
779  context.get_elem_jacobian() += Ke_transpose;
780  context.get_elem_jacobian() *= 0.5;
781  }
782 
783  if (apply_dof_constraints)
784  {
785  // Apply constraints, e.g. Dirichlet and periodic constraints
787  (context.get_elem_jacobian(), context.get_elem_residual(), context.get_dof_indices() );
788  }
789 
790  // Scale and add to global matrix and/or vector
791  context.get_elem_jacobian() *= scalar;
792  context.get_elem_residual() *= scalar;
793 
794  if (assemble_matrix)
795  {
796 
797  CouplingMatrix * coupling_matrix = get_dof_map()._dof_coupling;
798  if (!coupling_matrix)
799  {
800  // If we haven't defined a _dof_coupling matrix then just add
801  // the whole matrix
802  input_matrix->add_matrix (context.get_elem_jacobian(),
803  context.get_dof_indices() );
804  }
805  else
806  {
807  // Otherwise we should only add the relevant submatrices
808  for (unsigned int var1=0; var1<n_vars(); var1++)
809  {
810  ConstCouplingRow ccr(var1, *coupling_matrix);
813  ccr.begin(); it != end; ++it)
814  {
815  unsigned int var2 = *it;
816 
817  unsigned int sub_m = context.get_elem_jacobian( var1, var2 ).m();
818  unsigned int sub_n = context.get_elem_jacobian( var1, var2 ).n();
819  DenseMatrix<Number> sub_jac(sub_m, sub_n);
820  for (unsigned int row=0; row<sub_m; row++)
821  for (unsigned int col=0; col<sub_n; col++)
822  {
823  sub_jac(row,col) = context.get_elem_jacobian( var1, var2 ).el(row,col);
824  }
825  input_matrix->add_matrix (sub_jac,
826  context.get_dof_indices(var1),
827  context.get_dof_indices(var2) );
828  }
829  }
830  }
831 
832  }
833 
834  if (assemble_vector)
835  input_vector->add_vector (context.get_elem_residual(),
836  context.get_dof_indices() );
837  }
838 
839  if (assemble_matrix)
840  input_matrix->close();
841  if (assemble_vector)
842  input_vector->close();
843 }
OStreamProxy err
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
virtual void init_context(FEMContext &)
Initialize the FEMContext prior to performing an element loop.
bool impose_internal_fluxes
Boolean flag to indicate whether we impose "fluxes" (i.e.
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
const MeshBase & get_mesh() const
Definition: system.h:2014
const DofMap & get_dof_map() const
Definition: system.h:2030
dof_id_type numeric_index_type
Definition: id_types.h:92
void constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix and vector.
Definition: dof_map.h:1806
CouplingMatrix * _dof_coupling
Degree of freedom coupling.
Definition: dof_map.h:1240
FEGenericBase< Real > FEBase
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
const Parallel::Communicator & comm() const
static const bool value
Definition: xdr_io.C:108
unsigned int rank() const
Definition: parallel.h:724
unsigned int n_vars() const
Definition: system.h:2086
virtual void add(const numeric_index_type i, const T value)=0
Adds value to each entry of the vector.
ConstCouplingRowConstIterator const_iterator
virtual UniquePtr< DGFEMContext > build_context()
Builds a DGFEMContext object with enough information to do evaluations on each element.
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 1071 of file system.C.

References libMesh::System::add_vector().

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

1072 {
1073  std::ostringstream sensitivity_rhs_name;
1074  sensitivity_rhs_name << "sensitivity_rhs" << i;
1075 
1076  return this->add_vector(sensitivity_rhs_name.str(), false);
1077 }
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:681
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 926 of file system.C.

References libMesh::System::add_vector().

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

927 {
928  std::ostringstream sensitivity_name;
929  sensitivity_name << "sensitivity_solution" << i;
930 
931  return this->add_vector(sensitivity_name.str());
932 }
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:681
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 1101 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::libmesh_assert(), 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(), assemble_and_solve(), SolidSystem::init_data(), CurlCurlSystem::init_data(), HeatSystem::init_data(), main(), libMesh::ErrorVector::plot_error(), libMesh::System::project_solution_on_reinit(), ParsedFEMFunctionTest::setUp(), FETest< order, family, elem_type >::setUp(), SlitMeshRefinedSystemTest::setUp(), MeshfunctionDFEM::test_mesh_function_dfem(), MeshfunctionDFEM::test_mesh_function_dfem_grad(), DefaultCouplingTest::testCoupling(), PointNeighborCouplingTest::testCoupling(), DofMapTest::testDofOwner(), EquationSystemsTest::testPostInitAddElem(), EquationSystemsTest::testPostInitAddRealSystem(), SystemsTest::testProjectCube(), SystemsTest::testProjectCubeWithMeshFunction(), SystemsTest::testProjectLine(), SystemsTest::testProjectSquare(), EquationSystemsTest::testRefineThenReinitPreserveFlags(), BoundaryInfoTest::testShellFaceConstraints(), and WriteVecAndScalar::testWrite().

1104 {
1105  libmesh_assert(!this->is_initialized());
1106 
1107  // Make sure the variable isn't there already
1108  // or if it is, that it's the type we want
1109  for (unsigned int v=0; v<this->n_vars(); v++)
1110  if (this->variable_name(v) == var)
1111  {
1112  if (this->variable_type(v) == type)
1113  return _variables[v].number();
1114 
1115  libmesh_error_msg("ERROR: incompatible variable " << var << " has already been added for this system!");
1116  }
1117 
1118  // Optimize for VariableGroups here - if the user is adding multiple
1119  // variables of the same FEType and subdomain restriction, catch
1120  // that here and add them as members of the same VariableGroup.
1121  //
1122  // start by setting this flag to whatever the user has requested
1123  // and then consider the conditions which should negate it.
1124  bool should_be_in_vg = this->identify_variable_groups();
1125 
1126  // No variable groups, nothing to add to
1127  if (!this->n_variable_groups())
1128  should_be_in_vg = false;
1129 
1130  else
1131  {
1132  VariableGroup & vg(_variable_groups.back());
1133 
1134  // get a pointer to their subdomain restriction, if any.
1135  const std::set<subdomain_id_type> * const
1136  their_active_subdomains (vg.implicitly_active() ?
1137  libmesh_nullptr : &vg.active_subdomains());
1138 
1139  // Different types?
1140  if (vg.type() != type)
1141  should_be_in_vg = false;
1142 
1143  // they are restricted, we aren't?
1144  if (their_active_subdomains && !active_subdomains)
1145  should_be_in_vg = false;
1146 
1147  // they aren't restricted, we are?
1148  if (!their_active_subdomains && active_subdomains)
1149  should_be_in_vg = false;
1150 
1151  if (their_active_subdomains && active_subdomains)
1152  // restricted to different sets?
1153  if (*their_active_subdomains != *active_subdomains)
1154  should_be_in_vg = false;
1155 
1156  // OK, after all that, append the variable to the vg if none of the conditions
1157  // were violated
1158  if (should_be_in_vg)
1159  {
1160  const unsigned short curr_n_vars = cast_int<unsigned short>
1161  (this->n_vars());
1162 
1163  vg.append (var);
1164 
1165  _variables.push_back(vg(vg.n_variables()-1));
1166  _variable_numbers[var] = curr_n_vars;
1167  return curr_n_vars;
1168  }
1169  }
1170 
1171  // otherwise, fall back to adding a single variable group
1172  return this->add_variables (std::vector<std::string>(1, var),
1173  type,
1174  active_subdomains);
1175 }
std::map< std::string, unsigned short int > _variable_numbers
The variable numbers corresponding to user-specified names, useful for name-based lookups...
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
The Variable in this System.
Definition: system.h:1892
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
libmesh_assert(j)
unsigned int n_variable_groups() const
Definition: system.h:2094
std::vector< VariableGroup > _variable_groups
The VariableGroup in this System.
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)
Adds the variable var to the list of variables for this system.
Definition: system.C:1191
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 1179 of file system.C.

References libMesh::System::add_variable().

1183 {
1184  return this->add_variable(var,
1185  FEType(order, family),
1186  active_subdomains);
1187 }
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
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 1191 of file system.C.

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

1194 {
1195  libmesh_assert(!this->is_initialized());
1196 
1197  // Make sure the variable isn't there already
1198  // or if it is, that it's the type we want
1199  for (std::size_t ov=0; ov<vars.size(); ov++)
1200  for (unsigned int v=0; v<this->n_vars(); v++)
1201  if (this->variable_name(v) == vars[ov])
1202  {
1203  if (this->variable_type(v) == type)
1204  return _variables[v].number();
1205 
1206  libmesh_error_msg("ERROR: incompatible variable " << vars[ov] << " has already been added for this system!");
1207  }
1208 
1209  const unsigned short curr_n_vars = cast_int<unsigned short>
1210  (this->n_vars());
1211 
1212  const unsigned int next_first_component = this->n_components();
1213 
1214  // Add the variable group to the list
1215  _variable_groups.push_back((active_subdomains == libmesh_nullptr) ?
1216  VariableGroup(this, vars, curr_n_vars,
1217  next_first_component, type) :
1218  VariableGroup(this, vars, curr_n_vars,
1219  next_first_component, type, *active_subdomains));
1220 
1221  const VariableGroup & vg (_variable_groups.back());
1222 
1223  // Add each component of the group individually
1224  for (std::size_t v=0; v<vars.size(); v++)
1225  {
1226  _variables.push_back (vg(v));
1227  _variable_numbers[vars[v]] = cast_int<unsigned short>
1228  (curr_n_vars+v);
1229  }
1230 
1231  libmesh_assert_equal_to ((curr_n_vars+vars.size()), this->n_vars());
1232 
1233  // BSK - Defer this now to System::init_data() so we can detect
1234  // VariableGroups 12/28/2012
1235  // // Add the variable group to the _dof_map
1236  // _dof_map->add_variable_group (vg);
1237 
1238  // Return the number of the new variable
1239  return cast_int<unsigned int>(curr_n_vars+vars.size()-1);
1240 }
std::map< std::string, unsigned short int > _variable_numbers
The variable numbers corresponding to user-specified names, useful for name-based lookups...
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
The Variable in this System.
Definition: system.h:1892
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
libmesh_assert(j)
unsigned int n_components() const
Definition: system.h:2102
std::vector< VariableGroup > _variable_groups
The VariableGroup in this System.
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 1244 of file system.C.

References libMesh::System::add_variables().

1248 {
1249  return this->add_variables(vars,
1250  FEType(order, family),
1251  active_subdomains);
1252 }
unsigned int add_variables(const std::vector< std::string > &vars, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1191
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 681 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(), libMesh::ContinuationSystem::init_data(), main(), libMesh::NewmarkSystem::NewmarkSystem(), libMesh::FrequencySystem::set_frequencies(), libMesh::FrequencySystem::set_frequencies_by_range(), and libMesh::FrequencySystem::set_frequencies_by_steps().

684 {
685  // Return the vector if it is already there.
686  if (this->have_vector(vec_name))
687  return *(_vectors[vec_name]);
688 
689  // Otherwise build the vector
690  NumericVector<Number> * buf = NumericVector<Number>::build(this->comm()).release();
691  _vectors.insert (std::make_pair (vec_name, buf));
692  _vector_projections.insert (std::make_pair (vec_name, projections));
693 
694  _vector_types.insert (std::make_pair (vec_name, type));
695 
696  // Vectors are primal by default
697  _vector_is_adjoint.insert (std::make_pair (vec_name, -1));
698 
699  // Initialize it if necessary
700  if (_is_initialized)
701  {
702  if (type == GHOSTED)
703  {
704 #ifdef LIBMESH_ENABLE_GHOSTED
705  buf->init (this->n_dofs(), this->n_local_dofs(),
706  _dof_map->get_send_list(), false,
707  GHOSTED);
708 #else
709  libmesh_error_msg("Cannot initialize ghosted vectors when they are not enabled.");
710 #endif
711  }
712  else
713  buf->init (this->n_dofs(), this->n_local_dofs(), false, type);
714  }
715 
716  return *buf;
717 }
std::map< std::string, ParallelType > _vector_types
Holds the type of a vector.
Definition: system.h:1933
bool _is_initialized
true when additional vectors and variables do not require immediate initialization, false otherwise.
Definition: system.h:1952
UniquePtr< DofMap > _dof_map
Data structure describing the relationship between nodes, variables, etc...
Definition: system.h:1865
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
bool have_vector(const std::string &vec_name) const
Definition: system.h:2206
std::map< std::string, int > _vector_is_adjoint
Holds non-negative if a vector by that name should be projected using adjoint constraints/BCs, -1 if primal.
Definition: system.h:1928
dof_id_type n_local_dofs() const
Definition: system.C:185
std::map< std::string, NumericVector< Number > * > _vectors
Some systems need an arbitrary number of vectors.
Definition: system.h:1916
const Parallel::Communicator & comm() const
std::map< std::string, bool > _vector_projections
Holds true if a vector by that name should be projected onto a changed grid, false if it should be ze...
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 1009 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().

1010 {
1011  std::ostringstream adjoint_name;
1012  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1013 
1014  NumericVector<Number> & returnval = this->add_vector(adjoint_name.str());
1015  this->set_vector_as_adjoint(adjoint_name.str(), i);
1016  return returnval;
1017 }
void set_vector_as_adjoint(const std::string &vec_name, int qoi_num)
Allows one to set the QoI index controlling whether the vector identified by vec_name represents a so...
Definition: system.C:905
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:681
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 956 of file system.C.

References libMesh::System::add_vector().

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

957 {
958  return this->add_vector("weighted_sensitivity_solution");
959 }
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:681
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 708 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(), and main().

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

Assembles & solves the linear system (dR/du)^T*z = dq/du, for those quantities of interest q specified by qoi_indices.

Leave qoi_indices empty to solve all adjoint problems.

Returns
A pair with the total number of linear iterations performed and the (sum of the) final residual norms

Reimplemented from libMesh::System.

Reimplemented in libMesh::DifferentiableSystem.

Definition at line 379 of file implicit_system.C.

References libMesh::System::add_adjoint_solution(), libMesh::LinearSolver< T >::adjoint_solve(), libMesh::System::assemble_before_solve, libMesh::ExplicitSystem::assemble_qoi_derivative(), libMesh::ImplicitSystem::assembly(), libMesh::DofMap::enforce_adjoint_constraints_exactly(), libMesh::System::get_adjoint_rhs(), libMesh::System::get_adjoint_solution(), libMesh::System::get_dof_map(), libMesh::ImplicitSystem::get_linear_solve_parameters(), libMesh::ImplicitSystem::get_linear_solver(), libMesh::QoISet::has_index(), libMesh::ImplicitSystem::matrix, libMesh::System::qoi, and libMesh::ImplicitSystem::release_linear_solver().

Referenced by libMesh::ImplicitSystem::adjoint_qoi_parameter_sensitivity(), libMesh::DifferentiableSystem::adjoint_solve(), libMesh::ImplicitSystem::assembly(), libMesh::ImplicitSystem::qoi_parameter_hessian(), and libMesh::ImplicitSystem::qoi_parameter_hessian_vector_product().

380 {
381  // Log how long the linear solve takes.
382  LOG_SCOPE("adjoint_solve()", "ImplicitSystem");
383 
384  if (this->assemble_before_solve)
385  // Assemble the linear system
386  this->assembly (/* get_residual = */ false,
387  /* get_jacobian = */ true);
388 
389  // The adjoint problem is linear
390  LinearSolver<Number> * linear_solver = this->get_linear_solver();
391 
392  // Reset and build the RHS from the QOI derivative
393  this->assemble_qoi_derivative(qoi_indices,
394  /* include_liftfunc = */ false,
395  /* apply_constraints = */ true);
396 
397  // Our iteration counts and residuals will be sums of the individual
398  // results
399  std::pair<unsigned int, Real> solver_params =
401  std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);
402 
403  for (std::size_t i=0; i != this->qoi.size(); ++i)
404  if (qoi_indices.has_index(i))
405  {
406  const std::pair<unsigned int, Real> rval =
407  linear_solver->adjoint_solve (*matrix, this->add_adjoint_solution(i),
408  this->get_adjoint_rhs(i),
409  solver_params.second,
410  solver_params.first);
411 
412  totalrval.first += rval.first;
413  totalrval.second += rval.second;
414  }
415 
416  this->release_linear_solver(linear_solver);
417 
418  // The linear solver may not have fit our constraints exactly
419 #ifdef LIBMESH_ENABLE_CONSTRAINTS
420  for (std::size_t i=0; i != this->qoi.size(); ++i)
421  if (qoi_indices.has_index(i))
423  (this->get_adjoint_solution(i), i);
424 #endif
425 
426  return totalrval;
427 }
virtual LinearSolver< Number > * get_linear_solver() const
NumericVector< Number > & add_adjoint_solution(unsigned int i=0)
Definition: system.C:977
virtual std::pair< unsigned int, Real > get_linear_solve_parameters() const
virtual void assemble_qoi_derivative(const QoISet &qoi_indices=QoISet(), bool include_liftfunc=true, bool apply_constraints=true) libmesh_override
Prepares adjoint_rhs for quantity of interest derivative assembly, then calls user qoi derivative fun...
virtual void assembly(bool, bool, bool=false, bool=false)
Assembles a residual in rhs and/or a jacobian in matrix, as requested.
const DofMap & get_dof_map() const
Definition: system.h:2030
std::vector< Number > qoi
Values of the quantities of interest.
Definition: system.h:1553
SparseMatrix< Number > * matrix
The system matrix.
virtual void release_linear_solver(LinearSolver< Number > *) const
Releases a pointer to a linear solver acquired by this->get_linear_solver()
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
Definition: system.C:989
bool assemble_before_solve
Flag which tells the system to whether or not to call the user assembly function during each call to ...
Definition: system.h:1477
void enforce_adjoint_constraints_exactly(NumericVector< Number > &v, unsigned int q) const
Heterogenously constrains the numeric vector v, which represents an adjoint solution defined on the m...
Definition: dof_map.h:1820
NumericVector< Number > & get_adjoint_rhs(unsigned int i=0)
Definition: system.C:1051
void libMesh::RBConstruction::allocate_data_structures ( )
protectedvirtualinherited

Helper function that actually allocates all the data structures required by this class.

Reimplemented in libMesh::TransientRBConstruction.

Definition at line 492 of file rb_construction.C.

References libMesh::RBConstruction::Aq_vector, libMesh::DofMap::attach_matrix(), libMesh::SparseMatrix< T >::build(), libMesh::NumericVector< T >::build(), libMesh::ParallelObject::comm(), libMesh::RBConstruction::Fq_representor, libMesh::RBConstruction::Fq_representor_innerprods, libMesh::RBConstruction::Fq_vector, libMesh::System::get_dof_map(), libMesh::RBThetaExpansion::get_n_A_terms(), libMesh::RBThetaExpansion::get_n_F_terms(), libMesh::RBThetaExpansion::get_n_output_terms(), libMesh::RBThetaExpansion::get_n_outputs(), libMesh::RBConstruction::get_rb_theta_expansion(), libMesh::RBConstruction::inner_product_matrix, libMesh::System::n_dofs(), libMesh::System::n_local_dofs(), libMesh::RBConstruction::non_dirichlet_Aq_vector, libMesh::RBConstruction::non_dirichlet_Fq_vector, libMesh::RBConstruction::non_dirichlet_outputs_vector, libMesh::RBConstruction::output_dual_innerprods, libMesh::RBConstruction::outputs_vector, libMesh::PARALLEL, libMesh::RBConstruction::store_non_dirichlet_operators, and libMesh::RBConstruction::truth_outputs.

Referenced by libMesh::TransientRBConstruction::allocate_data_structures(), and libMesh::RBConstruction::initialize_rb_construction().

493 {
494  // Resize vectors for storing mesh-dependent data
495  Aq_vector.resize(get_rb_theta_expansion().get_n_A_terms());
496  Fq_vector.resize(get_rb_theta_expansion().get_n_F_terms());
497 
498  // Resize the Fq_representors and initialize each to NULL
499  // These are basis independent and hence stored here, whereas
500  // the Aq_representors are stored in RBEvaluation
501  Fq_representor.resize(get_rb_theta_expansion().get_n_F_terms());
502 
503  // Initialize vectors for the inner products of the Fq representors
504  // These are basis independent and therefore stored here.
505  unsigned int Q_f_hat = get_rb_theta_expansion().get_n_F_terms()*(get_rb_theta_expansion().get_n_F_terms()+1)/2;
506  Fq_representor_innerprods.resize(Q_f_hat);
507 
508  // Resize the output vectors
509  outputs_vector.resize(get_rb_theta_expansion().get_n_outputs());
510  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
512 
513  // Resize the output dual norm vectors
514  output_dual_innerprods.resize(get_rb_theta_expansion().get_n_outputs());
515  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
516  {
518  output_dual_innerprods[n].resize(Q_l_hat);
519  }
520 
521  {
522  DofMap & dof_map = this->get_dof_map();
523 
524  dof_map.attach_matrix(*inner_product_matrix);
525  inner_product_matrix->init();
526  inner_product_matrix->zero();
527 
528  for (unsigned int q=0; q<get_rb_theta_expansion().get_n_A_terms(); q++)
529  {
530  // Initialize the memory for the matrices
531  Aq_vector[q] = SparseMatrix<Number>::build(this->comm()).release();
532  dof_map.attach_matrix(*Aq_vector[q]);
533  Aq_vector[q]->init();
534  Aq_vector[q]->zero();
535  }
536 
537  // We also need to initialize a second set of non-Dirichlet operators
539  {
540  non_dirichlet_Aq_vector.resize(get_rb_theta_expansion().get_n_A_terms());
541  for (unsigned int q=0; q<get_rb_theta_expansion().get_n_A_terms(); q++)
542  {
543  // Initialize the memory for the matrices
545  dof_map.attach_matrix(*non_dirichlet_Aq_vector[q]);
546  non_dirichlet_Aq_vector[q]->init();
547  non_dirichlet_Aq_vector[q]->zero();
548  }
549  }
550  }
551 
552  // Initialize the vectors
553  for (unsigned int q=0; q<get_rb_theta_expansion().get_n_F_terms(); q++)
554  {
555  // Initialize the memory for the vectors
556  Fq_vector[q] = NumericVector<Number>::build(this->comm()).release();
557  Fq_vector[q]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
558  }
559 
560  // We also need to initialize a second set of non-Dirichlet operators
562  {
563  non_dirichlet_Fq_vector.resize(get_rb_theta_expansion().get_n_F_terms());
564  for (unsigned int q=0; q<get_rb_theta_expansion().get_n_F_terms(); q++)
565  {
566  // Initialize the memory for the vectors
568  non_dirichlet_Fq_vector[q]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
569  }
570  }
571 
572  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
573  for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
574  {
575  // Initialize the memory for the truth output vectors
576  outputs_vector[n][q_l] = (NumericVector<Number>::build(this->comm()).release());
577  outputs_vector[n][q_l]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
578  }
579 
581  {
582  non_dirichlet_outputs_vector.resize(get_rb_theta_expansion().get_n_outputs());
583  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
584  {
585  non_dirichlet_outputs_vector[n].resize( get_rb_theta_expansion().get_n_output_terms(n) );
586  for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
587  {
588  // Initialize the memory for the truth output vectors
589  non_dirichlet_outputs_vector[n][q_l] = (NumericVector<Number>::build(this->comm()).release());
590  non_dirichlet_outputs_vector[n][q_l]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
591  }
592  }
593  }
594 
595  // Resize truth_outputs vector
596  truth_outputs.resize(this->get_rb_theta_expansion().get_n_outputs());
597 }
std::vector< SparseMatrix< Number > * > non_dirichlet_Aq_vector
We may also need a second set of matrices/vectors that do not have the Dirichlet boundary conditions ...
unsigned int get_n_A_terms() const
Get Q_a, the number of terms in the affine expansion for the bilinear form.
unsigned int get_n_output_terms(unsigned int output_index) const
Get the number of affine terms associated with the specified output.
std::vector< std::vector< Number > > output_dual_innerprods
The vector storing the dual norm inner product terms for each output.
static UniquePtr< SparseMatrix< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a SparseMatrix<T> using the linear solver package specified by solver_package.
std::vector< Number > Fq_representor_innerprods
Vectors storing the residual representor inner products to be used in computing the residuals online...
UniquePtr< SparseMatrix< Number > > inner_product_matrix
The inner product matrix.
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
RBThetaExpansion & get_rb_theta_expansion()
Get a reference to the RBThetaExpansion object that that belongs to rb_eval.
std::vector< std::vector< NumericVector< Number > * > > non_dirichlet_outputs_vector
const DofMap & get_dof_map() const
Definition: system.h:2030
unsigned int get_n_F_terms() const
Get Q_f, the number of terms in the affine expansion for the right-hand side.
std::vector< std::vector< NumericVector< Number > * > > outputs_vector
The libMesh vectors that define the output functionals.
std::vector< NumericVector< Number > * > Fq_representor
Vector storing the residual representors associated with the right-hand side.
dof_id_type n_local_dofs() const
Definition: system.C:185
unsigned int get_n_outputs() const
Get n_outputs, the number output functionals.
std::vector< SparseMatrix< Number > * > Aq_vector
Vector storing the Q_a matrices from the affine expansion.
bool store_non_dirichlet_operators
Boolean flag to indicate whether we store a second copy of each affine operator and vector which does...
const Parallel::Communicator & comm() const
std::vector< Number > truth_outputs
Vector storing the truth output values from the most recent truth solve.
dof_id_type n_dofs() const
Definition: system.C:148
std::vector< NumericVector< Number > * > non_dirichlet_Fq_vector
std::vector< NumericVector< Number > * > Fq_vector
Vector storing the Q_f vectors in the affine decomposition of the right-hand side.
virtual void libMesh::LinearImplicitSystem::assemble ( )
virtualinherited

Prepares matrix and _dof_map for matrix assembly.

Does not actually assemble anything. For matrix assembly, use the assemble() in derived classes. Should be overridden in derived classes.

Reimplemented from libMesh::ImplicitSystem.

Reimplemented in libMesh::FrequencySystem, and libMesh::NewmarkSystem.

Definition at line 104 of file linear_implicit_system.h.

References libMesh::ImplicitSystem::assemble(), libMesh::LinearImplicitSystem::assembly(), libMesh::LinearImplicitSystem::get_linear_solver(), libMesh::LinearImplicitSystem::release_linear_solver(), libMesh::LinearImplicitSystem::restrict_solve_to(), libMesh::LinearImplicitSystem::solve(), and libMesh::SUBSET_ZERO.

Referenced by libMesh::NewmarkSystem::assemble(), libMesh::FrequencySystem::assemble(), libMesh::LinearImplicitSystem::assembly(), and libMesh::LinearImplicitSystem::solve().

virtual void assemble() libmesh_override
Prepares matrix and rhs for system assembly, then calls user assembly function.
void libMesh::RBConstruction::assemble_affine_expansion ( bool  skip_matrix_assembly,
bool  skip_vector_assembly 
)
virtualinherited

Assemble the matrices and vectors for this system.

Optionally skip matrix or vector assembly (e.g. we may want to read data in from disk instead).

Reimplemented in libMesh::TransientRBConstruction.

Definition at line 474 of file rb_construction.C.

References libMesh::RBConstruction::assemble_all_affine_operators(), libMesh::RBConstruction::assemble_all_affine_vectors(), libMesh::RBConstruction::assemble_all_output_vectors(), and libMesh::RBConstruction::assemble_misc_matrices().

Referenced by libMesh::TransientRBConstruction::assemble_affine_expansion(), libMesh::RBConstruction::initialize_rb_construction(), and libMesh::RBConstruction::is_quiet().

476 {
477  if (!skip_matrix_assembly)
478  {
479  // Assemble and store all of the matrices
480  this->assemble_misc_matrices();
482  }
483 
484  if (!skip_vector_assembly)
485  {
486  // Assemble and store all of the vectors
489  }
490 }
virtual void assemble_all_output_vectors()
Assemble and store the output vectors.
virtual void assemble_misc_matrices()
Assemble and store all the inner-product matrix, the constraint matrix (for constrained problems) and...
virtual void assemble_all_affine_vectors()
Assemble and store the affine RHS vectors.
virtual void assemble_all_affine_operators()
Assemble and store all Q_a affine operators as well as the inner-product matrix.
void libMesh::RBConstruction::assemble_all_affine_operators ( )
protectedvirtualinherited

Assemble and store all Q_a affine operators as well as the inner-product matrix.

Reimplemented in libMesh::TransientRBConstruction.

Definition at line 949 of file rb_construction.C.

References libMesh::RBConstruction::assemble_Aq_matrix(), libMesh::RBConstruction::get_Aq(), libMesh::RBThetaExpansion::get_n_A_terms(), libMesh::RBConstruction::get_non_dirichlet_Aq(), libMesh::RBConstruction::get_rb_theta_expansion(), libMesh::out, and libMesh::RBConstruction::store_non_dirichlet_operators.

Referenced by libMesh::RBConstruction::assemble_affine_expansion(), and libMesh::TransientRBConstruction::assemble_all_affine_operators().

950 {
951  for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
952  {
953  libMesh::out << "Assembling affine operator " << (q_a+1) << " of "
954  << get_rb_theta_expansion().get_n_A_terms() << std::endl;
955  assemble_Aq_matrix(q_a, get_Aq(q_a));
956  }
957 
959  {
960  for (unsigned int q_a=0; q_a<get_rb_theta_expansion().get_n_A_terms(); q_a++)
961  {
962  libMesh::out << "Assembling non-Dirichlet affine operator " << (q_a+1) << " of "
963  << get_rb_theta_expansion().get_n_A_terms() << std::endl;
964  assemble_Aq_matrix(q_a, get_non_dirichlet_Aq(q_a), false);
965  }
966  }
967 }
unsigned int get_n_A_terms() const
Get Q_a, the number of terms in the affine expansion for the bilinear form.
SparseMatrix< Number > * get_Aq(unsigned int q)
Get a pointer to Aq.
RBThetaExpansion & get_rb_theta_expansion()
Get a reference to the RBThetaExpansion object that that belongs to rb_eval.
SparseMatrix< Number > * get_non_dirichlet_Aq(unsigned int q)
Get a pointer to non_dirichlet_Aq.
bool store_non_dirichlet_operators
Boolean flag to indicate whether we store a second copy of each affine operator and vector which does...
OStreamProxy out
void assemble_Aq_matrix(unsigned int q, SparseMatrix< Number > *input_matrix, bool apply_dof_constraints=true)
Assemble the q^th affine matrix and store it in input_matrix.
void libMesh::RBConstruction::assemble_all_affine_vectors ( )
protectedvirtualinherited

Assemble and store the affine RHS vectors.

Definition at line 969 of file rb_construction.C.

References libMesh::RBConstruction::assemble_Fq_vector(), libMesh::RBConstruction::get_Fq(), libMesh::RBThetaExpansion::get_n_F_terms(), libMesh::RBConstruction::get_non_dirichlet_Fq(), libMesh::RBConstruction::get_rb_theta_expansion(), libMesh::out, and libMesh::RBConstruction::store_non_dirichlet_operators.

Referenced by libMesh::RBConstruction::assemble_affine_expansion().

970 {
971  for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
972  {
973  libMesh::out << "Assembling affine vector " << (q_f+1) << " of "
974  << get_rb_theta_expansion().get_n_F_terms() << std::endl;
975  assemble_Fq_vector(q_f, get_Fq(q_f));
976  }
977 
979  {
980  for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
981  {
982  libMesh::out << "Assembling non-Dirichlet affine vector " << (q_f+1) << " of "
983  << get_rb_theta_expansion().get_n_F_terms() << std::endl;
984  assemble_Fq_vector(q_f, get_non_dirichlet_Fq(q_f), false);
985  }
986  }
987 
988 }
NumericVector< Number > * get_Fq(unsigned int q)
Get a pointer to Fq.
RBThetaExpansion & get_rb_theta_expansion()
Get a reference to the RBThetaExpansion object that that belongs to rb_eval.
unsigned int get_n_F_terms() const
Get Q_f, the number of terms in the affine expansion for the right-hand side.
void assemble_Fq_vector(unsigned int q, NumericVector< Number > *input_vector, bool apply_dof_constraints=true)
Assemble the q^th affine vector and store it in input_matrix.
NumericVector< Number > * get_non_dirichlet_Fq(unsigned int q)
Get a pointer to non-Dirichlet Fq.
bool store_non_dirichlet_operators
Boolean flag to indicate whether we store a second copy of each affine operator and vector which does...
OStreamProxy out
void libMesh::RBConstruction::assemble_all_output_vectors ( )
protectedvirtualinherited

Assemble and store the output vectors.

Definition at line 1007 of file rb_construction.C.

References libMesh::RBConstruction::add_scaled_matrix_and_vector(), libMesh::RBThetaExpansion::get_n_output_terms(), libMesh::RBThetaExpansion::get_n_outputs(), libMesh::RBConstruction::get_non_dirichlet_output_vector(), libMesh::RBAssemblyExpansion::get_output_assembly(), libMesh::RBConstruction::get_output_vector(), libMesh::RBConstruction::get_rb_theta_expansion(), libmesh_nullptr, libMesh::out, libMesh::RBConstruction::rb_assembly_expansion, libMesh::RBConstruction::store_non_dirichlet_operators, and libMesh::NumericVector< T >::zero().

Referenced by libMesh::RBConstruction::assemble_affine_expansion().

1008 {
1009  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1010  for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
1011  {
1012  libMesh::out << "Assembling output vector, (" << (n+1) << "," << (q_l+1)
1013  << ") of (" << get_rb_theta_expansion().get_n_outputs()
1014  << "," << get_rb_theta_expansion().get_n_output_terms(n) << ")"
1015  << std::endl;
1016  get_output_vector(n, q_l)->zero();
1019  get_output_vector(n,q_l),
1020  false, /* symmetrize */
1021  true /* apply_dof_constraints */);
1022  }
1023 
1025  {
1026  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1027  for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
1028  {
1029  libMesh::out << "Assembling non-Dirichlet output vector, (" << (n+1) << "," << (q_l+1)
1030  << ") of (" << get_rb_theta_expansion().get_n_outputs()
1031  << "," << get_rb_theta_expansion().get_n_output_terms(n) << ")"
1032  << std::endl;
1037  false, /* symmetrize */
1038  false /* apply_dof_constraints */);
1039  }
1040  }
1041 }
unsigned int get_n_output_terms(unsigned int output_index) const
Get the number of affine terms associated with the specified output.
const class libmesh_nullptr_t libmesh_nullptr
ElemAssembly & get_output_assembly(unsigned int output_index, unsigned int q_l)
Return a reference to the specified output assembly object.
virtual void zero()=0
Set all entries to zero.
RBThetaExpansion & get_rb_theta_expansion()
Get a reference to the RBThetaExpansion object that that belongs to rb_eval.
void add_scaled_matrix_and_vector(Number scalar, ElemAssembly *elem_assembly, SparseMatrix< Number > *input_matrix, NumericVector< Number > *input_vector, bool symmetrize=false, bool apply_dof_constraints=true)
This function loops over the mesh and applies the specified interior and/or boundary assembly routine...
unsigned int get_n_outputs() const
Get n_outputs, the number output functionals.
bool store_non_dirichlet_operators
Boolean flag to indicate whether we store a second copy of each affine operator and vector which does...
OStreamProxy out
RBAssemblyExpansion * rb_assembly_expansion
This member holds the (parameter independent) assembly functors that define the "affine expansion" of...
NumericVector< Number > * get_output_vector(unsigned int n, unsigned int q_l)
Get a pointer to the n^th output.
NumericVector< Number > * get_non_dirichlet_output_vector(unsigned int n, unsigned int q_l)
Get a pointer to non-Dirichlet output vector.
void libMesh::RBConstruction::assemble_Aq_matrix ( unsigned int  q,
SparseMatrix< Number > *  input_matrix,
bool  apply_dof_constraints = true 
)
inherited

Assemble the q^th affine matrix and store it in input_matrix.

Definition at line 901 of file rb_construction.C.

References libMesh::RBConstruction::add_scaled_matrix_and_vector(), libMesh::RBAssemblyExpansion::get_A_assembly(), libMesh::RBConstruction::get_rb_theta_expansion(), libmesh_nullptr, libMesh::RBConstruction::rb_assembly_expansion, and libMesh::SparseMatrix< T >::zero().

Referenced by libMesh::RBConstruction::assemble_all_affine_operators(), and libMesh::RBConstruction::is_quiet().

904 {
905  if (q >= get_rb_theta_expansion().get_n_A_terms())
906  libmesh_error_msg("Error: We must have q < Q_a in assemble_Aq_matrix.");
907 
908  input_matrix->zero();
909 
912  input_matrix,
914  false, /* symmetrize */
915  apply_dof_constraints);
916 }
const class libmesh_nullptr_t libmesh_nullptr
RBThetaExpansion & get_rb_theta_expansion()
Get a reference to the RBThetaExpansion object that that belongs to rb_eval.
virtual void zero()=0
Set all entries to 0.
void add_scaled_matrix_and_vector(Number scalar, ElemAssembly *elem_assembly, SparseMatrix< Number > *input_matrix, NumericVector< Number > *input_vector, bool symmetrize=false, bool apply_dof_constraints=true)
This function loops over the mesh and applies the specified interior and/or boundary assembly routine...
ElemAssembly & get_A_assembly(unsigned int q)
Return a reference to the specified A_assembly object.
RBAssemblyExpansion * rb_assembly_expansion
This member holds the (parameter independent) assembly functors that define the "affine expansion" of...
void libMesh::RBConstruction::assemble_Fq_vector ( unsigned int  q,
NumericVector< Number > *  input_vector,
bool  apply_dof_constraints = true 
)
inherited

Assemble the q^th affine vector and store it in input_matrix.

Definition at line 990 of file rb_construction.C.

References libMesh::RBConstruction::add_scaled_matrix_and_vector(), libMesh::RBAssemblyExpansion::get_F_assembly(), libMesh::RBConstruction::get_rb_theta_expansion(), libmesh_nullptr, libMesh::RBConstruction::rb_assembly_expansion, and libMesh::NumericVector< T >::zero().

Referenced by libMesh::RBConstruction::assemble_all_affine_vectors(), and libMesh::RBConstruction::is_quiet().

993 {
994  if (q >= get_rb_theta_expansion().get_n_F_terms())
995  libmesh_error_msg("Error: We must have q < Q_f in assemble_Fq_vector.");
996 
997  input_vector->zero();
998 
1002  input_vector,
1003  false, /* symmetrize */
1004  apply_dof_constraints /* apply_dof_constraints */);
1005 }
ElemAssembly & get_F_assembly(unsigned int q)
Return a reference to the specified F_assembly object.
const class libmesh_nullptr_t libmesh_nullptr
virtual void zero()=0
Set all entries to zero.
RBThetaExpansion & get_rb_theta_expansion()
Get a reference to the RBThetaExpansion object that that belongs to rb_eval.
void add_scaled_matrix_and_vector(Number scalar, ElemAssembly *elem_assembly, SparseMatrix< Number > *input_matrix, NumericVector< Number > *input_vector, bool symmetrize=false, bool apply_dof_constraints=true)
This function loops over the mesh and applies the specified interior and/or boundary assembly routine...
RBAssemblyExpansion * rb_assembly_expansion
This member holds the (parameter independent) assembly functors that define the "affine expansion" of...
void libMesh::RBConstruction::assemble_inner_product_matrix ( SparseMatrix< Number > *  input_matrix,
bool  apply_dof_constraints = true 
)
inherited

Assemble the inner product matrix and store it in input_matrix.

Definition at line 889 of file rb_construction.C.

References libMesh::RBConstruction::add_scaled_matrix_and_vector(), libMesh::RBConstruction::inner_product_assembly, libmesh_nullptr, and libMesh::SparseMatrix< T >::zero().

Referenced by libMesh::RBConstruction::assemble_misc_matrices(), and libMesh::RBConstruction::is_quiet().

891 {
892  input_matrix->zero();
895  input_matrix,
897  false, /* symmetrize */
898  apply_dof_constraints);
899 }
const class libmesh_nullptr_t libmesh_nullptr
virtual void zero()=0
Set all entries to 0.
void add_scaled_matrix_and_vector(Number scalar, ElemAssembly *elem_assembly, SparseMatrix< Number > *input_matrix, NumericVector< Number > *input_vector, bool symmetrize=false, bool apply_dof_constraints=true)
This function loops over the mesh and applies the specified interior and/or boundary assembly routine...
ElemAssembly * inner_product_assembly
Pointer to inner product assembly.
void libMesh::RBConstruction::assemble_misc_matrices ( )
protectedvirtualinherited

Assemble and store all the inner-product matrix, the constraint matrix (for constrained problems) and the mass matrix (for time-dependent problems).

Reimplemented in libMesh::TransientRBConstruction.

Definition at line 943 of file rb_construction.C.

References libMesh::RBConstruction::assemble_inner_product_matrix(), libMesh::RBConstruction::inner_product_matrix, and libMesh::out.

Referenced by libMesh::RBConstruction::assemble_affine_expansion(), and libMesh::TransientRBConstruction::assemble_misc_matrices().

944 {
945  libMesh::out << "Assembling inner product matrix" << std::endl;
947 }
UniquePtr< SparseMatrix< Number > > inner_product_matrix
The inner product matrix.
void assemble_inner_product_matrix(SparseMatrix< Number > *input_matrix, bool apply_dof_constraints=true)
Assemble the inner product matrix and store it in input_matrix.
OStreamProxy out
void libMesh::ExplicitSystem::assemble_qoi ( const QoISet qoi_indices = QoISet())
virtualinherited

Prepares qoi for quantity of interest assembly, then calls user qoi function.

Can be overridden in derived classes.

Reimplemented from libMesh::System.

Reimplemented in libMesh::FEMSystem.

Definition at line 68 of file explicit_system.C.

References libMesh::System::assemble_qoi(), libMesh::QoISet::has_index(), and libMesh::System::qoi.

Referenced by libMesh::ImplicitSystem::adjoint_qoi_parameter_sensitivity(), libMesh::ImplicitSystem::forward_qoi_parameter_sensitivity(), libMesh::ImplicitSystem::qoi_parameter_hessian(), libMesh::ImplicitSystem::qoi_parameter_hessian_vector_product(), and libMesh::ExplicitSystem::system().

69 {
70  // The user quantity of interest assembly gets to expect to
71  // accumulate on initially zero values
72  for (std::size_t i=0; i != qoi.size(); ++i)
73  if (qoi_indices.has_index(i))
74  qoi[i] = 0;
75 
76  Parent::assemble_qoi (qoi_indices);
77 }
std::vector< Number > qoi
Values of the quantities of interest.
Definition: system.h:1553
virtual void assemble_qoi(const QoISet &qoi_indices=QoISet())
Calls user qoi function.
Definition: system.C:490
void libMesh::ExplicitSystem::assemble_qoi_derivative ( const QoISet qoi_indices = QoISet(),
bool  include_liftfunc = true,
bool  apply_constraints = true 
)
virtualinherited

Prepares adjoint_rhs for quantity of interest derivative assembly, then calls user qoi derivative function.

Can be overridden in derived classes.

Reimplemented from libMesh::System.

Reimplemented in libMesh::FEMSystem.

Definition at line 81 of file explicit_system.C.

References libMesh::System::add_adjoint_rhs(), libMesh::System::assemble_qoi_derivative(), libMesh::QoISet::has_index(), libMesh::System::qoi, and libMesh::NumericVector< T >::zero().

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

84 {
85  // The user quantity of interest derivative assembly gets to expect
86  // to accumulate on initially zero vectors
87  for (std::size_t i=0; i != qoi.size(); ++i)
88  if (qoi_indices.has_index(i))
89  this->add_adjoint_rhs(i).zero();
90 
91  Parent::assemble_qoi_derivative (qoi_indices, include_liftfunc,
92  apply_constraints);
93 }
virtual void assemble_qoi_derivative(const QoISet &qoi_indices=QoISet(), bool include_liftfunc=true, bool apply_constraints=true)
Calls user qoi derivative function.
Definition: system.C:501
virtual void zero()=0
Set all entries to zero.
std::vector< Number > qoi
Values of the quantities of interest.
Definition: system.h:1553
NumericVector< Number > & add_adjoint_rhs(unsigned int i=0)
Definition: system.C:1041
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 665 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().

666 {
667  ParameterVector & parameters =
668  const_cast<ParameterVector &>(parameters_in);
669 
670  const unsigned int Np = cast_int<unsigned int>
671  (parameters.size());
672 
673  for (unsigned int p=0; p != Np; ++p)
674  {
675  NumericVector<Number> & sensitivity_rhs = this->add_sensitivity_rhs(p);
676 
677  // Approximate -(partial R / partial p) by
678  // (R(p-dp) - R(p+dp)) / (2*dp)
679 
680  Number old_parameter = *parameters[p];
681 
682  const Real delta_p =
683  TOLERANCE * std::max(std::abs(old_parameter), 1e-3);
684 
685  *parameters[p] -= delta_p;
686 
687  // this->assembly(true, false, true);
688  this->assembly(true, false, false);
689  this->rhs->close();
690  sensitivity_rhs = *this->rhs;
691 
692  *parameters[p] = old_parameter + delta_p;
693 
694  // this->assembly(true, false, true);
695  this->assembly(true, false, false);
696  this->rhs->close();
697 
698  sensitivity_rhs -= *this->rhs;
699  sensitivity_rhs /= (2*delta_p);
700  sensitivity_rhs.close();
701 
702  *parameters[p] = old_parameter;
703  }
704 }
double abs(double a)
NumericVector< Number > * rhs
The system matrix.
NumericVector< Number > & add_sensitivity_rhs(unsigned int i=0)
Definition: system.C:1071
static const Real TOLERANCE
long double max(long double a, double b)
virtual void assembly(bool, bool, bool=false, bool=false)
Assembles a residual in rhs and/or a jacobian in matrix, as requested.
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void libMesh::LinearImplicitSystem::assembly ( bool  get_residual,
bool  get_jacobian,
bool  apply_heterogeneous_constraints = false,
bool  apply_no_constraints = false 
)
virtualinherited

Assembles a residual in rhs and/or a jacobian in matrix, as requested.

Reimplemented from libMesh::ImplicitSystem.

Definition at line 366 of file linear_implicit_system.C.

References libMesh::NumericVector< T >::add_vector(), libMesh::LinearImplicitSystem::assemble(), libMesh::NumericVector< T >::close(), libMesh::SparseMatrix< T >::close(), libMesh::ImplicitSystem::matrix, libMesh::ExplicitSystem::rhs, and libMesh::System::solution.

Referenced by libMesh::LinearImplicitSystem::assemble().

370 {
371  // Residual R(u(p),p) := A(p)*u(p) - b(p)
372  // partial R / partial u = A
373 
374  this->assemble();
375  this->rhs->close();
376  this->matrix->close();
377 
378  *(this->rhs) *= -1.0;
379  this->rhs->add_vector(*this->solution, *this->matrix);
380 }
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
NumericVector< Number > * rhs
The system matrix.
virtual void assemble() libmesh_override
Prepares matrix and _dof_map for matrix assembly.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
SparseMatrix< Number > * matrix
The system matrix.
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 1795 of file system.C.

References libMesh::System::_assemble_system_function, libMesh::System::_assemble_system_object, fptr(), libMesh::libmesh_assert(), libmesh_nullptr, and libMesh::out.

Referenced by assemble_and_solve(), main(), and libMesh::System::read_parallel_data().

1797 {
1799 
1801  {
1802  libmesh_here();
1803  libMesh::out << "WARNING: Cannot specify both assembly function and object!"
1804  << std::endl;
1805 
1807  }
1808 
1810 }
Assembly * _assemble_system_object
Object that assembles the system.
Definition: system.h:1822
const class libmesh_nullptr_t libmesh_nullptr
Number fptr(const Point &p, const Parameters &, const std::string &libmesh_dbg_var(sys_name), const std::string &unknown_name)
Definition: projection.C:78
libmesh_assert(j)
OStreamProxy out
void(* _assemble_system_function)(EquationSystems &es, const std::string &name)
Function that assembles the system.
Definition: system.h:1816
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 1814 of file system.C.

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

Referenced by main(), and libMesh::System::read_parallel_data().

1815 {
1817  {
1818  libmesh_here();
1819  libMesh::out << "WARNING: Cannot specify both assembly object and function!"
1820  << std::endl;
1821 
1823  }
1824 
1825  _assemble_system_object = &assemble_in;
1826 }
Assembly * _assemble_system_object
Object that assembles the system.
Definition: system.h:1822
const class libmesh_nullptr_t libmesh_nullptr
OStreamProxy out
void(* _assemble_system_function)(EquationSystems &es, const std::string &name)
Function that assembles the system.
Definition: system.h:1816
void libMesh::System::attach_constraint_function ( void   fptrEquationSystems &es, const std::string &name)
inherited

Register a user function for imposing constraints.

Definition at line 1830 of file system.C.

References libMesh::System::_constrain_system_function, libMesh::System::_constrain_system_object, fptr(), libMesh::libmesh_assert(), libmesh_nullptr, and libMesh::out.

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

1832 {
1834 
1836  {
1837  libmesh_here();
1838  libMesh::out << "WARNING: Cannot specify both constraint function and object!"
1839  << std::endl;
1840 
1842  }
1843 
1845 }
void(* _constrain_system_function)(EquationSystems &es, const std::string &name)
Function to impose constraints.
Definition: system.h:1827
Constraint * _constrain_system_object
Object that constrains the system.
Definition: system.h:1833
const class libmesh_nullptr_t libmesh_nullptr
Number fptr(const Point &p, const Parameters &, const std::string &libmesh_dbg_var(sys_name), const std::string &unknown_name)
Definition: projection.C:78
libmesh_assert(j)
OStreamProxy out
void libMesh::System::attach_constraint_object ( System::Constraint constrain)
inherited

Register a user object for imposing constraints.

Definition at line 1849 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().

1850 {
1852  {
1853  libmesh_here();
1854  libMesh::out << "WARNING: Cannot specify both constraint object and function!"
1855  << std::endl;
1856 
1858  }
1859 
1860  _constrain_system_object = &constrain;
1861 }
void(* _constrain_system_function)(EquationSystems &es, const std::string &name)
Function to impose constraints.
Definition: system.h:1827
Constraint * _constrain_system_object
Object that constrains the system.
Definition: system.h:1833
const class libmesh_nullptr_t libmesh_nullptr
OStreamProxy out
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 1760 of file system.C.

References libMesh::System::_init_system_function, libMesh::System::_init_system_object, fptr(), libMesh::libmesh_assert(), libmesh_nullptr, and libMesh::out.

Referenced by main(), and libMesh::System::read_parallel_data().

1762 {
1764 
1766  {
1767  libmesh_here();
1768  libMesh::out << "WARNING: Cannot specify both initialization function and object!"
1769  << std::endl;
1770 
1772  }
1773 
1775 }
const class libmesh_nullptr_t libmesh_nullptr
Number fptr(const Point &p, const Parameters &, const std::string &libmesh_dbg_var(sys_name), const std::string &unknown_name)
Definition: projection.C:78
libmesh_assert(j)
Initialization * _init_system_object
Object that initializes the system.
Definition: system.h:1811
void(* _init_system_function)(EquationSystems &es, const std::string &name)
Function that initializes the system.
Definition: system.h:1805
OStreamProxy out
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 1779 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().

1780 {
1782  {
1783  libmesh_here();
1784  libMesh::out << "WARNING: Cannot specify both initialization object and function!"
1785  << std::endl;
1786 
1788  }
1789 
1790  _init_system_object = &init_in;
1791 }
const class libmesh_nullptr_t libmesh_nullptr
Initialization * _init_system_object
Object that initializes the system.
Definition: system.h:1811
void(* _init_system_function)(EquationSystems &es, const std::string &name)
Function that initializes the system.
Definition: system.h:1805
OStreamProxy out
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 1901 of file system.C.

References libMesh::System::_qoi_evaluate_derivative_function, libMesh::System::_qoi_evaluate_derivative_object, fptr(), libMesh::libmesh_assert(), libmesh_nullptr, and libMesh::out.

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

1903 {
1905 
1907  {
1908  libmesh_here();
1909  libMesh::out << "WARNING: Cannot specify both QOI derivative function and object!"
1910  << std::endl;
1911 
1913  }
1914 
1916 }
const class libmesh_nullptr_t libmesh_nullptr
Number fptr(const Point &p, const Parameters &, const std::string &libmesh_dbg_var(sys_name), const std::string &unknown_name)
Definition: projection.C:78
libmesh_assert(j)
QOIDerivative * _qoi_evaluate_derivative_object
Object to compute derivatives of quantities of interest.
Definition: system.h:1859
OStreamProxy out
void(* _qoi_evaluate_derivative_function)(EquationSystems &es, const std::string &name, const QoISet &qoi_indices, bool include_liftfunc, bool apply_constraints)
Function to evaluate quantity of interest derivative.
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 1920 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().

1921 {
1923  {
1924  libmesh_here();
1925  libMesh::out << "WARNING: Cannot specify both QOI derivative object and function!"
1926  << std::endl;
1927 
1929  }
1930 
1931  _qoi_evaluate_derivative_object = &qoi_derivative;
1932 }
const class libmesh_nullptr_t libmesh_nullptr
QOIDerivative * _qoi_evaluate_derivative_object
Object to compute derivatives of quantities of interest.
Definition: system.h:1859
OStreamProxy out
void(* _qoi_evaluate_derivative_function)(EquationSystems &es, const std::string &name, const QoISet &qoi_indices, bool include_liftfunc, bool apply_constraints)
Function to evaluate quantity of interest derivative.
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 1865 of file system.C.

References libMesh::System::_qoi_evaluate_function, libMesh::System::_qoi_evaluate_object, fptr(), libMesh::libmesh_assert(), libmesh_nullptr, and libMesh::out.

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

1868 {
1870 
1872  {
1873  libmesh_here();
1874  libMesh::out << "WARNING: Cannot specify both QOI function and object!"
1875  << std::endl;
1876 
1878  }
1879 
1881 }
const class libmesh_nullptr_t libmesh_nullptr
void(* _qoi_evaluate_function)(EquationSystems &es, const std::string &name, const QoISet &qoi_indices)
Function to evaluate quantity of interest.
Definition: system.h:1838
Number fptr(const Point &p, const Parameters &, const std::string &libmesh_dbg_var(sys_name), const std::string &unknown_name)
Definition: projection.C:78
libmesh_assert(j)
QOI * _qoi_evaluate_object
Object to compute quantities of interest.
Definition: system.h:1845
OStreamProxy out
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 1885 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().

1886 {
1888  {
1889  libmesh_here();
1890  libMesh::out << "WARNING: Cannot specify both QOI object and function!"
1891  << std::endl;
1892 
1894  }
1895 
1896  _qoi_evaluate_object = &qoi_in;
1897 }
const class libmesh_nullptr_t libmesh_nullptr
void(* _qoi_evaluate_function)(EquationSystems &es, const std::string &name, const QoISet &qoi_indices)
Function to evaluate quantity of interest.
Definition: system.h:1838
QOI * _qoi_evaluate_object
Object to compute quantities of interest.
Definition: system.h:1845
OStreamProxy out
void libMesh::LinearImplicitSystem::attach_shell_matrix ( ShellMatrix< Number > *  shell_matrix)
inherited

This function enables the user to provide a shell matrix, i.e.

a matrix that is not stored element-wise, but as a function. When you register your shell matrix using this function, calling solve() will no longer use the matrix member but the registered shell matrix instead. You can reset this behaviour to its original state by supplying a NULL pointer to this function.

Definition at line 158 of file linear_implicit_system.C.

References libMesh::LinearImplicitSystem::_shell_matrix.

Referenced by libMesh::LinearImplicitSystem::detach_shell_matrix(), libMesh::LinearImplicitSystem::final_linear_residual(), and main().

159 {
160  _shell_matrix = shell_matrix;
161 }
ShellMatrix< Number > * _shell_matrix
User supplies shell matrix or NULL if no shell matrix is used.
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.

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

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.

Definition at line 999 of file system_projection.C.

References libMesh::System::boundary_project_vector(), fptr(), and gptr().

Referenced by libMesh::System::project_vector(), and libMesh::System::system_type().

1003 {
1004  this->boundary_project_vector(b, variables, *solution, f, g);
1005 
1006  solution->localize(*current_local_solution);
1007 }
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
Projects arbitrary boundary functions onto a vector of degree of freedom values for the current syste...
UniquePtr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1535
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
void libMesh::System::boundary_project_solution ( const std::set< boundary_id_type > &  b,
const std::vector< unsigned int > &  variables,
Number   fptrconst Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name,
Gradient   gptrconst Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name,
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.

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.

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

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.

Definition at line 1041 of file system_projection.C.

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

Referenced by libMesh::System::boundary_project_solution(), and libMesh::System::system_type().

1047 {
1048  LOG_SCOPE ("boundary_project_vector()", "System");
1049 
1051  (ConstElemRange (this->get_mesh().active_local_elements_begin(),
1052  this->get_mesh().active_local_elements_end() ),
1053  BoundaryProjectSolution(b, variables, *this, f, g,
1054  this->get_equation_systems().parameters,
1055  new_vector)
1056  );
1057 
1058  // We don't do SCALAR dofs when just projecting the boundary, so
1059  // we're done here.
1060 
1061  new_vector.close();
1062 
1063 #ifdef LIBMESH_ENABLE_CONSTRAINTS
1064  if (is_adjoint == -1)
1065  this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
1066  else if (is_adjoint >= 0)
1068  is_adjoint);
1069 #endif
1070 }
void parallel_for(const Range &range, const Body &body)
Execute the provided function object in parallel on the specified range.
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
Constrains the numeric vector v, which represents a solution defined on the mesh. ...
Definition: dof_map.h:1816
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
const EquationSystems & get_equation_systems() const
Definition: system.h:712
void enforce_adjoint_constraints_exactly(NumericVector< Number > &v, unsigned int q) const
Heterogenously constrains the numeric vector v, which represents an adjoint solution defined on the m...
Definition: dof_map.h:1820
void libMesh::System::boundary_project_vector ( const std::set< boundary_id_type > &  b,
const std::vector< unsigned int > &  variables,
Number   fptrconst Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name,
Gradient   gptrconst Point &p, const Parameters &parameters, const std::string &sys_name, const std::string &unknown_name,
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.

void libMesh::RBConstructionBase< LinearImplicitSystem >::broadcast_parameters ( unsigned int  proc_id)
inherited

Broadcasts parameters on processor proc_id to all processors.

Referenced by libMesh::RBConstruction::compute_max_error_bound().

UniquePtr< DGFEMContext > libMesh::RBConstruction::build_context ( )
protectedvirtualinherited

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

We use DGFEMContext since it allows for both DG and continuous Galerkin formulations.

Definition at line 599 of file rb_construction.C.

Referenced by libMesh::RBConstruction::add_scaled_matrix_and_vector().

600 {
601  return UniquePtr<DGFEMContext>(new DGFEMContext(*this));
602 }
UniquePtr< DirichletBoundary > libMesh::RBConstruction::build_zero_dirichlet_boundary_object ( )
staticinherited

It's helpful to be able to generate a DirichletBoundary that stores a ZeroFunction in order to impose Dirichlet boundary conditions.

Definition at line 1946 of file rb_construction.C.

Referenced by libMesh::RBConstruction::get_delta_N().

1947 {
1948  ZeroFunction<> zf;
1949 
1950  std::set<boundary_id_type> dirichlet_ids;
1951  std::vector<unsigned int> variables;
1952 
1953  // The DirichletBoundary constructor clones zf, so it's OK that zf is only in local scope
1954  return UniquePtr<DirichletBoundary> (new DirichletBoundary(dirichlet_ids, variables, &zf));
1955 }
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 1383 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(), main(), output_norms(), and libMesh::System::project_solution_on_reinit().

1387 {
1388  //short circuit to save time
1389  if (norm_type == DISCRETE_L1 ||
1390  norm_type == DISCRETE_L2 ||
1391  norm_type == DISCRETE_L_INF)
1392  return discrete_var_norm(v,var,norm_type);
1393 
1394  // Not a discrete norm
1395  std::vector<FEMNormType> norms(this->n_vars(), L2);
1396  std::vector<Real> weights(this->n_vars(), 0.0);
1397  norms[var] = norm_type;
1398  weights[var] = 1.0;
1399  Real val = this->calculate_norm(v, SystemNorm(norms, weights), skip_dimensions);
1400  return val;
1401 }
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:1383
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
Finds the discrete norm for the entries in the vector corresponding to Dofs associated with var...
Definition: system.C:1364
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 1405 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(), dim, 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::libmesh_assert(), 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().

1408 {
1409  // This function must be run on all processors at once
1410  parallel_object_only();
1411 
1412  LOG_SCOPE ("calculate_norm()", "System");
1413 
1414  // Zero the norm before summation
1415  Real v_norm = 0.;
1416 
1417  if (norm.is_discrete())
1418  {
1419  //Check to see if all weights are 1.0 and all types are equal
1420  FEMNormType norm_type0 = norm.type(0);
1421  unsigned int check_var = 0;
1422  for (; check_var != this->n_vars(); ++check_var)
1423  if ((norm.weight(check_var) != 1.0) || (norm.type(check_var) != norm_type0))
1424  break;
1425 
1426  //All weights were 1.0 so just do the full vector discrete norm
1427  if (check_var == this->n_vars())
1428  {
1429  if (norm_type0 == DISCRETE_L1)
1430  return v.l1_norm();
1431  if (norm_type0 == DISCRETE_L2)
1432  return v.l2_norm();
1433  if (norm_type0 == DISCRETE_L_INF)
1434  return v.linfty_norm();
1435  else
1436  libmesh_error_msg("Invalid norm_type0 = " << norm_type0);
1437  }
1438 
1439  for (unsigned int var=0; var != this->n_vars(); ++var)
1440  {
1441  // Skip any variables we don't need to integrate
1442  if (norm.weight(var) == 0.0)
1443  continue;
1444 
1445  v_norm += norm.weight(var) * discrete_var_norm(v, var, norm.type(var));
1446  }
1447 
1448  return v_norm;
1449  }
1450 
1451  // Localize the potentially parallel vector
1452  UniquePtr<NumericVector<Number>> local_v = NumericVector<Number>::build(this->comm());
1453  local_v->init(v.size(), true, SERIAL);
1454  v.localize (*local_v, _dof_map->get_send_list());
1455 
1456  // I'm not sure how best to mix Hilbert norms on some variables (for
1457  // which we'll want to square then sum then square root) with norms
1458  // like L_inf (for which we'll just want to take an absolute value
1459  // and then sum).
1460  bool using_hilbert_norm = true,
1461  using_nonhilbert_norm = true;
1462 
1463  // Loop over all variables
1464  for (unsigned int var=0; var != this->n_vars(); ++var)
1465  {
1466  // Skip any variables we don't need to integrate
1467  Real norm_weight_sq = norm.weight_sq(var);
1468  if (norm_weight_sq == 0.0)
1469  continue;
1470  Real norm_weight = norm.weight(var);
1471 
1472  // Check for unimplemented norms (rather than just returning 0).
1473  FEMNormType norm_type = norm.type(var);
1474  if ((norm_type==H1) ||
1475  (norm_type==H2) ||
1476  (norm_type==L2) ||
1477  (norm_type==H1_SEMINORM) ||
1478  (norm_type==H2_SEMINORM))
1479  {
1480  if (!using_hilbert_norm)
1481  libmesh_not_implemented();
1482  using_nonhilbert_norm = false;
1483  }
1484  else if ((norm_type==L1) ||
1485  (norm_type==L_INF) ||
1486  (norm_type==W1_INF_SEMINORM) ||
1487  (norm_type==W2_INF_SEMINORM))
1488  {
1489  if (!using_nonhilbert_norm)
1490  libmesh_not_implemented();
1491  using_hilbert_norm = false;
1492  }
1493  else
1494  libmesh_not_implemented();
1495 
1496  const FEType & fe_type = this->get_dof_map().variable_type(var);
1497 
1498  // Allow space for dims 0-3, even if we don't use them all
1499  std::vector<FEBase *> fe_ptrs(4,libmesh_nullptr);
1500  std::vector<QBase *> q_rules(4,libmesh_nullptr);
1501 
1502  const std::set<unsigned char> & elem_dims = _mesh.elem_dimensions();
1503 
1504  // Prepare finite elements for each dimension present in the mesh
1505  for (std::set<unsigned char>::const_iterator d_it = elem_dims.begin();
1506  d_it != elem_dims.end(); ++d_it)
1507  {
1508  if (skip_dimensions && skip_dimensions->find(*d_it) != skip_dimensions->end())
1509  continue;
1510 
1511  q_rules[*d_it] =
1512  fe_type.default_quadrature_rule (*d_it).release();
1513 
1514  // Construct finite element object
1515 
1516  fe_ptrs[*d_it] = FEBase::build(*d_it, fe_type).release();
1517 
1518  // Attach quadrature rule to FE object
1519  fe_ptrs[*d_it]->attach_quadrature_rule (q_rules[*d_it]);
1520  }
1521 
1522  std::vector<dof_id_type> dof_indices;
1523 
1524  // Begin the loop over the elements
1525  for (const auto & elem : this->get_mesh().active_local_element_ptr_range())
1526  {
1527  const unsigned int dim = elem->dim();
1528 
1529 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1530 
1531  // One way for implementing this would be to exchange the fe with the FEInterface- class.
1532  // However, it needs to be discussed whether integral-norms make sense for infinite elements.
1533  // or in which sense they could make sense.
1534  if (elem->infinite() )
1535  libmesh_not_implemented();
1536 
1537 #endif
1538 
1539  if (skip_dimensions && skip_dimensions->find(dim) != skip_dimensions->end())
1540  continue;
1541 
1542  FEBase * fe = fe_ptrs[dim];
1543  QBase * qrule = q_rules[dim];
1544  libmesh_assert(fe);
1545  libmesh_assert(qrule);
1546 
1547  const std::vector<Real> & JxW = fe->get_JxW();
1548  const std::vector<std::vector<Real>> * phi = libmesh_nullptr;
1549  if (norm_type == H1 ||
1550  norm_type == H2 ||
1551  norm_type == L2 ||
1552  norm_type == L1 ||
1553  norm_type == L_INF)
1554  phi = &(fe->get_phi());
1555 
1556  const std::vector<std::vector<RealGradient>> * dphi = libmesh_nullptr;
1557  if (norm_type == H1 ||
1558  norm_type == H2 ||
1559  norm_type == H1_SEMINORM ||
1560  norm_type == W1_INF_SEMINORM)
1561  dphi = &(fe->get_dphi());
1562 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1563  const std::vector<std::vector<RealTensor>> * d2phi = libmesh_nullptr;
1564  if (norm_type == H2 ||
1565  norm_type == H2_SEMINORM ||
1566  norm_type == W2_INF_SEMINORM)
1567  d2phi = &(fe->get_d2phi());
1568 #endif
1569 
1570  fe->reinit (elem);
1571 
1572  this->get_dof_map().dof_indices (elem, dof_indices, var);
1573 
1574  const unsigned int n_qp = qrule->n_points();
1575 
1576  const unsigned int n_sf = cast_int<unsigned int>
1577  (dof_indices.size());
1578 
1579  // Begin the loop over the Quadrature points.
1580  for (unsigned int qp=0; qp<n_qp; qp++)
1581  {
1582  if (norm_type == L1)
1583  {
1584  Number u_h = 0.;
1585  for (unsigned int i=0; i != n_sf; ++i)
1586  u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1587  v_norm += norm_weight *
1588  JxW[qp] * std::abs(u_h);
1589  }
1590 
1591  if (norm_type == L_INF)
1592  {
1593  Number u_h = 0.;
1594  for (unsigned int i=0; i != n_sf; ++i)
1595  u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1596  v_norm = std::max(v_norm, norm_weight * std::abs(u_h));
1597  }
1598 
1599  if (norm_type == H1 ||
1600  norm_type == H2 ||
1601  norm_type == L2)
1602  {
1603  Number u_h = 0.;
1604  for (unsigned int i=0; i != n_sf; ++i)
1605  u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1606  v_norm += norm_weight_sq *
1607  JxW[qp] * TensorTools::norm_sq(u_h);
1608  }
1609 
1610  if (norm_type == H1 ||
1611  norm_type == H2 ||
1612  norm_type == H1_SEMINORM)
1613  {
1614  Gradient grad_u_h;
1615  for (unsigned int i=0; i != n_sf; ++i)
1616  grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
1617  v_norm += norm_weight_sq *
1618  JxW[qp] * grad_u_h.norm_sq();
1619  }
1620 
1621  if (norm_type == W1_INF_SEMINORM)
1622  {
1623  Gradient grad_u_h;
1624  for (unsigned int i=0; i != n_sf; ++i)
1625  grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
1626  v_norm = std::max(v_norm, norm_weight * grad_u_h.norm());
1627  }
1628 
1629 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1630  if (norm_type == H2 ||
1631  norm_type == H2_SEMINORM)
1632  {
1633  Tensor hess_u_h;
1634  for (unsigned int i=0; i != n_sf; ++i)
1635  hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
1636  v_norm += norm_weight_sq *
1637  JxW[qp] * hess_u_h.norm_sq();
1638  }
1639 
1640  if (norm_type == W2_INF_SEMINORM)
1641  {
1642  Tensor hess_u_h;
1643  for (unsigned int i=0; i != n_sf; ++i)
1644  hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
1645  v_norm = std::max(v_norm, norm_weight * hess_u_h.norm());
1646  }
1647 #endif
1648  }
1649  }
1650 
1651  // Need to delete the FE and quadrature objects to prevent a memory leak
1652  for (std::size_t i=0; i<fe_ptrs.size(); i++)
1653  if (fe_ptrs[i])
1654  delete fe_ptrs[i];
1655 
1656  for (std::size_t i=0; i<q_rules.size(); i++)
1657  if (q_rules[i])
1658  delete q_rules[i];
1659  }
1660 
1661  if (using_hilbert_norm)
1662  {
1663  this->comm().sum(v_norm);
1664  v_norm = std::sqrt(v_norm);
1665  }
1666  else
1667  {
1668  this->comm().max(v_norm);
1669  }
1670 
1671  return v_norm;
1672 }
double abs(double a)
virtual numeric_index_type size() const =0
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
void max(T &r) const
Take a local variable and replace it with the maximum of it&#39;s values on all processors.
void add_scaled(const TypeTensor< T2 > &, const T)
Add a scaled tensor to this tensor without creating a temporary.
Definition: type_tensor.h:777
void add_scaled(const TypeVector< T2 > &, const T)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:624
virtual Real l2_norm() const =0
unsigned int dim
const class libmesh_nullptr_t libmesh_nullptr
virtual Real l1_norm() const =0
long double max(long double a, double b)
UniquePtr< DofMap > _dof_map
Data structure describing the relationship between nodes, variables, etc...
Definition: system.h:1865
libmesh_assert(j)
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
const MeshBase & get_mesh() const
Definition: system.h:2014
const DofMap & get_dof_map() const
Definition: system.h:2030
NumberVectorValue Gradient
const std::set< unsigned char > & elem_dimensions() const
Definition: mesh_base.h:206
FEGenericBase< Real > FEBase
virtual Real linfty_norm() const =0
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
void sum(T &r) const
Take a local variable and replace it with the sum of it&#39;s values on all processors.
MeshBase & _mesh
Constant reference to the mesh data structure used for the simulation.
Definition: system.h:1877
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.
Real discrete_var_norm(const NumericVector< Number > &v, unsigned int var, FEMNormType norm_type) const
Finds the discrete norm for the entries in the vector corresponding to Dofs associated with var...
Definition: system.C:1364
FEMNormType
defines an enum for norms defined on vectors of finite element coefficients
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
void libMesh::RBConstruction::check_convergence ( LinearSolver< Number > &  input_solver)
protectedinherited

Check if the linear solver reports convergence.

Throw an error when that is not the case.

Definition at line 2159 of file rb_construction.C.

References libMesh::LinearSolver< T >::get_converged_reason().

Referenced by libMesh::RBConstruction::compute_Fq_representor_innerprods(), libMesh::RBConstruction::compute_output_dual_innerprods(), libMesh::RBConstruction::init_context(), libMesh::TransientRBConstruction::truth_solve(), libMesh::RBEIMConstruction::truth_solve(), libMesh::RBConstruction::truth_solve(), libMesh::TransientRBConstruction::update_residual_terms(), and libMesh::RBConstruction::update_residual_terms().

2160 {
2162 
2163  conv_flag = input_solver.get_converged_reason();
2164 
2165  if (conv_flag < 0)
2166  {
2167  std::stringstream err_msg;
2168  err_msg << "Convergence error. Error id: " << conv_flag;
2169  libmesh_error_msg(err_msg.str());
2170  }
2171 }
virtual LinearConvergenceReason get_converged_reason() const =0
LinearConvergenceReason
Linear solver convergence flags (taken from the PETSc flags)
void libMesh::RBConstruction::clear ( )
virtualinherited

Clear all the data structures associated with the system.

Reimplemented from libMesh::RBConstructionBase< LinearImplicitSystem >.

Reimplemented in libMesh::TransientSystem< RBConstruction >, libMesh::RBEIMConstruction, and libMesh::TransientRBConstruction.

Definition at line 96 of file rb_construction.C.

References libMesh::RBConstruction::Aq_vector, libMesh::RBConstructionBase< LinearImplicitSystem >::clear(), libMesh::RBConstruction::Fq_representor, libMesh::RBConstruction::Fq_representor_innerprods_computed, libMesh::RBConstruction::Fq_vector, libmesh_nullptr, libMesh::RBConstruction::non_dirichlet_Aq_vector, libMesh::RBConstruction::non_dirichlet_Fq_vector, libMesh::RBConstruction::non_dirichlet_outputs_vector, libMesh::RBConstruction::outputs_vector, and libMesh::RBConstruction::store_non_dirichlet_operators.

Referenced by libMesh::RBEIMConstruction::clear(), and libMesh::RBConstruction::~RBConstruction().

97 {
98  LOG_SCOPE("clear()", "RBConstruction");
99 
100  Parent::clear();
101 
102  for (std::size_t q=0; q<Aq_vector.size(); q++)
103  {
104  delete Aq_vector[q];
106  }
107 
108  for (std::size_t q=0; q<Fq_vector.size(); q++)
109  {
110  delete Fq_vector[q];
112  }
113 
114  for (std::size_t i=0; i<outputs_vector.size(); i++)
115  for (std::size_t q_l=0; q_l<outputs_vector[i].size(); q_l++)
116  {
117  delete outputs_vector[i][q_l];
119  }
120 
122  {
123  for (std::size_t q=0; q<non_dirichlet_Aq_vector.size(); q++)
124  {
125  delete non_dirichlet_Aq_vector[q];
127  }
128 
129  for (std::size_t q=0; q<non_dirichlet_Fq_vector.size(); q++)
130  {
131  delete non_dirichlet_Fq_vector[q];
133  }
134 
135  for (std::size_t i=0; i<non_dirichlet_outputs_vector.size(); i++)
136  for (std::size_t q_l=0; q_l<non_dirichlet_outputs_vector[i].size(); q_l++)
137  {
138  delete non_dirichlet_outputs_vector[i][q_l];
140  }
141  }
142 
143  // Also delete the Fq representors
144  for (std::size_t q_f=0; q_f<Fq_representor.size(); q_f++)
145  {
146  delete Fq_representor[q_f];
148  }
149  // Set Fq_representor_innerprods_computed flag to false now
150  // that we've cleared the Fq representors
152 }
std::vector< SparseMatrix< Number > * > non_dirichlet_Aq_vector
We may also need a second set of matrices/vectors that do not have the Dirichlet boundary conditions ...
const class libmesh_nullptr_t libmesh_nullptr
bool Fq_representor_innerprods_computed
A boolean flag to indicate whether or not the Fq representor norms have already been computed — used...
std::vector< std::vector< NumericVector< Number > * > > non_dirichlet_outputs_vector
std::vector< std::vector< NumericVector< Number > * > > outputs_vector
The libMesh vectors that define the output functionals.
std::vector< NumericVector< Number > * > Fq_representor
Vector storing the residual representors associated with the right-hand side.
std::vector< SparseMatrix< Number > * > Aq_vector
Vector storing the Q_a matrices from the affine expansion.
bool store_non_dirichlet_operators
Boolean flag to indicate whether we store a second copy of each affine operator and vector which does...
std::vector< NumericVector< Number > * > non_dirichlet_Fq_vector
std::vector< NumericVector< Number > * > Fq_vector
Vector storing the Q_f vectors in the affine decomposition of the right-hand side.
virtual void clear()
Clear all the data structures associated with the system.
const Parallel::Communicator& libMesh::ParallelObject::comm ( ) const
inherited
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_monitor(), libMesh::__libmesh_petsc_diff_solver_residual(), libMesh::__libmesh_petsc_preconditioner_apply(), libMesh::__libmesh_petsc_snes_jacobian(), 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(), add_cube_convex_hull_to_mesh(), libMesh::TransientRBConstruction::add_IC_to_RB_space(), libMesh::ImplicitSystem::add_matrix(), libMesh::RBConstruction::add_scaled_matrix_and_vector(), libMesh::System::add_vector(), libMesh::EigenSparseLinearSolver< T >::adjoint_solve(), libMesh::UnstructuredMesh::all_second_order(), libMesh::MeshTools::Modification::all_tri(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::TransientRBConstruction::allocate_data_structures(), libMesh::RBConstruction::allocate_data_structures(), libMesh::TransientRBConstruction::assemble_affine_expansion(), 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::Parallel::Sort< KeyType, IdxType >::communicate_bins(), libMesh::RBConstruction::compute_Fq_representor_innerprods(), libMesh::RBConstruction::compute_max_error_bound(), 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::RBConstruction::compute_output_dual_innerprods(), libMesh::RBSCMConstruction::compute_SCM_bounds_on_training_set(), 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::RBEIMConstruction::enrich_RB_space(), libMesh::TransientRBConstruction::enrich_RB_space(), libMesh::RBConstruction::enrich_RB_space(), AssembleOptimization::equality_constraints(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::RBEIMConstruction::evaluate_mesh_function(), 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::MeshfreeInterpolation::gather_remote_data(), libMesh::CondensedEigenSystem::get_eigenpair(), libMesh::DofMap::get_info(), libMesh::ImplicitSystem::get_linear_solver(), libMesh::EquationSystems::get_solution(), AssembleOptimization::inequality_constraints(), AssembleOptimization::inequality_constraints_jacobian(), libMesh::LocationMap< T >::init(), libMesh::TopologyMap::init(), libMesh::PetscDiffSolver::init(), libMesh::TimeSolver::init(), libMesh::TaoOptimizationSolver< T >::init(), libMesh::PetscNonlinearSolver< 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::RBEIMConstruction::initialize_rb_construction(), integrate_function(), 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::MeshRefinement::limit_level_mismatch_at_edge(), libMesh::MeshRefinement::limit_level_mismatch_at_node(), libMesh::MeshRefinement::limit_overrefined_boundary(), libMesh::MeshRefinement::limit_underrefined_boundary(), main(), 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::TransientRBConstruction::mass_matrix_scaled_matvec(), libMesh::FEMSystem::mesh_position_set(), libMesh::MeshSerializer::MeshSerializer(), LinearElasticityWithContact::move_mesh(), 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(), new_function_base(), 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::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(), FEMParameters::read(), libMesh::Nemesis_IO::read(), libMesh::CheckpointIO::read(), libMesh::XdrIO::read(), libMesh::CheckpointIO::read_header(), libMesh::XdrIO::read_header(), libMesh::RBEvaluation::read_in_vectors_from_multiple_files(), libMesh::TransientRBConstruction::read_riesz_representors_from_files(), libMesh::RBConstruction::read_riesz_representors_from_files(), 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(), LinearElasticityWithContact::residual_and_jacobian(), libMesh::MeshCommunication::send_coarse_ghosts(), libMesh::TransientRBConstruction::set_error_temporal_data(), libMesh::RBEIMConstruction::set_explicit_sys_subvector(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::Partitioner::set_parent_processor_ids(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::Parallel::Sort< KeyType, IdxType >::sort(), 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(), MeshfunctionDFEM::test_mesh_function_dfem(), MeshfunctionDFEM::test_mesh_function_dfem_grad(), libMesh::MeshRefinement::test_unflagged(), PointLocatorTest::testLocator(), BoundaryInfoTest::testMesh(), SystemsTest::testProjectCubeWithMeshFunction(), libMesh::MeshTools::total_weight(), libMesh::MeshFunctionSolutionTransfer::transfer(), libMesh::MeshfreeSolutionTransfer::transfer(), libMesh::TransientRBConstruction::truth_assembly(), libMesh::RBConstruction::truth_assembly(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::TransientRBConstruction::update_RB_initial_condition_all_N(), libMesh::RBEIMConstruction::update_RB_system_matrices(), libMesh::TransientRBConstruction::update_RB_system_matrices(), libMesh::RBConstruction::update_RB_system_matrices(), libMesh::TransientRBConstruction::update_residual_terms(), libMesh::RBConstruction::update_residual_terms(), libMesh::NameBasedIO::write(), libMesh::XdrIO::write(), libMesh::RBEvaluation::write_out_vectors(), libMesh::TransientRBConstruction::write_riesz_representors_to_files(), libMesh::RBConstruction::write_riesz_representors_to_files(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), libMesh::RBDataSerialization::RBEvaluationSerialization::write_to_file(), libMesh::RBDataSerialization::TransientRBEvaluationSerialization::write_to_file(), libMesh::RBDataSerialization::RBEIMEvaluationSerialization::write_to_file(), and libMesh::RBDataSerialization::RBSCMEvaluationSerialization::write_to_file().

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 531 of file system.C.

References libMesh::System::_is_initialized, libMesh::System::_sys_name, libMesh::System::_vectors, libMesh::System::get_vector(), libMesh::libmesh_assert(), 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().

534 {
535  // we do not care for matrices, but for vectors
537  libmesh_assert (other_system._is_initialized);
538 
539  if (verbose)
540  {
541  libMesh::out << " Systems \"" << _sys_name << "\"" << std::endl;
542  libMesh::out << " comparing matrices not supported." << std::endl;
543  libMesh::out << " comparing names...";
544  }
545 
546  // compare the name: 0 means identical
547  const int name_result = _sys_name.compare(other_system.name());
548  if (verbose)
549  {
550  if (name_result == 0)
551  libMesh::out << " identical." << std::endl;
552  else
553  libMesh::out << " names not identical." << std::endl;
554  libMesh::out << " comparing solution vector...";
555  }
556 
557 
558  // compare the solution: -1 means identical
559  const int solu_result = solution->compare (*other_system.solution.get(),
560  threshold);
561 
562  if (verbose)
563  {
564  if (solu_result == -1)
565  libMesh::out << " identical up to threshold." << std::endl;
566  else
567  libMesh::out << " first difference occurred at index = "
568  << solu_result << "." << std::endl;
569  }
570 
571 
572  // safety check, whether we handle at least the same number
573  // of vectors
574  std::vector<int> ov_result;
575 
576  if (this->n_vectors() != other_system.n_vectors())
577  {
578  if (verbose)
579  {
580  libMesh::out << " Fatal difference. This system handles "
581  << this->n_vectors() << " add'l vectors," << std::endl
582  << " while the other system handles "
583  << other_system.n_vectors()
584  << " add'l vectors." << std::endl
585  << " Aborting comparison." << std::endl;
586  }
587  return false;
588  }
589  else if (this->n_vectors() == 0)
590  {
591  // there are no additional vectors...
592  ov_result.clear ();
593  }
594  else
595  {
596  // compare other vectors
597  for (const_vectors_iterator pos = _vectors.begin();
598  pos != _vectors.end(); ++pos)
599  {
600  if (verbose)
601  libMesh::out << " comparing vector \""
602  << pos->first << "\" ...";
603 
604  // assume they have the same name
605  const NumericVector<Number> & other_system_vector =
606  other_system.get_vector(pos->first);
607 
608  ov_result.push_back(pos->second->compare (other_system_vector,
609  threshold));
610 
611  if (verbose)
612  {
613  if (ov_result[ov_result.size()-1] == -1)
614  libMesh::out << " identical up to threshold." << std::endl;
615  else
616  libMesh::out << " first difference occurred at" << std::endl
617  << " index = " << ov_result[ov_result.size()-1] << "." << std::endl;
618  }
619 
620  }
621 
622  } // finished comparing additional vectors
623 
624 
625  bool overall_result;
626 
627  // sum up the results
628  if ((name_result==0) && (solu_result==-1))
629  {
630  if (ov_result.size()==0)
631  overall_result = true;
632  else
633  {
634  bool ov_identical;
635  unsigned int n = 0;
636  do
637  {
638  ov_identical = (ov_result[n]==-1);
639  n++;
640  }
641  while (ov_identical && n<ov_result.size());
642  overall_result = ov_identical;
643  }
644  }
645  else
646  overall_result = false;
647 
648  if (verbose)
649  {
650  libMesh::out << " finished comparisons, ";
651  if (overall_result)
652  libMesh::out << "found no differences." << std::endl << std::endl;
653  else
654  libMesh::out << "found differences." << std::endl << std::endl;
655  }
656 
657  return overall_result;
658 }
bool _is_initialized
true when additional vectors and variables do not require immediate initialization, false otherwise.
Definition: system.h:1952
std::map< std::string, NumericVector< Number > * >::const_iterator const_vectors_iterator
Definition: system.h:749
libmesh_assert(j)
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
std::map< std::string, NumericVector< Number > * > _vectors
Some systems need an arbitrary number of vectors.
Definition: system.h:1916
OStreamProxy out
const std::string _sys_name
A name associated with this system.
Definition: system.h:1882
unsigned int n_vectors() const
Definition: system.h:2214
void libMesh::RBConstruction::compute_Fq_representor_innerprods ( bool  compute_inner_products = true)
protectedvirtualinherited

Compute the terms that are combined `online' to determine the dual norm of the residual.

Here we compute the terms associated with the right-hand side. These terms are basis independent, hence we separate them from the rest of the calculations that are done in update_residual_terms. By default, inner product terms are also computed, but you can turn this feature off e.g. if you are already reading in that data from files.

Definition at line 1688 of file rb_construction.C.

References libMesh::NumericVector< T >::add(), libMesh::RBConstruction::assert_convergence, libMesh::NumericVector< T >::build(), libMesh::RBConstruction::check_convergence(), libMesh::ParallelObject::comm(), libMesh::LinearImplicitSystem::final_linear_residual(), libMesh::RBConstruction::Fq_representor, libMesh::RBEvaluation::Fq_representor_innerprods, libMesh::RBConstruction::Fq_representor_innerprods, libMesh::RBConstruction::Fq_representor_innerprods_computed, libMesh::RBConstruction::get_Fq(), libMesh::RBThetaExpansion::get_n_F_terms(), libMesh::RBConstruction::get_rb_evaluation(), libMesh::RBConstruction::get_rb_theta_expansion(), libMesh::Utility::get_timestamp(), libMesh::RBConstruction::inner_product_matrix, libMesh::RBConstruction::inner_product_solver, libMesh::RBConstructionBase< LinearImplicitSystem >::inner_product_storage_vector, libMesh::RBConstruction::is_quiet(), libMesh::libmesh_assert(), libMesh::System::n_dofs(), libMesh::LinearImplicitSystem::n_linear_iterations(), libMesh::System::n_local_dofs(), libMesh::out, libMesh::PARALLEL, libMesh::ExplicitSystem::rhs, libMesh::System::solution, libMesh::RBConstruction::solve_for_matrix_and_rhs(), and libMesh::NumericVector< T >::zero().

Referenced by libMesh::RBConstruction::recompute_all_residual_terms(), and libMesh::RBConstruction::train_reduced_basis().

1689 {
1690 
1691  // Skip calculations if we've already computed the Fq_representors
1693  {
1694  // Only log if we get to here
1695  LOG_SCOPE("compute_Fq_representor_innerprods()", "RBConstruction");
1696 
1697  for (unsigned int q_f=0; q_f<get_rb_theta_expansion().get_n_F_terms(); q_f++)
1698  {
1699  if (!Fq_representor[q_f])
1700  {
1701  Fq_representor[q_f] = (NumericVector<Number>::build(this->comm()).release());
1702  Fq_representor[q_f]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1703  }
1704 
1705  libmesh_assert(Fq_representor[q_f]->size() == this->n_dofs() &&
1706  Fq_representor[q_f]->local_size() == this->n_local_dofs() );
1707 
1708  rhs->zero();
1709  rhs->add(1., *get_Fq(q_f));
1710 
1711  if (!is_quiet())
1712  libMesh::out << "Starting solve q_f=" << q_f
1713  << " in RBConstruction::update_residual_terms() at "
1714  << Utility::get_timestamp() << std::endl;
1715 
1717 
1718  if (assert_convergence)
1720 
1721  if (!is_quiet())
1722  {
1723  libMesh::out << "Finished solve q_f=" << q_f
1724  << " in RBConstruction::update_residual_terms() at "
1725  << Utility::get_timestamp() << std::endl;
1726 
1728  << " iterations, final residual "
1729  << this->final_linear_residual() << std::endl;
1730  }
1731 
1732  *Fq_representor[q_f] = *solution;
1733  }
1734 
1735  if (compute_inner_products)
1736  {
1737  unsigned int q=0;
1738 
1739  for (unsigned int q_f1=0; q_f1<get_rb_theta_expansion().get_n_F_terms(); q_f1++)
1740  {
1742 
1743  for (unsigned int q_f2=q_f1; q_f2<get_rb_theta_expansion().get_n_F_terms(); q_f2++)
1744  {
1746 
1747  q++;
1748  }
1749  }
1750  } // end if (compute_inner_products)
1751 
1753  }
1754 
1756 }
unsigned int n_linear_iterations() const
std::vector< Number > Fq_representor_innerprods
Vectors storing the residual representor inner products to be used in computing the residuals online...
UniquePtr< NumericVector< Number > > inner_product_storage_vector
We keep an extra temporary vector that is useful for performing inner products (avoids unnecessary me...
bool assert_convergence
A boolean flag to indicate whether to check for proper convergence after each solve.
std::string get_timestamp()
Definition: timestamp.C:37
NumericVector< Number > * rhs
The system matrix.
NumericVector< Number > * get_Fq(unsigned int q)
Get a pointer to Fq.
std::vector< Number > Fq_representor_innerprods
Vectors storing the residual representor inner products to be used in computing the residuals online...
UniquePtr< SparseMatrix< Number > > inner_product_matrix
The inner product matrix.
UniquePtr< LinearSolver< Number > > inner_product_solver
We store an extra linear solver object which we can optionally use for solving all systems in which t...
virtual void zero()=0
Set all entries to zero.
libmesh_assert(j)
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
virtual void solve_for_matrix_and_rhs(LinearSolver< Number > &input_solver, SparseMatrix< Number > &input_matrix, NumericVector< Number > &input_rhs)
Assembles & solves the linear system A*x=b for the specified matrix input_matrix and right-hand side ...
bool Fq_representor_innerprods_computed
A boolean flag to indicate whether or not the Fq representor norms have already been computed — used...
RBThetaExpansion & get_rb_theta_expansion()
Get a reference to the RBThetaExpansion object that that belongs to rb_eval.
unsigned int get_n_F_terms() const
Get Q_f, the number of terms in the affine expansion for the right-hand side.
std::vector< NumericVector< Number > * > Fq_representor
Vector storing the residual representors associated with the right-hand side.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
dof_id_type n_local_dofs() const
Definition: system.C:185
bool is_quiet() const
Is the system in quiet mode?
void check_convergence(LinearSolver< Number > &input_solver)
Check if the linear solver reports convergence.
const Parallel::Communicator & comm() const
OStreamProxy out
dof_id_type n_dofs() const
Definition: system.C:148
RBEvaluation & get_rb_evaluation()
Get a reference to the RBEvaluation object.
virtual void add(const numeric_index_type i, const T value)=0
Adds value to each entry of the vector.
Real libMesh::RBConstruction::compute_max_error_bound ( )
virtualinherited

(i) Compute the a posteriori error bound for each set of parameters in the training set, (ii) set current_parameters to the parameters that maximize the error bound, and (iii) return the maximum error bound.

Definition at line 1345 of file rb_construction.C.

References libMesh::RBConstructionBase< LinearImplicitSystem >::broadcast_parameters(), libMesh::ParallelObject::comm(), libMesh::RBConstructionBase< LinearImplicitSystem >::get_first_local_training_index(), libMesh::RBConstructionBase< LinearImplicitSystem >::get_global_max_error_pair(), libMesh::RBConstructionBase< LinearImplicitSystem >::get_last_local_training_index(), libMesh::RBConstructionBase< LinearImplicitSystem >::get_local_n_training_samples(), libMesh::RBParametrized::get_n_params(), libMesh::RBConstruction::get_RB_error_bound(), libMesh::RBConstruction::get_rb_evaluation(), std::max(), libMesh::ParallelObject::processor_id(), libMesh::Real, libMesh::RBConstructionBase< LinearImplicitSystem >::serial_training_set, libMesh::RBConstructionBase< LinearImplicitSystem >::set_params_from_training_set(), libMesh::Parallel::Communicator::sum(), and libMesh::RBConstruction::training_error_bounds.

Referenced by libMesh::RBConstruction::train_reduced_basis().

1346 {
1347  LOG_SCOPE("compute_max_error_bound()", "RBConstruction");
1348 
1349  // Treat the case with no parameters in a special way
1350  if (get_n_params() == 0)
1351  {
1352  Real max_val;
1353  if (std::numeric_limits<Real>::has_infinity)
1354  {
1355  max_val = std::numeric_limits<Real>::infinity();
1356  }
1357  else
1358  {
1359  max_val = std::numeric_limits<Real>::max();
1360  }
1361 
1362  // Make sure we do at least one solve, but otherwise return a zero error bound
1363  // when we have no parameters
1364  return (get_rb_evaluation().get_n_basis_functions() == 0) ? max_val : 0.;
1365  }
1366 
1368 
1369  // keep track of the maximum error
1370  unsigned int max_err_index = 0;
1371  Real max_err = 0.;
1372 
1374  for (unsigned int i=0; i<get_local_n_training_samples(); i++)
1375  {
1376  // Load training parameter i, this is only loaded
1377  // locally since the RB solves are local.
1378  set_params_from_training_set( first_index+i );
1379 
1381 
1382  if (training_error_bounds[i] > max_err)
1383  {
1384  max_err_index = i;
1385  max_err = training_error_bounds[i];
1386  }
1387  }
1388 
1389  std::pair<numeric_index_type, Real> error_pair(first_index+max_err_index, max_err);
1390  get_global_max_error_pair(this->comm(),error_pair);
1391 
1392  // If we have a serial training set (i.e. a training set that is the same on all processors)
1393  // just set the parameters on all processors
1394  if (serial_training_set)
1395  {
1396  set_params_from_training_set( error_pair.first );
1397  }
1398  // otherwise, broadcast the parameter that produced the maximum error
1399  else
1400  {
1401  unsigned int root_id=0;
1402  if ((get_first_local_training_index() <= error_pair.first) &&
1403  (error_pair.first < get_last_local_training_index()))
1404  {
1405  set_params_from_training_set( error_pair.first );
1406  root_id = this->processor_id();
1407  }
1408 
1409  this->comm().sum(root_id); // root_id is only non-zero on one processor
1410  broadcast_parameters(root_id);
1411  }
1412 
1413  return error_pair.second;
1414 }
numeric_index_type get_local_n_training_samples() const
Get the total number of training samples local to this processor.
std::vector< Real > training_error_bounds
Vector storing the values of the error bound for each parameter in the training set — the parameter ...
void broadcast_parameters(unsigned int proc_id)
Broadcasts parameters on processor proc_id to all processors.
unsigned int get_n_params() const
Get the number of parameters.
static void get_global_max_error_pair(const Parallel::Communicator &communicator, std::pair< numeric_index_type, Real > &error_pair)
Static function to return the error pair (index,error) that is corresponds to the largest error on al...
long double max(long double a, double b)
dof_id_type numeric_index_type
Definition: id_types.h:92
void set_params_from_training_set(unsigned int index)
Set parameters to the RBParameters stored in index index of the training set.
virtual Real get_RB_error_bound()
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
numeric_index_type get_first_local_training_index() const
Get the first local index of the training parameters.
const Parallel::Communicator & comm() const
numeric_index_type get_last_local_training_index() const
Get the last local index of the training parameters.
bool serial_training_set
This boolean flag indicates whether or not the training set should be the same on all processors...
RBEvaluation & get_rb_evaluation()
Get a reference to the RBEvaluation object.
void sum(T &r) const
Take a local variable and replace it with the sum of it&#39;s values on all processors.
processor_id_type processor_id() const
void libMesh::RBConstruction::compute_output_dual_innerprods ( )
protectedvirtualinherited

Compute and store the dual norm of each output functional.

Definition at line 1583 of file rb_construction.C.

References libMesh::NumericVector< T >::add(), libMesh::RBConstruction::assert_convergence, libMesh::NumericVector< T >::build(), libMesh::RBConstruction::check_convergence(), libMesh::ParallelObject::comm(), libMesh::LinearImplicitSystem::final_linear_residual(), libMesh::LinearImplicitSystem::get_linear_solver(), libMesh::RBConstruction::get_matrix_for_output_dual_solves(), libMesh::RBThetaExpansion::get_n_output_terms(), libMesh::RBThetaExpansion::get_n_outputs(), libMesh::RBConstruction::get_output_vector(), libMesh::RBConstruction::get_rb_evaluation(), libMesh::RBConstruction::get_rb_theta_expansion(), libMesh::Utility::get_timestamp(), libMesh::RBConstructionBase< LinearImplicitSystem >::inner_product_storage_vector, libMesh::RBConstruction::is_quiet(), libmesh_nullptr, libMesh::LinearImplicitSystem::linear_solver, libMesh::System::n_dofs(), libMesh::LinearImplicitSystem::n_linear_iterations(), libMesh::System::n_local_dofs(), libMesh::out, libMesh::RBEvaluation::output_dual_innerprods, libMesh::RBConstruction::output_dual_innerprods, libMesh::RBConstruction::output_dual_innerprods_computed, libMesh::PARALLEL, libMesh::ExplicitSystem::rhs, libMesh::System::solution, libMesh::RBConstruction::solve_for_matrix_and_rhs(), libMesh::SparseMatrix< T >::vector_mult(), and libMesh::NumericVector< T >::zero().

Referenced by libMesh::RBConstruction::train_reduced_basis().

1584 {
1585  // Skip calculations if we've already computed the output dual norms
1587  {
1588  // Short circuit if we don't have any outputs
1589  if (get_rb_theta_expansion().get_n_outputs() == 0)
1590  {
1592  return;
1593  }
1594 
1595  // Only log if we get to here
1596  LOG_SCOPE("compute_output_dual_innerprods()", "RBConstruction");
1597 
1598  libMesh::out << "Compute output dual inner products" << std::endl;
1599 
1600  // Find out the largest value of Q_l
1601  unsigned int max_Q_l = 0;
1602  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1603  max_Q_l = (get_rb_theta_expansion().get_n_output_terms(n) > max_Q_l) ? get_rb_theta_expansion().get_n_output_terms(n) : max_Q_l;
1604 
1605  std::vector<NumericVector<Number> *> L_q_representor(max_Q_l);
1606  for (unsigned int q=0; q<max_Q_l; q++)
1607  {
1608  L_q_representor[q] = (NumericVector<Number>::build(this->comm()).release());
1609  L_q_representor[q]->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1610  }
1611 
1612  for (unsigned int n=0; n<get_rb_theta_expansion().get_n_outputs(); n++)
1613  {
1614  for (unsigned int q_l=0; q_l<get_rb_theta_expansion().get_n_output_terms(n); q_l++)
1615  {
1616  rhs->zero();
1617  rhs->add(1., *get_output_vector(n,q_l));
1618 
1619  if (!is_quiet())
1620  libMesh::out << "Starting solve n=" << n << ", q_l=" << q_l
1621  << " in RBConstruction::compute_output_dual_innerprods() at "
1622  << Utility::get_timestamp() << std::endl;
1623 
1624  // Use the main linear solver here instead of the inner_product solver, since
1625  // get_matrix_for_output_dual_solves() may not return the inner product matrix.
1627 
1628  // We possibly perform multiple solves here with the same matrix, hence
1629  // set reuse_preconditioner(true) (and set it back to false again below
1630  // at the end of this function).
1631  linear_solver->reuse_preconditioner(true);
1632 
1633  if (assert_convergence)
1635 
1636  if (!is_quiet())
1637  {
1638  libMesh::out << "Finished solve n=" << n << ", q_l=" << q_l
1639  << " in RBConstruction::compute_output_dual_innerprods() at "
1640  << Utility::get_timestamp() << std::endl;
1641 
1643  << " iterations, final residual "
1644  << this->final_linear_residual() << std::endl;
1645  }
1646 
1647  *L_q_representor[q_l] = *solution;
1648  }
1649 
1650  unsigned int q=0;
1651  for (unsigned int q_l1=0; q_l1<get_rb_theta_expansion().get_n_output_terms(n); q_l1++)
1652  {
1654 
1655  for (unsigned int q_l2=q_l1; q_l2<get_rb_theta_expansion().get_n_output_terms(n); q_l2++)
1656  {
1657  output_dual_innerprods[n][q] = L_q_representor[q_l2]->dot(*inner_product_storage_vector);
1658  libMesh::out << "output_dual_innerprods[" << n << "][" << q << "] = " << output_dual_innerprods[n][q] << std::endl;
1659 
1660  q++;
1661  }
1662  }
1663  }
1664 
1665  // Finally clear the L_q_representor vectors
1666  for (unsigned int q=0; q<max_Q_l; q++)
1667  {
1668  if (L_q_representor[q])
1669  {
1670  delete L_q_representor[q];
1671  L_q_representor[q] = libmesh_nullptr;
1672  }
1673  }
1674 
1675  // We may not need to use linear_solver again (e.g. this would happen if we use
1676  // extra_linear_solver for the truth_solves). As a result, let's clear linear_solver
1677  // to release any memory it may be taking up. If we do need it again, it will
1678  // be initialized when necessary.
1679  linear_solver->clear();
1680  linear_solver->reuse_preconditioner(false);
1681 
1683  }
1684 
1686 }
unsigned int n_linear_iterations() const
UniquePtr< NumericVector< Number > > inner_product_storage_vector
We keep an extra temporary vector that is useful for performing inner products (avoids unnecessary me...
bool assert_convergence
A boolean flag to indicate whether to check for proper convergence after each solve.
std::string get_timestamp()
Definition: timestamp.C:37
unsigned int get_n_output_terms(unsigned int output_index) const
Get the number of affine terms associated with the specified output.
std::vector< std::vector< Number > > output_dual_innerprods
The vector storing the dual norm inner product terms for each output.
void vector_mult(NumericVector< T > &dest, const NumericVector< T > &arg) const
Multiplies the matrix by the NumericVector arg and stores the result in NumericVector dest...
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
std::vector< std::vector< Number > > output_dual_innerprods
The vector storing the dual norm inner product terms for each output.
virtual void zero()=0
Set all entries to zero.
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
virtual void solve_for_matrix_and_rhs(LinearSolver< Number > &input_solver, SparseMatrix< Number > &input_matrix, NumericVector< Number > &input_rhs)
Assembles & solves the linear system A*x=b for the specified matrix input_matrix and right-hand side ...
RBThetaExpansion & get_rb_theta_expansion()
Get a reference to the RBThetaExpansion object that that belongs to rb_eval.
virtual LinearSolver< Number > * get_linear_solver() const libmesh_override
virtual SparseMatrix< Number > & get_matrix_for_output_dual_solves()
Return the matrix for the output residual dual norm solves.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
dof_id_type n_local_dofs() const
Definition: system.C:185
unsigned int get_n_outputs() const
Get n_outputs, the number output functionals.
UniquePtr< LinearSolver< Number > > linear_solver
The LinearSolver defines the interface used to solve the linear_implicit system.
bool is_quiet() const
Is the system in quiet mode?
void check_convergence(LinearSolver< Number > &input_solver)
Check if the linear solver reports convergence.
const Parallel::Communicator & comm() const
OStreamProxy out
NumericVector< Number > * get_output_vector(unsigned int n, unsigned int q_l)
Get a pointer to the n^th output.
dof_id_type n_dofs() const
Definition: system.C:148
RBEvaluation & get_rb_evaluation()
Get a reference to the RBEvaluation object.
virtual void add(const numeric_index_type i, const T value)=0
Adds value to each entry of the vector.
bool output_dual_innerprods_computed
A boolean flag to indicate whether or not the output dual norms have already been computed — used to...
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(), compute_stresses(), LinearElasticityWithContact::compute_stresses(), LinearElasticity::compute_stresses(), LargeDeformationElasticity::compute_stresses(), libMesh::ExactErrorEstimator::estimate_error(), main(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::HPCoarsenTest::select_refinement(), SolidSystem::side_time_derivative(), 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 }
UniquePtr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1535
UniquePtr< DofMap > _dof_map
Data structure describing the relationship between nodes, variables, etc...
Definition: system.h:1865
void libMesh::System::deactivate ( )
inherited

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 }
bool _active
Flag stating if the system is active or not.
Definition: system.h:1908
void libMesh::LinearImplicitSystem::detach_shell_matrix ( )
inherited

Detaches a shell matrix.

Same as attach_shell_matrix(libmesh_nullptr).

Definition at line 179 of file linear_implicit_system.h.

References libMesh::LinearImplicitSystem::attach_shell_matrix(), and libmesh_nullptr.

Referenced by main().

void attach_shell_matrix(ShellMatrix< Number > *shell_matrix)
This function enables the user to provide a shell matrix, i.e.
const class libmesh_nullptr_t libmesh_nullptr
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)
Set the same_preconditioner flag, which indicates if we reuse the same preconditioner for subsequent ...
bool assemble_before_solve
Flag which tells the system to whether or not to call the user assembly function during each call to ...
Definition: system.h:1477
void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited

Definition at line 107 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

Referenced by libMesh::LibMeshInit::LibMeshInit(), and libMesh::ReferenceCounter::n_objects().

108 {
109  _enable_print_counter = false;
110  return;
111 }
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called...
void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited

Definition at line 107 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

Referenced by libMesh::LibMeshInit::LibMeshInit(), and libMesh::ReferenceCounter::n_objects().

108 {
109  _enable_print_counter = false;
110  return;
111 }
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called...
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 }
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called...
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 }
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called...
void libMesh::RBConstruction::enrich_RB_space ( )
protectedvirtualinherited

Add a new basis function to the RB space.

This is called during train_reduced_basis.

Reimplemented in libMesh::TransientRBConstruction, and libMesh::RBEIMConstruction.

Definition at line 1262 of file rb_construction.C.

References libMesh::NumericVector< T >::add(), libMesh::RBEvaluation::basis_functions, libMesh::NumericVector< T >::build(), libMesh::ParallelObject::comm(), libMesh::RBEvaluation::get_n_basis_functions(), libMesh::RBConstruction::get_rb_evaluation(), libMesh::NumericVector< T >::init(), libMesh::RBConstruction::inner_product_matrix, libMesh::RBConstructionBase< LinearImplicitSystem >::inner_product_storage_vector, libMesh::System::n_dofs(), libMesh::System::n_local_dofs(), libMesh::PARALLEL, libMesh::NumericVector< T >::scale(), libMesh::System::solution, and libMesh::NumericVector< T >::zero().

Referenced by libMesh::RBConstruction::train_reduced_basis().

1263 {
1264  LOG_SCOPE("enrich_RB_space()", "RBConstruction");
1265 
1266  NumericVector<Number> * new_bf = NumericVector<Number>::build(this->comm()).release();
1267  new_bf->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
1268  *new_bf = *solution;
1269 
1270  for (unsigned int index=0; index<get_rb_evaluation().get_n_basis_functions(); index++)
1271  {
1272  inner_product_matrix->vector_mult(*inner_product_storage_vector, *new_bf);
1273 
1274  Number scalar =
1275  inner_product_storage_vector->dot(get_rb_evaluation().get_basis_function(index));
1276  new_bf->add(-scalar, get_rb_evaluation().get_basis_function(index));
1277  }
1278 
1279  // Normalize new_bf
1280  inner_product_matrix->vector_mult(*inner_product_storage_vector, *new_bf);
1281  Number new_bf_norm = std::sqrt( inner_product_storage_vector->dot(*new_bf) );
1282 
1283  if (new_bf_norm == 0.)
1284  {
1285  new_bf->zero(); // avoid potential nan's
1286  }
1287  else
1288  {
1289  new_bf->scale(1./new_bf_norm);
1290  }
1291 
1292  // load the new basis function into the basis_functions vector.
1293  get_rb_evaluation().basis_functions.push_back( new_bf );
1294 }
std::vector< NumericVector< Number > * > basis_functions
The libMesh vectors storing the finite element coefficients of the RB basis functions.
UniquePtr< NumericVector< Number > > inner_product_storage_vector
We keep an extra temporary vector that is useful for performing inner products (avoids unnecessary me...
UniquePtr< SparseMatrix< Number > > inner_product_matrix
The inner product matrix.
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
dof_id_type n_local_dofs() const
Definition: system.C:185
const Parallel::Communicator & comm() const
dof_id_type n_dofs() const
Definition: system.C:148
RBEvaluation & get_rb_evaluation()
Get a reference to the RBEvaluation object.
virtual unsigned int get_n_basis_functions() const
Get the current number of basis functions.
Real libMesh::LinearImplicitSystem::final_linear_residual ( ) const
inherited
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 819 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().

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