VCycle

allscale::api::user::algorithm::VCycle

Defined in header "allscale/api/user/algorithm/vcycle.h"

template<template<typename M,unsigned L> class StageBody, typename Mesh>
class VCycle;

A generic V-cycle implementation enabling the creation of a V-cycle solver by providing an implementation of a V-cycle stage body. For an implementation using the V-cycle operator take a look at the Meshes tutorial.

Member types

Member type	Definition
mesh_type	Mesh

Member functions

Member function	Description
VCycle(const mesh_type& mesh)	Construct a VCycle
vcycle_reference run(std::size_t numCycles = 1)	run numCycles iterations of the V-cycle
template<unsigned Level = 0> const StageBody<Mesh,Level>& getStageBody() const	Read access to the stage body on Level
template<unsigned Level = 0> StageBody<Mesh,Level>& getStageBody()	Mutable access to the stage body on Level

StageBody

The VCycle class takes the StageBody as template argument.

StageBody<Mesh,Level>

The StageBody has to implement the following methods:

Member function	Description
StageBody(const Mesh& mesh)	Construct a StageBody with mesh
void computeFineToCoarse()	compute stage while going up the hierarchy
void computeCoarseToFine()	compute stage while going down the hierarchy
void restrictFrom(StageBody<Mesh,Level-1>& childStage)	go from child stage to this stage (going up the hierarchy)
void prolongateTo(StageBody<Mesh,Level-1>& childStage)	go from this stage to child stage (going down the hierarchy)

Examples

The code below uses the TestStage as StageBody. The constructor of the struct initializes the updateCounters by settings it to 0 for all Cells on that layer.

The VCycle starts by calling computeFineToCoarse on layer zero which will increase the counter for each Cell on that layer by one. The restrictFrom method will be called afterwards and just copy the value of the Cells of layer zero to the Cells of layer one. This will be repeated until computeFineToCoarse is called for the top layer. In the example code the Cells of the top layer will have a value of three.

When going down the hierarchy of the mesh restrictFrom copies the value of the higher layer to the lower layer. Afterwards the computeCoarseToFine increases the value of each Cell by one. This is repeated until the bottom layer is reached. In the case of the example code the value of each Cell at the bottom layer is five. The VCycle gets executed 10 times which will result in the resulting value of 50 for each Cell at the bottom layer.

#include "allscale/api/user/data/mesh.h"
#include "allscale/api/user/algorithm/vcycle.h"

using namespace allscale::api::user;

// define 'object' types
struct Cell {};

// define 'relations'
struct Cell2Cell : public data::edge<Cell,Cell> {};

// and 'hierarchies'
struct Cell2Child : public data::hierarchy<Cell,Cell> {};

template<
    typename Mesh,
    unsigned Level
>
struct TestStage {

    template<typename NodeKind,typename ValueType>
    using attribute = typename Mesh::template mesh_data_type<NodeKind,ValueType,Level>;

    // capture mesh
    const Mesh& mesh;

    // counts the number of updates per node
    attribute<Cell,int> updateCounters;


    TestStage(const Mesh& mesh)
        : mesh(mesh),
            updateCounters(mesh.template createNodeData<Cell,int,Level>()) {

        mesh.template pforAll<Cell,Level>([&](const auto& cur){
            updateCounters[cur] = 0;
        });

    }

    void computeFineToCoarse() {
        mesh.template pforAll<Cell,Level>([&](const auto& cur) {
            updateCounters[cur]++;
        });
    }

    void computeCoarseToFine() {
        mesh.template pforAll<Cell,Level>([&](const auto& cur) {
            updateCounters[cur]++;
        });
    }


    void restrictFrom(const TestStage<Mesh,Level-1>& childStage) {

        // reduction by taking average of children
        mesh.template pforAll<Cell,Level>([&](auto cur){
            // check that all children have the same number of updates
            const auto& children = mesh.template getChildren<Cell2Child>(cur);
            assert_false(children.empty());

            int numUpdates = childStage.updateCounters[children.front()];
            for(const auto& child : children) {
                assert_eq(numUpdates,childStage.updateCounters[child]);
            }

            assert_gt(numUpdates,updateCounters[cur]);
            updateCounters[cur] = numUpdates;
        });

    }

    void prolongateTo(TestStage<Mesh,Level-1>& childStage) {

        // prolongation of results to children
        mesh.template pforAll<Cell,Level>([&](auto cur){
            // check that all children have the same number of updates
            const auto& children = mesh.template getChildren<Cell2Child>(cur);
            assert_false(children.empty());

            int numUpdates = childStage.updateCounters[children.front()];
            for(const auto& child : children) {
                assert_eq(numUpdates,childStage.updateCounters[child]);
            }

            assert_gt(updateCounters[cur],numUpdates);

            // propagate update counters to child stages
            for(const auto& child : children) {
                childStage.updateCounters[child] = updateCounters[cur];
            }

        });
    }

};

data::MeshBuilder<
    data::nodes<Cell>,
    data::edges<Cell2Cell>,
    data::hierarchies<Cell2Child>,
    3            // 3 layers
> mesh_builder;

// create four cells on layer 0
auto l0a = mesh_builder.create<Cell>();
auto l0b = mesh_builder.create<Cell>();
auto l0c = mesh_builder.create<Cell>();
auto l0d = mesh_builder.create<Cell>();

// create two cells on layer 1
auto l1a = mesh_builder.create<Cell,1>();
auto l1b = mesh_builder.create<Cell,1>();

// create a final cell on layer 2
auto l2a = mesh_builder.create<Cell,2>();

// link first layer
mesh_builder.link<Cell2Cell>(l0a,l0b);
mesh_builder.link<Cell2Cell>(l0b,l0a);
mesh_builder.link<Cell2Cell>(l0b,l0c);
mesh_builder.link<Cell2Cell>(l0c,l0b);
mesh_builder.link<Cell2Cell>(l0c,l0d);
mesh_builder.link<Cell2Cell>(l0d,l0c);

// link second layer
mesh_builder.link<Cell2Cell>(l1a,l1b);
mesh_builder.link<Cell2Cell>(l1b,l1a);

// link hierarchies
mesh_builder.link<Cell2Child>(l1a,l0a);
mesh_builder.link<Cell2Child>(l1a,l0b);
mesh_builder.link<Cell2Child>(l1b,l0c);
mesh_builder.link<Cell2Child>(l1b,l0d);

mesh_builder.link<Cell2Child>(l2a,l1a);
mesh_builder.link<Cell2Child>(l2a,l1b);

// build mesh
auto mesh = mesh_builder.build();

using vcycle_type = algorithm::VCycle<
    TestStage,
    decltype(mesh)
>;

// create vcycle instance
vcycle_type vcycle(mesh);
auto& counts = vcycle.getStageBody().updateCounters;

// counters should all be initially 0
for(std::size_t i = 0; i < mesh.getNumNodes<Cell>(); ++i) {
    assert_eq(0, counts[data::NodeRef<Cell>((data::node_index_t)i)]);
}
vcycle.run(10);

// now each element should be updated 50x
for(std::size_t i = 0; i < mesh.getNumNodes<Cell>(); ++i) {
    assert_eq(50, counts[data::NodeRef<Cell>((data::node_index_t)i)]);
}