ContinuationSystem Class Reference

#include <continuation_system.h>

Inheritance diagram for ContinuationSystem:

[legend]

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< DiffContext >	build_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_type &	system ()
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 MeshBase &	get_mesh () const
MeshBase &	get_mesh ()
const DofMap &	get_dof_map () const
DofMap &	get_dof_map ()
const EquationSystems &	get_equation_systems () const
EquationSystems &	get_equation_systems ()
bool	active () const
void	activate ()
void	deactivate ()
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 Variable &	variable (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 FEType &	variable_type (const unsigned int i) const
const FEType &	variable_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
Real *	continuation_parameter
bool	quiet
Real	continuation_parameter_tolerance
Real	solution_tolerance
Real	initial_newton_tolerance
Real	old_continuation_parameter
Real	min_continuation_parameter
Real	max_continuation_parameter
Real	Theta
Real	Theta_LOCA
unsigned int	n_backtrack_steps
unsigned int	n_arclength_reductions
Real	ds_min
Predictor	predictor
Real	newton_stepgrowth_aggressiveness
bool	newton_progress_check
bool	fe_reinit_during_postprocess
int	extra_quadrature_order
Real	numerical_jacobian_h
Real	verify_analytic_jacobians
bool	compute_internal_sides
bool	postprocess_sides
bool	assemble_qoi_sides
AutoPtr< TimeSolver >	time_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< Number >	qoi
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
NewtonSolver *	newton_solver
Real	dlambda_ds
Real	ds
Real	ds_current
Real	previous_dlambda_ds
Real	previous_ds
unsigned int	newton_step

---

Detailed Description

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

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

Author:

John W. Peterson 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.

typedef FEMSystem ContinuationSystem::Parent

The type of the parent.

Reimplemented from FEMSystem.

Definition at line 80 of file continuation_system.h.

typedef ContinuationSystem ContinuationSystem::sys_type

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

enum ContinuationSystem::Predictor

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

Enumerator:

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

Definition at line 222 of file continuation_system.h.

                 {
    Euler,

    AB2,

    Invalid_Predictor
  };

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.

                {Residual,
                 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.

  : Parent(es, name, number),
    continuation_parameter(NULL),
    quiet(true),
    continuation_parameter_tolerance(1.e-6),
    solution_tolerance(1.e-6),
    initial_newton_tolerance(0.01),
    old_continuation_parameter(0.),
    min_continuation_parameter(0.),
    max_continuation_parameter(0.),
    Theta(1.),
    Theta_LOCA(1.),
    //tau(1.),
    n_backtrack_steps(5),
    n_arclength_reductions(5),
    ds_min(1.e-8),
    predictor(Euler),
    newton_stepgrowth_aggressiveness(1.),
    newton_progress_check(true),
    rhs_mode(Residual),
    linear_solver(LinearSolver<Number>::build()),
    tangent_initialized(false),
    newton_solver(NULL),
    dlambda_ds(0.707),
    ds(0.1),
    ds_current(0.1),
    previous_dlambda_ds(0.),
    previous_ds(0.),
    newton_step(0)
{
  // Warn about using untested code
  libmesh_experimental();
}

ContinuationSystem::~ContinuationSystem

(

)

[virtual]

Destructor.

Definition at line 66 of file continuation_system.C.

References clear().

{
  this->clear();
}

---

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.

{
  _active = true;
}

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.

{
  return _active;  
}

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

{
  OStringStream adjoint_rhs_name;
  adjoint_rhs_name << "adjoint_rhs" << i;

  return this->add_vector(adjoint_rhs_name.str());
}

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

{
  OStringStream adjoint_name;
  adjoint_name << "adjoint_solution" << i;

  return this->add_vector(adjoint_name.str());
}

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

{
  // only add matrices before initializing...
  if (!_can_add_matrices)
    {
      std::cerr << "ERROR: Too late.  Cannot add matrices to the system after initialization"
                << std::endl
                << " any more.  You should have done this earlier."
                << std::endl;
      libmesh_error();
    }

  // Return the matrix if it is already there.
  if (this->have_matrix(mat_name))
    return *(_matrices[mat_name]);

  // Otherwise build the matrix and return it.
  SparseMatrix<Number>* buf = SparseMatrix<Number>::build().release();
  _matrices.insert (std::make_pair (mat_name, buf));

  return *buf;
}

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

{
  OStringStream sensitivity_rhs_name;
  sensitivity_rhs_name << "sensitivity_rhs" << i;

  return this->add_vector(sensitivity_rhs_name.str());
}

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

{
  OStringStream sensitivity_name;
  sensitivity_name << "sensitivity_solution" << i;

  return this->add_vector(sensitivity_name.str());
}

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

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

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

{  
  // Make sure the variable isn't there already
  // or if it is, that it's the type we want
  for (unsigned int v=0; v<this->n_vars(); v++)
    if (this->variable_name(v) == var)
      {
        if (this->variable_type(v) == type)
          return _variables[v].number();

        std::cerr << "ERROR: incompatible variable "
                  << var
                  << " has already been added for this system!"
                  << std::endl;
        libmesh_error();
      }

  const unsigned int curr_n_vars = this->n_vars();                                              

  // Add the variable to the list
  _variables.push_back((active_subdomains == NULL) ?
                       Variable(var, curr_n_vars, type) :
                       Variable(var, curr_n_vars, type, *active_subdomains));

  libmesh_assert ((curr_n_vars+1) == this->n_vars());

  _variable_numbers[var] = curr_n_vars;

  // Add the variable to the _dof_map
  _dof_map->add_variable (_variables.back());

  // Return the number of the new variable
  return curr_n_vars;
}

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

{
  // Return the vector if it is already there.
  if (this->have_vector(vec_name))
    return *(_vectors[vec_name]);

  // Otherwise build the vector
  NumericVector<Number>* buf = NumericVector<Number>::build().release();
  _vectors.insert (std::make_pair (vec_name, buf));
  _vector_projections.insert (std::make_pair (vec_name, projections));

  // Initialize it if necessary
  if (!_can_add_vectors)
    buf->init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);

  return *buf;
}

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

{
  OStringStream adjoint_name;
  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;

  return this->get_vector(adjoint_name.str());
}

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

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

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

{
  const unsigned int Np = parameters.size();
  const unsigned int Nq = qoi.size();

  // We currently get partial derivatives via central differencing
  const Real delta_p = TOLERANCE;

  // An introduction to the problem:
  //
  // Residual R(u(p),p) = 0
  // partial R / partial u = J = system matrix
  //
  // This implies that:
  // d/dp(R) = 0
  // (partial R / partial p) + 
  // (partial R / partial u) * (partial u / partial p) = 0

  // We first do an adjoint solve:
  // J^T * z = (partial q / partial u)

  this->adjoint_solve(qoi_indices);

  // Get ready to fill in senstivities:
  sensitivities.allocate_data(qoi_indices, *this, parameters);

  // We use the identities:
  // dq/dp = (partial q / partial p) + (partial q / partial u) *
  //         (partial u / partial p)
  // dq/dp = (partial q / partial p) + (J^T * z) *
  //         (partial u / partial p)
  // dq/dp = (partial q / partial p) + z * J *
  //         (partial u / partial p)
 
  // Leading to our final formula:
  // dq/dp = (partial q / partial p) - z * (partial R / partial p)

  for (unsigned int j=0; j != Np; ++j)
    {
      // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)
      // (partial R / partial p) ~= (rhs(p+dp) - rhs(p-dp))/(2*dp)

      Number old_parameter = *parameters[j];
      // Number old_qoi = this->qoi;

      *parameters[j] = old_parameter - delta_p;
      this->assemble_qoi(qoi_indices);
      std::vector<Number> qoi_minus = this->qoi;

      this->assembly(true, false);
      this->rhs->close();
      AutoPtr<NumericVector<Number> > partialR_partialp = this->rhs->clone();
      *partialR_partialp *= -1;

      *parameters[j] = old_parameter + delta_p;
      this->assemble_qoi(qoi_indices);
      std::vector<Number>& qoi_plus = this->qoi;

      std::vector<Number> partialq_partialp(Nq, 0);
      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);

      this->assembly(true, false);
      this->rhs->close();
      *partialR_partialp += *this->rhs;
      *partialR_partialp /= (2.*delta_p);

      // Don't leave the parameter changed
      *parameters[j] = old_parameter;

      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          sensitivities[i][j] = partialq_partialp[i] -
            partialR_partialp->dot(this->get_adjoint_solution(i));
    }

  // All parameters have been reset.
  // We didn't cache the original rhs or matrix for memory reasons,
  // but we can restore them to a state consistent solution -
  // principle of least surprise.
  this->assembly(true, true);
  this->rhs->close();
  this->matrix->close();
  this->assemble_qoi(qoi_indices);
}

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

{
  // Log how long the linear solve takes.
  START_LOG("adjoint_solve()", "ImplicitSystem");

  // The forward system should now already be solved.
  // Now assemble it's adjoint.

  if (this->assemble_before_solve)
    {
      // Build the Jacobian
      this->assembly(false, true);
      this->matrix->close();

      // Take the discrete adjoint
      matrix->get_transpose(*matrix);

      // Reset and build the RHS from the QOI derivative
      this->assemble_qoi_derivative(qoi_indices);
    }

  // The adjoint problem is linear
  LinearSolver<Number> *linear_solver = this->get_linear_solver();

  // Our iteration counts and residuals will be sums of the individual
  // results
  std::pair<unsigned int, Real> solver_params =
    this->get_linear_solve_parameters();
  std::pair<unsigned int, Real> totalrval = std::make_pair(0,0.0);

  for (unsigned int i=0; i != this->qoi.size(); ++i)
    if (qoi_indices.has_index(i))
      {
        const std::pair<unsigned int, Real> rval =
          linear_solver->solve (*matrix, this->add_adjoint_solution(i),
                                 this->get_adjoint_rhs(i),
                                 solver_params.second,
                                 solver_params.first);

        totalrval.first  += rval.first;
        totalrval.second += rval.second;
      }

  this->release_linear_solver(linear_solver);

  // The linear solver may not have fit our constraints exactly
#ifdef LIBMESH_ENABLE_AMR
  for (unsigned int i=0; i != this->qoi.size(); ++i)
    if (qoi_indices.has_index(i))
      this->get_dof_map().enforce_constraints_exactly
        (*this, &this->get_adjoint_solution(i));
#endif

  // Stop logging the nonlinear solve
  STOP_LOG("adjoint_solve()", "ImplicitSystem");

  return totalrval;
}

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

