Commit 5d6a1f99 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Rewriting adaptive grid

parent c03faa46
......@@ -26,20 +26,35 @@ dune-2.6 ubuntu-18.04 clang-6-17:
- dunecontrol --current make -j4 examples
dune-git debian-10 gcc-8-17:
image: registry.dune-project.org/docker/ci/dune:git-debian-10-gcc-8-17
.dune-git:
before_script:
- . /duneci/bin/duneci-init-job
- duneci-install-module https://gitlab.dune-project.org/core/dune-common.git
- duneci-install-module https://gitlab.dune-project.org/core/dune-geometry.git
- duneci-install-module https://gitlab.dune-project.org/core/dune-localfunctions.git
- duneci-install-module https://gitlab.dune-project.org/staging/dune-uggrid.git
- duneci-install-module https://gitlab.dune-project.org/core/dune-grid.git
- duneci-install-module https://gitlab.dune-project.org/core/dune-istl.git
- duneci-install-module https://gitlab.dune-project.org/staging/dune-typetree.git
- duneci-install-module https://gitlab.dune-project.org/staging/dune-functions.git
script:
- duneci-standard-test
- dunecontrol --current make -j4 examples
dune-git debian-10 gcc-8-17:
extends: .dune-git
image: registry.dune-project.org/docker/ci/debian:10
variables:
DUNECI_TOOLCHAIN: gcc-8-17
dune-git debian-9 gcc-6-14:
image: registry.dune-project.org/docker/ci/dune:git-debian-9-gcc-6-14
script:
- duneci-standard-test
- dunecontrol --current make -j4 examples
extends: .dune-git
image: registry.dune-project.org/docker/ci/debian:9
variables:
DUNECI_TOOLCHAIN: gcc-6-14
dune-git ubuntu-18.04 clang-6-17:
image: registry.dune-project.org/docker/ci/dune:git-ubuntu-18.04-clang-6-17
script:
- duneci-standard-test
- dunecontrol --current make -j4 examples
extends: .dune-git
image: registry.dune-project.org/docker/ci/ubuntu:18.04
variables:
DUNECI_TOOLCHAIN: clang-6-17
......@@ -8,103 +8,132 @@
#include <type_traits>
#include <utility>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/common/timer.hh>
#include <dune/common/version.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/geometry/type.hh>
#include <dune/grid/common/backuprestore.hh>
#include <dune/grid/common/capabilities.hh>
#include <dune/grid/common/grid.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <amdis/common/ConceptsBase.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/DefaultGridView.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)
{
grid.loadBalance(FWD(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.
// forward declaration
template <class HostGrid>
class AdaptiveGridFamily;
/// \brief Wrapper class for Dune-grids that allows automatic signalling of events
/// during grid adaptation.
/**
* 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>
*
* \tparam HG Host grid to be wrapped. Must implement the dune grid interface.
**/
template <class HG>
class AdaptiveGrid
: public Signals<event::preAdapt, event::adapt, event::postAdapt>
: public Dune::GridDefaultImplementation<
HG::dimension, HG::dimensionworld, typename HG::ctype, AdaptiveGridFamily<HG> >
, public Signals<event::preAdapt, event::adapt, event::postAdapt>
{
using Self = AdaptiveGrid<Grid>;
using Element = typename Grid::template Codim<0>::Entity;
struct HiddenStruct {};
using Self = AdaptiveGrid<HG>;
public:
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;
using HostGrid = HG;
using GridFamily = AdaptiveGridFamily<HG>;
using Traits = typename GridFamily::Traits;
using Element = typename Traits::template Codim<0>::Entity;
/// Constructor that may only be called indirectly via the instance method
explicit AdaptiveGrid(std::shared_ptr<HostGrid> grid, HiddenStruct)
: hostGrid_(std::move(grid))
public:
template <class HostGrid_,
Dune::disableCopyMove<Self, HostGrid_> = 0>
explicit AdaptiveGrid(HostGrid_&& hostGrid)
: hostGrid_(wrap_or_share(FWD(hostGrid)))
{}
/// 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.");
}
/// Return the underlying grid
std::shared_ptr<HostGrid> const& hostGrid() const { return hostGrid_; }
public:
/// Grid iterators
/// @{
/// Iterator to first entity of given codim on level
template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition>
auto lbegin(int level) const { return hostGrid_->levelGridView(level).template begin<codim,pt>(); }
/// one past the end on this level
template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition>
auto lend(int level) const { return hostGrid_->levelGridView(level).template end<codim,pt>(); }
/// Iterator to first leaf entity of given codim
template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition>
auto leafbegin() const { return hostGrid_->leafGridView().template begin<codim,pt>(); }
/// One past the end of the sequence of leaf entities
template <int codim, Dune::PartitionIteratorType pt = Dune::All_Partition>
auto leafend() const { return hostGrid_->leafGridView().template end<codim,pt>(); }
/// Obtain begin intersection iterator with respect to the level GridView
auto ilevelbegin(Element const& e) const { return hostGrid_->levelGridView(e.level()).ibegin(e); }
/// Obtain end intersection iterator with respect to the level GridView
auto ilevelend(Element const& e) const { return hostGrid_->levelGridView(e.level()).iend(e); }
/// Obtain begin intersection iterator with respect to the leaf GridView
auto ileafbegin(Element const& e) const { return hostGrid_->leafGridView().ibegin(e); }
/// Obtain end intersection iterator with respect to the leaf GridView
auto ileafend(Element const& e) const { return hostGrid_->leafGridView().iend(e); }
/// @}
/// Size methods
/// @{
/// Return maximum level defined in this grid.
int maxLevel() const { return hostGrid_->maxLevel(); }
/// Number of grid entities per level and codim
int size(int level, int codim) const { return hostGrid_->size(level, codim); }
/// Return number of leaf entities of a given codim in this process
......@@ -116,24 +145,32 @@ namespace AMDiS
/// 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(); }
/// Returns the number of boundary segments within the macro grid
std::size_t 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(); }
/// Access to index and id sets
/// @{
/// return const reference to the host grids local id set
LocalIdSet const& localIdSet() const { return hostGrid_->localIdSet(); }
/// Return const reference to the grids global id set
auto const& globalIdSet() const { return hostGrid_->globalIdSet(); }
/// 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 local id set
auto const& localIdSet() const { return hostGrid_->localIdSet(); }
/// return const reference to the host grids leaf index set
LeafIndexSet const& leafIndexSet() const { return hostGrid_->leafIndexSet(); }
/// Return const reference to the host grids level index set for level level
auto const& levelIndexSet(int level) const { return hostGrid_->levelIndexSet(level); }
/// Return const reference to the host grids leaf index set
auto const& leafIndexSet() const { return hostGrid_->leafIndexSet(); }
/// @}
/// Adaptivity and grid refinement
/// @{
/// Refines all grid elements refCount times.
void globalRefine(int refCount)
......@@ -148,7 +185,7 @@ namespace AMDiS
}
}
/// Mark entity for refinement
/// Marks an entity to be refined/coarsened in a subsequent adapt
bool mark(int refCount, Element const& e) { return hostGrid_->mark(refCount, e); }
/// Return refinement mark for entity
......@@ -161,7 +198,7 @@ namespace AMDiS
mightCoarsen_ = hostGrid_->preAdapt();
Dune::MPIHelper::getCollectiveCommunication().max(&mightCoarsen_, 1);
this->notify(event::preAdapt{mightCoarsen_});
info(2,"preAdapt needed {} seconds", t.elapsed());
info(2,"AdaptiveGrid::preAdapt needed {} seconds", t.elapsed());
return mightCoarsen_;
}
......@@ -172,6 +209,7 @@ namespace AMDiS
refined_ = hostGrid_->adapt();
Dune::MPIHelper::getCollectiveCommunication().max(&refined_, 1);
this->notify(event::adapt{mightCoarsen_ || refined_});
info(2,"AdaptiveGrid::adapt needed {} seconds", t.elapsed());
return refined_;
}
......@@ -182,126 +220,340 @@ namespace AMDiS
hostGrid_->postAdapt();
this->notify(event::postAdapt{});
changeIndex_++;
info(2,"postAdapt needed {} seconds", t.elapsed());
info(2,"AdaptiveGrid::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()
/// Returns the grid change index, see \ref changeIndex.
unsigned long changeIndex() const
{
return Impl::loadBalanceImpl(*hostGrid_);
return changeIndex_;
}
/* 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.
*/
// @}
/// Parallel data distribution and communication
/// @{
/// Return const reference to a collective communication object.
auto const& comm() const { return hostGrid_->comm(); }
/// Communicate data of level gridView
template <class DataHandle>
bool loadBalance (DataHandle& handle)
void communicate(DataHandle& handle, Dune::InterfaceType iftype,
Dune::CommunicationDirection dir, int level) const
{
return Impl::loadBalanceImpl(*hostGrid_, handle);
hostGrid_->levelGridView(level).communicate(handle,iftype,dir);
}
/// Returns the grid change index, see \ref changeIndex.
unsigned long changeIndex() const
/// Communicate data of leaf gridView
template <class DataHandle>
void communicate(DataHandle& handle, Dune::InterfaceType iftype,
Dune::CommunicationDirection dir) const
{
return changeIndex_;
hostGrid_->leafGridView().communicate(handle,iftype,dir);
}
private:
template <class Grid_>
static std::shared_ptr<Self>
instanceImpl(std::shared_ptr<Grid_> grid)
/// \brief 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()
{
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);
write_lock.unlock();
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())
{
test_exit(!std::is_const<Grid_>::value,
"No existing AdaptiveGrid found and no mutable grid passed to create a new instance");
auto ptr = std::make_shared<Self>(std::move(grid), HiddenStruct{});
read_lock.unlock();
write_lock.lock();
adaptiveGrids_.emplace_back(ptr);
return std::move(ptr);
}
else
{
return it->lock();
}
return Dune::Hybrid::ifElse(
/* if */ std::is_convertible<decltype(std::declval<HG>().loadBalance()), bool>{},
[&](auto id) { return id(hostGrid_)->loadBalance(); },
[&](auto id) { id(hostGrid_)->loadBalance(); return true; });
}
// Do-nothing specialization if argument is already an AdaptiveGrid
static std::shared_ptr<Self> instanceImpl(std::shared_ptr<Self> grid) { return grid; }
public:
/** 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)
/// \brief 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 instanceImpl(wrap_or_share(FWD(grid)));
return Dune::Hybrid::ifElse(
/* if */ std::is_convertible<decltype(std::declval<HG>().loadBalance(handle)), bool>{},
[&](auto id) { return id(hostGrid_)->loadBalance(handle); },
[&](auto id) { id(hostGrid_)->loadBalance(handle); return true; });
}
// 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_; }
/// Return size of the overlap region for a given codim on the level grid view.
int overlapSize(int level, int codim) const { return hostGrid_->levelGridView(level).overlapSize(codim); }
/// Return size of the overlap region for a given codim on the leaf grid view.
int overlapSize(int codim) const { return hostGrid_->leafGridView().overlapSize(codim); }
/// Return size of the ghost region for a given codim on the level grid view.
int ghostSize(int level, int codim) const { return hostGrid_->levelGridView(level).ghostSize(codim); }
/// Return size of the ghost region for a given codim on the leaf grid view.
int ghostSize(int codim) const { return hostGrid_->leafGridView().ghostSize(codim); }
/// @}
/// Obtain Entity from EntitySeed of the HostGrid.
template <class EntitySeed>
auto entity(EntitySeed const& seed) const { return hostGrid_->entity(seed); }
private:
/// 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()
/// 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.
/// 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_;
template <class HostGrid>
class AdaptiveGridFamily
{
public:
struct Traits : HostGrid::Traits
{
using Grid = AdaptiveGrid<HostGrid>;
using LeafGridView = Dune::GridView< AMDiS::DefaultLeafGridViewTraits<const Grid> >;
using LevelGridView = Dune::GridView< AMDiS::DefaultLevelGridViewTraits<const Grid> >;
};
};
/// Return a change index of a grid that is not an \ref AdaptiveGrid
template <class HostGrid>
unsigned long changeIndex(HostGrid const& /*hostGrid*/)
{
return 0;
}
/// Return a change index of an \ref AdaptiveGrid
template <class HostGrid>
unsigned long changeIndex(AdaptiveGrid<HostGrid> const& grid)
{
return grid.changeIndex();
}
namespace Impl
{
template <class HostGrid>
struct AdaptiveGridImpl
{
using type = AdaptiveGrid<HostGrid>;
};
template <class HostGrid>
struct AdaptiveGridImpl<AdaptiveGrid<HostGrid>>
{
using type = AdaptiveGrid<HostGrid>;
};
}
/// Always returning an \ref AdaptiveGrid. Returns the grid itself if it is
/// already an AdaptiveGrid.
template <class HostGrid>
using AdaptiveGrid_t = typename Impl::AdaptiveGridImpl<HostGrid>::type;
} // end namespace AMDiS
namespace Dune
{
/// Specialization of a \ref GridFactory to \ref AdaptiveGrid.
/// Provide a generic factory class for unstructured grids.
template <class HostGrid>
class GridFactory<AMDiS::AdaptiveGrid<HostGrid> >
: public GridFactoryInterface<AMDiS::AdaptiveGrid<HostGrid> >
{
using Self = GridFactory;
using GridType = AMDiS::AdaptiveGrid<HostGrid>;
using HostGridFactory = GridFactory<HostGrid>;
public:
using ctype = typename GridType::ctype;
enum { dim = GridType::dimension };
enum { dimworld = GridType::dimensionworld };
template <class... Args,
Dune::disableCopyMove<Self, Args...> = 0>
GridFactory(Args&&... args)
: hostFactory_(FWD(args)...)
{}
/// Insert a vertex into the coarse grid
void insertVertex(FieldVector<ctype,dimworld> const& pos) override
{
hostFactory_.insertVertex(pos);
}
template <class F, class... Args>
using HasInsertElement = decltype(std::declval<F>().insertElement(std::declval<Args>()...));
/// Insert an element into the coarse grid
void insertElement(GeometryType const& type,
std::vector<unsigned int> const& vertices) override
{
hostFactory_.insertElement(type, vertices);
}
using ElementParametrizationType = std::shared_ptr<VirtualFunction<FieldVector<ctype,dim>,FieldVector<ctype,dimworld> > >;
/// Insert a parametrized element into the coarse grid
void insertElement(GeometryType const& type,
std::vector<unsigned int> const& vertices,
ElementParametrizationType const& elementParametrization) override
{
using A0 = GeometryType;
using A1 = std::vector<unsigned int>;
using A2 = ElementParametrizationType;
Hybrid::ifElse(Std::is_detected<HasInsertElement, HostGridFactory, A0,A1,A2>{},
[&](auto id) { id(hostFactory_).insertElement(type, vertices, elementParametrization); },
[&](auto id) { AMDiS::error_exit("insertElement() not implemented for HostGrid type."); });
}
template <class F, class... Args>
using HasInsertBoundarySegment = decltype(std::declval<F>().insertBoundarySegment(std::declval<Args>()...));