ContinuationSystem Class Reference

#include <continuation_system.h>

Inheritance diagram for ContinuationSystem:

List of all members.

Public Types

enum  Predictor { Euler, AB2, Invalid_Predictor }
typedef ContinuationSystem sys_type
typedef FEMSystem Parent
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)
virtual void time_evolving (unsigned int var)
virtual void mesh_x_position (unsigned int sysnum, unsigned int var)
virtual void mesh_y_position (unsigned int sysnum, unsigned int var)
virtual void mesh_z_position (unsigned int sysnum, unsigned int var)
void mesh_position_get ()
void mesh_position_set ()
virtual bool eulerian_residual (bool request_jacobian, DiffContext &context)
virtual bool mass_residual (bool request_jacobian, DiffContext &context)
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())
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 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 element_postprocess (DiffContext &)
virtual void side_postprocess (DiffContext &)
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 bool side_mass_residual (bool request_jacobian, DiffContext &)
sys_typesystem ()
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
adjoint_solve (const QoISet &qoi_indices=QoISet())
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_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)
unsigned int n_matrices () const
void init ()
virtual void update ()
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 (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), Parameters &parameters) 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), 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 ()
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)
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)
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_dofs () const
unsigned int n_active_dofs () const
unsigned int n_constrained_dofs () const
unsigned int 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)
const Variablevariable (unsigned int var) 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
const FETypevariable_type (const unsigned int i) const
const FETypevariable_type (const std::string &var) const
Real calculate_norm (NumericVector< Number > &v, unsigned int var=0, FEMNormType norm_type=L2) const
Real calculate_norm (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)
void read_serialized_data (Xdr &io, const bool read_additional_data=true)
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
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_assemble_function (void fptr(EquationSystems &es, const std::string &name))
void attach_constraint_function (void fptr(EquationSystems &es, const std::string &name))
void attach_QOI_function (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))
void attach_QOI_derivative (void fptr(EquationSystems &es, const std::string &name, const QoISet &qoi_indices))
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 unsigned int global_dof_number) const
void local_dof_indices (const unsigned int var, std::set< unsigned int > &var_indices) const
void zero_variable (NumericVector< Number > &v, unsigned int var_num) const

Static Public Member Functions

static std::string get_info ()
static void print_info ()
static unsigned int n_objects ()

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
int extra_quadrature_order
Real numerical_jacobian_h
Real verify_analytic_jacobians
bool compute_internal_sides
bool postprocess_sides
bool assemble_qoi_sides
AutoPtr< TimeSolvertime_solver
Real time
Real deltat
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
bool use_fixed_solution
SparseMatrix< Number > * matrix
NumericVector< Number > * rhs
bool assemble_before_solve
AutoPtr< NumericVector< Number > > solution
AutoPtr< NumericVector< Number > > current_local_solution
std::vector< Numberqoi

Protected Types

enum  RHS_Mode { Residual, G_Lambda }
typedef bool(TimeSolver::* TimeSolverResPtr )(bool, DiffContext &)
typedef std::map< std::string,
std::pair< unsigned int,
unsigned int > > 
Counts

Protected Member Functions

virtual void init_data ()
void numerical_jacobian (TimeSolverResPtr res, FEMContext &context)
void numerical_elem_jacobian (FEMContext &context)
void numerical_side_jacobian (FEMContext &context)
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
unsigned int _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

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 56 of file continuation_system.h.


Member Typedef Documentation

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

Definition at line 256 of file implicit_system.h.

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

Definition at line 405 of file system.h.

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > ReferenceCounter::Counts [protected, inherited]

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

Definition at line 105 of file reference_counter.h.

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

Matrix iterator typedefs.

Definition at line 255 of file implicit_system.h.

The type of the parent.

Reimplemented from FEMSystem.

Definition at line 80 of file continuation_system.h.

The type of system.

Reimplemented from FEMSystem.

Definition at line 75 of file continuation_system.h.

typedef bool(TimeSolver::* FEMSystem::TimeSolverResPtr)(bool, DiffContext &) [protected, inherited]

Syntax sugar to make numerical_jacobian() declaration easier.

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

Vector iterator typedefs.

Definition at line 404 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 222 of file continuation_system.h.

00222                  {
00226     Euler,
00227 
00231     AB2,
00232 
00236     Invalid_Predictor
00237   };

enum ContinuationSystem::RHS_Mode [protected]

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

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

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

Enumerator:
Residual 
G_Lambda 

Definition at line 287 of file continuation_system.h.

00287                 {Residual,
00288                  G_Lambda};


Constructor & Destructor Documentation

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

Constructor. Optionally initializes required data structures.

Definition at line 27 of file continuation_system.C.

00030   : Parent(es, name, number),
00031     continuation_parameter(NULL),
00032     quiet(true),
00033     continuation_parameter_tolerance(1.e-6),
00034     solution_tolerance(1.e-6),
00035     initial_newton_tolerance(0.01),
00036     old_continuation_parameter(0.),
00037     min_continuation_parameter(0.),
00038     max_continuation_parameter(0.),
00039     Theta(1.),
00040     Theta_LOCA(1.),
00041     //tau(1.),
00042     n_backtrack_steps(5),
00043     n_arclength_reductions(5),
00044     ds_min(1.e-8),
00045     predictor(Euler),
00046     newton_stepgrowth_aggressiveness(1.),
00047     newton_progress_check(true),
00048     rhs_mode(Residual),
00049     linear_solver(LinearSolver<Number>::build()),
00050     tangent_initialized(false),
00051     newton_solver(NULL),
00052     dlambda_ds(0.707),
00053     ds(0.1),
00054     ds_current(0.1),
00055     previous_dlambda_ds(0.),
00056     previous_ds(0.),
00057     newton_step(0)
00058 {
00059   // Warn about using untested code
00060   libmesh_experimental();
00061 }

ContinuationSystem::~ContinuationSystem (  )  [virtual]

Destructor.

Definition at line 66 of file continuation_system.C.

References clear().

00067 {
00068   this->clear();
00069 }


Member Function Documentation

void System::activate (  )  [inline, inherited]

Activates the system. Only active systems are solved.

Definition at line 1374 of file system.h.

References System::_active.

01375 {
01376   _active = true;
01377 }

bool System::active (  )  const [inline, inherited]

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

Definition at line 1366 of file system.h.

References System::_active.

01367 {
01368   return _active;  
01369 }

NumericVector< Number > & 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 824 of file system.C.

References System::add_vector().

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

00825 {
00826   OStringStream adjoint_rhs_name;
00827   adjoint_rhs_name << "adjoint_rhs" << i;
00828 
00829   return this->add_vector(adjoint_rhs_name.str());
00830 }

NumericVector< Number > & 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 764 of file system.C.

References System::add_vector().

Referenced by ImplicitSystem::adjoint_solve().

00765 {
00766   OStringStream adjoint_name;
00767   adjoint_name << "adjoint_solution" << i;
00768 
00769   return this->add_vector(adjoint_name.str());
00770 }

SparseMatrix< Number > & 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 175 of file implicit_system.C.

References ImplicitSystem::_can_add_matrices, ImplicitSystem::_matrices, and ImplicitSystem::have_matrix().

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

00176 {
00177   // only add matrices before initializing...
00178   if (!_can_add_matrices)
00179     {
00180       std::cerr << "ERROR: Too late.  Cannot add matrices to the system after initialization"
00181                 << std::endl
00182                 << " any more.  You should have done this earlier."
00183                 << std::endl;
00184       libmesh_error();
00185     }
00186 
00187   // Return the matrix if it is already there.
00188   if (this->have_matrix(mat_name))
00189     return *(_matrices[mat_name]);
00190 
00191   // Otherwise build the matrix and return it.
00192   SparseMatrix<Number>* buf = SparseMatrix<Number>::build().release();
00193   _matrices.insert (std::make_pair (mat_name, buf));
00194 
00195   return *buf;
00196 }

NumericVector< Number > & 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 854 of file system.C.

References System::add_vector().

Referenced by ImplicitSystem::assemble_residual_derivatives().

00855 {
00856   OStringStream sensitivity_rhs_name;
00857   sensitivity_rhs_name << "sensitivity_rhs" << i;
00858 
00859   return this->add_vector(sensitivity_rhs_name.str());
00860 }

NumericVector< Number > & 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 714 of file system.C.

References System::add_vector().

00715 {
00716   OStringStream sensitivity_name;
00717   sensitivity_name << "sensitivity_solution" << i;
00718 
00719   return this->add_vector(sensitivity_name.str());
00720 }

unsigned int 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 923 of file system.C.

References System::add_variable().

00927 {
00928   return this->add_variable(var, 
00929                             FEType(order, family), 
00930                             active_subdomains);
00931 }

unsigned int 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 884 of file system.C.

References System::_dof_map, System::_variable_numbers, System::_variables, System::n_vars(), System::number(), System::variable_name(), and System::variable_type().

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

00887 {  
00888   // Make sure the variable isn't there already
00889   // or if it is, that it's the type we want
00890   for (unsigned int v=0; v<this->n_vars(); v++)
00891     if (this->variable_name(v) == var)
00892       {
00893         if (this->variable_type(v) == type)
00894           return _variables[v].number();
00895 
00896         std::cerr << "ERROR: incompatible variable "
00897                   << var
00898                   << " has already been added for this system!"
00899                   << std::endl;
00900         libmesh_error();
00901       }
00902 
00903   const unsigned int curr_n_vars = this->n_vars();                                              
00904 
00905   // Add the variable to the list
00906   _variables.push_back((active_subdomains == NULL) ?
00907                        Variable(var, curr_n_vars, type) :
00908                        Variable(var, curr_n_vars, type, *active_subdomains));
00909 
00910   libmesh_assert ((curr_n_vars+1) == this->n_vars());
00911 
00912   _variable_numbers[var] = curr_n_vars;
00913 
00914   // Add the variable to the _dof_map
00915   _dof_map->add_variable (_variables.back());
00916 
00917   // Return the number of the new variable
00918   return curr_n_vars;
00919 }

NumericVector< Number > & System::add_vector ( const std::string &  vec_name,
const bool  projections = true 
) [inherited]

Adds the additional vector vec_name to this system. Only allowed prior to init(). 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 549 of file system.C.

References System::_can_add_vectors, System::_vector_projections, System::_vectors, System::have_vector(), NumericVector< T >::init(), System::n_dofs(), System::n_local_dofs(), and libMeshEnums::PARALLEL.

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

00551 {
00552   // Return the vector if it is already there.
00553   if (this->have_vector(vec_name))
00554     return *(_vectors[vec_name]);
00555 
00556   // Otherwise build the vector
00557   NumericVector<Number>* buf = NumericVector<Number>::build().release();
00558   _vectors.insert (std::make_pair (vec_name, buf));
00559   _vector_projections.insert (std::make_pair (vec_name, projections));
00560 
00561   // Initialize it if necessary
00562   if (!_can_add_vectors)
00563     buf->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
00564 
00565   return *buf;
00566 }

NumericVector< Number > & 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 794 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::weighted_sensitivity_adjoint_solve().

00795 {
00796   OStringStream adjoint_name;
00797   adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
00798 
00799   return this->get_vector(adjoint_name.str());
00800 }

NumericVector< Number > & 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 743 of file system.C.

References System::add_vector().

Referenced by ImplicitSystem::weighted_sensitivity_solve().

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

void ImplicitSystem::adjoint_qoi_parameter_sensitivity ( const QoISet qoi_indices,
const ParameterVector parameters,
SensitivityData sensitivities 
) [virtual, inherited]

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 System.

Definition at line 657 of file implicit_system.C.

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

00660 {
00661   const unsigned int Np = parameters.size();
00662   const unsigned int Nq = qoi.size();
00663 
00664   // We currently get partial derivatives via central differencing
00665   const Real delta_p = TOLERANCE;
00666 
00667   // An introduction to the problem:
00668   //
00669   // Residual R(u(p),p) = 0
00670   // partial R / partial u = J = system matrix
00671   //
00672   // This implies that:
00673   // d/dp(R) = 0
00674   // (partial R / partial p) + 
00675   // (partial R / partial u) * (partial u / partial p) = 0
00676 
00677   // We first do an adjoint solve:
00678   // J^T * z = (partial q / partial u)
00679 
00680   this->adjoint_solve(qoi_indices);
00681 
00682   // Get ready to fill in senstivities:
00683   sensitivities.allocate_data(qoi_indices, *this, parameters);
00684 
00685   // We use the identities:
00686   // dq/dp = (partial q / partial p) + (partial q / partial u) *
00687   //         (partial u / partial p)
00688   // dq/dp = (partial q / partial p) + (J^T * z) *
00689   //         (partial u / partial p)
00690   // dq/dp = (partial q / partial p) + z * J *
00691   //         (partial u / partial p)
00692  
00693   // Leading to our final formula:
00694   // dq/dp = (partial q / partial p) - z * (partial R / partial p)
00695 
00696   for (unsigned int j=0; j != Np; ++j)
00697     {
00698       // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
00699       // (partial R / partial p) ~= (rhs(p+dp) - rhs(p-dp))/(2*dp)
00700 
00701       Number old_parameter = *parameters[j];
00702       // Number old_qoi = this->qoi;
00703 
00704       *parameters[j] = old_parameter - delta_p;
00705       this->assemble_qoi(qoi_indices);
00706       std::vector<Number> qoi_minus = this->qoi;
00707 
00708       this->assembly(true, false);
00709       this->rhs->close();
00710       AutoPtr<NumericVector<Number> > partialR_partialp = this->rhs->clone();
00711       *partialR_partialp *= -1;
00712 
00713       *parameters[j] = old_parameter + delta_p;
00714       this->assemble_qoi(qoi_indices);
00715       std::vector<Number>& qoi_plus = this->qoi;
00716 
00717       std::vector<Number> partialq_partialp(Nq, 0);
00718       for (unsigned int i=0; i != Nq; ++i)
00719         if (qoi_indices.has_index(i))
00720           partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
00721 
00722       this->assembly(true, false);
00723       this->rhs->close();
00724       *partialR_partialp += *this->rhs;
00725       *partialR_partialp /= (2.*delta_p);
00726 
00727       // Don't leave the parameter changed
00728       *parameters[j] = old_parameter;
00729 
00730       for (unsigned int i=0; i != Nq; ++i)
00731         if (qoi_indices.has_index(i))
00732           sensitivities[i][j] = partialq_partialp[i] -
00733             partialR_partialp->dot(this->get_adjoint_solution(i));
00734     }
00735 
00736   // All parameters have been reset.
00737   // We didn't cache the original rhs or matrix for memory reasons,
00738   // but we can restore them to a state consistent solution -
00739   // principle of least surprise.
00740   this->assembly(true, true);
00741   this->rhs->close();
00742   this->matrix->close();
00743   this->assemble_qoi(qoi_indices);
00744 }