{
  // Solve for the updated tangent du1/ds, d(lambda1)/ds
  solve_tangent();

  // Advance the solution and the parameter to the next value.
  update_solution();
}

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

{
  if (predictor == Euler)
    {
      // 1.) Euler Predictor
      // Predict next the solution
      solution->add(ds_current, *du_ds);
      solution->close();
  
      // Predict next parameter value
      *continuation_parameter += ds_current*dlambda_ds;
    }


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

  else
    {
      // Unknown predictor
      libmesh_error();
    }
  
}

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

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

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

{
  START_LOG("assemble_qoi()", "FEMSystem");

  const MeshBase& mesh = this->get_mesh();

  this->update();

  AutoPtr<DiffContext> con = this->build_context();
  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);

  // the quantity of interest is assumed to be a sum of element and
  // side terms
  for (unsigned int i=0; i != qoi.size(); ++i)
    if (qoi_indices.has_index(i))
      qoi[i] = 0;

  // Loop over every active mesh element on this processor
  MeshBase::const_element_iterator el =
    mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el =
    mesh.active_local_elements_end();

  for ( ; el != end_el; ++el)
    {
      _femcontext.reinit(*this, *el);

      this->element_qoi(_femcontext, qoi_indices);

      for (_femcontext.side = 0;
           _femcontext.side != _femcontext.elem->n_sides();
           ++_femcontext.side)
        {
          // Don't compute on non-boundary sides unless requested
          if (!assemble_qoi_sides ||
              (!compute_internal_sides &&
               _femcontext.elem->neighbor(_femcontext.side) != NULL))
            continue;

          std::map<FEType, FEBase *>::iterator fe_end =
            _femcontext.side_fe.end();
          for (std::map<FEType, FEBase *>::iterator i =
            _femcontext.side_fe.begin();
               i != fe_end; ++i)
            {
              i->second->reinit(_femcontext.elem, _femcontext.side);
            }

          this->side_qoi(_femcontext, qoi_indices);
        }
    }

  Parallel::sum(_femcontext.elem_qoi);

  this->qoi = _femcontext.elem_qoi;

  STOP_LOG("assemble_qoi()", "FEMSystem");
}

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

{
  START_LOG("assemble_qoi_derivative()", "FEMSystem");

  const MeshBase& mesh = this->get_mesh();

  this->update();

  AutoPtr<DiffContext> con = this->build_context();
  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);

  // The quantity of interest derivative assembly accumulates on
  // initially zero vectors
  for (unsigned int i=0; i != qoi.size(); ++i)
    if (qoi_indices.has_index(i))
      this->add_adjoint_rhs(i).zero();

  // Loop over every active mesh element on this processor
  MeshBase::const_element_iterator el =
    mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el =
    mesh.active_local_elements_end();

  for ( ; el != end_el; ++el)
    {
      _femcontext.reinit(*this, *el);

      this->element_qoi_derivative(_femcontext, qoi_indices);

      for (_femcontext.side = 0;
           _femcontext.side != _femcontext.elem->n_sides();
           ++_femcontext.side)
        {
          // Don't compute on non-boundary sides unless requested
          if (!assemble_qoi_sides ||
              (!compute_internal_sides &&
               _femcontext.elem->neighbor(_femcontext.side) != NULL))
            continue;

          std::map<FEType, FEBase *>::iterator fe_end =
            _femcontext.side_fe.end();
          for (std::map<FEType, FEBase *>::iterator i =
            _femcontext.side_fe.begin();
               i != fe_end; ++i)
            {
              i->second->reinit(_femcontext.elem, _femcontext.side);
            }

          this->side_qoi_derivative(_femcontext, qoi_indices);
        }

      // We need some unmodified indices to use for constraining
      // multiple vector
      // FIXME - there should be a DofMap::constrain_element_vectors
      // to do this more efficiently
      std::vector<unsigned int> original_dofs = _femcontext.dof_indices;

      for (unsigned int i=0; i != qoi.size(); ++i)
        if (qoi_indices.has_index(i))
          {
            _femcontext.dof_indices = original_dofs;
            this->get_dof_map().constrain_element_vector
              (_femcontext.elem_qoi_derivative[i], _femcontext.dof_indices, false);

            this->get_adjoint_rhs(i).add_vector
              (_femcontext.elem_qoi_derivative[i], _femcontext.dof_indices);
          }
    }

  STOP_LOG("assemble_qoi_derivative()", "FEMSystem");
}

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

{
  const unsigned int Np = parameters.size();
  Real deltap = TOLERANCE;

  for (unsigned int p=0; p != Np; ++p)
    {
      NumericVector<Number> &sensitivity_rhs = this->add_sensitivity_rhs(p);

      // Approximate -(partial R / partial p) by
      // (R(p-dp) - R(p+dp)) / (2*dp)

      Number old_parameter = *parameters[p];
      *parameters[p] -= deltap;

      this->assembly(true, false);
      this->rhs->close();
      sensitivity_rhs = *this->rhs;

      *parameters[p] = old_parameter + deltap;

      this->assembly(true, false);
      this->rhs->close();

      sensitivity_rhs -= *this->rhs;
      sensitivity_rhs /= (2*deltap);
      sensitivity_rhs.close();

      *parameters[p] = old_parameter;
    }
}

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

