oomph::PitchForkHandler Class Reference

#include <assembly_handler.h>

Inheritance diagram for oomph::PitchForkHandler:

Public Member Functions
	PitchForkHandler (Problem *const &problem_pt, AssemblyHandler *const &assembly_handler_pt, double *const &parameter_pt, const DoubleVector &symmetry_vector)
	Constructor, initialise the systems. More...
	~PitchForkHandler ()
	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)
	Get the derivatives of the residuals with respect to a parameter. More...
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)
void	solve_augmented_block_system ()
	Set to solve the augmented block system. More...
void	solve_block_system ()
	Set to solve the block 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 Types
enum	{ Full_augmented , Block_J , Block_augmented_J }

Private Member Functions
unsigned	global_eqn_number (const unsigned &i)

Private Attributes
unsigned	Solve_which_system
Problem *	Problem_pt
	Pointer to the problem. More...
AssemblyHandler *	Assembly_handler_pt
	Pointer to the underlying (original) assembly handler. More...
unsigned	Ndof
LinearAlgebraDistribution *	Dof_distribution_pt
	Store the original dof distribution. More...
LinearAlgebraDistribution *	Augmented_dof_distribution_pt
	The augmented distribution. More...
double	Sigma
DoubleVectorWithHaloEntries	Y
	Storage for the null vector. More...
DoubleVectorWithHaloEntries	Psi
	A constant vector that is specifies the symmetry being broken. More...
DoubleVectorWithHaloEntries	C
DoubleVectorWithHaloEntries	Count
Vector< unsigned >	Global_eqn_number
double *	Parameter_pt
	Storage for the pointer to the parameter. More...
unsigned	Nelement

Friends
class	BlockPitchForkLinearSolver
class	AugmentedBlockPitchForkLinearSolver

Detailed Description

A class that is used to assemble the augmented system that defines a pitchfork (symmetry-breaking) 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 pitchfork bifurcation may be specified by the augmented system of size \(2N+2\).

\[ R(u,\lambda) + \sigma \psi = 0,\]

\[ Jy = 0,\]

\[ \langle u, \phi \rangle = 0, \]

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

In the abovem \(J\) is the usual Jacobian matrix, \(dR/du\) and \(\phi\) is a constant vector that is chosen to ensure that the null vector, \(y\), is not trivial. Here \(\sigma \) is a slack variable that is used to enforce the constraint \( \langle u, \phi \rangle = 0 \) — the inner product of the solution with the chosen symmetry vector is zero. At the bifurcation \( \sigma \) should be very close to zero. Note that this formulation means that any odd symmetry in the problem must be about zero. Moreover, we use the dot product of two vectors to calculate the inner product, rather than an integral over the domain and so the mesh must be symmetric in the appropriate directions.

Member Enumeration Documentation

anonymous enum

anonymous enum

private

A little private enum to determine whether we are solving the block system or not

Enumerator
Full_augmented
Block_J
Block_augmented_J

    {
      Full_augmented,
      Block_J,
      Block_augmented_J
    };

Constructor & Destructor Documentation

PitchForkHandler()

oomph::PitchForkHandler::PitchForkHandler	(	Problem *const &	problem_pt,
		AssemblyHandler *const &	assembly_handler_pt,
		double *const &	parameter_pt,
		const DoubleVector &	symmetry_vector
	)

Constructor, initialise the systems.