std::pair< unsigned int, Real > ImplicitSystem::adjoint_solve ( const QoISet qoi_indices = QoISet()  )  [virtual, inherited]

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

Leave qoi_indices empty to solve all adjoint problems.

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

Reimplemented from System.

Definition at line 342 of file implicit_system.C.

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

00343 {
00344   // Log how long the linear solve takes.
00345   START_LOG("adjoint_solve()", "ImplicitSystem");
00346 
00347   // The forward system should now already be solved.
00348   // Now assemble it's adjoint.
00349 
00350   if (this->assemble_before_solve)
00351     {
00352       // Build the Jacobian
00353       this->assembly(false, true);
00354       this->matrix->close();
00355 
00356       // Take the discrete adjoint
00357       matrix->get_transpose(*matrix);
00358 
00359       // Reset and build the RHS from the QOI derivative
00360       this->assemble_qoi_derivative(qoi_indices);
00361     }
00362 
00363   // The adjoint problem is linear
00364   LinearSolver<Number> *linear_solver = this->get_linear_solver();
00365 
00366   // Our iteration counts and residuals will be sums of the individual
00367   // results
00368   std::pair<unsigned int, Real> solver_params =
00369     this->get_linear_solve_parameters();
00370   std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);
00371 
00372   for (unsigned int i=0; i != this->qoi.size(); ++i)
00373     if (qoi_indices.has_index(i))
00374       {
00375         const std::pair<unsigned int, Real> rval =
00376           linear_solver->solve (*matrix, this->add_adjoint_solution(i),
00377                                  this->get_adjoint_rhs(i),
00378                                  solver_params.second,
00379                                  solver_params.first);
00380 
00381         totalrval.first  += rval.first;
00382         totalrval.second += rval.second;
00383       }
00384 
00385   this->release_linear_solver(linear_solver);
00386 
00387   // The linear solver may not have fit our constraints exactly
00388 #ifdef LIBMESH_ENABLE_AMR
00389   for (unsigned int i=0; i != this->qoi.size(); ++i)
00390     if (qoi_indices.has_index(i))
00391       this->get_dof_map().enforce_constraints_exactly
00392         (*this, &this->get_adjoint_solution(i));
00393 #endif
00394 
00395   // Stop logging the nonlinear solve
00396   STOP_LOG("adjoint_solve()", "ImplicitSystem");
00397 
00398   return totalrval;
00399 }

void ContinuationSystem::advance_arcstep (  ) 

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

Definition at line 945 of file continuation_system.C.

References solve_tangent(), and update_solution().

00946 {
00947   // Solve for the updated tangent du1/ds, d(lambda1)/ds
00948   solve_tangent();
00949 
00950   // Advance the solution and the parameter to the next value.
00951   update_solution();
00952 }

void ContinuationSystem::apply_predictor (  )  [private]

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

Definition at line 1384 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, and System::solution.

Referenced by continuation_solve(), and update_solution().

01385 {
01386   if (predictor == Euler)
01387     {
01388       // 1.) Euler Predictor
01389       // Predict next the solution
01390       solution->add(ds_current, *du_ds);
01391       solution->close();
01392   
01393       // Predict next parameter value
01394       *continuation_parameter += ds_current*dlambda_ds;
01395     }
01396 
01397 
01398   else if (predictor == AB2)
01399     {
01400       // 2.) 2nd-order explicit AB predictor
01401       libmesh_assert(previous_ds != 0.0);
01402       const Real stepratio = ds_current/previous_ds;
01403       
01404       // Build up next solution value.
01405       solution->add( 0.5*ds_current*(2.+stepratio), *du_ds);
01406       solution->add(-0.5*ds_current*stepratio     , *previous_du_ds);
01407       solution->close();
01408       
01409       // Next parameter value
01410       *continuation_parameter +=
01411         0.5*ds_current*((2.+stepratio)*dlambda_ds -
01412                         stepratio*previous_dlambda_ds);
01413     }
01414 
01415   else
01416     {
01417       // Unknown predictor
01418       libmesh_error();
01419     }
01420   
01421 }

void DifferentiableSystem::assemble (  )  [virtual, inherited]

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

Reimplemented from ImplicitSystem.

Definition at line 102 of file diff_system.C.

References DifferentiableSystem::assembly().

00103 {
00104   this->assembly(true, true);
00105 }

void FEMSystem::assemble_qoi ( const QoISet indices = QoISet()  )  [virtual, inherited]

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 for quantities of interest that are not expressible as a sum of element qois.

Reimplemented from ExplicitSystem.

Definition at line 486 of file fem_system.C.

References MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), DifferentiableSystem::assemble_qoi_sides, FEMSystem::build_context(), DifferentiableSystem::compute_internal_sides, FEMContext::elem, DiffContext::elem_qoi, DifferentiableSystem::element_qoi(), System::get_mesh(), QoISet::has_index(), mesh, Elem::n_sides(), Elem::neighbor(), System::qoi, FEMContext::reinit(), FEMContext::side, FEMContext::side_fe, DifferentiableSystem::side_qoi(), and System::update().

00487 {
00488   START_LOG("assemble_qoi()", "FEMSystem");
00489 
00490   const MeshBase& mesh = this->get_mesh();
00491 
00492   this->update();
00493 
00494   AutoPtr<DiffContext> con = this->build_context();
00495   FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
00496 
00497   // the quantity of interest is assumed to be a sum of element and
00498   // side terms
00499   for (unsigned int i=0; i != qoi.size(); ++i)
00500     if (qoi_indices.has_index(i))
00501       qoi[i] = 0;
00502 
00503   // Loop over every active mesh element on this processor
00504   MeshBase::const_element_iterator el =
00505     mesh.active_local_elements_begin();
00506   const MeshBase::const_element_iterator end_el =
00507     mesh.active_local_elements_end();
00508 
00509   for ( ; el != end_el; ++el)
00510     {
00511       _femcontext.reinit(*this, *el);
00512 
00513       this->element_qoi(_femcontext, qoi_indices);
00514 
00515       for (_femcontext.side = 0;
00516            _femcontext.side != _femcontext.elem->n_sides();
00517            ++_femcontext.side)
00518         {
00519           // Don't compute on non-boundary sides unless requested
00520           if (!assemble_qoi_sides ||
00521               (!compute_internal_sides &&
00522                _femcontext.elem->neighbor(_femcontext.side) != NULL))
00523             continue;
00524 
00525           std::map<FEType, FEBase *>::iterator fe_end =
00526             _femcontext.side_fe.end();
00527           for (std::map<FEType, FEBase *>::iterator i =
00528             _femcontext.side_fe.begin();
00529                i != fe_end; ++i)
00530             {
00531               i->second->reinit(_femcontext.elem, _femcontext.side);
00532             }
00533 
00534           this->side_qoi(_femcontext, qoi_indices);
00535         }
00536     }
00537 
00538   Parallel::sum(_femcontext.elem_qoi);
00539 
00540   this->qoi = _femcontext.elem_qoi;
00541 
00542   STOP_LOG("assemble_qoi()", "FEMSystem");
00543 }

void FEMSystem::assemble_qoi_derivative ( const QoISet indices = QoISet()  )  [virtual, inherited]

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 ExplicitSystem.

Definition at line 547 of file fem_system.C.

References MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), System::add_adjoint_rhs(), DifferentiableSystem::assemble_qoi_sides, FEMSystem::build_context(), DifferentiableSystem::compute_internal_sides, DofMap::constrain_element_vector(), DiffContext::dof_indices, FEMContext::elem, DiffContext::elem_qoi_derivative, DifferentiableSystem::element_qoi_derivative(), System::get_adjoint_rhs(), System::get_dof_map(), System::get_mesh(), QoISet::has_index(), mesh, Elem::n_sides(), Elem::neighbor(), System::qoi, FEMContext::reinit(), FEMContext::side, FEMContext::side_fe, DifferentiableSystem::side_qoi_derivative(), and System::update().

00548 {
00549   START_LOG("assemble_qoi_derivative()", "FEMSystem");
00550 
00551   const MeshBase& mesh = this->get_mesh();
00552 
00553   this->update();
00554 
00555   AutoPtr<DiffContext> con = this->build_context();
00556   FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
00557 
00558   // The quantity of interest derivative assembly accumulates on
00559   // initially zero vectors
00560   for (unsigned int i=0; i != qoi.size(); ++i)
00561     if (qoi_indices.has_index(i))
00562       this->add_adjoint_rhs(i).zero();
00563 
00564   // Loop over every active mesh element on this processor
00565   MeshBase::const_element_iterator el =
00566     mesh.active_local_elements_begin();
00567   const MeshBase::const_element_iterator end_el =
00568     mesh.active_local_elements_end();
00569 
00570   for ( ; el != end_el; ++el)
00571     {
00572       _femcontext.reinit(*this, *el);
00573 
00574       this->element_qoi_derivative(_femcontext, qoi_indices);
00575 
00576       for (_femcontext.side = 0;
00577            _femcontext.side != _femcontext.elem->n_sides();
00578            ++_femcontext.side)
00579         {
00580           // Don't compute on non-boundary sides unless requested
00581           if (!assemble_qoi_sides ||
00582               (!compute_internal_sides &&
00583                _femcontext.elem->neighbor(_femcontext.side) != NULL))
00584             continue;
00585 
00586           std::map<FEType, FEBase *>::iterator fe_end =
00587             _femcontext.side_fe.end();
00588           for (std::map<FEType, FEBase *>::iterator i =
00589             _femcontext.side_fe.begin();
00590                i != fe_end; ++i)
00591             {
00592               i->second->reinit(_femcontext.elem, _femcontext.side);
00593             }
00594 
00595           this->side_qoi_derivative(_femcontext, qoi_indices);
00596         }
00597 
00598       // We need some unmodified indices to use for constraining
00599       // multiple vector
00600       // FIXME - there should be a DofMap::constrain_element_vectors
00601       // to do this more efficiently
00602       std::vector<unsigned int> original_dofs = _femcontext.dof_indices;
00603 
00604       for (unsigned int i=0; i != qoi.size(); ++i)
00605         if (qoi_indices.has_index(i))
00606           {
00607             _femcontext.dof_indices = original_dofs;
00608             this->get_dof_map().constrain_element_vector
00609               (_femcontext.elem_qoi_derivative[i], _femcontext.dof_indices, false);
00610 
00611             this->get_adjoint_rhs(i).add_vector
00612               (_femcontext.elem_qoi_derivative[i], _femcontext.dof_indices);
00613           }
00614     }
00615 
00616   STOP_LOG("assemble_qoi_derivative()", "FEMSystem");
00617 }

void ImplicitSystem::assemble_residual_derivatives ( const ParameterVector parameters  )  [virtual, inherited]

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 System.

Definition at line 622 of file implicit_system.C.

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

Referenced by ImplicitSystem::sensitivity_solve().

00623 {
00624   const unsigned int Np = parameters.size();
00625   Real deltap = TOLERANCE;
00626 
00627   for (unsigned int p=0; p != Np; ++p)
00628     {
00629       NumericVector<Number> &sensitivity_rhs = this->add_sensitivity_rhs(p);
00630 
00631       // Approximate -(partial R / partial p) by
00632       // (R(p-dp) - R(p+dp)) / (2*dp)
00633 
00634       Number old_parameter = *parameters[p];
00635       *parameters[p] -= deltap;
00636 
00637       this->assembly(true, false);
00638       this->rhs->close();
00639       sensitivity_rhs = *this->rhs;
00640 
00641       *parameters[p] = old_parameter + deltap;
00642 
00643       this->assembly(true, false);
00644       this->rhs->close();
00645 
00646       sensitivity_rhs -= *this->rhs;
00647       sensitivity_rhs /= (2*deltap);
00648       sensitivity_rhs.close();
00649 
00650       *parameters[p] = old_parameter;
00651     }
00652 }

void FEMSystem::assembly ( bool  get_residual,
bool  get_jacobian 
) [virtual, inherited]

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 DifferentiableSystem.

Definition at line 78 of file fem_system.C.

References MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), DenseMatrix< T >::add(), FEMSystem::build_context(), DifferentiableSystem::compute_internal_sides, DofMap::constrain_element_matrix(), DofMap::constrain_element_matrix_and_vector(), DofMap::constrain_element_vector(), DiffContext::dof_indices, FEMContext::elem, DiffContext::elem_jacobian, DiffContext::elem_residual, FEMContext::elem_side_fe_reinit(), AutoPtr< Tp >::get(), System::get_dof_map(), System::get_mesh(), DofObject::id(), DenseMatrix< T >::l1_norm(), ImplicitSystem::matrix, std::max(), mesh, Elem::n_sides(), Elem::neighbor(), FEMSystem::numerical_elem_jacobian(), FEMSystem::numerical_side_jacobian(), DifferentiableSystem::print_element_jacobians, DifferentiableSystem::print_jacobian_norms, DifferentiableSystem::print_jacobians, DifferentiableSystem::print_residual_norms, DifferentiableSystem::print_residuals, DifferentiableSystem::print_solution_norms, DifferentiableSystem::print_solutions, FEMContext::reinit(), ExplicitSystem::rhs, FEMContext::side, System::solution, DenseMatrix< T >::swap(), DifferentiableSystem::time_solver, System::update(), and FEMSystem::verify_analytic_jacobians.

Referenced by continuation_solve(), and solve_tangent().

