Solve. More...

Public Types
using	geometry_type = Geometry

using	matrix_type = Matrix

using	container_type = Container

using	value_type = get_value_type< Container >

Public Member Functions
	MultigridCG2d ()=default
	Allocate nothing, Call `construct` method before usage. More...

template<class ... ContainerParams>
	MultigridCG2d (const Geometry &grid, const unsigned stages, ContainerParams &&... ps)
	Construct the grids and the interpolation/projection operators. More...

template<class ... Params>
void	construct (Params &&...ps)
	Perfect forward parameters to one of the constructors. More...

template<class ContainerType0 >
void	project (const ContainerType0 &src, std::vector< ContainerType0 > &out) const
	Project vector to all involved grids. More...

template<class ContainerType0 >
std::vector< ContainerType0 >	project (const ContainerType0 &src) const
	Project vector to all involved grids (allocate memory version) More...

unsigned	stages () const

unsigned	num_stages () const

const Geometry &	grid (unsigned stage) const
	return the grid at given stage More...

unsigned	max_iter () const

void	set_max_iter (unsigned new_max)
	Set the maximum number of iterations allowed at stage 0. More...

void	set_benchmark (bool benchmark, std::string message="Nested Iterations")
	Set or unset performance timings during iterations. More...

const Container &	copyable () const
	Return an object of same size as the object used for construction on the finest grid. More...

template<class MatrixType , class ContainerType0 , class ContainerType1 >
std::vector< unsigned >	solve (std::vector< MatrixType > &ops, ContainerType0 &x, const ContainerType1 &b, value_type eps)
	Nested iterations. More...

template<class MatrixType , class ContainerType0 , class ContainerType1 >
std::vector< unsigned >	solve (std::vector< MatrixType > &ops, ContainerType0 &x, const ContainerType1 &b, std::vector< value_type > eps)
	Nested iterations. More...

Detailed Description

template<class Geometry, class Matrix, class Container>
struct dg::MultigridCG2d< Geometry, Matrix, Container >

Solve.

\[ \hat O \phi = \rho \]

for self-adjoint \(\hat O\)

using dg::nested_iterations with dg::PCG solvers for any operator \(\hat O\) that is self-adjoint in appropriate weights \(W\) We refine the grids in the first two dimensions (2d / x and y)

Note: The dg::Elliptic and dg::Helmholtz classes are self-adjoint so these are the intended target operators.; The preconditioner and weights for the dg::PCG solver are taken from the precond() and weights() method in the MatrixType class

    dg::Timer t;
    t.tic();
 
    const dg::CartesianGrid2d grid( 0, lx, 0, ly, n, Nx, Ny, bcx, bcy);
 
    const unsigned stages = 3;
 
    dg::MultigridCG2d<dg::aGeometry2d, dg::DMatrix, dg::DVec > multigrid( grid, stages);
 
    const dg::DVec chi =  dg::evaluate( pol, grid);
 
    const std::vector<dg::DVec> multi_chi = multigrid.project( chi);
 
    std::vector<dg::Elliptic<dg::aGeometry2d, dg::DMatrix, dg::DVec> > multi_pol( stages);
 
    for(unsigned u=0; u<stages; u++)
    {
        multi_pol[u].construct( multigrid.grid(u), dg::centered, jfactor);
        //this tests if elliptic can recover from NaN in the preconditioner
        multi_pol[u].set_chi(0.);
        // here we test if we can set the tensor part in elliptic
        multi_pol[u].set_chi( multigrid.grid(u).metric());
        // now set the actual scalar part in chi
        multi_pol[u].set_chi( multi_chi[u]);
        multi_pol[u].set_jump_weighting(jump_weight);
    }
 
    t.toc();
 
    std::cout << "Creation of multigrid took: "<<t.diff()<<"s\n";
    const dg::DVec b =    dg::evaluate( rhs,     grid);
    dg::DVec x       =    dg::evaluate( initial, grid);
    t.tic();
    std::vector<unsigned> number = multigrid.solve(multi_pol, x, b, {eps, 1.5*eps, 1.5*eps});
    t.toc();
    std::cout << "Solution took "<< t.diff() <<"s\n";
    for( unsigned u=0; u<number.size(); u++)
        std::cout << " # iterations stage "<< number.size()-1-u << " " << number[number.size()-1-u] << " \n";

Template Parameters

