Commit 707a73ff authored by Müller, Felix's avatar Müller, Felix
Browse files

Redesign of automatic basis, communicator and DOFVector updates

Major changes:
  Add wrapper class AdaptiveGrid<Grid> to be used instead of a
    dune-grid within AMDiS. This provides most of the grid interface
    and supports automatic updates of bases.
  Change: MeshCreator now returns an AdaptiveGrid
  Change: ProblemStat uses AdaptiveGrid internally and returns an
    object of that type
  Add class ParallelGlobalBasis<Basis> as a replacement for a
    dune-functions basis. This implements the complete interface and
    can be used as a replacement for a regular basis in all contexts.
    It supports automatic updates.
  Add free function makeGlobalBasis to make a ParallelGlobalBasis
  Change: AMDiS basis creators now return a ParallelGlobalBasis
  Change: ProblemStat uses ParallelGlobalBasis internally and returns
    an object of that type
  Add classes Observer/Signals. Classes can derive from those to take
    part in the automatic update feature.
  Add update(Event const& e) methods to classes using the feature
    above. They implement the action to be done when the registered
    Event is triggered.
  Remove GridTransfer[-Manager] and replace with the new
    Observer/Signals implementation
  Change: linear algebra communicator is now a part of
  Add wrapper for DataTransfer objects that can be used instead of a
    base class pointer and performs deep-copy operations on copy
    construction and assignment
  Change: DOFVector DataTransfer member is now stored using the new
    wrapper class DataTransferWrapper
  Change: preAdapt, adapt, postAdapt callbacks are now called after
    calling the respective grid function. Classes that previously used
    any of them had their methods changed accordingly.

Minor changes:
  Remove Comm argument for ctor of [Bi-]LinearForm, DOFVector
  Change ctor argument from Comm to Basis in
    linearalgebra/<impl>/[Matrix-, Vector-]Backend.hpp
  Simplify DOFVector copy/move by defaulting
  Remove DOFVector::resize override
  Remove DOFVectorInterface.hpp
  Change: linearalgebra/Traits.hpp includes the backend-specific
    Traits file. This is required to break an include loop.
  Add global 'solver category' initfile parameter used by the ISTL
    backend communicator when no parameter  at prefix + '->category'
    is found
  Fix commented out parts in DOFVectorTest, DiscreteFunctionTest
  Add unit test for Observer functionality
  Add helper function wrap_or_share (see amdis/common/SharedPtr.hpp)
  Fix backend = ISTL check for ISTLCommTest in CMakeLists
