27template<
class Geometry,
class Matrix,
class Container>
35 NestedGrids(): m_stages(0), m_grids(0), m_inter(0), m_project(0){}
44 template<
class ...ContainerParams>
52 m_grids[0].reset(
grid);
55 for(
unsigned u=1; u<
stages; u++)
57 m_grids[u] = m_grids[u-1];
58 m_grids[u]->multiplyCellNumbers(0.5, 0.5);
63 m_project.resize(
stages-1);
64 for(
unsigned u=0; u<
stages-1; u++)
69 2, 2), std::forward<ContainerParams>(ps)...);
71 1, 2, 2), std::forward<ContainerParams>(ps)...);
73 for(
unsigned u=0; u<m_stages; u++)
75 *m_grids[u]), std::forward<ContainerParams>(ps)...);
76 m_w = m_r = m_b = m_x;
81 template<
class ...Params>
93 template<
class ContainerType0>
94 void project(
const ContainerType0& src, std::vector<ContainerType0>& out)
const
97 for(
unsigned u=0; u<m_stages-1; u++)
106 template<
class ContainerType0>
107 std::vector<ContainerType0>
project(
const ContainerType0& src)
const
110 std::vector<Container> out( m_x);
126 const Geometry&
grid(
unsigned stage)
const {
127 return *(m_grids[stage]);
134 return m_inter[stage];
140 return m_project[stage];
142 Container&
x(
unsigned stage){
return m_x[stage];}
143 const Container&
x(
unsigned stage)
const{
return m_x[stage];}
144 Container&
r(
unsigned stage){
return m_r[stage];}
145 const Container&
r(
unsigned stage)
const{
return m_r[stage];}
146 Container&
b(
unsigned stage){
return m_b[stage];}
147 const Container&
b(
unsigned stage)
const{
return m_b[stage];}
148 Container&
w(
unsigned stage){
return m_w[stage];}
149 const Container&
w(
unsigned stage)
const{
return m_w[stage];}
153 std::vector< dg::ClonePtr< Geometry> > m_grids;
154 std::vector< MultiMatrix<Matrix, Container> > m_inter;
155 std::vector< MultiMatrix<Matrix, Container> > m_project;
156 std::vector< Container> m_x, m_r, m_b, m_w;
196template<
class MatrixType0,
class ContainerType0,
class ContainerType1,
class MatrixType1,
class NestedGr
ids>
198 std::vector<MatrixType0>& ops, ContainerType0& x,
const ContainerType1& b,
199 std::vector<MatrixType1>& inverse_ops,
NestedGrids& nested_grids)
207 for(
unsigned u=0; u<nested.
stages()-1; u++)
218 for(
unsigned u=nested.
stages()-1; u>0; u--)
221 dg::apply( inverse_ops[u], nested.
b(u), nested.
x(u));
275template<
class NestedGr
ids,
class MatrixType0,
class MatrixType1,
class MatrixType2>
277 std::vector<MatrixType0>& ops,
278 std::vector<MatrixType1>& inverse_ops_down,
279 std::vector<MatrixType2>& inverse_ops_up,
280 NestedGrids& nested_grids,
unsigned gamma,
unsigned p)
296 dg::apply( inverse_ops_down[p], nested.
b(p), nested.
x(p));
311 if( p+1 == nested.
stages()-1)
314 dg::apply( inverse_ops_up[p+1], nested.
b(p+1), nested.
x(p+1));
323 for(
unsigned u=0; u<gamma; u++)
335 dg::apply(inverse_ops_up[p], nested.
b(p), nested.
x(p));
364template<
class MatrixType0,
class MatrixType1,
class MatrixType2,
class NestedGr
ids,
class ContainerType0,
class ContainerType1>
366 std::vector<MatrixType0>& ops, ContainerType0& x,
const ContainerType1& b,
367 std::vector<MatrixType1>& inverse_ops_down,
368 std::vector<MatrixType2>& inverse_ops_up,
369 NestedGrids& nested_grids,
unsigned gamma,
unsigned mu)
378 for(
unsigned u=0; u<nested.
stages()-1; u++)
389 unsigned s = nested.
stages()-1;
391 dg::apply( inverse_ops_up[s], nested.
b(s), nested.
x(s));
400 for(
int p=nested.
stages()-2; p>=1; p--)
402 for(
unsigned u=0; u<mu; u++)
403 multigrid_cycle( ops, inverse_ops_down, inverse_ops_up, nested, gamma, p);
409 for(
unsigned u=0; u<mu; u++)
410 multigrid_cycle( ops, inverse_ops_down, inverse_ops_up, nested, gamma, 0);
438template<
class NestedGrids,
class MatrixType0,
class MatrixType1,
class MatrixType2,
439 class ContainerType0,
class ContainerType1,
class ContainerType2>
441 std::vector<MatrixType0>& ops,
442 ContainerType0& x,
const ContainerType1& b,
443 std::vector<MatrixType1>& inverse_ops_down,
444 std::vector<MatrixType2>& inverse_ops_up,
446 const ContainerType2& weights,
double eps,
unsigned gamma)
453 full_multigrid( ops,
x, b, inverse_ops_down, inverse_ops_up, nested_grids, gamma, 1);
461 double error = sqrt(
blas2::dot(weights,nested_grids.
r(0)) );
463 while ( error > eps*(nrmb + 1))
469 full_multigrid( ops,
x, b, inverse_ops_down, inverse_ops_up, nested_grids, gamma, 1);
477 error = sqrt(
blas2::dot(weights,nested_grids.
r(0)) );
499template<
class Geometry,
class Matrix,
class Container>
516 template<
class ...ContainerParams>
518 ContainerParams&& ... ps):
522 for (
unsigned u = 0; u <
stages; u++)
523 m_pcg[u].construct(m_nested.x(u), m_nested.grid(u).size());
532 template<
class ...Params>
540 template<
class ContainerType0>
541 void project(
const ContainerType0& src, std::vector<ContainerType0>& out)
const
543 m_nested.project( src, out);
547 template<
class ContainerType0>
548 std::vector<ContainerType0>
project(
const ContainerType0& src)
const
550 return m_nested.project( src);
553 unsigned stages()
const{
return m_nested.stages();}
559 const Geometry&
grid(
unsigned stage)
const {
560 return m_nested.grid(stage);
566 unsigned max_iter()
const{
return m_pcg[0].get_max();}
583 void set_benchmark(
bool benchmark, std::string message =
"Nested Iterations"){
584 m_benchmark = benchmark;
590 const Container&
copyable()
const {
return m_nested.copyable();}
617 template<
class MatrixType,
class ContainerType0,
class ContainerType1>
618 std::vector<unsigned>
solve( std::vector<MatrixType>& ops, ContainerType0& x,
const ContainerType1& b,
value_type eps)
620 std::vector<value_type> v_eps( m_stages, eps);
621 for(
unsigned u=m_stages-1; u>0; u--)
623 return solve( ops,
x, b, v_eps);
626 template<
class MatrixType,
class ContainerType0,
class ContainerType1>
627 std::vector<unsigned>
solve( std::vector<MatrixType>& ops, ContainerType0& x,
const ContainerType1& b, std::vector<value_type> eps)
631 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
633 std::vector<unsigned> number(m_stages);
634 std::vector<std::function<void(
const ContainerType1&, ContainerType0&)> >
635 multi_inv_pol(m_stages);
636 for(
unsigned u=0; u<m_stages; u++)
638 multi_inv_pol[u] = [&, u, &pcg = m_pcg[u], &pol = ops[u]](
639 const auto&
y,
auto&
x)
644 number[u] = pcg.solve( pol,
x,
y, pol.precond(),
645 pol.weights(), eps[u], 1, 1);
647 number[u] = pcg.solve( pol,
x,
y, pol.precond(),
648 pol.weights(), eps[u], 1, 10);
651 DG_RANK0 std::cout <<
"# `"<<m_message<<
"` stage: " << u <<
", iter: " << number[u] <<
", took "<<t.
diff()<<
"s\n";
661 std::vector< PCG<Container> > m_pcg;
663 bool m_benchmark =
true;
664 std::string m_message =
"Nested Iterations";
class intended for the use in throw statements
Definition exceptions.h:83
void append_line(const Message &message)
Appends a newline and a message verbatim to the what string.
Definition exceptions.h:108
small class holding a stringstream
Definition exceptions.h:29
Error classes or the dg library.
#define _ping_
Definition exceptions.h:12
A matrix type for fast interpolations/projections.
DG_DEVICE T zero(T x, Ts ...xs)
This enum can be used in dg::evaluate.
Definition functions.h:19
void copy(const ContainerTypeIn &source, ContainerTypeOut &target)
Definition blas1.h:243
void axpby(value_type alpha, const ContainerType1 &x, value_type1 beta, ContainerType &y)
Definition blas1.h:306
ContainerType construct(const from_ContainerType &from, Params &&... ps)
Generic way to construct an object of ContainerType given a from_ContainerType object and optional ad...
Definition blas1.h:792
void gemv(get_value_type< ContainerType1 > alpha, MatrixType &&M, const ContainerType1 &x, get_value_type< ContainerType1 > beta, ContainerType2 &y)
Alias for blas2::symv ;.
Definition blas2.h:339
auto dot(const ContainerType1 &x, const MatrixType &m, const ContainerType2 &y)
; Binary reproducible general dot product
Definition blas2.h:94
void apply(get_value_type< ContainerType1 > alpha, MatrixType &&M, const ContainerType1 &x, get_value_type< ContainerType1 > beta, ContainerType2 &y)
Alias for dg::blas2::symv ;.
Definition blas2.h:490
void symv(MatrixType &&M, const ContainerType1 &x, ContainerType2 &y)
Definition blas2.h:325
@ forward
forward derivative (cell to the right and current cell)
Definition enums.h:98
typename TensorTraits< std::decay_t< Vector > >::value_type get_value_type
Definition tensor_traits.h:45
auto evaluate(Functor &&f, const Topology &g)
Evaluate a function on grid coordinates
Definition evaluation.h:74
EllSparseBlockMat< real_type, thrust::host_vector > fast_interpolation(unsigned coord, const aRealTopology< real_type, Nd > &t, unsigned multiplyn, unsigned multiplyNx)
Create interpolation matrix for integer multipliers.
Definition fast_interpolation.h:302
EllSparseBlockMat< real_type, thrust::host_vector > fast_projection(unsigned coord, const aRealTopology< real_type, Nd > &t, unsigned dividen, unsigned divideNx)
Create projecton matrix for integer dividers.
Definition fast_interpolation.h:314
void multigrid_cycle(std::vector< MatrixType0 > &ops, std::vector< MatrixType1 > &inverse_ops_down, std::vector< MatrixType2 > &inverse_ops_up, NestedGrids &nested_grids, unsigned gamma, unsigned p)
EXPERIMENTAL Full approximation multigrid cycle.
Definition multigrid.h:276
void fmg_solve(std::vector< MatrixType0 > &ops, ContainerType0 &x, const ContainerType1 &b, std::vector< MatrixType1 > &inverse_ops_down, std::vector< MatrixType2 > &inverse_ops_up, NestedGrids &nested_grids, const ContainerType2 &weights, double eps, unsigned gamma)
EXPERIMENTAL Full multigrid cycles.
Definition multigrid.h:440
void full_multigrid(std::vector< MatrixType0 > &ops, ContainerType0 &x, const ContainerType1 &b, std::vector< MatrixType1 > &inverse_ops_down, std::vector< MatrixType2 > &inverse_ops_up, NestedGrids &nested_grids, unsigned gamma, unsigned mu)
EXPERIMENTAL One Full multigrid cycle.
Definition multigrid.h:365
void nested_iterations(std::vector< MatrixType0 > &ops, ContainerType0 &x, const ContainerType1 &b, std::vector< MatrixType1 > &inverse_ops, NestedGrids &nested_grids)
Full approximation nested iterations.
Definition multigrid.h:197
Interpolation matrix creation functions.
Useful MPI typedefs and overloads of interpolation and projection.
This is the namespace for all functions and classes defined and used by the discontinuous Galerkin li...
Struct that applies given matrices one after the other.
Definition fast_interpolation.h:37
Solve.
Definition multigrid.h:501
std::vector< unsigned > solve(std::vector< MatrixType > &ops, ContainerType0 &x, const ContainerType1 &b, std::vector< value_type > eps)
Nested iterations.
Definition multigrid.h:627
get_value_type< Container > value_type
Definition multigrid.h:505
std::vector< unsigned > solve(std::vector< MatrixType > &ops, ContainerType0 &x, const ContainerType1 &b, value_type eps)
Nested iterations.
Definition multigrid.h:618
unsigned stages() const
Definition multigrid.h:553
const Container & copyable() const
Return an object of same size as the object used for construction on the finest grid.
Definition multigrid.h:590
void project(const ContainerType0 &src, std::vector< ContainerType0 > &out) const
Project vector to all involved grids.
Definition multigrid.h:541
void construct(Params &&...ps)
Perfect forward parameters to one of the constructors.
Definition multigrid.h:533
MultigridCG2d(const Geometry &grid, const unsigned stages, ContainerParams &&... ps)
Construct the grids and the interpolation/projection operators.
Definition multigrid.h:517
const Geometry & grid(unsigned stage) const
return the grid at given stage
Definition multigrid.h:559
Geometry geometry_type
Definition multigrid.h:502
void set_max_iter(unsigned new_max)
Set the maximum number of iterations allowed at stage 0.
Definition multigrid.h:575
MultigridCG2d()=default
Allocate nothing, Call construct method before usage.
unsigned max_iter() const
Definition multigrid.h:566
void set_benchmark(bool benchmark, std::string message="Nested Iterations")
Set or unset performance timings during iterations.
Definition multigrid.h:583
std::vector< ContainerType0 > project(const ContainerType0 &src) const
Project vector to all involved grids (allocate memory version)
Definition multigrid.h:548
unsigned num_stages() const
Definition multigrid.h:555
Container container_type
Definition multigrid.h:504
Matrix matrix_type
Definition multigrid.h:503
Hold nested grids and provide dg fast interpolation and projection matrices.
Definition multigrid.h:29
const Container & w(unsigned stage) const
Definition multigrid.h:149
const Container & copyable() const
Return an object of same size as the object used for construction on the finest grid.
Definition multigrid.h:117
Container & b(unsigned stage)
Definition multigrid.h:146
const Container & b(unsigned stage) const
Definition multigrid.h:147
const MultiMatrix< Matrix, Container > & interpolation(unsigned stage) const
return the interpolation matrix at given stage
Definition multigrid.h:132
Container container_type
Definition multigrid.h:32
void project(const ContainerType0 &src, std::vector< ContainerType0 > &out) const
Project vector to all involved grids.
Definition multigrid.h:94
const MultiMatrix< Matrix, Container > & projection(unsigned stage) const
return the projection matrix at given stage
Definition multigrid.h:138
void construct(Params &&...ps)
Perfect forward parameters to one of the constructors.
Definition multigrid.h:82
std::vector< ContainerType0 > project(const ContainerType0 &src) const
Project vector to all involved grids (allocate memory version)
Definition multigrid.h:107
NestedGrids()
Allocate nothing, Call construct method before usage.
Definition multigrid.h:35
Container & r(unsigned stage)
Definition multigrid.h:144
const Container & x(unsigned stage) const
Definition multigrid.h:143
unsigned stages() const
Definition multigrid.h:120
Geometry geometry_type
Definition multigrid.h:30
const Container & r(unsigned stage) const
Definition multigrid.h:145
Container & w(unsigned stage)
Definition multigrid.h:148
const Geometry & grid(unsigned stage) const
return the grid at given stage
Definition multigrid.h:126
unsigned num_stages() const
Definition multigrid.h:122
NestedGrids(const Geometry &grid, const unsigned stages, ContainerParams &&...ps)
Construct the grids and the interpolation/projection operators.
Definition multigrid.h:45
Container & x(unsigned stage)
Definition multigrid.h:142
get_value_type< Container > value_type
Definition multigrid.h:33
Matrix matrix_type
Definition multigrid.h:31
Simple tool for performance measuring.
Definition dg_doc.h:352
double diff() const
Return time in seconds elapsed between tic and toc.
#define DG_RANK0
Definition typedefs.h:146