Geometry	A type that is or derives from one of the abstract geometry base classes ( `aGeometry2d`, `aGeometry3d`, `aMPIGeometry2d`, ...). `Geometry` determines which `Matrix` and `Container` types can be used:
Matrix	A class for which the `dg::blas2::symv` functions are callable in connection with the `Container` class and to which the return type of `dg::create::dx()` can be converted using `dg::blas2::transfer`. The `Matrix` type can be one of: `dg::HMatrix` with `dg::HVec` and one of the shared memory geometries `dg::DMatrix` with `dg::DVec` and one of the shared memory geometries `dg::MHMatrix` with `dg::MHVec` and one of the MPI geometries `dg::MDMatrix` with `dg::MDVec` and one of the MPI geometries
Container	A data container class for which the `blas1` functionality is overloaded and to which the return type of `blas1::subroutine()` can be converted using `dg::assign`. We assume that `Container` is copyable/assignable and has a swap member function. In connection with `Geometry` this is one of `dg::HVec`, `dg::DVec` when `Geometry` is a shared memory geometry `dg::MHVec` or `dg::MDVec` when `Geometry` is one of the MPI geometries

See also: Extrapolation to generate an initial guess

Member Typedef Documentation

◆ container_type

template<class Geometry , class Matrix , class Container >

using dg::MultigridCG2d< Geometry, Matrix, Container >::container_type = Container

◆ geometry_type

template<class Geometry , class Matrix , class Container >

using dg::MultigridCG2d< Geometry, Matrix, Container >::geometry_type = Geometry

◆ matrix_type

template<class Geometry , class Matrix , class Container >

using dg::MultigridCG2d< Geometry, Matrix, Container >::matrix_type = Matrix

◆ value_type

template<class Geometry , class Matrix , class Container >

using dg::MultigridCG2d< Geometry, Matrix, Container >::value_type = get_value_type<Container>

Constructor & Destructor Documentation

◆ MultigridCG2d() [1/2]

template<class Geometry , class Matrix , class Container >

dg::MultigridCG2d< Geometry, Matrix, Container >::MultigridCG2d ( )

default

Allocate nothing, Call construct method before usage.

◆ MultigridCG2d() [2/2]

template<class Geometry , class Matrix , class Container >

template<class ... ContainerParams>

dg::MultigridCG2d< Geometry, Matrix, Container >::MultigridCG2d	(	const Geometry &	grid,
		const unsigned	stages,
		ContainerParams &&...	ps
	)

inline

Construct the grids and the interpolation/projection operators.

Parameters

