libMesh::ContinuationSystem Class Reference

#include <continuation_system.h>

Inheritance diagram for libMesh::ContinuationSystem:

Public Types

enum  Predictor { Euler, AB2, Invalid_Predictor }
 
typedef ContinuationSystem sys_type
 
typedef FEMSystem Parent
 
typedef bool(TimeSolver::* TimeSolverResPtr )(bool, DiffContext &)
 
typedef std::map< std::string,
SparseMatrix< Number >
* >::iterator 
matrices_iterator
 
typedef std::map< std::string,
SparseMatrix< Number >
* >::const_iterator 
const_matrices_iterator
 
typedef std::map< std::string,
NumericVector< Number >
* >::iterator 
vectors_iterator
 
typedef std::map< std::string,
NumericVector< Number >
* >::const_iterator 
const_vectors_iterator
 

Public Member Functions

 ContinuationSystem (EquationSystems &es, const std::string &name, const unsigned int number)
 
virtual ~ContinuationSystem ()
 
virtual void clear ()
 
virtual void solve ()
 
void continuation_solve ()
 
void advance_arcstep ()
 
void set_max_arclength_stepsize (Real maxds)
 
void save_current_solution ()
 
virtual void assembly (bool get_residual, bool get_jacobian)
 
void mesh_position_get ()
 
void mesh_position_set ()
 
virtual AutoPtr< DiffContextbuild_context ()
 
virtual void init_context (DiffContext &)
 
virtual void postprocess ()
 
virtual void assemble_qoi (const QoISet &indices=QoISet())
 
virtual void assemble_qoi_derivative (const QoISet &indices=QoISet())
 
void numerical_jacobian (TimeSolverResPtr res, FEMContext &context) const
 
void numerical_elem_jacobian (FEMContext &context) const
 
void numerical_side_jacobian (FEMContext &context) const
 
virtual void reinit ()
 
virtual void assemble ()
 
virtual LinearSolver< Number > * get_linear_solver () const
 
virtual std::pair< unsigned
int, Real
get_linear_solve_parameters () const
 
virtual void release_linear_solver (LinearSolver< Number > *) const
 
virtual std::pair< unsigned
int, Real
adjoint_solve (const QoISet &qoi_indices=QoISet())
 
virtual AutoPtr
< DifferentiablePhysics
clone_physics ()
 
virtual AutoPtr
< DifferentiableQoI
clone ()
 
const DifferentiablePhysicsget_physics () const
 
DifferentiablePhysicsget_physics ()
 
void attach_physics (DifferentiablePhysics *physics_in)
 
const DifferentiableQoIget_qoi () const
 
DifferentiableQoIget_qoi ()
 
void attach_qoi (DifferentiableQoI *qoi_in)
 
void set_time_solver (AutoPtr< TimeSolver > _time_solver)
 
TimeSolverget_time_solver ()
 
const TimeSolverget_time_solver () const
 
virtual void element_postprocess (DiffContext &)
 
virtual void side_postprocess (DiffContext &)
 
sys_typesystem ()
 
virtual void disable_cache ()
 
virtual std::string system_type () const
 
virtual void assemble_residual_derivatives (const ParameterVector &parameters)
 
virtual std::pair< unsigned
int, Real
sensitivity_solve (const ParameterVector &parameters)
 
virtual std::pair< unsigned
int, Real
weighted_sensitivity_solve (const ParameterVector &parameters, const ParameterVector &weights)
 
virtual std::pair< unsigned
int, Real
weighted_sensitivity_adjoint_solve (const ParameterVector &parameters, const ParameterVector &weights, const QoISet &qoi_indices=QoISet())
 
virtual void adjoint_qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
 
virtual void forward_qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
 
virtual void qoi_parameter_hessian (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &hessian)
 
virtual void qoi_parameter_hessian_vector_product (const QoISet &qoi_indices, const ParameterVector &parameters, const ParameterVector &vector, SensitivityData &product)
 
SparseMatrix< Number > & add_matrix (const std::string &mat_name)
 
bool have_matrix (const std::string &mat_name) const
 
const SparseMatrix< Number > * request_matrix (const std::string &mat_name) const
 
SparseMatrix< Number > * request_matrix (const std::string &mat_name)
 
const SparseMatrix< Number > & get_matrix (const std::string &mat_name) const
 
SparseMatrix< Number > & get_matrix (const std::string &mat_name)
 
virtual unsigned int n_matrices () const
 
void init ()
 
virtual void update ()
 
virtual void restrict_solve_to (const SystemSubset *subset, const SubsetSolveMode subset_solve_mode=SUBSET_ZERO)
 
bool is_adjoint_already_solved () const
 
void set_adjoint_already_solved (bool setting)
 
virtual void qoi_parameter_sensitivity (const QoISet &qoi_indices, const ParameterVector &parameters, SensitivityData &sensitivities)
 
virtual bool compare (const System &other_system, const Real threshold, const bool verbose) const
 
const std::string & name () const
 
void project_solution (FunctionBase< Number > *f, FunctionBase< Gradient > *g=NULL) const
 
void project_solution (FEMFunctionBase< Number > *f, FEMFunctionBase< Gradient > *g=NULL) const
 
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
 
void project_vector (NumericVector< Number > &new_vector, FunctionBase< Number > *f, FunctionBase< Gradient > *g=NULL) const
 
void project_vector (NumericVector< Number > &new_vector, FEMFunctionBase< Number > *f, FEMFunctionBase< Gradient > *g=NULL) const
 
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) const
 
void boundary_project_solution (const std::set< boundary_id_type > &b, const std::vector< unsigned int > &variables, FunctionBase< Number > *f, FunctionBase< Gradient > *g=NULL)
 
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)
 
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=NULL) const
 
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) const
 
unsigned int number () const
 
void update_global_solution (std::vector< Number > &global_soln) const
 
void update_global_solution (std::vector< Number > &global_soln, const unsigned int dest_proc) const
 
const MeshBaseget_mesh () const
 
MeshBaseget_mesh ()
 
const DofMapget_dof_map () const
 
DofMapget_dof_map ()
 
const EquationSystemsget_equation_systems () const
 
EquationSystemsget_equation_systems ()
 
bool active () const
 
void activate ()
 
void deactivate ()
 
void set_basic_system_only ()
 
vectors_iterator vectors_begin ()
 
const_vectors_iterator vectors_begin () const
 
vectors_iterator vectors_end ()
 
const_vectors_iterator vectors_end () const
 
NumericVector< Number > & add_vector (const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
 
void remove_vector (const std::string &vec_name)
 
bool & project_solution_on_reinit (void)
 
bool have_vector (const std::string &vec_name) const
 
const NumericVector< Number > * request_vector (const std::string &vec_name) const
 
NumericVector< Number > * request_vector (const std::string &vec_name)
 
const NumericVector< Number > * request_vector (const unsigned int vec_num) const
 
NumericVector< Number > * request_vector (const unsigned int vec_num)
 
const NumericVector< Number > & get_vector (const std::string &vec_name) const
 
NumericVector< Number > & get_vector (const std::string &vec_name)
 
const NumericVector< Number > & get_vector (const unsigned int vec_num) const
 
NumericVector< Number > & get_vector (const unsigned int vec_num)
 
const std::string & vector_name (const unsigned int vec_num) const
 
const std::string & vector_name (const NumericVector< Number > &vec_reference) const
 
void set_vector_preservation (const std::string &vec_name, bool preserve)
 
bool vector_preservation (const std::string &vec_name) const
 
NumericVector< Number > & add_adjoint_solution (unsigned int i=0)
 
NumericVector< Number > & get_adjoint_solution (unsigned int i=0)
 
const NumericVector< Number > & get_adjoint_solution (unsigned int i=0) const
 
NumericVector< Number > & add_sensitivity_solution (unsigned int i=0)
 
NumericVector< Number > & get_sensitivity_solution (unsigned int i=0)
 
const NumericVector< Number > & get_sensitivity_solution (unsigned int i=0) const
 
NumericVector< Number > & add_weighted_sensitivity_adjoint_solution (unsigned int i=0)
 
NumericVector< Number > & get_weighted_sensitivity_adjoint_solution (unsigned int i=0)
 
const NumericVector< Number > & get_weighted_sensitivity_adjoint_solution (unsigned int i=0) const
 
NumericVector< Number > & add_weighted_sensitivity_solution ()
 
NumericVector< Number > & get_weighted_sensitivity_solution ()
 
const NumericVector< Number > & get_weighted_sensitivity_solution () const
 
NumericVector< Number > & add_adjoint_rhs (unsigned int i=0)
 
NumericVector< Number > & get_adjoint_rhs (unsigned int i=0)
 
const NumericVector< Number > & get_adjoint_rhs (unsigned int i=0) const
 
NumericVector< Number > & add_sensitivity_rhs (unsigned int i=0)
 
NumericVector< Number > & get_sensitivity_rhs (unsigned int i=0)
 
const NumericVector< Number > & get_sensitivity_rhs (unsigned int i=0) const
 
unsigned int n_vectors () const
 
unsigned int n_vars () const
 
unsigned int n_variable_groups () const
 
unsigned int n_components () const
 
dof_id_type n_dofs () const
 
dof_id_type n_active_dofs () const
 
dof_id_type n_constrained_dofs () const
 
dof_id_type n_local_constrained_dofs () const
 
dof_id_type n_local_dofs () const
 
unsigned int add_variable (const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=NULL)
 
unsigned int add_variable (const std::string &var, const Order order=FIRST, const FEFamily=LAGRANGE, const std::set< subdomain_id_type > *const active_subdomains=NULL)
 
unsigned int add_variables (const std::vector< std::string > &vars, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=NULL)
 
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=NULL)
 
const Variablevariable (unsigned int var) const
 
const VariableGroupvariable_group (unsigned int vg) const
 
bool has_variable (const std::string &var) const
 
const std::string & variable_name (const unsigned int i) const
 
unsigned short int variable_number (const std::string &var) const
 
void get_all_variable_numbers (std::vector< unsigned int > &all_variable_numbers) const
 
unsigned int variable_scalar_number (const std::string &var, unsigned int component) const
 
unsigned int variable_scalar_number (unsigned int var_num, unsigned int component) const
 
const FETypevariable_type (const unsigned int i) const
 
const FETypevariable_type (const std::string &var) const
 
bool identify_variable_groups () const
 
void identify_variable_groups (const bool)
 
Real calculate_norm (const NumericVector< Number > &v, unsigned int var=0, FEMNormType norm_type=L2) const
 
Real calculate_norm (const NumericVector< Number > &v, const SystemNorm &norm) const
 
void read_header (Xdr &io, const std::string &version, const bool read_header=true, const bool read_additional_data=true, const bool read_legacy_format=false)
 
void read_legacy_data (Xdr &io, const bool read_additional_data=true)
 
template<typename ValType >
void read_serialized_data (Xdr &io, const bool read_additional_data=true)
 
void read_serialized_data (Xdr &io, const bool read_additional_data=true)
 
template<typename InValType >
std::size_t read_serialized_vectors (Xdr &io, const std::vector< NumericVector< Number > * > &vectors) const
 
std::size_t read_serialized_vectors (Xdr &io, const std::vector< NumericVector< Number > * > &vectors) const
 
template<typename InValType >
void read_parallel_data (Xdr &io, const bool read_additional_data)
 
void read_parallel_data (Xdr &io, const bool read_additional_data)
 
void write_header (Xdr &io, const std::string &version, const bool write_additional_data) const
 
void write_serialized_data (Xdr &io, const bool write_additional_data=true) const
 
dof_id_type write_serialized_vectors (Xdr &io, const std::vector< const NumericVector< Number > * > &vectors) const
 
void write_parallel_data (Xdr &io, const bool write_additional_data) const
 
std::string get_info () const
 
void attach_init_function (void fptr(EquationSystems &es, const std::string &name))
 
void attach_init_object (Initialization &init)
 
void attach_assemble_function (void fptr(EquationSystems &es, const std::string &name))
 
void attach_assemble_object (Assembly &assemble)
 
void attach_constraint_function (void fptr(EquationSystems &es, const std::string &name))
 
void attach_constraint_object (Constraint &constrain)
 
void attach_QOI_function (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))
 
void attach_QOI_object (QOI &qoi)
 
void attach_QOI_derivative (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))
 
void attach_QOI_derivative_object (QOIDerivative &qoi_derivative)
 
virtual void user_initialization ()
 
virtual void user_assembly ()
 
virtual void user_constrain ()
 
virtual void user_QOI (const QoISet &qoi_indices)
 
virtual void user_QOI_derivative (const QoISet &qoi_indices)
 
virtual void re_update ()
 
virtual void restrict_vectors ()
 
virtual void prolong_vectors ()
 
Number current_solution (const dof_id_type global_dof_number) const
 
Number point_value (unsigned int var, const Point &p, const bool insist_on_success=true) const
 
Number point_value (unsigned int var, const Point &p, const Elem &e) const
 
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
 
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
 
void local_dof_indices (const unsigned int var, std::set< dof_id_type > &var_indices) const
 
void zero_variable (NumericVector< Number > &v, unsigned int var_num) const
 
const Parallel::Communicatorcomm () const
 
processor_id_type n_processors () const
 
processor_id_type processor_id () const
 
virtual void clear_physics ()
 
virtual void init_physics (const System &sys)
 
virtual bool element_time_derivative (bool request_jacobian, DiffContext &)
 
virtual bool element_constraint (bool request_jacobian, DiffContext &)
 
virtual bool side_time_derivative (bool request_jacobian, DiffContext &)
 
virtual bool side_constraint (bool request_jacobian, DiffContext &)
 
virtual void time_evolving (unsigned int var)
 
bool is_time_evolving (unsigned int var) const
 
virtual bool eulerian_residual (bool request_jacobian, DiffContext &)
 
virtual bool eulerian_residual (bool request_jacobian, DiffContext &context)
 
virtual bool mass_residual (bool request_jacobian, DiffContext &)
 
virtual bool mass_residual (bool request_jacobian, DiffContext &)
 
virtual bool side_mass_residual (bool request_jacobian, DiffContext &)
 
virtual void set_mesh_system (System *sys)
 
const Systemget_mesh_system () const
 
Systemget_mesh_system ()
 
virtual void set_mesh_x_var (unsigned int var)
 
unsigned int get_mesh_x_var () const
 
virtual void set_mesh_y_var (unsigned int var)
 
unsigned int get_mesh_y_var () const
 
virtual void set_mesh_z_var (unsigned int var)
 
unsigned int get_mesh_z_var () const
 
virtual void init_qoi (std::vector< Number > &)
 
virtual void clear_qoi ()
 
virtual void element_qoi (DiffContext &, const QoISet &)
 
virtual void element_qoi_derivative (DiffContext &, const QoISet &)
 
virtual void side_qoi (DiffContext &, const QoISet &)
 
virtual void side_qoi_derivative (DiffContext &, const QoISet &)
 
virtual void thread_join (std::vector< Number > &qoi, const std::vector< Number > &other_qoi, const QoISet &qoi_indices)
 
virtual void parallel_op (const Parallel::Communicator &communicator, std::vector< Number > &sys_qoi, std::vector< Number > &local_qoi, const QoISet &qoi_indices)
 

Static Public Member Functions

static std::string get_info ()
 
static void print_info (std::ostream &out=libMesh::out)
 
static unsigned int n_objects ()
 
static void enable_print_counter_info ()
 
static void disable_print_counter_info ()
 

Public Attributes

Realcontinuation_parameter
 
bool quiet
 
Real continuation_parameter_tolerance
 
Real solution_tolerance
 
Real initial_newton_tolerance
 
Real old_continuation_parameter
 
Real min_continuation_parameter
 
Real max_continuation_parameter
 
Real Theta
 
Real Theta_LOCA
 
unsigned int n_backtrack_steps
 
unsigned int n_arclength_reductions
 
Real ds_min
 
Predictor predictor
 
Real newton_stepgrowth_aggressiveness
 
bool newton_progress_check
 
bool fe_reinit_during_postprocess
 
Real numerical_jacobian_h
 
Real verify_analytic_jacobians
 
AutoPtr< TimeSolvertime_solver
 
Real deltat
 
bool postprocess_sides
 
bool print_solution_norms
 
bool print_solutions
 
bool print_residual_norms
 
bool print_residuals
 
bool print_jacobian_norms
 
bool print_jacobians
 
bool print_element_jacobians
 
SparseMatrix< Number > * matrix
 
NumericVector< Number > * rhs
 
bool assemble_before_solve
 
bool use_fixed_solution
 
int extra_quadrature_order
 
AutoPtr< NumericVector< Number > > solution
 
AutoPtr< NumericVector< Number > > current_local_solution
 
Real time
 
std::vector< Numberqoi
 
bool compute_internal_sides
 
bool assemble_qoi_sides
 
bool assemble_qoi_internal_sides
 
bool assemble_qoi_elements
 

Protected Types

enum  RHS_Mode { Residual, G_Lambda }
 
typedef std::map< std::string,
std::pair< unsigned int,
unsigned int > > 
Counts
 

Protected Member Functions

virtual void init_data ()
 
virtual void init_matrices ()
 
void project_vector (NumericVector< Number > &) const
 
void project_vector (const NumericVector< Number > &, NumericVector< Number > &) const
 
void increment_constructor_count (const std::string &name)
 
void increment_destructor_count (const std::string &name)
 

Protected Attributes

RHS_Mode rhs_mode
 
DifferentiablePhysics_diff_physics
 
DifferentiableQoIdiff_qoi
 
const Parallel::Communicator_communicator
 
System_mesh_sys
 
unsigned int _mesh_x_var
 
unsigned int _mesh_y_var
 
unsigned int _mesh_z_var
 
std::vector< bool > _time_evolving
 

Static Protected Attributes

static Counts _counts
 
static Threads::atomic
< unsigned int > 
_n_objects
 
static Threads::spin_mutex _mutex
 
static bool _enable_print_counter = true
 

Private Member Functions

void initialize_tangent ()
 
void solve_tangent ()
 
void update_solution ()
 
void set_Theta ()
 
void set_Theta_LOCA ()
 
void apply_predictor ()
 

Private Attributes

NumericVector< Number > * du_ds
 
NumericVector< Number > * previous_du_ds
 
NumericVector< Number > * previous_u
 
NumericVector< Number > * y
 
NumericVector< Number > * y_old
 
NumericVector< Number > * z
 
NumericVector< Number > * delta_u
 
AutoPtr< LinearSolver< Number > > linear_solver
 
bool tangent_initialized
 
NewtonSolvernewton_solver
 
Real dlambda_ds
 
Real ds
 
Real ds_current
 
Real previous_dlambda_ds
 
Real previous_ds
 
unsigned int newton_step
 

Detailed Description

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

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

Definition at line 57 of file continuation_system.h.

Member Typedef Documentation

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

Definition at line 277 of file implicit_system.h.

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