parent af7a257f
...@@ -65,7 +65,7 @@ void run(Grid& grid) ...@@ -65,7 +65,7 @@ void run(Grid& grid)
using BC = PeriodicBC<FieldVector<double,2>, typename Traits::GlobalBasis::MultiIndex>; using BC = PeriodicBC<FieldVector<double,2>, typename Traits::GlobalBasis::MultiIndex>;
BC periodicBC(prob.boundaryManagerPtr(),-1,{{{1.0,0.0}, {0.0,1.0}}, {1.0, 0.0}}); BC periodicBC(prob.boundaryManagerPtr(),-1,{{{1.0,0.0}, {0.0,1.0}}, {1.0, 0.0}});
periodicBC.init(prob.globalBasis(), prob.globalBasis()); periodicBC.init(*prob.globalBasis(), *prob.globalBasis());
std::cout << "periodicNodes:\n"; std::cout << "periodicNodes:\n";
for (std::size_t i = 0; i < periodicBC.periodicNodes().size(); ++i) for (std::size_t i = 0; i < periodicBC.periodicNodes().size(); ++i)
std::cout << i << ": " << periodicBC.periodicNodes()[i] << "\n"; std::cout << i << ": " << periodicBC.periodicNodes()[i] << "\n";
#pragma once
#include <algorithm>
#include <list>
#include <memory>
#include <mutex>
#include <shared_mutex>
#include <type_traits>
#include <utility>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/geometry/type.hh>
#include <amdis/common/ConceptsBase.hpp>
#include <amdis/common/SharedPtr.hpp>
#include <amdis/Observer.hpp>
#include <amdis/Output.hpp>
namespace AMDiS
namespace Impl
template <class Grid, class... Args,
REQUIRES(std::is_convertible<decltype(std::declval<Grid>().loadBalance()), bool>::value)>
bool loadBalanceImpl(Grid& grid, Args&&... args)
return grid.loadBalance(FWD(args)...);
// Workaround for grids not adhering to the correct interface
template <class Grid, class... Args,
REQUIRES(not std::is_convertible<decltype(std::declval<Grid>().loadBalance()), bool>::value)>
bool loadBalanceImpl(Grid& grid, Args&&... args)
return true;
namespace event
/** Event generated from an AdaptiveGrid when calling preAdapt(). It contains the return value
* of preAdapt() as its 'mightCoarsen' member and is passed to registered observers after
* calling preAdapt on the underlying grid.
struct preAdapt { bool mightCoarsen; };
/** Event generated from an AdaptiveGrid when calling adapt(). Its 'adapted' member contains the
* value true if either preAdapt() or adapt() returned true. This event is passed to registered
* observers after calling adapt on the underlying grid.
struct adapt { bool adapted; };
/** Event generated from an AdaptiveGrid when calling postAdapt().This event is passed to
* registered observers after calling postAdapt on the underlying grid.
struct postAdapt {};
/** Wrapper class for Dune-grids that allows automatic signalling of events during grid adaptation
* This class needs to be created after construction of the associated grid by calling the
* instance method.
* Calls to grid.preAdapt(), grid.adapt() and grid.postAdapt() need to be replaced by calls to
* the corresponding methods of this class to use the automatic update functionality.
template <class Grid>
class AdaptiveGrid
: public Signals<event::preAdapt, event::adapt, event::postAdapt>
using Self = AdaptiveGrid<Grid>;
using Element = typename Grid::template Codim<0>::Entity;
struct HiddenStruct {};
using HostGrid = Grid;
enum { dimension = HostGrid::dimension };
enum { dimensionworld = HostGrid::dimensionworld };
using LeafGridView = typename HostGrid::LeafGridView;
using LevelGridView = typename HostGrid::LevelGridView;
template <int cd>
struct Codim
using Geometry = typename HostGrid::template Codim<cd>::Geometry;
using Entity = typename HostGrid::template Codim<cd>::Entity;
using GlobalIdSet = typename HostGrid::GlobalIdSet;
using LocalIdSet = typename HostGrid::LocalIdSet;
using LevelIndexSet = typename HostGrid::LevelIndexSet;
using LeafIndexSet = typename HostGrid::LeafIndexSet;
/// Constructor that may only be called indirectly via the instance method
explicit AdaptiveGrid(std::shared_ptr<HostGrid> grid, HiddenStruct)
: hostGrid_(std::move(grid))
/// Unreachable constructor required to compile the unreachable branch in instance method
explicit AdaptiveGrid(std::shared_ptr<HostGrid const> grid, HiddenStruct)
error_exit("Cannot create AdaptiveGrid from const Grid.");
int maxLevel() const { return hostGrid_->maxLevel(); }
int size(int level, int codim) const { return hostGrid_->size(level, codim); }
/// Return number of leaf entities of a given codim in this process
int size(int codim) const { return hostGrid_->size(codim); }
/// Return number of entities per level and geometry type in this process
int size(int level, Dune::GeometryType type) const { return hostGrid_->size(level, type); }
/// Return number of leaf entities per geometry type in this process
int size(Dune::GeometryType type) const { return hostGrid_->size(type); }
auto numBoundarySegments () const { return hostGrid_->numBoundarySegments(); }
/// View for a grid level for All_Partition
LevelGridView levelGridView(int l) const { return hostGrid_->levelGridView(l); }
/// View for the leaf grid for All_Partition
LeafGridView leafGridView() const { return hostGrid_->leafGridView(); }
GlobalIdSet const& globalIdSet() const { return hostGrid_->globalIdSet(); }
/// return const reference to the host grids local id set
LocalIdSet const& localIdSet() const { return hostGrid_->localIdSet(); }
/// return const reference to the host grids level index set for level level
LevelIndexSet const& levelIndexSet(int level) const { return hostGrid_->levelIndexSet(level); }
/// return const reference to the host grids leaf index set
LeafIndexSet const& leafIndexSet() const { return hostGrid_->leafIndexSet(); }
/// Refines all grid elements refCount times.
void globalRefine(int refCount)
for (int i = 0; i < refCount; ++i) {
// mark all entities for grid refinement
for (const auto& element : elements(hostGrid_->leafGridView()))
hostGrid_->mark(1, element);
/// Mark entity for refinement
bool mark(int refCount, Element const& e) { return hostGrid_->mark(refCount, e); }
/// Return refinement mark for entity
int getMark(Element const& e) const { return hostGrid_->getMark(e); }
/// Prepare the grid for adaptation and notify observers of the preAdapt event
bool preAdapt()
Dune::Timer t;
mightCoarsen_ = hostGrid_->preAdapt();
Dune::MPIHelper::getCollectiveCommunication().max(&mightCoarsen_, 1);
info(2,"preAdapt needed {} seconds", t.elapsed());
return mightCoarsen_;
/// Adapt the grid and notify observers of the adapt event
bool adapt()
Dune::Timer t;
refined_ = hostGrid_->adapt();
Dune::MPIHelper::getCollectiveCommunication().max(&refined_, 1);
this->notify(event::adapt{mightCoarsen_ || refined_});
return refined_;
/// Perform cleanup after grid adaptation and notify observers of the postAdapt event
void postAdapt()
Dune::Timer t;
info(2,"postAdapt needed {} seconds", t.elapsed());
/** Calls loadBalance on the underlying grid.
* Re-balances the load each process has to handle for a parallel grid.
* \return True if the grid has changed.
bool loadBalance()
return Impl::loadBalanceImpl(*hostGrid_);
/* Calls loadBalance(handle) on the underlying grid.
* Re-balances the load each process has to handle for a parallel grid and moves the data.
* \param data A data handle telling the method which data should be communicated
* and how. Has to adhere to the interface described by Dune::CommDataHandleIf.
* \return True if the grid has changed.
template <class DataHandle>
bool loadBalance (DataHandle& handle)
return Impl::loadBalanceImpl(*hostGrid_, handle);
/// Returns the grid change index, see \ref changeIndex.
unsigned long changeIndex() const
return changeIndex_;
template <class Grid_>
static std::shared_ptr<Self>
instanceImpl(std::shared_ptr<Grid_> grid)
using mutex_type = std::shared_timed_mutex;
static mutex_type access_mutex;
// 1. Cleanup. Since the iterators may be invalidated only one accessor is allowed. Erases all
// expired weak pointers (all pointers that no longer have valid instances attached).
std::unique_lock<mutex_type> write_lock(access_mutex);
for (auto it = adaptiveGrids_.begin(); it != adaptiveGrids_.end(); ++it)
if (it->expired())
it = adaptiveGrids_.erase(it);
std::shared_lock<mutex_type> read_lock(access_mutex);
// 2. Find matching observed grid object. We obtain a lock here to avoid complications with
// insertions made by other threads in step 3.
auto it = find_if(adaptiveGrids_.begin(), adaptiveGrids_.end(),
[&](std::weak_ptr<Self>& wPtr) {
auto ag = wPtr.lock();
return ag->hostGrid_ == grid;
// 3. Register new object or return existing. In case of inserting a new instance we obtain a
// write lock.
if (it == adaptiveGrids_.end())
"No existing AdaptiveGrid found and no mutable grid passed to create a new instance");
auto ptr = std::make_shared<Self>(std::move(grid), HiddenStruct{});
return std::move(ptr);
return it->lock();
// Do-nothing specialization if argument is already an AdaptiveGrid
static std::shared_ptr<Self> instanceImpl(std::shared_ptr<Self> grid) { return grid; }
/** Returns the AdaptiveGrid associated to the grid passed.
* If no AdaptiveGrid exists, this method creates a new one if the passed grid is mutable,
* otherwise the call fails with an error.
template <class Grid_>
static std::shared_ptr<Self> instance(Grid_&& grid)
return instanceImpl(wrap_or_share(FWD(grid)));
// Test for equality of the grid pointers
bool operator==(AdaptiveGrid<HostGrid> const& that) const
return (hostGrid_ == that.hostGrid_);
/// Return the underlying grid
std::shared_ptr<HostGrid> hostGrid() const { return hostGrid_; }
/// List of previously created instances, indexed by address of the HostGrid.
* For each grid type Grid we maintain a list of instances handed out by the instance method.
* We use weak pointers here that are valid as long as there is at least one other place where
* the shared pointer to the instance is used. When the instance is no longer used the weak
* pointers here do not hinder the object's destruction.
static std::list<std::weak_ptr<Self>> adaptiveGrids_;
/// The underlying grid implementation
std::shared_ptr<HostGrid> hostGrid_;
/// Flag set during \ref preAdapt(), indicating whether any element might be coarsened in \ref adapt()
bool mightCoarsen_ = false;
/// Flag set during \ref adapt() indicating that at least one entity was refined
bool refined_ = false;
/// This index is incremented every time the grid is changed, e.g. by refinement or coarsening.
unsigned long changeIndex_ = 0;
template <class Grid>
std::list<std::weak_ptr<AdaptiveGrid<Grid>>> AdaptiveGrid<Grid>::adaptiveGrids_;
} // end namespace AMDiS
...@@ -45,12 +45,11 @@ namespace AMDiS ...@@ -45,12 +45,11 @@ namespace AMDiS
using ElementMatrix = FlatMatrix<CoefficientType>; using ElementMatrix = FlatMatrix<CoefficientType>;
public: public:
/// Constructor. Wraps the reference into a non-destroying shared_ptr or moves /// Constructor. Wraps the reference into a non-destroying shared_ptr or moves the basis into
/// the basis into a new shared_ptr. The additional arguments can be the communication /// a new shared_ptr.
/// object. If not provided, a new communication is created. template <class RB_, class CB_>
template <class RB_, class CB_, class... Args> BiLinearForm(RB_&& rowBasis, CB_&& colBasis)
BiLinearForm(RB_&& rowBasis, CB_&& colBasis, Args&&... args) : Super(FWD(rowBasis), FWD(colBasis))
: Super(FWD(rowBasis), FWD(colBasis), FWD(args)...)
{ {
operators_.init(*Super::rowBasis(), *Super::colBasis()); operators_.init(*Super::rowBasis(), *Super::colBasis());
} }
...@@ -123,16 +122,16 @@ namespace AMDiS ...@@ -123,16 +122,16 @@ namespace AMDiS
template <class RB, class CB, class... Args> template <class RB, class CB>
BiLinearForm(RB&&, CB&&, Args&&...) BiLinearForm(RB&&, CB&&)
-> BiLinearForm<Underlying_t<RB>, Underlying_t<CB>>; -> BiLinearForm<Underlying_t<RB>, Underlying_t<CB>>;
#endif #endif
template <class RB, class CB, class... Args> template <class RB, class CB, class... Args>
auto makeBiLinearForm(RB&& rowBasis, CB&& colBasis, Args&&... args) auto makeBiLinearForm(RB&& rowBasis, CB&& colBasis)
{ {
using BLF = BiLinearForm<Underlying_t<RB>, Underlying_t<CB>>; using BLF = BiLinearForm<Underlying_t<RB>, Underlying_t<CB>>;
return BLF{FWD(rowBasis), FWD(colBasis), FWD(args)...}; return BLF{FWD(rowBasis), FWD(colBasis)};
} }
} // end namespace AMDiS } // end namespace AMDiS
...@@ -18,6 +18,7 @@ install(FILES ...@@ -18,6 +18,7 @@ install(FILES
AdaptInfo.hpp AdaptInfo.hpp
AdaptInstationary.hpp AdaptInstationary.hpp
AdaptionInterface.hpp AdaptionInterface.hpp
AdaptStationary.hpp AdaptStationary.hpp
AMDiS.hpp AMDiS.hpp
Assembler.hpp Assembler.hpp
...@@ -37,13 +38,10 @@ install(FILES ...@@ -37,13 +38,10 @@ install(FILES
DOFVector.hpp DOFVector.hpp
Environment.hpp Environment.hpp
Flag.hpp Flag.hpp
GridFunctionOperator.hpp GridFunctionOperator.hpp
GridFunctions.hpp GridFunctions.hpp
Initfile.hpp Initfile.hpp
InitfileParser.hpp InitfileParser.hpp
Integrate.hpp Integrate.hpp
...@@ -55,6 +53,7 @@ install(FILES ...@@ -55,6 +53,7 @@ install(FILES
Marker.hpp Marker.hpp
MeshCreator.hpp MeshCreator.hpp
Operations.hpp Operations.hpp
OperatorList.hpp OperatorList.hpp
Output.hpp Output.hpp
...@@ -6,16 +6,35 @@ ...@@ -6,16 +6,35 @@
#include <dune/common/shared_ptr.hh> #include <dune/common/shared_ptr.hh>
#include <amdis/DataTransfer.hpp> #include <amdis/DataTransfer.hpp>
#include <amdis/DOFVectorInterface.hpp>
#include <amdis/GridTransferManager.hpp>
#include <amdis/LinearAlgebra.hpp> #include <amdis/LinearAlgebra.hpp>
#include <amdis/Observer.hpp>
#include <amdis/common/Concepts.hpp> #include <amdis/common/Concepts.hpp>
#include <amdis/common/TypeTraits.hpp> #include <amdis/common/TypeTraits.hpp>
#include <amdis/gridfunctions/GridFunction.hpp> #include <amdis/gridfunctions/GridFunction.hpp>
#include <amdis/typetree/TreePath.hpp> #include <amdis/typetree/TreePath.hpp>
namespace Dune
namespace Functions
template <class PB>
class DefaultGlobalBasis;
namespace AMDiS namespace AMDiS
{ {
// Forward declarations
template <class PB>
class ParallelGlobalBasis;
namespace event
struct preAdapt;
struct postAdapt;
template <class B> struct basisUpdate;
/// \brief The basic container that stores a base vector and a corresponding basis /// \brief The basic container that stores a base vector and a corresponding basis
/** /**
* \tparam GB Basis of the vector * \tparam GB Basis of the vector
...@@ -24,14 +43,14 @@ namespace AMDiS ...@@ -24,14 +43,14 @@ namespace AMDiS
template <class GB, class T = double> template <class GB, class T = double>
class DOFVector class DOFVector
: public VectorBase<GB, VectorBackend<BackendTraits<GB,T>>> : public VectorBase<GB, VectorBackend<BackendTraits<GB,T>>>
, public DOFVectorInterface , public Observer<GB, event::preAdapt, event::basisUpdate<GB>, event::postAdapt>
{ {
using Self = DOFVector; using Self = DOFVector;
using Super = VectorBase<GB, VectorBackend<BackendTraits<GB,T>>>; using Super = VectorBase<GB, VectorBackend<BackendTraits<GB,T>>>;
using Obs = Observer<GB, event::preAdapt, event::basisUpdate<GB>, event::postAdapt>;
public: public:
using Backend = VectorBackend<BackendTraits<GB,T>>; using Backend = VectorBackend<BackendTraits<GB,T>>;
using Comm = typename Backend::Traits::Comm;
using GlobalBasis = GB; using GlobalBasis = GB;
...@@ -41,89 +60,32 @@ namespace AMDiS ...@@ -41,89 +60,32 @@ namespace AMDiS
/// The type of the elements of the DOFVector /// The type of the elements of the DOFVector
using value_type = typename Backend::value_type; using value_type = typename Backend::value_type;
/// Defines an interface to transfer the data during grid adaption /// Wrapper for the implementation of the transfer of data during grid adaption
using DataTransfer = DataTransferInterface<Self>; using DataTransfer = DataTransferWrapper<Self>;
/// A creator for a concrete data transfer object, depending on \ref DataTransferOperation
using DataTransferFactory = AMDiS::DataTransferFactory<Self>;
public: public:
/// Constructor. Stores the shared_ptr of the basis and creates a new DataTransfer. /// Constructor. Stores the shared_ptr of the basis and creates a new DataTransfer.
DOFVector(std::shared_ptr<GlobalBasis> basis, std::shared_ptr<Comm> comm, DataTransferOperation op = DataTransferOperation::INTERPOLATE) DOFVector(std::shared_ptr<GB> basis,
: Super(basis, std::move(comm)) DataTransferOperation op = DataTransferOperation::INTERPOLATE)
, dataTransfer_(DataTransferFactory::create(op, basis)) : Super(basis)
{ , Obs(basis)
attachToGridTransfer(); , dataTransfer_(op, basis)
} , basis_(basis)
/// Forwarding constructor that wraps the arguments into shared_ptr
template <class GB_, class Comm_,
DOFVector(GB_&& basis, Comm_&& comm, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
: DOFVector(Dune::wrap_or_move(FWD(basis)), Dune::wrap_or_move(FWD(comm)), op)
{} {}
/// Constructor creating a new basis from a Dune::Functions::DefaultGlobalBasis.
// TODO(FM): Replace explicit type with concept check
DOFVector(Dune::Functions::DefaultGlobalBasis<typename GB::PreBasis>&& basis,
DataTransferOperation op = DataTransferOperation::INTERPOLATE)
: DOFVector(std::make_shared<GB>(std::move(basis)), op)
/// Construct the DOFVector from a basis by first creating a comm object. This constructor
/// Registers the basis and comm object into the \ref GridTransferManager so that it gets updated
/// on grid changes.
DOFVector(std::shared_ptr<GlobalBasis> basis, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
: DOFVector(basis, CommunicationCreator<Comm>::create(*basis), op)
GridTransferManager::attach(this->basis()->gridView().grid(), *this->comm(), [c=this->comm(), b=this->basis()]() -> void {
/// Forwarding constructor that wraps the arguments into shared_ptr
template <class GB_, template <class GB_,
REQUIRES(Concepts::Similar<GB_,GB>)> REQUIRES(Concepts::Similar<GB_,GB>)>
DOFVector(GB_&& basis, DataTransferOperation op = DataTransferOperation::INTERPOLATE) DOFVector(GB_&& basis, DataTransferOperation op = DataTransferOperation::INTERPOLATE)
: DOFVector(Dune::wrap_or_move(FWD(basis)), op) : DOFVector(Dune::wrap_or_move(FWD(basis)), op)
{} {}
/// Copy constructor
DOFVector(Self const& that)
: Super(static_cast<Super const&>(that))
, dataTransfer_(that.dataTransfer_)
/// Move constructor
DOFVector(Self&& that)
: Super(static_cast<Super&&>(that))
, dataTransfer_(std::move(that.dataTransfer_))
/// Destructor
~DOFVector() override
/// Copy assignment operator
Self& operator=(Self const& that)
dataTransfer_ = that.dataTransfer_;
return *this;
/// Move assignment
Self& operator=(Self&& that) = default;
void resize() override
template <class TreePath = RootTreePath> template <class TreePath = RootTreePath>
auto child(TreePath const& path = {}) auto child(TreePath const& path = {})
...@@ -175,71 +137,68 @@ namespace AMDiS ...@@ -175,71 +137,68 @@ namespace AMDiS
/// Return the associated DataTransfer object /// Return the associated DataTransfer object
std::shared_ptr<DataTransfer const> dataTransfer() const