00079 {
00080   libmesh_assert(get_residual || get_jacobian);
00081   std::string log_name;
00082   if (get_residual && get_jacobian)
00083     log_name = "assembly()";
00084   else if (get_residual)
00085     log_name = "assembly(get_residual)";
00086   else
00087     log_name = "assembly(get_jacobian)";
00088 
00089   START_LOG(log_name, "FEMSystem");
00090 
00091   const MeshBase& mesh = this->get_mesh();
00092 
00093 //  this->get_vector("_nonlinear_solution").localize
00094 //    (*current_local_nonlinear_solution,
00095 //     dof_map.get_send_list());
00096   this->update();
00097 
00098   if (print_solution_norms)
00099     {
00100 //      this->get_vector("_nonlinear_solution").close();
00101       this->solution->close();
00102       std::cout << "|U| = "
00103 //                << this->get_vector("_nonlinear_solution").l1_norm()
00104                 << this->solution->l1_norm()
00105                 << std::endl;
00106     }
00107   if (print_solutions)
00108     {
00109       unsigned int old_precision = std::cout.precision();
00110       std::cout.precision(16);
00111 //      std::cout << "U = [" << this->get_vector("_nonlinear_solution")
00112       std::cout << "U = [" << *(this->solution)
00113                 << "];" << std::endl;
00114       std::cout.precision(old_precision);
00115     }
00116 
00117   // Is this definitely necessary? [RHS]
00118   if (get_jacobian)
00119     matrix->zero();
00120   if (get_residual)
00121     rhs->zero();
00122 
00123   // Stupid C++ lets you set *Real* verify_analytic_jacobians = true!
00124   if (verify_analytic_jacobians > 0.5)
00125     {
00126       std::cerr << "WARNING!  verify_analytic_jacobians was set "
00127                 << "to absurdly large value of "
00128                 << verify_analytic_jacobians << std::endl;
00129       std::cerr << "Resetting to 1e-6!" << std::endl;
00130       verify_analytic_jacobians = 1e-6;
00131     }
00132 
00133   // In time-dependent problems, the nonlinear function we're trying
00134   // to solve at each timestep may depend on the particular solver
00135   // we're using
00136   libmesh_assert (time_solver.get() != NULL);
00137 
00138   AutoPtr<DiffContext> con = this->build_context();
00139   FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);
00140 
00141   // Build the residual and jacobian contributions on every active
00142   // mesh element on this processor
00143   MeshBase::const_element_iterator el =
00144     mesh.active_local_elements_begin();
00145   const MeshBase::const_element_iterator end_el =
00146     mesh.active_local_elements_end();
00147 
00148   for ( ; el != end_el; ++el)
00149     {
00150       _femcontext.reinit(*this, *el);
00151 
00152       bool jacobian_computed =
00153         time_solver->element_residual(get_jacobian, _femcontext);
00154 
00155       // Compute a numeric jacobian if we have to
00156       if (get_jacobian && !jacobian_computed)
00157         {
00158           // Make sure we didn't compute a jacobian and lie about it
00159           libmesh_assert(_femcontext.elem_jacobian.l1_norm() == 0.0);
00160           // Logging of numerical jacobians is done separately
00161           this->numerical_elem_jacobian(_femcontext);
00162         }
00163 
00164       // Compute a numeric jacobian if we're asked to verify the
00165       // analytic jacobian we got
00166       if (get_jacobian && jacobian_computed &&
00167           verify_analytic_jacobians != 0.0)
00168         {
00169           DenseMatrix<Number> analytic_jacobian(_femcontext.elem_jacobian);
00170 
00171           _femcontext.elem_jacobian.zero();
00172           // Logging of numerical jacobians is done separately
00173           this->numerical_elem_jacobian(_femcontext);
00174 
00175           Real analytic_norm = analytic_jacobian.l1_norm();
00176           Real numerical_norm = _femcontext.elem_jacobian.l1_norm();
00177 
00178           // If we can continue, we'll probably prefer the analytic jacobian
00179           analytic_jacobian.swap(_femcontext.elem_jacobian);
00180 
00181           // The matrix "analytic_jacobian" will now hold the error matrix
00182           analytic_jacobian.add(-1.0, _femcontext.elem_jacobian);
00183           Real error_norm = analytic_jacobian.l1_norm();
00184 
00185           Real relative_error = error_norm /
00186                                 std::max(analytic_norm, numerical_norm);
00187 
00188           if (relative_error > verify_analytic_jacobians)
00189             {
00190               std::cerr << "Relative error " << relative_error
00191                         << " detected in analytic jacobian on element "
00192                         << _femcontext.elem->id() << '!' << std::endl;
00193 
00194               unsigned int old_precision = std::cout.precision();
00195               std::cout.precision(16);
00196               std::cout << "J_analytic " << _femcontext.elem->id() << " = "
00197                         << _femcontext.elem_jacobian << std::endl;
00198               analytic_jacobian.add(1.0, _femcontext.elem_jacobian);
00199               std::cout << "J_numeric " << _femcontext.elem->id() << " = "
00200                         << analytic_jacobian << std::endl;
00201 
00202               std::cout.precision(old_precision);
00203 
00204               libmesh_error();
00205             }
00206         }
00207 
00208       for (_femcontext.side = 0;
00209            _femcontext.side != _femcontext.elem->n_sides();
00210            ++_femcontext.side)
00211         {
00212           // Don't compute on non-boundary sides unless requested
00213           if (!compute_internal_sides && 
00214               _femcontext.elem->neighbor(_femcontext.side) != NULL)
00215             continue;
00216 
00217           // Any mesh movement has already been done (and restored,
00218           // if the TimeSolver isn't broken), but
00219           // reinitializing the side FE objects is still necessary
00220           _femcontext.elem_side_fe_reinit();
00221 
00222           DenseMatrix<Number> old_jacobian;
00223           // If we're in DEBUG mode, we should always verify that the
00224           // user's side_residual function doesn't alter our existing
00225           // jacobian and then lie about it
00226 #ifndef DEBUG
00227           // Even if we're not in DEBUG mode, when we're verifying
00228           // analytic jacobians we'll want to verify each side's
00229           // jacobian contribution separately
00230           if (verify_analytic_jacobians != 0.0 && get_jacobian)
00231 #endif // ifndef DEBUG
00232             {
00233               old_jacobian = _femcontext.elem_jacobian;
00234               _femcontext.elem_jacobian.zero();
00235             }
00236           jacobian_computed =
00237             time_solver->side_residual(get_jacobian, _femcontext);
00238 
00239           // Compute a numeric jacobian if we have to
00240           if (get_jacobian && !jacobian_computed)
00241             {
00242               // In DEBUG mode, we've already set elem_jacobian == 0,
00243               // so we can make sure side_residual didn't compute a
00244               // jacobian and lie about it
00245 #ifdef DEBUG
00246               libmesh_assert(_femcontext.elem_jacobian.l1_norm() == 0.0);
00247 #endif
00248               // Logging of numerical jacobians is done separately
00249               this->numerical_side_jacobian(_femcontext);
00250 
00251               // If we're in DEBUG mode or if
00252               // verify_analytic_jacobians is on, we've moved
00253               // elem_jacobian's accumulated values into old_jacobian.
00254               // Now let's add them back.
00255 #ifndef DEBUG
00256               if (verify_analytic_jacobians != 0.0)
00257 #endif // ifndef DEBUG
00258                 _femcontext.elem_jacobian += old_jacobian;
00259             }
00260 
00261           // Compute a numeric jacobian if we're asked to verify the
00262           // analytic jacobian we got
00263           if (get_jacobian && jacobian_computed &&
00264               verify_analytic_jacobians != 0.0)
00265             {
00266               DenseMatrix<Number> analytic_jacobian(_femcontext.elem_jacobian);
00267 
00268               _femcontext.elem_jacobian.zero();
00269               // Logging of numerical jacobians is done separately
00270               this->numerical_side_jacobian(_femcontext);
00271 
00272               Real analytic_norm = analytic_jacobian.l1_norm();
00273               Real numerical_norm = _femcontext.elem_jacobian.l1_norm();
00274 
00275               // If we can continue, we'll probably prefer the analytic jacobian
00276               analytic_jacobian.swap(_femcontext.elem_jacobian);
00277 
00278               // The matrix "analytic_jacobian" will now hold the error matrix
00279               analytic_jacobian.add(-1.0, _femcontext.elem_jacobian);
00280               Real error_norm = analytic_jacobian.l1_norm();
00281 
00282               Real relative_error = error_norm /
00283                                     std::max(analytic_norm, numerical_norm);
00284 
00285               if (relative_error > verify_analytic_jacobians)
00286                 {
00287                   std::cerr << "Relative error " << relative_error
00288                             << " detected in analytic jacobian on element "
00289                             << _femcontext.elem->id()
00290                             << ", side "
00291                             << _femcontext.side << '!' << std::endl;
00292 
00293                   unsigned int old_precision = std::cout.precision();
00294                   std::cout.precision(16);
00295                   std::cout << "J_analytic " << _femcontext.elem->id() << " = "
00296                             << _femcontext.elem_jacobian << std::endl;
00297                   analytic_jacobian.add(1.0, _femcontext.elem_jacobian);
00298                   std::cout << "J_numeric " << _femcontext.elem->id() << " = "
00299                             << analytic_jacobian << std::endl;
00300                   std::cout.precision(old_precision);
00301 
00302                   libmesh_error();
00303                 }
00304               // Once we've verified a side, we'll want to add back the
00305               // rest of the accumulated jacobian
00306               _femcontext.elem_jacobian += old_jacobian;
00307             }
00308           // In DEBUG mode, we've set elem_jacobian == 0, and we
00309           // may still need to add the old jacobian back
00310 #ifdef DEBUG
00311           if (get_jacobian && jacobian_computed &&
00312               verify_analytic_jacobians == 0.0)
00313             {
00314               _femcontext.elem_jacobian += old_jacobian;
00315             }
00316 #endif // ifdef DEBUG
00317         }
00318 
00319 #ifdef LIBMESH_ENABLE_AMR
00320       // We turn off the asymmetric constraint application;
00321       // enforce_constraints_exactly() should be called in the solver
00322       if (get_residual && get_jacobian)
00323         this->get_dof_map().constrain_element_matrix_and_vector
00324           (_femcontext.elem_jacobian, _femcontext.elem_residual,
00325            _femcontext.dof_indices, false);
00326       else if (get_residual)
00327         this->get_dof_map().constrain_element_vector
00328           (_femcontext.elem_residual, _femcontext.dof_indices, false);
00329       else if (get_jacobian)
00330         this->get_dof_map().constrain_element_matrix
00331           (_femcontext.elem_jacobian, _femcontext.dof_indices, false);
00332 #endif // #ifdef LIBMESH_ENABLE_AMR
00333 
00334       if (get_jacobian && print_element_jacobians)
00335         {
00336           unsigned int old_precision = std::cout.precision();
00337           std::cout.precision(16);
00338           std::cout << "J_elem " << _femcontext.elem->id() << " = "
00339                     << _femcontext.elem_jacobian << std::endl;
00340           std::cout.precision(old_precision);
00341         }
00342 
00343       if (get_jacobian)
00344         this->matrix->add_matrix (_femcontext.elem_jacobian,
00345                                   _femcontext.dof_indices);
00346       if (get_residual)
00347         this->rhs->add_vector (_femcontext.elem_residual,
00348                                _femcontext.dof_indices);
00349     }
00350 
00351 
00352   if (get_residual && print_residual_norms)
00353     {
00354       this->rhs->close();
00355       std::cout << "|F| = " << this->rhs->l1_norm() << std::endl;
00356     }
00357   if (get_residual && print_residuals)
00358     {
00359       this->rhs->close();
00360       unsigned int old_precision = std::cout.precision();
00361       std::cout.precision(16);
00362       std::cout << "F = [" << *(this->rhs) << "];" << std::endl;
00363       std::cout.precision(old_precision);
00364     }
00365   if (get_jacobian && print_jacobian_norms)
00366     {
00367       this->matrix->close();
00368       std::cout << "|J| = " << this->matrix->l1_norm() << std::endl;
00369     }
00370   if (get_jacobian && print_jacobians)
00371     {
00372       this->matrix->close();
00373       unsigned int old_precision = std::cout.precision();
00374       std::cout.precision(16);
00375       std::cout << "J = [" << *(this->matrix) << "];" << std::endl;
00376       std::cout.precision(old_precision);
00377     }
00378   STOP_LOG(log_name, "FEMSystem");
00379 }

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

References System::_assemble_system.

01318 {
01319   libmesh_assert (fptr != NULL);
01320   
01321   _assemble_system = fptr;  
01322 }

void System::attach_constraint_function ( void   fptrEquationSystems &es,const std::string &name  )  [inherited]

Register a user function for imposing constraints.

Definition at line 1326 of file system.C.

References System::_constrain_system.

01328 {
01329   libmesh_assert (fptr != NULL);
01330   
01331   _constrain_system = fptr;  
01332 }

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

References System::_init_system.

01308 {
01309   libmesh_assert (fptr != NULL);
01310   
01311   _init_system = fptr;
01312 }

void 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

void 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

AutoPtr< DiffContext > FEMSystem::build_context (  )  [virtual, inherited]

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 DifferentiableSystem.

Definition at line 742 of file fem_system.C.

Referenced by FEMSystem::assemble_qoi(), FEMSystem::assemble_qoi_derivative(), FEMSystem::assembly(), FEMSystem::mesh_position_get(), FEMSystem::mesh_position_set(), and FEMSystem::postprocess().

00743 {
00744   return AutoPtr<DiffContext>(new FEMContext(*this));
00745 }