Definition at line 718 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 113 of file reference_counter.h.

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

Matrix iterator typedefs.

Definition at line 276 of file implicit_system.h.

The type of the parent.

Definition at line 81 of file continuation_system.h.

The type of system.

Definition at line 76 of file continuation_system.h.

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

Syntax sugar to make numerical_jacobian() declaration easier.

Definition at line 196 of file fem_system.h.

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

Vector iterator typedefs.

Definition at line 717 of file system.h.

Member Enumeration Documentation

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

Enumerator
Euler 

First-order Euler predictor

AB2 

Second-order explicit Adams-Bashforth predictor

Invalid_Predictor 

Invalid predictor

Definition at line 223 of file continuation_system.h.

223  {
227  Euler,
228 
232  AB2,
233 
238  };

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

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

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

Enumerator
Residual 
G_Lambda 

Definition at line 288 of file continuation_system.h.

288  {Residual,
289  G_Lambda};

Constructor & Destructor Documentation

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

Constructor. Optionally initializes required data structures.

Definition at line 29 of file continuation_system.C.

32  : Parent(es, name_in, number_in),
34  quiet(true),
36  solution_tolerance(1.e-6),
41  Theta(1.),
42  Theta_LOCA(1.),
43  //tau(1.),
46  ds_min(1.e-8),
52  tangent_initialized(false),
53  newton_solver(NULL),
54  dlambda_ds(0.707),
55  ds(0.1),
56  ds_current(0.1),
58  previous_ds(0.),
59  newton_step(0)
60 {
61  // Warn about using untested code
62  libmesh_experimental();
63 }
libMesh::ContinuationSystem::~ContinuationSystem ( )
virtual

Destructor.

Definition at line 68 of file continuation_system.C.

References clear().

69 {
70  this->clear();
71 }

Member Function Documentation

void libMesh::System::activate ( )
inlineinherited

Activates the system. Only active systems are solved.

Definition at line 1927 of file system.h.

References libMesh::System::_active.

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

Definition at line 1919 of file system.h.

References libMesh::System::_active.

1920 {
1921  return _active;
1922 }
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 1017 of file system.C.

References libMesh::System::add_vector().

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

1018 {
1019  std::ostringstream adjoint_rhs_name;
1020  adjoint_rhs_name << "adjoint_rhs" << i;
1021 
1022  return this->add_vector(adjoint_rhs_name.str(), false);
1023 }
NumericVector< Number > & libMesh::System::add_adjoint_solution ( unsigned int  i = 0)
inherited
Returns
a reference to one of the system's adjoint solution vectors, by default the one corresponding to the first qoi. Creates the vector if it doesn't already exist.

Definition at line 957 of file system.C.

References libMesh::System::add_vector().

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

958 {
959  std::ostringstream adjoint_name;
960  adjoint_name << "adjoint_solution" << i;
961 
962  return this->add_vector(adjoint_name.str());
963 }
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 207 of file implicit_system.C.

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

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

208 {
209  // only add matrices before initializing...
210  if (!_can_add_matrices)
211  {
212  libMesh::err << "ERROR: Too late. Cannot add matrices to the system after initialization"
213  << std::endl
214  << " any more. You should have done this earlier."
215  << std::endl;
216  libmesh_error();
217  }
218 
219  // Return the matrix if it is already there.
220  if (this->have_matrix(mat_name))
221  return *(_matrices[mat_name]);
222 
223  // Otherwise build the matrix and return it.
224  SparseMatrix<Number>* buf = SparseMatrix<Number>::build(this->comm()).release();
225  _matrices.insert (std::make_pair (mat_name, buf));
226 
227  return *buf;
228 }
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 1047 of file system.C.

References libMesh::System::add_vector().

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

1048 {
1049  std::ostringstream sensitivity_rhs_name;
1050  sensitivity_rhs_name << "sensitivity_rhs" << i;
1051 
1052  return this->add_vector(sensitivity_rhs_name.str(), false);
1053 }
NumericVector< Number > & libMesh::System::add_sensitivity_solution ( unsigned int  i = 0)
inherited
Returns
a reference to one of the system's solution sensitivity vectors, by default the one corresponding to the first parameter. Creates the vector if it doesn't already exist.

Definition at line 906 of file system.C.

References libMesh::System::add_vector().

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

907 {
908  std::ostringstream sensitivity_name;
909  sensitivity_name << "sensitivity_solution" << i;
910 
911  return this->add_vector(sensitivity_name.str());
912 }
unsigned int libMesh::System::add_variable ( const std::string &  var,
const FEType type,
const std::set< subdomain_id_type > *const  active_subdomains = NULL 
)
inherited

Adds the variable var to the list of variables for this system. Returns the index number for the new variable.

Definition at line 1077 of file system.C.

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

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

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

References libMesh::System::add_variable().

1160 {
1161  return this->add_variable(var,
1162  FEType(order, family),
1163  active_subdomains);
1164 }
unsigned int libMesh::System::add_variables ( const std::vector< std::string > &  vars,
const FEType type,
const std::set< subdomain_id_type > *const  active_subdomains = NULL 
)
inherited

Adds the variable var to the list of variables for this system. Returns the index number for the new variable.

Definition at line 1168 of file system.C.

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

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

1171 {
1172  // Make sure the variable isn't there already
1173  // or if it is, that it's the type we want
1174  for (unsigned int ov=0; ov<vars.size(); ov++)
1175  for (unsigned int v=0; v<this->n_vars(); v++)
1176  if (this->variable_name(v) == vars[ov])
1177  {
1178  if (this->variable_type(v) == type)
1179  return _variables[v].number();
1180 
1181  libMesh::err << "ERROR: incompatible variable "
1182  << vars[ov]
1183  << " has already been added for this system!"
1184  << std::endl;
1185  libmesh_error();
1186  }
1187 
1188  const unsigned int curr_n_vars = this->n_vars();
1189 
1190  const unsigned int next_first_component = this->n_components();
1191 
1192  // Add the variable group to the list
1193  _variable_groups.push_back((active_subdomains == NULL) ?
1194  VariableGroup(this, vars, curr_n_vars,
1195  next_first_component, type) :
1196  VariableGroup(this, vars, curr_n_vars,
1197  next_first_component, type, *active_subdomains));
1198 
1199  const VariableGroup &vg (_variable_groups.back());
1200 
1201  // Add each component of the group individually
1202  for (unsigned int v=0; v<vars.size(); v++)
1203  {
1204  _variables.push_back (vg(v));
1205  _variable_numbers[vars[v]] = curr_n_vars+v;
1206  }
1207 
1208  libmesh_assert_equal_to ((curr_n_vars+vars.size()), this->n_vars());
1209 
1210  // BSK - Defer this now to System::init_data() so we can detect
1211  // VariableGroups 12/28/2012
1212  // // Add the variable group to the _dof_map
1213  // _dof_map->add_variable_group (vg);
1214 
1215  // Return the number of the new variable
1216  return curr_n_vars+vars.size()-1;
1217 }
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 = NULL 
)
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 1221 of file system.C.

References libMesh::System::add_variables().

1225 {
1226  return this->add_variables(vars,
1227  FEType(order, family),
1228  active_subdomains);
1229 }
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 inititialized 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 676 of file system.C.

References libMesh::System::_can_add_vectors, libMesh::System::_dof_map, libMesh::System::_vector_projections, libMesh::System::_vector_types, libMesh::System::_vectors, libMesh::NumericVector< T >::build(), libMesh::ParallelObject::comm(), libMesh::err, libMeshEnums::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::UnsteadySolver::init(), init_data(), libMesh::NewmarkSystem::NewmarkSystem(), libMesh::System::read_header(), libMesh::FrequencySystem::set_frequencies(), libMesh::FrequencySystem::set_frequencies_by_range(), and libMesh::FrequencySystem::set_frequencies_by_steps().

679 {
680  // Return the vector if it is already there.
681  if (this->have_vector(vec_name))
682  return *(_vectors[vec_name]);
683 
684  // Otherwise build the vector
685  NumericVector<Number>* buf = NumericVector<Number>::build(this->comm()).release();
686  _vectors.insert (std::make_pair (vec_name, buf));
687  _vector_projections.insert (std::make_pair (vec_name, projections));
688 
689  _vector_types.insert (std::make_pair (vec_name, type));
690 
691  // Initialize it if necessary
692  if (!_can_add_vectors)
693  {
694  if(type == GHOSTED)
695  {
696 #ifdef LIBMESH_ENABLE_GHOSTED
697  buf->init (this->n_dofs(), this->n_local_dofs(),
698  _dof_map->get_send_list(), false,
699  GHOSTED);
700 #else
701  libMesh::err << "Cannot initialize ghosted vectors when they are not enabled." << std::endl;
702  libmesh_error();
703 #endif
704  }
705  else
706  buf->init (this->n_dofs(), this->n_local_dofs(), false, type);
707  }
708 
709  return *buf;
710 }
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 987 of file system.C.

References libMesh::System::add_vector().

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

988 {
989  std::ostringstream adjoint_name;
990  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
991 
992  return this->add_vector(adjoint_name.str());
993 }
NumericVector< Number > & libMesh::System::add_weighted_sensitivity_solution ( )
inherited
Returns
a reference to the solution of the last weighted sensitivity solve Creates the vector if it doesn't already exist.

Definition at line 936 of file system.C.

References libMesh::System::add_vector().

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

937 {
938  return this->add_vector("weighted_sensitivity_solution");
939 }
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 701 of file implicit_system.C.

References libMesh::SensitivityData::allocate_data(), libMesh::AutoPtr< Tp >::get(), libMesh::QoISet::has_index(), libMesh::Real, libMesh::ParameterVector::size(), and libMesh::TOLERANCE.

704 {
705  const unsigned int Np = libmesh_cast_int<unsigned int>
706  (parameters.size());
707  const unsigned int Nq = libmesh_cast_int<unsigned int>
708  (qoi.size());
709 
710  // We currently get partial derivatives via central differencing
711  const Real delta_p = TOLERANCE;
712 
713  // An introduction to the problem:
714  //
715  // Residual R(u(p),p) = 0
716  // partial R / partial u = J = system matrix
717  //
718  // This implies that:
719  // d/dp(R) = 0
720  // (partial R / partial p) +
721  // (partial R / partial u) * (partial u / partial p) = 0
722 
723  // We first do an adjoint solve:
724  // J^T * z = (partial q / partial u)
725  // if we havent already or dont have an initial condition for the adjoint
726  if (!this->is_adjoint_already_solved())
727  {
728  this->adjoint_solve(qoi_indices);
729  }
730 
731  // Get ready to fill in senstivities:
732  sensitivities.allocate_data(qoi_indices, *this, parameters);
733 
734  // We use the identities:
735  // dq/dp = (partial q / partial p) + (partial q / partial u) *
736  // (partial u / partial p)
737  // dq/dp = (partial q / partial p) + (J^T * z) *
738  // (partial u / partial p)
739  // dq/dp = (partial q / partial p) + z * J *
740  // (partial u / partial p)
741 
742  // Leading to our final formula:
743  // dq/dp = (partial q / partial p) - z * (partial R / partial p)
744 
745  // In the case of adjoints with heterogenous Dirichlet boundary
746  // function phi, where
747  // q := R(u,phi) + S(u)
748  // the final formula works out to:
749  // dq/dp = (partial S / partial p) - z * (partial R / partial p)
750  // Because we currently have no direct access to
751  // (partial S / partial p), we use the identity
752  // (partial S / partial p) = (partial q / partial p) -
753  // phi * (partial R / partial p)
754  // to derive an equivalent equation:
755  // dq/dp = (partial q / partial p) - (z+phi) * (partial R / partial p)
756 
757  for (unsigned int j=0; j != Np; ++j)
758  {
759  // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
760  // (partial R / partial p) ~= (rhs(p+dp) - rhs(p-dp))/(2*dp)
761 
762  Number old_parameter = *parameters[j];
763  // Number old_qoi = this->qoi;
764 
765  *parameters[j] = old_parameter - delta_p;
766  this->assemble_qoi(qoi_indices);
767  std::vector<Number> qoi_minus = this->qoi;
768 
769  this->assembly(true, false);
770  this->rhs->close();
771 
772  // FIXME - this can and should be optimized to avoid the clone()
773  AutoPtr<NumericVector<Number> > partialR_partialp = this->rhs->clone();
774  *partialR_partialp *= -1;
775 
776  *parameters[j] = old_parameter + delta_p;
777  this->assemble_qoi(qoi_indices);
778  std::vector<Number>& qoi_plus = this->qoi;
779 
780  std::vector<Number> partialq_partialp(Nq, 0);
781  for (unsigned int i=0; i != Nq; ++i)
782  if (qoi_indices.has_index(i))
783  partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
784 
785  this->assembly(true, false);
786  this->rhs->close();
787  *partialR_partialp += *this->rhs;
788  *partialR_partialp /= (2.*delta_p);
789 
790  // Don't leave the parameter changed
791  *parameters[j] = old_parameter;
792 
793  for (unsigned int i=0; i != Nq; ++i)
794  if (qoi_indices.has_index(i))
795  {
796 
798  {
799  AutoPtr<NumericVector<Number> > lift_func =
800  this->get_adjoint_solution(i).zero_clone();
802  (*this, lift_func.get(),
803  /* homogeneous = */ false);
804  sensitivities[i][j] = partialq_partialp[i] -
805  partialR_partialp->dot(*lift_func) -
806  partialR_partialp->dot(this->get_adjoint_solution(i));
807  }
808  else
809  sensitivities[i][j] = partialq_partialp[i] -
810  partialR_partialp->dot(this->get_adjoint_solution(i));
811  }
812  }
813 
814  // All parameters have been reset.
815  // We didn't cache the original rhs or matrix for memory reasons,
816  // but we can restore them to a state consistent solution -
817  // principle of least surprise.
818  this->assembly(true, true);
819  this->rhs->close();
820  this->matrix->close();
821  this->assemble_qoi(qoi_indices);
822 }
std::pair< unsigned int, Real > libMesh::DifferentiableSystem::adjoint_solve ( const QoISet qoi_indices = QoISet())
virtualinherited

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

Reimplemented from libMesh::ImplicitSystem.

Definition at line 144 of file diff_system.C.

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

145 {
146  // Get the time solver object associated with the system, and tell it that
147  // we are solving the adjoint problem
148  this->get_time_solver().set_is_adjoint(true);
149 
150  return this->ImplicitSystem::adjoint_solve(qoi_indices);
151 }
void libMesh::ContinuationSystem::advance_arcstep ( )

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

Definition at line 947 of file continuation_system.C.

References solve_tangent(), and update_solution().

948 {
949  // Solve for the updated tangent du1/ds, d(lambda1)/ds
950  solve_tangent();
951 
952  // Advance the solution and the parameter to the next value.
953  update_solution();
954 }
void libMesh::ContinuationSystem::apply_predictor ( )
private

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

Definition at line 1386 of file continuation_system.C.

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

Referenced by continuation_solve(), and update_solution().

1387 {
1388  if (predictor == Euler)
1389  {
1390  // 1.) Euler Predictor
1391  // Predict next the solution
1392  solution->add(ds_current, *du_ds);
1393  solution->close();
1394 
1395  // Predict next parameter value
1397  }
1398 
1399 
1400  else if (predictor == AB2)
1401  {
1402  // 2.) 2nd-order explicit AB predictor
1403  libmesh_assert_not_equal_to (previous_ds, 0.0);
1404  const Real stepratio = ds_current/previous_ds;
1405 
1406  // Build up next solution value.
1407  solution->add( 0.5*ds_current*(2.+stepratio), *du_ds);
1408  solution->add(-0.5*ds_current*stepratio , *previous_du_ds);
1409  solution->close();
1410 
1411  // Next parameter value
1413  0.5*ds_current*((2.+stepratio)*dlambda_ds -
1414  stepratio*previous_dlambda_ds);
1415  }
1416 
1417  else
1418  {
1419  // Unknown predictor
1420  libmesh_error();
1421  }
1422 
1423 }
void libMesh::DifferentiableSystem::assemble ( )
virtualinherited

Prepares matrix and rhs for matrix assembly. Users should not reimplement this

Reimplemented from libMesh::ImplicitSystem.

Definition at line 125 of file diff_system.C.

References libMesh::DifferentiableSystem::assembly().

126 {
127  this->assembly(true, true);
128 }
void libMesh::FEMSystem::assemble_qoi ( const QoISet indices = QoISet())
virtualinherited

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

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

Reimplemented from libMesh::ExplicitSystem.

Definition at line 773 of file fem_system.C.

References libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::ParallelObject::comm(), libMesh::DifferentiableSystem::diff_qoi, libMesh::System::get_mesh(), libMesh::QoISet::has_index(), mesh, libMesh::DifferentiableQoI::parallel_op(), libMesh::Threads::parallel_reduce(), libMesh::System::qoi, libMesh::START_LOG(), libMesh::STOP_LOG(), and libMesh::System::update().

774 {
775  START_LOG("assemble_qoi()", "FEMSystem");
776 
777  const MeshBase& mesh = this->get_mesh();
778 
779  this->update();
780 
781  const unsigned int Nq = libmesh_cast_int<unsigned int>(qoi.size());
782 
783  // the quantity of interest is assumed to be a sum of element and
784  // side terms
785  for (unsigned int i=0; i != Nq; ++i)
786  if (qoi_indices.has_index(i))
787  qoi[i] = 0;
788 
789  // Create a non-temporary qoi_contributions object, so we can query
790  // its results after the reduction
791  QoIContributions qoi_contributions(*this, *(this->diff_qoi), qoi_indices);
792 
793  // Loop over every active mesh element on this processor
796  qoi_contributions);
797 
798  this->diff_qoi->parallel_op( this->comm(), this->qoi, qoi_contributions.qoi, qoi_indices );
799 
800  STOP_LOG("assemble_qoi()", "FEMSystem");
801 }
void libMesh::FEMSystem::assemble_qoi_derivative ( const QoISet indices = QoISet())
virtualinherited

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

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

Reimplemented from libMesh::ExplicitSystem.

Definition at line 805 of file fem_system.C.

References libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::System::add_adjoint_rhs(), libMesh::DifferentiableSystem::diff_qoi, libMesh::System::get_mesh(), libMesh::QoISet::has_index(), mesh, libMesh::Threads::parallel_for(), libMesh::System::qoi, libMesh::START_LOG(), libMesh::STOP_LOG(), libMesh::System::update(), and libMesh::NumericVector< T >::zero().

