oomph::BlockHopfLinearSolver Class Reference

#include <assembly_handler.h>

Inheritance diagram for oomph::BlockHopfLinearSolver:

Public Member Functions
	BlockHopfLinearSolver (LinearSolver *const linear_solver_pt)
	Constructor, inherits the original linear solver. More...
	~BlockHopfLinearSolver ()
	Destructor: clean up the allocated memory. More...
void	solve_for_two_rhs (Problem *const &problem_pt, DoubleVector &result, const DoubleVector &rhs2, DoubleVector &result2)
	Solve for two right hand sides. More...
void	solve (Problem *const &problem_pt, DoubleVector &result)
	The solve function uses the block factorisation. More...
void	solve (DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)
void	solve (DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
void	resolve (const DoubleVector &rhs, DoubleVector &result)
	The resolve function also uses the block factorisation. More...
LinearSolver *	linear_solver_pt () const
	Access function to the original linear solver. More...
Public Member Functions inherited from oomph::LinearSolver
	LinearSolver ()
	Empty constructor, initialise the member data. More...
	LinearSolver (const LinearSolver &dummy)=delete
	Broken copy constructor. More...
void	operator= (const LinearSolver &)=delete
	Broken assignment operator. More...
virtual	~LinearSolver ()
	Empty virtual destructor. More...
void	enable_doc_time ()
	Enable documentation of solve times. More...
void	disable_doc_time ()
	Disable documentation of solve times. More...
bool	is_doc_time_enabled () const
	Is documentation of solve times enabled? More...
bool	is_resolve_enabled () const
	Boolean flag indicating if resolves are enabled. More...
virtual void	enable_resolve ()
virtual void	disable_resolve ()
virtual void	solve_transpose (Problem *const &problem_pt, DoubleVector &result)
virtual void	solve_transpose (DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)
virtual void	solve_transpose (DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
virtual void	resolve_transpose (const DoubleVector &rhs, DoubleVector &result)
virtual void	clean_up_memory ()
virtual double	jacobian_setup_time () const
virtual double	linear_solver_solution_time () const
virtual void	enable_computation_of_gradient ()
void	disable_computation_of_gradient ()
void	reset_gradient ()
void	get_gradient (DoubleVector &gradient)
	function to access the gradient, provided it has been computed More...
Public Member Functions inherited from oomph::DistributableLinearAlgebraObject
	DistributableLinearAlgebraObject ()
	Default constructor - create a distribution. More...
	DistributableLinearAlgebraObject (const DistributableLinearAlgebraObject &matrix)=delete
	Broken copy constructor. More...
void	operator= (const DistributableLinearAlgebraObject &)=delete
	Broken assignment operator. More...
virtual	~DistributableLinearAlgebraObject ()
	Destructor. More...
LinearAlgebraDistribution *	distribution_pt () const
	access to the LinearAlgebraDistribution More...
unsigned	nrow () const
	access function to the number of global rows. More...
unsigned	nrow_local () const
	access function for the num of local rows on this processor. More...
unsigned	nrow_local (const unsigned &p) const
	access function for the num of local rows on this processor. More...
unsigned	first_row () const
	access function for the first row on this processor More...
unsigned	first_row (const unsigned &p) const
	access function for the first row on this processor More...
bool	distributed () const
	distribution is serial or distributed More...
bool	distribution_built () const
void	build_distribution (const LinearAlgebraDistribution *const dist_pt)
void	build_distribution (const LinearAlgebraDistribution &dist)

Private Attributes
LinearSolver *	Linear_solver_pt
	Pointer to the original linear solver. More...
Problem *	Problem_pt
	Pointer to the problem, used in the resolve. More...
DoubleVector *	A_pt
	Pointer to the storage for the vector a. More...
DoubleVector *	E_pt
	Pointer to the storage for the vector e (0 to n-1) More...
DoubleVector *	G_pt
	Pointer to the storage for the vector g (0 to n-1) More...

Additional Inherited Members
Protected Member Functions inherited from oomph::DistributableLinearAlgebraObject
void	clear_distribution ()
Protected Attributes inherited from oomph::LinearSolver
bool	Enable_resolve
bool	Doc_time
	Boolean flag that indicates whether the time taken. More...
bool	Compute_gradient
bool	Gradient_has_been_computed
	flag that indicates whether the gradient was computed or not More...
DoubleVector	Gradient_for_glob_conv_newton_solve

Detailed Description

A custom linear solver class that is used to solve a block-factorised version of the Hopf bifurcation detection problem.

Constructor & Destructor Documentation

BlockHopfLinearSolver()

oomph::BlockHopfLinearSolver::BlockHopfLinearSolver ( LinearSolver *const  linear_solver_pt )

inline

Constructor, inherits the original linear solver.

      : Linear_solver_pt(linear_solver_pt),
        Problem_pt(0),
        A_pt(0),
        E_pt(0),
        G_pt(0)
    {
    }

~BlockHopfLinearSolver()

oomph::BlockHopfLinearSolver::~BlockHopfLinearSolver

(

)

Destructor: clean up the allocated memory.

Clean up the memory that may have been allocated by the solver.

  {
    if (A_pt != 0)
    {
      delete A_pt;
    }
    if (E_pt != 0)
    {
      delete E_pt;
    }
    if (G_pt != 0)
    {
      delete G_pt;
    }
  }

References A_pt, E_pt, and G_pt.

Member Function Documentation

linear_solver_pt()

LinearSolver* oomph::BlockHopfLinearSolver::linear_solver_pt ( ) const

inline

Access function to the original linear solver.

    {
      return Linear_solver_pt;
    }

References Linear_solver_pt.

Referenced by oomph::HopfHandler::~HopfHandler().

resolve()

void oomph::BlockHopfLinearSolver::resolve ( const DoubleVector &  rhs, DoubleVector &  result  )

virtual

The resolve function also uses the block factorisation.

Reimplemented from oomph::LinearSolver.

  {
    throw OomphLibError("resolve() is not implemented for this solver",
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
  }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

solve() [1/3]

void oomph::BlockHopfLinearSolver::solve ( DoubleMatrixBase *const &  matrix_pt, const DoubleVector &  rhs, DoubleVector &  result  )

inlinevirtual

The linear-algebra-type solver does not make sense. The interface is deliberately broken

Reimplemented from oomph::LinearSolver.

    {
      throw OomphLibError(
        "Linear-algebra interface does not make sense for this linear solver\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

solve() [2/3]

void oomph::BlockHopfLinearSolver::solve ( DoubleMatrixBase *const &  matrix_pt, const Vector< double > &  rhs, Vector< double > &  result  )

inlinevirtual

The linear-algebra-type solver does not make sense. The interface is deliberately broken

Reimplemented from oomph::LinearSolver.

    {
      throw OomphLibError(
        "Linear-algebra interface does not make sense for this linear solver\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

solve() [3/3]

void oomph::BlockHopfLinearSolver::solve ( Problem *const &  problem_pt, DoubleVector &  result  )

virtual

The solve function uses the block factorisation.

Use a block factorisation to solve the augmented system associated with a Hopf bifurcation.

Implements oomph::LinearSolver.

  {
    // if the result is setup then it should not be distributed
#ifdef PARANOID
    if (result.built())
    {
      if (result.distributed())
      {
        throw OomphLibError("The result vector must not be distributed",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // Locally cache the pointer to the handler.
    HopfHandler* handler_pt =
      static_cast<HopfHandler*>(problem_pt->assembly_handler_pt());
 
    // Find the number of dofs in the augmented problem
    unsigned n_dof = problem_pt->ndof();
 
    // create the linear algebra distribution for this solver
    // currently only global (non-distributed) distributions are allowed
    LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
    this->build_distribution(dist);
 
    // Firstly, let's calculate the derivative of the residuals wrt
    // the parameter
    DoubleVector dRdparam(this->distribution_pt(), 0.0);
 
    const double FD_step = 1.0e-8;
    {
      // Cache pointer to the parameter (second last entry in the vector
      double* param_pt = &problem_pt->dof(n_dof - 2);
      // Backup the parameter
      double old_var = *param_pt;
      ;
      // Increment the parameter
      *param_pt += FD_step;
      problem_pt->actions_after_change_in_bifurcation_parameter();
 
      // Now calculate the new residuals
      problem_pt->get_residuals(dRdparam);
 
      // Now calculate the difference assume original residuals in resul
      for (unsigned n = 0; n < n_dof; n++)
      {
        dRdparam[n] = (dRdparam[n] - result[n]) / FD_step;
      }
 
      // Reset the parameter
      *param_pt = old_var;
      problem_pt->actions_after_change_in_bifurcation_parameter();
    }
 
    // Switch the handler to "standard" mode
    handler_pt->solve_standard_system();
 
    // Find out the number of dofs in the non-standard problem
    n_dof = problem_pt->ndof();
 
    // update the distribution pt
    dist.build(problem_pt->communicator_pt(), n_dof, false);
    this->build_distribution(dist);
 
    // Setup storage for temporary vector
    DoubleVector y1(this->distribution_pt(), 0.0),
      alpha(this->distribution_pt(), 0.0);
 
    // Allocate storage for A which can be used in the resolve
    if (A_pt != 0)
    {
      delete A_pt;
    }
    A_pt = new DoubleVector(this->distribution_pt(), 0.0);
 
    // We are going to do resolves using the underlying linear solver
    Linear_solver_pt->enable_resolve();
 
    // Solve the first system Jy1 = R
    Linear_solver_pt->solve(problem_pt, y1);
 
    // This should have set the sign of the jacobian, store it
    int sign_of_jacobian = problem_pt->sign_of_jacobian();
 
    // Now let's get the appropriate bit of alpha
    for (unsigned n = 0; n < n_dof; n++)
    {
      alpha[n] = dRdparam[n];
    }
 
    // Resolve to find A
    Linear_solver_pt->resolve(alpha, *A_pt);
 
    // Now set to the complex system
    handler_pt->solve_complex_system();
 
    // update the distribution
    dist.build(problem_pt->communicator_pt(), n_dof * 2, false);
    this->build_distribution(dist);
 
    // Resize the stored vector G
    if (G_pt != 0)
    {
      delete G_pt;
    }
    G_pt = new DoubleVector(this->distribution_pt(), 0.0);
 
    // Solve the first Zg = (Mz -My)
    Linear_solver_pt->solve(problem_pt, *G_pt);
 
    // This should have set the sign of the complex matrix's determinant,
    // multiply
    sign_of_jacobian *= problem_pt->sign_of_jacobian();
 
    // We can now construct our multipliers
    // Prepare to scale
    double dof_length = 0.0, a_length = 0.0, y1_length = 0.0;
    // Loop over the standard number of dofs
    for (unsigned n = 0; n < n_dof; n++)
    {
      if (std::fabs(problem_pt->dof(n)) > dof_length)
      {
        dof_length = std::fabs(problem_pt->dof(n));
      }
      if (std::fabs((*A_pt)[n]) > a_length)
      {
        a_length = std::fabs((*A_pt)[n]);
      }
      if (std::fabs(y1[n]) > y1_length)
      {
        y1_length = std::fabs(y1[n]);
      }
    }
 
    double a_mult = dof_length / a_length;
    double y1_mult = dof_length / y1_length;
    a_mult += FD_step;
    y1_mult += FD_step;
    a_mult *= FD_step;
    y1_mult *= FD_step;
 
    // Local storage for the product terms
    Vector<double> Jprod_a(2 * n_dof, 0.0), Jprod_y1(2 * n_dof, 0.0);
 
    // Temporary storage
    Vector<double> rhs(2 * n_dof);
 
    // find the number of elements
    unsigned long n_element = problem_pt->mesh_pt()->nelement();
 
    // Calculate the product of the jacobian matrices, etc
    for (unsigned long e = 0; e < n_element; e++)
    {
      GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
      // Loop over the ndofs in each element
      unsigned n_var = elem_pt->ndof();
      // Get the jacobian matrices
      DenseMatrix<double> jac(n_var), jac_a(n_var), jac_y1(n_var);
      DenseMatrix<double> M(n_var), M_a(n_var), M_y1(n_var);
      // Get unperturbed jacobian
      elem_pt->get_jacobian_and_mass_matrix(rhs, jac, M);
 
      // Backup the dofs
      Vector<double> dof_bac(n_var);
      // Perturb the original dofs
      for (unsigned n = 0; n < n_var; n++)
      {
        unsigned eqn_number = elem_pt->eqn_number(n);
        dof_bac[n] = problem_pt->dof(eqn_number);
        // Pertub by vector A
        problem_pt->dof(eqn_number) += a_mult * (*A_pt)[eqn_number];
      }
 
      // Now get the new jacobian and mass matrix
      elem_pt->get_jacobian_and_mass_matrix(rhs, jac_a, M_a);
 
      // Perturb the dofs
      for (unsigned n = 0; n < n_var; n++)
      {
        unsigned eqn_number = elem_pt->eqn_number(n);
        problem_pt->dof(eqn_number) = dof_bac[n];
        // Pertub by vector y1
        problem_pt->dof(eqn_number) += y1_mult * y1[eqn_number];
      }
 
      // Now get the new jacobian and mass matrix
      elem_pt->get_jacobian_and_mass_matrix(rhs, jac_y1, M_y1);
 
      // Reset the dofs
      for (unsigned n = 0; n < n_var; n++)
      {
        unsigned eqn_number = elem_pt->eqn_number(n);
        problem_pt->dof(eqn_number) = dof_bac[n];
      }
 
      // OK, now work out the products
      for (unsigned n = 0; n < n_var; n++)
      {
        unsigned eqn_number = elem_pt->eqn_number(n);
        double prod_a1 = 0.0, prod_y11 = 0.0;
        double prod_a2 = 0.0, prod_y12 = 0.0;
        for (unsigned m = 0; m < n_var; m++)
        {
          const unsigned unknown = elem_pt->eqn_number(m);
          const double y = handler_pt->Phi[unknown];
          const double z = handler_pt->Psi[unknown];
          const double omega = handler_pt->Omega;
          // Real part (first line)
          prod_a1 +=
            (jac_a(n, m) - jac(n, m)) * y + omega * (M_a(n, m) - M(n, m)) * z;
          prod_y11 +=
            (jac_y1(n, m) - jac(n, m)) * y + omega * (M_y1(n, m) - M(n, m)) * z;
          // Imag par (second line)
          prod_a2 +=
            (jac_a(n, m) - jac(n, m)) * z - omega * (M_a(n, m) - M(n, m)) * y;
          prod_y12 +=
            (jac_y1(n, m) - jac(n, m)) * z - omega * (M_y1(n, m) - M(n, m)) * y;
        }
        Jprod_a[eqn_number] += prod_a1 / a_mult;
        Jprod_y1[eqn_number] += prod_y11 / y1_mult;
        Jprod_a[n_dof + eqn_number] += prod_a2 / a_mult;
        Jprod_y1[n_dof + eqn_number] += prod_y12 / y1_mult;
      }
    }
 
    // The assumption here is still that the result has been set to the
    // residuals.
    for (unsigned n = 0; n < 2 * n_dof; n++)
    {
      rhs[n] = result[n_dof + n] - Jprod_y1[n];
    }
 
    // Temporary storage
    DoubleVector y2(this->distribution_pt(), 0.0);
 
    // DoubleVector for rhs
    DoubleVector rhs2(this->distribution_pt(), 0.0);
    for (unsigned i = 0; i < 2 * n_dof; i++)
    {
      rhs2[i] = rhs[i];
    }
 
    // Solve it
    Linear_solver_pt->resolve(rhs2, y2);
 
    // Assemble the next RHS
    for (unsigned n = 0; n < 2 * n_dof; n++)
    {
      rhs[n] = dRdparam[n_dof + n] - Jprod_a[n];
    }
 
    // Resive the storage
    if (E_pt != 0)
    {
      delete E_pt;
    }
    E_pt = new DoubleVector(this->distribution_pt(), 0.0);
 
    // Solve for the next RHS
    for (unsigned i = 0; i < 2 * n_dof; i++)
    {
      rhs2[i] = rhs[i];
    }
    Linear_solver_pt->resolve(rhs2, *E_pt);
 
    // We can now calculate the final corrections
    // We need to work out a large number of dot products
    double dot_c = 0.0, dot_d = 0.0, dot_e = 0.0, dot_f = 0.0, dot_g = 0.0;
    double dot_h = 0.0;
 
    for (unsigned n = 0; n < n_dof; n++)
    {
      // Get the appopriate entry
      const double Cn = handler_pt->C[n];
      dot_c += Cn * y2[n];
      dot_d += Cn * y2[n_dof + n];
      dot_e += Cn * (*E_pt)[n];
      dot_f += Cn * (*E_pt)[n_dof + n];
      dot_g += Cn * (*G_pt)[n];
      dot_h += Cn * (*G_pt)[n_dof + n];
    }
 
    // Now we should be able to work out the corrections
    double denom = dot_e * dot_h - dot_g * dot_f;
 
    // Copy the previous residuals
    double R31 = result[3 * n_dof], R32 = result[3 * n_dof + 1];
    // Delta parameter
    const double delta_param =
      ((R32 - dot_d) * dot_g - (R31 - dot_c) * dot_h) / denom;
    // Delta frequency
    const double delta_w = -((R32 - dot_d) + dot_f * delta_param) / (dot_h);
 
    // Load into the result vector
    result[3 * n_dof] = delta_param;
    result[3 * n_dof + 1] = delta_w;
 
    // The corrections to the null vector
    for (unsigned n = 0; n < 2 * n_dof; n++)
    {
      result[n_dof + n] =
        y2[n] - (*E_pt)[n] * delta_param - (*G_pt)[n] * delta_w;
    }
 
    // Finally add the corrections to the unknowns
    for (unsigned n = 0; n < n_dof; n++)
    {
      result[n] = y1[n] - (*A_pt)[n] * delta_param;
    }
 
    // The sign of the jacobian is the previous signs multiplied by the
    // sign of the denominator
    problem_pt->sign_of_jacobian() =
      sign_of_jacobian * static_cast<int>(std::fabs(denom) / denom);
 
    // Switch things to our full solver
    handler_pt->solve_full_system();
 
    // If we are not storing things, clear the results
    if (!Enable_resolve)
    {
      // We no longer need to store the matrix
      Linear_solver_pt->disable_resolve();
      delete A_pt;
      A_pt = 0;
      delete E_pt;
      E_pt = 0;
      delete G_pt;
      G_pt = 0;
    }
    // Otherwise, also store the problem pointer
    else
    {
      Problem_pt = problem_pt;
    }
  }

References A_pt, oomph::Problem::actions_after_change_in_bifurcation_parameter(), alpha, oomph::Problem::assembly_handler_pt(), oomph::LinearAlgebraDistribution::build(), oomph::DistributableLinearAlgebraObject::build_distribution(), oomph::DoubleVector::built(), oomph::HopfHandler::C, oomph::Problem::communicator_pt(), oomph::LinearSolver::disable_resolve(), oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::Problem::dof(), e(), E_pt, oomph::Mesh::element_pt(), oomph::LinearSolver::Enable_resolve, oomph::LinearSolver::enable_resolve(), oomph::GeneralisedElement::eqn_number(), boost::multiprecision::fabs(), oomph::BlackBoxFDNewtonSolver::FD_step, G_pt, oomph::GeneralisedElement::get_jacobian_and_mass_matrix(), oomph::Problem::get_residuals(), i, Linear_solver_pt, m, oomph::Problem::mesh_pt(), n, oomph::GeneralisedElement::ndof(), oomph::Problem::ndof(), oomph::Mesh::nelement(), Eigen::internal::omega(), oomph::HopfHandler::Omega, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::HopfHandler::Phi, Problem_pt, oomph::HopfHandler::Psi, oomph::LinearSolver::resolve(), oomph::Problem::sign_of_jacobian(), oomph::LinearSolver::solve(), oomph::HopfHandler::solve_complex_system(), oomph::HopfHandler::solve_full_system(), oomph::HopfHandler::solve_standard_system(), and y.

solve_for_two_rhs()

void oomph::BlockHopfLinearSolver::solve_for_two_rhs	(	Problem *const &	problem_pt,
		DoubleVector &	result,
		const DoubleVector &	rhs2,
		DoubleVector &	result2
	)

Solve for two right hand sides.

Solve for two right hand sides, required for (efficient) continuation because otherwise we have to store the inverse of the jacobian and the complex equivalent ... nasty.

  {
    // if the result is setup then it should not be distributed
#ifdef PARANOID
    if (result.built())
    {
      if (result.distributed())
      {
        throw OomphLibError("The result vector must not be distributed",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
    if (result2.built())
    {
      if (result2.distributed())
      {
        throw OomphLibError("The result2 vector must not be distributed",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // Locally cache the pointer to the handler.
    HopfHandler* handler_pt =
      static_cast<HopfHandler*>(problem_pt->assembly_handler_pt());
 
    // Find the number of dofs in the augmented problem
    unsigned n_dof = problem_pt->ndof();
 
    // create the linear algebra distribution for this solver
    // currently only global (non-distributed) distributions are allowed
    LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n_dof, false);
    this->build_distribution(dist);
 
    // if the result vector is not setup then rebuild with distribution = global
    if (!result.built())
    {
      result.build(this->distribution_pt(), 0.0);
    }
    if (!result2.built())
    {
      result2.build(this->distribution_pt(), 0.0);
    }
 
    // Firstly, let's calculate the derivative of the residuals wrt
    // the parameter
    DoubleVector dRdparam(this->distribution_pt(), 0.0);
 
    const double FD_step = 1.0e-8;
    {
      // Cache pointer to the parameter (second last entry in the vector
      double* param_pt = &problem_pt->dof(n_dof - 2);
      // Backup the parameter
      double old_var = *param_pt;
      ;
      // Increment the parameter
      *param_pt += FD_step;
      problem_pt->actions_after_change_in_bifurcation_parameter();
 
      // Now calculate the new residuals
      problem_pt->get_residuals(dRdparam);
 
      // Now calculate the difference assume original residuals in resul
      for (unsigned n = 0; n < n_dof; n++)
      {
        dRdparam[n] = (dRdparam[n] - result[n]) / FD_step;
      }
 
      // Reset the parameter
      *param_pt = old_var;
      problem_pt->actions_after_change_in_bifurcation_parameter();
    }
 
    // Switch the handler to "standard" mode
    handler_pt->solve_standard_system();
 
    // Find out the number of dofs in the non-standard problem
    n_dof = problem_pt->ndof();
 
    //
    dist.build(problem_pt->communicator_pt(), n_dof, false);
    this->build_distribution(dist);
 
    // Setup storage for temporary vector
    DoubleVector y1(this->distribution_pt(), 0.0),
      alpha(this->distribution_pt(), 0.0),
      y1_resolve(this->distribution_pt(), 0.0);
 
    // Allocate storage for A which can be used in the resolve
    if (A_pt != 0)
    {
      delete A_pt;
    }
    A_pt = new DoubleVector(this->distribution_pt(), 0.0);
 
    // We are going to do resolves using the underlying linear solver
    Linear_solver_pt->enable_resolve();
 
    // Solve the first system Jy1 =
    Linear_solver_pt->solve(problem_pt, y1);
 
    // This should have set the sign of the jacobian, store it
    int sign_of_jacobian = problem_pt->sign_of_jacobian();
 
    // Now let's get the appropriate bit of alpha
    for (unsigned n = 0; n < n_dof; n++)
    {
      alpha[n] = dRdparam[n];
    }
 
    // Resolve to find A
    Linear_solver_pt->resolve(alpha, *A_pt);
 
    // Get the solution for the second rhs
    for (unsigned n = 0; n < n_dof; n++)
    {
      alpha[n] = rhs2[n];
    }
 
    // Resolve to find y1_resolve
    Linear_solver_pt->resolve(alpha, y1_resolve);
 
    // Now set to the complex system
    handler_pt->solve_complex_system();
 
    // rebuild the Distribution
    dist.build(problem_pt->communicator_pt(), n_dof * 2, false);
    this->build_distribution(dist);
 
    // Resize the stored vector G
    if (G_pt != 0)
    {
      delete G_pt;
    }
    G_pt = new DoubleVector(this->distribution_pt(), 0.0);
 
    // Solve the first Zg = (Mz -My)
    Linear_solver_pt->solve(problem_pt, *G_pt);
 
    // This should have set the sign of the complex matrix's determinant,
    // multiply
    sign_of_jacobian *= problem_pt->sign_of_jacobian();
 
    // We can now construct our multipliers
    // Prepare to scale
    double dof_length = 0.0, a_length = 0.0, y1_length = 0.0,
           y1_resolve_length = 0.0;
    // Loop over the standard number of dofs
    for (unsigned n = 0; n < n_dof; n++)
    {
      if (std::fabs(problem_pt->dof(n)) > dof_length)
      {
        dof_length = std::fabs(problem_pt->dof(n));
      }
      if (std::fabs((*A_pt)[n]) > a_length)
      {
        a_length = std::fabs((*A_pt)[n]);
      }
      if (std::fabs(y1[n]) > y1_length)
      {
        y1_length = std::fabs(y1[n]);
      }
      if (std::fabs(y1_resolve[n]) > y1_resolve_length)
      {
        y1_resolve_length = std::fabs(y1[n]);
      }
    }
 
 
    double a_mult = dof_length / a_length;
    double y1_mult = dof_length / y1_length;
    double y1_resolve_mult = dof_length / y1_resolve_length;
    a_mult += FD_step;
    y1_mult += FD_step;
    y1_resolve_mult += FD_step;
    a_mult *= FD_step;
    y1_mult *= FD_step;
    y1_resolve_mult *= FD_step;
 
    // Local storage for the product terms
    Vector<double> Jprod_a(2 * n_dof, 0.0), Jprod_y1(2 * n_dof, 0.0);
    Vector<double> Jprod_y1_resolve(2 * n_dof, 0.0);
 
    // Temporary storage
    Vector<double> rhs(2 * n_dof);
 
    // find the number of elements
    unsigned long n_element = problem_pt->mesh_pt()->nelement();
 
    // Calculate the product of the jacobian matrices, etc
    for (unsigned long e = 0; e < n_element; e++)
    {
      GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);
      // Loop over the ndofs in each element
      unsigned n_var = elem_pt->ndof();
      // Get the jacobian matrices
      DenseMatrix<double> jac(n_var), jac_a(n_var), jac_y1(n_var),
        jac_y1_resolve(n_var);
      DenseMatrix<double> M(n_var), M_a(n_var), M_y1(n_var),
        M_y1_resolve(n_var);
      // Get unperturbed jacobian
      elem_pt->get_jacobian_and_mass_matrix(rhs, jac, M);
 
      // Backup the dofs
      Vector<double> dof_bac(n_var);
      // Perturb the original dofs
      for (unsigned n = 0; n < n_var; n++)
      {
        unsigned eqn_number = elem_pt->eqn_number(n);
        dof_bac[n] = problem_pt->dof(eqn_number);
        // Pertub by vector A
        problem_pt->dof(eqn_number) += a_mult * (*A_pt)[eqn_number];
      }
 
      // Now get the new jacobian and mass matrix
      elem_pt->get_jacobian_and_mass_matrix(rhs, jac_a, M_a);
 
      // Perturb the dofs
      for (unsigned n = 0; n < n_var; n++)
      {
        unsigned eqn_number = elem_pt->eqn_number(n);
        problem_pt->dof(eqn_number) = dof_bac[n];
        // Pertub by vector y1
        problem_pt->dof(eqn_number) += y1_mult * y1[eqn_number];
      }
 
      // Now get the new jacobian and mass matrix
      elem_pt->get_jacobian_and_mass_matrix(rhs, jac_y1, M_y1);
 
      // Perturb the dofs
      for (unsigned n = 0; n < n_var; n++)
      {
        unsigned eqn_number = elem_pt->eqn_number(n);
        problem_pt->dof(eqn_number) = dof_bac[n];
        // Pertub by vector y1
        problem_pt->dof(eqn_number) += y1_resolve_mult * y1_resolve[eqn_number];
      }
 
      // Now get the new jacobian and mass matrix
      elem_pt->get_jacobian_and_mass_matrix(rhs, jac_y1_resolve, M_y1_resolve);
 
      // Reset the dofs
      for (unsigned n = 0; n < n_var; n++)
      {
        unsigned eqn_number = elem_pt->eqn_number(n);
        problem_pt->dof(eqn_number) = dof_bac[n];
      }
 
      // OK, now work out the products
      for (unsigned n = 0; n < n_var; n++)
      {
        unsigned eqn_number = elem_pt->eqn_number(n);
        double prod_a1 = 0.0, prod_y11 = 0.0, prod_y1_resolve1 = 0.0;
        double prod_a2 = 0.0, prod_y12 = 0.0, prod_y1_resolve2 = 0.0;
        for (unsigned m = 0; m < n_var; m++)
        {
          const unsigned unknown = elem_pt->eqn_number(m);
          const double y = handler_pt->Phi[unknown];
          const double z = handler_pt->Psi[unknown];
          const double omega = handler_pt->Omega;
          // Real part (first line)
          prod_a1 +=
            (jac_a(n, m) - jac(n, m)) * y + omega * (M_a(n, m) - M(n, m)) * z;
          prod_y11 +=
            (jac_y1(n, m) - jac(n, m)) * y + omega * (M_y1(n, m) - M(n, m)) * z;
          prod_y1_resolve1 += (jac_y1_resolve(n, m) - jac(n, m)) * y +
                              omega * (M_y1_resolve(n, m) - M(n, m)) * z;
          // Imag par (second line)
          prod_a2 +=
            (jac_a(n, m) - jac(n, m)) * z - omega * (M_a(n, m) - M(n, m)) * y;
          prod_y12 +=
            (jac_y1(n, m) - jac(n, m)) * z - omega * (M_y1(n, m) - M(n, m)) * y;
          prod_y1_resolve2 += (jac_y1_resolve(n, m) - jac(n, m)) * z -
                              omega * (M_y1_resolve(n, m) - M(n, m)) * y;
        }
        Jprod_a[eqn_number] += prod_a1 / a_mult;
        Jprod_y1[eqn_number] += prod_y11 / y1_mult;
        Jprod_y1_resolve[eqn_number] += prod_y1_resolve1 / y1_resolve_mult;
        Jprod_a[n_dof + eqn_number] += prod_a2 / a_mult;
        Jprod_y1[n_dof + eqn_number] += prod_y12 / y1_mult;
        Jprod_y1_resolve[n_dof + eqn_number] +=
          prod_y1_resolve2 / y1_resolve_mult;
      }
    }
 
    // The assumption here is still that the result has been set to the
    // residuals.
    for (unsigned n = 0; n < 2 * n_dof; n++)
    {
      rhs[n] = result[n_dof + n] - Jprod_y1[n];
    }
 
    // Temporary storage
    DoubleVector y2(this->distribution_pt(), 0.0);
 
    // Solve it
    DoubleVector temp_rhs(this->distribution_pt(), 0.0);
    for (unsigned i = 0; i < 2 * n_dof; i++)
    {
      temp_rhs[i] = rhs[i];
    }
    Linear_solver_pt->resolve(temp_rhs, y2);
 
    // Assemble the next RHS
    for (unsigned n = 0; n < 2 * n_dof; n++)
    {
      rhs[n] = dRdparam[n_dof + n] - Jprod_a[n];
    }
 
    // Resive the storage
    if (E_pt != 0)
    {
      delete E_pt;
    }
    E_pt = new DoubleVector(this->distribution_pt(), 0.0);
 
    // Solve for the next RHS
    for (unsigned i = 0; i < 2 * n_dof; i++)
    {
      temp_rhs[i] = rhs[i];
    }
    Linear_solver_pt->resolve(temp_rhs, *E_pt);
 
    // Assemble the next RHS
    for (unsigned n = 0; n < 2 * n_dof; n++)
    {
      rhs[n] = rhs2[n_dof + n] - Jprod_y1_resolve[n];
    }
 
    DoubleVector y2_resolve(this->distribution_pt(), 0.0);
    for (unsigned i = 0; i < 2 * n_dof; i++)
    {
      temp_rhs[i] = rhs[i];
    }
    Linear_solver_pt->resolve(temp_rhs, y2_resolve);
 
 
    // We can now calculate the final corrections
    // We need to work out a large number of dot products
    double dot_c = 0.0, dot_d = 0.0, dot_e = 0.0, dot_f = 0.0, dot_g = 0.0;
    double dot_h = 0.0;
 
    double dot_c_resolve = 0.0, dot_d_resolve = 0.0;
 
    for (unsigned n = 0; n < n_dof; n++)
    {
      // Get the appopriate entry
      const double Cn = handler_pt->C[n];
      dot_c += Cn * y2[n];
      dot_d += Cn * y2[n_dof + n];
      dot_c_resolve += Cn * y2_resolve[n];
      dot_d_resolve += Cn * y2_resolve[n_dof + n];
      dot_e += Cn * (*E_pt)[n];
      dot_f += Cn * (*E_pt)[n_dof + n];
      dot_g += Cn * (*G_pt)[n];
      dot_h += Cn * (*G_pt)[n_dof + n];
    }
 
    // Now we should be able to work out the corrections
    double denom = dot_e * dot_h - dot_g * dot_f;
 
    // Copy the previous residuals
    double R31 = result[3 * n_dof], R32 = result[3 * n_dof + 1];
    // Delta parameter
    const double delta_param =
      ((R32 - dot_d) * dot_g - (R31 - dot_c) * dot_h) / denom;
    // Delta frequency
    const double delta_w = -((R32 - dot_d) + dot_f * delta_param) / (dot_h);
 
    // Corrections
    double R31_resolve = rhs2[3 * n_dof], R32_resolve = rhs2[3 * n_dof + 1];
    // Delta parameter
    const double delta_param_resolve = ((R32_resolve - dot_d_resolve) * dot_g -
                                        (R31_resolve - dot_c_resolve) * dot_h) /
                                       denom;
    // Delta frequency
    const double delta_w_resolve =
      -((R32_resolve - dot_d_resolve) + dot_f * delta_param_resolve) / (dot_h);
 
 
    // Load into the result vector
    result[3 * n_dof] = delta_param;
    result[3 * n_dof + 1] = delta_w;
 
    // The corrections to the null vector
    for (unsigned n = 0; n < 2 * n_dof; n++)
    {
      result[n_dof + n] =
        y2[n] - (*E_pt)[n] * delta_param - (*G_pt)[n] * delta_w;
    }
 
    // Finally add the corrections to the unknowns
    for (unsigned n = 0; n < n_dof; n++)
    {
      result[n] = y1[n] - (*A_pt)[n] * delta_param;
    }
 
    // Load into the result vector
    result2[3 * n_dof] = delta_param_resolve;
    result2[3 * n_dof + 1] = delta_w_resolve;
 
    // The corrections to the null vector
    for (unsigned n = 0; n < 2 * n_dof; n++)
    {
      result2[n_dof + n] = y2_resolve[n] - (*E_pt)[n] * delta_param_resolve -
                           (*G_pt)[n] * delta_w_resolve;
    }
 
    // Finally add the corrections to the unknowns
    for (unsigned n = 0; n < n_dof; n++)
    {
      result2[n] = y1_resolve[n] - (*A_pt)[n] * delta_param_resolve;
    }
 
 
    // The sign of the jacobian is the previous signs multiplied by the
    // sign of the denominator
    problem_pt->sign_of_jacobian() =
      sign_of_jacobian * static_cast<int>(std::fabs(denom) / denom);
 
    // Switch things to our full solver
    handler_pt->solve_full_system();
 
    // If we are not storing things, clear the results
    if (!Enable_resolve)
    {
      // We no longer need to store the matrix
      Linear_solver_pt->disable_resolve();
      delete A_pt;
      A_pt = 0;
      delete E_pt;
      E_pt = 0;
      delete G_pt;
      G_pt = 0;
    }
    // Otherwise, also store the problem pointer
    else
    {
      Problem_pt = problem_pt;
    }
  }

References A_pt, oomph::Problem::actions_after_change_in_bifurcation_parameter(), alpha, oomph::Problem::assembly_handler_pt(), oomph::DoubleVector::build(), oomph::LinearAlgebraDistribution::build(), oomph::DistributableLinearAlgebraObject::build_distribution(), oomph::DoubleVector::built(), oomph::HopfHandler::C, oomph::Problem::communicator_pt(), oomph::LinearSolver::disable_resolve(), oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::Problem::dof(), e(), E_pt, oomph::Mesh::element_pt(), oomph::LinearSolver::Enable_resolve, oomph::LinearSolver::enable_resolve(), oomph::GeneralisedElement::eqn_number(), boost::multiprecision::fabs(), oomph::BlackBoxFDNewtonSolver::FD_step, G_pt, oomph::GeneralisedElement::get_jacobian_and_mass_matrix(), oomph::Problem::get_residuals(), i, Linear_solver_pt, m, oomph::Problem::mesh_pt(), n, oomph::GeneralisedElement::ndof(), oomph::Problem::ndof(), oomph::Mesh::nelement(), Eigen::internal::omega(), oomph::HopfHandler::Omega, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::HopfHandler::Phi, Problem_pt, oomph::HopfHandler::Psi, oomph::LinearSolver::resolve(), oomph::Problem::sign_of_jacobian(), oomph::LinearSolver::solve(), oomph::HopfHandler::solve_complex_system(), oomph::HopfHandler::solve_full_system(), oomph::HopfHandler::solve_standard_system(), and y.

Member Data Documentation

A_pt

DoubleVector* oomph::BlockHopfLinearSolver::A_pt

private

Pointer to the storage for the vector a.

Referenced by solve(), solve_for_two_rhs(), and ~BlockHopfLinearSolver().

E_pt

DoubleVector* oomph::BlockHopfLinearSolver::E_pt

private

Pointer to the storage for the vector e (0 to n-1)

Referenced by solve(), solve_for_two_rhs(), and ~BlockHopfLinearSolver().

G_pt

DoubleVector* oomph::BlockHopfLinearSolver::G_pt

private

Pointer to the storage for the vector g (0 to n-1)

Referenced by solve(), solve_for_two_rhs(), and ~BlockHopfLinearSolver().

Linear_solver_pt

LinearSolver* oomph::BlockHopfLinearSolver::Linear_solver_pt

private

Pointer to the original linear solver.

Referenced by linear_solver_pt(), solve(), and solve_for_two_rhs().

Problem_pt

Problem* oomph::BlockHopfLinearSolver::Problem_pt

private

Pointer to the problem, used in the resolve.

Referenced by solve(), and solve_for_two_rhs().

---

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

• assembly_handler.h
• assembly_handler.cc