{
  libmesh_assert(get_residual || get_jacobian);
  std::string log_name;
  if (get_residual && get_jacobian)
    log_name = "assembly()";
  else if (get_residual)
    log_name = "assembly(get_residual)";
  else
    log_name = "assembly(get_jacobian)";

  START_LOG(log_name, "FEMSystem");

  const MeshBase& mesh = this->get_mesh();

//  this->get_vector("_nonlinear_solution").localize
//    (*current_local_nonlinear_solution,
//     dof_map.get_send_list());
  this->update();

  if (print_solution_norms)
    {
//      this->get_vector("_nonlinear_solution").close();
      this->solution->close();
      std::cout << "|U| = "
//                << this->get_vector("_nonlinear_solution").l1_norm()
                << this->solution->l1_norm()
                << std::endl;
    }
  if (print_solutions)
    {
      unsigned int old_precision = std::cout.precision();
      std::cout.precision(16);
//      std::cout << "U = [" << this->get_vector("_nonlinear_solution")
      std::cout << "U = [" << *(this->solution)
                << "];" << std::endl;
      std::cout.precision(old_precision);
    }

  // Is this definitely necessary? [RHS]
  if (get_jacobian)
    matrix->zero();
  if (get_residual)
    rhs->zero();

  // Stupid C++ lets you set *Real* verify_analytic_jacobians = true!
  if (verify_analytic_jacobians > 0.5)
    {
      std::cerr << "WARNING!  verify_analytic_jacobians was set "
                << "to absurdly large value of "
                << verify_analytic_jacobians << std::endl;
      std::cerr << "Resetting to 1e-6!" << std::endl;
      verify_analytic_jacobians = 1e-6;
    }

  // In time-dependent problems, the nonlinear function we're trying
  // to solve at each timestep may depend on the particular solver
  // we're using
  libmesh_assert (time_solver.get() != NULL);

  AutoPtr<DiffContext> con = this->build_context();
  FEMContext &_femcontext = libmesh_cast_ref<FEMContext&>(*con);

  // Build the residual and jacobian contributions on every active
  // mesh element on this processor
  MeshBase::const_element_iterator el =
    mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el =
    mesh.active_local_elements_end();

  for ( ; el != end_el; ++el)
    {
      _femcontext.reinit(*this, *el);

      bool jacobian_computed =
        time_solver->element_residual(get_jacobian, _femcontext);

      // Compute a numeric jacobian if we have to
      if (get_jacobian && !jacobian_computed)
        {
          // Make sure we didn't compute a jacobian and lie about it
          libmesh_assert(_femcontext.elem_jacobian.l1_norm() == 0.0);
          // Logging of numerical jacobians is done separately
          this->numerical_elem_jacobian(_femcontext);
        }

      // Compute a numeric jacobian if we're asked to verify the
      // analytic jacobian we got
      if (get_jacobian && jacobian_computed &&
          verify_analytic_jacobians != 0.0)
        {
          DenseMatrix<Number> analytic_jacobian(_femcontext.elem_jacobian);

          _femcontext.elem_jacobian.zero();
          // Logging of numerical jacobians is done separately
          this->numerical_elem_jacobian(_femcontext);

          Real analytic_norm = analytic_jacobian.l1_norm();
          Real numerical_norm = _femcontext.elem_jacobian.l1_norm();

          // If we can continue, we'll probably prefer the analytic jacobian
          analytic_jacobian.swap(_femcontext.elem_jacobian);

          // The matrix "analytic_jacobian" will now hold the error matrix
          analytic_jacobian.add(-1.0, _femcontext.elem_jacobian);
          Real error_norm = analytic_jacobian.l1_norm();

          Real relative_error = error_norm /
                                std::max(analytic_norm, numerical_norm);

          if (relative_error > verify_analytic_jacobians)
            {
              std::cerr << "Relative error " << relative_error
                        << " detected in analytic jacobian on element "
                        << _femcontext.elem->id() << '!' << std::endl;

              unsigned int old_precision = std::cout.precision();
              std::cout.precision(16);
              std::cout << "J_analytic " << _femcontext.elem->id() << " = "
                        << _femcontext.elem_jacobian << std::endl;
              analytic_jacobian.add(1.0, _femcontext.elem_jacobian);
              std::cout << "J_numeric " << _femcontext.elem->id() << " = "
                        << analytic_jacobian << std::endl;

              std::cout.precision(old_precision);

              libmesh_error();
            }
        }

      for (_femcontext.side = 0;
           _femcontext.side != _femcontext.elem->n_sides();
           ++_femcontext.side)
        {
          // Don't compute on non-boundary sides unless requested
          if (!compute_internal_sides && 
              _femcontext.elem->neighbor(_femcontext.side) != NULL)
            continue;

          // Any mesh movement has already been done (and restored,
          // if the TimeSolver isn't broken), but
          // reinitializing the side FE objects is still necessary
          _femcontext.elem_side_fe_reinit();

          DenseMatrix<Number> old_jacobian;
          // If we're in DEBUG mode, we should always verify that the
          // user's side_residual function doesn't alter our existing
          // jacobian and then lie about it
#ifndef DEBUG
          // Even if we're not in DEBUG mode, when we're verifying
          // analytic jacobians we'll want to verify each side's
          // jacobian contribution separately
          if (verify_analytic_jacobians != 0.0 && get_jacobian)
#endif // ifndef DEBUG
            {
              old_jacobian = _femcontext.elem_jacobian;
              _femcontext.elem_jacobian.zero();
            }
          jacobian_computed =
            time_solver->side_residual(get_jacobian, _femcontext);

          // Compute a numeric jacobian if we have to
          if (get_jacobian && !jacobian_computed)
            {
              // In DEBUG mode, we've already set elem_jacobian == 0,
              // so we can make sure side_residual didn't compute a
              // jacobian and lie about it
#ifdef DEBUG
              libmesh_assert(_femcontext.elem_jacobian.l1_norm() == 0.0);
#endif
              // Logging of numerical jacobians is done separately
              this->numerical_side_jacobian(_femcontext);

              // If we're in DEBUG mode or if
              // verify_analytic_jacobians is on, we've moved
              // elem_jacobian's accumulated values into old_jacobian.
              // Now let's add them back.
#ifndef DEBUG
              if (verify_analytic_jacobians != 0.0)
#endif // ifndef DEBUG
                _femcontext.elem_jacobian += old_jacobian;
            }

          // Compute a numeric jacobian if we're asked to verify the
          // analytic jacobian we got
          if (get_jacobian && jacobian_computed &&
              verify_analytic_jacobians != 0.0)
            {
              DenseMatrix<Number> analytic_jacobian(_femcontext.elem_jacobian);

              _femcontext.elem_jacobian.zero();
              // Logging of numerical jacobians is done separately
              this->numerical_side_jacobian(_femcontext);

              Real analytic_norm = analytic_jacobian.l1_norm();
              Real numerical_norm = _femcontext.elem_jacobian.l1_norm();

              // If we can continue, we'll probably prefer the analytic jacobian
              analytic_jacobian.swap(_femcontext.elem_jacobian);

              // The matrix "analytic_jacobian" will now hold the error matrix
              analytic_jacobian.add(-1.0, _femcontext.elem_jacobian);
              Real error_norm = analytic_jacobian.l1_norm();

              Real relative_error = error_norm /
                                    std::max(analytic_norm, numerical_norm);

              if (relative_error > verify_analytic_jacobians)
                {
                  std::cerr << "Relative error " << relative_error
                            << " detected in analytic jacobian on element "
                            << _femcontext.elem->id()
                            << ", side "
                            << _femcontext.side << '!' << std::endl;

                  unsigned int old_precision = std::cout.precision();
                  std::cout.precision(16);
                  std::cout << "J_analytic " << _femcontext.elem->id() << " = "
                            << _femcontext.elem_jacobian << std::endl;
                  analytic_jacobian.add(1.0, _femcontext.elem_jacobian);
                  std::cout << "J_numeric " << _femcontext.elem->id() << " = "
                            << analytic_jacobian << std::endl;
                  std::cout.precision(old_precision);

                  libmesh_error();
                }
              // Once we've verified a side, we'll want to add back the
              // rest of the accumulated jacobian
              _femcontext.elem_jacobian += old_jacobian;
            }
          // In DEBUG mode, we've set elem_jacobian == 0, and we
          // may still need to add the old jacobian back
#ifdef DEBUG
          if (get_jacobian && jacobian_computed &&
              verify_analytic_jacobians == 0.0)
            {
              _femcontext.elem_jacobian += old_jacobian;
            }
#endif // ifdef DEBUG
        }

#ifdef LIBMESH_ENABLE_AMR
      // We turn off the asymmetric constraint application;
      // enforce_constraints_exactly() should be called in the solver
      if (get_residual && get_jacobian)
        this->get_dof_map().constrain_element_matrix_and_vector
          (_femcontext.elem_jacobian, _femcontext.elem_residual,
           _femcontext.dof_indices, false);
      else if (get_residual)
        this->get_dof_map().constrain_element_vector
          (_femcontext.elem_residual, _femcontext.dof_indices, false);
      else if (get_jacobian)
        this->get_dof_map().constrain_element_matrix
          (_femcontext.elem_jacobian, _femcontext.dof_indices, false);
#endif // #ifdef LIBMESH_ENABLE_AMR

      if (get_jacobian && print_element_jacobians)
        {
          unsigned int old_precision = std::cout.precision();
          std::cout.precision(16);
          std::cout << "J_elem " << _femcontext.elem->id() << " = "
                    << _femcontext.elem_jacobian << std::endl;
          std::cout.precision(old_precision);
        }

      if (get_jacobian)
        this->matrix->add_matrix (_femcontext.elem_jacobian,
                                  _femcontext.dof_indices);
      if (get_residual)
        this->rhs->add_vector (_femcontext.elem_residual,
                               _femcontext.dof_indices);
    }


  if (get_residual && print_residual_norms)
    {
      this->rhs->close();
      std::cout << "|F| = " << this->rhs->l1_norm() << std::endl;
    }
  if (get_residual && print_residuals)
    {
      this->rhs->close();
      unsigned int old_precision = std::cout.precision();
      std::cout.precision(16);
      std::cout << "F = [" << *(this->rhs) << "];" << std::endl;
      std::cout.precision(old_precision);
    }
  if (get_jacobian && print_jacobian_norms)
    {
      this->matrix->close();
      std::cout << "|J| = " << this->matrix->l1_norm() << std::endl;
    }
  if (get_jacobian && print_jacobians)
    {
      this->matrix->close();
      unsigned int old_precision = std::cout.precision();
      std::cout.precision(16);
      std::cout << "J = [" << *(this->matrix) << "];" << std::endl;
      std::cout.precision(old_precision);
    }
  STOP_LOG(log_name, "FEMSystem");
}

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.

{
  libmesh_assert (fptr != NULL);
  
  _assemble_system = fptr;  
}

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.

{
  libmesh_assert (fptr != NULL);
  
  _constrain_system = fptr;  
}

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.

{
  libmesh_assert (fptr != NULL);
  
  _init_system = fptr;
}

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

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

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

{
  // This function must be run on all processors at once
  parallel_only();

  START_LOG ("calculate_norm()", "System");

  // Zero the norm before summation
  Real v_norm = 0.;

  if (norm.is_discrete())
    {
      STOP_LOG ("calculate_norm()", "System");
      //Check to see if all weights are 1.0
      unsigned int check_var = 0;
      for (; check_var != this->n_vars(); ++check_var)
        if(norm.weight(check_var) != 1.0)
          break;

      //All weights were 1.0 so just do the full vector discrete norm
      if(check_var == this->n_vars())
        {
          FEMNormType norm_type = norm.type(0);
          
          if(norm_type == DISCRETE_L1)
            return v.l1_norm();
          if(norm_type == DISCRETE_L2)
            return v.l2_norm();
          if(norm_type == DISCRETE_L_INF)
            return v.linfty_norm();
          else
            libmesh_error();
        }

      for (unsigned int var=0; var != this->n_vars(); ++var)
        {
          // Skip any variables we don't need to integrate
          if (norm.weight(var) == 0.0)
            continue;

          v_norm += norm.weight(var) * discrete_var_norm(v, var, norm.type(var));
        }

      return v_norm;
    }

  // Localize the potentially parallel vector
  AutoPtr<NumericVector<Number> > local_v = NumericVector<Number>::build();
  local_v->init(v.size(), true, SERIAL);
  v.localize (*local_v, _dof_map->get_send_list()); 

  unsigned int dim = this->get_mesh().mesh_dimension();

  // Loop over all variables
  for (unsigned int var=0; var != this->n_vars(); ++var)
    {
      // Skip any variables we don't need to integrate
      if (norm.weight(var) == 0.0)
        continue;

      const FEType& fe_type = this->get_dof_map().variable_type(var);
      AutoPtr<QBase> qrule =
        fe_type.default_quadrature_rule (dim);
      AutoPtr<FEBase> fe
        (FEBase::build(dim, fe_type));
      fe->attach_quadrature_rule (qrule.get());

      const std::vector<Real>&               JxW = fe->get_JxW();
      const std::vector<std::vector<Real> >* phi = NULL;
      if (norm.type(var) == H1 ||
          norm.type(var) == H2 ||
          norm.type(var) == L2)
        phi = &(fe->get_phi());

      const std::vector<std::vector<RealGradient> >* dphi = NULL;
      if (norm.type(var) == H1 ||
          norm.type(var) == H2 ||
          norm.type(var) == H1_SEMINORM)
        dphi = &(fe->get_dphi());
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
      const std::vector<std::vector<RealTensor> >*   d2phi = NULL;
      if (norm.type(var) == H2 ||
          norm.type(var) == H2_SEMINORM)
        d2phi = &(fe->get_d2phi());
#endif

      std::vector<unsigned int> dof_indices;

      // Begin the loop over the elements
      MeshBase::const_element_iterator       el     =
        this->get_mesh().active_local_elements_begin();
      const MeshBase::const_element_iterator end_el =
        this->get_mesh().active_local_elements_end();

      for ( ; el != end_el; ++el)
        {
          const Elem* elem = *el;

          fe->reinit (elem);

          this->get_dof_map().dof_indices (elem, dof_indices, var);

          const unsigned int n_qp = qrule->n_points();

          const unsigned int n_sf = dof_indices.size();

          // Begin the loop over the Quadrature points.
          for (unsigned int qp=0; qp<n_qp; qp++)
            {
              if (norm.type(var) == H1 ||
                  norm.type(var) == H2 ||
                  norm.type(var) == L2)
                {
                  Number u_h = 0.;
                  for (unsigned int i=0; i != n_sf; ++i)
                    u_h += (*phi)[i][qp] * (*local_v)(dof_indices[i]);
                  v_norm += norm.weight_sq(var) *
                            JxW[qp] * libmesh_norm(u_h);
                }

              if (norm.type(var) == H1 ||
                  norm.type(var) == H2 ||
                  norm.type(var) == H1_SEMINORM)
                {
                  Gradient grad_u_h;
                  for (unsigned int i=0; i != n_sf; ++i)
                    grad_u_h.add_scaled((*dphi)[i][qp], (*local_v)(dof_indices[i]));
                  v_norm += norm.weight_sq(var) *
                            JxW[qp] * grad_u_h.size_sq();
                }

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
              if (norm.type(var) == H2 ||
                  norm.type(var) == H2_SEMINORM)
                {
                  Tensor hess_u_h;
                  for (unsigned int i=0; i != n_sf; ++i)
                    hess_u_h.add_scaled((*d2phi)[i][qp], (*local_v)(dof_indices[i]));
                  v_norm += norm.weight_sq(var) *
                            JxW[qp] * hess_u_h.size_sq();
                }
#endif
            }
        }
    }

  Parallel::sum(v_norm);

  STOP_LOG ("calculate_norm()", "System");

  return std::sqrt(v_norm);
}

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