806 {
807  START_LOG("assemble_qoi_derivative()", "FEMSystem");
808 
809  const MeshBase& mesh = this->get_mesh();
810 
811  this->update();
812 
813  // The quantity of interest derivative assembly accumulates on
814  // initially zero vectors
815  for (unsigned int i=0; i != qoi.size(); ++i)
816  if (qoi_indices.has_index(i))
817  this->add_adjoint_rhs(i).zero();
818 
819  // Loop over every active mesh element on this processor
820  Threads::parallel_for(elem_range.reset(mesh.active_local_elements_begin(),
822  QoIDerivativeContributions(*this, qoi_indices,
823  *(this->diff_qoi)));
824 
825  STOP_LOG("assemble_qoi_derivative()", "FEMSystem");
826 }
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 overloaded in derived classes.

Reimplemented from libMesh::System.

Definition at line 665 of file implicit_system.C.

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

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

666 {
667  const unsigned int Np = libmesh_cast_int<unsigned int>
668  (parameters.size());
669  Real deltap = TOLERANCE;
670 
671  for (unsigned int p=0; p != Np; ++p)
672  {
673  NumericVector<Number> &sensitivity_rhs = this->add_sensitivity_rhs(p);
674 
675  // Approximate -(partial R / partial p) by
676  // (R(p-dp) - R(p+dp)) / (2*dp)
677 
678  Number old_parameter = *parameters[p];
679  *parameters[p] -= deltap;
680 
681  this->assembly(true, false);
682  this->rhs->close();
683  sensitivity_rhs = *this->rhs;
684 
685  *parameters[p] = old_parameter + deltap;
686 
687  this->assembly(true, false);
688  this->rhs->close();
689 
690  sensitivity_rhs -= *this->rhs;
691  sensitivity_rhs /= (2*deltap);
692  sensitivity_rhs.close();
693 
694  *parameters[p] = old_parameter;
695  }
696 }
void libMesh::FEMSystem::assembly ( bool  get_residual,
bool  get_jacobian 
)
virtualinherited

Prepares matrix or rhs for matrix assembly. Users may reimplement this to add pre- or post-assembly code before or after calling FEMSystem::assembly()

Implements libMesh::DifferentiableSystem.

Definition at line 585 of file fem_system.C.

References libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::NumericVector< T >::close(), libMesh::SparseMatrix< T >::close(), libMesh::err, libMesh::System::get_mesh(), libMesh::NumericVector< T >::l1_norm(), libMesh::SparseMatrix< T >::l1_norm(), libMesh::libmesh_assert(), libMesh::ImplicitSystem::matrix, mesh, libMesh::out, libMesh::Threads::parallel_for(), libMesh::BasicOStreamProxy< charT, traits >::precision(), libMesh::DifferentiableSystem::print_jacobian_norms, libMesh::DifferentiableSystem::print_jacobians, libMesh::DifferentiableSystem::print_residual_norms, libMesh::DifferentiableSystem::print_residuals, libMesh::DifferentiableSystem::print_solution_norms, libMesh::DifferentiableSystem::print_solutions, libMesh::ExplicitSystem::rhs, libMesh::System::solution, libMesh::START_LOG(), libMesh::STOP_LOG(), libMesh::DifferentiableSystem::time_solver, libMesh::System::update(), libMesh::FEMSystem::verify_analytic_jacobians, libMesh::NumericVector< T >::zero(), and libMesh::SparseMatrix< T >::zero().

Referenced by continuation_solve(), and solve_tangent().

586 {
587  libmesh_assert(get_residual || get_jacobian);
588  std::string log_name;
589  if (get_residual && get_jacobian)
590  log_name = "assembly()";
591  else if (get_residual)
592  log_name = "assembly(get_residual)";
593  else
594  log_name = "assembly(get_jacobian)";
595 
596  START_LOG(log_name, "FEMSystem");
597 
598  const MeshBase& mesh = this->get_mesh();
599 
600 // this->get_vector("_nonlinear_solution").localize
601 // (*current_local_nonlinear_solution,
602 // dof_map.get_send_list());
603  this->update();
604 
606  {
607 // this->get_vector("_nonlinear_solution").close();
608  this->solution->close();
609 
610  std::streamsize old_precision = libMesh::out.precision();
612  libMesh::out << "|U| = "
613 // << this->get_vector("_nonlinear_solution").l1_norm()
614  << this->solution->l1_norm()
615  << std::endl;
616  libMesh::out.precision(old_precision);
617  }
618  if (print_solutions)
619  {
620  std::streamsize old_precision = libMesh::out.precision();
622 // libMesh::out << "U = [" << this->get_vector("_nonlinear_solution")
623  libMesh::out << "U = [" << *(this->solution)
624  << "];" << std::endl;
625  libMesh::out.precision(old_precision);
626  }
627 
628  // Is this definitely necessary? [RHS]
629  // Yes. [RHS 2012]
630  if (get_jacobian)
631  matrix->zero();
632  if (get_residual)
633  rhs->zero();
634 
635  // Stupid C++ lets you set *Real* verify_analytic_jacobians = true!
636  if (verify_analytic_jacobians > 0.5)
637  {
638  libMesh::err << "WARNING! verify_analytic_jacobians was set "
639  << "to absurdly large value of "
640  << verify_analytic_jacobians << std::endl;
641  libMesh::err << "Resetting to 1e-6!" << std::endl;
643  }
644 
645  // In time-dependent problems, the nonlinear function we're trying
646  // to solve at each timestep may depend on the particular solver
647  // we're using
649 
650  // Build the residual and jacobian contributions on every active
651  // mesh element on this processor
652  Threads::parallel_for(elem_range.reset(mesh.active_local_elements_begin(),
654  AssemblyContributions(*this, get_residual, get_jacobian));
655 
656 
657  if (get_residual && (print_residual_norms || print_residuals))
658  this->rhs->close();
659  if (get_residual && print_residual_norms)
660  {
661  std::streamsize old_precision = libMesh::out.precision();
663  libMesh::out << "|F| = " << this->rhs->l1_norm() << std::endl;
664  libMesh::out.precision(old_precision);
665  }
666  if (get_residual && print_residuals)
667  {
668  std::streamsize old_precision = libMesh::out.precision();
670  libMesh::out << "F = [" << *(this->rhs) << "];" << std::endl;
671  libMesh::out.precision(old_precision);
672  }
673 
674  if (get_jacobian && (print_jacobian_norms || print_jacobians))
675  this->matrix->close();
676  if (get_jacobian && print_jacobian_norms)
677  {
678  std::streamsize old_precision = libMesh::out.precision();
680  libMesh::out << "|J| = " << this->matrix->l1_norm() << std::endl;
681  libMesh::out.precision(old_precision);
682  }
683  if (get_jacobian && print_jacobians)
684  {
685  std::streamsize old_precision = libMesh::out.precision();
687  libMesh::out << "J = [" << *(this->matrix) << "];" << std::endl;
688  libMesh::out.precision(old_precision);
689  }
690  STOP_LOG(log_name, "FEMSystem");
691 }
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 1754 of file system.C.

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

1756 {
1757  libmesh_assert(fptr);
1758 
1759  if (_assemble_system_object != NULL)
1760  {
1761  libmesh_here();
1762  libMesh::out << "WARNING: Cannot specify both assembly function and object!"
1763  << std::endl;
1764 
1765  _assemble_system_object = NULL;
1766  }
1767 
1769 }
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 1773 of file system.C.

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

1774 {
1775  if (_assemble_system_function != NULL)
1776  {
1777  libmesh_here();
1778  libMesh::out << "WARNING: Cannot specify both assembly object and function!"
1779  << std::endl;
1780 
1782  }
1783 
1784  _assemble_system_object = &assemble_in;
1785 }
void libMesh::System::attach_constraint_function ( void   fptrEquationSystems &es,const std::string &name)
inherited

Register a user function for imposing constraints.

Definition at line 1789 of file system.C.

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

1791 {
1792  libmesh_assert(fptr);
1793 
1794  if (_constrain_system_object != NULL)
1795  {
1796  libmesh_here();
1797  libMesh::out << "WARNING: Cannot specify both constraint function and object!"
1798  << std::endl;
1799 
1800  _constrain_system_object = NULL;
1801  }
1802 
1804 }
void libMesh::System::attach_constraint_object ( System::Constraint constrain)
inherited

Register a user object for imposing constraints.

Definition at line 1808 of file system.C.

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

1809 {
1810  if (_constrain_system_function != NULL)
1811  {
1812  libmesh_here();
1813  libMesh::out << "WARNING: Cannot specify both constraint object and function!"
1814  << std::endl;
1815 
1817  }
1818 
1819  _constrain_system_object = &constrain;
1820 }
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 1719 of file system.C.

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

1721 {
1722  libmesh_assert(fptr);
1723 
1724  if (_init_system_object != NULL)
1725  {
1726  libmesh_here();
1727  libMesh::out << "WARNING: Cannot specify both initialization function and object!"
1728  << std::endl;
1729 
1730  _init_system_object = NULL;
1731  }
1732 
1733  _init_system_function = fptr;
1734 }
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 1738 of file system.C.

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

1739 {
1740  if (_init_system_function != NULL)
1741  {
1742  libmesh_here();
1743  libMesh::out << "WARNING: Cannot specify both initialization object and function!"
1744  << std::endl;
1745 
1746  _init_system_function = NULL;
1747  }
1748 
1749  _init_system_object = &init_in;
1750 }
void libMesh::DifferentiableSystem::attach_physics ( DifferentiablePhysics physics_in)
inlineinherited

Attach external Physics object.

Definition at line 176 of file diff_system.h.

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

177  { this->_diff_physics = (physics_in->clone_physics()).release();
178  this->_diff_physics->init_physics(*this);}
void libMesh::DifferentiableSystem::attach_qoi ( DifferentiableQoI qoi_in)
inlineinherited

Attach external QoI object.

Definition at line 197 of file diff_system.h.

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

198  { this->diff_qoi = (qoi_in->clone()).release();
199  // User needs to resize qoi system qoi accordingly
200  this->diff_qoi->init_qoi( this->qoi );}
void libMesh::System::attach_QOI_derivative ( void   fptrEquationSystems &es,const std::string &name,const QoISet &qoi_indices)
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 1860 of file system.C.

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

1863 {
1864  libmesh_assert(fptr);
1865 
1866  if (_qoi_evaluate_derivative_object != NULL)
1867  {
1868  libmesh_here();
1869  libMesh::out << "WARNING: Cannot specify both QOI derivative function and object!"
1870  << std::endl;
1871 
1873  }
1874 
1876 }
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 1880 of file system.C.

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

1881 {
1883  {
1884  libmesh_here();
1885  libMesh::out << "WARNING: Cannot specify both QOI derivative object and function!"
1886  << std::endl;
1887 
1889  }
1890 
1891  _qoi_evaluate_derivative_object = &qoi_derivative;
1892 }
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 1824 of file system.C.

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

1827 {
1828  libmesh_assert(fptr);
1829 
1830  if (_qoi_evaluate_object != NULL)
1831  {
1832  libmesh_here();
1833  libMesh::out << "WARNING: Cannot specify both QOI function and object!"
1834  << std::endl;
1835 
1836  _qoi_evaluate_object = NULL;
1837  }
1838 
1839  _qoi_evaluate_function = fptr;
1840 }
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 1844 of file system.C.

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

1845 {
1846  if (_qoi_evaluate_function != NULL)
1847  {
1848  libmesh_here();
1849  libMesh::out << "WARNING: Cannot specify both QOI object and function!"
1850  << std::endl;
1851 
1852  _qoi_evaluate_function = NULL;
1853  }
1854 
1855  _qoi_evaluate_object = &qoi_in;
1856 }
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 = NULL 
)
inherited

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

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

Definition at line 665 of file system_projection.C.

669 {
670  this->boundary_project_vector(b, variables, *solution, f, g);
671 
672  solution->localize(*current_local_solution);
673 }
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.

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

Definition at line 641 of file system_projection.C.

652 {
653  WrappedFunction<Number> f(*this, fptr, &parameters);
654  WrappedFunction<Gradient> g(*this, gptr, &parameters);
655  this->boundary_project_solution(b, variables, &f, &g);
656 }
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 = NULL 
) const
inherited

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

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

Definition at line 707 of file system_projection.C.

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

712 {
713  START_LOG ("boundary_project_vector()", "System");
714 
716  (ConstElemRange (this->get_mesh().active_local_elements_begin(),
717  this->get_mesh().active_local_elements_end() ),
718  BoundaryProjectSolution(b, variables, *this, f, g,
719  this->get_equation_systems().parameters,
720  new_vector)
721  );
722 
723  // We don't do SCALAR dofs when just projecting the boundary, so
724  // we're done here.
725 
726  new_vector.close();
727 
728 #ifdef LIBMESH_ENABLE_CONSTRAINTS
729  this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
730 #endif
731 
732  STOP_LOG("boundary_project_vector()", "System");
733 }
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 
) 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.

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

Definition at line 684 of file system_projection.C.

696 {
697  WrappedFunction<Number> f(*this, fptr, &parameters);
698  WrappedFunction<Gradient> g(*this, gptr, &parameters);
699  this->boundary_project_vector(b, variables, new_vector, &f, &g);
700 }
AutoPtr< DiffContext > libMesh::FEMSystem::build_context ( )
virtualinherited

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

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

Reimplemented from libMesh::DifferentiableSystem.

Definition at line 951 of file fem_system.C.

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

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

952 {
953  FEMContext* fc = new FEMContext(*this);
954 
955  AutoPtr<DiffContext> ap(fc);
956 
957  DifferentiablePhysics* phys = this->get_physics();
958 
959  libmesh_assert (phys);
960 
961  // If we are solving a moving mesh problem, tell that to the Context
962  fc->set_mesh_system(phys->get_mesh_system());
963  fc->set_mesh_x_var(phys->get_mesh_x_var());
964  fc->set_mesh_y_var(phys->get_mesh_y_var());
965  fc->set_mesh_z_var(phys->get_mesh_z_var());
966 
967  ap->set_deltat_pointer( &deltat );
968 
969  // If we are solving the adjoint problem, tell that to the Context
970  ap->is_adjoint() = this->get_time_solver().is_adjoint();
971 
972  return ap;
973 }
Real libMesh::System::calculate_norm ( const NumericVector< Number > &  v,
unsigned int  var = 0,
FEMNormType  norm_type = L2 
) 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 1380 of file system.C.

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

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

1383 {
1384  //short circuit to save time
1385  if(norm_type == DISCRETE_L1 ||
1386  norm_type == DISCRETE_L2 ||
1387  norm_type == DISCRETE_L_INF)
1388  return discrete_var_norm(v,var,norm_type);
1389 
1390  // Not a discrete norm
1391  std::vector<FEMNormType> norms(this->n_vars(), L2);
1392  std::vector<Real> weights(this->n_vars(), 0.0);
1393  norms[var] = norm_type;
1394  weights[var] = 1.0;
1395  Real val = this->calculate_norm(v, SystemNorm(norms, weights));
1396  return val;
1397 }
Real libMesh::System::calculate_norm ( const NumericVector< Number > &  v,
const SystemNorm norm 
) 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 1401 of file system.C.

References libMesh::System::_dof_map, std::abs(), libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::TypeVector< T >::add_scaled(), libMesh::TypeTensor< T >::add_scaled(), libMesh::NumericVector< T >::build(), libMesh::FEGenericBase< T >::build(), libMesh::ParallelObject::comm(), libMesh::FEType::default_quadrature_rule(), libMesh::dim, libMeshEnums::DISCRETE_L1, libMeshEnums::DISCRETE_L2, libMeshEnums::DISCRETE_L_INF, libMesh::System::discrete_var_norm(), libMesh::DofMap::dof_indices(), libMesh::AutoPtr< Tp >::get(), libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMeshEnums::H1, libMeshEnums::H1_SEMINORM, libMeshEnums::H2, libMeshEnums::H2_SEMINORM, libMesh::SystemNorm::is_discrete(), libMeshEnums::L1, libMesh::NumericVector< T >::l1_norm(), libMeshEnums::L2, libMesh::NumericVector< T >::l2_norm(), libMeshEnums::L_INF, libMesh::NumericVector< T >::linfty_norm(), libMesh::NumericVector< T >::localize(), std::max(), libMesh::Parallel::Communicator::max(), libMesh::MeshBase::mesh_dimension(), libMesh::System::n_vars(), libMesh::TensorTools::norm_sq(), libMesh::Real, libMeshEnums::SERIAL, libMesh::TypeVector< T >::size(), libMesh::TypeTensor< T >::size(), libMesh::NumericVector< T >::size(), libMesh::TypeVector< T >::size_sq(), libMesh::TypeTensor< T >::size_sq(), libMesh::START_LOG(), libMesh::STOP_LOG(), libMesh::Parallel::Communicator::sum(), libMesh::SystemNorm::type(), libMesh::DofMap::variable_type(), libMeshEnums::W1_INF_SEMINORM, libMeshEnums::W2_INF_SEMINORM, libMesh::SystemNorm::weight(), and libMesh::SystemNorm::weight_sq().

