#include <complex_smoother.h>

Inheritance diagram for oomph::ComplexDampedJacobi< MATRIX >:

Public Member Functions
	ComplexDampedJacobi (const double &omega=0.5)
	Constructor (empty) More...

	~ComplexDampedJacobi ()
	Empty destructor. More...

void	clean_up_memory ()
	Cleanup data that's stored for resolve (if any has been stored) More...

	ComplexDampedJacobi (const ComplexDampedJacobi &)=delete
	Broken copy constructor. More...

void	operator= (const ComplexDampedJacobi &)=delete
	Broken assignment operator. More...

void	calculate_omega (const double &k, const double &h)

double &	omega ()
	Get access to the value of Omega (lvalue) More...

void	complex_smoother_setup (Vector< CRDoubleMatrix * > helmholtz_matrix_pt)
	Setup: Pass pointer to the matrix and store in cast form. More...

void	complex_smoother_solve (const Vector< DoubleVector > &rhs, Vector< DoubleVector > &solution)

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

unsigned	iterations () const
	Number of iterations taken. More...

Public Member Functions inherited from oomph::HelmholtzSmoother
	HelmholtzSmoother ()
	Empty constructor. More...

virtual	~HelmholtzSmoother ()
	Virtual empty destructor. More...

void	complex_matrix_multiplication (Vector< CRDoubleMatrix * > matrices_pt, const Vector< DoubleVector > &x, Vector< DoubleVector > &soln)

template<typename MATRIX >
void	check_validity_of_solve_helper_inputs (CRDoubleMatrix const &real_matrix_pt, CRDoubleMatrix const &imag_matrix_pt, const Vector< DoubleVector > &rhs, Vector< DoubleVector > &solution, const double &n_dof)

Public Member Functions inherited from oomph::IterativeLinearSolver
	IterativeLinearSolver ()

	IterativeLinearSolver (const IterativeLinearSolver &)=delete
	Broken copy constructor. More...

void	operator= (const IterativeLinearSolver &)=delete
	Broken assignment operator. More...

virtual	~IterativeLinearSolver ()
	Destructor (empty) More...

Preconditioner *&	preconditioner_pt ()
	Access function to preconditioner. More...

Preconditioner *const &	preconditioner_pt () const
	Access function to preconditioner (const version) More...

double &	tolerance ()
	Access to convergence tolerance. More...

unsigned &	max_iter ()
	Access to max. number of iterations. More...

void	enable_doc_convergence_history ()
	Enable documentation of the convergence history. More...

void	disable_doc_convergence_history ()
	Disable documentation of the convergence history. More...

void	open_convergence_history_file_stream (const std::string &file_name, const std::string &zone_title="")

void	close_convergence_history_file_stream ()
	Close convergence history output stream. More...

double	jacobian_setup_time () const

double	linear_solver_solution_time () const
	return the time taken to solve the linear system More...

virtual double	preconditioner_setup_time () const
	returns the the time taken to setup the preconditioner More...

void	enable_setup_preconditioner_before_solve ()
	Setup the preconditioner before the solve. More...

void	disable_setup_preconditioner_before_solve ()
	Don't set up the preconditioner before the solve. More...

void	enable_error_after_max_iter ()
	Throw an error if we don't converge within max_iter. More...

void	disable_error_after_max_iter ()
	Don't throw an error if we don't converge within max_iter (default). More...

void	enable_iterative_solver_as_preconditioner ()

void	disable_iterative_solver_as_preconditioner ()

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 (DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)

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

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 (const DoubleVector &rhs, DoubleVector &result)

virtual void	resolve_transpose (const DoubleVector &rhs, DoubleVector &result)

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 Member Functions
void	complex_solve_helper (const Vector< DoubleVector > &rhs, Vector< DoubleVector > &solution)
	This is where the actual work is done. More...

Private Attributes
bool	Matrix_can_be_deleted

CRDoubleMatrix *	Matrix_real_pt
	Pointer to the real part of the system matrix. More...