grid	the original grid (Nx() and Ny() must be evenly divisable by pow(2, stages-1)
stages	number of grids in total (The second grid contains half the points of the original grids, The third grid contains half of the second grid ...). Must be >= 1. A good number to start is 3.
ps	parameters necessary for `dg::construct` to construct a `Container` from a `dg::HVec`

Member Function Documentation

◆ construct()

template<class Geometry , class Matrix , class Container >

template<class ... Params>

void dg::MultigridCG2d< Geometry, Matrix, Container >::construct ( Params &&... ps )

inline

Perfect forward parameters to one of the constructors.

Template Parameters

Params deduced by the compiler

Parameters

ps	parameters forwarded to constructors

◆ copyable()

template<class Geometry , class Matrix , class Container >

const Container & dg::MultigridCG2d< Geometry, Matrix, Container >::copyable ( ) const

inline

Return an object of same size as the object used for construction on the finest grid.

Returns: A copyable object; what it contains is undefined, its size is important

◆ grid()

template<class Geometry , class Matrix , class Container >

const Geometry & dg::MultigridCG2d< Geometry, Matrix, Container >::grid ( unsigned stage ) const

inline

return the grid at given stage

Parameters

stage must fulfill 0 <= stage < stages()

◆ max_iter()

template<class Geometry , class Matrix , class Container >

unsigned dg::MultigridCG2d< Geometry, Matrix, Container >::max_iter ( ) const

inline

The maximum number of iterations allowed at stage 0 (if the solution method returns this number, failure is indicated)

◆ num_stages()

template<class Geometry , class Matrix , class Container >

unsigned dg::MultigridCG2d< Geometry, Matrix, Container >::num_stages ( ) const

inline

Returns: number of stages (same as stages)

◆ project() [1/2]

template<class Geometry , class Matrix , class Container >

template<class ContainerType0 >

std::vector< ContainerType0 > dg::MultigridCG2d< Geometry, Matrix, Container >::project ( const ContainerType0 & src ) const

inline

Project vector to all involved grids (allocate memory version)

Parameters

src	the input vector

Returns: the input vector projected to all grids ( index 0 contains a copy of src, 1 is the projetion to the first coarse grid, 2 is the next coarser grid, ...)

◆ project() [2/2]

template<class Geometry , class Matrix , class Container >

template<class ContainerType0 >

void dg::MultigridCG2d< Geometry, Matrix, Container >::project	(	const ContainerType0 &	src,
		std::vector< ContainerType0 > &	out
	)		const

inline

Project vector to all involved grids.

Parameters

src	the input vector (may alias first element of out)
out	the input vector projected to all grids ( index 0 contains a copy of src, 1 is the projetion to the first coarse grid, 2 is the next coarser grid, ...)

Note: out is not resized

◆ set_benchmark()

template<class Geometry , class Matrix , class Container >

void dg::MultigridCG2d< Geometry, Matrix, Container >::set_benchmark	(	bool	benchmark,
		std::string	message = `"Nested Iterations"`
	)

inline

Set or unset performance timings during iterations.

Parameters

benchmark	If true, additional output will be written to `std::cout` during solution
message	An optional identifier that is printed together with the benchmark (intended use is to distinguish different messages in the output)

◆ set_max_iter()

template<class Geometry , class Matrix , class Container >

void dg::MultigridCG2d< Geometry, Matrix, Container >::set_max_iter ( unsigned new_max )

inline

Set the maximum number of iterations allowed at stage 0.

By default this number is the grid size. However, for large simulations you may want to prevent the solver from iterating to that number in case of failure.

Parameters

new_max new maximum number of iterations allowed at stage 0

◆ solve() [1/2]

template<class Geometry , class Matrix , class Container >

template<class MatrixType , class ContainerType0 , class ContainerType1 >

std::vector< unsigned > dg::MultigridCG2d< Geometry, Matrix, Container >::solve	(	std::vector< MatrixType > &	ops,
		ContainerType0 &	x,
		const ContainerType1 &	b,
		std::vector< value_type >	eps
	)

inline

Nested iterations.

Equivalent to the following

Compute residual with given initial guess.
Project residual down to the coarsest grid.
Solve equation on the coarse grid.

interpolate solution up to next finer grid and repeat 3 and 4 until the original grid is reached.

See also: dg::nested_iterations

Note: The weights and preconditioner for the dg::PCG solver is taken from the weights() and precond() method in the MatrixType class

Parameters

ops	Index 0 is the `MatrixType` on the original grid, 1 on the half grid, 2 on the quarter grid, ... `ops`[u].precond() and `ops`[u].weights() need to be callable!
x	(read/write) contains initial guess on input and the solution on output (if the initial guess is good enough the solve may return immediately)
b	The right hand side
eps	the accuracy: iteration stops if \( \|\|b - Ax\|\| < \epsilon( \|\|b\|\| + 1) \). If needed (and it is recommended to tune these values) the accuracy can be set for each stage separately. Per default the same accuracy is used at all stages but \( \epsilon_i = 0.5\epsilon_0\) for i > 0 may be a good value as well.

Returns: the number of iterations in each of the stages beginning with the finest grid

Note: the convergence test on the coarse grids is only evaluated every 10th iteration. This effectively saves one dot product per iteration. The dot product is the main performance bottleneck on the coarse grids.

Template Parameters

MatrixType

Any class for which a specialization of TensorTraits exists and defines a tensor_category derived from AnyMatrixTag. Furthermore, any functor/lambda type with signature void operator()( const ContainerType0&, ContainerType1&) . For example

scalar or Container: Scalars and containers act as diagonal matrices.
std::function<void( const ContainerType0&, ContainerType1&)>
dg::HMatrix and dg::IHMatrix with dg::HVec or std::vector<dg::HVec>
dg::DMatrix and dg::IDMatrix with dg::DVec or std::vector<dg::DVec>
dg::MHMatrix with dg::MHVec or std::vector<dg::MHVec>
dg::MDMatrix with dg::MDVec or std::vector<dg::MDVec>

In case of SelfMadeMatrixTag only those blas2 functions that have a corresponding member function in the Matrix class (e.g. symv( const ContainerType0&, ContainerType1&); ) can be called. If a Container has the RecursiveVectorTag, then the matrix is applied to each of the elements unless the type has the SelfMadeMatrixTag or is a Functor type.

ContainerType

Any class for which a specialization of TensorTraits exists and which fulfills the requirements of the there defined data and execution policies derived from AnyVectorTag and AnyPolicyTag. Among others

dg::HVec (serial), dg::DVec (cuda / omp), dg::MHVec (mpi + serial) or dg::MDVec (mpi + cuda / omp)
std::vector<dg::DVec> (vector of shared device vectors), std::array<double, 4> (array of 4 doubles) or std::map < std::string, dg::DVec> ( a map of named vectors)
double (scalar) and other primitive types ...

If there are several ContainerTypes in the argument list, then TensorTraits must exist for all of them

See also: See The dg dispatch system for a detailed explanation of our type dispatch system

◆ solve() [2/2]

template<class Geometry , class Matrix , class Container >

template<class MatrixType , class ContainerType0 , class ContainerType1 >

std::vector< unsigned > dg::MultigridCG2d< Geometry, Matrix, Container >::solve	(	std::vector< MatrixType > &	ops,
		ContainerType0 &	x,
		const ContainerType1 &	b,
		value_type	eps
	)

inline

Nested iterations.

Equivalent to the following

Compute residual with given initial guess.
Project residual down to the coarsest grid.
Solve equation on the coarse grid.

interpolate solution up to next finer grid and repeat 3 and 4 until the original grid is reached.

See also: dg::nested_iterations

Note: The weights and preconditioner for the dg::PCG solver is taken from the weights() and precond() method in the MatrixType class

Parameters

ops	Index 0 is the `MatrixType` on the original grid, 1 on the half grid, 2 on the quarter grid, ... `ops`[u].precond() and `ops`[u].weights() need to be callable!
x	(read/write) contains initial guess on input and the solution on output (if the initial guess is good enough the solve may return immediately)
b	The right hand side
eps	the accuracy: iteration stops if \( \|\|b - Ax\|\| < \epsilon( \|\|b\|\| + 1) \). If needed (and it is recommended to tune these values) the accuracy can be set for each stage separately. Per default the same accuracy is used at all stages but \( \epsilon_i = 0.5\epsilon_0\) for i > 0 may be a good value as well.

Returns: the number of iterations in each of the stages beginning with the finest grid

Note: the convergence test on the coarse grids is only evaluated every 10th iteration. This effectively saves one dot product per iteration. The dot product is the main performance bottleneck on the coarse grids.

Template Parameters

MatrixType

Any class for which a specialization of TensorTraits exists and defines a tensor_category derived from AnyMatrixTag. Furthermore, any functor/lambda type with signature void operator()( const ContainerType0&, ContainerType1&) . For example

scalar or Container: Scalars and containers act as diagonal matrices.
std::function<void( const ContainerType0&, ContainerType1&)>
dg::HMatrix and dg::IHMatrix with dg::HVec or std::vector<dg::HVec>
dg::DMatrix and dg::IDMatrix with dg::DVec or std::vector<dg::DVec>
dg::MHMatrix with dg::MHVec or std::vector<dg::MHVec>
dg::MDMatrix with dg::MDVec or std::vector<dg::MDVec>

In case of SelfMadeMatrixTag only those blas2 functions that have a corresponding member function in the Matrix class (e.g. symv( const ContainerType0&, ContainerType1&); ) can be called. If a Container has the RecursiveVectorTag, then the matrix is applied to each of the elements unless the type has the SelfMadeMatrixTag or is a Functor type.

ContainerType

Any class for which a specialization of TensorTraits exists and which fulfills the requirements of the there defined data and execution policies derived from AnyVectorTag and AnyPolicyTag. Among others

dg::HVec (serial), dg::DVec (cuda / omp), dg::MHVec (mpi + serial) or dg::MDVec (mpi + cuda / omp)
std::vector<dg::DVec> (vector of shared device vectors), std::array<double, 4> (array of 4 doubles) or std::map < std::string, dg::DVec> ( a map of named vectors)
double (scalar) and other primitive types ...

If there are several ContainerTypes in the argument list, then TensorTraits must exist for all of them

See also: See The dg dispatch system for a detailed explanation of our type dispatch system

◆ stages()

template<class Geometry , class Matrix , class Container >

unsigned dg::MultigridCG2d< Geometry, Matrix, Container >::stages ( ) const

inline

Returns: number of stages (same as num_stages)

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

multigrid.h

Public Types

Public Member Functions

Detailed Description

Member Typedef Documentation

◆ container_type

◆ geometry_type

◆ matrix_type

◆ value_type

Constructor & Destructor Documentation

◆ MultigridCG2d() [1/2]

◆ MultigridCG2d() [2/2]

Member Function Documentation

◆ construct()

◆ copyable()

◆ grid()

◆ max_iter()

◆ num_stages()

◆ project() [1/2]

◆ project() [2/2]

◆ set_benchmark()

◆ set_max_iter()

◆ solve() [1/2]

◆ solve() [2/2]

◆ stages()