Real System::calculate_norm ( 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 1077 of file system.C.

References System::_dof_map, MeshBase::active_local_elements_begin(), MeshBase::active_local_elements_end(), TypeTensor< T >::add_scaled(), TypeVector< T >::add_scaled(), FEBase::build(), FEType::default_quadrature_rule(), dim, libMeshEnums::DISCRETE_L1, libMeshEnums::DISCRETE_L2, libMeshEnums::DISCRETE_L_INF, System::discrete_var_norm(), DofMap::dof_indices(), AutoPtr< Tp >::get(), System::get_dof_map(), System::get_mesh(), libMeshEnums::H1, libMeshEnums::H1_SEMINORM, libMeshEnums::H2, libMeshEnums::H2_SEMINORM, SystemNorm::is_discrete(), NumericVector< T >::l1_norm(), libMeshEnums::L2, NumericVector< T >::l2_norm(), libmesh_norm(), NumericVector< T >::linfty_norm(), NumericVector< T >::localize(), MeshBase::mesh_dimension(), System::n_vars(), libMeshEnums::SERIAL, NumericVector< T >::size(), TypeTensor< T >::size_sq(), TypeVector< T >::size_sq(), SystemNorm::type(), DofMap::variable_type(), SystemNorm::weight(), and SystemNorm::weight_sq().

01079 {
01080   // This function must be run on all processors at once
01081   parallel_only();
01082 
01083   START_LOG ("calculate_norm()", "System");
01084 
01085   // Zero the norm before summation
01086   Real v_norm = 0.;
01087 
01088   if (norm.is_discrete())
01089     {
01090       STOP_LOG ("calculate_norm()", "System");
01091       //Check to see if all weights are 1.0
01092       unsigned int check_var = 0;
01093       for (; check_var != this->n_vars(); ++check_var)
01094         if(norm.weight(check_var) != 1.0)
01095           break;
01096 
01097       //All weights were 1.0 so just do the full vector discrete norm
01098       if(check_var == this->n_vars())
01099         {
01100           FEMNormType norm_type = norm.type(0);
01101           
01102           if(norm_type == DISCRETE_L1)
01103             return v.l1_norm();
01104           if(norm_type == DISCRETE_L2)
01105             return v.l2_norm();
01106           if(norm_type == DISCRETE_L_INF)
01107             return v.linfty_norm();
01108           else
01109             libmesh_error();
01110         }
01111 
01112       for (unsigned int var=0; var != this->n_vars(); ++var)
01113         {
01114           // Skip any variables we don't need to integrate
01115           if (norm.weight(var) == 0.0)
01116             continue;
01117 
01118           v_norm += norm.weight(var) * discrete_var_norm(v, var, norm.type(var));
01119         }
01120 
01121       return v_norm;
01122     }
01123 
01124   // Localize the potentially parallel vector
01125   AutoPtr<NumericVector<Number> > local_v = NumericVector<Number>::build();
01126   local_v->init(v.size(), true, SERIAL);
01127   v.localize (*local_v, _dof_map->get_send_list()); 
01128 
01129   unsigned int dim = this->get_mesh().mesh_dimension();
01130 
01131   // Loop over all variables
01132   for (unsigned int var=0; var != this->n_vars(); ++var)
01133     {
01134       // Skip any variables we don't need to integrate
01135       if (norm.weight(var) == 0.0)
01136         continue;
01137 
01138       const FEType& fe_type = this->get_dof_map().variable_type(var);
01139       AutoPtr<QBase> qrule =
01140         fe_type.default_quadrature_rule (dim);
01141       AutoPtr<FEBase> fe
01142         (FEBase::build(dim, fe_type));
01143       fe->attach_quadrature_rule (qrule.get());
01144 
01145       const std::vector<Real>&               JxW = fe->get_JxW();
01146       const std::vector<std::vector<Real> >* phi = NULL;
01147       if (norm.type(var) == H1 ||
01148           norm.type(var) == H2 ||
01149           norm.type(var) == L2)
01150         phi = &(fe->get_phi());
01151 
01152       const std::vector<std::vector<RealGradient> >* dphi = NULL;
01153       if (norm.type(var) == H1 ||
01154           norm.type(var) == H2 ||
01155           norm.type(var) == H1_SEMINORM)
01156         dphi = &(fe->get_dphi());
01157 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
01158       const std::vector<std::vector<RealTensor> >*   d2phi = NULL;
01159       if (norm.type(var) == H2 ||
01160           norm.type(var) == H2_SEMINORM)
01161         d2phi = &(fe->get_d2phi());
01162 #endif
01163 
01164       std::vector<unsigned int> dof_indices;
01165 
01166       // Begin the loop over the elements
01167       MeshBase::const_element_iterator       el     =
01168         this->get_mesh().active_local_elements_begin();
01169       const MeshBase::const_element_iterator end_el =
01170         this->get_mesh().active_local_elements_end();
01171 
01172       for ( ; el != end_el; ++el)
01173         {
01174           const Elem* elem = *el;
01175 
01176           fe->reinit (elem);
01177 
01178           this->get_dof_map().dof_indices (elem, dof_indices, var);
01179 
01180           const unsigned int n_qp = qrule->n_points();
01181 
01182           const unsigned int n_sf = dof_indices.size();
01183 
01184           // Begin the loop over the Quadrature points.
01185           for (unsigned int qp=0; qp<n_qp; qp++)
01186             {
01187               if (norm.type(var) == H1 ||
01188                   norm.type(var) == H2 ||
01189                   norm.type(var) == L2)
01190                 {
01191                   Number u_h = 0.;
01192                   for (unsigned int i=0; i != n_sf; ++i)
01193                     u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
01194                   v_norm += norm.weight_sq(var) *
01195                             JxW[qp] * libmesh_norm(u_h);
01196                 }
01197 
01198               if (norm.type(var) == H1 ||
01199                   norm.type(var) == H2 ||
01200                   norm.type(var) == H1_SEMINORM)
01201                 {
01202                   Gradient grad_u_h;
01203                   for (unsigned int i=0; i != n_sf; ++i)
01204                     grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
01205                   v_norm += norm.weight_sq(var) *
01206                             JxW[qp] * grad_u_h.size_sq();
01207                 }
01208 
01209 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
01210               if (norm.type(var) == H2 ||
01211                   norm.type(var) == H2_SEMINORM)
01212                 {
01213                   Tensor hess_u_h;
01214                   for (unsigned int i=0; i != n_sf; ++i)
01215                     hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
01216                   v_norm += norm.weight_sq(var) *
01217                             JxW[qp] * hess_u_h.size_sq();
01218                 }
01219 #endif
01220             }
01221         }
01222     }
01223 
01224   Parallel::sum(v_norm);
01225 
01226   STOP_LOG ("calculate_norm()", "System");
01227 
01228   return std::sqrt(v_norm);
01229 }

Real System::calculate_norm ( 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, H0, H1)

Definition at line 1056 of file system.C.

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

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

01059 {
01060   //short circuit to save time
01061   if(norm_type == DISCRETE_L1 ||
01062      norm_type == DISCRETE_L2 ||
01063      norm_type == DISCRETE_L_INF)
01064     return discrete_var_norm(v,var,norm_type);
01065 
01066   // Not a discrete norm
01067   std::vector<FEMNormType> norms(this->n_vars(), L2);
01068   std::vector<Real> weights(this->n_vars(), 0.0);
01069   norms[var] = norm_type;
01070   weights[var] = 1.0;
01071   Real val = this->calculate_norm(v, SystemNorm(norms, weights));
01072   return val;
01073 }

void ContinuationSystem::clear (  )  [virtual]

Clear all the data structures associated with the system.

Reimplemented from FEMSystem.

Definition at line 74 of file continuation_system.C.

References FEMSystem::clear().

Referenced by ~ContinuationSystem().

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

bool System::compare ( const System other_system,
const Real  threshold,
const bool  verbose 
) const [virtual, inherited]

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

References System::_can_add_vectors, System::_sys_name, System::_vectors, System::get_vector(), System::n_vectors(), System::name(), and System::solution.

Referenced by EquationSystems::compare().

00402 {
00403   // we do not care for matrices, but for vectors
00404   libmesh_assert (!_can_add_vectors);
00405   libmesh_assert (!other_system._can_add_vectors);
00406 
00407   if (verbose)
00408     {
00409       std::cout << "  Systems \"" << _sys_name << "\"" << std::endl;
00410       std::cout << "   comparing matrices not supported." << std::endl;
00411       std::cout << "   comparing names...";
00412     }
00413 
00414   // compare the name: 0 means identical
00415   const int name_result = _sys_name.compare(other_system.name());
00416   if (verbose)
00417     {
00418       if (name_result == 0)
00419         std::cout << " identical." << std::endl;
00420       else
00421         std::cout << "  names not identical." << std::endl;
00422       std::cout << "   comparing solution vector...";
00423     }
00424 
00425 
00426   // compare the solution: -1 means identical
00427   const int solu_result = solution->compare (*other_system.solution.get(),
00428                                              threshold);
00429 
00430   if (verbose)
00431     {
00432       if (solu_result == -1)
00433         std::cout << " identical up to threshold." << std::endl;
00434       else
00435         std::cout << "  first difference occured at index = " 
00436                   << solu_result << "." << std::endl;
00437     }
00438 
00439 
00440   // safety check, whether we handle at least the same number
00441   // of vectors
00442   std::vector<int> ov_result;
00443 
00444   if (this->n_vectors() != other_system.n_vectors())
00445     {
00446       if (verbose)
00447         {
00448           std::cout << "   Fatal difference. This system handles " 
00449                     << this->n_vectors() << " add'l vectors," << std::endl
00450                     << "   while the other system handles "
00451                     << other_system.n_vectors() 
00452                     << " add'l vectors." << std::endl
00453                     << "   Aborting comparison." << std::endl;
00454         }
00455       return false;
00456     }
00457   else if (this->n_vectors() == 0)
00458     {
00459       // there are no additional vectors...
00460       ov_result.clear ();
00461     }
00462   else
00463     {
00464       // compare other vectors
00465       for (const_vectors_iterator pos = _vectors.begin();
00466            pos != _vectors.end(); ++pos)
00467         {
00468           if (verbose)
00469               std::cout << "   comparing vector \""
00470                         << pos->first << "\" ...";
00471 
00472           // assume they have the same name
00473           const NumericVector<Number>& other_system_vector = 
00474               other_system.get_vector(pos->first);
00475 
00476           ov_result.push_back(pos->second->compare (other_system_vector,
00477                                                     threshold));
00478 
00479           if (verbose)
00480             {
00481               if (ov_result[ov_result.size()-1] == -1)
00482                 std::cout << " identical up to threshold." << std::endl;
00483               else
00484                 std::cout << " first difference occured at" << std::endl
00485                           << "   index = " << ov_result[ov_result.size()-1] << "." << std::endl;
00486             }
00487 
00488         }
00489 
00490     } // finished comparing additional vectors
00491 
00492 
00493   bool overall_result;
00494        
00495   // sum up the results
00496   if ((name_result==0) && (solu_result==-1))
00497     {
00498       if (ov_result.size()==0)
00499         overall_result = true;
00500       else
00501         {
00502           bool ov_identical;
00503           unsigned int n    = 0;
00504           do
00505             {
00506               ov_identical = (ov_result[n]==-1);
00507               n++;
00508             }
00509           while (ov_identical && n<ov_result.size());
00510           overall_result = ov_identical;
00511         }
00512     }
00513   else
00514     overall_result = false;
00515 
00516   if (verbose)
00517     {
00518       std::cout << "   finished comparisons, ";
00519       if (overall_result)
00520         std::cout << "found no differences." << std::endl << std::endl;
00521       else 
00522         std::cout << "found differences." << std::endl << std::endl;
00523     }
00524           
00525   return overall_result;
00526 }

void 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 352 of file continuation_system.C.

References apply_predictor(), FEMSystem::assembly(), continuation_parameter, continuation_parameter_tolerance, delta_u, dlambda_ds, ds_current, du_ds, G_Lambda, AutoPtr< Tp >::get(), initial_newton_tolerance, initialize_tangent(), libmesh_real(), linear_solver, ImplicitSystem::matrix, max_continuation_parameter, DiffSolver::max_linear_iterations, 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, Utility::pow(), previous_u, quiet, Residual, ExplicitSystem::rhs, rhs_mode, System::solution, solution_tolerance, tangent_initialized, Theta, Theta_LOCA, DifferentiableSystem::time_solver, y, y_old, and z.

00353 {
00354   // Be sure the user has set the continuation parameter pointer
00355   if (!continuation_parameter)
00356     {
00357       std::cerr << "You must set the continuation_parameter pointer "
00358                 << "to a member variable of the derived class, preferably in the "
00359                 << "Derived class's init_data function.  This is how the ContinuationSystem "
00360                 << "updates the continuation parameter."
00361                 << std::endl;
00362       
00363       libmesh_error();
00364     }
00365 
00366   // Use extra precision for all the numbers printed in this function.
00367   unsigned int old_precision = std::cout.precision();
00368   std::cout.precision(16);
00369   std::cout.setf(std::ios_base::scientific);
00370   
00371   // We can't start solving the augmented PDE system unless the tangent
00372   // vectors have been initialized.  This only needs to occur once.
00373   if (!tangent_initialized)
00374     initialize_tangent();
00375 
00376   // Save the old value of -du/dlambda.  This will be used after the Newton iterations
00377   // to compute the angle between previous tangent vectors.  This cosine of this angle is
00378   //
00379   // tau := abs( (du/d(lambda)_i , du/d(lambda)_{i-1}) / (||du/d(lambda)_i|| * ||du/d(lambda)_{i-1}||) )
00380   //
00381   // The scaling factor tau (which should vary between 0 and 1) is used to shrink the step-size ds
00382   // when we are approaching a turning point.  Note that it can only shrink the step size.  
00383   *y_old = *y;
00384   
00385   // Set pointer to underlying Newton solver
00386   if (!newton_solver)
00387     newton_solver = libmesh_cast_ptr<NewtonSolver*> (this->time_solver->diff_solver().get());
00388   
00389   // A pair for catching return values from linear system solves.
00390   std::pair<unsigned int, Real> rval;
00391 
00392   // Convergence flag for the entire arcstep
00393   bool arcstep_converged = false;
00394 
00395   // Begin loop over arcstep reductions.
00396   for (unsigned int ns=0; ns<n_arclength_reductions; ++ns)
00397     {
00398       if (!quiet)
00399         {
00400           std::cout << "Current arclength stepsize, ds_current=" << ds_current << std::endl;
00401           std::cout << "Current parameter value, lambda=" << *continuation_parameter << std::endl;
00402         }
00403       
00404       // Upon exit from the nonlinear loop, the newton_converged flag
00405       // will tell us the convergence status of Newton's method.
00406       bool newton_converged = false;
00407 
00408       // The nonlinear residual before *any* nonlinear steps have been taken.
00409       Real nonlinear_residual_firststep = 0.;
00410 
00411       // The nonlinear residual from the current "k" Newton step, before the Newton step
00412       Real nonlinear_residual_beforestep = 0.;
00413 
00414       // The nonlinear residual from the current "k" Newton step, after the Newton step
00415       Real nonlinear_residual_afterstep = 0.;
00416 
00417       // The linear solver tolerance, can be updated dynamically at each Newton step.
00418       Real current_linear_tolerance = 0.;
00419       
00420       // The nonlinear loop
00421       for (newton_step=0; newton_step<newton_solver->max_nonlinear_iterations; ++newton_step)
00422         {
00423           std::cout << "\n === Starting Newton step " << newton_step << " ===" << std::endl;
00424           
00425           // Set the linear system solver tolerance
00426 //        // 1.) Set the current linear tolerance based as a multiple of the current residual of the system.
00427 //        const Real residual_multiple = 1.e-4;
00428 //        Real current_linear_tolerance = residual_multiple*nonlinear_residual_beforestep;
00429 
00430 //        // But if the current residual isn't small, don't let the solver exit with zero iterations!
00431 //        if (current_linear_tolerance > 1.)
00432 //          current_linear_tolerance = residual_multiple;
00433 
00434           // 2.) Set the current linear tolerance based on the method based on technique of Eisenstat & Walker.
00435           if (newton_step==0)
00436             {
00437               // At first step, only try reducing the residual by a small amount
00438               current_linear_tolerance = initial_newton_tolerance;//0.01;
00439             }
00440 
00441           else
00442             {
00443               // The new tolerance is based on the ratio of the most recent tolerances
00444               const Real alp=0.5*(1.+std::sqrt(5.));
00445               const Real gam=0.9;
00446 
00447               libmesh_assert (nonlinear_residual_beforestep != 0.0);
00448               libmesh_assert (nonlinear_residual_afterstep != 0.0);
00449 
00450               current_linear_tolerance = std::min(gam*std::pow(nonlinear_residual_afterstep/nonlinear_residual_beforestep, alp),
00451                                                   current_linear_tolerance*current_linear_tolerance
00452                                                   );
00453 
00454               // Don't let it get ridiculously small!!
00455               if (current_linear_tolerance < 1.e-12)
00456                 current_linear_tolerance = 1.e-12;
00457             }
00458 
00459           if (!quiet)
00460             std::cout << "Using current_linear_tolerance=" << current_linear_tolerance << std::endl;
00461 
00462           
00463           // Assemble the residual (and Jacobian).
00464           rhs_mode = Residual;
00465           assembly(true,   // Residual
00466                    true); // Jacobian
00467           rhs->close();
00468 
00469           // Save the current nonlinear residual.  We don't need to recompute the residual unless
00470           // this is the first step, since it was already computed as part of the convergence check
00471           // at the end of the last loop iteration.
00472           if (newton_step==0)
00473             {
00474               nonlinear_residual_beforestep = rhs->l2_norm();
00475 
00476               // Store the residual before any steps have been taken.  This will *not*
00477               // be updated at each step, and can be used to see if any progress has
00478               // been made from the initial residual at later steps.
00479               nonlinear_residual_firststep = nonlinear_residual_beforestep;
00480 
00481               const Real old_norm_u = solution->l2_norm();
00482               std::cout << "  (before step) ||R||_{L2} = " << nonlinear_residual_beforestep << std::endl;
00483               std::cout << "  (before step) ||R||_{L2}/||u|| = " << nonlinear_residual_beforestep / old_norm_u << std::endl;
00484 
00485               // In rare cases (very small arcsteps), it's possible that the residual is
00486               // already below our absolute linear tolerance.
00487               if (nonlinear_residual_beforestep  < solution_tolerance)
00488                 {
00489                   if (!quiet)
00490                     std::cout << "Initial guess satisfied linear tolerance, exiting with zero Newton iterations!" << std::endl;
00491 
00492                   // Since we go straight from here to the solve of the next tangent, we
00493                   // have to close the matrix before it can be assembled again.
00494                   matrix->close();
00495                   newton_converged=true;
00496                   break; // out of Newton iterations, with newton_converged=true
00497                 }
00498             }
00499 
00500           else
00501             {
00502               nonlinear_residual_beforestep = nonlinear_residual_afterstep;
00503             }
00504 
00505           
00506           // Solve the linear system G_u*z = G
00507           // Initial guess?
00508           z->zero(); // It seems to be extremely important to zero z here, otherwise the solver quits early.
00509           z->close();
00510           
00511           // It's possible that we have selected the current_linear_tolerance so large that
00512           // a guess of z=zero yields a linear system residual |Az + R| small enough that the
00513           // linear solver exits in zero iterations.  If this happens, we will reduce the
00514           // current_linear_tolerance until the linear solver does at least 1 iteration.
00515           do
00516             {
00517               rval =
00518                 linear_solver->solve(*matrix,
00519                                      *z,
00520                                      *rhs,
00521                                      //1.e-12,
00522                                      current_linear_tolerance,
00523                                      newton_solver->max_linear_iterations);   // max linear iterations
00524               
00525               if (rval.first==0)
00526                 {
00527                   if (newton_step==0)
00528                     {
00529                       std::cout << "Repeating initial solve with smaller linear tolerance!" << std::endl;
00530                       current_linear_tolerance *= initial_newton_tolerance; // reduce the linear tolerance to force the solver to do some work
00531                     }
00532                   else
00533                     {
00534                       // We shouldn't get here ... it means the linear solver did no work on a Newton
00535                       // step other than the first one.  If this happens, we need to think more about our
00536                       // tolerance selection.
00537                       libmesh_error();
00538                     }
00539                 }
00540               
00541             } while (rval.first==0);
00542 
00543           
00544           if (!quiet)
00545             std::cout << "  G_u*z = G solver converged at step "
00546                       << rval.first
00547                       << " linear tolerance = "
00548                       << rval.second
00549                       << "."
00550                       << std::endl;
00551 
00552           // Sometimes (I am not sure why) the linear solver exits after zero iterations.
00553           // Perhaps it is hitting PETSc's divergence tolerance dtol???  If this occurs,
00554           // we should break out of the Newton iteration loop because nothing further is
00555           // going to happen...  Of course if the tolerance is already small enough after
00556           // zero iterations (how can this happen?!) we should not quit.
00557           if ((rval.first == 0) && (rval.second > current_linear_tolerance*nonlinear_residual_beforestep))
00558             {
00559               if (!quiet)
00560                 std::cout << "Linear solver exited in zero iterations!" << std::endl;
00561 
00562               // Try to find out the reason for convergence/divergence
00563               linear_solver->print_converged_reason();
00564               
00565               break; // out of Newton iterations
00566             }
00567           
00568           // Note: need to scale z by -1 since our code always solves Jx=R
00569           // instead of Jx=-R.
00570           z->scale(-1.);
00571           z->close();
00572 
00573 
00574 
00575 
00576 
00577           
00578           // Assemble the G_Lambda vector, skip residual.
00579           rhs_mode = G_Lambda;
00580 
00581           // Assemble both rhs and Jacobian
00582           assembly(true,  // Residual
00583                    false); // Jacobian
00584 
00585           // Not sure if this is really necessary
00586           rhs->close();
00587           const Real yrhsnorm=rhs->l2_norm();
00588           if (yrhsnorm == 0.0)
00589             {
00590               std::cout << "||G_Lambda|| = 0" << std::endl;
00591               libmesh_error();
00592             }
00593 
00594           // We select a tolerance for the y-system which is based on the inexact Newton
00595           // tolerance but scaled by an extra term proportional to the RHS (which is not -> 0 in this case)
00596           const Real ysystemtol=current_linear_tolerance*(nonlinear_residual_beforestep/yrhsnorm);
00597           if (!quiet)
00598             std::cout << "ysystemtol=" << ysystemtol << std::endl;
00599           
00600           // Solve G_u*y = G_{\lambda}
00601           // FIXME: Initial guess?  This is really a solve for -du/dlambda so we could try
00602           // initializing it with the latest approximation to that... du/dlambda ~ du/ds * ds/dlambda
00603           //*y = *solution;
00604           //y->add(-1., *previous_u);
00605           //y->scale(-1. / (*continuation_parameter - old_continuation_parameter)); // Be careful of divide by zero...
00606           //y->close();
00607           
00608           //      const unsigned int max_attempts=1;
00609           // unsigned int attempt=0;
00610           //      do
00611           //        {
00612           //          if (!quiet)
00613           //            std::cout << "Trying to solve tangent system, attempt " << attempt << std::endl;
00614               
00615               rval =
00616                 linear_solver->solve(*matrix,
00617                                      *y,
00618                                      *rhs,
00619                                      //1.e-12, 
00620                                      ysystemtol,
00621                                      newton_solver->max_linear_iterations);   // max linear iterations
00622 
00623               if (!quiet)
00624                 std::cout << "  G_u*y = G_{lambda} solver converged at step "
00625                           << rval.first
00626                           << ", linear tolerance = "
00627                           << rval.second
00628                           << "."
00629                           << std::endl;
00630 
00631               // Sometimes (I am not sure why) the linear solver exits after zero iterations.
00632               // Perhaps it is hitting PETSc's divergence tolerance dtol???  If this occurs,
00633               // we should break out of the Newton iteration loop because nothing further is
00634               // going to happen...
00635               if ((rval.first == 0) && (rval.second > ysystemtol))
00636                 {
00637                   if (!quiet)
00638                     std::cout << "Linear solver exited in zero iterations!" << std::endl;
00639 
00640                   break; // out of Newton iterations
00641                 }
00642               
00643 //            ++attempt;
00644 //          } while ((attempt<max_attempts) && (rval.first==newton_solver->max_linear_iterations));
00645             
00646       
00647       
00648       
00649       
00650           // Compute N, the residual of the arclength constraint eqn.
00651           // Note 1: N(u,lambda,s) := (u-u_{old}, du_ds) + (lambda-lambda_{old}, dlambda_ds) - _ds 
00652           // We temporarily use the delta_u vector as a temporary vector for this calculation.
00653           *delta_u = *solution;
00654           delta_u->add(-1., *previous_u);
00655 
00656           // First part of the arclength constraint
00657           const Number N1 = Theta_LOCA*Theta_LOCA*Theta*delta_u->dot(*du_ds);
00658           const Number N2 = ((*continuation_parameter) - old_continuation_parameter)*dlambda_ds;
00659           const Number N3 = ds_current;
00660 
00661           if (!quiet)
00662             {
00663               std::cout << "  N1=" << N1 << std::endl;
00664               std::cout << "  N2=" << N2 << std::endl;
00665               std::cout << "  N3=" << N3 << std::endl;
00666             }
00667       
00668           // The arclength constraint value 
00669           const Number N = N1+N2-N3;
00670       
00671           if (!quiet)
00672             std::cout << "  N=" << N << std::endl;
00673 
00674           const Number duds_dot_z = du_ds->dot(*z);
00675           const Number duds_dot_y = du_ds->dot(*y);
00676 
00677           //std::cout << "duds_dot_z=" << duds_dot_z << std::endl;
00678           //std::cout << "duds_dot_y=" << duds_dot_y << std::endl;
00679           //std::cout << "dlambda_ds=" << dlambda_ds << std::endl;
00680 
00681           const Number delta_lambda_numerator   = -(N          + Theta_LOCA*Theta_LOCA*Theta*duds_dot_z);
00682           const Number delta_lambda_denominator =  (dlambda_ds - Theta_LOCA*Theta_LOCA*Theta*duds_dot_y);
00683 
00684           libmesh_assert (delta_lambda_denominator != 0.0);
00685 
00686           // Now, we are ready to compute the step delta_lambda
00687           const Number delta_lambda_comp = delta_lambda_numerator /
00688                                            delta_lambda_denominator;
00689           // Lambda is real-valued
00690           const Real delta_lambda = libmesh_real(delta_lambda_comp);
00691 
00692           // Knowing delta_lambda, we are ready to update delta_u
00693           // delta_u = z - delta_lambda*y
00694           delta_u->zero();
00695           delta_u->add(1., *z);
00696           delta_u->add(-delta_lambda, *y);
00697           delta_u->close();
00698 
00699           // Update the system solution and the continuation parameter.
00700           solution->add(1., *delta_u);
00701           solution->close();
00702           *continuation_parameter += delta_lambda;
00703 
00704           // Did the Newton step actually reduce the residual?
00705           rhs_mode = Residual;
00706           assembly(true,   // Residual
00707                    false); // Jacobian
00708           rhs->close();
00709           nonlinear_residual_afterstep = rhs->l2_norm();
00710 
00711           
00712           // In a "normal" Newton step, ||du||/||R|| > 1 since the most recent
00713           // step is where you "just were" and the current residual is where
00714           // you are now.  It can occur that ||du||/||R|| < 1, but these are
00715           // likely not good cases to attempt backtracking (?).
00716           const Real norm_du_norm_R = delta_u->l2_norm() / nonlinear_residual_afterstep;
00717           if (!quiet)
00718             std::cout << "  norm_du_norm_R=" << norm_du_norm_R << std::endl;
00719       
00720       
00721           // Factor to decrease the stepsize by for backtracking
00722           Real newton_stepfactor = 1.;
00723 
00724           const bool attempt_backtracking =
00725             (nonlinear_residual_afterstep > solution_tolerance)
00726             && (nonlinear_residual_afterstep > nonlinear_residual_beforestep)
00727             && (n_backtrack_steps>0)
00728             && (norm_du_norm_R > 1.)
00729             ;
00730     
00731           // If residual is not reduced, do Newton back tracking.  
00732           if (attempt_backtracking)
00733             {
00734               if (!quiet)
00735                 std::cout << "Newton step did not reduce residual." << std::endl;
00736 
00737               // back off the previous step.
00738               solution->add(-1., *delta_u);
00739               solution->close();
00740               *continuation_parameter -= delta_lambda;
00741 
00742               // Backtracking: start cutting the Newton stepsize by halves until
00743               // the new residual is actually smaller...
00744               for (unsigned int backtrack_step=0; backtrack_step<n_backtrack_steps; ++backtrack_step)
00745                 {
00746                   newton_stepfactor *= 0.5;
00747 
00748                   if (!quiet)
00749                     std::cout << "Shrinking step size by " << newton_stepfactor << std::endl;
00750               
00751                   // Take fractional step
00752                   solution->add(newton_stepfactor, *delta_u);
00753                   solution->close();
00754                   *continuation_parameter += newton_stepfactor*delta_lambda;
00755               
00756                   rhs_mode = Residual;
00757                   assembly(true,   // Residual
00758                            false); // Jacobian
00759                   rhs->close();
00760                   nonlinear_residual_afterstep = rhs->l2_norm();
00761 
00762                   if (!quiet)
00763                     std::cout << "At shrink step "
00764                               << backtrack_step
00765                               << ", nonlinear_residual_afterstep="
00766                               << nonlinear_residual_afterstep
00767                               << std::endl;
00768 
00769                   if (nonlinear_residual_afterstep < nonlinear_residual_beforestep)
00770                     {
00771                       if (!quiet)
00772                         std::cout << "Backtracking succeeded!" << std::endl;
00773                       
00774                       break; // out of backtracking loop
00775                     }
00776                   
00777                   else
00778                     {
00779                       // Back off that step
00780                       solution->add(-newton_stepfactor, *delta_u);
00781                       solution->close();
00782                       *continuation_parameter -= newton_stepfactor*delta_lambda;
00783                     }
00784               
00785                   // Save a copy of the solution from before the Newton step.
00786                   //AutoPtr<NumericVector<Number> > prior_iterate = solution->clone();
00787                 }
00788             } // end if (attempte_backtracking)
00789 
00790 
00791           // If we tried backtracking but the residual is still not reduced, print message.
00792           if ((attempt_backtracking) && (nonlinear_residual_afterstep > nonlinear_residual_beforestep))
00793             {
00794               //std::cerr << "Backtracking failed." << std::endl;
00795               std::cout << "Backtracking failed." << std::endl;
00796           
00797               // 1.) Quit, exit program.
00798               //libmesh_error();
00799 
00800               // 2.) Continue with last newton_stepfactor
00801               if (newton_step<3)
00802                 {
00803                   solution->add(newton_stepfactor, *delta_u);
00804                   solution->close();
00805                   *continuation_parameter += newton_stepfactor*delta_lambda;
00806                   if (!quiet)
00807                     std::cout << "Backtracking could not reduce residual ... continuing anyway!" << std::endl;
00808                 }
00809               
00810               // 3.) Break out of Newton iteration loop with newton_converged = false,
00811               //     reduce the arclength stepsize, and try again.
00812               else
00813                 {
00814                   break; // out of Newton iteration loop, with newton_converged=false
00815                 }
00816             }
00817 
00818           // Another type of convergence check: suppose the residual has not been reduced
00819           // from its initial value after half of the allowed Newton steps have occurred.
00820           // In our experience, this typically means that it isn't going to converge and
00821           // we could probably save time by dropping out of the Newton iteration loop and
00822           // trying a smaller arcstep.
00823           if (this->newton_progress_check)
00824             {
00825               if ((nonlinear_residual_afterstep > nonlinear_residual_firststep) &&
00826                   (newton_step+1 > static_cast<unsigned int>(0.5*newton_solver->max_nonlinear_iterations)))
00827                 {
00828                   std::cout << "Progress check failed: the current residual: "
00829                             << nonlinear_residual_afterstep
00830                             << ", is\n"
00831                             << "larger than the initial residual, and half of the allowed\n"
00832                             << "number of Newton iterations have elapsed.\n"
00833                             << "Exiting Newton iterations with converged==false." << std::endl;
00834               
00835                   break; // out of Newton iteration loop, newton_converged = false
00836                 }
00837             }
00838           
00839           // Safety check: Check the current continuation parameter against user-provided min-allowable parameter value
00840           if (*continuation_parameter < min_continuation_parameter)
00841             {
00842               std::cout << "Continuation parameter fell below min-allowable value." << std::endl;
00843               // libmesh_error();
00844               break; // out of Newton iteration loop, newton_converged = false
00845             }
00846           
00847           // Safety check: Check the current continuation parameter against user-provided max-allowable parameter value
00848           if ( (max_continuation_parameter != 0.0) &&
00849                (*continuation_parameter > max_continuation_parameter) )
00850             {
00851               std::cout << "Current continuation parameter value: "
00852                         << *continuation_parameter
00853                         << " exceeded max-allowable value."
00854                         << std::endl;
00855               // libmesh_error();
00856               break; // out of Newton iteration loop, newton_converged = false
00857             }
00858           
00859       
00860           // Check the convergence of the parameter and the solution.  If they are small
00861           // enough, we can break out of the Newton iteration loop.  
00862           const Real norm_delta_u = delta_u->l2_norm();
00863           const Real norm_u = solution->l2_norm();
00864           std::cout << "  delta_lambda                   = " << delta_lambda << std::endl;
00865           std::cout << "  newton_stepfactor*delta_lambda = " << newton_stepfactor*delta_lambda << std::endl;
00866           std::cout << "  lambda_current                 = " << *continuation_parameter << std::endl;
00867           std::cout << "  ||delta_u||                    = " << norm_delta_u << std::endl;
00868           std::cout << "  ||delta_u||/||u||              = " << norm_delta_u / norm_u << std::endl;
00869       
00870 
00871           // Evaluate the residual at the current Newton iterate.  We don't want to detect
00872           // convergence due to a small Newton step when the residual is still not small.
00873           rhs_mode = Residual;
00874           assembly(true,   // Residual
00875                    false); // Jacobian
00876           rhs->close();
00877           const Real norm_residual = rhs->l2_norm();
00878           std::cout << "  ||R||_{L2} = " << norm_residual << std::endl;
00879           std::cout << "  ||R||_{L2}/||u|| = " << norm_residual / norm_u << std::endl;
00880       
00881       
00882           // FIXME: The norm_delta_u tolerance (at least) should be relative.
00883           // It doesn't make sense to converge a solution whose size is ~ 10^5 to
00884           // a tolerance of 1.e-6.  Oh, and we should also probably check the
00885           // (relative) size of the residual as well, instead of just the step.
00886           if ((std::abs(delta_lambda) < continuation_parameter_tolerance) &&
00887               //(norm_delta_u       < solution_tolerance)               && // This is a *very* strict criterion we can probably skip
00888               (norm_residual      < solution_tolerance))
00889             {
00890               if (!quiet)
00891                 std::cout << "Newton iterations converged!" << std::endl;
00892 
00893               newton_converged = true;
00894               break; // out of Newton iterations
00895             }
00896         } // end nonlinear loop
00897 
00898       if (!newton_converged)
00899         {
00900           std::cout << "Newton iterations of augmented system did not converge!" << std::endl;
00901           
00902           // Reduce ds_current, recompute the solution and parameter, and continue to next
00903           // arcstep, if there is one.
00904           ds_current *= 0.5;
00905 
00906           // Go back to previous solution and parameter value.
00907           *solution = *previous_u;
00908           *continuation_parameter = old_continuation_parameter;
00909 
00910           // Compute new predictor with smaller ds
00911           apply_predictor();
00912         }
00913       else
00914         {
00915           // Set step convergence and break out
00916           arcstep_converged=true;
00917           break; // out of arclength reduction loop
00918         }
00919       
00920     } // end loop over arclength reductions
00921 
00922   // Check for convergence of the whole arcstep.  If not converged at this
00923   // point, we have no choice but to quit.
00924   if (!arcstep_converged)
00925     {
00926       std::cout << "Arcstep failed to converge after max number of reductions! Exiting..." << std::endl;
00927       libmesh_error();
00928     }
00929   
00930   // Print converged solution control parameter and max value.
00931   std::cout << "lambda_current=" << *continuation_parameter << std::endl;
00932   //std::cout << "u_max=" << solution->max() << std::endl;
00933 
00934   // Reset old stream precision and flags.
00935   std::cout.precision(old_precision);
00936   std::cout.unsetf(std::ios_base::scientific);
00937 
00938   // Note: we don't want to go on to the next guess yet, since the user may
00939   // want to post-process this data.  It's up to the user to call advance_arcstep()
00940   // when they are ready to go on.
00941 }

Number System::current_solution ( const unsigned int  global_dof_number  )  const [inherited]

Returns:
the current solution for the specified global DOF.

Definition at line 117 of file system.C.

References System::current_local_solution, and System::n_dofs().

Referenced by ExactSolution::_compute_error(), HPCoarsenTest::add_projection(), JumpErrorEstimator::estimate_error(), ExactErrorEstimator::estimate_error(), FEMSystem::eulerian_residual(), PatchRecoveryErrorEstimator::EstimateError::operator()(), FEMContext::reinit(), HPCoarsenTest::select_refinement(), VTKIO::solution_to_vtk(), EnsightIO::write_scalar_ascii(), and EnsightIO::write_vector_ascii().

00118 {
00119   // Check the sizes
00120   libmesh_assert (global_dof_number < _dof_map->n_dofs());
00121   libmesh_assert (global_dof_number < current_local_solution->size());
00122    
00123   return (*current_local_solution)(global_dof_number);
00124 }

void System::deactivate (  )  [inline, inherited]

Deactivates the system. Only active systems are solved.

Definition at line 1382 of file system.h.

References System::_active.

01383 {
01384   _active = false;
01385 }

virtual bool DifferentiableSystem::element_constraint ( bool  request_jacobian,
DiffContext  
) [inline, virtual, inherited]

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.

Definition at line 157 of file diff_system.h.

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

00158                                                   {
00159     return request_jacobian;
00160   }

virtual void DifferentiableSystem::element_postprocess ( DiffContext  )  [inline, virtual, inherited]

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

Definition at line 225 of file diff_system.h.

Referenced by FEMSystem::postprocess().

00225 {}

virtual void DifferentiableSystem::element_qoi ( DiffContext ,
const QoISet  
) [inline, virtual, inherited]

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 240 of file diff_system.h.

Referenced by FEMSystem::assemble_qoi().

00242     {}

virtual void DifferentiableSystem::element_qoi_derivative ( DiffContext ,
const QoISet  
) [inline, virtual, inherited]

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 252 of file diff_system.h.

Referenced by FEMSystem::assemble_qoi_derivative().

00254     {}

virtual bool DifferentiableSystem::element_time_derivative ( bool  request_jacobian,
DiffContext  
) [inline, virtual, inherited]

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.

Definition at line 140 of file diff_system.h.

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

00141                                                        {
00142     return request_jacobian;
00143   }

bool FEMSystem::eulerian_residual ( bool  request_jacobian,
DiffContext context 
) [virtual, inherited]

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 DifferentiableSystem.

Definition at line 864 of file fem_system.C.

References FEMSystem::_mesh_sys, FEMSystem::_mesh_x_var, FEMSystem::_mesh_y_var, FEMSystem::_mesh_z_var, DifferentiableSystem::_time_evolving, System::current_solution(), DifferentiableSystem::deltat, DiffContext::dof_indices_var, DiffContext::elem_subjacobians, DiffContext::elem_subresiduals, FEMContext::element_fe_var, FEMContext::element_qrule, AutoPtr< Tp >::get(), FEMContext::interior_gradient(), libMesh::invalid_uint, libmesh_real(), QBase::n_points(), System::n_vars(), System::number(), UnsteadySolver::old_nonlinear_solution(), and DifferentiableSystem::time_solver.

00866 {
00867   // Only calculate a mesh movement residual if it's necessary
00868   if (_mesh_sys == libMesh::invalid_uint)
00869     return request_jacobian;
00870 
00871   FEMContext &context = libmesh_cast_ref<FEMContext&>(c);
00872 
00873   // This function only supports fully coupled mesh motion for now
00874   libmesh_assert(_mesh_sys == this->number());
00875 
00876   unsigned int n_qpoints = context.element_qrule->n_points();
00877 
00878   const unsigned int n_x_dofs = (_mesh_x_var == libMesh::invalid_uint) ?
00879                                 0 : context.dof_indices_var[_mesh_x_var].size();
00880   const unsigned int n_y_dofs = (_mesh_y_var == libMesh::invalid_uint) ?
00881                                 0 : context.dof_indices_var[_mesh_y_var].size();
00882   const unsigned int n_z_dofs = (_mesh_z_var == libMesh::invalid_uint) ?
00883                                 0 : context.dof_indices_var[_mesh_z_var].size();
00884 
00885   const unsigned int mesh_xyz_var = n_x_dofs ? _mesh_x_var :
00886                                    (n_y_dofs ? _mesh_y_var :
00887                                    (n_z_dofs ? _mesh_z_var :
00888                                     libMesh::invalid_uint));
00889 
00890   // If we're our own _mesh_sys, we'd better be in charge of
00891   // at least one coordinate, and we'd better have the same
00892   // FE type for all coordinates we are in charge of
00893   libmesh_assert(mesh_xyz_var != libMesh::invalid_uint);
00894   libmesh_assert(!n_x_dofs || context.element_fe_var[_mesh_x_var] ==
00895                               context.element_fe_var[mesh_xyz_var]);
00896   libmesh_assert(!n_y_dofs || context.element_fe_var[_mesh_y_var] ==
00897                               context.element_fe_var[mesh_xyz_var]);
00898   libmesh_assert(!n_z_dofs || context.element_fe_var[_mesh_z_var] ==
00899                               context.element_fe_var[mesh_xyz_var]);
00900 
00901   const std::vector<std::vector<Real> >     &psi =
00902     context.element_fe_var[mesh_xyz_var]->get_phi();
00903 
00904   for (unsigned int var = 0; var != this->n_vars(); ++var)
00905     {
00906       // Mesh motion only affects time-evolving variables
00907       if (!_time_evolving[var])
00908         continue;
00909 
00910       // The mesh coordinate variables themselves are Lagrangian,
00911       // not Eulerian, and no convective term is desired.
00912       if (_mesh_sys == this->number() &&
00913           (var == _mesh_x_var ||
00914            var == _mesh_y_var ||
00915            var == _mesh_z_var))
00916         continue;
00917 
00918       // Some of this code currently relies on the assumption that
00919       // we can pull mesh coordinate data from our own system
00920       if (_mesh_sys != this->number())
00921         libmesh_not_implemented();
00922 
00923       // This residual should only be called by unsteady solvers:
00924       // if the mesh is steady, there's no mesh convection term!
00925       UnsteadySolver *unsteady =
00926         dynamic_cast<UnsteadySolver *>(this->time_solver.get());
00927       if (!unsteady)
00928         return request_jacobian;
00929 
00930       const std::vector<Real> &JxW = 
00931         context.element_fe_var[var]->get_JxW();
00932 
00933       const std::vector<std::vector<Real> >     &phi =
00934         context.element_fe_var[var]->get_phi();
00935 
00936       const std::vector<std::vector<RealGradient> > &dphi =
00937         context.element_fe_var[var]->get_dphi();
00938 
00939       const unsigned int n_u_dofs = context.dof_indices_var[var].size();
00940 
00941       DenseSubVector<Number> &Fu = *context.elem_subresiduals[var];
00942       DenseSubMatrix<Number> &Kuu = *context.elem_subjacobians[var][var];
00943 
00944       DenseSubMatrix<Number> *Kux = n_x_dofs ?
00945         context.elem_subjacobians[var][_mesh_x_var] : NULL;
00946       DenseSubMatrix<Number> *Kuy = n_y_dofs ?
00947         context.elem_subjacobians[var][_mesh_y_var] : NULL;
00948       DenseSubMatrix<Number> *Kuz = n_z_dofs ?
00949         context.elem_subjacobians[var][_mesh_z_var] : NULL;
00950 
00951       std::vector<Real> delta_x(n_x_dofs, 0.);
00952       std::vector<Real> delta_y(n_y_dofs, 0.);
00953       std::vector<Real> delta_z(n_z_dofs, 0.);
00954 
00955       for (unsigned int i = 0; i != n_x_dofs; ++i)
00956         {
00957           unsigned int j = context.dof_indices_var[_mesh_x_var][i];
00958           delta_x[i] = libmesh_real(this->current_solution(j)) -
00959                        libmesh_real(unsteady->old_nonlinear_solution(j));
00960         }
00961 
00962       for (unsigned int i = 0; i != n_y_dofs; ++i)
00963         {
00964           unsigned int j = context.dof_indices_var[_mesh_y_var][i];
00965           delta_y[i] = libmesh_real(this->current_solution(j)) -
00966                        libmesh_real(unsteady->old_nonlinear_solution(j));
00967         }
00968 
00969       for (unsigned int i = 0; i != n_z_dofs; ++i)
00970         {
00971           unsigned int j = context.dof_indices_var[_mesh_z_var][i];
00972           delta_z[i] = libmesh_real(this->current_solution(j)) -
00973                        libmesh_real(unsteady->old_nonlinear_solution(j));
00974         }
00975 
00976       for (unsigned int qp = 0; qp != n_qpoints; ++qp)
00977         {
00978           Gradient grad_u = context.interior_gradient(var, qp);
00979           RealGradient convection(0.);
00980 
00981           for (unsigned int i = 0; i != n_x_dofs; ++i)
00982             convection(0) += delta_x[i] * psi[i][qp];
00983           for (unsigned int i = 0; i != n_y_dofs; ++i)
00984             convection(1) += delta_y[i] * psi[i][qp];
00985           for (unsigned int i = 0; i != n_z_dofs; ++i)
00986             convection(2) += delta_z[i] * psi[i][qp];
00987 
00988           for (unsigned int i = 0; i != n_u_dofs; ++i)
00989             {
00990               Number JxWxPhiI = JxW[qp] * phi[i][qp];
00991               Fu(i) += (convection * grad_u) * JxWxPhiI;
00992               if (request_jacobian)
00993                 {
00994                   Number JxWxPhiI = JxW[qp] * phi[i][qp];
00995                   for (unsigned int j = 0; j != n_u_dofs; ++j)
00996                     Kuu(i,j) += JxWxPhiI * (convection * dphi[j][qp]);
00997 
00998                   Number JxWxPhiIoverDT = JxWxPhiI/this->deltat;
00999 
01000                   Number JxWxPhiIxDUDXoverDT = JxWxPhiIoverDT * grad_u(0);
01001                   for (unsigned int j = 0; j != n_x_dofs; ++j)
01002                     (*Kux)(i,j) += JxWxPhiIxDUDXoverDT * psi[j][qp];
01003 
01004                   Number JxWxPhiIxDUDYoverDT = JxWxPhiIoverDT * grad_u(1);
01005                   for (unsigned int j = 0; j != n_y_dofs; ++j)
01006                     (*Kuy)(i,j) += JxWxPhiIxDUDYoverDT * psi[j][qp];
01007 
01008                   Number JxWxPhiIxDUDZoverDT = JxWxPhiIoverDT * grad_u(2);
01009                   for (unsigned int j = 0; j != n_z_dofs; ++j)
01010                     (*Kuz)(i,j) += JxWxPhiIxDUDZoverDT * psi[j][qp];
01011                 }
01012             }
01013         }
01014     }
01015 
01016   return request_jacobian;
01017 }

void ImplicitSystem::forward_qoi_parameter_sensitivity ( const QoISet qoi_indices,
const ParameterVector parameters,
SensitivityData sensitivities 
) [virtual, inherited]

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 System.

Definition at line 749 of file implicit_system.C.

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

00752 {
00753   const unsigned int Np = parameters.size();
00754   const unsigned int Nq = qoi.size();
00755 
00756   // We currently get partial derivatives via central differencing
00757   const Real delta_p = TOLERANCE;
00758 
00759   // An introduction to the problem:
00760   //
00761   // Residual R(u(p),p) = 0
00762   // partial R / partial u = J = system matrix
00763   //
00764   // This implies that:
00765   // d/dp(R) = 0
00766   // (partial R / partial p) + 
00767   // (partial R / partial u) * (partial u / partial p) = 0
00768 
00769   // We first solve for (partial u / partial p) for each parameter:
00770   // J * (partial u / partial p) = - (partial R / partial p)
00771 
00772   this->sensitivity_solve(parameters);
00773 
00774   // Get ready to fill in senstivities:
00775   sensitivities.allocate_data(qoi_indices, *this, parameters);
00776 
00777   // We use the identity:
00778   // dq/dp = (partial q / partial p) + (partial q / partial u) *
00779   //         (partial u / partial p)
00780  
00781   // We get (partial q / partial u) from the user
00782   this->assemble_qoi_derivative(qoi_indices);
00783 
00784   for (unsigned int j=0; j != Np; ++j)
00785     {
00786       // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
00787 
00788       Number old_parameter = *parameters[j];
00789 
00790       *parameters[j] = old_parameter - delta_p;
00791       this->assemble_qoi();
00792       std::vector<Number> qoi_minus = this->qoi;
00793 
00794       *parameters[j] = old_parameter + delta_p;
00795       this->assemble_qoi();
00796       std::vector<Number>& qoi_plus = this->qoi;
00797 
00798       std::vector<Number> partialq_partialp(Nq, 0);
00799       for (unsigned int i=0; i != Nq; ++i)
00800         if (qoi_indices.has_index(i))
00801           partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);
00802 
00803       // Don't leave the parameter changed
00804       *parameters[j] = old_parameter;
00805 
00806       for (unsigned int i=0; i != Nq; ++i)
00807         if (qoi_indices.has_index(i))
00808           sensitivities[i][j] = partialq_partialp[i] +
00809             this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(i));
00810     }
00811 
00812   // All parameters have been reset.
00813   // We didn't cache the original rhs or matrix for memory reasons,
00814   // but we can restore them to a state consistent solution -
00815   // principle of least surprise.
00816   this->assembly(true, true);
00817   this->rhs->close();
00818   this->matrix->close();
00819   this->assemble_qoi(qoi_indices);
00820 }