CRDoubleMatrix *	Matrix_imag_pt
	Pointer to the real part of the system matrix. More...

Vector< double >	Matrix_diagonal_real
	Vector containing the diagonal entries of A_r (real(A)) More...

Vector< double >	Matrix_diagonal_imag
	Vector containing the diagonal entries of A_c (imag(A)) More...

Vector< double >	Matrix_diagonal_inv_sq
	Vector whose j'th entry is given by 1/(A_r(j,j)^2 + A_c(j,j)^2) More...

unsigned	Iterations
	Number of iterations taken. More...

double	Omega
	Damping factor. More...

Additional Inherited Members
Protected Member Functions inherited from oomph::DistributableLinearAlgebraObject
void	clear_distribution ()

Protected Attributes inherited from oomph::HelmholtzSmoother
bool	Use_as_smoother

Protected Attributes inherited from oomph::IterativeLinearSolver
bool	Doc_convergence_history

std::ofstream	Output_file_stream
	Output file stream for convergence history. More...

double	Tolerance
	Convergence tolerance. More...

unsigned	Max_iter
	Maximum number of iterations. More...

Preconditioner *	Preconditioner_pt
	Pointer to the preconditioner. More...

double	Jacobian_setup_time
	Jacobian setup time. More...

double	Solution_time
	linear solver solution time More...

double	Preconditioner_setup_time
	Preconditioner setup time. More...

bool	Setup_preconditioner_before_solve

bool	Throw_error_after_max_iter

bool	Use_iterative_solver_as_preconditioner
	Use the iterative solver as preconditioner. More...

bool	First_time_solve_when_used_as_preconditioner

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

Static Protected Attributes inherited from oomph::IterativeLinearSolver
static IdentityPreconditioner	Default_preconditioner

Detailed Description

template<typename MATRIX>
class oomph::ComplexDampedJacobi< MATRIX >

//////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////// Damped Jacobi "solver" templated by matrix type. The "solver" exists in many different incarnations: It's an IterativeLinearSolver, and a Smoother, all of which use the same basic iteration.

Constructor & Destructor Documentation

◆ ComplexDampedJacobi() [1/2]

template<typename MATRIX >

oomph::ComplexDampedJacobi< MATRIX >::ComplexDampedJacobi ( const double & omega = 0.5 )

inline

Constructor (empty)

       : Matrix_can_be_deleted(true),
         Matrix_real_pt(0),
         Matrix_imag_pt(0),
         Omega(omega){};

◆ ~ComplexDampedJacobi()

template<typename MATRIX >

oomph::ComplexDampedJacobi< MATRIX >::~ComplexDampedJacobi ( )

inline

Empty destructor.

     {
       // Run clean_up_memory
       clean_up_memory();
     } // End of ~ComplexDampedJacobi

References oomph::ComplexDampedJacobi< MATRIX >::clean_up_memory().

◆ ComplexDampedJacobi() [2/2]

template<typename MATRIX >

oomph::ComplexDampedJacobi< MATRIX >::ComplexDampedJacobi ( const ComplexDampedJacobi< MATRIX > & )

delete

Broken copy constructor.

Member Function Documentation

◆ calculate_omega()

template<typename MATRIX >

void oomph::ComplexDampedJacobi< MATRIX >::calculate_omega	(	const double &	k,
		const double &	h
	)

inline