1403 {
1404  // This function must be run on all processors at once
1405  parallel_object_only();
1406 
1407  START_LOG ("calculate_norm()", "System");
1408 
1409  // Zero the norm before summation
1410  Real v_norm = 0.;
1411 
1412  if (norm.is_discrete())
1413  {
1414  STOP_LOG ("calculate_norm()", "System");
1415  //Check to see if all weights are 1.0 and all types are equal
1416  FEMNormType norm_type0 = norm.type(0);
1417  unsigned int check_var = 0;
1418  for (; check_var != this->n_vars(); ++check_var)
1419  if((norm.weight(check_var) != 1.0) || (norm.type(check_var) != norm_type0))
1420  break;
1421 
1422  //All weights were 1.0 so just do the full vector discrete norm
1423  if(check_var == this->n_vars())
1424  {
1425  if(norm_type0 == DISCRETE_L1)
1426  return v.l1_norm();
1427  if(norm_type0 == DISCRETE_L2)
1428  return v.l2_norm();
1429  if(norm_type0 == DISCRETE_L_INF)
1430  return v.linfty_norm();
1431  else
1432  libmesh_error();
1433  }
1434 
1435  for (unsigned int var=0; var != this->n_vars(); ++var)
1436  {
1437  // Skip any variables we don't need to integrate
1438  if (norm.weight(var) == 0.0)
1439  continue;
1440 
1441  v_norm += norm.weight(var) * discrete_var_norm(v, var, norm.type(var));
1442  }
1443 
1444  return v_norm;
1445  }
1446 
1447  // Localize the potentially parallel vector
1448  AutoPtr<NumericVector<Number> > local_v = NumericVector<Number>::build(this->comm());
1449  local_v->init(v.size(), true, SERIAL);
1450  v.localize (*local_v, _dof_map->get_send_list());
1451 
1452  unsigned int dim = this->get_mesh().mesh_dimension();
1453 
1454  // I'm not sure how best to mix Hilbert norms on some variables (for
1455  // which we'll want to square then sum then square root) with norms
1456  // like L_inf (for which we'll just want to take an absolute value
1457  // and then sum).
1458  bool using_hilbert_norm = true,
1459  using_nonhilbert_norm = true;
1460 
1461  // Loop over all variables
1462  for (unsigned int var=0; var != this->n_vars(); ++var)
1463  {
1464  // Skip any variables we don't need to integrate
1465  Real norm_weight_sq = norm.weight_sq(var);
1466  if (norm_weight_sq == 0.0)
1467  continue;
1468  Real norm_weight = norm.weight(var);
1469 
1470  // Check for unimplemented norms (rather than just returning 0).
1471  FEMNormType norm_type = norm.type(var);
1472  if((norm_type==H1) ||
1473  (norm_type==H2) ||
1474  (norm_type==L2) ||
1475  (norm_type==H1_SEMINORM) ||
1476  (norm_type==H2_SEMINORM))
1477  {
1478  if (!using_hilbert_norm)
1479  libmesh_not_implemented();
1480  using_nonhilbert_norm = false;
1481  }
1482  else if ((norm_type==L1) ||
1483  (norm_type==L_INF) ||
1484  (norm_type==W1_INF_SEMINORM) ||
1485  (norm_type==W2_INF_SEMINORM))
1486  {
1487  if (!using_nonhilbert_norm)
1488  libmesh_not_implemented();
1489  using_hilbert_norm = false;
1490  }
1491  else
1492  libmesh_not_implemented();
1493 
1494  const FEType& fe_type = this->get_dof_map().variable_type(var);
1495  AutoPtr<QBase> qrule =
1496  fe_type.default_quadrature_rule (dim);
1497  AutoPtr<FEBase> fe
1498  (FEBase::build(dim, fe_type));
1499  fe->attach_quadrature_rule (qrule.get());
1500 
1501  const std::vector<Real>& JxW = fe->get_JxW();
1502  const std::vector<std::vector<Real> >* phi = NULL;
1503  if (norm_type == H1 ||
1504  norm_type == H2 ||
1505  norm_type == L2 ||
1506  norm_type == L1 ||
1507  norm_type == L_INF)
1508  phi = &(fe->get_phi());
1509 
1510  const std::vector<std::vector<RealGradient> >* dphi = NULL;
1511  if (norm_type == H1 ||
1512  norm_type == H2 ||
1513  norm_type == H1_SEMINORM ||
1514  norm_type == W1_INF_SEMINORM)
1515  dphi = &(fe->get_dphi());
1516 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1517  const std::vector<std::vector<RealTensor> >* d2phi = NULL;
1518  if (norm_type == H2 ||
1519  norm_type == H2_SEMINORM ||
1520  norm_type == W2_INF_SEMINORM)
1521  d2phi = &(fe->get_d2phi());
1522 #endif
1523 
1524  std::vector<dof_id_type> dof_indices;
1525 
1526  // Begin the loop over the elements
1527  MeshBase::const_element_iterator el =
1529  const MeshBase::const_element_iterator end_el =
1531 
1532  for ( ; el != end_el; ++el)
1533  {
1534  const Elem* elem = *el;
1535 
1536  fe->reinit (elem);
1537 
1538  this->get_dof_map().dof_indices (elem, dof_indices, var);
1539 
1540  const unsigned int n_qp = qrule->n_points();
1541 
1542  const unsigned int n_sf = libmesh_cast_int<unsigned int>
1543  (dof_indices.size());
1544 
1545  // Begin the loop over the Quadrature points.
1546  for (unsigned int qp=0; qp<n_qp; qp++)
1547  {
1548  if (norm_type == L1)
1549  {
1550  Number u_h = 0.;
1551  for (unsigned int i=0; i != n_sf; ++i)
1552  u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1553  v_norm += norm_weight *
1554  JxW[qp] * std::abs(u_h);
1555  }
1556 
1557  if (norm_type == L_INF)
1558  {
1559  Number u_h = 0.;
1560  for (unsigned int i=0; i != n_sf; ++i)
1561  u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1562  v_norm = std::max(v_norm, norm_weight * std::abs(u_h));
1563  }
1564 
1565  if (norm_type == H1 ||
1566  norm_type == H2 ||
1567  norm_type == L2)
1568  {
1569  Number u_h = 0.;
1570  for (unsigned int i=0; i != n_sf; ++i)
1571  u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
1572  v_norm += norm_weight_sq *
1573  JxW[qp] * TensorTools::norm_sq(u_h);
1574  }
1575 
1576  if (norm_type == H1 ||
1577  norm_type == H2 ||
1578  norm_type == H1_SEMINORM)
1579  {
1580  Gradient grad_u_h;
1581  for (unsigned int i=0; i != n_sf; ++i)
1582  grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
1583  v_norm += norm_weight_sq *
1584  JxW[qp] * grad_u_h.size_sq();
1585  }
1586 
1587  if (norm_type == W1_INF_SEMINORM)
1588  {
1589  Gradient grad_u_h;
1590  for (unsigned int i=0; i != n_sf; ++i)
1591  grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
1592  v_norm = std::max(v_norm, norm_weight * grad_u_h.size());
1593  }
1594 
1595 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1596  if (norm_type == H2 ||
1597  norm_type == H2_SEMINORM)
1598  {
1599  Tensor hess_u_h;
1600  for (unsigned int i=0; i != n_sf; ++i)
1601  hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
1602  v_norm += norm_weight_sq *
1603  JxW[qp] * hess_u_h.size_sq();
1604  }
1605 
1606  if (norm_type == W2_INF_SEMINORM)
1607  {
1608  Tensor hess_u_h;
1609  for (unsigned int i=0; i != n_sf; ++i)
1610  hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
1611  v_norm = std::max(v_norm, norm_weight * hess_u_h.size());
1612  }
1613 #endif
1614  }
1615  }
1616  }
1617 
1618  if (using_hilbert_norm)
1619  {
1620  this->comm().sum(v_norm);
1621  v_norm = std::sqrt(v_norm);
1622  }
1623  else
1624  {
1625  this->comm().max(v_norm);
1626  }
1627 
1628  STOP_LOG ("calculate_norm()", "System");
1629 
1630  return v_norm;
1631 }
void libMesh::ContinuationSystem::clear ( )
virtual

Clear all the data structures associated with the system.

Reimplemented from libMesh::FEMSystem.

Definition at line 76 of file continuation_system.C.

References libMesh::FEMSystem::clear().

Referenced by ~ContinuationSystem().

77 {
78  // FIXME: Do anything here, e.g. zero vectors, etc?
79 
80  // Call the Parent's clear function
81  Parent::clear();
82 }
virtual void libMesh::DifferentiablePhysics::clear_physics ( )
virtualinherited

Clear any data structures associated with the physics.

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

virtual void libMesh::DifferentiableQoI::clear_qoi ( )
inlinevirtualinherited

Clear all the data structures associated with the QoI.

Definition at line 75 of file diff_qoi.h.

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

75 {}
virtual AutoPtr<DifferentiableQoI> libMesh::DifferentiableSystem::clone ( )
inlinevirtualinherited

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

Implements libMesh::DifferentiableQoI.

Definition at line 153 of file diff_system.h.

154  { libmesh_error();
155  // dummy
156  return AutoPtr<DifferentiableQoI>(this); }
virtual AutoPtr<DifferentiablePhysics> libMesh::DifferentiableSystem::clone_physics ( )
inlinevirtualinherited

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

Implements libMesh::DifferentiablePhysics.

Definition at line 145 of file diff_system.h.

146  { libmesh_error();
147  // dummy
148  return AutoPtr<DifferentiablePhysics>(this); }
const Parallel::Communicator& libMesh::ParallelObject::comm ( ) const
inlineinherited
Returns
a reference to the Parallel::Communicator object used by this mesh.

Definition at line 86 of file parallel_object.h.