{
  //short circuit to save time
  if(norm_type == DISCRETE_L1 ||
     norm_type == DISCRETE_L2 ||
     norm_type == DISCRETE_L_INF)
    return discrete_var_norm(v,var,norm_type);

  // Not a discrete norm
  std::vector<FEMNormType> norms(this->n_vars(), L2);
  std::vector<Real> weights(this->n_vars(), 0.0);
  norms[var] = norm_type;
  weights[var] = 1.0;
  Real val = this->calculate_norm(v, SystemNorm(norms, weights));
  return val;
}

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

{
  // FIXME: Do anything here, e.g. zero vectors, etc?

  // Call the Parent's clear function
  Parent::clear();
}

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

{
  // we do not care for matrices, but for vectors
  libmesh_assert (!_can_add_vectors);
  libmesh_assert (!other_system._can_add_vectors);

  if (verbose)
    {
      std::cout << "  Systems \"" << _sys_name << "\"" << std::endl;
      std::cout << "   comparing matrices not supported." << std::endl;
      std::cout << "   comparing names...";
    }

  // compare the name: 0 means identical
  const int name_result = _sys_name.compare(other_system.name());
  if (verbose)
    {
      if (name_result == 0)
        std::cout << " identical." << std::endl;
      else
        std::cout << "  names not identical." << std::endl;
      std::cout << "   comparing solution vector...";
    }


  // compare the solution: -1 means identical
  const int solu_result = solution->compare (*other_system.solution.get(),
                                             threshold);

  if (verbose)
    {
      if (solu_result == -1)
        std::cout << " identical up to threshold." << std::endl;
      else
        std::cout << "  first difference occured at index = " 
                  << solu_result << "." << std::endl;
    }


  // safety check, whether we handle at least the same number
  // of vectors
  std::vector<int> ov_result;

  if (this->n_vectors() != other_system.n_vectors())
    {
      if (verbose)
        {
          std::cout << "   Fatal difference. This system handles " 
                    << this->n_vectors() << " add'l vectors," << std::endl
                    << "   while the other system handles "
                    << other_system.n_vectors() 
                    << " add'l vectors." << std::endl
                    << "   Aborting comparison." << std::endl;
        }
      return false;
    }
  else if (this->n_vectors() == 0)
    {
      // there are no additional vectors...
      ov_result.clear ();
    }
  else
    {
      // compare other vectors
      for (const_vectors_iterator pos = _vectors.begin();
           pos != _vectors.end(); ++pos)
        {
          if (verbose)
              std::cout << "   comparing vector \""
                        << pos->first << "\" ...";

          // assume they have the same name
          const NumericVector<Number>& other_system_vector = 
              other_system.get_vector(pos->first);

          ov_result.push_back(pos->second->compare (other_system_vector,
                                                    threshold));

          if (verbose)
            {
              if (ov_result[ov_result.size()-1] == -1)
                std::cout << " identical up to threshold." << std::endl;
              else
                std::cout << " first difference occured at" << std::endl
                          << "   index = " << ov_result[ov_result.size()-1] << "." << std::endl;
            }

        }

    } // finished comparing additional vectors


  bool overall_result;
       
  // sum up the results
  if ((name_result==0) && (solu_result==-1))
    {
      if (ov_result.size()==0)
        overall_result = true;
      else
        {
          bool ov_identical;
          unsigned int n    = 0;
          do
            {
              ov_identical = (ov_result[n]==-1);
              n++;
            }
          while (ov_identical && n<ov_result.size());
          overall_result = ov_identical;
        }
    }
  else
    overall_result = false;

  if (verbose)
    {
      std::cout << "   finished comparisons, ";
      if (overall_result)
        std::cout << "found no differences." << std::endl << std::endl;
      else 
        std::cout << "found differences." << std::endl << std::endl;
    }
          
  return overall_result;
}

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.