Function to calculate the value of Omega by passing in the value of k and h [see Elman et al. "A multigrid method enhanced by Krylov subspace iteration for discrete Helmholtz equations", p.1303]

     {
       // Create storage for the parameter kh
       double kh = k * h;
  
       // Store the value of pi
       double pi = MathematicalConstants::Pi;
  
       // Calculate the value of Omega
       double omega = (12.0 - 4.0 * pow(kh, 2.0)) / (18.0 - 3.0 * pow(kh, 2.0));
  
       // Set the tolerance for how close omega can be to 0
       double tolerance = 1.0e-03;
  
       // Only store the value of omega if it lies in the range (0,1]. If it
       // isn't it will produce a divergent scheme. Note, to avoid letting
       // omega become too small we make sure it is greater than some tolerance
       if ((omega > tolerance) && !(omega > 1))
       {
         // Since omega lies in the specified range, store it
         Omega = omega;
       }
       // On the coarsest grids: if the wavenumber is greater than pi and
       // kh>2cos(pi*h/2) then we can choose omega from the range (0,omega_2)
       // where omega_2 is defined in [Elman et al."A multigrid method
       // enhanced by Krylov subspace iteration for discrete Helmholtz
       // equations", p.1295]
       else if ((k > pi) && (kh > 2 * cos(pi * h / 2)))
       {
         // Calculate the numerator of omega_2
         double omega_2 = (2.0 - pow(kh, 2.0));
  
         // Divide by the denominator of the fraction
         omega_2 /= (2.0 * pow(sin(pi * h / 2), 2.0) - 0.5 * pow(kh, 2.0));
  
         // If 2/3 lies in the range (0,omega_2), use it
         if (omega_2 > (2.0 / 3.0))
         {
           // Assign Omega the damping factor used for the Poisson equation
           Omega = 2.0 / 3.0;
         }
         // If omega_2 is less than 2/3 use the midpoint of (tolerance,omega_2)
         else
         {
           // Set the value of Omega
           Omega = 0.5 * (tolerance + omega_2);
         }
       }
       // Since the value of kh must be fairly large, make the value of
       // omega small to compensate
       else
       {
         // Since kh doesn't lie in the chosen range, set it to some small value
         Omega = 0.2;
       }
     } // End of calculate_omega

References cos(), k, oomph::ComplexDampedJacobi< MATRIX >::omega(), oomph::ComplexDampedJacobi< MATRIX >::Omega, constants::pi, oomph::MathematicalConstants::Pi, Eigen::bfloat16_impl::pow(), sin(), and oomph::IterativeLinearSolver::tolerance().

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::setup_smoothers().

◆ clean_up_memory()

template<typename MATRIX >

void oomph::ComplexDampedJacobi< MATRIX >::clean_up_memory ( )

inlinevirtual

Cleanup data that's stored for resolve (if any has been stored)

Reimplemented from oomph::LinearSolver.

     {
       // If the real matrix pointer isn't null AND we're allowed to delete
       // the matrix which is only when we create the matrix ourselves
       if ((Matrix_real_pt != 0) && (Matrix_can_be_deleted))
       {
         // Delete the matrix
         delete Matrix_real_pt;
  
         // Assign the associated pointer the value NULL
         Matrix_real_pt = 0;
       }
  
       // If the real matrix pointer isn't null AND we're allowed to delete
       // the matrix which is only when we create the matrix ourselves
       if ((Matrix_imag_pt != 0) && (Matrix_can_be_deleted))
       {
         // Delete the matrix
         delete Matrix_imag_pt;
  
         // Assign the associated pointer the value NULL
         Matrix_imag_pt = 0;
       }
     } // End of clean_up_memory

References oomph::ComplexDampedJacobi< MATRIX >::Matrix_can_be_deleted, oomph::ComplexDampedJacobi< MATRIX >::Matrix_imag_pt, and oomph::ComplexDampedJacobi< MATRIX >::Matrix_real_pt.

Referenced by oomph::ComplexDampedJacobi< MATRIX >::~ComplexDampedJacobi().

◆ complex_smoother_setup()

template<typename MATRIX >

void oomph::ComplexDampedJacobi< MATRIX >::complex_smoother_setup ( Vector< CRDoubleMatrix * > helmholtz_matrix_pt )

inlinevirtual

Setup: Pass pointer to the matrix and store in cast form.