const NumericVector< Number > & 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 844 of file system.C.

References System::get_vector().

00845 {
00846   OStringStream adjoint_rhs_name;
00847   adjoint_rhs_name << "adjoint_rhs" << i;
00848 
00849   return this->get_vector(adjoint_rhs_name.str());
00850 }

NumericVector< Number > & 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 834 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::adjoint_solve(), FEMSystem::assemble_qoi_derivative(), and ImplicitSystem::weighted_sensitivity_adjoint_solve().

00835 {
00836   OStringStream adjoint_rhs_name;
00837   adjoint_rhs_name << "adjoint_rhs" << i;
00838 
00839   return this->get_vector(adjoint_rhs_name.str());
00840 }

const NumericVector< Number > & 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 784 of file system.C.

References System::get_vector().

00785 {
00786   OStringStream adjoint_name;
00787   adjoint_name << "adjoint_solution" << i;
00788 
00789   return this->get_vector(adjoint_name.str());
00790 }

NumericVector< Number > & 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 774 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::adjoint_solve(), AdjointResidualErrorEstimator::estimate_error(), and ImplicitSystem::weighted_sensitivity_adjoint_solve().

00775 {
00776   OStringStream adjoint_name;
00777   adjoint_name << "adjoint_solution" << i;
00778 
00779   return this->get_vector(adjoint_name.str());
00780 }