{
  // Be sure the user has set the continuation parameter pointer
  if (!continuation_parameter)
    {
      std::cerr << "You must set the continuation_parameter pointer "
                << "to a member variable of the derived class, preferably in the "
                << "Derived class's init_data function.  This is how the ContinuationSystem "
                << "updates the continuation parameter."
                << std::endl;
      
      libmesh_error();
    }

  // Use extra precision for all the numbers printed in this function.
  unsigned int old_precision = std::cout.precision();
  std::cout.precision(16);
  std::cout.setf(std::ios_base::scientific);
  
  // We can't start solving the augmented PDE system unless the tangent
  // vectors have been initialized.  This only needs to occur once.
  if (!tangent_initialized)
    initialize_tangent();

  // Save the old value of -du/dlambda.  This will be used after the Newton iterations
  // to compute the angle between previous tangent vectors.  This cosine of this angle is
  //
  // tau := abs( (du/d(lambda)_i , du/d(lambda)_{i-1}) / (||du/d(lambda)_i|| * ||du/d(lambda)_{i-1}||) )
  //
  // The scaling factor tau (which should vary between 0 and 1) is used to shrink the step-size ds
  // when we are approaching a turning point.  Note that it can only shrink the step size.  
  *y_old = *y;
  
  // Set pointer to underlying Newton solver
  if (!newton_solver)
    newton_solver = libmesh_cast_ptr<NewtonSolver*> (this->time_solver->diff_solver().get());
  
  // A pair for catching return values from linear system solves.
  std::pair<unsigned int, Real> rval;

  // Convergence flag for the entire arcstep
  bool arcstep_converged = false;

  // Begin loop over arcstep reductions.
  for (unsigned int ns=0; ns<n_arclength_reductions; ++ns)
    {
      if (!quiet)
        {
          std::cout << "Current arclength stepsize, ds_current=" << ds_current << std::endl;
          std::cout << "Current parameter value, lambda=" << *continuation_parameter << std::endl;
        }
      
      // Upon exit from the nonlinear loop, the newton_converged flag
      // will tell us the convergence status of Newton's method.
      bool newton_converged = false;

      // The nonlinear residual before *any* nonlinear steps have been taken.
      Real nonlinear_residual_firststep = 0.;

      // The nonlinear residual from the current "k" Newton step, before the Newton step
      Real nonlinear_residual_beforestep = 0.;

      // The nonlinear residual from the current "k" Newton step, after the Newton step
      Real nonlinear_residual_afterstep = 0.;

      // The linear solver tolerance, can be updated dynamically at each Newton step.
      Real current_linear_tolerance = 0.;
      
      // The nonlinear loop
      for (newton_step=0; newton_step<newton_solver->max_nonlinear_iterations; ++newton_step)
        {
          std::cout << "\n === Starting Newton step " << newton_step << " ===" << std::endl;
          
          // Set the linear system solver tolerance
//        // 1.) Set the current linear tolerance based as a multiple of the current residual of the system.
//        const Real residual_multiple = 1.e-4;
//        Real current_linear_tolerance = residual_multiple*nonlinear_residual_beforestep;

//        // But if the current residual isn't small, don't let the solver exit with zero iterations!
//        if (current_linear_tolerance > 1.)
//          current_linear_tolerance = residual_multiple;

          // 2.) Set the current linear tolerance based on the method based on technique of Eisenstat & Walker.
          if (newton_step==0)
            {
              // At first step, only try reducing the residual by a small amount
              current_linear_tolerance = initial_newton_tolerance;//0.01;
            }

          else
            {
              // The new tolerance is based on the ratio of the most recent tolerances
              const Real alp=0.5*(1.+std::sqrt(5.));
              const Real gam=0.9;

              libmesh_assert (nonlinear_residual_beforestep != 0.0);
              libmesh_assert (nonlinear_residual_afterstep != 0.0);

              current_linear_tolerance = std::min(gam*std::pow(nonlinear_residual_afterstep/nonlinear_residual_beforestep, alp),
                                                  current_linear_tolerance*current_linear_tolerance
                                                  );

              // Don't let it get ridiculously small!!
              if (current_linear_tolerance < 1.e-12)
                current_linear_tolerance = 1.e-12;
            }

          if (!quiet)
            std::cout << "Using current_linear_tolerance=" << current_linear_tolerance << std::endl;

          
          // Assemble the residual (and Jacobian).
          rhs_mode = Residual;
          assembly(true,   // Residual
                   true); // Jacobian
          rhs->close();

          // Save the current nonlinear residual.  We don't need to recompute the residual unless
          // this is the first step, since it was already computed as part of the convergence check
          // at the end of the last loop iteration.
          if (newton_step==0)
            {
              nonlinear_residual_beforestep = rhs->l2_norm();

              // Store the residual before any steps have been taken.  This will *not*
              // be updated at each step, and can be used to see if any progress has
              // been made from the initial residual at later steps.
              nonlinear_residual_firststep = nonlinear_residual_beforestep;

              const Real old_norm_u = solution->l2_norm();
              std::cout << "  (before step) ||R||_{L2} = " << nonlinear_residual_beforestep << std::endl;
              std::cout << "  (before step) ||R||_{L2}/||u|| = " << nonlinear_residual_beforestep / old_norm_u << std::endl;

              // In rare cases (very small arcsteps), it's possible that the residual is
              // already below our absolute linear tolerance.
              if (nonlinear_residual_beforestep  < solution_tolerance)
                {
                  if (!quiet)
                    std::cout << "Initial guess satisfied linear tolerance, exiting with zero Newton iterations!" << std::endl;

                  // Since we go straight from here to the solve of the next tangent, we
                  // have to close the matrix before it can be assembled again.
                  matrix->close();
                  newton_converged=true;
                  break; // out of Newton iterations, with newton_converged=true
                }
            }

          else
            {
              nonlinear_residual_beforestep = nonlinear_residual_afterstep;
            }

          
          // Solve the linear system G_u*z = G
          // Initial guess?
          z->zero(); // It seems to be extremely important to zero z here, otherwise the solver quits early.
          z->close();
          
          // It's possible that we have selected the current_linear_tolerance so large that
          // a guess of z=zero yields a linear system residual |Az + R| small enough that the
          // linear solver exits in zero iterations.  If this happens, we will reduce the
          // current_linear_tolerance until the linear solver does at least 1 iteration.
          do
            {
              rval =
                linear_solver->solve(*matrix,
                                     *z,
                                     *rhs,
                                     //1.e-12,
                                     current_linear_tolerance,
                                     newton_solver->max_linear_iterations);   // max linear iterations
              
              if (rval.first==0)
                {
                  if (newton_step==0)
                    {
                      std::cout << "Repeating initial solve with smaller linear tolerance!" << std::endl;
                      current_linear_tolerance *= initial_newton_tolerance; // reduce the linear tolerance to force the solver to do some work
                    }
                  else
                    {
                      // We shouldn't get here ... it means the linear solver did no work on a Newton
                      // step other than the first one.  If this happens, we need to think more about our
                      // tolerance selection.
                      libmesh_error();
                    }
                }
              
            } while (rval.first==0);

          
          if (!quiet)
            std::cout << "  G_u*z = G solver converged at step "
                      << rval.first
                      << " linear tolerance = "
                      << rval.second
                      << "."
                      << std::endl;

          // Sometimes (I am not sure why) the linear solver exits after zero iterations.
          // Perhaps it is hitting PETSc's divergence tolerance dtol???  If this occurs,
          // we should break out of the Newton iteration loop because nothing further is
          // going to happen...  Of course if the tolerance is already small enough after
          // zero iterations (how can this happen?!) we should not quit.
          if ((rval.first == 0) && (rval.second > current_linear_tolerance*nonlinear_residual_beforestep))
            {
              if (!quiet)
                std::cout << "Linear solver exited in zero iterations!" << std::endl;

              // Try to find out the reason for convergence/divergence
              linear_solver->print_converged_reason();
              
              break; // out of Newton iterations
            }
          
          // Note: need to scale z by -1 since our code always solves Jx=R
          // instead of Jx=-R.
          z->scale(-1.);
          z->close();





          
          // Assemble the G_Lambda vector, skip residual.
          rhs_mode = G_Lambda;

          // Assemble both rhs and Jacobian
          assembly(true,  // Residual
                   false); // Jacobian

          // Not sure if this is really necessary
          rhs->close();
          const Real yrhsnorm=rhs->l2_norm();
          if (yrhsnorm == 0.0)
            {
              std::cout << "||G_Lambda|| = 0" << std::endl;
              libmesh_error();
            }

          // We select a tolerance for the y-system which is based on the inexact Newton
          // tolerance but scaled by an extra term proportional to the RHS (which is not -> 0 in this case)
          const Real ysystemtol=current_linear_tolerance*(nonlinear_residual_beforestep/yrhsnorm);
          if (!quiet)
            std::cout << "ysystemtol=" << ysystemtol << std::endl;
          
          // Solve G_u*y = G_{\lambda}
          // FIXME: Initial guess?  This is really a solve for -du/dlambda so we could try
          // initializing it with the latest approximation to that... du/dlambda ~ du/ds * ds/dlambda
          //*y = *solution;
          //y->add(-1., *previous_u);
          //y->scale(-1. / (*continuation_parameter - old_continuation_parameter)); // Be careful of divide by zero...
          //y->close();
          
          //      const unsigned int max_attempts=1;
          // unsigned int attempt=0;
          //      do
          //        {
          //          if (!quiet)
          //            std::cout << "Trying to solve tangent system, attempt " << attempt << std::endl;
              
              rval =
                linear_solver->solve(*matrix,
                                     *y,
                                     *rhs,
                                     //1.e-12, 
                                     ysystemtol,
                                     newton_solver->max_linear_iterations);   // max linear iterations

              if (!quiet)
                std::cout << "  G_u*y = G_{lambda} solver converged at step "
                          << rval.first
                          << ", linear tolerance = "
                          << rval.second
                          << "."
                          << std::endl;

              // Sometimes (I am not sure why) the linear solver exits after zero iterations.
              // Perhaps it is hitting PETSc's divergence tolerance dtol???  If this occurs,
              // we should break out of the Newton iteration loop because nothing further is
              // going to happen...
              if ((rval.first == 0) && (rval.second > ysystemtol))
                {
                  if (!quiet)
                    std::cout << "Linear solver exited in zero iterations!" << std::endl;

                  break; // out of Newton iterations
                }
              
//            ++attempt;
//          } while ((attempt<max_attempts) && (rval.first==newton_solver->max_linear_iterations));
            
      
      
      
      
          // Compute N, the residual of the arclength constraint eqn.
          // Note 1: N(u,lambda,s) := (u-u_{old}, du_ds) + (lambda-lambda_{old}, dlambda_ds) - _ds 
          // We temporarily use the delta_u vector as a temporary vector for this calculation.
          *delta_u = *solution;
          delta_u->add(-1., *previous_u);

          // First part of the arclength constraint
          const Number N1 = Theta_LOCA*Theta_LOCA*Theta*delta_u->dot(*du_ds);
          const Number N2 = ((*continuation_parameter) - old_continuation_parameter)*dlambda_ds;
          const Number N3 = ds_current;

          if (!quiet)
            {
              std::cout << "  N1=" << N1 << std::endl;
              std::cout << "  N2=" << N2 << std::endl;
              std::cout << "  N3=" << N3 << std::endl;
            }
      
          // The arclength constraint value 
          const Number N = N1+N2-N3;
      
          if (!quiet)
            std::cout << "  N=" << N << std::endl;

          const Number duds_dot_z = du_ds->dot(*z);
          const Number duds_dot_y = du_ds->dot(*y);

          //std::cout << "duds_dot_z=" << duds_dot_z << std::endl;
          //std::cout << "duds_dot_y=" << duds_dot_y << std::endl;
          //std::cout << "dlambda_ds=" << dlambda_ds << std::endl;

          const Number delta_lambda_numerator   = -(N          + Theta_LOCA*Theta_LOCA*Theta*duds_dot_z);
          const Number delta_lambda_denominator =  (dlambda_ds - Theta_LOCA*Theta_LOCA*Theta*duds_dot_y);

          libmesh_assert (delta_lambda_denominator != 0.0);

          // Now, we are ready to compute the step delta_lambda
          const Number delta_lambda_comp = delta_lambda_numerator /
                                           delta_lambda_denominator;
          // Lambda is real-valued
          const Real delta_lambda = libmesh_real(delta_lambda_comp);

          // Knowing delta_lambda, we are ready to update delta_u
          // delta_u = z - delta_lambda*y
          delta_u->zero();
          delta_u->add(1., *z);
          delta_u->add(-delta_lambda, *y);
          delta_u->close();

          // Update the system solution and the continuation parameter.
          solution->add(1., *delta_u);
          solution->close();
          *continuation_parameter += delta_lambda;

          // Did the Newton step actually reduce the residual?
          rhs_mode = Residual;
          assembly(true,   // Residual
                   false); // Jacobian
          rhs->close();
          nonlinear_residual_afterstep = rhs->l2_norm();

          
          // In a "normal" Newton step, ||du||/||R|| > 1 since the most recent
          // step is where you "just were" and the current residual is where
          // you are now.  It can occur that ||du||/||R|| < 1, but these are
          // likely not good cases to attempt backtracking (?).
          const Real norm_du_norm_R = delta_u->l2_norm() / nonlinear_residual_afterstep;
          if (!quiet)
            std::cout << "  norm_du_norm_R=" << norm_du_norm_R << std::endl;
      
      
          // Factor to decrease the stepsize by for backtracking
          Real newton_stepfactor = 1.;

          const bool attempt_backtracking =
            (nonlinear_residual_afterstep > solution_tolerance)
            && (nonlinear_residual_afterstep > nonlinear_residual_beforestep)
            && (n_backtrack_steps>0)
            && (norm_du_norm_R > 1.)
            ;
    
          // If residual is not reduced, do Newton back tracking.  
          if (attempt_backtracking)
            {
              if (!quiet)
                std::cout << "Newton step did not reduce residual." << std::endl;

              // back off the previous step.
              solution->add(-1., *delta_u);
              solution->close();
              *continuation_parameter -= delta_lambda;

              // Backtracking: start cutting the Newton stepsize by halves until
              // the new residual is actually smaller...
              for (unsigned int backtrack_step=0; backtrack_step<n_backtrack_steps; ++backtrack_step)
                {
                  newton_stepfactor *= 0.5;

                  if (!quiet)
                    std::cout << "Shrinking step size by " << newton_stepfactor << std::endl;
              
                  // Take fractional step
                  solution->add(newton_stepfactor, *delta_u);
                  solution->close();
                  *continuation_parameter += newton_stepfactor*delta_lambda;
              
                  rhs_mode = Residual;
                  assembly(true,   // Residual
                           false); // Jacobian
                  rhs->close();
                  nonlinear_residual_afterstep = rhs->l2_norm();

                  if (!quiet)
                    std::cout << "At shrink step "
                              << backtrack_step
                              << ", nonlinear_residual_afterstep="
                              << nonlinear_residual_afterstep
                              << std::endl;

                  if (nonlinear_residual_afterstep < nonlinear_residual_beforestep)
                    {
                      if (!quiet)
                        std::cout << "Backtracking succeeded!" << std::endl;
                      
                      break; // out of backtracking loop
                    }
                  
                  else
                    {
                      // Back off that step
                      solution->add(-newton_stepfactor, *delta_u);
                      solution->close();
                      *continuation_parameter -= newton_stepfactor*delta_lambda;
                    }
              
                  // Save a copy of the solution from before the Newton step.
                  //AutoPtr<NumericVector<Number> > prior_iterate = solution->clone();
                }
            } // end if (attempte_backtracking)


          // If we tried backtracking but the residual is still not reduced, print message.
          if ((attempt_backtracking) && (nonlinear_residual_afterstep > nonlinear_residual_beforestep))
            {
              //std::cerr << "Backtracking failed." << std::endl;
              std::cout << "Backtracking failed." << std::endl;
          
              // 1.) Quit, exit program.
              //libmesh_error();

              // 2.) Continue with last newton_stepfactor
              if (newton_step<3)
                {
                  solution->add(newton_stepfactor, *delta_u);
                  solution->close();
                  *continuation_parameter += newton_stepfactor*delta_lambda;
                  if (!quiet)
                    std::cout << "Backtracking could not reduce residual ... continuing anyway!" << std::endl;
                }
              
              // 3.) Break out of Newton iteration loop with newton_converged = false,
              //     reduce the arclength stepsize, and try again.
              else
                {
                  break; // out of Newton iteration loop, with newton_converged=false
                }
            }

          // Another type of convergence check: suppose the residual has not been reduced
          // from its initial value after half of the allowed Newton steps have occurred.
          // In our experience, this typically means that it isn't going to converge and
          // we could probably save time by dropping out of the Newton iteration loop and
          // trying a smaller arcstep.
          if (this->newton_progress_check)
            {
              if ((nonlinear_residual_afterstep > nonlinear_residual_firststep) &&
                  (newton_step+1 > static_cast<unsigned int>(0.5*newton_solver->max_nonlinear_iterations)))
                {
                  std::cout << "Progress check failed: the current residual: "
                            << nonlinear_residual_afterstep
                            << ", is\n"
                            << "larger than the initial residual, and half of the allowed\n"
                            << "number of Newton iterations have elapsed.\n"
                            << "Exiting Newton iterations with converged==false." << std::endl;
              
                  break; // out of Newton iteration loop, newton_converged = false
                }
            }
          
          // Safety check: Check the current continuation parameter against user-provided min-allowable parameter value
          if (*continuation_parameter < min_continuation_parameter)
            {
              std::cout << "Continuation parameter fell below min-allowable value." << std::endl;
              // libmesh_error();
              break; // out of Newton iteration loop, newton_converged = false
            }
          
          // Safety check: Check the current continuation parameter against user-provided max-allowable parameter value
          if ( (max_continuation_parameter != 0.0) &&
               (*continuation_parameter > max_continuation_parameter) )
            {
              std::cout << "Current continuation parameter value: "
                        << *continuation_parameter
                        << " exceeded max-allowable value."
                        << std::endl;
              // libmesh_error();
              break; // out of Newton iteration loop, newton_converged = false
            }
          
      
          // Check the convergence of the parameter and the solution.  If they are small
          // enough, we can break out of the Newton iteration loop.  
          const Real norm_delta_u = delta_u->l2_norm();
          const Real norm_u = solution->l2_norm();
          std::cout << "  delta_lambda                   = " << delta_lambda << std::endl;
          std::cout << "  newton_stepfactor*delta_lambda = " << newton_stepfactor*delta_lambda << std::endl;
          std::cout << "  lambda_current                 = " << *continuation_parameter << std::endl;
          std::cout << "  ||delta_u||                    = " << norm_delta_u << std::endl;
          std::cout << "  ||delta_u||/||u||              = " << norm_delta_u / norm_u << std::endl;
      

          // Evaluate the residual at the current Newton iterate.  We don't want to detect
          // convergence due to a small Newton step when the residual is still not small.
          rhs_mode = Residual;
          assembly(true,   // Residual
                   false); // Jacobian
          rhs->close();
          const Real norm_residual = rhs->l2_norm();
          std::cout << "  ||R||_{L2} = " << norm_residual << std::endl;
          std::cout << "  ||R||_{L2}/||u|| = " << norm_residual / norm_u << std::endl;
      
      
          // FIXME: The norm_delta_u tolerance (at least) should be relative.
          // It doesn't make sense to converge a solution whose size is ~ 10^5 to
          // a tolerance of 1.e-6.  Oh, and we should also probably check the
          // (relative) size of the residual as well, instead of just the step.
          if ((std::abs(delta_lambda) < continuation_parameter_tolerance) &&
              //(norm_delta_u       < solution_tolerance)               && // This is a *very* strict criterion we can probably skip
              (norm_residual      < solution_tolerance))
            {
              if (!quiet)
                std::cout << "Newton iterations converged!" << std::endl;

              newton_converged = true;
              break; // out of Newton iterations
            }
        } // end nonlinear loop

      if (!newton_converged)
        {
          std::cout << "Newton iterations of augmented system did not converge!" << std::endl;
          
          // Reduce ds_current, recompute the solution and parameter, and continue to next
          // arcstep, if there is one.
          ds_current *= 0.5;

          // Go back to previous solution and parameter value.
          *solution = *previous_u;
          *continuation_parameter = old_continuation_parameter;

          // Compute new predictor with smaller ds
          apply_predictor();
        }
      else
        {
          // Set step convergence and break out
          arcstep_converged=true;
          break; // out of arclength reduction loop
        }
      
    } // end loop over arclength reductions

  // Check for convergence of the whole arcstep.  If not converged at this
  // point, we have no choice but to quit.
  if (!arcstep_converged)
    {
      std::cout << "Arcstep failed to converge after max number of reductions! Exiting..." << std::endl;
      libmesh_error();
    }
  
  // Print converged solution control parameter and max value.
  std::cout << "lambda_current=" << *continuation_parameter << std::endl;
  //std::cout << "u_max=" << solution->max() << std::endl;

  // Reset old stream precision and flags.
  std::cout.precision(old_precision);
  std::cout.unsetf(std::ios_base::scientific);

  // Note: we don't want to go on to the next guess yet, since the user may
  // want to post-process this data.  It's up to the user to call advance_arcstep()
  // when they are ready to go on.
}

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

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