Constructor: Initialise the PitchForkHandler by setting intial guesses for Sigma, Y, specifying C and Psi and calculating count.

    : Solve_which_system(Full_augmented), Sigma(0.0), Parameter_pt(parameter_pt)
  {
    // Set the problem pointer
    Problem_pt = problem_pt;
    // Set the assembly handler
    Assembly_handler_pt = assembly_handler_pt;
    // Set the number of degrees of freedom
    Ndof = Problem_pt->ndof();
    // Backup the distribution
    Dof_distribution_pt = Problem_pt->Dof_distribution_pt;
#ifdef OOMPH_HAS_MPI
    // Set the distribution flag
    Distributed = Problem_pt->distributed();
#endif
 
    // Find the elements in the problem
    unsigned n_element = Problem_pt->mesh_pt()->nelement();
 
#ifdef OOMPH_HAS_MPI
 
    // Work out all the global equations to which this processor
    // contributes and store the halo data in the problem
    Problem_pt->setup_dof_halo_scheme();
 
#endif
 
 
    // Now use the dof distribution for all double vectors
    Psi.build(Dof_distribution_pt);
    Y.build(Dof_distribution_pt);
    C.build(Dof_distribution_pt);
    Count.build(Dof_distribution_pt);
#ifdef OOMPH_HAS_MPI
    Count.build_halo_scheme(Problem_pt->Halo_scheme_pt);
#endif
 
    // Loop over the elements and count the entries
    // and number of (non-halo) elements
    unsigned n_non_halo_element_local = 0;
    for (unsigned e = 0; e < n_element; e++)
    {
      GeneralisedElement* elem_pt = Problem_pt->mesh_pt()->element_pt(e);
#ifdef OOMPH_HAS_MPI
      // Ignore halo elements
      if (!elem_pt->is_halo())
      {
#endif
        // Increment the number of non halo elements
        ++n_non_halo_element_local;
        // Now count the number of times the element contributes to a value
        unsigned n_var = assembly_handler_pt->ndof(elem_pt);
        for (unsigned n = 0; n < n_var; n++)
        {
          ++Count.global_value(assembly_handler_pt->eqn_number(elem_pt, n));
        }
#ifdef OOMPH_HAS_MPI
      }
#endif
    }
 
    // Add together all the counts
#ifdef OOMPH_HAS_MPI
    Count.sum_all_halo_and_haloed_values();
 
    // If distributed, find the total number of elements in the problem
    if (Distributed)
    {
      // Need to gather the total number of non halo elements
      MPI_Allreduce(&n_non_halo_element_local,
                    &Nelement,
                    1,
                    MPI_UNSIGNED,
                    MPI_SUM,
                    Problem_pt->communicator_pt()->mpi_comm());
    }
    // Otherwise the total number is the same on each processor
    else
#endif
    {
      Nelement = n_non_halo_element_local;
    }
 
#ifdef OOMPH_HAS_MPI
    // Only add the parameter to the first processor, if distributed
    int my_rank = Problem_pt->communicator_pt()->my_rank();
    if ((!Distributed) || (my_rank == 0))
#endif
    {
      // Add the parameter to the problem
      Problem_pt->Dof_pt.push_back(parameter_pt);
    }
 
    // Find length of the symmetry vector
    double length = symmetry_vector.norm();
 
    // Add the unknowns for the null vector to the problem
    // Normalise the symmetry vector and initialise the null vector to the
    // symmetry vector and set the constant vector, c, to the
    // normalised symmetry vector.
 
    // Need to redistribute the symmetry vector into the (natural)
    // distribution of the Dofs
#ifdef OOMPH_HAS_MPI
    // Check that the symmetry vector has the correct distribution
    if (!(*symmetry_vector.distribution_pt() == *Dof_distribution_pt))
    {
      throw OomphLibError(
        "The symmetry vector must have the same distribution as the dofs\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Only loop over the local unknowns
    const unsigned n_dof_local = Dof_distribution_pt->nrow_local();
    for (unsigned n = 0; n < n_dof_local; n++)
    {
      Problem_pt->Dof_pt.push_back(&Y[n]);
      // Psi[n] = symmetry_vector[n];
      Psi[n] = Y[n] = C[n] = symmetry_vector[n] / length;
    }
 
 
#ifdef OOMPH_HAS_MPI
    // Set up the required halo schemes (which also synchronises)
    Psi.build_halo_scheme(Problem_pt->Halo_scheme_pt);
    Y.build_halo_scheme(Problem_pt->Halo_scheme_pt);
    C.build_halo_scheme(Problem_pt->Halo_scheme_pt);
 
 
    if ((!Distributed) || (my_rank == 0))
#endif
    // Add the slack parameter to the problem on the first processor
    // if distributed
    {
      Problem_pt->Dof_pt.push_back(&Sigma);
    }
 
#ifdef OOMPH_HAS_MPI
    if (Distributed)
    {
      unsigned augmented_first_row = 0;
      unsigned augmented_n_row_local = 0;
 
      // Set up the translation scheme on every processor
      Global_eqn_number.resize(2 * Ndof + 2);
      int n_proc = Problem_pt->communicator_pt()->nproc();
      unsigned global_eqn_count = 0;
      for (int d = 0; d < n_proc; d++)
      {
        // Find out the first row of the current processor
        if (my_rank == d)
        {
          augmented_first_row = global_eqn_count;
        }
 
        const unsigned n_row_local = Dof_distribution_pt->nrow_local(d);
        const unsigned first_row = Dof_distribution_pt->first_row(d);
        // Add the basic equations
        for (unsigned n = 0; n < n_row_local; n++)
        {
          Global_eqn_number[first_row + n] = global_eqn_count;
          ++global_eqn_count;
        }
        // If on the first processor add the pointer to the parameter
        if (d == 0)
        {
          Global_eqn_number[Ndof] = global_eqn_count;
          ++global_eqn_count;
        }
        // Add the eigenfunction
        for (unsigned n = 0; n < n_row_local; n++)
        {
          Global_eqn_number[Ndof + 1 + first_row + n] = global_eqn_count;
          ++global_eqn_count;
        }
        // Finally add the slack parameter
        if (d == 0)
        {
          Global_eqn_number[2 * Ndof + 1] = global_eqn_count;
          ++global_eqn_count;
        }
        // Find out the number of rows of the current processor
        if (my_rank == d)
        {
          augmented_n_row_local = global_eqn_count - augmented_first_row;
        }
      }
 
      // Make a new linear algebra distribution
      Augmented_dof_distribution_pt =
        new LinearAlgebraDistribution(Problem_pt->communicator_pt(),
                                      augmented_first_row,
                                      augmented_n_row_local);
    }
    else
#endif
    {
      Augmented_dof_distribution_pt = new LinearAlgebraDistribution(
        Problem_pt->communicator_pt(), 2 * Ndof + 2, false);
    }
 
    // resize
    // Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
    //                                       2*Ndof+2,false);
 
    Problem_pt->Dof_distribution_pt = Augmented_dof_distribution_pt;
 
    // Remove all previous sparse storage used during Jacobian assembly
    Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
  }

References Assembly_handler_pt, Augmented_dof_distribution_pt, oomph::DoubleVector::build(), oomph::DoubleVectorWithHaloEntries::build_halo_scheme(), oomph::Problem::communicator_pt(), Count, oomph::Problem::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), Dof_distribution_pt, oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, e(), oomph::Mesh::element_pt(), oomph::AssemblyHandler::eqn_number(), oomph::LinearAlgebraDistribution::first_row(), Global_eqn_number, oomph::DoubleVectorWithHaloEntries::global_value(), oomph::Problem::mesh_pt(), oomph::OomphCommunicator::my_rank(), n, Ndof, oomph::Problem::ndof(), oomph::AssemblyHandler::ndof(), Nelement, oomph::Mesh::nelement(), oomph::DoubleVector::norm(), oomph::OomphCommunicator::nproc(), oomph::LinearAlgebraDistribution::nrow_local(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Problem_pt, Psi, Sigma, oomph::Problem::Sparse_assemble_with_arrays_previous_allocation, oomph::DoubleVectorWithHaloEntries::sum_all_halo_and_haloed_values(), and Y.

~PitchForkHandler()

oomph::PitchForkHandler::~PitchForkHandler

(

)

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

References Dof_distribution_pt, oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, oomph::Problem::linear_solver_pt(), oomph::BlockPitchForkLinearSolver::linear_solver_pt(), oomph::AugmentedBlockPitchForkLinearSolver::linear_solver_pt(), oomph::LinearAlgebraDistribution::nrow_local(), Problem_pt, and oomph::Problem::Sparse_assemble_with_arrays_previous_allocation.

Member Function Documentation

bifurcation_parameter_pt()

double* oomph::PitchForkHandler::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::PitchForkHandler::bifurcation_type ( ) const

inlinevirtual

Indicate that we are tracking a pitchfork bifurcation by returning 2

Reimplemented from oomph::AssemblyHandler.

    {
      return 2;
    }

eqn_number()

unsigned long oomph::PitchForkHandler::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 = 0;
 
    // Which system are we solving
    switch (Solve_which_system)
    {
        // In the block case, it's just the standard numbering
      case Block_J:
        global_eqn = elem_pt->eqn_number(ieqn_local);
        break;
 
        // Block augmented not done properly yet
      case Block_augmented_J:
        // Not done the distributed case
#ifdef OOMPH_HAS_MPI
        if (Distributed)
        {
          throw OomphLibError(
            "Block Augmented solver not implemented for distributed case\n",
            OOMPH_CURRENT_FUNCTION,
            OOMPH_EXCEPTION_LOCATION);
        }
#endif
 
        // Full case
      case Full_augmented:
        // The usual equations
        if (ieqn_local < raw_ndof)
        {
          global_eqn = this->global_eqn_number(elem_pt->eqn_number(ieqn_local));
        }
        // The bifurcation parameter equation
        else if (ieqn_local == raw_ndof)
        {
          global_eqn = this->global_eqn_number(Ndof);
        }
        // If we are assembling the full system we also have
        // The components of the null vector
        else if (ieqn_local < (2 * raw_ndof + 1))
        {
          global_eqn = this->global_eqn_number(
            Ndof + 1 + elem_pt->eqn_number(ieqn_local - 1 - raw_ndof));
        }
        // The slack parameter
        else
        {
          global_eqn = this->global_eqn_number(2 * Ndof + 1);
        }
        break;
    } // End of switch
    return global_eqn;
  }

References Block_augmented_J, Block_J, oomph::GeneralisedElement::eqn_number(), Full_augmented, global_eqn_number(), Ndof, oomph::GeneralisedElement::ndof(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and Solve_which_system.

Referenced by oomph::AugmentedBlockPitchForkLinearSolver::resolve(), and oomph::AugmentedBlockPitchForkLinearSolver::solve().

get_djacobian_dparameter()

void oomph::PitchForkHandler::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::PitchForkHandler::get_dresiduals_dparameter ( GeneralisedElement *const &  elem_pt, double *const &  parameter_pt, Vector< double > &  dres_dparam  )

virtual

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

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

Reimplemented from oomph::AssemblyHandler.

  {
    // Need to get raw residuals and jacobian
    unsigned raw_ndof = elem_pt->ndof();
    Problem_pt->actions_before_newton_convergence_check();
    // Find out which system we are solving
    switch (Solve_which_system)
    {
        // If we are solving the original system
      case Block_J:
      {
        // get the basic residual derivatives
        elem_pt->get_dresiduals_dparameter(parameter_pt, dres_dparam);
        // Slack parameter term does not depened explicitly on the parameter
        // so is not added
      }
      break;
 
      // If we are solving the augmented-by-one system
      case Block_augmented_J:
      {
        // Get the basic residual derivatives
        elem_pt->get_dresiduals_dparameter(parameter_pt, dres_dparam);
 
        // Zero the final residuals derivative
        dres_dparam[raw_ndof] = 0.0;
 
        // Other terms must not depend on the parameter
      }
      break;
 
      // Otherwise we are solving the fully augemented system
      case Full_augmented:
      {
        DenseMatrix<double> djac_dparam(raw_ndof);
 
        // Get the basic residuals and jacobian derivatives
        elem_pt->get_djacobian_dparameter(
          parameter_pt, dres_dparam, djac_dparam);
 
        // The "final" residual derivatives do not depend on the parameter
        dres_dparam[raw_ndof] = 0.0;
        dres_dparam[2 * raw_ndof + 1] = 0.0;
 
        // Now multiply to fill in the residuals associated
        // with the null vector condition
        for (unsigned i = 0; i < raw_ndof; i++)
        {
          dres_dparam[raw_ndof + 1 + i] = 0.0;
          for (unsigned j = 0; j < raw_ndof; j++)
          {
            unsigned local_unknown = elem_pt->eqn_number(j);
            dres_dparam[raw_ndof + 1 + i] +=
              djac_dparam(i, j) * Y.global_value(local_unknown);
          }
        }
      }
      break;
 
      default:
        std::ostringstream error_stream;
        error_stream
          << "The Solve_which_system flag can only take values 0, 1, 2"
          << " not " << Solve_which_system << "\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  }

References oomph::Problem::actions_before_newton_convergence_check(), Block_augmented_J, Block_J, oomph::GeneralisedElement::eqn_number(), Full_augmented, oomph::GeneralisedElement::get_djacobian_dparameter(), oomph::GeneralisedElement::get_dresiduals_dparameter(), oomph::DoubleVectorWithHaloEntries::global_value(), i, j, oomph::GeneralisedElement::ndof(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Problem_pt, Solve_which_system, and Y.

get_eigenfunction()

void oomph::PitchForkHandler::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 only one (real) null vector
    eigenfunction.resize(1);
    // Rebuild the vector
    eigenfunction[0].build(this->Dof_distribution_pt, 0.0);
    // Set the value to be the null vector
    const unsigned n_row_local = eigenfunction[0].nrow_local();
    for (unsigned n = 0; n < n_row_local; n++)
    {
      eigenfunction[0][n] = Y[n];
    }
  }

References Dof_distribution_pt, n, and Y.

get_hessian_vector_products()

void oomph::PitchForkHandler::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 calls the underlying element's hessian function.

Reimplemented from oomph::AssemblyHandler.

  {
    elem_pt->get_hessian_vector_products(Y, C, product);
  }

References oomph::GeneralisedElement::get_hessian_vector_products(), product(), and Y.

get_jacobian()

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

virtual

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

ALICE:

ALICE

Reimplemented from oomph::AssemblyHandler.

  {
    unsigned augmented_ndof = ndof(elem_pt);
    unsigned raw_ndof = elem_pt->ndof();
 
    // Which system are we solving
    switch (Solve_which_system)
    {
        // If we are solving the original system
      case Block_J:
      {
        // get the raw residuals and jacobian
        elem_pt->get_jacobian(residuals, jacobian);
 
        // Now multiply to fill in the residuals for the final term
        for (unsigned i = 0; i < raw_ndof; i++)
        {
          unsigned local_eqn = elem_pt->eqn_number(i);
          // Add the slack parameter to the final residuals
          residuals[i] +=
            Sigma * Psi.global_value(local_eqn) / Count.global_value(local_eqn);
        }
      }
      break;
 
      // If we are solving the augmented-by-one system
      case Block_augmented_J:
      {
        // Get the full residuals, we need them
        get_residuals(elem_pt, residuals);
 
        // Need to get the raw jacobian (and raw residuals)
        Vector<double> newres(augmented_ndof);
        elem_pt->get_jacobian(newres, jacobian);
 
        // Now do finite differencing stuff
        const double FD_step = 1.0e-8;
        // Fill in the first lot of finite differences
        {
          // increase the global parameter
          double* unknown_pt = Parameter_pt;
          double init = *unknown_pt;
          *unknown_pt += FD_step;
 
          // Not do any possible updates
          Problem_pt->actions_after_change_in_bifurcation_parameter();
 
          // Get the new (modified) residuals
          get_residuals(elem_pt, newres);
 
          // The final column  is given by the difference
          // between the residuals
          for (unsigned n = 0; n < raw_ndof; n++)
          {
            jacobian(n, augmented_ndof - 1) =
              (newres[n] - residuals[n]) / FD_step;
          }
          // Reset the global parameter
          *unknown_pt = init;
 
          // Now do any possible updates
          Problem_pt->actions_after_change_in_bifurcation_parameter();
        }
 
        // Fill in the bottom row
        for (unsigned n = 0; n < raw_ndof; n++)
        {
          unsigned local_eqn = elem_pt->eqn_number(n);
          jacobian(augmented_ndof - 1, n) =
            Psi.global_value(local_eqn) / Count.global_value(local_eqn);
        }
      }
      break;
 
      // Otherwise solving the full system
      case Full_augmented:
      {
        Problem_pt->actions_before_newton_convergence_check();
        // Get the basic jacobian and residuals
        elem_pt->get_jacobian(residuals, jacobian);
        // get the proper residuals
        get_residuals(elem_pt, residuals);
 
        // Now fill in the next diagonal jacobian entry
        for (unsigned n = 0; n < raw_ndof; n++)
        {
          for (unsigned m = 0; m < raw_ndof; m++)
          {
            jacobian(raw_ndof + 1 + n, raw_ndof + 1 + m) = jacobian(n, m);
          }
          unsigned local_eqn = elem_pt->eqn_number(n);
          // Add in the sigma contribution
          jacobian(n, 2 * raw_ndof + 1) =
            Psi.global_value(local_eqn) / Count.global_value(local_eqn);
          // Symmetry constraint
          jacobian(raw_ndof, n) =
            Psi.global_value(local_eqn) / Count.global_value(local_eqn);
          // Non-zero constraint
          jacobian(2 * raw_ndof + 1, raw_ndof + 1 + n) =
            C.global_value(local_eqn) / Count.global_value(local_eqn);
        }
 
        // Finite difference the remaining blocks
        const double FD_step = 1.0e-8;
 
        Vector<double> newres_p(augmented_ndof);
 
        // Loop over the ndofs only
        for (unsigned n = 0; n < raw_ndof; ++n)
        {
          // Get the original (unaugmented) global equation number
          unsigned long global_eqn =
            Assembly_handler_pt->eqn_number(elem_pt, n);
          double* unknown_pt = Problem_pt->global_dof_pt(global_eqn);
          double init = *unknown_pt;
          *unknown_pt += FD_step;
          Problem_pt->actions_before_newton_convergence_check();
          // Get the new residuals
          get_residuals(elem_pt, newres_p);
          // Fill in the entries in the block d(Jy)/dx
          for (unsigned m = 0; m < raw_ndof; m++)
          {
            jacobian(raw_ndof + 1 + m, n) =
              (newres_p[raw_ndof + 1 + m] - residuals[raw_ndof + 1 + m]) /
              (FD_step);
          }
 
          // Reset the unknown
          *unknown_pt = init;
          Problem_pt->actions_before_newton_convergence_check();
        }
 
        {
          // Now increase the global parameter
          double* unknown_pt = Parameter_pt;
          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);
 
          // Add in the first block d/dx
          for (unsigned m = 0; m < raw_ndof; m++)
          {
            jacobian(m, raw_ndof) = (newres_p[m] - residuals[m]) / FD_step;
          }
          // Add in the second block d/dy
          for (unsigned m = raw_ndof + 1; m < augmented_ndof - 1; m++)
          {
            jacobian(m, raw_ndof) = (newres_p[m] - residuals[m]) / FD_step;
          }
 
          // Reset the unknown
          *unknown_pt = init;
          Problem_pt->actions_after_change_in_bifurcation_parameter();
        }
      }
      break;
 
      default:
        std::ostringstream error_stream;
        error_stream
          << "The Solve_which_system flag can only take values 0, 1, 2"
          << " not " << Solve_which_system << "\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    Problem_pt->actions_after_change_in_bifurcation_parameter();
  }

References oomph::Problem::actions_after_change_in_bifurcation_parameter(), oomph::Problem::actions_before_newton_convergence_check(), Assembly_handler_pt, Block_augmented_J, Block_J, Count, oomph::GeneralisedElement::eqn_number(), oomph::AssemblyHandler::eqn_number(), oomph::BlackBoxFDNewtonSolver::FD_step, Full_augmented, oomph::GeneralisedElement::get_jacobian(), get_residuals(), oomph::Problem::global_dof_pt(), oomph::DoubleVectorWithHaloEntries::global_value(), i, m, n, oomph::GeneralisedElement::ndof(), ndof(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Parameter_pt, Problem_pt, Psi, Sigma, and Solve_which_system.

Referenced by oomph::AugmentedBlockPitchForkLinearSolver::resolve(), and oomph::AugmentedBlockPitchForkLinearSolver::solve().

get_residuals()

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

virtual

Get the residuals.

Reimplemented from oomph::AssemblyHandler.

  {
    // Need to get raw residuals and jacobian
    unsigned raw_ndof = elem_pt->ndof();
 
    // Find out which system we are solving
    switch (Solve_which_system)
    {
        // If we are solving the original system
      case Block_J:
      {
        // get the basic residuals
        elem_pt->get_residuals(residuals);
        // Now multiply to fill in the residuals for the final term
        for (unsigned i = 0; i < raw_ndof; i++)
        {
          unsigned local_eqn = elem_pt->eqn_number(i);
          // Add the slack parameter to the final residuals
          residuals[i] +=
            Sigma * Psi.global_value(local_eqn) / Count.global_value(local_eqn);
        }
      }
      break;
 
      // If we are solving the augmented-by-one system
      case Block_augmented_J:
      {
        // Get the basic residuals
        elem_pt->get_residuals(residuals);
 
        // Zero the final residual
        residuals[raw_ndof] = 0.0;
        // Now multiply to fill in the residuals for the final term
        for (unsigned i = 0; i < raw_ndof; i++)
        {
          unsigned local_eqn = elem_pt->eqn_number(i);
          // Add the slack parameter to the final residuals
          residuals[i] +=
            Sigma * Psi.global_value(local_eqn) / Count.global_value(local_eqn);
          // Final term that specifies the symmetry
          residuals[raw_ndof] += ((*Problem_pt->global_dof_pt(local_eqn)) *
                                  Psi.global_value(local_eqn)) /
                                 Count.global_value(local_eqn);
        }
      }
      break;
 
      // Otherwise we are solving the fully augemented system
      case Full_augmented:
      {
        DenseMatrix<double> jacobian(raw_ndof);
 
        // Get the basic residuals and jacobian
        elem_pt->get_jacobian(residuals, jacobian);
 
        // Initialise the final residuals
        residuals[raw_ndof] = 0.0;
        residuals[2 * raw_ndof + 1] = -1.0 / Nelement;
 
        // Now multiply to fill in the residuals associated
        // with the null vector condition
        for (unsigned i = 0; i < raw_ndof; i++)
        {
          unsigned local_eqn = elem_pt->eqn_number(i);
          residuals[raw_ndof + 1 + i] = 0.0;
          for (unsigned j = 0; j < raw_ndof; j++)
          {
            unsigned local_unknown = elem_pt->eqn_number(j);
            residuals[raw_ndof + 1 + i] +=
              jacobian(i, j) * Y.global_value(local_unknown);
          }
 
          // Add the slack parameter to the governing equations
          residuals[i] +=
            Sigma * Psi.global_value(local_eqn) / Count.global_value(local_eqn);
 
          // Specify the symmetry
          residuals[raw_ndof] += ((*Problem_pt->global_dof_pt(local_eqn)) *
                                  Psi.global_value(local_eqn)) /
                                 Count.global_value(local_eqn);
          // Specify the normalisation
          residuals[2 * raw_ndof + 1] +=
            (Y.global_value(local_eqn) * C.global_value(local_eqn)) /
            Count.global_value(local_eqn);
        }
      }
      break;
 
      default:
        std::ostringstream error_stream;
        error_stream
          << "The Solve_which_system flag can only take values 0, 1, 2"
          << " not " << Solve_which_system << "\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  }

References Block_augmented_J, Block_J, Count, oomph::GeneralisedElement::eqn_number(), Full_augmented, oomph::GeneralisedElement::get_jacobian(), oomph::GeneralisedElement::get_residuals(), oomph::Problem::global_dof_pt(), oomph::DoubleVectorWithHaloEntries::global_value(), i, j, oomph::GeneralisedElement::ndof(), Nelement, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Problem_pt, Psi, Sigma, Solve_which_system, and Y.

Referenced by get_jacobian().

global_eqn_number()

unsigned oomph::PitchForkHandler::global_eqn_number ( const unsigned &  i )

inlineprivate

Function that is used to return map the global equations using the simplistic numbering scheme into the actual distributed scheme

    {
#ifdef OOMPH_HAS_MPI
      // If the problem is distributed I have to do something
      if (Distributed)
      {
        return Global_eqn_number[i];
      }
      // Otherwise it's just i
      else
#endif
      {
        return i;
      }
    }

References Global_eqn_number, and i.

Referenced by eqn_number().

ndof()

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

virtual

Get the number of elemental degrees of freedom.

Reimplemented from oomph::AssemblyHandler.

  {
    unsigned raw_ndof = elem_pt->ndof();
 
    // Return different values depending on the type of block decomposition
    switch (Solve_which_system)
    {
      case Full_augmented:
        return (2 * raw_ndof + 2);
        break;
 
      case Block_augmented_J:
        return (raw_ndof + 1);
        break;
 
      case Block_J:
        return raw_ndof;
        break;
 
      default:
        std::ostringstream error_stream;
        error_stream
          << "The Solve_which_system flag can only take values 0, 1, 2"
          << " not " << Solve_which_system << "\n";
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  }

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

Referenced by get_jacobian(), oomph::AugmentedBlockPitchForkLinearSolver::resolve(), and oomph::AugmentedBlockPitchForkLinearSolver::solve().

solve_augmented_block_system()

void oomph::PitchForkHandler::solve_augmented_block_system

(

)

Set to solve the augmented block system.

Set to solve the augmented-by-one block system.

  {
    // Only bother to do anything if we haven't already set the flag
    if (Solve_which_system != Block_augmented_J)
    {
      // If we were solving the system with the original jacobian add the
      // parameter back
      if (Solve_which_system == Block_J)
      {
        Problem_pt->Dof_pt.push_back(Parameter_pt);
      }
 
      // Restrict the problem to the standard variables and
      // the bifurcation parameter only
      Problem_pt->Dof_distribution_pt = new LinearAlgebraDistribution(
        Problem_pt->communicator_pt(), Ndof + 1, false);
      Problem_pt->Dof_pt.resize(Ndof + 1);
 
      // Problem_pt->Dof_pt.resize(Ndof+1);
      //    Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
      //                                         Ndof+1,false);
 
      // Remove all previous sparse storage used during Jacobian assembly
      Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
 
      // Set the solve flag
      Solve_which_system = Block_augmented_J;
    }
  }

References Block_augmented_J, Block_J, oomph::Problem::communicator_pt(), oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, Ndof, Parameter_pt, Problem_pt, Solve_which_system, and oomph::Problem::Sparse_assemble_with_arrays_previous_allocation.

Referenced by oomph::AugmentedBlockPitchForkLinearSolver::resolve(), and oomph::AugmentedBlockPitchForkLinearSolver::solve().

solve_block_system()

void oomph::PitchForkHandler::solve_block_system

(

)

Set to solve the block system.

Set to solve the block system. The boolean flag is used to specify whether the block decomposition should use exactly the same jacobian as the original system.

  {
    // Only bother to do anything if we haven't already set the flag
    if (Solve_which_system != Block_J)
    {
      // Restrict the problem to the standard variables
      Problem_pt->Dof_pt.resize(this->Dof_distribution_pt->nrow_local());
      Problem_pt->Dof_distribution_pt = this->Dof_distribution_pt;
 
      // 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);
 
      // Set the solve flag
      Solve_which_system = Block_J;
    }
  }

References Block_J, Dof_distribution_pt, oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, oomph::LinearAlgebraDistribution::nrow_local(), Problem_pt, Solve_which_system, and oomph::Problem::Sparse_assemble_with_arrays_previous_allocation.

Referenced by oomph::BlockPitchForkLinearSolver::resolve(), and oomph::BlockPitchForkLinearSolver::solve().

solve_full_system()

void oomph::PitchForkHandler::solve_full_system

(

)

Solve non-block system.

  {
    // Only do something if we are not solving the full system
    if (Solve_which_system != Full_augmented)
    {
#ifdef OOMPH_HAS_MPI
      int my_rank = Problem_pt->communicator_pt()->my_rank();
#endif
      // If we were solving the problem without any augementation
      // add the parameter again
      if (Solve_which_system == Block_J)
      {
#ifdef OOMPH_HAS_MPI
 
        if ((!Distributed) || (my_rank == 0))
#endif
        {
          Problem_pt->Dof_pt.push_back(Parameter_pt);
        }
      }
 
      // Always add the additional unknowns back into the problem
      const unsigned n_dof_local = Dof_distribution_pt->nrow_local();
      for (unsigned n = 0; n < n_dof_local; n++)
      {
        Problem_pt->Dof_pt.push_back(&Y[n]);
      }
 
 
#ifdef OOMPH_HAS_MPI
      if ((!Distributed) || (my_rank == 0))
#endif
      // Add the slack parameter to the problem on the first processor
      // if distributed
      {
        Problem_pt->Dof_pt.push_back(&Sigma);
      }
 
 
      // Delete the distribtion created in the augmented block solve
      if (Problem_pt->Dof_distribution_pt != this->Dof_distribution_pt)
      {
        delete Problem_pt->Dof_distribution_pt;
      }
 
      // update the Dof distribution
      Problem_pt->Dof_distribution_pt = this->Augmented_dof_distribution_pt;
      // Problem_pt->Dof_distribution_pt->build(Problem_pt->communicator_pt(),
      //                                       Ndof*2+2,false);
      // Remove all previous sparse storage used during Jacobian assembly
      Problem_pt->Sparse_assemble_with_arrays_previous_allocation.resize(0);
 
      // Set the solve flag
      Solve_which_system = Full_augmented;
    }
  }

References Augmented_dof_distribution_pt, Block_J, oomph::Problem::communicator_pt(), Dof_distribution_pt, oomph::Problem::Dof_distribution_pt, oomph::Problem::Dof_pt, Full_augmented, oomph::OomphCommunicator::my_rank(), n, oomph::LinearAlgebraDistribution::nrow_local(), Parameter_pt, Problem_pt, Sigma, Solve_which_system, oomph::Problem::Sparse_assemble_with_arrays_previous_allocation, and Y.

Referenced by oomph::BlockPitchForkLinearSolver::resolve(), oomph::AugmentedBlockPitchForkLinearSolver::resolve(), oomph::BlockPitchForkLinearSolver::solve(), and oomph::AugmentedBlockPitchForkLinearSolver::solve().

Friends And Related Function Documentation

AugmentedBlockPitchForkLinearSolver

friend class AugmentedBlockPitchForkLinearSolver

friend

BlockPitchForkLinearSolver

friend class BlockPitchForkLinearSolver

friend

Member Data Documentation

Assembly_handler_pt

AssemblyHandler* oomph::PitchForkHandler::Assembly_handler_pt

private

Pointer to the underlying (original) assembly handler.

Referenced by get_jacobian(), and PitchForkHandler().

Augmented_dof_distribution_pt

LinearAlgebraDistribution* oomph::PitchForkHandler::Augmented_dof_distribution_pt

private

The augmented distribution.

Referenced by PitchForkHandler(), oomph::BlockPitchForkLinearSolver::resolve(), oomph::BlockPitchForkLinearSolver::solve(), and solve_full_system().

C

DoubleVectorWithHaloEntries oomph::PitchForkHandler::C

private

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

Count

DoubleVectorWithHaloEntries oomph::PitchForkHandler::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. This should really be an integer, but its double so that the distribution can be used

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

Dof_distribution_pt

LinearAlgebraDistribution* oomph::PitchForkHandler::Dof_distribution_pt

private

Store the original dof distribution.

Referenced by get_eigenfunction(), PitchForkHandler(), oomph::BlockPitchForkLinearSolver::resolve(), oomph::BlockPitchForkLinearSolver::solve(), solve_block_system(), solve_full_system(), and ~PitchForkHandler().

Global_eqn_number

Vector<unsigned> oomph::PitchForkHandler::Global_eqn_number

private

A vector that is used to map the global equations to their actual location in a distributed problem

Referenced by global_eqn_number(), and PitchForkHandler().

Ndof

unsigned oomph::PitchForkHandler::Ndof

private

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

Referenced by eqn_number(), PitchForkHandler(), and solve_augmented_block_system().

Nelement

unsigned oomph::PitchForkHandler::Nelement

private

Referenced by get_residuals(), and PitchForkHandler().

Parameter_pt

double* oomph::PitchForkHandler::Parameter_pt

private

Storage for the pointer to the parameter.

Referenced by bifurcation_parameter_pt(), get_jacobian(), oomph::BlockPitchForkLinearSolver::solve(), solve_augmented_block_system(), and solve_full_system().

Problem_pt

Problem* oomph::PitchForkHandler::Problem_pt

private

Pointer to the problem.

Referenced by get_dresiduals_dparameter(), get_jacobian(), get_residuals(), PitchForkHandler(), solve_augmented_block_system(), solve_block_system(), solve_full_system(), and ~PitchForkHandler().

Psi

DoubleVectorWithHaloEntries oomph::PitchForkHandler::Psi

private

A constant vector that is specifies the symmetry being broken.

Referenced by get_jacobian(), get_residuals(), PitchForkHandler(), oomph::BlockPitchForkLinearSolver::resolve(), oomph::BlockPitchForkLinearSolver::solve(), and oomph::AugmentedBlockPitchForkLinearSolver::solve().

Sigma

double oomph::PitchForkHandler::Sigma

private

A slack variable used to specify the amount of antisymmetry in the solution.

Referenced by get_jacobian(), get_residuals(), PitchForkHandler(), and solve_full_system().

Solve_which_system

unsigned oomph::PitchForkHandler::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), a system in which the jacobian is augmented by 1 row and column to ensure that it is non-singular (2). See the enum above

Referenced by eqn_number(), get_dresiduals_dparameter(), get_jacobian(), get_residuals(), ndof(), solve_augmented_block_system(), solve_block_system(), and solve_full_system().

Y

DoubleVectorWithHaloEntries oomph::PitchForkHandler::Y

private

Storage for the null vector.

Referenced by get_dresiduals_dparameter(), get_eigenfunction(), get_hessian_vector_products(), get_residuals(), PitchForkHandler(), oomph::BlockPitchForkLinearSolver::resolve(), oomph::AugmentedBlockPitchForkLinearSolver::resolve(), oomph::BlockPitchForkLinearSolver::solve(), oomph::AugmentedBlockPitchForkLinearSolver::solve(), and solve_full_system().

---

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

• assembly_handler.h
• assembly_handler.cc