DofMap & System::get_dof_map (  )  [inline, inherited]

Returns:
a writeable reference to this system's _dof_map.

Definition at line 1358 of file system.h.

References System::_dof_map.

01359 {
01360   return *_dof_map;
01361 }

EquationSystems& System::get_equation_systems (  )  [inline, inherited]

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

Definition at line 383 of file system.h.

References System::_equation_systems.

00383 { return _equation_systems; }

std::string ReferenceCounter::get_info (  )  [static, inherited]

Gets a string containing the reference information.

Definition at line 45 of file reference_counter.C.

References ReferenceCounter::_counts, and QuadratureRules::name().

Referenced by ReferenceCounter::print_info().

00046 {
00047 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
00048 
00049   std::ostringstream out;
00050   
00051   out << '\n'
00052       << " ---------------------------------------------------------------------------- \n"
00053       << "| Reference count information                                                |\n"
00054       << " ---------------------------------------------------------------------------- \n";
00055   
00056   for (Counts::iterator it = _counts.begin();
00057        it != _counts.end(); ++it)
00058     {
00059       const std::string name(it->first);
00060       const unsigned int creations    = it->second.first;
00061       const unsigned int destructions = it->second.second;
00062 
00063       out << "| " << name << " reference count information:\n"
00064           << "|  Creations:    " << creations    << '\n'
00065           << "|  Destructions: " << destructions << '\n';
00066     }
00067   
00068   out << " ---------------------------------------------------------------------------- \n";
00069 
00070   return out.str();
00071 
00072 #else
00073 
00074   return "";
00075   
00076 #endif
00077 }