void System::deactivate

(

)

[inline, inherited]

Deactivates the system. Only active systems are solved.

Definition at line 1382 of file system.h.

References System::_active.

{
  _active = false;
}

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

                                                  {
    return request_jacobian;
  }

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

{}

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

    {}

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

    {}

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

                                                       {
    return request_jacobian;
  }

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.

{
  // Only calculate a mesh movement residual if it's necessary
  if (_mesh_sys == libMesh::invalid_uint)
    return request_jacobian;

  FEMContext &context = libmesh_cast_ref<FEMContext&>(c);

  // This function only supports fully coupled mesh motion for now
  libmesh_assert(_mesh_sys == this->number());

  unsigned int n_qpoints = context.element_qrule->n_points();

  const unsigned int n_x_dofs = (_mesh_x_var == libMesh::invalid_uint) ?
                                0 : context.dof_indices_var[_mesh_x_var].size();
  const unsigned int n_y_dofs = (_mesh_y_var == libMesh::invalid_uint) ?
                                0 : context.dof_indices_var[_mesh_y_var].size();
  const unsigned int n_z_dofs = (_mesh_z_var == libMesh::invalid_uint) ?
                                0 : context.dof_indices_var[_mesh_z_var].size();

  const unsigned int mesh_xyz_var = n_x_dofs ? _mesh_x_var :
                                   (n_y_dofs ? _mesh_y_var :
                                   (n_z_dofs ? _mesh_z_var :
                                    libMesh::invalid_uint));

  // If we're our own _mesh_sys, we'd better be in charge of
  // at least one coordinate, and we'd better have the same
  // FE type for all coordinates we are in charge of
  libmesh_assert(mesh_xyz_var != libMesh::invalid_uint);
  libmesh_assert(!n_x_dofs || context.element_fe_var[_mesh_x_var] ==
                              context.element_fe_var[mesh_xyz_var]);
  libmesh_assert(!n_y_dofs || context.element_fe_var[_mesh_y_var] ==
                              context.element_fe_var[mesh_xyz_var]);
  libmesh_assert(!n_z_dofs || context.element_fe_var[_mesh_z_var] ==
                              context.element_fe_var[mesh_xyz_var]);

  const std::vector<std::vector<Real> >     &psi =
    context.element_fe_var[mesh_xyz_var]->get_phi();

  for (unsigned int var = 0; var != this->n_vars(); ++var)
    {
      // Mesh motion only affects time-evolving variables
      if (!_time_evolving[var])
        continue;

      // The mesh coordinate variables themselves are Lagrangian,
      // not Eulerian, and no convective term is desired.
      if (_mesh_sys == this->number() &&
          (var == _mesh_x_var ||
           var == _mesh_y_var ||
           var == _mesh_z_var))
        continue;

      // Some of this code currently relies on the assumption that
      // we can pull mesh coordinate data from our own system
      if (_mesh_sys != this->number())
        libmesh_not_implemented();

      // This residual should only be called by unsteady solvers:
      // if the mesh is steady, there's no mesh convection term!
      UnsteadySolver *unsteady =
        dynamic_cast<UnsteadySolver *>(this->time_solver.get());
      if (!unsteady)
        return request_jacobian;

      const std::vector<Real> &JxW = 
        context.element_fe_var[var]->get_JxW();

      const std::vector<std::vector<Real> >     &phi =
        context.element_fe_var[var]->get_phi();

      const std::vector<std::vector<RealGradient> > &dphi =
        context.element_fe_var[var]->get_dphi();

      const unsigned int n_u_dofs = context.dof_indices_var[var].size();

      DenseSubVector<Number> &Fu = *context.elem_subresiduals[var];
      DenseSubMatrix<Number> &Kuu = *context.elem_subjacobians[var][var];

      DenseSubMatrix<Number> *Kux = n_x_dofs ?
        context.elem_subjacobians[var][_mesh_x_var] : NULL;
      DenseSubMatrix<Number> *Kuy = n_y_dofs ?
        context.elem_subjacobians[var][_mesh_y_var] : NULL;
      DenseSubMatrix<Number> *Kuz = n_z_dofs ?
        context.elem_subjacobians[var][_mesh_z_var] : NULL;

      std::vector<Real> delta_x(n_x_dofs, 0.);
      std::vector<Real> delta_y(n_y_dofs, 0.);
      std::vector<Real> delta_z(n_z_dofs, 0.);

      for (unsigned int i = 0; i != n_x_dofs; ++i)
        {
          unsigned int j = context.dof_indices_var[_mesh_x_var][i];
          delta_x[i] = libmesh_real(this->current_solution(j)) -
                       libmesh_real(unsteady->old_nonlinear_solution(j));
        }

      for (unsigned int i = 0; i != n_y_dofs; ++i)
        {
          unsigned int j = context.dof_indices_var[_mesh_y_var][i];
          delta_y[i] = libmesh_real(this->current_solution(j)) -
                       libmesh_real(unsteady->old_nonlinear_solution(j));
        }

      for (unsigned int i = 0; i != n_z_dofs; ++i)
        {
          unsigned int j = context.dof_indices_var[_mesh_z_var][i];
          delta_z[i] = libmesh_real(this->current_solution(j)) -
                       libmesh_real(unsteady->old_nonlinear_solution(j));
        }

      for (unsigned int qp = 0; qp != n_qpoints; ++qp)
        {
          Gradient grad_u = context.interior_gradient(var, qp);
          RealGradient convection(0.);

          for (unsigned int i = 0; i != n_x_dofs; ++i)
            convection(0) += delta_x[i] * psi[i][qp];
          for (unsigned int i = 0; i != n_y_dofs; ++i)
            convection(1) += delta_y[i] * psi[i][qp];
          for (unsigned int i = 0; i != n_z_dofs; ++i)
            convection(2) += delta_z[i] * psi[i][qp];

          for (unsigned int i = 0; i != n_u_dofs; ++i)
            {
              Number JxWxPhiI = JxW[qp] * phi[i][qp];
              Fu(i) += (convection * grad_u) * JxWxPhiI;
              if (request_jacobian)
                {
                  Number JxWxPhiI = JxW[qp] * phi[i][qp];
                  for (unsigned int j = 0; j != n_u_dofs; ++j)
                    Kuu(i,j) += JxWxPhiI * (convection * dphi[j][qp]);

                  Number JxWxPhiIoverDT = JxWxPhiI/this->deltat;

                  Number JxWxPhiIxDUDXoverDT = JxWxPhiIoverDT * grad_u(0);
                  for (unsigned int j = 0; j != n_x_dofs; ++j)
                    (*Kux)(i,j) += JxWxPhiIxDUDXoverDT * psi[j][qp];

                  Number JxWxPhiIxDUDYoverDT = JxWxPhiIoverDT * grad_u(1);
                  for (unsigned int j = 0; j != n_y_dofs; ++j)
                    (*Kuy)(i,j) += JxWxPhiIxDUDYoverDT * psi[j][qp];

                  Number JxWxPhiIxDUDZoverDT = JxWxPhiIoverDT * grad_u(2);
                  for (unsigned int j = 0; j != n_z_dofs; ++j)
                    (*Kuz)(i,j) += JxWxPhiIxDUDZoverDT * psi[j][qp];
                }
            }
        }
    }

  return request_jacobian;
}

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