Implements oomph::HelmholtzSmoother.

     {
       // Assume the matrices have been passed in from the outside so we must
       // not delete it. This is needed to avoid pre- and post-smoothers
       // deleting the same matrix in the MG solver. If this was originally
       // set to TRUE then this will be sorted out in the other functions
       // from which this was called
       Matrix_can_be_deleted = false;
  
       // Assign the real matrix pointers
       Matrix_real_pt = helmholtz_matrix_pt[0];
  
       // Assign the imaginary matrix pointers
       Matrix_imag_pt = helmholtz_matrix_pt[1];
  
 #ifdef PARANOID
       // PARANOID check that if the matrix is distributable. If it is not then
       // it should not be distributed
       if (Matrix_real_pt->nrow() != Matrix_imag_pt->nrow())
       {
         std::ostringstream error_message_stream;
         error_message_stream << "Incompatible real and complex matrix sizes.";
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
 #endif
  
       //--------------------------------------------------------------------
       // We need the matrix diagonals to calculate inv(D) (where D is the
       // diagonal of A) as it remains the same for each use of the iterative
       // scheme. To avoid unnecessary computations we store it now so it can
       // simply be called in each iteration.
       //--------------------------------------------------------------------
  
       // Grab the diagonal entries of the real part of the system matrix
       Matrix_diagonal_real = Matrix_real_pt->diagonal_entries();
  
       // Grab the diagonal entries of the imaginary part of the system matrix
       Matrix_diagonal_imag = Matrix_imag_pt->diagonal_entries();
  
       // Find the number of entries in the vectors containing the diagonal
       // entries of both the real and the imaginary matrix
       unsigned n_dof = Matrix_diagonal_real.size();
  
       // Make a dummy vector to store the entries of Matrix_diagonal_inv_sq
       Matrix_diagonal_inv_sq.resize(n_dof, 0.0);
  
       // Create a helper variable to store A_r(j,j), for some j
       double dummy_r;
  
       // Create a helper variable to store A_c(j,j), for some j
       double dummy_c;
  
       // Loop over the entries of Matrix_diagonal_inv_sq and set the i-th
       // entry to 1/(A_r(i,i)^2 + A_c(i,i)^2)
       for (unsigned i = 0; i < n_dof; i++)
       {
         // Store the value of A_r(i,i)
         dummy_r = Matrix_diagonal_real[i];
  
         // Store the value of A_c(i,i)
         dummy_c = Matrix_diagonal_imag[i];
  
         // Store the value of 1/(A_r(i,i)^2 + A_c(i,i)^2)
         Matrix_diagonal_inv_sq[i] =
           1.0 / (dummy_r * dummy_r + dummy_c * dummy_c);
       }
     } // End of complex_smoother_setup

References oomph::CRDoubleMatrix::diagonal_entries(), i, oomph::ComplexDampedJacobi< MATRIX >::Matrix_can_be_deleted, oomph::ComplexDampedJacobi< MATRIX >::Matrix_diagonal_imag, oomph::ComplexDampedJacobi< MATRIX >::Matrix_diagonal_inv_sq, oomph::ComplexDampedJacobi< MATRIX >::Matrix_diagonal_real, oomph::ComplexDampedJacobi< MATRIX >::Matrix_imag_pt, oomph::ComplexDampedJacobi< MATRIX >::Matrix_real_pt, oomph::CRDoubleMatrix::nrow(), OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ complex_smoother_solve()

template<typename MATRIX >

void oomph::ComplexDampedJacobi< MATRIX >::complex_smoother_solve	(	const Vector< DoubleVector > &	rhs,
		Vector< DoubleVector > &	solution
	)

inlinevirtual

The smoother_solve function performs fixed number of iterations on the system A*result=rhs. The number of (smoothing) iterations is the same as the max. number of iterations in the underlying IterativeLinearSolver class.

Implements oomph::HelmholtzSmoother.

     {
       // If you use a smoother but you don't want to calculate the residual
       Use_as_smoother = true;
  
       // Call the helper function
       complex_solve_helper(rhs, solution);
     } // End of complex_smoother_solve

References oomph::ComplexDampedJacobi< MATRIX >::complex_solve_helper(), BiharmonicTestFunctions1::solution(), and oomph::HelmholtzSmoother::Use_as_smoother.

◆ complex_solve_helper()

template<typename MATRIX >

void oomph::ComplexDampedJacobi< MATRIX >::complex_solve_helper	(	const Vector< DoubleVector > &	rhs,
		Vector< DoubleVector > &	solution
	)

private

This is where the actual work is done.

\[ \|a\|_2^2=\sum_{i=0}^{i=n-1}\real(a[i])^2+\imag(a[i])^2. \]

\[ \|a\|_2^2=\|\real(a)\|_2^2+\|\imag(a)\|_2^2. \]

   {
     // Get number of dofs
     unsigned n_dof = Matrix_real_pt->nrow();
  
 #ifdef PARANOID
     // Upcast the matrix to the appropriate type
     CRDoubleMatrix* tmp_rmatrix_pt =
       dynamic_cast<CRDoubleMatrix*>(Matrix_real_pt);
  
     // Upcast the matrix to the appropriate type
     CRDoubleMatrix* tmp_imatrix_pt =
       dynamic_cast<CRDoubleMatrix*>(Matrix_imag_pt);
  
     // PARANOID Run the self-tests to check the inputs are correct
     this->check_validity_of_solve_helper_inputs<MATRIX>(
       tmp_rmatrix_pt, tmp_imatrix_pt, rhs, solution, n_dof);
  
     // We don't need the real matrix pointer any more
     tmp_rmatrix_pt = 0;
  
     // We don't need the imaginary matrix pointer any more
     tmp_imatrix_pt = 0;
 #endif
  
     // Setup the solution if it is not
     if ((!solution[0].distribution_pt()->built()) ||
         (!solution[1].distribution_pt()->built()))
     {
       solution[0].build(rhs[0].distribution_pt(), 0.0);
       solution[1].build(rhs[1].distribution_pt(), 0.0);
     }
  
     // Initialise timer
     double t_start = TimingHelpers::timer();
  
     // Copy the real and imaginary part of the solution vector
     DoubleVector x_real(solution[0]);
     DoubleVector x_imag(solution[1]);
  
     // Copy the real and imaginary part of the RHS vector
     DoubleVector rhs_real(rhs[0]);
     DoubleVector rhs_imag(rhs[1]);
  
     // Create storage for the real and imaginary part of the constant term
     DoubleVector constant_term_real(this->distribution_pt(), 0.0);
     DoubleVector constant_term_imag(this->distribution_pt(), 0.0);
  
     // Create storage for the real and imaginary part of the residual vector.
     // These aren't used/built if we're inside the multigrid solver
     DoubleVector residual_real;
     DoubleVector residual_imag;
  
     // Create storage for the norm of the real and imaginary parts of the
     // residual vector. These aren't used if we're inside the multigrid
     // solver
     double res_real_norm = 0.0;
     double res_imag_norm = 0.0;
  
     // Variable to hold the current residual norm. This isn't used if
     // we're inside the multigrid solver
     double norm_res = 0.0;
  
     // Variables to hold the initial residual norm. This isn't used if
     // we're inside the multigrid solver
     double norm_f = 0.0;
  
     // Initialise the value of Iterations
     Iterations = 0;
  
     //------------------------------------------------------------------------
     // Initial guess isn't necessarily zero (restricted solution from finer
     // grids) therefore both x and the residual need to be assigned values
     // from inputs. The constant term at the end of the stationary iteration
     // is also calculated here since it does not change at all throughout the
     // iterative process:
     //------------------------------------------------------------------------
  
     // Loop over all the entries of each vector
     for (unsigned i = 0; i < n_dof; i++)
     {
       // Calculate the numerator of the i'th entry in the real component of
       // the constant term
       constant_term_real[i] = (Matrix_diagonal_real[i] * rhs_real[i] +
                                Matrix_diagonal_imag[i] * rhs_imag[i]);
  
       // Divide by the denominator
       constant_term_real[i] *= Matrix_diagonal_inv_sq[i];
  
       // Scale by the damping factor
       constant_term_real[i] *= Omega;
  
       // Calculate the numerator of the i'th entry in the imaginary component of
       // the constant term
       constant_term_imag[i] = (Matrix_diagonal_real[i] * rhs_imag[i] -
                                Matrix_diagonal_imag[i] * rhs_real[i]);
  
       // Divide by the denominator
       constant_term_imag[i] *= Matrix_diagonal_inv_sq[i];
  
       // Scale by the damping factor
       constant_term_imag[i] *= Omega;
     }
  
     // Create 4 temporary vectors to store the various matrix-vector products
     // required. The appropriate combination has been appended to the name. For
     // instance, the product A_r*x_c corresponds to temp_vec_rc (real*imag)
     DoubleVector temp_vec_rr(this->distribution_pt(), 0.0);
     DoubleVector temp_vec_cc(this->distribution_pt(), 0.0);
     DoubleVector temp_vec_rc(this->distribution_pt(), 0.0);
     DoubleVector temp_vec_cr(this->distribution_pt(), 0.0);
  
     // Calculate the different combinations of A*x (or A*e depending on the
     // level in the heirarchy) in the complex case (i.e. A_r*x_r, A_c*x_c,
     // A_r*x_c and A_c*x_r)
     Matrix_real_pt->multiply(x_real, temp_vec_rr);
     Matrix_imag_pt->multiply(x_imag, temp_vec_cc);
     Matrix_real_pt->multiply(x_imag, temp_vec_rc);
     Matrix_imag_pt->multiply(x_real, temp_vec_cr);
  
     //---------------------------------------------------------------------------
     // Calculating the residual r=b-Ax in the complex case requires more care
     // than the real case. The correct calculation can be derived rather easily.
     // Consider the splitting of A, x and b into their complex components:
     //           r = b - A*x
     //             = (b_r + i*b_c) - (A_r + i*A_c)*(x_r + i*x_c)
     //             = [b_r - A_r*x_r + A_c*x_c] + i*[b_c - A_r*x_c - A_c*x_r]
     // ==> real(r) = b_r - A_r*x_r + A_c*x_c
     //   & imag(r) = b_c - A_r*x_c - A_c*x_r.
     //---------------------------------------------------------------------------
  
     // Calculate the residual only if we're not inside the multigrid solver
     if (!Use_as_smoother)
     {
       // Create storage for the real and imaginary part of the residual vector
       residual_real.build(this->distribution_pt(), 0.0);
       residual_imag.build(this->distribution_pt(), 0.0);
  
       // Real part of the residual:
       residual_real = rhs_real;
       residual_real -= temp_vec_rr;
       residual_real += temp_vec_cc;
  
       // Imaginary part of the residual:
       residual_imag = rhs_imag;
       residual_imag -= temp_vec_rc;
       residual_imag -= temp_vec_cr;
  
       // Calculate the 2-norm (noting that the 2-norm of a complex vector a of
       // length n is simply the square root of the sum of the squares of each
       // real and imaginary component). That is:
       // can be written as:
       double res_real_norm = residual_real.norm();
       double res_imag_norm = residual_imag.norm();
       double norm_res =
         sqrt(res_real_norm * res_real_norm + res_imag_norm * res_imag_norm);
  
       // If required will document convergence history to screen
       // or file (if stream is open)
       if (Doc_convergence_history)
       {
         if (!Output_file_stream.is_open())
         {
           oomph_info << Iterations << " " << norm_res << std::endl;
         }
         else
         {
           Output_file_stream << Iterations << " " << norm_res << std::endl;
         }
       } // if (Doc_convergence_history)
     } // if (!Use_as_smoother)
  
     // Two temporary vectors to store the value of A_r*x_r - A_c*x_c and
     // A_c*x_r + A_r*x_c in each iteration
     DoubleVector temp_vec_real(this->distribution_pt(), 0.0);
     DoubleVector temp_vec_imag(this->distribution_pt(), 0.0);
  
     // Calculate A_r*x_r - A_c*x_c
     temp_vec_real = temp_vec_rr;
     temp_vec_real -= temp_vec_cc;
  
     // Calculate A_c*x_r + A_r*x_c
     temp_vec_imag = temp_vec_cr;
     temp_vec_imag += temp_vec_rc;
  
     // Outermost loop: Run up to Max_iter times
     for (unsigned iter_num = 0; iter_num < Max_iter; iter_num++)
     {
       // Loop over each degree of freedom and update
       // the current approximation
       for (unsigned i = 0; i < n_dof; i++)
       {
         // Make a couple of dummy variables to help with calculations
         double dummy_r = 0.0;
         double dummy_c = 0.0;
  
         // Calculate one part of the correction term (real)
         dummy_r = (Matrix_diagonal_real[i] * temp_vec_real[i] +
                    Matrix_diagonal_imag[i] * temp_vec_imag[i]);
  
         // Calculate one part of the correction term (imaginary)
         dummy_c = (Matrix_diagonal_real[i] * temp_vec_imag[i] -
                    Matrix_diagonal_imag[i] * temp_vec_real[i]);
  
         // Scale by Omega/([A(i,i)_r]^2+[A(i,i)_c]^2)
         dummy_r *= Omega * Matrix_diagonal_inv_sq[i];
         dummy_c *= Omega * Matrix_diagonal_inv_sq[i];
  
         // Update the value of x_real
         x_real[i] -= dummy_r;
         x_imag[i] -= dummy_c;
       }
  
       // Update the value of x (or e depending on the level in the heirarchy)
       x_real += constant_term_real;
       x_imag += constant_term_imag;
  
       // Calculate the different combinations of A*x (or A*e depending on the
       // level in the heirarchy) in the complex case (i.e. A_r*x_r, A_c*x_c,
       // A_r*x_c and A_c*x_r)
       Matrix_real_pt->multiply(x_real, temp_vec_rr);
       Matrix_imag_pt->multiply(x_imag, temp_vec_cc);
       Matrix_real_pt->multiply(x_imag, temp_vec_rc);
       Matrix_imag_pt->multiply(x_real, temp_vec_cr);
  
       // Calculate A_r*x_r - A_c*x_c
       temp_vec_real = temp_vec_rr;
       temp_vec_real -= temp_vec_cc;
  
       // Calculate A_c*x_r + A_r*x_c
       temp_vec_imag = temp_vec_cr;
       temp_vec_imag += temp_vec_rc;
  
       // Calculate the residual only if we're not inside the multigrid solver
       if (!Use_as_smoother)
       {
         // Calculate the residual
         residual_real = rhs_real;
         residual_real -= temp_vec_rr;
         residual_real += temp_vec_cc;
  
         // Calculate the imaginary part of the residual vector
         residual_imag = rhs_imag;
         residual_imag -= temp_vec_rc;
         residual_imag -= temp_vec_cr;
  
         // Calculate the 2-norm using the method mentioned previously
         res_real_norm = residual_real.norm();
         res_imag_norm = residual_imag.norm();
         norm_res =
           sqrt(res_real_norm * res_real_norm + res_imag_norm * res_imag_norm) /
           norm_f;
  
         // If required, this will document convergence history to
         // screen or file (if the stream is open)
         if (Doc_convergence_history)
         {
           if (!Output_file_stream.is_open())
           {
             oomph_info << Iterations << " " << norm_res << std::endl;
           }
           else
           {
             Output_file_stream << Iterations << " " << norm_res << std::endl;
           }
         } // if (Doc_convergence_history)
  
         // Check the tolerance only if the residual norm is being computed
         if (norm_res < Tolerance)
         {
           // Break out of the for-loop
           break;
         }
       } // if (!Use_as_smoother)
     } // for (unsigned iter_num=0;iter_num<Max_iter;iter_num++)
  
     // Calculate the residual only if we're not inside the multigrid solver
     if (!Use_as_smoother)
     {
       // If time documentation is enabled
       if (Doc_time)
       {
         oomph_info << "\n(Complex) damped Jacobi converged. Residual norm: "
                    << norm_res
                    << "\nNumber of iterations to convergence: " << Iterations
                    << "\n"
                    << std::endl;
       }
     } // if (!Use_as_smoother)
  
     // Copy the solution into the solution vector
     solution[0] = x_real;
     solution[1] = x_imag;
  
     // Doc time for solver
     double t_end = TimingHelpers::timer();
     Solution_time = t_end - t_start;
     if (Doc_time)
     {
       oomph_info << "Time for solve with (complex) damped Jacobi [sec]: "
                  << Solution_time << std::endl;
     }
  
     // If the solver failed to converge and the user asked for an error if
     // this happened
     if ((Iterations > Max_iter - 1) && (Throw_error_after_max_iter))
     {
       std::string error_message =
         "Solver failed to converge and you requested ";
       error_message += "an error on convergence failures.";
       throw OomphLibError(
         error_message, OOMPH_EXCEPTION_LOCATION, OOMPH_CURRENT_FUNCTION);
     }
   } // End of complex_solve_helper function

References oomph::DoubleVector::build(), i, Global_Variables::Iterations, oomph::BlackBoxFDNewtonSolver::Max_iter, oomph::DoubleVector::norm(), oomph::SarahBL::Omega, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, BiharmonicTestFunctions1::solution(), sqrt(), oomph::Global_string_for_annotation::string(), and oomph::TimingHelpers::timer().

Referenced by oomph::ComplexDampedJacobi< MATRIX >::complex_smoother_solve().

◆ iterations()

template<typename MATRIX >

unsigned oomph::ComplexDampedJacobi< MATRIX >::iterations ( ) const

inlinevirtual

Number of iterations taken.

Implements oomph::IterativeLinearSolver.

     {
       return Iterations;
     }

References oomph::ComplexDampedJacobi< MATRIX >::Iterations.

◆ omega()

template<typename MATRIX >

double& oomph::ComplexDampedJacobi< MATRIX >::omega ( )

inline

Get access to the value of Omega (lvalue)

     {
       // Return the value of Omega
       return Omega;
     } // End of omega

References oomph::ComplexDampedJacobi< MATRIX >::Omega.

Referenced by oomph::ComplexDampedJacobi< MATRIX >::calculate_omega().

◆ operator=()

template<typename MATRIX >

void oomph::ComplexDampedJacobi< MATRIX >::operator= ( const ComplexDampedJacobi< MATRIX > & )

delete

Broken assignment operator.

◆ solve()

template<typename MATRIX >

void oomph::ComplexDampedJacobi< MATRIX >::solve	(	Problem *const &	problem_pt,
		DoubleVector &	result
	)

inlinevirtual

Use damped Jacobi iteration as an IterativeLinearSolver: This obtains the Jacobian matrix J and the residual vector r (needed for the Newton method) from the problem's get_jacobian function and returns the result of Jx=r.

Implements oomph::LinearSolver.

     {
       BrokenCopy::broken_assign("ComplexDampedJacobi");
     }

Damping factor.

Referenced by oomph::ComplexDampedJacobi< MATRIX >::calculate_omega(), and oomph::ComplexDampedJacobi< MATRIX >::omega().

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

complex_smoother.h

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

template<typename MATRIX> class oomph::ComplexDampedJacobi< MATRIX >

Constructor & Destructor Documentation

◆ ComplexDampedJacobi() [1/2]

◆ ~ComplexDampedJacobi()

◆ ComplexDampedJacobi() [2/2]

Member Function Documentation

◆ calculate_omega()

◆ clean_up_memory()

◆ complex_smoother_setup()

◆ complex_smoother_solve()

◆ complex_solve_helper()

◆ iterations()

◆ omega()

◆ operator=()

◆ solve()

Member Data Documentation

◆ Iterations

◆ Matrix_can_be_deleted

◆ Matrix_diagonal_imag

◆ Matrix_diagonal_inv_sq

◆ Matrix_diagonal_real

◆ Matrix_imag_pt

◆ Matrix_real_pt

◆ Omega

template<typename MATRIX>
class oomph::ComplexDampedJacobi< MATRIX >