From cafdbe777bcae7e64975ee75862af2c22e2a3ecc Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Sun, 31 Mar 2019 11:19:29 +0200
Subject: [PATCH] backup and restore for DOFVectors
---
src/amdis/BackupRestore.hpp | 177 ++++++++++++++++++
src/amdis/ProblemStat.inc.hpp | 17 +-
src/amdis/linearalgebra/DOFVectorBase.hpp | 6 +
src/amdis/linearalgebra/DOFVectorBase.inc.hpp | 76 +++++++-
test/BackupRestoreTest.cpp | 53 +++++-
test/CMakeLists.txt | 12 +-
6 files changed, 329 insertions(+), 12 deletions(-)
create mode 100644 src/amdis/BackupRestore.hpp
diff --git a/src/amdis/BackupRestore.hpp b/src/amdis/BackupRestore.hpp
new file mode 100644
index 00000000..33fad315
--- /dev/null
+++ b/src/amdis/BackupRestore.hpp
@@ -0,0 +1,177 @@
+#pragma once
+
+#include <cstdint>
+#include <fstream>
+#include <memory>
+#include <vector>
+
+#include <dune/grid/common/backuprestore.hh>
+#include <dune/grid/common/gridfactory.hh>
+
+namespace AMDiS
+{
+ // Fallback Implementation for Grids not supporting direct BackupRestore
+ template <class Grid>
+ class BackupRestoreByGridFactory
+ {
+ template <class GridView>
+ class Writer
+ {
+ enum { dim = GridView::dimension };
+ enum { dow = GridView::dimensionworld };
+
+ using ct = typename GridView::ctype;
+
+ public:
+ Writer(GridView const& gv)
+ : gridView_(gv)
+ {
+ std::int64_t num = 0;
+ indexMap_.resize(gridView_.size(dim));
+ auto const& indexSet = gridView_.indexSet();
+ for (auto const& vertex : vertices(gridView_))
+ indexMap_[indexSet.index(vertex)] = std::int64_t(num++);
+ }
+
+ void writeVertices(std::ostream& out) const
+ {
+ for (auto const& vertex : vertices(gridView_)) {
+ auto v = vertex.geometry().center();
+ out.write((char*)&v, dow*sizeof(ct));
+ }
+ }
+
+ void writeElements(std::ostream& out) const
+ {
+ auto const& indexSet = gridView_.indexSet();
+
+ std::vector<std::int64_t> connectivity;
+ connectivity.reserve(8);
+ for (auto const& e : elements(gridView_)) {
+ unsigned int id = e.type().id();
+ out.write((char*)&id, sizeof(unsigned int));
+
+ connectivity.clear();
+ for (unsigned int j = 0; j < e.subEntities(dim); ++j)
+ connectivity.emplace_back(indexMap_[indexSet.subIndex(e,j,dim)]);
+
+ out.write((char*)connectivity.data(), connectivity.size()*sizeof(std::int64_t));
+ }
+ }
+
+ private:
+ GridView gridView_;
+ std::vector<std::int64_t> indexMap_;
+ };
+
+ class Reader
+ {
+ enum { dim = Grid::dimension };
+ enum { dow = Grid::dimensionworld };
+
+ using ct = typename Grid::ctype;
+ using Factory = Dune::GridFactory<Grid>;
+ using GlobalCoordinates = Dune::FieldVector<ct,dow>;
+
+ public:
+ Reader(Factory& factory, std::istream& in)
+ : factory_(factory)
+ {
+ in.read((char*)&numElements_, sizeof(std::int64_t));
+ in.read((char*)&numVertices_, sizeof(std::int64_t));
+ }
+
+ void readVertices(std::istream& in) const
+ {
+ GlobalCoordinates p;
+ for (std::int64_t i = 0; i < numVertices_; ++i) {
+ in.read((char*)&p[0], dow*sizeof(ct));
+ factory_.insertVertex(p);
+ }
+ }
+
+ void readElements(std::istream& in) const
+ {
+ std::vector<std::int64_t> connectivity(8);
+ std::vector<unsigned int> vertices;
+ vertices.reserve(8);
+
+ for (std::int64_t i = 0; i < numElements_; ++i) {
+ unsigned int id = 0;
+ in.read((char*)&id, sizeof(unsigned int));
+
+ Dune::GeometryType type(id,dim);
+ auto refElem = Dune::referenceElement<ct,dim>(type);
+
+ in.read((char*)connectivity.data(), refElem.size(dim)*sizeof(std::int64_t));
+ vertices.clear();
+ std::copy_n(connectivity.begin(),refElem.size(dim),std::back_inserter(vertices));
+
+ factory_.insertElement(type, vertices);
+ }
+ }
+
+ private:
+ Factory& factory_;
+ std::int64_t numElements_ = 0;
+ std::int64_t numVertices_ = 0;
+ };
+
+ public:
+ /// \brief Write a hierarchic grid to disk
+ template <class GridView>
+ static void backup(GridView const& gv, std::string const& filename)
+ {
+ std::ofstream out(filename, std::ios::binary);
+ backup(gv,out);
+ }
+
+ /// \brief write a hierarchic grid into a stream
+ template <class GridView>
+ static void backup(GridView const& gv, std::ostream& out)
+ {
+ std::int32_t dim = GridView::dimension;
+ std::int32_t dow = GridView::dimensionworld;
+ std::int64_t num_elements = gv.size(0);
+ std::int64_t num_vertices = gv.size(dim);
+
+ out.write((char*)&dim, sizeof(std::int32_t));
+ out.write((char*)&dow, sizeof(std::int32_t));
+ out.write((char*)&num_elements, sizeof(std::int64_t));
+ out.write((char*)&num_vertices, sizeof(std::int64_t));
+
+ Writer writer(gv);
+ writer.writeVertices(out);
+ writer.writeElements(out);
+ }
+
+ /// \brief read a hierarchic grid from disk
+ static Grid* restore(std::string const& filename)
+ {
+ std::ifstream in(filename, std::ios::binary);
+ return restore(in);
+ }
+
+ /// \brief read a hierarchic grid from a stream
+ static Grid* restore(std::istream& in)
+ {
+ Dune::GridFactory<Grid> factory;
+
+ std::int32_t dim = 0, dow = 0;
+
+ in.read((char*)&dim, sizeof(std::int32_t));
+ in.read((char*)&dow, sizeof(std::int32_t));
+
+ assert(Grid::dimension == dim);
+ assert(Grid::dimensionworld == dow);
+
+ Reader reader(factory, in);
+ reader.readVertices(in);
+ reader.readElements(in);
+
+ std::unique_ptr<Grid> ptr(factory.createGrid());
+ return ptr.release();
+ }
+ };
+
+} // end namespace AMDiS
diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp
index 2cfa566e..142148cd 100644
--- a/src/amdis/ProblemStat.inc.hpp
+++ b/src/amdis/ProblemStat.inc.hpp
@@ -9,6 +9,7 @@
#include <dune/typetree/childextraction.hh>
#include <amdis/AdaptInfo.hpp>
+#include <amdis/BackupRestore.hpp>
#include <amdis/FileWriter.hpp>
#include <amdis/Assembler.hpp>
#include <amdis/GridFunctionOperator.hpp>
@@ -491,7 +492,12 @@ backup(AdaptInfo& adaptInfo)
auto param = Parameters::get<std::string>(name_ + "->backup filename");
std::string filename = param ? *param : name_ + "_" + std::to_string(adaptInfo.timestepNumber()) + ".backup";
- Dune::BackupRestoreFacility<Grid>::backup(*grid_, filename);
+ try {
+ Dune::BackupRestoreFacility<Grid>::backup(*grid_, filename);
+ } catch (Dune::NotImplemented const&) {
+ msg("Falling back to backup of gridview.");
+ BackupRestoreByGridFactory<Grid>::backup(gridView(), filename);
+ }
// TODO: backup data
msg("Problem backuped to file '{}'.", filename);
@@ -506,8 +512,13 @@ restore(Flag initFlag)
test_exit(filesystem::exists(filename), "Restore file '{}' not found.", filename);
// restore grid from file
- std::shared_ptr<Grid> grid(Dune::BackupRestoreFacility<Grid>::restore(filename));
- adoptGrid(grid);
+ std::shared_ptr<Grid> grid;
+ try {
+ grid.reset(Dune::BackupRestoreFacility<Grid>::restore(filename));
+ } catch (Dune::NotImplemented const&) {
+ grid.reset(BackupRestoreByGridFactory<Grid>::restore(filename));
+ }
+ adoptGrid(std::move(grid));
// create fespace
if (initFlag.isSet(INIT_FE_SPACE) || initFlag.isSet(INIT_SYSTEM))
diff --git a/src/amdis/linearalgebra/DOFVectorBase.hpp b/src/amdis/linearalgebra/DOFVectorBase.hpp
index 91e1c57d..fc293862 100644
--- a/src/amdis/linearalgebra/DOFVectorBase.hpp
+++ b/src/amdis/linearalgebra/DOFVectorBase.hpp
@@ -189,6 +189,12 @@ namespace AMDiS
/// Assemble all vector operators
void assemble();
+ /// Write DOFVector to file
+ void backup(std::string const& filename);
+
+ /// Read backup data from file
+ void restore(std::string const& filename);
+
/// Return the associated DataTransfer object
std::shared_ptr<DataTransfer const> dataTransfer() const
{
diff --git a/src/amdis/linearalgebra/DOFVectorBase.inc.hpp b/src/amdis/linearalgebra/DOFVectorBase.inc.hpp
index 4795c768..610b54b9 100644
--- a/src/amdis/linearalgebra/DOFVectorBase.inc.hpp
+++ b/src/amdis/linearalgebra/DOFVectorBase.inc.hpp
@@ -1,5 +1,9 @@
#pragma once
+#include <cstdint>
+#include <fstream>
+#include <functional>
+
#include <amdis/Assembler.hpp>
#include <amdis/LocalOperator.hpp>
#include <amdis/typetree/Traversal.hpp>
@@ -89,9 +93,79 @@ assemble()
for (auto const& element : elements(basis_->gridView())) {
localView.bind(element);
assemble(localView);
- localView.unbind(element);
+ localView.unbind();
}
finish(true);
}
+
+template <class GB, class B>
+void DOFVectorBase<GB,B>::
+backup(std::string const& filename)
+{
+ std::ofstream out(filename, std::ios::binary);
+
+ std::int64_t numElements = basis_->gridView().size(0);
+ out.write((char*)&numElements, sizeof(std::int64_t));
+
+ auto localView = basis_->localView();
+ std::vector<value_type> data;
+ for (auto const& element : elements(basis_->gridView()))
+ {
+ localView.bind(element);
+ data.clear();
+ for_each_leaf_node(localView.tree(), [&](auto const& node, auto) -> void {
+ auto const& fe = node.finiteElement();
+ std::size_t size = fe.size();
+
+ for (std::size_t i = 0; i < size; ++i)
+ data.push_back((*this)[localView.index(node.localIndex(i))]);
+ });
+
+ std::uint64_t len = data.size();
+ out.write((char*)&len, sizeof(std::uint64_t));
+ out.write((char*)data.data(), len*sizeof(value_type));
+
+ localView.unbind();
+ }
+}
+
+template <class GB, class B>
+void DOFVectorBase<GB,B>::
+restore(std::string const& filename)
+{
+ std::ifstream in(filename, std::ios::binary);
+
+ std::int64_t numElements = 0;
+ in.read((char*)&numElements, sizeof(std::int64_t));
+ assert(numElements == basis_->gridView().size(0));
+
+ // assume the order of element traversal is not changed
+ auto localView = basis_->localView();
+ std::vector<value_type> data;
+ for (auto const& element : elements(basis_->gridView()))
+ {
+ std::uint64_t len = 0;
+ in.read((char*)&len, sizeof(std::uint64_t));
+ data.resize(len);
+
+ in.read((char*)data.data(), len*sizeof(value_type));
+
+ localView.bind(element);
+ std::size_t shift = 0;
+ for_each_leaf_node(localView.tree(), [&](auto const& node, auto) -> void {
+ auto const& fe = node.finiteElement();
+ std::size_t size = fe.size();
+
+ assert(data.size() >= shift+size);
+ for (std::size_t i = 0; i < size; ++i)
+ (*this)[localView.index(node.localIndex(i))] = data[shift + i];
+
+ shift += size;
+ });
+
+ localView.unbind();
+ }
+}
+
} // end namespace AMDiS
diff --git a/test/BackupRestoreTest.cpp b/test/BackupRestoreTest.cpp
index 7471e87b..7e29077b 100644
--- a/test/BackupRestoreTest.cpp
+++ b/test/BackupRestoreTest.cpp
@@ -1,6 +1,17 @@
#include <amdis/AMDiS.hpp>
+#include <amdis/Integrate.hpp>
#include <amdis/ProblemStat.hpp>
+#if HAVE_DUNE_SPGRID
+#include <dune/grid/spgrid.hh>
+#endif
+#if HAVE_DUNE_ALUGRID
+#include <dune/alugrid/grid.hh>
+#endif
+#if HAVE_DUNE_FOAMGRID
+#include <dune/foamgrid/foamgrid.hh>
+#endif
+
#include "Tests.hpp"
using namespace AMDiS;
@@ -21,11 +32,15 @@ void test(Factory factory)
prob.initialize(INIT_ALL);
prob.globalRefine(2);
+ prob.solution(Dune::Indices::_0) << [](auto const& x) { return x; };
+ prob.solution(Dune::Indices::_1) << [](auto const& x) { return two_norm(x); };
+
num_elements = prob.grid()->size(0);
num_vertices = prob.grid()->size(Grid::dimension);
AdaptInfo adaptInfo("adapt");
prob.backup(adaptInfo);
+ prob.solutionVector()->backup("solution.backup");
}
// restore
@@ -35,6 +50,20 @@ void test(Factory factory)
AMDIS_TEST_EQ(num_elements, prob.grid()->size(0));
AMDIS_TEST_EQ(num_vertices, prob.grid()->size(Grid::dimension));
+
+ auto vec = *prob.solutionVector();
+ vec.restore("solution.backup");
+
+ prob.solution(Dune::Indices::_0) << [](auto const& x) { return x; };
+ prob.solution(Dune::Indices::_1) << [](auto const& x) { return two_norm(x); };
+
+ double error0 = std::sqrt(integrate(
+ unary_dot(prob.solution(Dune::Indices::_0) - makeDiscreteFunction(vec, Dune::Indices::_0)), prob.gridView()));
+ AMDIS_TEST(error0 < 1.e-10);
+
+ double error1 = std::sqrt(integrate(
+ pow<2>(prob.solution(Dune::Indices::_1) - makeDiscreteFunction(vec, Dune::Indices::_1)), prob.gridView()));
+ AMDIS_TEST(error1 < 1.e-10);
}
}
@@ -60,12 +89,26 @@ int main(int argc, char** argv)
Parameters::set("test->backup filename", filename);
Parameters::set("test->restore filename", filename);
+#if GRID_ID == 0
test_cube<Dune::YaspGrid<2>>();
- // test_cube<Dune::UGGrid<2>>();
-
- // test_simplex<Dune::AlbertaGrid<2,2>>();
- // test_simplex<Dune::UGGrid<2>>();
+#elif GRID_ID == 1 && HAVE_DUNE_UGGRID
+ test_cube<Dune::UGGrid<2>>();
+#elif GRID_ID == 2 && HAVE_DUNE_ALUGRID
+ test_cube<Dune::ALUGrid<2,2,Dune::cube,Dune::nonconforming>>();
+#elif GRID_ID == 3 && HAVE_DUNE_SPGRID
+ test<Dune::SPGrid<double,2>>([]() {
+ return std::unique_ptr<Dune::SPGrid<double,2>>(
+ new Dune::SPGrid<double,2>(Dune::SPDomain<double, 2>({0.0,0.0}, {1.0,1.0}), Dune::SPMultiIndex<2>({2,2}))); });
+#elif GRID_ID == 4 && HAVE_DUNE_UGGRID
+ test_simplex<Dune::UGGrid<2>>();
+#elif GRID_ID == 5 && HAVE_DUNE_ALUGRID
+ test_simplex<Dune::ALUGrid<2,2,Dune::simplex,Dune::nonconforming>>();
+#elif GRID_ID == 6 && HAVE_DUNE_ALUGRID
+ test_simplex<Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming>>();
+#endif
+ // test_simplex<Dune::AlbertaGrid<2,2>>(); // Segmentation fault
+ // test_simplex<Dune::FoamGrid<2,2>>(); // Segmentation fault
AMDiS::finalize();
- return 0;
+ return report_errors();
}
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index f5ef52b2..515a07b5 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -2,9 +2,15 @@
dune_add_test(SOURCES AdaptInfoTest.cpp
LINK_LIBRARIES amdis)
-dune_add_test(SOURCES BackupRestoreTest.cpp
- LINK_LIBRARIES amdis)
-add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 BackupRestoreTest)
+foreach(_GRID RANGE 6)
+ dune_add_test(NAME "BackupRestoreTest_${_GRID}"
+ SOURCES BackupRestoreTest.cpp
+ COMPILE_DEFINITIONS "GRID_ID=${_GRID}"
+ LABELS "BackupRestore"
+ LINK_LIBRARIES amdis)
+ add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 "BackupRestoreTest_${_GRID}")
+endforeach()
+unset(_GRID)
dune_add_test(SOURCES ConceptsTest.cpp
LINK_LIBRARIES amdis)
--
GitLab