{
  const unsigned int Np = parameters.size();
  const unsigned int Nq = qoi.size();

  // We currently get partial derivatives via central differencing
  const Real delta_p = TOLERANCE;

  // An introduction to the problem:
  //
  // Residual R(u(p),p) = 0
  // partial R / partial u = J = system matrix
  //
  // This implies that:
  // d/dp(R) = 0
  // (partial R / partial p) + 
  // (partial R / partial u) * (partial u / partial p) = 0

  // We first solve for (partial u / partial p) for each parameter:
  // J * (partial u / partial p) = - (partial R / partial p)

  this->sensitivity_solve(parameters);

  // Get ready to fill in senstivities:
  sensitivities.allocate_data(qoi_indices, *this, parameters);

  // We use the identity:
  // dq/dp = (partial q / partial p) + (partial q / partial u) *
  //         (partial u / partial p)
 
  // We get (partial q / partial u) from the user
  this->assemble_qoi_derivative(qoi_indices);

  for (unsigned int j=0; j != Np; ++j)
    {
      // (partial q / partial p) ~= (q(p+dp)-q(p-dp))/(2*dp)

      Number old_parameter = *parameters[j];

      *parameters[j] = old_parameter - delta_p;
      this->assemble_qoi();
      std::vector<Number> qoi_minus = this->qoi;

      *parameters[j] = old_parameter + delta_p;
      this->assemble_qoi();
      std::vector<Number>& qoi_plus = this->qoi;

      std::vector<Number> partialq_partialp(Nq, 0);
      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          partialq_partialp[i] = (qoi_plus[i] - qoi_minus[i]) / (2.*delta_p);

      // Don't leave the parameter changed
      *parameters[j] = old_parameter;

      for (unsigned int i=0; i != Nq; ++i)
        if (qoi_indices.has_index(i))
          sensitivities[i][j] = partialq_partialp[i] +
            this->get_adjoint_rhs(i).dot(this->get_sensitivity_solution(i));
    }

  // All parameters have been reset.
  // We didn't cache the original rhs or matrix for memory reasons,
  // but we can restore them to a state consistent solution -
  // principle of least surprise.
  this->assembly(true, true);
  this->rhs->close();
  this->matrix->close();
  this->assemble_qoi(qoi_indices);
}

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

{
  OStringStream adjoint_rhs_name;
  adjoint_rhs_name << "adjoint_rhs" << i;

  return this->get_vector(adjoint_rhs_name.str());
}

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

{
  OStringStream adjoint_rhs_name;
  adjoint_rhs_name << "adjoint_rhs" << i;

  return this->get_vector(adjoint_rhs_name.str());
}

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

{
  OStringStream adjoint_name;
  adjoint_name << "adjoint_solution" << i;

  return this->get_vector(adjoint_name.str());
}

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

{
  OStringStream adjoint_name;
  adjoint_name << "adjoint_solution" << i;

  return this->get_vector(adjoint_name.str());
}

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.

{
  return *_dof_map;
}

const DofMap & System::get_dof_map

(

)

const [inline, inherited]

Returns:

a constant reference to this system's _dof_map.

Definition at line 1350 of file system.h.

References System::_dof_map.

Referenced by __libmesh_petsc_diff_solver_jacobian(), __libmesh_petsc_diff_solver_residual(), ExactSolution::_compute_error(), HPCoarsenTest::add_projection(), ImplicitSystem::adjoint_solve(), FEMSystem::assemble_qoi_derivative(), FEMSystem::assembly(), EquationSystems::build_discontinuous_solution_vector(), EquationSystems::build_solution_vector(), System::calculate_norm(), DofMap::enforce_constraints_exactly(), JumpErrorEstimator::estimate_error(), ExactErrorEstimator::estimate_error(), System::get_info(), EigenSystem::init_data(), ImplicitSystem::init_matrices(), System::local_dof_indices(), DofMap::max_constraint_error(), FEMSystem::mesh_position_get(), UnsteadySolver::old_nonlinear_solution(), System::ProjectVector::operator()(), PatchRecoveryErrorEstimator::EstimateError::operator()(), ErrorVector::plot_error(), System::project_vector(), FEMContext::reinit(), EquationSystems::reinit(), HPCoarsenTest::select_refinement(), ImplicitSystem::sensitivity_solve(), UnsteadySolver::solve(), PetscDiffSolver::solve(), NewtonSolver::solve(), ImplicitSystem::weighted_sensitivity_adjoint_solve(), ImplicitSystem::weighted_sensitivity_solve(), EnsightIO::write_scalar_ascii(), and EnsightIO::write_vector_ascii().

{
  return *_dof_map;
}

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.

{ return _equation_systems; }

const EquationSystems& System::get_equation_systems

(

)

const [inline, inherited]

Returns:

a constant reference to this system's parent EquationSystems object.

Definition at line 378 of file system.h.

References System::_equation_systems.

Referenced by NewmarkSystem::clear(), FrequencySystem::clear_all(), ExactErrorEstimator::find_squared_element_error(), ImplicitSystem::get_linear_solve_parameters(), FrequencySystem::init_data(), FrequencySystem::n_frequencies(), FrequencySystem::set_current_frequency(), FrequencySystem::set_frequencies(), FrequencySystem::set_frequencies_by_range(), FrequencySystem::set_frequencies_by_steps(), NewmarkSystem::set_newmark_parameters(), NonlinearImplicitSystem::set_solver_parameters(), LinearImplicitSystem::solve(), FrequencySystem::solve(), and EigenSystem::solve().

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

{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)

  std::ostringstream out;
  
  out << '\n'
      << " ---------------------------------------------------------------------------- \n"
      << "| Reference count information                                                |\n"
      << " ---------------------------------------------------------------------------- \n";
  
  for (Counts::iterator it = _counts.begin();
       it != _counts.end(); ++it)
    {
      const std::string name(it->first);
      const unsigned int creations    = it->second.first;
      const unsigned int destructions = it->second.second;

      out << "| " << name << " reference count information:\n"
          << "|  Creations:    " << creations    << '\n'
          << "|  Destructions: " << destructions << '\n';
    }
  
  out << " ---------------------------------------------------------------------------- \n";

  return out.str();

#else

  return "";
  
#endif
}

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

{
  std::ostringstream out;

  
  const std::string& sys_name = this->name();
      
  out << "   System \"" << sys_name << "\"\n"
      << "    Type \""  << this->system_type() << "\"\n"
      << "    Variables=";
  
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
      out << "\"" << this->variable_name(vn) << "\" ";
     
  out << '\n';

  out << "    Finite Element Types=";
#ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
    out << "\""
        << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_type(vn).family)
        << "\" ";  
#else
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
    {
      out << "\""
          << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_type(vn).family)
          << "\", \""
          << Utility::enum_to_string<FEFamily>(this->get_dof_map().variable_type(vn).radial_family)
          << "\" ";
    }

  out << '\n' << "    Infinite Element Mapping=";
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
    out << "\""
        << Utility::enum_to_string<InfMapType>(this->get_dof_map().variable_type(vn).inf_map)
        << "\" ";
#endif      

  out << '\n';
      
  out << "    Approximation Orders=";
  for (unsigned int vn=0; vn<this->n_vars(); vn++)
    {
#ifndef LIBMESH_ENABLE_INFINITE_ELEMENTS
      out << "\""
          << Utility::enum_to_string<Order>(this->get_dof_map().variable_type(vn).order)
          << "\" ";
#else
      out << "\""
          << Utility::enum_to_string<Order>(this->get_dof_map().variable_type(vn).order)
          << "\", \""
          << Utility::enum_to_string<Order>(this->get_dof_map().variable_type(vn).radial_order)
          << "\" ";
#endif
    }

  out << '\n';
      
  out << "    n_dofs()="             << this->n_dofs()             << '\n';
  out << "    n_local_dofs()="       << this->n_local_dofs()       << '\n';
#ifdef LIBMESH_ENABLE_AMR
  out << "    n_constrained_dofs()=" << this->n_constrained_dofs() << '\n';
#endif

  out << "    " << "n_vectors()="  << this->n_vectors()  << '\n';
//   out << "    " << "n_additional_matrices()=" << this->n_additional_matrices() << '\n';
  
  return out.str();
}

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.

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

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.

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

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.

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

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

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

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.

{
  return _mesh;
}

const MeshBase & System::get_mesh

(

)

const [inline, inherited]

Returns:

a constant reference to this systems's _mesh.

Definition at line 1334 of file system.h.

References System::_mesh.

Referenced by ExactSolution::_compute_error(), HPCoarsenTest::add_projection(), FEMSystem::assemble_qoi(), FEMSystem::assemble_qoi_derivative(), FEMSystem::assembly(), System::calculate_norm(), PatchRecoveryErrorEstimator::estimate_error(), JumpErrorEstimator::estimate_error(), ExactErrorEstimator::estimate_error(), AdjointResidualErrorEstimator::estimate_error(), System::init_data(), EigenSystem::init_data(), ImplicitSystem::init_matrices(), System::local_dof_indices(), DofMap::max_constraint_error(), FEMSystem::mesh_position_get(), FEMSystem::mesh_position_set(), System::ProjectVector::operator()(), PatchRecoveryErrorEstimator::EstimateError::operator()(), FEMSystem::postprocess(), System::project_vector(), System::read_header(), System::read_legacy_data(), System::read_parallel_data(), System::read_serialized_vector(), HPSingularity::select_refinement(), HPCoarsenTest::select_refinement(), System::write_header(), System::write_parallel_data(), System::write_serialized_vector(), and System::zero_variable().