std::string System::get_info (  )  const [inherited]

Returns:
a string containing information about the system.

Definition at line 1233 of file system.C.

References Utility::enum_to_string< FEFamily >(), Utility::enum_to_string< InfMapType >(), Utility::enum_to_string< Order >(), FEType::family, System::get_dof_map(), FEType::inf_map, System::n_constrained_dofs(), System::n_dofs(), System::n_local_dofs(), System::n_vars(), System::n_vectors(), System::name(), FEType::order, FEType::radial_family, FEType::radial_order, System::system_type(), System::variable_name(), and DofMap::variable_type().

01234 {
01235   std::ostringstream out;
01236 
01237   
01238   const std::string& sys_name = this->name();
01239       
01240   out << "   System \"" << sys_name << "\"\n"
01241       << "    Type \""  << this->system_type() << "\"\n"
01242       << "    Variables=";
01243   
01244   for (unsigned int vn=0; vn<this->n_vars(); vn++)
01245       out << "\"" << this->variable_name(vn) << "\" ";
01246      
01247   out << '\n';
01248 
01249   out << "    Finite Element Types=";
01250 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
01251   for (unsigned int vn=0; vn<this->n_vars(); vn++)
01252     out << "\""
01253         << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_type(vn).family)
01254         << "\" ";  
01255 #else
01256   for (unsigned int vn=0; vn<this->n_vars(); vn++)
01257     {
01258       out << "\""
01259           << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_type(vn).family)
01260           << "\", \""
01261           << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_type(vn).radial_family)
01262           << "\" ";
01263     }
01264 
01265   out << '\n' << "    Infinite Element Mapping=";
01266   for (unsigned int vn=0; vn<this->n_vars(); vn++)
01267     out << "\""
01268         << Utility::enum_to_string<InfMapType>(this->get_dof_map().variable_type(vn).inf_map)
01269         << "\" ";
01270 #endif      
01271 
01272   out << '\n';
01273       
01274   out << "    Approximation Orders=";
01275   for (unsigned int vn=0; vn<this->n_vars(); vn++)
01276     {
01277 #ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
01278       out << "\""
01279           << Utility::enum_to_string<Order>(this->get_dof_map().variable_type(vn).order)
01280           << "\" ";
01281 #else
01282       out << "\""
01283           << Utility::enum_to_string<Order>(this->get_dof_map().variable_type(vn).order)
01284           << "\", \""
01285           << Utility::enum_to_string<Order>(this->get_dof_map().variable_type(vn).radial_order)
01286           << "\" ";
01287 #endif
01288     }
01289 
01290   out << '\n';
01291       
01292   out << "    n_dofs()="             << this->n_dofs()             << '\n';
01293   out << "    n_local_dofs()="       << this->n_local_dofs()       << '\n';
01294 #ifdef LIBMESH_ENABLE_AMR
01295   out << "    n_constrained_dofs()=" << this->n_constrained_dofs() << '\n';
01296 #endif
01297 
01298   out << "    " << "n_vectors()="  << this->n_vectors()  << '\n';
01299 //   out << "    " << "n_additional_matrices()=" << this->n_additional_matrices() << '\n';
01300   
01301   return out.str();
01302 }

std::pair< unsigned int, Real > DifferentiableSystem::get_linear_solve_parameters (  )  const [virtual, inherited]

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 ImplicitSystem.

Definition at line 123 of file diff_system.C.

References DifferentiableSystem::time_solver.

00124 {
00125   return std::make_pair(this->time_solver->diff_solver()->max_linear_iterations,
00126                         this->time_solver->diff_solver()->relative_residual_tolerance);
00127 }

LinearSolver< Number > * DifferentiableSystem::get_linear_solver (  )  const [virtual, inherited]

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

Reimplemented from ImplicitSystem.

Definition at line 116 of file diff_system.C.

References AutoPtr< Tp >::get(), and DifferentiableSystem::time_solver.

00117 {
00118   return this->time_solver->linear_solver().get();
00119 }

SparseMatrix< Number > & 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 245 of file implicit_system.C.

References ImplicitSystem::_matrices.

00246 {
00247   // Make sure the matrix exists
00248   matrices_iterator pos = _matrices.find (mat_name);
00249   
00250   if (pos == _matrices.end())
00251     {
00252       std::cerr << "ERROR: matrix "
00253                 << mat_name
00254                 << " does not exist in this system!"
00255                 << std::endl;      
00256       libmesh_error();
00257     }
00258   
00259   return *(pos->second);
00260 }

const SparseMatrix< Number > & 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 226 of file implicit_system.C.

References ImplicitSystem::_matrices.

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

00227 {
00228   // Make sure the matrix exists
00229   const_matrices_iterator pos = _matrices.find (mat_name);
00230   
00231   if (pos == _matrices.end())
00232     {
00233       std::cerr << "ERROR: matrix "
00234                 << mat_name
00235                 << " does not exist in this system!"
00236                 << std::endl;      
00237       libmesh_error();
00238     }
00239   
00240   return *(pos->second);
00241 }

MeshBase & System::get_mesh (  )  [inline, inherited]

Returns:
a reference to this systems's _mesh.

Definition at line 1342 of file system.h.

References System::_mesh.

01343 {
01344   return _mesh;
01345 }

const NumericVector< Number > & 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 874 of file system.C.

References System::get_vector().

00875 {
00876   OStringStream sensitivity_rhs_name;
00877   sensitivity_rhs_name << "sensitivity_rhs" << i;
00878 
00879   return this->get_vector(sensitivity_rhs_name.str());
00880 }

NumericVector< Number > & 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 864 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::sensitivity_solve().

00865 {
00866   OStringStream sensitivity_rhs_name;
00867   sensitivity_rhs_name << "sensitivity_rhs" << i;
00868 
00869   return this->get_vector(sensitivity_rhs_name.str());
00870 }

const NumericVector< Number > & 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 733 of file system.C.

References System::get_vector().

00734 {
00735   OStringStream sensitivity_name;
00736   sensitivity_name << "sensitivity_solution" << i;
00737 
00738   return this->get_vector(sensitivity_name.str());
00739 }

NumericVector< Number > & 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 723 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::sensitivity_solve().

00724 {
00725   OStringStream sensitivity_name;
00726   sensitivity_name << "sensitivity_solution" << i;
00727 
00728   return this->get_vector(sensitivity_name.str());
00729 }

NumericVector< Number > & 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 682 of file system.C.

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

00683 {
00684   vectors_iterator v = vectors_begin();
00685   vectors_iterator v_end = vectors_end();
00686   unsigned int num = 0;
00687   while((num<vec_num) && (v!=v_end))
00688     {
00689       num++;
00690       ++v;
00691     }
00692   libmesh_assert(v!=v_end);
00693   return *(v->second);
00694 }

const NumericVector< Number > & 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 666 of file system.C.

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

00667 {
00668   const_vectors_iterator v = vectors_begin();
00669   const_vectors_iterator v_end = vectors_end();
00670   unsigned int num = 0;
00671   while((num<vec_num) && (v!=v_end))
00672     {
00673       num++;
00674       ++v;
00675     }
00676   libmesh_assert(v!=v_end);
00677   return *(v->second);
00678 }

NumericVector< Number > & 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 647 of file system.C.

References System::_vectors.

00648 {
00649   // Make sure the vector exists
00650   vectors_iterator pos = _vectors.find(vec_name);
00651   
00652   if (pos == _vectors.end())
00653     {
00654       std::cerr << "ERROR: vector "
00655                 << vec_name
00656                 << " does not exist in this system!"
00657                 << std::endl;      
00658       libmesh_error();
00659     }
00660   
00661   return *(pos->second);
00662 }