References libMesh::ParallelObject::_communicator.

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), libMesh::__libmesh_petsc_diff_solver_monitor(), libMesh::__libmesh_petsc_diff_solver_residual(), libMesh::__libmesh_petsc_snes_jacobian(), libMesh::__libmesh_petsc_snes_residual(), libMesh::MeshRefinement::_coarsen_elements(), libMesh::ExactSolution::_compute_error(), libMesh::MetisPartitioner::_do_partition(), libMesh::ParmetisPartitioner::_do_repartition(), libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_mult(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult_add(), libMesh::EquationSystems::_read_impl(), libMesh::MeshRefinement::_refine_elements(), libMesh::ParallelMesh::add_elem(), libMesh::ImplicitSystem::add_matrix(), libMesh::ParallelMesh::add_node(), libMesh::System::add_vector(), libMesh::UnstructuredMesh::all_second_order(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::FEMSystem::assemble_qoi(), libMesh::MeshCommunication::assign_global_indices(), libMesh::ParmetisPartitioner::assign_partitioning(), libMesh::DofMap::attach_matrix(), libMesh::MeshTools::bounding_box(), libMesh::System::calculate_norm(), libMesh::MeshRefinement::coarsen_elements(), libMesh::Nemesis_IO_Helper::compute_num_global_elem_blocks(), libMesh::Nemesis_IO_Helper::compute_num_global_nodesets(), libMesh::Nemesis_IO_Helper::compute_num_global_sidesets(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), libMesh::MeshTools::correct_node_proc_ids(), libMesh::MeshCommunication::delete_remote_elements(), libMesh::DofMap::distribute_dofs(), DMlibMeshFunction(), DMlibMeshJacobian(), DMLibMeshSetSystem(), DMVariableBounds_libMesh(), libMesh::MeshRefinement::eliminate_unrefined_patches(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::MeshRefinement::flag_elements_by_elem_fraction(), libMesh::MeshRefinement::flag_elements_by_error_fraction(), libMesh::MeshRefinement::flag_elements_by_nelem_target(), libMesh::for(), libMesh::CondensedEigenSystem::get_eigenpair(), libMesh::ImplicitSystem::get_linear_solver(), libMesh::LocationMap< T >::init(), libMesh::PetscDiffSolver::init(), libMesh::TimeSolver::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::EigenSystem::init_data(), libMesh::EigenSystem::init_matrices(), libMesh::ParmetisPartitioner::initialize(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), libMesh::ParallelMesh::libmesh_assert_valid_parallel_flags(), libMesh::MeshTools::libmesh_assert_valid_procids< Elem >(), libMesh::MeshTools::libmesh_assert_valid_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_refinement_flags(), libMesh::MeshRefinement::limit_level_mismatch_at_edge(), libMesh::MeshRefinement::limit_level_mismatch_at_node(), libMesh::MeshRefinement::make_coarsening_compatible(), libMesh::MeshCommunication::make_elems_parallel_consistent(), libMesh::MeshRefinement::make_flags_parallel_consistent(), libMesh::MeshCommunication::make_node_ids_parallel_consistent(), libMesh::MeshCommunication::make_node_proc_ids_parallel_consistent(), libMesh::MeshCommunication::make_nodes_parallel_consistent(), libMesh::MeshRefinement::make_refinement_compatible(), libMesh::FEMSystem::mesh_position_set(), libMesh::MeshSerializer::MeshSerializer(), libMesh::ParallelMesh::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::ParallelMesh::parallel_max_elem_id(), libMesh::ParallelMesh::parallel_max_node_id(), libMesh::ParallelMesh::parallel_n_elem(), libMesh::ParallelMesh::parallel_n_nodes(), libMesh::Partitioner::partition(), libMesh::Partitioner::partition_unpartitioned_elements(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::MeshBase::prepare_for_use(), libMesh::System::project_vector(), libMesh::Nemesis_IO::read(), libMesh::XdrIO::read(), libMesh::System::read_header(), libMesh::System::read_legacy_data(), libMesh::System::read_SCALAR_dofs(), libMesh::XdrIO::read_serialized_bc_names(), libMesh::XdrIO::read_serialized_bcs(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::System::read_serialized_vector(), libMesh::MeshBase::recalculate_n_partitions(), libMesh::MeshRefinement::refine_and_coarsen_elements(), libMesh::MeshRefinement::refine_elements(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::MeshBase::subdomain_ids(), libMesh::BoundaryInfo::sync(), libMesh::Parallel::sync_element_data_by_parent_id(), libMesh::MeshRefinement::test_level_one(), libMesh::MeshRefinement::test_unflagged(), libMesh::MeshTools::total_weight(), libMesh::CheckpointIO::write(), libMesh::XdrIO::write(), libMesh::UnstructuredMesh::write(), libMesh::LegacyXdrIO::write_mesh(), libMesh::System::write_SCALAR_dofs(), libMesh::XdrIO::write_serialized_bcs(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), and libMesh::DivaIO::write_stream().

87  { return _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 526 of file system.C.

References libMesh::System::_can_add_vectors, 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().

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

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

Definition at line 354 of file continuation_system.C.

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

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

Deactivates the system. Only active systems are solved.

Definition at line 1935 of file system.h.

References libMesh::System::_active.

1936 {
1937  _active = false;
1938 }
void libMesh::ImplicitSystem::disable_cache ( )
virtualinherited

Assembles & solves the linear system Ax=b. Avoids use of any cached data that might affect any solve result. Should be overloaded 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 }
void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited

Definition at line 106 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

107 {
108  _enable_print_counter = false;
109  return;
110 }
virtual bool libMesh::DifferentiablePhysics::element_constraint ( bool  request_jacobian,
DiffContext  
)
inlinevirtualinherited

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

Users may need to reimplement this for their particular PDE.

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

Definition at line 123 of file diff_physics.h.

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

124  {
125  return request_jacobian;
126  }
virtual void libMesh::DifferentiableSystem::element_postprocess ( DiffContext )
inlinevirtualinherited

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

Definition at line 251 of file diff_system.h.

251 {}
virtual void libMesh::DifferentiableQoI::element_qoi ( DiffContext ,
const QoISet  
)
inlinevirtualinherited

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

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

Definition at line 108 of file diff_qoi.h.

110  {}
virtual void libMesh::DifferentiableQoI::element_qoi_derivative ( DiffContext ,
const QoISet  
)
inlinevirtualinherited

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

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

Definition at line 120 of file diff_qoi.h.

122  {}
virtual bool libMesh::DifferentiablePhysics::element_time_derivative ( bool  request_jacobian,
DiffContext  
)
inlinevirtualinherited

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

Users need to reimplement this for their particular PDE.

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

Definition at line 105 of file diff_physics.h.

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

106  {
107  return request_jacobian;
108  }
void libMesh::ReferenceCounter::enable_print_counter_info ( )
staticinherited

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

Definition at line 100 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

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

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

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

Reimplemented from libMesh::DifferentiablePhysics.

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

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

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

Reimplemented in libMesh::FEMPhysics.

Definition at line 213 of file diff_physics.h.

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

214  {
215  return request_jacobian;
216  }
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 827 of file implicit_system.C.

References libMesh::SensitivityData::allocate_data(), libMesh::QoISet::has_index(), libMesh::Real, libMesh::ParameterVector::size(), and libMesh::TOLERANCE.

830 {
831  const unsigned int Np = libmesh_cast_int<unsigned int>
832  (parameters.size());
833  const unsigned int Nq = libmesh_cast_int<unsigned int>
834  (qoi.size());
835 
836  // We currently get partial derivatives via central differencing
837  const Real delta_p = TOLERANCE;
838 
839  // An introduction to the problem:
840  //
841  // Residual R(u(p),p) = 0
842  // partial R / partial u = J = system matrix
843  //
844  // This implies that:
845  // d/dp(R) = 0
846  // (partial R / partial p) +
847  // (partial R / partial u) * (partial u / partial p) = 0
848 
849  // We first solve for (partial u / partial p) for each parameter:
850  // J * (partial u / partial p) = - (partial R / partial p)
851 
852  this->sensitivity_solve(parameters);
853 
854  // Get ready to fill in senstivities:
855  sensitivities.allocate_data(qoi_indices, *this, parameters);
856 
857  // We use the identity:
858  // dq/dp = (partial q / partial p) + (partial q / partial u) *
859  // (partial u / partial p)
860 
861  // We get (partial q / partial u) from the user
862  this->assemble_qoi_derivative(qoi_indices);
863 
864  // FIXME: what do we do with adjoint boundary conditions here?
865 
866  // We don't need these to be closed() in this function, but libMesh
867  // standard practice is to have them closed() by the time the
868  // function exits
869  for (unsigned int i=0; i != this->qoi.size(); ++i)
870  if (qoi_indices.has_index(i))
871  this->get_adjoint_rhs(i).close();
872 
873  for (unsigned int j=0; j != Np; ++j)
874  {
875  // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
876 
877  Number old_parameter = *parameters[j];
878 
879  *parameters[j] = old_parameter - delta_p;
880  this->assemble_qoi();
881  std::vector<Number> qoi_minus = this->qoi;
882 
883  *parameters[j] = old_parameter + delta_p;
884  this->assemble_qoi();
885  std::vector<Number>& qoi_plus = this->qoi;
886 
887  std::vector<Number> partialq_partialp(Nq, 0);
888  for (unsigned int i=0; i != Nq; ++i)
889  if (qoi_indices.has_index(i))
890  partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
891 
892  // Don't leave the parameter changed
893  *parameters[j] = old_parameter;
894 
895  for (unsigned int i=0; i != Nq; ++i)
896  if (qoi_indices.has_index(i))
897  sensitivities[i][j] = partialq_partialp[i] +
898  this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(j));
899  }
900 
901  // All parameters have been reset.
902  // We didn't cache the original rhs or matrix for memory reasons,
903  // but we can restore them to a state consistent solution -
904  // principle of least surprise.
905  this->assembly(true, true);
906  this->rhs->close();
907  this->matrix->close();
908  this->assemble_qoi(qoi_indices);
909 }
NumericVector< Number > & libMesh::System::get_adjoint_rhs ( unsigned int  i = 0)
inherited
Returns
a reference to one of the system's adjoint rhs vectors, by default the one corresponding to the first qoi. This what the user's QoI derivative code should assemble when setting up an adjoint problem

Definition at line 1027 of file system.C.

References libMesh::System::get_vector().

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

1028 {
1029  std::ostringstream adjoint_rhs_name;
1030  adjoint_rhs_name << "adjoint_rhs" << i;
1031 
1032  return this->get_vector(adjoint_rhs_name.str());
1033 }
const NumericVector< Number > & libMesh::System::get_adjoint_rhs ( unsigned int  i = 0) const
inherited
Returns
a reference to one of the system's adjoint rhs vectors, by default the one corresponding to the first qoi.

Definition at line 1037 of file system.C.

References libMesh::System::get_vector().

1038 {
1039  std::ostringstream adjoint_rhs_name;
1040  adjoint_rhs_name << "adjoint_rhs" << i;
1041 
1042  return this->get_vector(adjoint_rhs_name.str());
1043 }
NumericVector< Number > & libMesh::System::get_adjoint_solution ( unsigned int  i = 0)
inherited
Returns
a reference to one of the system's adjoint solution vectors, by default the one corresponding to the first qoi.

Definition at line 967 of file system.C.

References libMesh::System::get_vector().

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

968 {
969  std::ostringstream adjoint_name;
970  adjoint_name << "adjoint_solution" << i;
971 
972  return this->get_vector(adjoint_name.str());
973 }
const NumericVector< Number > & libMesh::System::get_adjoint_solution ( unsigned int  i = 0) const
inherited
Returns
a reference to one of the system's adjoint solution vectors, by default the one corresponding to the first qoi.

Definition at line 977 of file system.C.

References libMesh::System::get_vector().

978 {
979  std::ostringstream adjoint_name;
980  adjoint_name << "adjoint_solution" << i;
981 
982  return this->get_vector(adjoint_name.str());
983 }
void libMesh::System::get_all_variable_numbers ( std::vector< unsigned int > &  all_variable_numbers) const
inherited

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

Definition at line 1260 of file system.C.

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

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

Definition at line 1903 of file system.h.

References libMesh::System::_dof_map.

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), libMesh::__libmesh_petsc_diff_solver_residual(), libMesh::__libmesh_petsc_snes_jacobian(), libMesh::__libmesh_petsc_snes_residual(), libMesh::ExactSolution::_compute_error(), libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::HPCoarsenTest::add_projection(), libMesh::UnsteadySolver::adjoint_advance_timestep(), libMesh::ImplicitSystem::adjoint_solve(), libMesh::UnsteadySolver::advance_timestep(), libMesh::EquationSystems::allgather(), libMesh::EquationSystems::build_discontinuous_solution_vector(), libMesh::EquationSystems::build_solution_vector(), libMesh::System::calculate_norm(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), DMCreateDomainDecomposition_libMesh(), DMCreateFieldDecomposition_libMesh(), DMlibMeshFunction(), DMlibMeshJacobian(), DMLibMeshSetSystem(), libMesh::DofMap::enforce_constraints_exactly(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::for(), libMesh::System::get_info(), libMesh::EquationSystems::get_solution(), libMesh::SystemSubsetBySubdomain::init(), libMesh::UnsteadySolver::init_data(), libMesh::EigenSystem::init_matrices(), libMesh::ImplicitSystem::init_matrices(), libMesh::CondensedEigenSystem::initialize_condensed_dofs(), libMesh::System::local_dof_indices(), libMesh::DofMap::max_constraint_error(), libMesh::DGFEMContext::neighbor_side_fe_reinit(), libMesh::UnsteadySolver::old_nonlinear_solution(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::ProjectSolution::operator()(), libMesh::ProjectFEMSolution::operator()(), libMesh::BoundaryProjectSolution::operator()(), libMesh::ErrorVector::plot_error(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::FEMContext::pre_fe_reinit(), libMesh::System::project_vector(), libMesh::System::re_update(), libMesh::System::read_parallel_data(), libMesh::System::read_SCALAR_dofs(), libMesh::UnsteadySolver::reinit(), libMesh::ImplicitSystem::reinit(), libMesh::EigenSystem::reinit(), libMesh::EquationSystems::reinit(), libMesh::UnsteadySolver::retrieve_timestep(), libMesh::HPCoarsenTest::select_refinement(), libMesh::ImplicitSystem::sensitivity_solve(), libMesh::NewtonSolver::solve(), libMesh::PetscDiffSolver::solve(), libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve(), libMesh::ImplicitSystem::weighted_sensitivity_solve(), libMesh::System::write_parallel_data(), libMesh::EnsightIO::write_scalar_ascii(), libMesh::System::write_SCALAR_dofs(), and libMesh::EnsightIO::write_vector_ascii().

1904 {
1905  return *_dof_map;
1906 }
DofMap & libMesh::System::get_dof_map ( )
inlineinherited
Returns
a writeable reference to this system's _dof_map.

Definition at line 1911 of file system.h.

References libMesh::System::_dof_map.

1912 {
1913  return *_dof_map;
1914 }
EquationSystems& libMesh::System::get_equation_systems ( )
inlineinherited
Returns
a reference to this system's parent EquationSystems object.

Definition at line 686 of file system.h.

References libMesh::System::_equation_systems.

686 { return _equation_systems; }
std::string libMesh::ReferenceCounter::get_info ( )
staticinherited

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

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

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

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

Definition at line 1635 of file system.C.

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

1636 {
1637  std::ostringstream oss;
1638 
1639 
1640  const std::string& sys_name = this->name();
1641 
1642  oss << " System #" << this->number() << ", \"" << sys_name << "\"\n"
1643  << " Type \"" << this->system_type() << "\"\n"
1644  << " Variables=";
1645 
1646  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1647  {
1648  const VariableGroup &vg_description (this->variable_group(vg));
1649 
1650  if (vg_description.n_variables() > 1) oss << "{ ";
1651  for (unsigned int vn=0; vn<vg_description.n_variables(); vn++)
1652  oss << "\"" << vg_description.name(vn) << "\" ";
1653  if (vg_description.n_variables() > 1) oss << "} ";
1654  }
1655 
1656  oss << '\n';
1657 
1658  oss << " Finite Element Types=";
1659 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
1660  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1661  oss << "\""
1662  << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
1663  << "\" ";
1664 #else
1665  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1666  {
1667  oss << "\""
1668  << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().family)
1669  << "\", \""
1670  << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_group(vg).type().radial_family)
1671  << "\" ";
1672  }
1673 
1674  oss << '\n' << " Infinite Element Mapping=";
1675  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1676  oss << "\""
1677  << Utility::enum_to_string<InfMapType>(this->get_dof_map().variable_group(vg).type().inf_map)
1678  << "\" ";
1679 #endif
1680 
1681  oss << '\n';
1682 
1683  oss << " Approximation Orders=";
1684  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
1685  {
1686 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
1687  oss << "\""
1688  << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
1689  << "\" ";
1690 #else
1691  oss << "\""
1692  << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().order)
1693  << "\", \""
1694  << Utility::enum_to_string<Order>(this->get_dof_map().variable_group(vg).type().radial_order)
1695  << "\" ";
1696 #endif
1697  }
1698 
1699  oss << '\n';
1700 
1701  oss << " n_dofs()=" << this->n_dofs() << '\n';
1702  oss << " n_local_dofs()=" << this->n_local_dofs() << '\n';
1703 #ifdef LIBMESH_ENABLE_CONSTRAINTS
1704  oss << " n_constrained_dofs()=" << this->n_constrained_dofs() << '\n';
1705  oss << " n_local_constrained_dofs()=" << this->n_local_constrained_dofs() << '\n';
1706 #endif
1707 
1708  oss << " " << "n_vectors()=" << this->n_vectors() << '\n';
1709  oss << " " << "n_matrices()=" << this->n_matrices() << '\n';
1710 // oss << " " << "n_additional_matrices()=" << this->n_additional_matrices() << '\n';
1711 
1712  oss << this->get_dof_map().get_info();
1713 
1714  return oss.str();
1715 }
std::pair< unsigned int, Real > libMesh::DifferentiableSystem::get_linear_solve_parameters ( ) const
virtualinherited

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

Reimplemented from libMesh::ImplicitSystem.

Definition at line 164 of file diff_system.C.

References libMesh::libmesh_assert(), and libMesh::DifferentiableSystem::time_solver.

165 {
167  libmesh_assert_equal_to (&(time_solver->system()), this);
168  return std::make_pair(this->time_solver->diff_solver()->max_linear_iterations,
169  this->time_solver->diff_solver()->relative_residual_tolerance);
170 }
LinearSolver< Number > * libMesh::DifferentiableSystem::get_linear_solver ( ) const
virtualinherited

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

Reimplemented from libMesh::ImplicitSystem.

Definition at line 155 of file diff_system.C.

References libMesh::libmesh_assert(), and libMesh::DifferentiableSystem::time_solver.

156 {
158  libmesh_assert_equal_to (&(time_solver->system()), this);
159  return this->time_solver->linear_solver().get();
160 }
const SparseMatrix< Number > & libMesh::ImplicitSystem::get_matrix ( const std::string &  mat_name) const
inherited
Returns
a const reference to this system's additional matrix named mat_name. None of these matrices is involved in the solution process. Access is only granted when the matrix is already properly initialized.

Definition at line 258 of file implicit_system.C.

References libMesh::ImplicitSystem::_matrices, and libMesh::err.

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

259 {
260  // Make sure the matrix exists
261  const_matrices_iterator pos = _matrices.find (mat_name);
262 
263  if (pos == _matrices.end())
264  {
265  libMesh::err << "ERROR: matrix "
266  << mat_name
267  << " does not exist in this system!"
268  << std::endl;
269  libmesh_error();
270  }
271 
272  return *(pos->second);
273 }
SparseMatrix< Number > & libMesh::ImplicitSystem::get_matrix ( const std::string &  mat_name)
inherited
Returns
a writeable reference to this system's additional matrix named mat_name. None of these matrices is involved in the solution process. Access is only granted when the matrix is already properly initialized.

Definition at line 277 of file implicit_system.C.

References libMesh::ImplicitSystem::_matrices, and libMesh::err.

278 {
279  // Make sure the matrix exists
280  matrices_iterator pos = _matrices.find (mat_name);
281 
282  if (pos == _matrices.end())
283  {
284  libMesh::err << "ERROR: matrix "
285  << mat_name
286  << " does not exist in this system!"
287  << std::endl;
288  libmesh_error();
289  }
290 
291  return *(pos->second);
292 }
const MeshBase & libMesh::System::get_mesh ( ) const
inlineinherited
Returns
a constant reference to this systems's _mesh.

Definition at line 1887 of file system.h.

References libMesh::System::_mesh.

Referenced by libMesh::ExactSolution::_compute_error(), libMesh::HPCoarsenTest::add_projection(), libMesh::FEMSystem::assemble_qoi(), libMesh::FEMSystem::assemble_qoi_derivative(), libMesh::FEMSystem::assembly(), libMesh::System::calculate_norm(), DMCreateDomainDecomposition_libMesh(), DMCreateFieldDecomposition_libMesh(), DMLibMeshSetSystem(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointResidualErrorEstimator::estimate_error(), libMesh::SystemSubsetBySubdomain::init(), libMesh::System::init_data(), libMesh::EigenSystem::init_matrices(), libMesh::ImplicitSystem::init_matrices(), libMesh::System::local_dof_indices(), libMesh::DofMap::max_constraint_error(), libMesh::FEMSystem::mesh_position_get(), libMesh::FEMSystem::mesh_position_set(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::ProjectSolution::operator()(), libMesh::BoundaryProjectSolution::operator()(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::FEMSystem::postprocess(), libMesh::System::project_vector(), libMesh::System::read_header(), libMesh::System::read_legacy_data(), libMesh::System::read_parallel_data(), libMesh::System::read_serialized_vector(), libMesh::System::read_serialized_vectors(), libMesh::ImplicitSystem::reinit(), libMesh::EigenSystem::reinit(), libMesh::HPSingularity::select_refinement(), libMesh::HPCoarsenTest::select_refinement(), libMesh::System::write_header(), libMesh::System::write_parallel_data(), libMesh::System::write_serialized_vector(), libMesh::System::write_serialized_vectors(), and libMesh::System::zero_variable().

1888 {
1889  return _mesh;
1890 }
MeshBase & libMesh::System::get_mesh ( )
inlineinherited
Returns
a reference to this systems's _mesh.

Definition at line 1895 of file system.h.

References libMesh::System::_mesh.

1896 {
1897  return _mesh;
1898 }
const System * libMesh::DifferentiablePhysics::get_mesh_system ( ) const
inlineinherited

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

Definition at line 414 of file diff_physics.h.

References libMesh::DifferentiablePhysics::_mesh_sys.

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

415 {
416  return _mesh_sys;
417 }
System * libMesh::DifferentiablePhysics::get_mesh_system ( )
inlineinherited

Returns a reference to the system with variables corresponding to mesh nodal coordinates, or NULL if the mesh is fixed.

Definition at line 420 of file diff_physics.h.

References libMesh::DifferentiablePhysics::_mesh_sys.

421 {
422  return _mesh_sys;
423 }
unsigned int libMesh::DifferentiablePhysics::get_mesh_x_var ( ) const
inlineinherited

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

Definition at line 426 of file diff_physics.h.

References libMesh::DifferentiablePhysics::_mesh_x_var.

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

427 {
428  return _mesh_x_var;
429 }
unsigned int libMesh::DifferentiablePhysics::get_mesh_y_var ( ) const
inlineinherited

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

Definition at line 432 of file diff_physics.h.

References libMesh::DifferentiablePhysics::_mesh_y_var.

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

433 {
434  return _mesh_y_var;
435 }
unsigned int libMesh::DifferentiablePhysics::get_mesh_z_var ( ) const
inlineinherited

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

Definition at line 438 of file diff_physics.h.

References libMesh::DifferentiablePhysics::_mesh_z_var.

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

439 {
440  return _mesh_z_var;
441 }
const DifferentiablePhysics* libMesh::DifferentiableSystem::get_physics ( ) const
inlineinherited

Returns const reference to DifferentiablePhysics object. Note that if no external Physics object is attached, the default is this.

Definition at line 163 of file diff_system.h.

References libMesh::DifferentiableSystem::_diff_physics.

Referenced by libMesh::FEMSystem::build_context(), and libMesh::FEMSystem::init_context().

164  { return this->_diff_physics; }
DifferentiablePhysics* libMesh::DifferentiableSystem::get_physics ( )
inlineinherited

Returns reference to DifferentiablePhysics object. Note that if no external Physics object is attached, the default is this.

Definition at line 170 of file diff_system.h.

References libMesh::DifferentiableSystem::_diff_physics.

171  { return this->_diff_physics; }
const DifferentiableQoI* libMesh::DifferentiableSystem::get_qoi ( ) const
inlineinherited

Returns const reference to DifferentiableQoI object. Note that if no external QoI object is attached, the default is this.

Definition at line 184 of file diff_system.h.

References libMesh::DifferentiableSystem::diff_qoi.

185  { return this->diff_qoi; }
DifferentiableQoI* libMesh::DifferentiableSystem::get_qoi ( )
inlineinherited

Returns reference to DifferentiableQoI object. Note that if no external QoI object is attached, the default is this.

Definition at line 191 of file diff_system.h.

References libMesh::DifferentiableSystem::diff_qoi.

192  { return this->diff_qoi; }
NumericVector< Number > & libMesh::System::get_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. By default these vectors are built by the library, using finite differences, when assemble_residual_derivatives() is called.

When assembled, this vector should hold -(partial R / partial p_i)

Definition at line 1057 of file system.C.

References libMesh::System::get_vector().

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

1058 {
1059  std::ostringstream sensitivity_rhs_name;
1060  sensitivity_rhs_name << "sensitivity_rhs" << i;
1061 
1062  return this->get_vector(sensitivity_rhs_name.str());
1063 }
const NumericVector< Number > & libMesh::System::get_sensitivity_rhs ( unsigned int  i = 0) const
inherited
Returns
a reference to one of the system's sensitivity rhs vectors, by default the one corresponding to the first parameter.

Definition at line 1067 of file system.C.

References libMesh::System::get_vector().

1068 {
1069  std::ostringstream sensitivity_rhs_name;
1070  sensitivity_rhs_name << "sensitivity_rhs" << i;
1071 
1072  return this->get_vector(sensitivity_rhs_name.str());
1073 }
NumericVector< Number > & libMesh::System::get_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.

Definition at line 916 of file system.C.

References libMesh::System::get_vector().

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

917 {
918  std::ostringstream sensitivity_name;
919  sensitivity_name << "sensitivity_solution" << i;
920 
921  return this->get_vector(sensitivity_name.str());
922 }
const NumericVector< Number > & libMesh::System::get_sensitivity_solution ( unsigned int  i = 0) const
inherited
Returns
a reference to one of the system's solution sensitivity vectors, by default the one corresponding to the first parameter.

Definition at line 926 of file system.C.

References libMesh::System::get_vector().

927 {
928  std::ostringstream sensitivity_name;
929  sensitivity_name << "sensitivity_solution" << i;
930 
931  return this->get_vector(sensitivity_name.str());
932 }
TimeSolver & libMesh::DifferentiableSystem::get_time_solver ( )
inlineinherited

Returns a pointer to the time solver attached to the calling system

Definition at line 328 of file diff_system.h.

References libMesh::libmesh_assert(), and libMesh::DifferentiableSystem::time_solver.

Referenced by libMesh::DifferentiableSystem::adjoint_solve(), libMesh::FEMSystem::build_context(), libMesh::FEMSystem::postprocess(), and libMesh::DifferentiableSystem::solve().

329 {
331  libmesh_assert_equal_to (&(time_solver->system()), this);
332  return *time_solver;
333 }
const TimeSolver & libMesh::DifferentiableSystem::get_time_solver ( ) const
inlineinherited

Non-const version of the above

Definition at line 336 of file diff_system.h.

References libMesh::libmesh_assert(), and libMesh::DifferentiableSystem::time_solver.

337 {
339  libmesh_assert_equal_to (&(time_solver->system()), this);
340  return *time_solver;
341 }
const NumericVector< Number > & libMesh::System::get_vector ( const std::string &  vec_name) const
inherited
Returns
a const reference to this system's additional vector named vec_name. Access is only granted when the vector is already properly initialized.

Definition at line 785 of file system.C.

References libMesh::System::_vectors, and libMesh::err.

Referenced by libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::UnsteadySolver::adjoint_advance_timestep(), libMesh::AdaptiveTimeSolver::advance_timestep(), libMesh::UnsteadySolver::advance_timestep(), libMesh::System::compare(), libMesh::UnsteadySolver::du(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::System::get_adjoint_rhs(), libMesh::System::get_adjoint_solution(), libMesh::System::get_sensitivity_rhs(), libMesh::System::get_sensitivity_solution(), libMesh::System::get_weighted_sensitivity_adjoint_solution(), libMesh::System::get_weighted_sensitivity_solution(), libMesh::NewmarkSystem::initial_conditions(), libMesh::UnsteadySolver::reinit(), libMesh::MemorySolutionHistory::retrieve(), libMesh::UnsteadySolver::retrieve_timestep(), libMesh::TwostepTimeSolver::solve(), libMesh::FrequencySystem::solve(), libMesh::NewmarkSystem::update_rhs(), and libMesh::NewmarkSystem::update_u_v_a().

786 {
787  // Make sure the vector exists
788  const_vectors_iterator pos = _vectors.find(vec_name);
789 
790  if (pos == _vectors.end())
791  {
792  libMesh::err << "ERROR: vector "
793  << vec_name
794  << " does not exist in this system!"
795  << std::endl;
796  libmesh_error();
797  }
798 
799  return *(pos->second);
800 }
NumericVector< Number > & libMesh::System::get_vector ( const std::string &  vec_name)
inherited
Returns
a writeable reference to this system's additional vector named vec_name. Access is only granted when the vector is already properly initialized.

Definition at line 804 of file system.C.

References libMesh::System::_vectors, and libMesh::err.

805 {
806  // Make sure the vector exists
807  vectors_iterator pos = _vectors.find(vec_name);
808 
809  if (pos == _vectors.end())
810  {
811  libMesh::err << "ERROR: vector "
812  << vec_name
813  << " does not exist in this system!"
814  << std::endl;
815  libmesh_error();
816  }
817 
818  return *(pos->second);
819 }
const NumericVector< Number > & libMesh::System::get_vector ( const unsigned int  vec_num) const
inherited
Returns
a const reference to this system's additional vector number vec_num (where the vectors are counted starting with 0).

Definition at line 823 of file system.C.

References libMesh::libmesh_assert(), libMesh::System::vectors_begin(), and libMesh::System::vectors_end().

824 {
827  unsigned int num = 0;
828  while((num<vec_num) && (v!=v_end))
829  {
830  num++;
831  ++v;
832  }
833  libmesh_assert (v != v_end);
834  return *(v->second);
835 }
NumericVector< Number > & libMesh::System::get_vector ( const unsigned int  vec_num)
inherited
Returns
a writeable reference to this system's additional vector number vec_num (where the vectors are counted starting with 0).

Definition at line 839 of file system.C.

References libMesh::libmesh_assert(), libMesh::System::vectors_begin(), and libMesh::System::vectors_end().

840 {
842  vectors_iterator v_end = vectors_end();
843  unsigned int num = 0;
844  while((num<vec_num) && (v!=v_end))
845  {
846  num++;
847  ++v;
848  }
849  libmesh_assert (v != v_end);
850  return *(v->second);
851 }
NumericVector< Number > & libMesh::System::get_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.

Definition at line 997 of file system.C.

References libMesh::System::get_vector().

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

998 {
999  std::ostringstream adjoint_name;
1000  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1001 
1002  return this->get_vector(adjoint_name.str());
1003 }
const NumericVector< Number > & libMesh::System::get_weighted_sensitivity_adjoint_solution ( unsigned int  i = 0) const
inherited
Returns
a reference to one of the system's weighted sensitivity adjoint solution vectors, by default the one corresponding to the first qoi.

Definition at line 1007 of file system.C.

References libMesh::System::get_vector().

1008 {
1009  std::ostringstream adjoint_name;
1010  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
1011 
1012  return this->get_vector(adjoint_name.str());
1013 }
NumericVector< Number > & libMesh::System::get_weighted_sensitivity_solution ( )
inherited
Returns
a reference to the solution of the last weighted sensitivity solve

Definition at line 943 of file system.C.

References libMesh::System::get_vector().

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

944 {
945  return this->get_vector("weighted_sensitivity_solution");
946 }
const NumericVector< Number > & libMesh::System::get_weighted_sensitivity_solution ( ) const
inherited
Returns
a reference to the solution of the last weighted sensitivity solve

Definition at line 950 of file system.C.

References libMesh::System::get_vector().

951 {
952  return this->get_vector("weighted_sensitivity_solution");
953 }
bool libMesh::System::has_variable ( const std::string &  var) const
inherited
Returns
true if a variable named var exists in this System

Definition at line 1233 of file system.C.

References libMesh::System::_variable_numbers.

Referenced by libMesh::GMVIO::copy_nodal_solution().

1234 {
1235  return _variable_numbers.count(var);
1236 }
bool libMesh::ImplicitSystem::have_matrix ( const std::string &  mat_name) const
inlineinherited
Returns
true if this System has a matrix associated with the given name, false otherwise.

Definition at line 378 of file implicit_system.h.

References libMesh::ImplicitSystem::_matrices.

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

379 {
380  return (_matrices.count(mat_name));
381 }
bool libMesh::System::have_vector ( const std::string &  vec_name) const
inlineinherited
Returns
true if this System has a vector associated with the given name, false otherwise.

Definition at line 2071 of file system.h.

References libMesh::System::_vectors.

Referenced by libMesh::System::add_vector(), and libMesh::System::remove_vector().

2072 {
2073  return (_vectors.count(vec_name));
2074 }
bool libMesh::System::identify_variable_groups ( ) const
inlineinherited
Returns
true when VariableGroup structures should be automatically identified, false otherwise.

Definition at line 2047 of file system.h.

References libMesh::System::_identify_variable_groups.

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

2048 {
2050 }
void libMesh::System::identify_variable_groups ( const bool  ivg)
inlineinherited

Toggle automatic VariableGroup identification.

Definition at line 2055 of file system.h.

References libMesh::System::_identify_variable_groups.

2056 {
2058 }
void libMesh::ReferenceCounter::increment_constructor_count ( const std::string &  name)
inlineprotectedinherited

Increments the construction counter. Should be called in the constructor of any derived class that will be reference counted.

Definition at line 163 of file reference_counter.h.

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

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::ReferenceCountedObject().

164 {
165  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
166  std::pair<unsigned int, unsigned int>& p = _counts[name];
167 
168  p.first++;
169 }
void libMesh::ReferenceCounter::increment_destructor_count ( const std::string &  name)
inlineprotectedinherited

Increments the destruction counter. Should be called in the destructor of any derived class that will be reference counted.

Definition at line 176 of file reference_counter.h.

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

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::~ReferenceCountedObject().

177 {
178  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
179  std::pair<unsigned int, unsigned int>& p = _counts[name];
180 
181  p.second++;
182 }
void libMesh::System::init ( )
inherited

Initializes degrees of freedom on the current mesh. Sets the

Definition at line 226 of file system.C.

References libMesh::System::_basic_system_only, libMesh::System::init_data(), libMesh::System::n_vars(), and libMesh::System::user_initialization().

227 {
228  // First initialize any required data:
229  // either only the basic System data
230  if (_basic_system_only)
232  // or all the derived class' data too
233  else
234  this->init_data();
235 
236  // If no variables have been added to this system
237  // don't do anything
238  if(!this->n_vars())
239  return;
240 
241  // Then call the user-provided intialization function
242  this->user_initialization();
243 }
void libMesh::FEMSystem::init_context ( DiffContext c)
virtualinherited

Reimplemented from libMesh::DifferentiablePhysics.

Definition at line 977 of file fem_system.C.

References libMesh::DifferentiableSystem::deltat, libMesh::FEInterface::field_type(), libMesh::FEMContext::get_element_fe(), libMesh::FEAbstract::get_JxW(), libMesh::FEGenericBase< T >::get_phi(), libMesh::DifferentiableSystem::get_physics(), libMesh::DifferentiablePhysics::is_time_evolving(), libMesh::System::n_vars(), libMesh::DiffContext::set_deltat_pointer(), libMeshEnums::TYPE_SCALAR, libMeshEnums::TYPE_VECTOR, and libMesh::System::variable_type().

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

978 {
979  // Parent::init_context(c); // may be a good idea in derived classes
980 
981  // Although we do this in DiffSystem::build_context() and
982  // FEMSystem::build_context() as well, we do it here just to be
983  // extra sure that the deltat pointer gets set. Since the
984  // intended behavior is for classes derived from FEMSystem to
985  // call Parent::init_context() in their own init_context()
986  // overloads, we can ensure that those classes get the correct
987  // deltat pointers even if they have different build_context()
988  // overloads.
989  c.set_deltat_pointer ( &deltat );
990 
991  FEMContext &context = libmesh_cast_ref<FEMContext&>(c);
992 
993  // Make sure we're prepared to do mass integration
994  for (unsigned int var = 0; var != this->n_vars(); ++var)
995  if (this->get_physics()->is_time_evolving(var))
996  {
997  // Request shape functions based on FEType
998  switch( FEInterface::field_type( this->variable_type(var) ) )
999  {
1000  case( TYPE_SCALAR ):
1001  {
1002  FEBase* elem_fe = NULL;
1003  context.get_element_fe(var, elem_fe);
1004  elem_fe->get_JxW();
1005  elem_fe->get_phi();
1006  }
1007  break;
1008  case( TYPE_VECTOR ):
1009  {
1010  FEGenericBase<RealGradient>* elem_fe = NULL;
1011  context.get_element_fe(var, elem_fe);
1012  elem_fe->get_JxW();
1013  elem_fe->get_phi();
1014  }
1015  break;
1016  default:
1017  {
1018  libmesh_error();
1019  }
1020  }
1021  }
1022 }
void libMesh::ContinuationSystem::init_data ( )
protectedvirtual

Initializes the member data fields associated with the system, so that, e.g., assemble() may be used.

Reimplemented from libMesh::FEMSystem.

Definition at line 86 of file continuation_system.C.

References libMesh::System::add_vector(), delta_u, du_ds, libMesh::FEMSystem::init_data(), previous_du_ds, previous_u, y, y_old, and z.

87 {
88  // Add a vector which stores the tangent "du/ds" to the system and save its pointer.
89  du_ds = &(add_vector("du_ds"));
90 
91  // Add a vector which stores the tangent "du/ds" to the system and save its pointer.
92  previous_du_ds = &(add_vector("previous_du_ds"));
93 
94  // Add a vector to keep track of the previous nonlinear solution
95  // at the old value of lambda.
96  previous_u = &(add_vector("previous_u"));
97 
98  // Add a vector to keep track of the temporary solution "y" of Ay=G_{\lambda}.
99  y = &(add_vector("y"));
100 
101  // Add a vector to keep track of the "old value" of "y" which is the solution of Ay=G_{\lambda}.
102  y_old = &(add_vector("y_old"));
103 
104  // Add a vector to keep track of the temporary solution "z" of Az=-G.
105  z = &(add_vector("z"));
106 
107  // Add a vector to keep track of the Newton update during the constrained PDE solves.
108  delta_u = &(add_vector("delta_u"));
109 
110  // Call the Parent's initialization routine.
112 }
void libMesh::ImplicitSystem::init_matrices ( )
protectedvirtualinherited

Initializes the matrices associated with this system.

Definition at line 105 of file implicit_system.C.

References libMesh::ImplicitSystem::_can_add_matrices, libMesh::ImplicitSystem::_matrices, libMesh::DofMap::attach_matrix(), libMesh::DofMap::compute_sparsity(), libMesh::dof_map, libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::SparseMatrix< T >::initialized(), libMesh::DofMap::is_attached(), libMesh::libmesh_assert(), and libMesh::ImplicitSystem::matrix.

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

106 {
108 
109  // Check for quick return in case the system matrix
110  // (and by extension all the matrices) has already
111  // been initialized
112  if (matrix->initialized())
113  return;
114 
115  // Get a reference to the DofMap
116  DofMap& dof_map = this->get_dof_map();
117 
118  // no chance to add other matrices
119  _can_add_matrices = false;
120 
121  // Tell the matrices about the dof map, and vice versa
122  for (matrices_iterator pos = _matrices.begin();
123  pos != _matrices.end(); ++pos)
124  {
125  SparseMatrix<Number> &m = *(pos->second);
126  libmesh_assert (!m.initialized());
127 
128  // We want to allow repeated init() on systems, but we don't
129  // want to attach the same matrix to the DofMap twice
130  if (!dof_map.is_attached(m))
131  dof_map.attach_matrix (m);
132  }
133 
134  // Compute the sparsity pattern for the current
135  // mesh and DOF distribution. This also updates
136  // additional matrices, \p DofMap now knows them
137  dof_map.compute_sparsity (this->get_mesh());
138 
139  // Initialize matrices
140  for (matrices_iterator pos = _matrices.begin();
141  pos != _matrices.end(); ++pos)
142  pos->second->init ();
143 
144  // Set the additional matrices to 0.
145  for (matrices_iterator pos = _matrices.begin();
146  pos != _matrices.end(); ++pos)
147  pos->second->zero ();
148 }
virtual void libMesh::DifferentiablePhysics::init_physics ( const System sys)
virtualinherited

Initialize any data structures associated with the physics.

Referenced by libMesh::DifferentiableSystem::attach_physics(), and libMesh::DifferentiableSystem::init_data().

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

Initialize system qoi. By default, does nothing in order to maintain backward compatibility for FEMSystem applications that control qoi.

Definition at line 69 of file diff_qoi.h.

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

69 {}
void libMesh::ContinuationSystem::initialize_tangent ( )
private

Before starting arclength continuation, we need at least 2 prior solutions (both solution and u_previous should be filled with meaningful values) And we need to initialize the tangent vector. This only needs to be called once.

Definition at line 127 of file continuation_system.C.

References libMesh::NumericVector< T >::add(), libMesh::NumericVector< T >::close(), continuation_parameter, dlambda_ds, ds_current, du_ds, libMesh::NumericVector< T >::l2_norm(), libMesh::libmesh_assert(), old_continuation_parameter, libMesh::out, previous_u, quiet, libMesh::Real, libMesh::NumericVector< T >::scale(), set_Theta(), libMesh::System::solution, solve_tangent(), tangent_initialized, Theta, Theta_LOCA, update_solution(), and y.

Referenced by continuation_solve().

128 {
129  // Be sure the tangent was not already initialized.
131 
132  // Compute delta_s_zero, the initial arclength travelled during the
133  // first step. Here we assume that previous_u and lambda_old store
134  // the previous solution and control parameter. You may need to
135  // read in an old solution (or solve the non-continuation system)
136  // first and call save_current_solution() before getting here.
137 
138  // 1.) Compute delta_s_zero as ||u|| - ||u_old|| + ...
139  // Compute norms of the current and previous solutions
140 // Real norm_u = solution->l2_norm();
141 // Real norm_previous_u = previous_u->l2_norm();
142 
143 // if (!quiet)
144 // {
145 // libMesh::out << "norm_u=" << norm_u << std::endl;
146 // libMesh::out << "norm_previous_u=" << norm_previous_u << std::endl;
147 // }
148 
149 // if (norm_u == norm_previous_u)
150 // {
151 // libMesh::err << "Warning, it appears u and previous_u are the "
152 // << "same, are you sure this is correct?"
153 // << "It's possible you forgot to set one or the other..."
154 // << std::endl;
155 // }
156 
157 // Real delta_s_zero = std::sqrt(
158 // (norm_u - norm_previous_u)*(norm_u - norm_previous_u) +
159 // (*continuation_parameter-old_continuation_parameter)*
160 // (*continuation_parameter-old_continuation_parameter)
161 // );
162 
163 // // 2.) Compute delta_s_zero as ||u -u_old|| + ...
164 // *delta_u = *solution;
165 // delta_u->add(-1., *previous_u);
166 // delta_u->close();
167 // Real norm_delta_u = delta_u->l2_norm();
168 // Real norm_u = solution->l2_norm();
169 // Real norm_previous_u = previous_u->l2_norm();
170 
171 // // Scale norm_delta_u by the bigger of either norm_u or norm_previous_u
172 // norm_delta_u /= std::max(norm_u, norm_previous_u);
173 
174 // if (!quiet)
175 // {
176 // libMesh::out << "norm_u=" << norm_u << std::endl;
177 // libMesh::out << "norm_previous_u=" << norm_previous_u << std::endl;
178 // //libMesh::out << "norm_delta_u=" << norm_delta_u << std::endl;
179 // libMesh::out << "norm_delta_u/max(|u|,|u_old|)=" << norm_delta_u << std::endl;
180 // libMesh::out << "|norm_u-norm_previous_u|=" << std::abs(norm_u - norm_previous_u) << std::endl;
181 // }
182 
183 // const Real dlambda = *continuation_parameter-old_continuation_parameter;
184 
185 // if (!quiet)
186 // libMesh::out << "dlambda=" << dlambda << std::endl;
187 
188 // Real delta_s_zero = std::sqrt(
189 // (norm_delta_u*norm_delta_u) +
190 // (dlambda*dlambda)
191 // );
192 
193 // if (!quiet)
194 // libMesh::out << "delta_s_zero=" << delta_s_zero << std::endl;
195 
196  // 1.) + 2.)
197 // // Now approximate the initial tangent d(lambda)/ds
198 // this->dlambda_ds = (*continuation_parameter-old_continuation_parameter) / delta_s_zero;
199 
200 
201 // // We can also approximate the deriv. wrt s by finite differences:
202 // // du/ds = (u1 - u0) / delta_s_zero.
203 // // FIXME: Use delta_u from above if we decide to keep that method.
204 // *du_ds = *solution;
205 // du_ds->add(-1., *previous_u);
206 // du_ds->scale(1./delta_s_zero);
207 // du_ds->close();
208 
209 
210  // 3.) Treating (u-previous_u)/(lambda - lambda_old) as an approximation to du/d(lambda),
211  // we follow the same technique as Carnes and Shadid.
212 // const Real dlambda = *continuation_parameter-old_continuation_parameter;
213 // libmesh_assert_greater (dlambda, 0.);
214 
215 // // Use delta_u for temporary calculation of du/d(lambda)
216 // *delta_u = *solution;
217 // delta_u->add(-1., *previous_u);
218 // delta_u->scale(1. / dlambda);
219 // delta_u->close();
220 
221 // // Determine initial normalization parameter
222 // const Real solution_size = std::max(solution->l2_norm(), previous_u->l2_norm());
223 // if (solution_size > 1.)
224 // {
225 // Theta = 1./solution_size;
226 
227 // if (!quiet)
228 // libMesh::out << "Setting Normalization Parameter Theta=" << Theta << std::endl;
229 // }
230 
231 // // Compute d(lambda)/ds
232 // // The correct sign of d(lambda)/ds should be positive, since we assume that (lambda > lambda_old)
233 // // but we could always double-check that as well.
234 // Real norm_delta_u = delta_u->l2_norm();
235 // this->dlambda_ds = 1. / std::sqrt(1. + Theta*Theta*norm_delta_u*norm_delta_u);
236 
237 // // Finally, compute du/ds = d(lambda)/ds * du/d(lambda)
238 // *du_ds = *delta_u;
239 // du_ds->scale(dlambda_ds);
240 // du_ds->close();
241 
242 
243  // 4.) Use normalized arclength formula to estimate delta_s_zero
244 // // Determine initial normalization parameter
245 // set_Theta();
246 
247 // // Compute (normalized) delta_s_zero
248 // *delta_u = *solution;
249 // delta_u->add(-1., *previous_u);
250 // delta_u->close();
251 // Real norm_delta_u = delta_u->l2_norm();
252 
253 // const Real dlambda = *continuation_parameter-old_continuation_parameter;
254 
255 // if (!quiet)
256 // libMesh::out << "dlambda=" << dlambda << std::endl;
257 
258 // Real delta_s_zero = std::sqrt(
259 // (Theta_LOCA*Theta_LOCA*Theta*norm_delta_u*norm_delta_u) +
260 // (dlambda*dlambda)
261 // );
262 // *du_ds = *delta_u;
263 // du_ds->scale(1./delta_s_zero);
264 // dlambda_ds = dlambda / delta_s_zero;
265 
266 // if (!quiet)
267 // {
268 // libMesh::out << "delta_s_zero=" << delta_s_zero << std::endl;
269 // libMesh::out << "initial d(lambda)/ds|_0 = " << dlambda_ds << std::endl;
270 // libMesh::out << "initial ||du_ds||_0 = " << du_ds->l2_norm() << std::endl;
271 // }
272 
273 // // FIXME: Also store the initial finite-differenced approximation to -du/dlambda as y.
274 // // We stick to the convention of storing negative y, since that is what we typically
275 // // solve for anyway.
276 // *y = *delta_u;
277 // y->scale(-1./dlambda);
278 // y->close();
279 
280 
281 
282  // 5.) Assume dlambda/ds_0 ~ 1/sqrt(2) and determine the value of Theta_LOCA which
283  // will satisfy this criterion
284 
285  // Initial change in parameter
287  libmesh_assert_not_equal_to (dlambda, 0.0);
288 
289  // Ideal initial value of dlambda_ds
290  dlambda_ds = 1. / std::sqrt(2.);
291  if (dlambda < 0.)
292  dlambda_ds *= -1.;
293 
294  // This also implies the initial value of ds
295  ds_current = dlambda / dlambda_ds;
296 
297  if (!quiet)
298  libMesh::out << "Setting ds_current|_0=" << ds_current << std::endl;
299 
300  // Set y = -du/dlambda using finite difference approximation
301  *y = *solution;
302  y->add(-1., *previous_u);
303  y->scale(-1./dlambda);
304  y->close();
305  const Real ynorm=y->l2_norm();
306 
307  // Finally, set the value of du_ds to be used in the upcoming
308  // tangent calculation. du/ds = du/dlambda * dlambda/ds
309  *du_ds = *y;
311  du_ds->close();
312 
313  // Determine additional solution normalization parameter
314  // (Since we just set du/ds, it will be: ||du||*||du/ds||)
315  set_Theta();
316 
317  // The value of Theta_LOCA which makes dlambda_ds = 1/sqrt(2),
318  // assuming our Theta = ||du||^2.
319  // Theta_LOCA = std::abs(dlambda);
320 
321  // Assuming general Theta
322  Theta_LOCA = std::sqrt(1./Theta/ynorm/ynorm);
323 
324 
325  if (!quiet)
326  {
327  libMesh::out << "Setting initial Theta_LOCA = " << Theta_LOCA << std::endl;
328  libMesh::out << "Theta_LOCA^2*Theta = " << Theta_LOCA*Theta_LOCA*Theta << std::endl;
329  libMesh::out << "initial d(lambda)/ds|_0 = " << dlambda_ds << std::endl;
330  libMesh::out << "initial ||du_ds||_0 = " << du_ds->l2_norm() << std::endl;
331  }
332 
333 
334 
335  // OK, we estimated the tangent at point u0.
336  // Now, to estimate the tangent at point u1, we call the solve_tangent routine.
337 
338  // Set the flag which tells us the method has been initialized.
339  tangent_initialized = true;
340 
341  solve_tangent();
342 
343  // Advance the solution and the parameter to the next value.
344  update_solution();
345 }
bool libMesh::System::is_adjoint_already_solved ( ) const
inlineinherited

Accessor for the adjoint_already_solved boolean

Definition at line 361 of file system.h.

References libMesh::System::adjoint_already_solved.

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

362  { return adjoint_already_solved;}
bool libMesh::DifferentiablePhysics::is_time_evolving ( unsigned int  var) const
inlineinherited

Returns true iff variable var is evolving with respect to time. In general, the user's init() function should have set time_evolving() for any variables which behave like du/dt = F(u), and should not call time_evolving() for any variables which behave like 0 = G(u).

Definition at line 201 of file diff_physics.h.

References libMesh::DifferentiablePhysics::_time_evolving.

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

201  {
202  return _time_evolving[var];
203  }
void libMesh::System::local_dof_indices ( const unsigned int  var,
std::set< dof_id_type > &  var_indices 
) const
inherited

Fills the std::set with the degrees of freedom on the local processor corresponding the the variable number passed in.

Definition at line 1279 of file system.C.

References libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::DofMap::dof_indices(), libMesh::DofMap::end_dof(), libMesh::DofMap::first_dof(), libMesh::System::get_dof_map(), and libMesh::System::get_mesh().

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

1281 {
1282  // Make sure the set is clear
1283  var_indices.clear();
1284 
1285  std::vector<dof_id_type> dof_indices;
1286 
1287  // Begin the loop over the elements
1288  MeshBase::const_element_iterator el =
1290  const MeshBase::const_element_iterator end_el =
1292 
1293  const dof_id_type
1294  first_local = this->get_dof_map().first_dof(),
1295  end_local = this->get_dof_map().end_dof();
1296 
1297  for ( ; el != end_el; ++el)
1298  {
1299  const Elem* elem = *el;
1300  this->get_dof_map().dof_indices (elem, dof_indices, var);
1301 
1302  for(unsigned int i=0; i<dof_indices.size(); i++)
1303  {
1304  dof_id_type dof = dof_indices[i];
1305 
1306  //If the dof is owned by the local processor
1307  if(first_local <= dof && dof < end_local)
1308  var_indices.insert(dof_indices[i]);
1309  }
1310  }
1311 }
virtual bool libMesh::FEMPhysics::mass_residual ( bool  request_jacobian,
DiffContext  
)
virtualinherited

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

Most problems can use the reimplementation in FEMPhysics::mass_residual; few users will need to reimplement this themselves.

Reimplemented from libMesh::DifferentiablePhysics.

virtual bool libMesh::DifferentiablePhysics::mass_residual ( bool  request_jacobian,
DiffContext  
)
inlinevirtualinherited

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

Most problems can use the reimplementation in FEMSystem::mass_residual; few users will need to reimplement this themselves.

Reimplemented in libMesh::FEMPhysics.

Definition at line 229 of file diff_physics.h.

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

230  {
231  return request_jacobian;
232  }
void libMesh::FEMSystem::mesh_position_get ( )
inherited

Tells the FEMSystem to set the degree of freedom coefficients which should correspond to mesh nodal coordinates.

Definition at line 1026 of file fem_system.C.

References libMesh::DifferentiablePhysics::_mesh_sys, libMesh::DifferentiablePhysics::_mesh_x_var, libMesh::DifferentiablePhysics::_mesh_y_var, libMesh::DifferentiablePhysics::_mesh_z_var, libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::FEMSystem::build_context(), libMesh::FEMContext::elem_position_get(), libMesh::DiffContext::get_dof_indices(), libMesh::DiffContext::get_elem_solution(), libMesh::System::get_mesh(), libMesh::FEMSystem::init_context(), libMesh::invalid_uint, mesh, libMesh::FEMContext::pre_fe_reinit(), libMesh::System::solution, and libMesh::System::update().

1027 {
1028  // This function makes no sense unless we've already picked out some
1029  // variable(s) to reflect mesh position coordinates
1030  if (!_mesh_sys)
1031  libmesh_error();
1032 
1033  // We currently assume mesh variables are in our own system
1034  if (_mesh_sys != this)
1035  libmesh_not_implemented();
1036 
1037  // Loop over every active mesh element on this processor
1038  const MeshBase& mesh = this->get_mesh();
1039 
1042  const MeshBase::const_element_iterator end_el =
1044 
1045  AutoPtr<DiffContext> con = this->build_context();
1046  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
1047  this->init_context(_femcontext);
1048 
1049  // Get the solution's mesh variables from every element
1050  for ( ; el != end_el; ++el)
1051  {
1052  _femcontext.pre_fe_reinit(*this, *el);
1053 
1054  _femcontext.elem_position_get();
1055 
1057  this->solution->insert(_femcontext.get_elem_solution(_mesh_x_var),
1058  _femcontext.get_dof_indices(_mesh_x_var) );
1060  this->solution->insert(_femcontext.get_elem_solution(_mesh_y_var),
1061  _femcontext.get_dof_indices(_mesh_y_var));
1063  this->solution->insert(_femcontext.get_elem_solution(_mesh_z_var),
1064  _femcontext.get_dof_indices(_mesh_z_var));
1065  }
1066 
1067  this->solution->close();
1068 
1069  // And make sure the current_local_solution is up to date too
1070  this->System::update();
1071 }
void libMesh::FEMSystem::mesh_position_set ( )
inherited

Tells the FEMSystem to set the mesh nodal coordinates which should correspond to degree of freedom coefficients.

Definition at line 706 of file fem_system.C.

References libMesh::DifferentiablePhysics::_mesh_sys, libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), libMesh::FEMSystem::build_context(), libMesh::ParallelObject::comm(), libMesh::FEMContext::elem_fe_reinit(), libMesh::FEMContext::elem_position_set(), libMesh::FEMContext::get_elem(), libMesh::System::get_mesh(), libMesh::Elem::has_children(), libMesh::FEMSystem::init_context(), libMesh::libmesh_assert(), mesh, libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), libMesh::FEMContext::pre_fe_reinit(), and libMesh::Parallel::sync_dofobject_data_by_id().

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

707 {
708  // If we don't need to move the mesh, we're done
709  if (_mesh_sys != this)
710  return;
711 
712  MeshBase& mesh = this->get_mesh();
713 
714  AutoPtr<DiffContext> con = this->build_context();
715  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
716  this->init_context(_femcontext);
717 
718  // Move every mesh element we can
721  const MeshBase::const_element_iterator end_el =
723 
724  for ( ; el != end_el; ++el)
725  {
726  // We need the algebraic data
727  _femcontext.pre_fe_reinit(*this, *el);
728  // And when asserts are on, we also need the FE so
729  // we can assert that the mesh data is of the right type.
730 #ifndef NDEBUG
731  _femcontext.elem_fe_reinit();
732 #endif
733 
734  // This code won't handle moving subactive elements
735  libmesh_assert(!_femcontext.get_elem().has_children());
736 
737  _femcontext.elem_position_set(1.);
738  }
739 
740  // We've now got positions set on all local nodes (and some
741  // semilocal nodes); let's request positions for non-local nodes
742  // from their processors.
743 
744  SyncNodalPositions sync_object(mesh);
746  (this->comm(), mesh.nodes_begin(), mesh.nodes_end(), sync_object);
747 }
dof_id_type libMesh::System::n_active_dofs ( ) const
inlineinherited

Returns the number of active degrees of freedom for this System.

Definition at line 2063 of file system.h.

References libMesh::System::n_constrained_dofs(), and libMesh::System::n_dofs().

2064 {
2065  return this->n_dofs() - this->n_constrained_dofs();
2066 }
unsigned int libMesh::System::n_components ( ) const
inlineinherited
Returns
the total number of scalar components in the system's variables. This will equal n_vars() in the case of all scalar-valued variables.

Definition at line 1967 of file system.h.

References libMesh::System::_variables, libMesh::Variable::first_scalar_number(), and libMesh::Variable::n_components().

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

1968 {
1969  if (_variables.empty())
1970  return 0;
1971 
1972  const Variable& last = _variables.back();
1973  return last.first_scalar_number() + last.n_components();
1974 }
dof_id_type libMesh::System::n_constrained_dofs ( ) const
inherited
Returns
the total number of constrained degrees of freedom in the system.

Definition at line 147 of file system.C.

References libMesh::System::_dof_map.

Referenced by libMesh::System::get_info(), and libMesh::System::n_active_dofs().

148 {
149 #ifdef LIBMESH_ENABLE_CONSTRAINTS
150 
151  return _dof_map->n_constrained_dofs();
152 
153 #else
154 
155  return 0;
156 
157 #endif
158 }
dof_id_type libMesh::System::n_local_constrained_dofs ( ) const
inherited
Returns
the number of constrained degrees of freedom on this processor.

Definition at line 162 of file system.C.

References libMesh::System::_dof_map.

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

163 {
164 #ifdef LIBMESH_ENABLE_CONSTRAINTS
165 
166  return _dof_map->n_local_constrained_dofs();
167 
168 #else
169 
170  return 0;
171 
172 #endif
173 }
unsigned int libMesh::ImplicitSystem::n_matrices ( ) const
inlinevirtualinherited
Returns
the number of matrices handled by this system

Reimplemented from libMesh::System.

Definition at line 385 of file implicit_system.h.

References libMesh::ImplicitSystem::_matrices.

386 {
387  return libmesh_cast_int<unsigned int>(_matrices.size());
388 }
static unsigned int libMesh::ReferenceCounter::n_objects ( )
inlinestaticinherited

Prints the number of outstanding (created, but not yet destroyed) objects.

Definition at line 79 of file reference_counter.h.

References libMesh::ReferenceCounter::_n_objects.

80  { return _n_objects; }
processor_id_type libMesh::ParallelObject::n_processors ( ) const
inlineinherited
Returns
the number of processors in the group.

Definition at line 92 of file parallel_object.h.

References libMesh::ParallelObject::_communicator, and libMesh::Parallel::Communicator::size().

Referenced by libMesh::ParmetisPartitioner::_do_repartition(), libMesh::ParallelMesh::add_elem(), libMesh::ParallelMesh::add_node(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::ParmetisPartitioner::assign_partitioning(), libMesh::ParallelMesh::assign_unique_ids(), libMesh::AztecLinearSolver< T >::AztecLinearSolver(), libMesh::ParallelMesh::clear(), libMesh::Nemesis_IO_Helper::compute_border_node_ids(), libMesh::Nemesis_IO_Helper::construct_nemesis_filename(), libMesh::UnstructuredMesh::create_pid_mesh(), libMesh::DofMap::distribute_dofs(), libMesh::DofMap::distribute_local_dofs_node_major(), libMesh::DofMap::distribute_local_dofs_var_major(), libMesh::EnsightIO::EnsightIO(), libMesh::MeshBase::get_info(), libMesh::EquationSystems::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::ParmetisPartitioner::initialize(), libMesh::Nemesis_IO_Helper::initialize(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), libMesh::MeshTools::libmesh_assert_valid_procids< Elem >(), libMesh::MeshTools::libmesh_assert_valid_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_refinement_flags(), libMesh::MeshBase::n_active_elem_on_proc(), libMesh::MeshBase::n_elem_on_proc(), libMesh::MeshBase::n_nodes_on_proc(), libMesh::Partitioner::partition(), libMesh::MeshBase::partition(), libMesh::Partitioner::partition_unpartitioned_elements(), libMesh::PetscLinearSolver< T >::PetscLinearSolver(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::MeshTools::processor_bounding_box(), libMesh::System::project_vector(), libMesh::Nemesis_IO::read(), libMesh::CheckpointIO::read(), libMesh::UnstructuredMesh::read(), libMesh::System::read_parallel_data(), libMesh::System::read_SCALAR_dofs(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::System::read_serialized_vector(), libMesh::Partitioner::repartition(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::BoundaryInfo::sync(), libMesh::ParallelMesh::update_parallel_id_counts(), libMesh::CheckpointIO::write(), libMesh::GMVIO::write_binary(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::System::write_parallel_data(), libMesh::System::write_SCALAR_dofs(), libMesh::XdrIO::write_serialized_bcs(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), and libMesh::XdrIO::write_serialized_nodesets().

93  { return libmesh_cast_int<processor_id_type>(_communicator.size()); }
unsigned int libMesh::System::n_variable_groups ( ) const
inlineinherited
Returns
the number of VariableGroup variable groups in the system

Definition at line 1959 of file system.h.

References libMesh::System::_variable_groups.

Referenced by libMesh::System::add_variable(), libMesh::System::get_info(), and libMesh::System::init_data().

1960 {
1961  return libmesh_cast_int<unsigned int>(_variable_groups.size());
1962 }
unsigned int libMesh::System::n_vars ( ) const
inlineinherited
Returns
the number of variables in the system

Definition at line 1951 of file system.h.

References libMesh::System::_variables.

Referenced by libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::DiffContext::add_localized_vector(), libMesh::System::add_variable(), libMesh::System::add_variables(), libMesh::EquationSystems::build_discontinuous_solution_vector(), libMesh::EquationSystems::build_solution_vector(), libMesh::System::calculate_norm(), libMesh::DGFEMContext::DGFEMContext(), libMesh::DiffContext::DiffContext(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointResidualErrorEstimator::estimate_error(), libMesh::ErrorEstimator::estimate_errors(), libMesh::ExactSolution::ExactSolution(), libMesh::FEMContext::FEMContext(), libMesh::System::get_all_variable_numbers(), libMesh::EquationSystems::get_solution(), libMesh::System::init(), libMesh::FEMSystem::init_context(), libMesh::DGFEMContext::neighbor_side_fe_reinit(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::ProjectSolution::operator()(), libMesh::FEMContext::pre_fe_reinit(), libMesh::System::project_vector(), libMesh::System::re_update(), libMesh::System::read_legacy_data(), libMesh::System::read_parallel_data(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::System::read_serialized_vector(), libMesh::System::read_serialized_vectors(), libMesh::EquationSystems::reinit(), libMesh::System::reinit(), libMesh::HPCoarsenTest::select_refinement(), libMesh::SystemSubsetBySubdomain::set_var_nums(), libMesh::System::write_header(), libMesh::System::write_parallel_data(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::System::write_serialized_vector(), libMesh::System::write_serialized_vectors(), and libMesh::System::zero_variable().

1952 {
1953  return libmesh_cast_int<unsigned int>(_variables.size());
1954 }
unsigned int libMesh::System::n_vectors ( ) const
inlineinherited
Returns
the number of vectors (in addition to the solution) handled by this system This is the size of the _vectors map

Definition at line 2079 of file system.h.

References libMesh::System::_vectors.

Referenced by libMesh::ExplicitSystem::add_system_rhs(), libMesh::System::compare(), libMesh::System::get_info(), and libMesh::System::write_header().

2080 {
2081  return libmesh_cast_int<unsigned int>(_vectors.size());
2082 }
void libMesh::FEMSystem::numerical_elem_jacobian ( FEMContext context) const
inherited

Uses the results of multiple element_residual() calls to numerically differentiate the corresponding jacobian on an element.

Definition at line 933 of file fem_system.C.

References libMesh::TimeSolver::element_residual(), libMesh::FEMSystem::numerical_jacobian(), libMesh::START_LOG(), and libMesh::STOP_LOG().

934 {
935  START_LOG("numerical_elem_jacobian()", "FEMSystem");
937  STOP_LOG("numerical_elem_jacobian()", "FEMSystem");
938 }
void libMesh::FEMSystem::numerical_jacobian ( TimeSolverResPtr  res,
FEMContext context 
) const
inherited

Uses the results of multiple res calls to numerically differentiate the corresponding jacobian.

Definition at line 830 of file fem_system.C.

References libMesh::DifferentiablePhysics::_mesh_sys, libMesh::DifferentiablePhysics::_mesh_x_var, libMesh::DifferentiablePhysics::_mesh_y_var, libMesh::DifferentiablePhysics::_mesh_z_var, libMesh::DiffContext::get_dof_indices(), libMesh::FEMContext::get_elem(), libMesh::DiffContext::get_elem_jacobian(), libMesh::DiffContext::get_elem_residual(), libMesh::DiffContext::get_elem_solution(), libMesh::Elem::hmin(), libMesh::invalid_uint, libMesh::libmesh_real(), libMesh::FEMSystem::numerical_jacobian_h, libMesh::Elem::point(), libMesh::Real, and libMesh::DenseVector< T >::zero().

Referenced by libMesh::FEMSystem::numerical_elem_jacobian(), and libMesh::FEMSystem::numerical_side_jacobian().

832 {
833  // Logging is done by numerical_elem_jacobian
834  // or numerical_side_jacobian
835 
836  DenseVector<Number> original_residual(context.get_elem_residual());
837  DenseVector<Number> backwards_residual(context.get_elem_residual());
838  DenseMatrix<Number> numeric_jacobian(context.get_elem_jacobian());
839 #ifdef DEBUG
840  DenseMatrix<Number> old_jacobian(context.get_elem_jacobian());
841 #endif
842 
843  Real numerical_point_h = 0.;
844  if (_mesh_sys == this)
845  numerical_point_h = numerical_jacobian_h * context.get_elem().hmin();
846 
847  for (unsigned int j = 0; j != context.get_dof_indices().size(); ++j)
848  {
849  // Take the "minus" side of a central differenced first derivative
850  Number original_solution = context.get_elem_solution()(j);
851  context.get_elem_solution()(j) -= numerical_jacobian_h;
852 
853  // Make sure to catch any moving mesh terms
854  // FIXME - this could be less ugly
855  Real *coord = NULL;
856  if (_mesh_sys == this)
857  {
859  for (unsigned int k = 0;
860  k != context.get_dof_indices( _mesh_x_var ).size(); ++k)
861  if (context.get_dof_indices( _mesh_x_var )[k] ==
862  context.get_dof_indices()[j])
863  coord = &(context.get_elem().point(k)(0));
865  for (unsigned int k = 0;
866  k != context.get_dof_indices( _mesh_y_var ).size(); ++k)
867  if (context.get_dof_indices( _mesh_y_var )[k] ==
868  context.get_dof_indices()[j])
869  coord = &(context.get_elem().point(k)(1));
871  for (unsigned int k = 0;
872  k != context.get_dof_indices( _mesh_z_var ).size(); ++k)
873  if (context.get_dof_indices( _mesh_z_var )[k] ==
874  context.get_dof_indices()[j])
875  coord = &(context.get_elem().point(k)(2));
876  }
877  if (coord)
878  {
879  // We have enough information to scale the perturbations
880  // here appropriately
881  context.get_elem_solution()(j) = original_solution - numerical_point_h;
882  *coord = libmesh_real(context.get_elem_solution()(j));
883  }
884 
885  context.get_elem_residual().zero();
886  ((*time_solver).*(res))(false, context);
887 #ifdef DEBUG
888  libmesh_assert_equal_to (old_jacobian, context.get_elem_jacobian());
889 #endif
890  backwards_residual = context.get_elem_residual();
891 
892  // Take the "plus" side of a central differenced first derivative
893  context.get_elem_solution()(j) = original_solution + numerical_jacobian_h;
894  if (coord)
895  {
896  context.get_elem_solution()(j) = original_solution + numerical_point_h;
897  *coord = libmesh_real(context.get_elem_solution()(j));
898  }
899  context.get_elem_residual().zero();
900  ((*time_solver).*(res))(false, context);
901 #ifdef DEBUG
902  libmesh_assert_equal_to (old_jacobian, context.get_elem_jacobian());
903 #endif
904 
905  context.get_elem_solution()(j) = original_solution;
906  if (coord)
907  {
908  *coord = libmesh_real(context.get_elem_solution()(j));
909  for (unsigned int i = 0; i != context.get_dof_indices().size(); ++i)
910  {
911  numeric_jacobian(i,j) =
912  (context.get_elem_residual()(i) - backwards_residual(i)) /
913  2. / numerical_point_h;
914  }
915  }
916  else
917  {
918  for (unsigned int i = 0; i != context.get_dof_indices().size(); ++i)
919  {
920  numeric_jacobian(i,j) =
921  (context.get_elem_residual()(i) - backwards_residual(i)) /
923  }
924  }
925  }
926 
927  context.get_elem_residual() = original_residual;
928  context.get_elem_jacobian() = numeric_jacobian;
929 }
void libMesh::FEMSystem::numerical_side_jacobian ( FEMContext context) const
inherited

Uses the results of multiple side_residual() calls to numerically differentiate the corresponding jacobian on an element's side.

Definition at line 942 of file fem_system.C.

References libMesh::FEMSystem::numerical_jacobian(), libMesh::TimeSolver::side_residual(), libMesh::START_LOG(), and libMesh::STOP_LOG().

943 {
944  START_LOG("numerical_side_jacobian()", "FEMSystem");
946  STOP_LOG("numerical_side_jacobian()", "FEMSystem");
947 }
virtual void libMesh::DifferentiableQoI::parallel_op ( const Parallel::Communicator communicator,
std::vector< Number > &  sys_qoi,
std::vector< Number > &  local_qoi,
const QoISet qoi_indices 
)
virtualinherited

Method to populate system qoi data structure with process-local qoi. By default, simply sums process qois into system qoi.

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

Gradient libMesh::System::point_gradient ( unsigned int  var,
const Point p,
const bool  insist_on_success = true 
) const
inherited

Returns the gradient of the solution variable var at the physical point p in the mesh, similarly to point_value.

Definition at line 2063 of file system.C.

References libMesh::Parallel::Communicator::broadcast(), libMesh::ParallelObject::comm(), libMesh::PointLocatorBase::enable_out_of_mesh_mode(), libMesh::System::get_mesh(), libMesh::libmesh_assert(), mesh, libMesh::Parallel::Communicator::min(), libMesh::ParallelObject::n_processors(), libMesh::ParallelObject::processor_id(), libMesh::DofObject::processor_id(), libMesh::MeshBase::sub_point_locator(), and libMesh::Parallel::Communicator::verify().

2064 {
2065  // This function must be called on every processor; there's no
2066  // telling where in the partition p falls.
2067  parallel_object_only();
2068 
2069  // And every processor had better agree about which point we're
2070  // looking for
2071 #ifndef NDEBUG
2072  this->comm().verify(p);
2073 #endif // NDEBUG
2074 
2075  // Get a reference to the mesh object associated with the system object that calls this function
2076  const MeshBase &mesh = this->get_mesh();
2077 
2078  // Use an existing PointLocator or create a new one
2079  AutoPtr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
2080  PointLocatorBase& locator = *locator_ptr;
2081 
2082  if (!insist_on_success)
2083  locator.enable_out_of_mesh_mode();
2084 
2085  // Get a pointer to the element that contains P
2086  const Elem *e = locator(p);
2087 
2088  Gradient grad_u;
2089 
2090  if (e && e->processor_id() == this->processor_id())
2091  grad_u = point_gradient(var, p, *e);
2092 
2093  // If I have an element containing p, then let's let everyone know
2094  processor_id_type lowest_owner =
2095  (e && (e->processor_id() == this->processor_id())) ?
2096  this->processor_id() : this->n_processors();
2097  this->comm().min(lowest_owner);
2098 
2099  // Everybody should get their value from a processor that was able
2100  // to compute it.
2101  // If nobody admits owning the point, we may have a problem.
2102  if (lowest_owner != this->n_processors())
2103  this->comm().broadcast(grad_u, lowest_owner);
2104  else
2105  libmesh_assert(!insist_on_success);
2106 
2107  return grad_u;
2108 }
Gradient libMesh::System::point_gradient ( unsigned int  var,
const Point p,
const Elem e 
) const
inherited

Returns the gradient of the solution variable var at the physical point p in local Elem e in the mesh, similarly to point_value.

Definition at line 2111 of file system.C.

References libMesh::TypeVector< T >::add_scaled(), libMesh::FEGenericBase< T >::build(), libMesh::Elem::contains_point(), libMesh::System::current_solution(), libMesh::Elem::dim(), libMesh::dof_map, libMesh::System::get_dof_map(), libMesh::FEInterface::inverse_map(), libMesh::libmesh_assert(), libMesh::ParallelObject::processor_id(), and libMesh::DofObject::processor_id().

2112 {
2113  libmesh_assert_equal_to (e.processor_id(), this->processor_id());
2114 
2115  // Ensuring that the given point is really in the element is an
2116  // expensive assert, but as long as debugging is turned on we might
2117  // as well try to catch a particularly nasty potential error
2118  libmesh_assert (e.contains_point(p));
2119 
2120  // Get the dof map to get the proper indices for our computation
2121  const DofMap& dof_map = this->get_dof_map();
2122 
2123  // Need dof_indices for phi[i][j]
2124  std::vector<dof_id_type> dof_indices;
2125 
2126  // Fill in the dof_indices for our element
2127  dof_map.dof_indices (&e, dof_indices, var);
2128 
2129  // Get the no of dofs assciated with this point
2130  const unsigned int num_dofs = libmesh_cast_int<unsigned int>
2131  (dof_indices.size());
2132 
2133  FEType fe_type = dof_map.variable_type(var);
2134 
2135  // Build a FE again so we can calculate u(p)
2136  AutoPtr<FEBase> fe (FEBase::build(e.dim(), fe_type));
2137 
2138  // Map the physical co-ordinates to the master co-ordinates using the inverse_map from fe_interface.h
2139  // Build a vector of point co-ordinates to send to reinit
2140  std::vector<Point> coor(1, FEInterface::inverse_map(e.dim(), fe_type, &e, p));
2141 
2142  // Get the values of the shape function derivatives
2143  const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
2144 
2145  // Reinitialize the element and compute the shape function values at coor
2146  fe->reinit (&e, &coor);
2147 
2148  // Get ready to accumulate a gradient
2149  Gradient grad_u;
2150 
2151  for (unsigned int l=0; l<num_dofs; l++)
2152  {
2153  grad_u.add_scaled (dphi[l][0], this->current_solution (dof_indices[l]));
2154  }
2155 
2156  return grad_u;
2157 }
Tensor libMesh::System::point_hessian ( unsigned int  var,
const Point p,
const bool  insist_on_success = true 
) const
inherited

Returns the second derivative tensor of the solution variable var at the physical point p in the mesh, similarly to point_value.

Definition at line 2162 of file system.C.

References libMesh::Parallel::Communicator::broadcast(), libMesh::ParallelObject::comm(), libMesh::PointLocatorBase::enable_out_of_mesh_mode(), libMesh::System::get_mesh(), libMesh::libmesh_assert(), mesh, libMesh::Parallel::Communicator::min(), libMesh::ParallelObject::n_processors(), libMesh::ParallelObject::processor_id(), libMesh::DofObject::processor_id(), libMesh::MeshBase::sub_point_locator(), and libMesh::Parallel::Communicator::verify().

2163 {
2164  // This function must be called on every processor; there's no
2165  // telling where in the partition p falls.
2166  parallel_object_only();
2167 
2168  // And every processor had better agree about which point we're
2169  // looking for
2170 #ifndef NDEBUG
2171  this->comm().verify(p);
2172 #endif // NDEBUG
2173 
2174  // Get a reference to the mesh object associated with the system object that calls this function
2175  const MeshBase &mesh = this->get_mesh();
2176 
2177  // Use an existing PointLocator or create a new one
2178  AutoPtr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
2179  PointLocatorBase& locator = *locator_ptr;
2180 
2181  if (!insist_on_success)
2182  locator.enable_out_of_mesh_mode();
2183 
2184  // Get a pointer to the element that contains P
2185  const Elem *e = locator(p);
2186 
2187  Tensor hess_u;
2188 
2189  if (e && e->processor_id() == this->processor_id())
2190  hess_u = point_hessian(var, p, *e);
2191 
2192  // If I have an element containing p, then let's let everyone know
2193  processor_id_type lowest_owner =
2194  (e && (e->processor_id() == this->processor_id())) ?
2195  this->processor_id() : this->n_processors();
2196  this->comm().min(lowest_owner);
2197 
2198  // Everybody should get their value from a processor that was able
2199  // to compute it.
2200  // If nobody admits owning the point, we may have a problem.
2201  if (lowest_owner != this->n_processors())
2202  this->comm().broadcast(hess_u, lowest_owner);
2203  else
2204  libmesh_assert(!insist_on_success);
2205 
2206  return hess_u;
2207 }
Tensor libMesh::System::point_hessian ( unsigned int  var,
const Point p,
const Elem e 
) const
inherited

Returns the second derivative tensor of the solution variable var at the physical point p in local Elem e in the mesh, similarly to point_value.

Definition at line 2209 of file system.C.

References libMesh::TypeTensor< T >::add_scaled(), libMesh::FEGenericBase< T >::build(), libMesh::Elem::contains_point(), libMesh::System::current_solution(), libMesh::Elem::dim(), libMesh::dof_map, libMesh::System::get_dof_map(), libMesh::FEInterface::inverse_map(), libMesh::libmesh_assert(), libMesh::ParallelObject::processor_id(), and libMesh::DofObject::processor_id().

2210 {
2211  libmesh_assert_equal_to (e.processor_id(), this->processor_id());
2212 
2213  // Ensuring that the given point is really in the element is an
2214  // expensive assert, but as long as debugging is turned on we might
2215  // as well try to catch a particularly nasty potential error
2216  libmesh_assert (e.contains_point(p));
2217 
2218  // Get the dof map to get the proper indices for our computation
2219  const DofMap& dof_map = this->get_dof_map();
2220 
2221  // Need dof_indices for phi[i][j]
2222  std::vector<dof_id_type> dof_indices;
2223 
2224  // Fill in the dof_indices for our element
2225  dof_map.dof_indices (&e, dof_indices, var);
2226 
2227  // Get the no of dofs assciated with this point
2228  const unsigned int num_dofs = libmesh_cast_int<unsigned int>
2229  (dof_indices.size());
2230 
2231  FEType fe_type = dof_map.variable_type(var);
2232 
2233  // Build a FE again so we can calculate u(p)
2234  AutoPtr<FEBase> fe (FEBase::build(e.dim(), fe_type));
2235 
2236  // Map the physical co-ordinates to the master co-ordinates using the inverse_map from fe_interface.h
2237  // Build a vector of point co-ordinates to send to reinit
2238  std::vector<Point> coor(1, FEInterface::inverse_map(e.dim(), fe_type, &e, p));
2239 
2240  // Get the values of the shape function derivatives
2241  const std::vector<std::vector<RealTensor> >& d2phi = fe->get_d2phi();
2242 
2243  // Reinitialize the element and compute the shape function values at coor
2244  fe->reinit (&e, &coor);
2245 
2246  // Get ready to accumulate a hessian
2247  Tensor hess_u;
2248 
2249  for (unsigned int l=0; l<num_dofs; l++)
2250  {
2251  hess_u.add_scaled (d2phi[l][0], this->current_solution (dof_indices[l]));
2252  }
2253 
2254  return hess_u;
2255 }
Number libMesh::System::point_value ( unsigned int  var,
const Point p,
const bool  insist_on_success = true 
) const
inherited

Returns the value of the solution variable var at the physical point p in the mesh, without knowing a priori which element contains p.

Note that this function uses MeshBase::sub_point_locator(); users may or may not want to call MeshBase::clear_point_locator() afterward. Also, point_locator() is expensive (N log N for initial construction, log N for evaluations). Avoid using this function in any context where you are already looping over elements.

Because the element containing p may lie on any processor, this function is parallel-only.

By default this method expects the point to reside inside the domain and will abort if no element can be found which contains . The optional parameter insist_on_success can be set to false to allow the method to return 0 when the point is not located.

Definition at line 1966 of file system.C.

References libMesh::Parallel::Communicator::broadcast(), libMesh::ParallelObject::comm(), libMesh::PointLocatorBase::enable_out_of_mesh_mode(), libMesh::System::get_mesh(), libMesh::libmesh_assert(), mesh, libMesh::Parallel::Communicator::min(), libMesh::ParallelObject::n_processors(), libMesh::ParallelObject::processor_id(), libMesh::DofObject::processor_id(), libMesh::MeshBase::sub_point_locator(), and libMesh::Parallel::Communicator::verify().

1967 {
1968  // This function must be called on every processor; there's no
1969  // telling where in the partition p falls.
1970  parallel_object_only();
1971 
1972  // And every processor had better agree about which point we're
1973  // looking for
1974 #ifndef NDEBUG
1975  this->comm().verify(p);
1976 #endif // NDEBUG
1977 
1978  // Get a reference to the mesh object associated with the system object that calls this function
1979  const MeshBase &mesh = this->get_mesh();
1980 
1981  // Use an existing PointLocator or create a new one
1982  AutoPtr<PointLocatorBase> locator_ptr = mesh.sub_point_locator();
1983  PointLocatorBase& locator = *locator_ptr;
1984 
1985  if (!insist_on_success)
1986  locator.enable_out_of_mesh_mode();
1987 
1988  // Get a pointer to the element that contains P
1989  const Elem *e = locator(p);
1990 
1991  Number u = 0;
1992 
1993  if (e && e->processor_id() == this->processor_id())
1994  u = point_value(var, p, *e);
1995 
1996  // If I have an element containing p, then let's let everyone know
1997  processor_id_type lowest_owner =
1998  (e && (e->processor_id() == this->processor_id())) ?
1999  this->processor_id() : this->n_processors();
2000  this->comm().min(lowest_owner);
2001 
2002  // Everybody should get their value from a processor that was able
2003  // to compute it.
2004  // If nobody admits owning the point, we have a problem.
2005  if (lowest_owner != this->n_processors())
2006  this->comm().broadcast(u, lowest_owner);
2007  else
2008  libmesh_assert(!insist_on_success);
2009 
2010  return u;
2011 }
Number libMesh::System::point_value ( unsigned int  var,
const Point p,
const Elem e 
) const
inherited

Returns the value of the solution variable var at the physical point p contained in local Elem e

This version of point_value can be run in serial, but assumes e is in the local mesh partition.

Definition at line 2013 of file system.C.

References libMesh::FEGenericBase< T >::build(), libMesh::Elem::contains_point(), libMesh::System::current_solution(), libMesh::Elem::dim(), libMesh::dof_map, libMesh::System::get_dof_map(), libMesh::FEInterface::inverse_map(), libMesh::libmesh_assert(), libMesh::ParallelObject::processor_id(), and libMesh::DofObject::processor_id().

2014 {
2015  libmesh_assert_equal_to (e.processor_id(), this->processor_id());
2016 
2017  // Ensuring that the given point is really in the element is an
2018  // expensive assert, but as long as debugging is turned on we might
2019  // as well try to catch a particularly nasty potential error
2020  libmesh_assert (e.contains_point(p));
2021 
2022  // Get the dof map to get the proper indices for our computation
2023  const DofMap&