oomph::HopfHandler Class Reference

#include <assembly_handler.h>

Inheritance diagram for oomph::HopfHandler:

Public Member Functions
	HopfHandler (Problem *const &problem_pt, double *const &parameter_pt)
	Constructor. More...
	HopfHandler (Problem *const &problem_pt, double *const &paramter_pt, const double &omega, const DoubleVector &phi, const DoubleVector &psi)
	~HopfHandler ()
	Destructor return the problem to its original (non-augmented) state. More...
unsigned	ndof (GeneralisedElement *const &elem_pt)
	Get the number of elemental degrees of freedom. More...
unsigned long	eqn_number (GeneralisedElement *const &elem_pt, const unsigned &ieqn_local)
	Get the global equation number of the local unknown. More...
void	get_residuals (GeneralisedElement *const &elem_pt, Vector< double > &residuals)
	Get the residuals. More...
void	get_jacobian (GeneralisedElement *const &elem_pt, Vector< double > &residuals, DenseMatrix< double > &jacobian)
void	get_dresiduals_dparameter (GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam)
void	get_djacobian_dparameter (GeneralisedElement *const &elem_pt, double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
void	get_hessian_vector_products (GeneralisedElement *const &elem_pt, Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
int	bifurcation_type () const
double *	bifurcation_parameter_pt () const
void	get_eigenfunction (Vector< DoubleVector > &eigenfunction)
const double &	omega () const
	Return the frequency of the bifurcation. More...
void	solve_standard_system ()
	Set to solve the standard system. More...
void	solve_complex_system ()
	Set to solve the complex system. More...
void	solve_full_system ()
	Solve non-block system. More...
Public Member Functions inherited from oomph::AssemblyHandler
	AssemblyHandler ()
	Empty constructor. More...
virtual void	dof_vector (GeneralisedElement *const &elem_pt, const unsigned &t, Vector< double > &dof)
	Return vector of dofs at time level t in the element elem_pt. More...
virtual void	dof_pt_vector (GeneralisedElement *const &elem_pt, Vector< double * > &dof_pt)
	Return vector of pointers to dofs in the element elem_pt. More...
virtual double &	local_problem_dof (Problem *const &problem_pt, const unsigned &t, const unsigned &i)
virtual void	get_all_vectors_and_matrices (GeneralisedElement *const &elem_pt, Vector< Vector< double >> &vec, Vector< DenseMatrix< double >> &matrix)
virtual void	get_inner_products (GeneralisedElement *const &elem_pt, Vector< std::pair< unsigned, unsigned >> const &history_index, Vector< double > &inner_product)
virtual void	get_inner_product_vectors (GeneralisedElement *const &elem_pt, Vector< unsigned > const &history_index, Vector< Vector< double >> &inner_product_vector)
virtual	~AssemblyHandler ()
	Empty virtual destructor. More...

Private Attributes
unsigned	Solve_which_system
Problem *	Problem_pt
	Pointer to the problem. More...
double *	Parameter_pt
	Pointer to the parameter. More...
unsigned	Ndof
double	Omega
	The critical frequency of the bifurcation. More...
Vector< double >	Phi
	The real part of the null vector. More...
Vector< double >	Psi
	The imaginary part of the null vector. More...
Vector< double >	C
Vector< int >	Count

Friends
class	BlockHopfLinearSolver

Detailed Description

A class that is used to assemble the augmented system that defines a Hopf bifurcation. The "standard" problem must be a function of a global parameter \( \lambda \) and a solution is \( R(u,\lambda) = 0 \), where \( u \) are the unknowns in the problem. A Hopf bifurcation may be specified by the augmented system of size \( 3N+2 \).

\[ R(u,\lambda) = 0, \]

\[ J\phi + \omega M \psi = 0, \]

\[ J\psi - \omega M \phi = 0, \]

\[ c \cdot \phi = 1. \]

\[ c \cdot \psi = 0. \]

In the above \( J \) is the usual Jacobian matrix, \( dR/du \) and \( M \) is the mass matrix that multiplies the time derivative terms. \( \phi + i\psi \) is the (complex) null vector of the complex matrix \( J - i\omega M \), where \( \omega \) is the critical frequency. \( c \) is a constant vector that is used to ensure that the null vector is non-trivial.

Constructor & Destructor Documentation

HopfHandler() [1/2]

oomph::HopfHandler::HopfHandler	(	Problem *const &	problem_pt,
		double *const &	parameter_pt
	)

Constructor.

Constructor: Initialise the hopf handler, by setting initial guesses for Phi, Psi and calculating Count. If the system changes, a new handler must be constructed.

    : Solve_which_system(0), Parameter_pt(parameter_pt), Omega(0.0)
  {
    // Set the problem pointer
    Problem_pt = problem_pt;
    // Set the number of non-augmented degrees of freedom
    Ndof = problem_pt->ndof();
 
    // create the linear algebra distribution for this solver
    // currently only global (non-distributed) distributions are allowed
    LinearAlgebraDistribution* dist_pt =
      new LinearAlgebraDistribution(problem_pt->communicator_pt(), Ndof, false);
 
    // Resize the vectors of additional dofs
    Phi.resize(Ndof);
    Psi.resize(Ndof);
    C.resize(Ndof);
    Count.resize(Ndof, 0);
 
    // Loop over all the elements in the problem
    unsigned n_element = problem_pt->mesh_pt()->nelement();
    for (unsigned e = 0; e < n_element; e++)
    {
      GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
      // Loop over the local freedoms in an element
      unsigned n_var = elem_pt->ndof();
      for (unsigned n = 0; n < n_var; n++)
      {
        // Increase the associated global equation number counter
        ++Count[elem_pt->eqn_number(n)];
      }
    }
 
    // Calculate the value Phi by
    // solving the system JPhi = dF/dlambda
 
    // Locally cache the linear solver
    LinearSolver* const linear_solver_pt = problem_pt->linear_solver_pt();
 
    // Save the status before entry to this routine
    bool enable_resolve = linear_solver_pt->is_resolve_enabled();
 
    // We need to do a resolve
    linear_solver_pt->enable_resolve();
 
    // Storage for the solution
    DoubleVector x(dist_pt, 0.0);
 
    // Solve the standard problem, we only want to make sure that
    // we factorise the matrix, if it has not been factorised. We shall
    // ignore the return value of x.
    linear_solver_pt->solve(problem_pt, x);
 
    // Get the vector dresiduals/dparameter
    problem_pt->get_derivative_wrt_global_parameter(parameter_pt, x);
 
    // Copy rhs vector into local storage so it doesn't get overwritten
    // if the linear solver decides to initialise the solution vector, say,
    // which it's quite entitled to do!
    DoubleVector input_x(x);
 
    // Now resolve the system with the new RHS and overwrite the solution
    linear_solver_pt->resolve(input_x, x);
 
    // Restore the storage status of the linear solver
    if (enable_resolve)
    {
      linear_solver_pt->enable_resolve();
    }
    else
    {
      linear_solver_pt->disable_resolve();
    }
 
    // Normalise the solution x
    double length = 0.0;
    for (unsigned n = 0; n < Ndof; n++)
    {
      length += x[n] * x[n];
    }
    length = sqrt(length);
 
    // Now add the real part of the null space components to the problem
    // unknowns and initialise it all
    // This is dumb at the moment ... fix with eigensolver?
    for (unsigned n = 0; n < Ndof; n++)
    {
      problem_pt->Dof_pt.push_back(&Phi[n]);
      C[n] = Phi[n] = -x[n] / length;
    }
 
    // Set the imaginary part so that the appropriate residual is
    // zero initially (eigensolvers)
    for (unsigned n = 0; n < Ndof; n += 2)
    {
      // Make sure that we are not at the end of an array of odd length
      if (n != Ndof - 1)
      {
        Psi[n] = C[n + 1];
        Psi[n + 1] = -C[n];
      }
      // If it's odd set the final entry to zero
      else
      {
        Psi[n] = 0.0;
      }
    }
 
    // Next add the imaginary parts of the null space components to the problem
    for (unsigned n = 0; n < Ndof; n++)
    {
      problem_pt->Dof_pt.push_back(&Psi[n]);
    }
    // Now add the parameter
    problem_pt->Dof_pt.push_back(parameter_pt);
    // Finally add the frequency
    problem_pt->Dof_pt.push_back(&Omega);
 
    // rebuild the dof dist
    Problem_pt->Dof_distribution_pt->build(
      Problem_pt->communicator_pt(), Ndof * 3 + 2, false);
    // Remove all previous sparse storage used during Jacobian assembly
    Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
 
    // delete the dist_pt
    delete dist_pt;
  }

References oomph::LinearAlgebraDistribution::build(), oomph::Problem::communicator_pt(), Count, oomph::LinearSolver::disable_resolve(), oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, e(), oomph::Mesh::element_pt(), oomph::LinearSolver::enable_resolve(), oomph::GeneralisedElement::eqn_number(), oomph::Problem::get_derivative_wrt_global_parameter(), oomph::LinearSolver::is_resolve_enabled(), oomph::Problem::linear_solver_pt(), oomph::Problem::mesh_pt(), n, Ndof, oomph::GeneralisedElement::ndof(), oomph::Problem::ndof(), oomph::Mesh::nelement(), Omega, Phi, Problem_pt, Psi, oomph::LinearSolver::resolve(), oomph::LinearSolver::solve(), oomph::Problem::Sparse_assemble_with_arrays_previous_allocation, sqrt(), and plotDoE::x.

HopfHandler() [2/2]

oomph::HopfHandler::HopfHandler	(	Problem *const &	problem_pt,
		double *const &	parameter_pt,
		const double &	omega,
		const DoubleVector &	phi,
		const DoubleVector &	psi
	)

Constructor with initial guesses for the frequency and null vectors, such as might be provided by an eigensolver

Constructor: Initialise the hopf handler, by setting initial guesses for Phi, Psi, Omega and calculating Count. If the system changes, a new handler must be constructed.

    : Solve_which_system(0), Parameter_pt(parameter_pt), Omega(omega)
  {
    // Set the problem pointer
    Problem_pt = problem_pt;
    // Set the number of non-augmented degrees of freedom
    Ndof = problem_pt->ndof();
 
    // Resize the vectors of additional dofs
    Phi.resize(Ndof);
    Psi.resize(Ndof);
    C.resize(Ndof);
    Count.resize(Ndof, 0);
 
    // Loop over all the elements in the problem
    unsigned n_element = problem_pt->mesh_pt()->nelement();
    for (unsigned e = 0; e < n_element; e++)
    {
      GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
      // Loop over the local freedoms in an element
      unsigned n_var = elem_pt->ndof();
      for (unsigned n = 0; n < n_var; n++)
      {
        // Increase the associated global equation number counter
        ++Count[elem_pt->eqn_number(n)];
      }
    }
 
    // Normalise the guess for phi
    double length = 0.0;
    for (unsigned n = 0; n < Ndof; n++)
    {
      length += phi[n] * phi[n];
    }
    length = sqrt(length);
 
    // Now add the real part of the null space components to the problem
    // unknowns and initialise it all
    for (unsigned n = 0; n < Ndof; n++)
    {
      problem_pt->Dof_pt.push_back(&Phi[n]);
      C[n] = Phi[n] = phi[n] / length;
      Psi[n] = psi[n] / length;
    }
 
    // Next add the imaginary parts of the null space components to the problem
    for (unsigned n = 0; n < Ndof; n++)
    {
      problem_pt->Dof_pt.push_back(&Psi[n]);
    }
 
    // Now add the parameter
    problem_pt->Dof_pt.push_back(parameter_pt);
    // Finally add the frequency
    problem_pt->Dof_pt.push_back(&Omega);
 
    // rebuild the Dof_distribution_pt
    Problem_pt->Dof_distribution_pt->build(
      Problem_pt->communicator_pt(), Ndof * 3 + 2, false);
    // Remove all previous sparse storage used during Jacobian assembly
    Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
  }

References oomph::LinearAlgebraDistribution::build(), oomph::Problem::communicator_pt(), Count, oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, e(), oomph::Mesh::element_pt(), oomph::GeneralisedElement::eqn_number(), oomph::Problem::mesh_pt(), n, Ndof, oomph::GeneralisedElement::ndof(), oomph::Problem::ndof(), oomph::Mesh::nelement(), Omega, Phi, Problem_pt, Psi, oomph::Problem::Sparse_assemble_with_arrays_previous_allocation, and sqrt().

~HopfHandler()

oomph::HopfHandler::~HopfHandler

(

)

Destructor return the problem to its original (non-augmented) state.

Destructor, return the problem to its original state, before the augmented system was added

  {
    // If we are using the block solver reset the problem's linear solver
    // to the original one
    BlockHopfLinearSolver* block_hopf_solver_pt =
      dynamic_cast<BlockHopfLinearSolver*>(Problem_pt->linear_solver_pt());
    if (block_hopf_solver_pt)
    {
      // Reset the problem's linear solver
      Problem_pt->linear_solver_pt() = block_hopf_solver_pt->linear_solver_pt();
      // Delete the block solver
      delete block_hopf_solver_pt;
    }
    // Now return the problem to its original size
    Problem_pt->Dof_pt.resize(Ndof);
    Problem_pt->Dof_distribution_pt->build(
      Problem_pt->communicator_pt(), Ndof, false);
    // Remove all previous sparse storage used during Jacobian assembly
    Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
  }

References oomph::LinearAlgebraDistribution::build(), oomph::Problem::communicator_pt(), oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, oomph::Problem::linear_solver_pt(), oomph::BlockHopfLinearSolver::linear_solver_pt(), Ndof, Problem_pt, and oomph::Problem::Sparse_assemble_with_arrays_previous_allocation.

Member Function Documentation

bifurcation_parameter_pt()

double* oomph::HopfHandler::bifurcation_parameter_pt ( ) const

inlinevirtual

Return a pointer to the bifurcation parameter in bifurcation tracking problems

Reimplemented from oomph::AssemblyHandler.

    {
      return Parameter_pt;
    }

References Parameter_pt.

bifurcation_type()

int oomph::HopfHandler::bifurcation_type ( ) const

inlinevirtual

Indicate that we are tracking a Hopf bifurcation by returning 3

Reimplemented from oomph::AssemblyHandler.

    {
      return 3;
    }

eqn_number()

unsigned long oomph::HopfHandler::eqn_number ( GeneralisedElement *const &  elem_pt, const unsigned &  ieqn_local  )

virtual

Get the global equation number of the local unknown.

Reimplemented from oomph::AssemblyHandler.

  {
    // Get the raw value
    unsigned raw_ndof = elem_pt->ndof();
    unsigned long global_eqn;
    if (ieqn_local < raw_ndof)
    {
      global_eqn = elem_pt->eqn_number(ieqn_local);
    }
    else if (ieqn_local < 2 * raw_ndof)
    {
      global_eqn = Ndof + elem_pt->eqn_number(ieqn_local - raw_ndof);
    }
    else if (ieqn_local < 3 * raw_ndof)
    {
      global_eqn = 2 * Ndof + elem_pt->eqn_number(ieqn_local - 2 * raw_ndof);
    }
    else if (ieqn_local == 3 * raw_ndof)
    {
      global_eqn = 3 * Ndof;
    }
    else
    {
      global_eqn = 3 * Ndof + 1;
    }
    return global_eqn;
  }

References oomph::GeneralisedElement::eqn_number(), Ndof, and oomph::GeneralisedElement::ndof().

Referenced by get_jacobian().

get_djacobian_dparameter()

void oomph::HopfHandler::get_djacobian_dparameter ( GeneralisedElement *const &  elem_pt, double *const &  parameter_pt, Vector< double > &  dres_dparam, DenseMatrix< double > &  djac_dparam  )

virtual

Overload the derivative of the residuals and jacobian with respect to a parameter so that it breaks

Overload the derivative of the residuals and jacobian with respect to a parameter so that it breaks because it should not be required

Reimplemented from oomph::AssemblyHandler.

  {
    std::ostringstream error_stream;
    error_stream
      << "This function has not been implemented because it is not required\n";
    error_stream << "in standard problems.\n";
    error_stream
      << "If you find that you need it, you will have to implement it!\n\n";
 
    throw OomphLibError(
      error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
  }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

get_dresiduals_dparameter()

void oomph::HopfHandler::get_dresiduals_dparameter ( GeneralisedElement *const &  elem_pt, double *const &  parameter_pt, Vector< double > &  dres_dparam  )

virtual

Overload the derivatives of the residuals with respect to a parameter to apply to the augmented system

Get the derivatives of the augmented residuals with respect to a parameter

Reimplemented from oomph::AssemblyHandler.

  {
    // Should only call get residuals for the full system
    if (Solve_which_system == 0)
    {
      // Need to get raw residuals and jacobian
      unsigned raw_ndof = elem_pt->ndof();
 
      DenseMatrix<double> djac_dparam(raw_ndof), dM_dparam(raw_ndof);
      // Get the basic residuals, jacobian and mass matrix
      elem_pt->get_djacobian_and_dmass_matrix_dparameter(
        parameter_pt, dres_dparam, djac_dparam, dM_dparam);
 
      // Initialise the pen-ultimate residual, which does not
      // depend on the parameter
      dres_dparam[3 * raw_ndof] = 0.0;
      dres_dparam[3 * raw_ndof + 1] = 0.0;
 
      // Now multiply to fill in the residuals
      for (unsigned i = 0; i < raw_ndof; i++)
      {
        dres_dparam[raw_ndof + i] = 0.0;
        dres_dparam[2 * raw_ndof + i] = 0.0;
        for (unsigned j = 0; j < raw_ndof; j++)
        {
          unsigned global_unknown = elem_pt->eqn_number(j);
          // Real part
          dres_dparam[raw_ndof + i] +=
            djac_dparam(i, j) * Phi[global_unknown] +
            Omega * dM_dparam(i, j) * Psi[global_unknown];
          // Imaginary part
          dres_dparam[2 * raw_ndof + i] +=
            djac_dparam(i, j) * Psi[global_unknown] -
            Omega * dM_dparam(i, j) * Phi[global_unknown];
        }
      }
    }
    else
    {
      throw OomphLibError("Solve_which_system can only be 0",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
  }

References oomph::GeneralisedElement::eqn_number(), oomph::GeneralisedElement::get_djacobian_and_dmass_matrix_dparameter(), i, j, oomph::GeneralisedElement::ndof(), Omega, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Phi, Psi, and Solve_which_system.

get_eigenfunction()

void oomph::HopfHandler::get_eigenfunction ( Vector< DoubleVector > &  eigenfunction )

virtual

Return the eigenfunction(s) associated with the bifurcation that has been detected in bifurcation tracking problems

Reimplemented from oomph::AssemblyHandler.

  {
    // There is a real and imaginary part of the null vector
    eigenfunction.resize(2);
    LinearAlgebraDistribution dist(Problem_pt->communicator_pt(), Ndof, false);
    // Rebuild the vector
    eigenfunction[0].build(&dist, 0.0);
    eigenfunction[1].build(&dist, 0.0);
    // Set the value to be the null vector
    for (unsigned n = 0; n < Ndof; n++)
    {
      eigenfunction[0][n] = Phi[n];
      eigenfunction[1][n] = Psi[n];
    }
  }

References oomph::Problem::communicator_pt(), n, Ndof, Phi, Problem_pt, and Psi.

get_hessian_vector_products()

void oomph::HopfHandler::get_hessian_vector_products ( GeneralisedElement *const &  elem_pt, Vector< double > const &  Y, DenseMatrix< double > const &  C, DenseMatrix< double > &  product  )

virtual

Overload the hessian vector product function so that it breaks

Overload the hessian vector product function so that it breaks because it should not be required

Reimplemented from oomph::AssemblyHandler.

  {
    std::ostringstream error_stream;
    error_stream
      << "This function has not been implemented because it is not required\n";
    error_stream << "in standard problems.\n";
    error_stream
      << "If you find that you need it, you will have to implement it!\n\n";
 
    throw OomphLibError(
      error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
  }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

get_jacobian()

void oomph::HopfHandler::get_jacobian ( GeneralisedElement *const &  elem_pt, Vector< double > &  residuals, DenseMatrix< double > &  jacobian  )

virtual

Calculate the elemental Jacobian matrix "d equation / d variable".

Reimplemented from oomph::AssemblyHandler.

  {
    // The standard case
    if (Solve_which_system == 0)
    {
      unsigned augmented_ndof = ndof(elem_pt);
      unsigned raw_ndof = elem_pt->ndof();
 
      // Get the basic residuals and jacobian
      DenseMatrix<double> M(raw_ndof);
      elem_pt->get_jacobian_and_mass_matrix(residuals, jacobian, M);
      // Now fill in the actual residuals
      get_residuals(elem_pt, residuals);
 
      // Now the jacobian appears in other entries
      for (unsigned n = 0; n < raw_ndof; ++n)
      {
        for (unsigned m = 0; m < raw_ndof; ++m)
        {
          jacobian(raw_ndof + n, raw_ndof + m) = jacobian(n, m);
          jacobian(raw_ndof + n, 2 * raw_ndof + m) = Omega * M(n, m);
          jacobian(2 * raw_ndof + n, 2 * raw_ndof + m) = jacobian(n, m);
          jacobian(2 * raw_ndof + n, raw_ndof + m) = -Omega * M(n, m);
          unsigned global_eqn = elem_pt->eqn_number(m);
          jacobian(raw_ndof + n, 3 * raw_ndof + 1) += M(n, m) * Psi[global_eqn];
          jacobian(2 * raw_ndof + n, 3 * raw_ndof + 1) -=
            M(n, m) * Phi[global_eqn];
        }
 
        unsigned local_eqn = elem_pt->eqn_number(n);
        jacobian(3 * raw_ndof, raw_ndof + n) = C[local_eqn] / Count[local_eqn];
        jacobian(3 * raw_ndof + 1, 2 * raw_ndof + n) =
          C[local_eqn] / Count[local_eqn];
      }
 
      const double FD_step = 1.0e-8;
 
      Vector<double> newres_p(augmented_ndof), newres_m(augmented_ndof);
 
      // Loop over the dofs
      for (unsigned n = 0; n < raw_ndof; n++)
      {
        // Just do the x's
        unsigned long global_eqn = eqn_number(elem_pt, n);
        double* unknown_pt = Problem_pt->Dof_pt[global_eqn];
        double init = *unknown_pt;
        *unknown_pt += FD_step;
 
        // Get the new residuals
        get_residuals(elem_pt, newres_p);
 
        // Reset
        *unknown_pt = init;
 
        // Subtract
        *unknown_pt -= FD_step;
        get_residuals(elem_pt, newres_m);
 
        for (unsigned m = 0; m < raw_ndof; m++)
        {
          jacobian(raw_ndof + m, n) =
            (newres_p[raw_ndof + m] - residuals[raw_ndof + m]) / (FD_step);
          jacobian(2 * raw_ndof + m, n) =
            (newres_p[2 * raw_ndof + m] - residuals[2 * raw_ndof + m]) /
            (FD_step);
        }
        // Reset the unknown
        *unknown_pt = init;
      }
 
      {
        // Now do the global parameter
        double* unknown_pt = Problem_pt->Dof_pt[3 * Ndof];
        double init = *unknown_pt;
        *unknown_pt += FD_step;
 
        Problem_pt->actions_after_change_in_bifurcation_parameter();
        // Get the new residuals
        get_residuals(elem_pt, newres_p);
 
        // Reset
        *unknown_pt = init;
 
        // Subtract
        *unknown_pt -= FD_step;
        get_residuals(elem_pt, newres_m);
 
        for (unsigned m = 0; m < augmented_ndof - 2; m++)
        {
          jacobian(m, 3 * raw_ndof) = (newres_p[m] - residuals[m]) / FD_step;
        }
        // Reset the unknown
        *unknown_pt = init;
        Problem_pt->actions_after_change_in_bifurcation_parameter();
      }
    } // End of standard case
    // Normal case
    else if (Solve_which_system == 1)
    {
      // Just get the normal jacobian and residuals
      elem_pt->get_jacobian(residuals, jacobian);
    }
    // Otherwise the augmented complex case
    else if (Solve_which_system == 2)
    {
      unsigned raw_ndof = elem_pt->ndof();
 
      // Get the basic residuals and jacobian
      DenseMatrix<double> M(raw_ndof);
      elem_pt->get_jacobian_and_mass_matrix(residuals, jacobian, M);
 
      // We now need to fill in the other blocks
      for (unsigned n = 0; n < raw_ndof; n++)
      {
        for (unsigned m = 0; m < raw_ndof; m++)
        {
          jacobian(n, raw_ndof + m) = Omega * M(n, m);
          jacobian(raw_ndof + n, m) = -Omega * M(n, m);
          jacobian(raw_ndof + n, raw_ndof + m) = jacobian(n, m);
        }
      }
 
      // Now overwrite to fill in the residuals
      // The decision take is to solve for the mass matrix multiplied
      // terms in the residuals because they require no additional
      // information to assemble.
      for (unsigned n = 0; n < raw_ndof; n++)
      {
        residuals[n] = 0.0;
        residuals[raw_ndof + n] = 0.0;
        for (unsigned m = 0; m < raw_ndof; m++)
        {
          unsigned global_unknown = elem_pt->eqn_number(m);
          // Real part
          residuals[n] += M(n, m) * Psi[global_unknown];
          // Imaginary part
          residuals[raw_ndof + n] -= M(n, m) * Phi[global_unknown];
        }
      }
    } // End of complex augmented case
    else
    {
      throw OomphLibError("Solve_which_system can only be 0,1 or 2",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
  }

References oomph::Problem::actions_after_change_in_bifurcation_parameter(), Count, oomph::Problem::Dof_pt, oomph::GeneralisedElement::eqn_number(), eqn_number(), oomph::BlackBoxFDNewtonSolver::FD_step, oomph::GeneralisedElement::get_jacobian(), oomph::GeneralisedElement::get_jacobian_and_mass_matrix(), get_residuals(), m, n, Ndof, oomph::GeneralisedElement::ndof(), ndof(), Omega, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Phi, Problem_pt, Psi, and Solve_which_system.

get_residuals()

void oomph::HopfHandler::get_residuals ( GeneralisedElement *const &  elem_pt, Vector< double > &  residuals  )

virtual

Get the residuals.

Reimplemented from oomph::AssemblyHandler.

  {
    // Should only call get residuals for the full system
    if (Solve_which_system == 0)
    {
      // Need to get raw residuals and jacobian
      unsigned raw_ndof = elem_pt->ndof();
 
      DenseMatrix<double> jacobian(raw_ndof), M(raw_ndof);
      // Get the basic residuals, jacobian and mass matrix
      elem_pt->get_jacobian_and_mass_matrix(residuals, jacobian, M);
 
      // Initialise the pen-ultimate residual
      residuals[3 * raw_ndof] =
        -1.0 / (double)(Problem_pt->mesh_pt()->nelement());
      residuals[3 * raw_ndof + 1] = 0.0;
 
      // Now multiply to fill in the residuals
      for (unsigned i = 0; i < raw_ndof; i++)
      {
        residuals[raw_ndof + i] = 0.0;
        residuals[2 * raw_ndof + i] = 0.0;
        for (unsigned j = 0; j < raw_ndof; j++)
        {
          unsigned global_unknown = elem_pt->eqn_number(j);
          // Real part
          residuals[raw_ndof + i] += jacobian(i, j) * Phi[global_unknown] +
                                     Omega * M(i, j) * Psi[global_unknown];
          // Imaginary part
          residuals[2 * raw_ndof + i] += jacobian(i, j) * Psi[global_unknown] -
                                         Omega * M(i, j) * Phi[global_unknown];
        }
        // Get the global equation number
        unsigned global_eqn = elem_pt->eqn_number(i);
 
        // Real part
        residuals[3 * raw_ndof] +=
          (Phi[global_eqn] * C[global_eqn]) / Count[global_eqn];
        // Imaginary part
        residuals[3 * raw_ndof + 1] +=
          (Psi[global_eqn] * C[global_eqn]) / Count[global_eqn];
      }
    }
    else
    {
      throw OomphLibError("Solve_which_system can only be 0",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
  }

References Count, oomph::GeneralisedElement::eqn_number(), oomph::GeneralisedElement::get_jacobian_and_mass_matrix(), i, j, oomph::Problem::mesh_pt(), oomph::GeneralisedElement::ndof(), oomph::Mesh::nelement(), Omega, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Phi, Problem_pt, Psi, and Solve_which_system.

Referenced by get_jacobian().

ndof()

unsigned oomph::HopfHandler::ndof ( GeneralisedElement *const &  elem_pt )

virtual

Get the number of elemental degrees of freedom.

Reimplemented from oomph::AssemblyHandler.

  {
    unsigned raw_ndof = elem_pt->ndof();
    switch (Solve_which_system)
    {
        // Full augmented system
      case 0:
        return (3 * raw_ndof + 2);
        break;
        // Standard non-augmented system
      case 1:
        return raw_ndof;
        break;
        // Complex system
      case 2:
        return (2 * raw_ndof);
        break;
 
      default:
        throw OomphLibError("Solve_which_system can only be 0,1 or 2",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
    }
  }

References oomph::GeneralisedElement::ndof(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and Solve_which_system.

Referenced by get_jacobian().

omega()

const double& oomph::HopfHandler::omega ( ) const

inline

Return the frequency of the bifurcation.

    {
      return Omega;
    }

References Omega.

solve_complex_system()

void oomph::HopfHandler::solve_complex_system

(

)

Set to solve the complex system.

Set to solve the complex (jacobian and mass matrix) system.

  {
    // If we were not solving the complex system resize the unknowns
    // accordingly
    if (Solve_which_system != 2)
    {
      Solve_which_system = 2;
      // Resize to the first Ndofs (will work whichever system we were
      // solving before)
      Problem_pt->Dof_pt.resize(Ndof);
      // Add the first (real) part of the eigenfunction back into the problem
      for (unsigned n = 0; n < Ndof; n++)
      {
        Problem_pt->Dof_pt.push_back(&Phi[n]);
      }
      Problem_pt->Dof_distribution_pt->build(
        Problem_pt->communicator_pt(), Ndof * 2, false);
      // Remove all previous sparse storage used during Jacobian assembly
      Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
    }
  }

References oomph::LinearAlgebraDistribution::build(), oomph::Problem::communicator_pt(), oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, n, Ndof, Phi, Problem_pt, Solve_which_system, and oomph::Problem::Sparse_assemble_with_arrays_previous_allocation.

Referenced by oomph::BlockHopfLinearSolver::solve(), and oomph::BlockHopfLinearSolver::solve_for_two_rhs().

solve_full_system()

void oomph::HopfHandler::solve_full_system

(

)

Solve non-block system.

Set to Solve full system system.

  {
    // If we are starting from another system
    if (Solve_which_system)
    {
      Solve_which_system = 0;
 
      // Resize to the first Ndofs (will work whichever system we were
      // solving before)
      Problem_pt->Dof_pt.resize(Ndof);
      // Add the additional unknowns back into the problem
      for (unsigned n = 0; n < Ndof; n++)
      {
        Problem_pt->Dof_pt.push_back(&Phi[n]);
      }
      for (unsigned n = 0; n < Ndof; n++)
      {
        Problem_pt->Dof_pt.push_back(&Psi[n]);
      }
      // Now add the parameter
      Problem_pt->Dof_pt.push_back(Parameter_pt);
      // Finally add the frequency
      Problem_pt->Dof_pt.push_back(&Omega);
 
      //
      Problem_pt->Dof_distribution_pt->build(
        Problem_pt->communicator_pt(), 3 * Ndof + 2, false);
      // Remove all previous sparse storage used during Jacobian assembly
      Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
    }
  }

References oomph::LinearAlgebraDistribution::build(), oomph::Problem::communicator_pt(), oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, n, Ndof, Omega, Parameter_pt, Phi, Problem_pt, Psi, Solve_which_system, and oomph::Problem::Sparse_assemble_with_arrays_previous_allocation.

Referenced by oomph::BlockHopfLinearSolver::solve(), and oomph::BlockHopfLinearSolver::solve_for_two_rhs().

solve_standard_system()

void oomph::HopfHandler::solve_standard_system

(

)

Set to solve the standard system.

Set to solve the standard (underlying jacobian) system.

  {
    if (Solve_which_system != 1)
    {
      Solve_which_system = 1;
      // Restrict the problem to the standard variables only
      Problem_pt->Dof_pt.resize(Ndof);
      Problem_pt->Dof_distribution_pt->build(
        Problem_pt->communicator_pt(), Ndof, false);
      // Remove all previous sparse storage used during Jacobian assembly
      Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
    }
  }

References oomph::LinearAlgebraDistribution::build(), oomph::Problem::communicator_pt(), oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, Ndof, Problem_pt, Solve_which_system, and oomph::Problem::Sparse_assemble_with_arrays_previous_allocation.

Referenced by oomph::BlockHopfLinearSolver::solve(), and oomph::BlockHopfLinearSolver::solve_for_two_rhs().

Friends And Related Function Documentation

BlockHopfLinearSolver

friend class BlockHopfLinearSolver

friend

Member Data Documentation

C

Vector<double> oomph::HopfHandler::C

private

A constant vector used to ensure that the null vector is not trivial

Referenced by oomph::BlockHopfLinearSolver::solve(), and oomph::BlockHopfLinearSolver::solve_for_two_rhs().

Count

Vector<int> oomph::HopfHandler::Count

private

A vector that is used to determine how many elements contribute to a particular equation. It is used to ensure that the global system is correctly formulated.

Referenced by get_jacobian(), get_residuals(), and HopfHandler().

Ndof

unsigned oomph::HopfHandler::Ndof

private

Store the number of degrees of freedom in the non-augmented problem

Referenced by eqn_number(), get_eigenfunction(), get_jacobian(), HopfHandler(), solve_complex_system(), solve_full_system(), solve_standard_system(), and ~HopfHandler().

Omega

double oomph::HopfHandler::Omega

private

The critical frequency of the bifurcation.

Referenced by get_dresiduals_dparameter(), get_jacobian(), get_residuals(), HopfHandler(), omega(), oomph::BlockHopfLinearSolver::solve(), oomph::BlockHopfLinearSolver::solve_for_two_rhs(), and solve_full_system().

Parameter_pt

double* oomph::HopfHandler::Parameter_pt

private

Pointer to the parameter.

Referenced by bifurcation_parameter_pt(), and solve_full_system().

Phi

Vector<double> oomph::HopfHandler::Phi

private

The real part of the null vector.

Referenced by get_dresiduals_dparameter(), get_eigenfunction(), get_jacobian(), get_residuals(), HopfHandler(), oomph::BlockHopfLinearSolver::solve(), solve_complex_system(), oomph::BlockHopfLinearSolver::solve_for_two_rhs(), and solve_full_system().

Problem_pt

Problem* oomph::HopfHandler::Problem_pt

private

Pointer to the problem.

Referenced by get_eigenfunction(), get_jacobian(), get_residuals(), HopfHandler(), solve_complex_system(), solve_full_system(), solve_standard_system(), and ~HopfHandler().

Psi

Vector<double> oomph::HopfHandler::Psi

private

The imaginary part of the null vector.

Referenced by get_dresiduals_dparameter(), get_eigenfunction(), get_jacobian(), get_residuals(), HopfHandler(), oomph::BlockHopfLinearSolver::solve(), oomph::BlockHopfLinearSolver::solve_for_two_rhs(), and solve_full_system().

Solve_which_system

unsigned oomph::HopfHandler::Solve_which_system

private

Integer flag to indicate which system should be assembled. There are three possibilities. The full augmented system (0), the non-augmented jacobian system (1), and complex system (2), where the matrix is a combination of the jacobian and mass matrices.

Referenced by get_dresiduals_dparameter(), get_jacobian(), get_residuals(), ndof(), solve_complex_system(), solve_full_system(), and solve_standard_system().

---

The documentation for this class was generated from the following files:

• assembly_handler.h
• assembly_handler.cc