const NumericVector< Number > & 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 628 of file system.C.

References System::_vectors.

Referenced by System::add_weighted_sensitivity_adjoint_solution(), UnsteadySolver::advance_timestep(), AdaptiveTimeSolver::advance_timestep(), System::compare(), UnsteadySolver::du(), System::get_adjoint_rhs(), System::get_adjoint_solution(), System::get_sensitivity_rhs(), System::get_sensitivity_solution(), System::get_weighted_sensitivity_adjoint_solution(), System::get_weighted_sensitivity_solution(), NewmarkSystem::initial_conditions(), UnsteadySolver::solve(), TwostepTimeSolver::solve(), FrequencySystem::solve(), NewmarkSystem::update_rhs(), and NewmarkSystem::update_u_v_a().

00629 {
00630   // Make sure the vector exists
00631   const_vectors_iterator pos = _vectors.find(vec_name);
00632   
00633   if (pos == _vectors.end())
00634     {
00635       std::cerr << "ERROR: vector "
00636                 << vec_name
00637                 << " does not exist in this system!"
00638                 << std::endl;      
00639       libmesh_error();
00640     }
00641   
00642   return *(pos->second);
00643 }

const NumericVector< Number > & 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 814 of file system.C.

References System::get_vector().

00815 {
00816   OStringStream adjoint_name;
00817   adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
00818 
00819   return this->get_vector(adjoint_name.str());
00820 }

NumericVector< Number > & 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 804 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::weighted_sensitivity_adjoint_solve().

00805 {
00806   OStringStream adjoint_name;
00807   adjoint_name << "weighted_sensitivity_adjoint_solution" << i;
00808 
00809   return this->get_vector(adjoint_name.str());
00810 }

const NumericVector< Number > & System::get_weighted_sensitivity_solution (  )  const [inherited]

Returns:
a reference to the solution of the last weighted sensitivity solve

Definition at line 757 of file system.C.

References System::get_vector().

00758 {
00759   return this->get_vector("weighted_sensitivity_solution");
00760 }

NumericVector< Number > & System::get_weighted_sensitivity_solution (  )  [inherited]

Returns:
a reference to the solution of the last weighted sensitivity solve

Definition at line 750 of file system.C.

References System::get_vector().

Referenced by ImplicitSystem::weighted_sensitivity_solve().

00751 {
00752   return this->get_vector("weighted_sensitivity_solution");
00753 }

bool System::has_variable ( const std::string &  var  )  const [inherited]

Returns:
true if a variable named var exists in this System

Definition at line 935 of file system.C.

References System::_variable_numbers.

Referenced by GMVIO::copy_nodal_solution().

00936 {
00937   return _variable_numbers.count(var);
00938 }

bool ImplicitSystem::have_matrix ( const std::string &  mat_name  )  const [inline, inherited]

Returns:
true if this System has a matrix associated with the given name, false otherwise.

Definition at line 357 of file implicit_system.h.

References ImplicitSystem::_matrices.

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

00358 {
00359   return (_matrices.count(mat_name));
00360 }

bool System::have_vector ( const std::string &  vec_name  )  const [inline, inherited]

Returns:
true if this System has a vector associated with the given name, false otherwise.

Definition at line 1444 of file system.h.

References System::_vectors.

Referenced by System::add_vector().

01445 {
01446   return (_vectors.count(vec_name));
01447 }

void ReferenceCounter::increment_constructor_count ( const std::string &  name  )  [inline, protected, inherited]

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

Definition at line 149 of file reference_counter.h.

References ReferenceCounter::_counts, and Threads::spin_mtx.

Referenced by ReferenceCountedObject< SparseMatrix< T > >::ReferenceCountedObject().

00150 {
00151   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
00152   std::pair<unsigned int, unsigned int>& p = _counts[name];
00153 
00154   p.first++;
00155 }

void ReferenceCounter::increment_destructor_count ( const std::string &  name  )  [inline, protected, inherited]

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

Definition at line 167 of file reference_counter.h.

References ReferenceCounter::_counts, and Threads::spin_mtx.

Referenced by ReferenceCountedObject< SparseMatrix< T > >::~ReferenceCountedObject().

00168 {
00169   Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
00170   std::pair<unsigned int, unsigned int>& p = _counts[name];
00171 
00172   p.second++;
00173 }

void System::init (  )  [inherited]

Initializes degrees of freedom on the current mesh. Sets the

Definition at line 158 of file system.C.

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

00159 { 
00160   // First initialize any required data
00161   this->init_data();
00162   
00163   //If no variables have been added to this system
00164   //don't do anything
00165   if(!this->n_vars())
00166     return;
00167  
00168   // Then call the user-provided intialization function
00169   this->user_initialization();
00170 }

void FEMSystem::init_context ( DiffContext c  )  [virtual, inherited]

Reimplemented from DifferentiableSystem.

Definition at line 749 of file fem_system.C.

References DifferentiableSystem::_time_evolving, FEMContext::element_fe_var, DifferentiableSystem::init_context(), and System::n_vars().

Referenced by FEMSystem::mesh_position_get().

00750 {
00751   Parent::init_context(c);
00752 
00753   FEMContext &context = libmesh_cast_ref<FEMContext&>(c);
00754 
00755   // Make sure we're prepared to do mass integration
00756   for (unsigned int var = 0; var != this->n_vars(); ++var)
00757     if (_time_evolving[var])
00758       {
00759         context.element_fe_var[var]->get_JxW();
00760         context.element_fe_var[var]->get_phi();
00761       }
00762 }

void ContinuationSystem::init_data (  )  [protected, virtual]

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

Reimplemented from FEMSystem.

Definition at line 84 of file continuation_system.C.

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

00085 {
00086   // Add a vector which stores the tangent "du/ds" to the system and save its pointer.
00087   du_ds = &(add_vector("du_ds"));
00088 
00089   // Add a vector which stores the tangent "du/ds" to the system and save its pointer.
00090   previous_du_ds = &(add_vector("previous_du_ds"));
00091 
00092   // Add a vector to keep track of the previous nonlinear solution
00093   // at the old value of lambda.
00094   previous_u = &(add_vector("previous_u"));
00095 
00096   // Add a vector to keep track of the temporary solution "y" of Ay=G_{\lambda}.
00097   y = &(add_vector("y"));
00098 
00099   // Add a vector to keep track of the "old value" of "y" which is the solution of Ay=G_{\lambda}.
00100   y_old = &(add_vector("y_old"));
00101 
00102   // Add a vector to keep track of the temporary solution "z" of Az=-G.
00103   z = &(add_vector("z"));
00104 
00105   // Add a vector to keep track of the Newton update during the constrained PDE solves.
00106   delta_u = &(add_vector("delta_u"));
00107 
00108   // Call the Parent's initialization routine.
00109   Parent::init_data();
00110 }

void ImplicitSystem::init_matrices (  )  [protected, virtual, inherited]

Initializes the matrices associated with this system.

Definition at line 99 of file implicit_system.C.

References ImplicitSystem::_can_add_matrices, ImplicitSystem::_matrices, DofMap::attach_matrix(), DofMap::compute_sparsity(), System::get_dof_map(), System::get_mesh(), and ImplicitSystem::matrix.

Referenced by ImplicitSystem::init_data(), and ImplicitSystem::reinit().

00100 {
00101   libmesh_assert (matrix != NULL);
00102 
00103   // Check for quick return in case the system matrix
00104   // (and by extension all the matrices) has already
00105   // been initialized
00106   if (matrix->initialized())
00107     return;
00108 
00109   // Get a reference to the DofMap
00110   DofMap& dof_map = this->get_dof_map();
00111   
00112   // no chance to add other matrices
00113   _can_add_matrices = false;
00114   
00115   // Tell the matrices about the dof map, and vice versa
00116   for (matrices_iterator pos = _matrices.begin();
00117        pos != _matrices.end(); ++pos)
00118     {
00119       libmesh_assert (!pos->second->initialized());
00120       dof_map.attach_matrix (*(pos->second));
00121     }
00122   
00123   // Compute the sparsity pattern for the current
00124   // mesh and DOF distribution.  This also updates
00125   // additional matrices, \p DofMap now knows them
00126   dof_map.compute_sparsity (this->get_mesh());
00127   
00128   // Initialize matrices
00129   for (matrices_iterator pos = _matrices.begin(); 
00130        pos != _matrices.end(); ++pos)
00131     pos->second->init ();
00132   
00133   // Set the additional matrices to 0.
00134   for (matrices_iterator pos = _matrices.begin(); 
00135        pos != _matrices.end(); ++pos)
00136     pos->second->zero ();
00137 }

void 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 125 of file continuation_system.C.

References continuation_parameter, dlambda_ds, ds_current, du_ds, old_continuation_parameter, previous_u, quiet, set_Theta(), System::solution, solve_tangent(), tangent_initialized, Theta, Theta_LOCA, update_solution(), and y.

Referenced by continuation_solve().

00126 {
00127   // Be sure the tangent was not already initialized.
00128   libmesh_assert (!tangent_initialized);
00129   
00130   // Compute delta_s_zero, the initial arclength travelled during the
00131   // first step.  Here we assume that previous_u and lambda_old store
00132   // the previous solution and control parameter.  You may need to
00133   // read in an old solution (or solve the non-continuation system)
00134   // first and call save_current_solution() before getting here.
00135   
00136   // 1.) Compute delta_s_zero as ||u|| - ||u_old|| + ...
00137   // Compute norms of the current and previous solutions
00138 //   Real norm_u          = solution->l2_norm();
00139 //   Real norm_previous_u = previous_u->l2_norm();
00140   
00141 //   if (!quiet)
00142 //     {
00143 //       std::cout << "norm_u=" << norm_u << std::endl;
00144 //       std::cout << "norm_previous_u=" << norm_previous_u << std::endl;
00145 //     }
00146   
00147 //   if (norm_u == norm_previous_u)
00148 //     {
00149 //       std::cerr << "Warning, it appears u and previous_u are the "
00150 //              << "same, are you sure this is correct?"
00151 //              << "It's possible you forgot to set one or the other..."
00152 //              << std::endl;
00153 //     }
00154   
00155 //   Real delta_s_zero = std::sqrt(
00156 //                              (norm_u - norm_previous_u)*(norm_u - norm_previous_u) +
00157 //                              (*continuation_parameter-old_continuation_parameter)*
00158 //                              (*continuation_parameter-old_continuation_parameter)
00159 //                              );
00160 
00161 //   // 2.) Compute delta_s_zero as ||u -u_old|| + ...
00162 //   *delta_u = *solution;
00163 //   delta_u->add(-1., *previous_u);
00164 //   delta_u->close();
00165 //   Real norm_delta_u = delta_u->l2_norm();
00166 //   Real norm_u          = solution->l2_norm();
00167 //   Real norm_previous_u = previous_u->l2_norm();
00168 
00169 //   // Scale norm_delta_u by the bigger of either norm_u or norm_previous_u
00170 //   norm_delta_u /= std::max(norm_u, norm_previous_u);
00171   
00172 //   if (!quiet)
00173 //     {
00174 //       std::cout << "norm_u=" << norm_u << std::endl;
00175 //       std::cout << "norm_previous_u=" << norm_previous_u << std::endl;
00176 //       //std::cout << "norm_delta_u=" << norm_delta_u << std::endl;
00177 //       std::cout << "norm_delta_u/max(|u|,|u_old|)=" << norm_delta_u << std::endl;
00178 //       std::cout << "|norm_u-norm_previous_u|=" << std::abs(norm_u - norm_previous_u) << std::endl;
00179 //     }
00180   
00181 //   const Real dlambda = *continuation_parameter-old_continuation_parameter;
00182 
00183 //   if (!quiet)
00184 //     std::cout << "dlambda=" << dlambda << std::endl;
00185   
00186 //   Real delta_s_zero = std::sqrt(
00187 //                              (norm_delta_u*norm_delta_u) +
00188 //                              (dlambda*dlambda)
00189 //                              );
00190   
00191 //   if (!quiet)
00192 //     std::cout << "delta_s_zero=" << delta_s_zero << std::endl;
00193 
00194   // 1.) + 2.)
00195 //   // Now approximate the initial tangent d(lambda)/ds
00196 //   this->dlambda_ds = (*continuation_parameter-old_continuation_parameter) / delta_s_zero;
00197 
00198 
00199 //   // We can also approximate the deriv. wrt s by finite differences:
00200 //   // du/ds = (u1 - u0) / delta_s_zero.
00201 //   // FIXME: Use delta_u from above if we decide to keep that method.
00202 //   *du_ds = *solution;
00203 //   du_ds->add(-1., *previous_u);
00204 //   du_ds->scale(1./delta_s_zero);
00205 //   du_ds->close();
00206 
00207 
00208   // 3.) Treating (u-previous_u)/(lambda - lambda_old) as an approximation to du/d(lambda),
00209   // we follow the same technique as Carnes and Shadid.
00210 //   const Real dlambda = *continuation_parameter-old_continuation_parameter;
00211 //   libmesh_assert (dlambda > 0.);
00212   
00213 //   // Use delta_u for temporary calculation of du/d(lambda)
00214 //   *delta_u = *solution;
00215 //   delta_u->add(-1., *previous_u);
00216 //   delta_u->scale(1. / dlambda);
00217 //   delta_u->close();
00218 
00219 //   // Determine initial normalization parameter
00220 //   const Real solution_size = std::max(solution->l2_norm(), previous_u->l2_norm());
00221 //   if (solution_size > 1.)
00222 //     {
00223 //       Theta = 1./solution_size;
00224       
00225 //       if (!quiet)
00226 //      std::cout << "Setting Normalization Parameter Theta=" << Theta << std::endl;
00227 //     }
00228   
00229 //   // Compute d(lambda)/ds
00230 //   // The correct sign of d(lambda)/ds should be positive, since we assume that (lambda > lambda_old)
00231 //   // but we could always double-check that as well.
00232 //   Real norm_delta_u = delta_u->l2_norm();
00233 //   this->dlambda_ds = 1. / std::sqrt(1. + Theta*Theta*norm_delta_u*norm_delta_u);
00234 
00235 //   // Finally, compute du/ds = d(lambda)/ds * du/d(lambda)
00236 //   *du_ds = *delta_u;
00237 //   du_ds->scale(dlambda_ds);
00238 //   du_ds->close();
00239 
00240 
00241   // 4.) Use normalized arclength formula to estimate delta_s_zero
00242 //   // Determine initial normalization parameter
00243 //   set_Theta();
00244 
00245 //   // Compute (normalized) delta_s_zero
00246 //   *delta_u = *solution;
00247 //   delta_u->add(-1., *previous_u);
00248 //   delta_u->close();
00249 //   Real norm_delta_u = delta_u->l2_norm();
00250   
00251 //   const Real dlambda = *