{
  return _mesh;
}

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

{
  OStringStream sensitivity_rhs_name;
  sensitivity_rhs_name << "sensitivity_rhs" << i;

  return this->get_vector(sensitivity_rhs_name.str());
}

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

{
  OStringStream sensitivity_rhs_name;
  sensitivity_rhs_name << "sensitivity_rhs" << i;

  return this->get_vector(sensitivity_rhs_name.str());
}

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

{
  OStringStream sensitivity_name;
  sensitivity_name << "sensitivity_solution" << i;

  return this->get_vector(sensitivity_name.str());
}

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

{
  OStringStream sensitivity_name;
  sensitivity_name << "sensitivity_solution" << i;

  return this->get_vector(sensitivity_name.str());
}

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

{
  vectors_iterator v = vectors_begin();
  vectors_iterator v_end = vectors_end();
  unsigned int num = 0;
  while((num<vec_num) && (v!=v_end))
    {
      num++;
      ++v;
    }
  libmesh_assert(v!=v_end);
  return *(v->second);
}

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

{
  const_vectors_iterator v = vectors_begin();
  const_vectors_iterator v_end = vectors_end();
  unsigned int num = 0;
  while((num<vec_num) && (v!=v_end))
    {
      num++;
      ++v;
    }
  libmesh_assert(v!=v_end);
  return *(v->second);
}

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.

{
  // Make sure the vector exists
  vectors_iterator pos = _vectors.find(vec_name);
  
  if (pos == _vectors.end())
    {
      std::cerr << "ERROR: vector "
                << vec_name
                << " does not exist in this system!"
                << std::endl;      
      libmesh_error();
    }
  
  return *(pos->second);
}

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

{
  // Make sure the vector exists
  const_vectors_iterator pos = _vectors.find(vec_name);
  
  if (pos == _vectors.end())
    {
      std::cerr << "ERROR: vector "
                << vec_name
                << " does not exist in this system!"
                << std::endl;      
      libmesh_error();
    }
  
  return *(pos->second);
}

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

{
  OStringStream adjoint_name;
  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;

  return this->get_vector(adjoint_name.str());
}

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

{
  OStringStream adjoint_name;
  adjoint_name << "weighted_sensitivity_adjoint_solution" << i;

  return this->get_vector(adjoint_name.str());
}

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

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

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

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

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

{
  return _variable_numbers.count(var);
}

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

{
  return (_matrices.count(mat_name));
}

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

{
  return (_vectors.count(vec_name));
}

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

{
  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
  std::pair<unsigned int, unsigned int>& p = _counts[name];

  p.first++;
}

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

{
  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
  std::pair<unsigned int, unsigned int>& p = _counts[name];

  p.second++;
}

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

{ 
  // First initialize any required data
  this->init_data();
  
  //If no variables have been added to this system
  //don't do anything
  if(!this->n_vars())
    return;
 
  // Then call the user-provided intialization function
  this->user_initialization();
}

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

{
  Parent::init_context(c);

  FEMContext &context = libmesh_cast_ref<FEMContext&>(c);

  // Make sure we're prepared to do mass integration
  for (unsigned int var = 0; var != this->n_vars(); ++var)
    if (_time_evolving[var])
      {
        context.element_fe_var[var]->get_JxW();
        context.element_fe_var[var]->get_phi();
      }
}

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.

{
  // Add a vector which stores the tangent "du/ds" to the system and save its pointer.
  du_ds = &(add_vector("du_ds"));

  // Add a vector which stores the tangent "du/ds" to the system and save its pointer.
  previous_du_ds = &(add_vector("previous_du_ds"));

  // Add a vector to keep track of the previous nonlinear solution
  // at the old value of lambda.
  previous_u = &(add_vector("previous_u"));

  // Add a vector to keep track of the temporary solution "y" of Ay=G_{\lambda}.
  y = &(add_vector("y"));

  // Add a vector to keep track of the "old value" of "y" which is the solution of Ay=G_{\lambda}.
  y_old = &(add_vector("y_old"));

  // Add a vector to keep track of the temporary solution "z" of Az=-G.
  z = &(add_vector("z"));

  // Add a vector to keep track of the Newton update during the constrained PDE solves.
  delta_u = &(add_vector("delta_u"));

  // Call the Parent's initialization routine.
  Parent::init_data();
}

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

{
  libmesh_assert (matrix != NULL);

  // Check for quick return in case the system matrix
  // (and by extension all the matrices) has already
  // been initialized
  if (matrix->initialized())
    return;

  // Get a reference to the DofMap
  DofMap& dof_map = this->get_dof_map();
  
  // no chance to add other matrices
  _can_add_matrices = false;
  
  // Tell the matrices about the dof map, and vice versa
  for (matrices_iterator pos = _matrices.begin();
       pos != _matrices.end(); ++pos)
    {
      libmesh_assert (!pos->second->initialized());
      dof_map.attach_matrix (*(pos->second));
    }
  
  // Compute the sparsity pattern for the current
  // mesh and DOF distribution.  This also updates
  // additional matrices, \p DofMap now knows them
  dof_map.compute_sparsity (this->get_mesh());
  
  // Initialize matrices
  for (matrices_iterator pos = _matrices.begin(); 
       pos != _matrices.end(); ++pos)
    pos->second->init ();
  
  // Set the additional matrices to 0.
  for (matrices_iterator pos = _matrices.begin(); 
       pos != _matrices.end(); ++pos)
    pos->second->zero ();
}

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

{
  // Be sure the tangent was not already initialized.
  libmesh_assert (!tangent_initialized);
  
  // Compute delta_s_zero, the initial arclength travelled during the
  // first step.  Here we assume that previous_u and lambda_old store
  // the previous solution and control parameter.  You may need to
  // read in an old solution (or solve the non-continuation system)
  // first and call save_current_solution() before getting here.
  
  // 1.) Compute delta_s_zero as ||u|| - ||u_old|| + ...
  // Compute norms of the current and previous solutions
//   Real norm_u          = solution->l2_norm();
//   Real norm_previous_u = previous_u->l2_norm();
  
//   if (!quiet)
//     {
//       std::cout << "norm_u=" << norm_u << std::endl;
//       std::cout << "norm_previous_u=" << norm_previous_u << std::endl;
//     }
  
//   if (norm_u == norm_previous_u)
//     {
//       std::cerr << "Warning, it appears u and previous_u are the "
//              << "same, are you sure this is correct?"
//              << "It's possible you forgot to set one or the other..."
//              << std::endl;
//     }
  
//   Real delta_s_zero = std::sqrt(
//                              (norm_u - norm_previous_u)*(norm_u - norm_previous_u) +
//                              (*continuation_parameter-old_continuation_parameter)*
//                              (*continuation_parameter-old_continuation_parameter)
//                              );

//   // 2.) Compute delta_s_zero as ||u -u_old|| + ...
//   *delta_u = *solution;
//   delta_u->add(-1., *previous_u);
//   delta_u->close();
//   Real norm_delta_u = delta_u->l2_norm();
//   Real norm_u          = solution->l2_norm();
//   Real norm_previous_u = previous_u->l2_norm();

//   // Scale norm_delta_u by the bigger of either norm_u or norm_previous_u
//   norm_delta_u /= std::max(norm_u, norm_previous_u);
  
//   if (!quiet)
//     {
//       std::cout << "norm_u=" << norm_u << std::endl;
//       std::cout << "norm_previous_u=" << norm_previous_u << std::endl;
//       //std::cout << "norm_delta_u=" << norm_delta_u << std::endl;
//       std::cout << "norm_delta_u/max(|u|,|u_old|)=" << norm_delta_u << std::endl;
//       std::cout << "|norm_u-norm_previous_u|=" << std::abs(norm_u - norm_previous_u) << std::endl;
//     }
  
//   const Real dlambda = *continuation_parameter-old_continuation_parameter;

//   if (!quiet)
//     std::cout << "dlambda=" << dlambda << std::endl;
  
//   Real delta_s_zero = std::sqrt(
//                              (norm_delta_u*norm_delta_u) +
//                              (dlambda*dlambda)
//                              );
  
//   if (!quiet)
//     std::cout << "delta_s_zero=" << delta_s_zero << std::endl;

  // 1.) + 2.)
//   // Now approximate the initial tangent d(lambda)/ds
//   this->dlambda_ds = (*continuation_parameter-old_continuation_parameter) / delta_s_zero;


//   // We can also approximate the deriv. wrt s by finite differences:
//   // du/ds = (u1 - u0) / delta_s_zero.
//   // FIXME: Use delta_u from above if we decide to keep that method.
//   *du_ds = *solution;
//   du_ds->add(-1., *previous_u);
//   du_ds->scale(1./delta_s_zero);
//   du_ds->close();


  // 3.) Treating (u-previous_u)/(lambda - lambda_old) as an approximation to du/d(lambda),
  // we follow the same technique as Carnes and Shadid.
//   const Real dlambda = *continuation_parameter-old_continuation_parameter;
//   libmesh_assert (dlambda > 0.);
  
//   // Use delta_u for temporary calculation of du/d(lambda)
//   *delta_u = *solution;
//   delta_u->add(-1., *previous_u);
//   delta_u->scale(1. / dlambda);
//   delta_u->close();

//   // Determine initial normalization parameter
//   const Real solution_size = std::max(solution->l2_norm(), previous_u->l2_norm());
//   if (solution_size > 1.)
//     {
//       Theta = 1./solution_size;
      
//       if (!quiet)
//      std::cout << "Setting Normalization Parameter Theta=" << Theta << std::endl;
//     }
  
//   // Compute d(lambda)/ds
//   // The correct sign of d(lambda)/ds should be positive, since we assume that (lambda > lambda_old)
//   // but we could always double-check that as well.
//   Real norm_delta_u = delta_u->l2_norm();
//   this->dlambda_ds = 1. / std::sqrt(1. + Theta*Theta*norm_delta_u*norm_delta_u);

//   // Finally, compute du/ds = d(lambda)/ds * du/d(lambda)
//   *du_ds = *delta_u;
//   du_ds->scale(dlambda_ds);
//   du_ds->close();


  // 4.) Use normalized arclength formula to estimate delta_s_zero
//   // Determine initial normalization parameter
//   set_Theta();

//   // Compute (normalized) delta_s_zero
//   *delta_u = *solution;
//   delta_u->add(-1., *previous_u);
//   delta_u->close();
//   Real norm_delta_u = delta_u->l2_norm();
  
//   const Real dlambda = *