From f27ab991e117cda1f2b955b7c840e7276978d565 Mon Sep 17 00:00:00 2001 From: podlesny <podlesny@zedat.fu-berlin.de> Date: Fri, 29 May 2020 19:32:15 +0200 Subject: [PATCH] . --- .../faultfactories/CMakeLists.txt | 2 + .../sparsecantorfaultfactory.hh | 246 ++++++++++++ dune/faultnetworks/localproblem.hh | 26 +- .../patchfactories/CMakeLists.txt | 2 + .../patchfactories/cellpatchfactory.hh | 176 +++++++++ src/cantorfaultnetworks/CMakeLists.txt | 4 + .../results/sparse/plot.tex | 28 ++ .../sparsecantorfaultnetwork.cc | 354 ++++++++++++++++++ .../sparsecantorfaultnetwork.parset | 23 ++ 9 files changed, 856 insertions(+), 5 deletions(-) create mode 100644 dune/faultnetworks/faultfactories/sparsecantorfaultfactory.hh create mode 100644 dune/faultnetworks/patchfactories/cellpatchfactory.hh create mode 100644 src/cantorfaultnetworks/results/sparse/plot.tex create mode 100644 src/cantorfaultnetworks/sparsecantorfaultnetwork.cc create mode 100644 src/cantorfaultnetworks/sparsecantorfaultnetwork.parset diff --git a/dune/faultnetworks/faultfactories/CMakeLists.txt b/dune/faultnetworks/faultfactories/CMakeLists.txt index 2b32532..b1baf5e 100644 --- a/dune/faultnetworks/faultfactories/CMakeLists.txt +++ b/dune/faultnetworks/faultfactories/CMakeLists.txt @@ -4,6 +4,7 @@ add_custom_target(faultnetworks_faultfactories_src SOURCES geofaultfactory.hh cantorconvergencefaultfactory.hh cantorfaultfactory.hh + sparsecantorfaultfactory.hh ) install(FILES @@ -11,4 +12,5 @@ install(FILES geofaultfactory.hh cantorconvergencefaultfactory.hh cantorfaultfactory.hh + sparsecantorfaultfactory.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/faultnetworks/faultfactories) diff --git a/dune/faultnetworks/faultfactories/sparsecantorfaultfactory.hh b/dune/faultnetworks/faultfactories/sparsecantorfaultfactory.hh new file mode 100644 index 0000000..0a88095 --- /dev/null +++ b/dune/faultnetworks/faultfactories/sparsecantorfaultfactory.hh @@ -0,0 +1,246 @@ +#ifndef SPARSE_CANTOR_FAULT_FACTORY_HH +#define SPARSE_CANTOR_FAULT_FACTORY_HH + +#include <math.h> + +#include <dune/faultnetworks/utils/debugutils.hh> +#include <dune/faultnetworks/faultfactories/oscunitcube.hh> +#include <dune/faultnetworks/levelinterfacenetwork.hh> +#include <dune/faultnetworks/interfacenetwork.hh> + +#include <dune/grid/common/mcmgmapper.hh> + + +class SparseCantorIndexHierarchy { + public: + typedef std::map<std::pair<int,int>, bool> LevelCantorIndices; + + private: + const size_t maxLevel_; + std::vector<size_t> levelN_; + std::vector<LevelCantorIndices> cantorIndexHierarchy_; + + void prolongIndices(size_t currentLevel, size_t newLevel) { + const LevelCantorIndices& indices = cantorIndexHierarchy_[currentLevel]; + LevelCantorIndices& newIndices = cantorIndexHierarchy_[newLevel]; + + size_t offset = 2*levelN_[currentLevel]; + + std::map<std::pair<int,int>, bool>::const_iterator it = indices.begin(); + std::map<std::pair<int,int>, bool>::const_iterator endIt = indices.end(); + for (; it!=endIt; ++it) { + int xID = it->first.first; + int yID = it->first.second; + newIndices[std::make_pair(xID, yID)] = true; + newIndices[std::make_pair(xID+offset, yID)] = true; + newIndices[std::make_pair(xID, yID+offset)] = true; + } + } + + void print(const LevelCantorIndices& indices) const { + std::cout << "LevelCantorIndices: " << std::endl; + + std::map<std::pair<int,int>, bool>::const_iterator it = indices.begin(); + std::map<std::pair<int,int>, bool>::const_iterator endIt = indices.end(); + for (; it!=endIt; ++it) { + std::cout << "(" << it->first.first << ", " << it->first.second << ")"<< std::endl; + } + } + + public: + SparseCantorIndexHierarchy(size_t maxLevel) : + maxLevel_(maxLevel) + { + levelN_.resize(maxLevel_+1); + cantorIndexHierarchy_.resize(maxLevel_+1); + + // init levelCantorIndices on level 1 + LevelCantorIndices& initCantorIndices = cantorIndexHierarchy_[0]; + initCantorIndices[std::make_pair(2, 1)] = true; + initCantorIndices[std::make_pair(2, 3)] = true; + initCantorIndices[std::make_pair(2, 5)] = true; + initCantorIndices[std::make_pair(2, 7)] = true; + + initCantorIndices[std::make_pair(1, 2)] = true; + initCantorIndices[std::make_pair(3, 2)] = true; + initCantorIndices[std::make_pair(5, 2)] = true; + initCantorIndices[std::make_pair(7, 2)] = true; + + initCantorIndices[std::make_pair(1, 4)] = true; + initCantorIndices[std::make_pair(4, 1)] = true; + + for (size_t i=0; i<maxLevel_; i++) { + levelN_[i] = std::pow(4, i+1); + prolongIndices(i, i+1); + } + } + + const LevelCantorIndices& levelCantorIndices(size_t i) const { + return cantorIndexHierarchy_[i]; + } + + size_t levelN(const size_t i) const { + return levelN_[i]; + } +}; + + +template <class GridType> +class SparseCantorFaultFactory +{ + //! Parameter for mapper class + template<int dim> + struct FaceMapperLayout + { + bool contains (Dune::GeometryType gt) + { + return gt.dim() == dim-1; + } + }; + + protected: + typedef OscUnitCube<GridType, 2> GridOb; + static const int dimworld = GridType::dimensionworld ; + static const int dim = GridType::dimension; + typedef typename GridType::ctype ctype; + typedef typename GridType::LevelGridView GV; + + typedef typename Dune::MultipleCodimMultipleGeomTypeMapper<GV, FaceMapperLayout > FaceMapper; + + private: + const size_t minCantorLevel_; + const size_t maxCantorLevel_; + const SparseCantorIndexHierarchy cantorIndexHierarchy_; + + std::shared_ptr<GridType> grid_; + std::vector<double> levelResolutions_; + + std::shared_ptr<InterfaceNetwork<GridType>> interfaceNetwork_; + +private: + typename std::enable_if<!std::numeric_limits<ctype>::is_integer, bool>::type almost_equal(ctype x, ctype y, int ulp) const { + return std::abs(x-y) < std::numeric_limits<ctype>::epsilon() * std::abs(x+y) * ulp || std::abs(x-y) < std::numeric_limits<ctype>::min(); + } + + bool computeID(Dune::FieldVector<ctype, dimworld> vertex, const size_t gridN, std::pair<size_t, size_t>& IDs) const { + ctype xID = vertex[0]*gridN*2; + ctype yID = vertex[1]*gridN*2; + + ctype x = 0; + ctype y = 0; + + bool xIsInt = almost_equal(std::modf(xID, &x), 0.0, 2); + bool yIsInt = almost_equal(std::modf(yID, &y), 0.0, 2); + + if (xIsInt && yIsInt) { + IDs = std::make_pair((size_t) x, (size_t) y); + return true; + } + + if (xIsInt) { + y = std::ceil(yID); + if ((size_t) y % 2==0) + y = y-1; + IDs = std::make_pair((size_t) x, (size_t) y); + return true; + } + + if (yIsInt) { + x = std::ceil(xID); + if ((size_t) x % 2==0) + x = x-1; + IDs = std::make_pair((size_t) x, (size_t) y); + return true; + } + + return false; + } + +public: + //setup + SparseCantorFaultFactory(const size_t minCantorLevel, const size_t maxCantorLevel) : + minCantorLevel_(minCantorLevel), + maxCantorLevel_(maxCantorLevel), + cantorIndexHierarchy_(maxCantorLevel_-1), + interfaceNetwork_(nullptr) + { + Dune::UGGrid<dim>::setDefaultHeapSize(4000); + GridOb unitCube(std::pow(4, minCantorLevel_)); + grid_ = unitCube.grid(); + + size_t maxLevel = 2*(maxCantorLevel_-minCantorLevel_); + + grid_->globalRefine(maxLevel); + + // init interface network + interfaceNetwork_ = std::make_shared<InterfaceNetwork<GridType>>(*grid_); + + for (size_t i=minCantorLevel_; i<=maxCantorLevel_; i++) { + if (i>minCantorLevel_) { + interfaceNetwork_->prolongLevel(2*(i-minCantorLevel_)-2, 2*(i-minCantorLevel_)-1); + interfaceNetwork_->prolongLevel(2*(i-minCantorLevel_)-1, 2*(i-minCantorLevel_)); + } + + const size_t cantorResolution = cantorIndexHierarchy_.levelN(i-1); + const typename SparseCantorIndexHierarchy::LevelCantorIndices& levelCantorIndices = cantorIndexHierarchy_.levelCantorIndices(i-1); + std::shared_ptr<LevelInterfaceNetwork<GV>> levelInterfaceNetwork = interfaceNetwork_->levelInterfaceNetworkPtr(2*(i-minCantorLevel_)); + + const GV& gridView = levelInterfaceNetwork->levelGridView(); + FaceMapper intersectionMapper(gridView); + std::vector<bool> intersectionHandled(intersectionMapper.size(),false); + + for (const auto& elem:elements(gridView)) { + for (const auto& isect:intersections(gridView, elem)) { + + if (intersectionHandled[intersectionMapper.subIndex(elem, isect.indexInInside(),1)]) + continue; + + intersectionHandled[intersectionMapper.subIndex(elem, isect.indexInInside(),1)]=true; + + if (isect.boundary()) + continue; + + std::pair<size_t, size_t> intersectionID; + bool isAdmissibleID = computeID(isect.geometry().center(), cantorResolution, intersectionID); + + if (isAdmissibleID && levelCantorIndices.count(intersectionID)) { + levelInterfaceNetwork->addIntersection(isect, i); + } + } + } + } + } + + /* + void prolongToAll() { + // prolong all faults to all subsequent levels + for (int i=maxLevel_-1; i>=0; i--) { + if (interfaceNetwork_->size(i)>0) { + std::set<int> toLevels; + for (size_t j=i+1; j<=maxLevel_; j++) { + toLevels.insert(j); + } + + interfaceNetwork_->prolongLevelInterfaces(i, toLevels); + } + } + interfaceNetwork_->build(); + }*/ + + /*void prolongToAll() { + // prolong all faults to all subsequent levels + for (int i=interfaceNetwork_->size()-1; i>=0; i--) { + interfaceNetwork_->prolongLevelInterfaces(i, maxLevel_); + } + interfaceNetwork_->build(); + }*/ + + const GridType& grid() const { + return *grid_; + } + + InterfaceNetwork<GridType>& interfaceNetwork() { + return *interfaceNetwork_; + } +}; +#endif diff --git a/dune/faultnetworks/localproblem.hh b/dune/faultnetworks/localproblem.hh index 3da7063..7ae7a56 100644 --- a/dune/faultnetworks/localproblem.hh +++ b/dune/faultnetworks/localproblem.hh @@ -7,6 +7,7 @@ #include <dune/common/timer.hh> #include <dune/istl/matrixindexset.hh> +#include <dune/matrix-vector/axpy.hh> #include <dune/istl/umfpack.hh> @@ -113,12 +114,12 @@ public: return localMat; } - /* void setRhs(const RangeType& d){ - localRhs.resize(d.size()); + /* void setRhs(const RangeType& d){ + localRhs.resize(d.size()); - for (size_t i=0; i<localRhs.size(); ++i) { - localRhs[i] = d[i]; - } + for (size_t i=0; i<localRhs.size(); ++i) { + localRhs[i] = d[i]; + } }*/ void updateBoundary(const RangeType& boundary){ @@ -128,6 +129,15 @@ public: } } + void updateVector(const RangeType& update, RangeType& vec){ + assert(update.size()==dofs_.size()); + + for (size_t i=0; i<localRhs.size(); i++) { + if (!boundaryDofs_[i][0]) + vec[dofs_[i]] = update[i]; + } + } + void updateRhs(const RangeType& rhs){ for (size_t i=0; i<localRhs.size(); i++) { if (!boundaryDofs_[i][0]) @@ -180,6 +190,12 @@ public: } } + + void localResidual(const DomainType& x, RangeType& res){ + res = localRhs; + Dune::MatrixVector::subtractProduct(res, localMat, x); + } + }; #endif diff --git a/dune/faultnetworks/patchfactories/CMakeLists.txt b/dune/faultnetworks/patchfactories/CMakeLists.txt index 146bebe..b416642 100644 --- a/dune/faultnetworks/patchfactories/CMakeLists.txt +++ b/dune/faultnetworks/patchfactories/CMakeLists.txt @@ -1,9 +1,11 @@ #add_subdirectory("test" EXCLUDE_FROM_ALL) add_custom_target(faultnetworks_patchfactories_src SOURCES + cellpatchfactory.hh supportpatchfactory.hh ) install(FILES + cellpatchfactory.hh supportpatchfactory.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/faultnetworks/patchfactories) diff --git a/dune/faultnetworks/patchfactories/cellpatchfactory.hh b/dune/faultnetworks/patchfactories/cellpatchfactory.hh new file mode 100644 index 0000000..530c462 --- /dev/null +++ b/dune/faultnetworks/patchfactories/cellpatchfactory.hh @@ -0,0 +1,176 @@ +#ifndef CELL_PATCH_FACTORY_HH +#define CELL_PATCH_FACTORY_HH + +#include <set> +#include <queue> + +#include <dune/common/bitsetvector.hh> +#include <dune/common/fvector.hh> + +#include <dune/fufem/referenceelementhelper.hh> +#include <dune/grid/common/mcmgmapper.hh> + +#include <dune/faultnetworks/hierarchicleveliterator.hh> + +template <class BasisType> +class CellPatchFactory +{ + protected: + typedef typename BasisType::GridView GV; + typedef typename BasisType::LocalFiniteElement LFE; + + typedef typename GV::Grid GridType; + typedef typename GridType::ctype ctype; + static const int dim = GridType::dimension; + typedef typename GridType::template Codim<0>::Entity Element; + typedef typename GridType::template Codim<dim>::Entity Vertex; + + typedef HierarchicLevelIterator<GridType> HierarchicLevelIteratorType; + typedef typename Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType, Dune::MCMGElementLayout > ElementMapper; + + private: + const BasisType& basis_; + const int level_; + + const GridType& grid; + const GV& gridView; + + std::vector<std::vector<int>> cellLocalToGlobal; + std::vector<Dune::BitSetVector<1>> cellBoundaryDofs; + std::vector<std::vector<Element>> cellElements; + + ElementMapper mapper; + +private: + void print(const Dune::BitSetVector<1>& x, std::string message){ + std::cout << message << std::endl; + for (size_t i=0; i<x.size(); i++) { + std::cout << x[i][0] << " "; + } + std::cout << std::endl << std::endl; + } + + bool containsInsideSubentity(const Element& elem, const typename GridType::LevelIntersection& intersection, int subEntity, int codim) { + return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, codim); + } + +public: + //setup + CellPatchFactory(const BasisType& basis, const int level) : + basis_(basis), + level_(level), + grid(basis_.getGridView().grid()), + gridView(basis_.getGridView()), + mapper(grid, level_) { + + const auto& interfaces = basis_.faultNetwork(); + + cellElements.resize(0); + + Dune::BitSetVector<1> visited(mapper.size()); + visited.unsetAll(); + + std::queue<Element> cellSeeds; + cellSeeds.push(*gridView. template begin <0>()); + size_t cellCount = 0; + + while (!cellSeeds.empty()) { + const auto cellSeed = cellSeeds.front(); + cellSeeds.pop(); + + if (visited[mapper.index(cellSeed)][0] == true) + continue; + + cellCount++; + + cellElements.resize(cellCount); + auto& elements = cellElements.back(); + + std::queue<Element> elementQueue; + elementQueue.push(cellSeed); + + + cellLocalToGlobal.resize(cellCount); + cellBoundaryDofs.resize(cellCount); + std::set<int> dofs; + std::set<int> boundaryDofs; + + while (!elementQueue.empty()) { + const auto elem = elementQueue.front(); + const auto elemID = mapper.index(elem); + elementQueue.pop(); + + if (visited[elemID][0] == true) + continue; + + elements.push_back(elem); + visited[elemID][0] = true; + + const LFE& lfe = basis_.getLocalFiniteElement(elem); + // insert dofs + for (size_t i=0; i<lfe.localBasis().size(); ++i) { + int dofIndex = basis_.index(elem, i); + dofs.insert(dofIndex); + } + + for (const auto& is : intersections(gridView, elem)) { + if (is.neighbor()) { + const auto neighbor = is.outside(); + const auto neighborId = mapper.index(neighbor); + + if (!(visited[neighborId][0])) { + if (interfaces.isInterfaceIntersection(is)) { + cellSeeds.push(neighbor); + } else { + elementQueue.push(neighbor); + } + } + } + + // set boundaryDofs + if (is.boundary()) { + const auto& localCoefficients = lfe.localCoefficients(); + + for (size_t i=0; i<localCoefficients.size(); i++) { + auto entity = localCoefficients.localKey(i).subEntity(); + auto codim = localCoefficients.localKey(i).codim(); + + if (containsInsideSubentity(elem, is, entity, codim)) + boundaryDofs.insert(basis_.index(elem, i)); + } + } + } + } + + auto& localToGlobal = cellLocalToGlobal.back(); + auto& localBoundaryDofs = cellBoundaryDofs.back(); + + localToGlobal.resize(dofs.size()); + localBoundaryDofs.resize(dofs.size()); + localBoundaryDofs.unsetAll(); + + size_t i=0; + for (const auto& dof : dofs) { + localToGlobal[i] = dof; + + if (boundaryDofs.count(dof)) { + localBoundaryDofs[i][0] = true; + } + i++; + } + } + } + + auto& getLocalToGlobal() { + return cellLocalToGlobal; + } + + auto& getRegionElements() { + return cellElements; + } + + auto& getBoundaryDofs() { + return cellBoundaryDofs; + } +}; +#endif diff --git a/src/cantorfaultnetworks/CMakeLists.txt b/src/cantorfaultnetworks/CMakeLists.txt index a025875..6435f75 100644 --- a/src/cantorfaultnetworks/CMakeLists.txt +++ b/src/cantorfaultnetworks/CMakeLists.txt @@ -1,10 +1,14 @@ add_executable("cantorfaultnetwork" cantorfaultnetwork.cc) target_link_dune_default_libraries("cantorfaultnetwork") +add_executable("sparsecantorfaultnetwork" sparsecantorfaultnetwork.cc) +target_link_dune_default_libraries("sparsecantorfaultnetwork") + add_executable("cantorconvergence" cantorconvergence.cc) target_link_dune_default_libraries("cantorconvergence") add_custom_target(cantorfaultnetworks_srcs SOURCES cantorfaultnetwork.parset +sparsecantorfaultnetwork.parset cantorconvergence.parset ) diff --git a/src/cantorfaultnetworks/results/sparse/plot.tex b/src/cantorfaultnetworks/results/sparse/plot.tex new file mode 100644 index 0000000..5c038a7 --- /dev/null +++ b/src/cantorfaultnetworks/results/sparse/plot.tex @@ -0,0 +1,28 @@ +\documentclass[a4paper]{article} + +\usepackage{tikz} + + +\begin{document} + + \begin{figure} + \def\scale{5} + \input{levelinterfacenetwork_0.tikz} + \end{figure} + + \begin{figure} + \def\scale{5} + \input{levelinterfacenetwork_2.tikz} + \end{figure} + + \begin{figure} + \def\scale{5} + \input{levelinterfacenetwork_4.tikz} + \end{figure} + +% \begin{figure} +% \def\scale{5} +% \input{levelinterfacenetwork_6.tikz} +% \end{figure} + +\end{document} \ No newline at end of file diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc new file mode 100644 index 0000000..3441ee8 --- /dev/null +++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.cc @@ -0,0 +1,354 @@ +#include <config.h> + +#include <cmath> +#include <iostream> +#include <sstream> +#include <fstream> +#include <string> + +/* +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wunused-parameter" +#include <dune/grid/uggrid/ugwrapper.hh> +#pragma GCC diagnostic pop +*/ + +// dune-common includes +#include <dune/common/fvector.hh> +#include <dune/common/fmatrix.hh> +#include <dune/common/shared_ptr.hh> +#include <dune/common/parametertree.hh> +#include <dune/common/parametertreeparser.hh> +#include <dune/common/timer.hh> + +// dune-istl includes +#include <dune/istl/matrix.hh> +#include <dune/istl/bcrsmatrix.hh> + +// dune-fufem includes +#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh> + +// dune-grid includes +#include <dune/grid/uggrid.hh> +#include <dune/grid/io/file/vtk/vtkwriter.hh> + +// dune-solver includes +#include <dune/solvers/norms/energynorm.hh> +#include <dune/solvers/norms/twonorm.hh> + +// dune-faultnetworks includes +#include <dune/faultnetworks/utils/debugutils.hh> +#include <dune/faultnetworks/utils/matrixreader.hh> +#include <dune/faultnetworks/utils/vectorreader.hh> +#include <dune/faultnetworks/utils/levelinterfacenetworkwriter.hh> + +#include <dune/faultnetworks/preconditioners/multilevelpatchpreconditioner.hh> +#include <dune/faultnetworks/solvers/osccgsolver.hh> + +#include <dune/faultnetworks/assemblers/intersectionjumpmassassembler.hh> +#include <dune/faultnetworks/assemblers/osclocalassembler.hh> + +#include <dune/faultnetworks/oscrhs.hh> +#include <dune/faultnetworks/oscdata.hh> +#include <dune/faultnetworks/oscdatahierarchy.hh> + +#include <dune/faultnetworks/faultp1nodalbasis.hh> +#include <dune/faultnetworks/interfacenetwork.hh> +#include <dune/faultnetworks/levelinterfacenetwork.hh> +#include <dune/faultnetworks/levelinterfacenetworkproblem.hh> +#include <dune/faultnetworks/dgmgtransfer.hh> + +#include <dune/faultnetworks/faultfactories/sparsecantorfaultfactory.hh> + +int main(int argc, char** argv) { try +{ + MPIHelper::instance(argc, argv); + + const int dim = 2; + //typedef double ctype; + + + typedef typename Dune::FieldVector<double, dim> WorldVectorType; + typedef typename Dune::BlockVector<Dune::FieldVector<double,1>> VectorType; + typedef typename Dune::BitSetVector<1> BitVector; + typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>> MatrixType; + + + #if HAVE_UG + typedef typename Dune::UGGrid<dim> GridType; + #else + #error No UG! + #endif + + typedef typename GridType::LevelGridView GV; + typedef FaultP1NodalBasis<GV, double> DGBasis; + typedef typename DGBasis::LocalFiniteElement DGFE; + + typedef OscLocalAssembler<GridType, DGFE, DGFE> LocalOscAssembler; + typedef IntersectionJumpMassAssembler<GV, DGFE, DGFE> LocalInterfaceAssembler; + + // parse parameter file + ParameterTree parameterSet; + if (argc==2) + ParameterTreeParser::readINITree(argv[1], parameterSet); + else + ParameterTreeParser::readINITree("sparsecantorfaultnetwork.parset", parameterSet); + + const std::string path = parameterSet.get<std::string>("path"); + const std::string resultPath = parameterSet.get<std::string>("resultPath"); + + int problemCount = 0; + while (parameterSet.hasSub("problem" + std::to_string(problemCount))) { + + ParameterTree& problemParameters = parameterSet.sub("problem" + std::to_string(problemCount)); + + const std::string oscDataFile = problemParameters.get<std::string>("oscDataFile"); + + const size_t minCantorLevel = problemParameters.get<size_t>("coarseCantorLevel"); + const size_t fineCantorLevel = problemParameters.get<size_t>("fineCantorLevel"); + const size_t maxCantorLevel = problemParameters.get<size_t>("maxCantorLevel"); + + const double c = problemParameters.get<double>("penaltyFactor"); + const int patchDepth = problemParameters.get<int>("patchDepth"); + + const int maxIterations = problemParameters.get<int>("maxIterations"); + const double solverTolerance = problemParameters.get<double>("tol"); + const bool nestedIteration = problemParameters.get<bool>("nestedIteration"); + + const int fineGridN = std::pow(4, fineCantorLevel); + + std::ofstream out(resultPath + oscDataFile + "_" + std::to_string(minCantorLevel) + "_" + std::to_string(fineCantorLevel) + "_" + std::to_string(maxCantorLevel) + "_mult.log"); + std::streambuf *coutbuf = std::cout.rdbuf(); //save old buf + std::cout.rdbuf(out.rdbuf()); //redirect std::cout to out.txt + + + std::cout << std::endl; + std::cout << "Parameter file read successfully!" << std::endl; + + typedef Dune::Matrix<Dune::FieldMatrix<double,1,1>, std::allocator<Dune::FieldMatrix<double,1,1>>> OscMatrixType; + OscMatrixType matrix; + + // interface integral jump penalty, B: + OscMatrixType B(2,2); + B[0][0] = 1; + B[1][1] = 1; + B[0][1] = -1; + B[1][0] = -1; + + MatrixReader<OscMatrixType> matrixReader(path + oscDataFile); + matrixReader.read(matrix); + + const int oscGridN = matrix.N(); + const int oscGridLevelIdx = std::round(std::log2(oscGridN)) - minCantorLevel; + + if (oscGridN>fineGridN) + DUNE_THROW(Dune::Exception, "Provided oscData too large!"); + + std::cout << "OscData file read successfully!" << std::endl << std::endl; + + const int fineLevelIdx = 2*(fineCantorLevel - minCantorLevel); + const int exactLevelIdx = 2*(maxCantorLevel - minCantorLevel); + + // build multilevel cantor fault network + SparseCantorFaultFactory<GridType> faultFactory(minCantorLevel, maxCantorLevel); + const auto& interfaceNetwork = faultFactory.interfaceNetwork(); + + if (problemCount==0) { + std::vector<int> writeLevels {0, 2, 4, 6}; + for (size_t i=0; i<writeLevels.size(); i++) { + int writeLevel = writeLevels[i]; + bool writeGrid = !(writeLevel==8); + + LevelInterfaceNetworkWriter networkWriter(resultPath + "levelinterfacenetwork_" + std::to_string(writeLevel) + ".tikz"); + networkWriter.write(interfaceNetwork.levelInterfaceNetwork(writeLevel), writeGrid); + } + } + + const GridType& grid = faultFactory.grid(); + + OscData<GridType, OscMatrixType> oscData(grid, oscGridLevelIdx); + oscData.set(matrix); + OscDataHierarchy<GridType, OscMatrixType> oscDataHierarchy(oscData); + + Timer totalTimer; + totalTimer.reset(); + totalTimer.start(); + + + // ---------------------------------------------------------------- + // --- compute initial iterate --- + // ---------------------------------------------------------------- + VectorType coarseSol; + std::shared_ptr<DGBasis> coarseBasis; + + OscRhs<WorldVectorType, Dune::FieldVector<double,1>> f; + L2FunctionalAssembler<GridType, DGFE> l2FunctionalAssembler(f); + + std::cout << std::endl; + std::cout << "Computing initial iterate!" << std::endl; + + if (!nestedIteration) { + const LevelInterfaceNetwork<GV>& coarseLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(0); + LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> coarseGlobalAssembler(coarseLevelInterfaceNetwork); + LocalOscAssembler coarseLocalAssembler(oscDataHierarchy[0]->data(), oscDataHierarchy[0]->mapper()); + + std::vector<std::shared_ptr<LocalInterfaceAssembler>> coarseLocalInterfaceAssemblers(coarseLevelInterfaceNetwork.size()); + for (size_t i=0; i<coarseLocalInterfaceAssemblers.size(); i++) { + const int k = minCantorLevel + coarseLevelInterfaceNetwork.getIntersectionsLevels().at(i); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1); + coarseLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); + } + + coarseGlobalAssembler.assemble(coarseLocalAssembler, coarseLocalInterfaceAssemblers, l2FunctionalAssembler); + coarseGlobalAssembler.solve(); + + coarseSol = coarseGlobalAssembler.getSol(); + coarseBasis = coarseGlobalAssembler.basis(); + } else { + //TODO + const LevelInterfaceNetwork<GV>& nLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(fineLevelIdx-1); + coarseBasis = std::make_shared<DGBasis>(nLevelInterfaceNetwork); + + coarseSol.resize(coarseBasis->size()); + + std::stringstream solFilename; + solFilename << resultPath << "solutions/level_" << std::to_string(fineCantorLevel-1); + + std::ifstream file(solFilename.str().c_str(), std::ios::in|std::ios::binary); + if (not(file)) + DUNE_THROW(SolverError, "Couldn't open file " << solFilename.str() << " for reading"); + + Dune::MatrixVector::Generic::readBinary(file, coarseSol); + + file.close(); + } + + // ---------------------------------------------------------------- + // --- compute exact solution --- + // ---------------------------------------------------------------- + std::cout << std::endl; + std::cout << "Computing exact solution!" << std::endl; + + const LevelInterfaceNetwork<GV>& exactLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(exactLevelIdx); + LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> exactGlobalAssembler(exactLevelInterfaceNetwork); + LocalOscAssembler exactLocalAssembler(oscDataHierarchy[exactLevelIdx]->data(), oscDataHierarchy[exactLevelIdx]->mapper()); + + std::vector<std::shared_ptr<LocalInterfaceAssembler>> exactLocalInterfaceAssemblers(exactLevelInterfaceNetwork.size()); + for (size_t i=0; i<exactLocalInterfaceAssemblers.size(); i++) { + const int k = minCantorLevel + exactLevelInterfaceNetwork.getIntersectionsLevels().at(i); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1); + exactLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); + } + + exactGlobalAssembler.assemble(exactLocalAssembler, exactLocalInterfaceAssemblers, l2FunctionalAssembler); + exactGlobalAssembler.solve(); + + const VectorType& exactSol = exactGlobalAssembler.getSol(); + std::shared_ptr<DGBasis> exactBasis = exactGlobalAssembler.basis(); + + EnergyNorm<MatrixType, VectorType> exactEnergyNorm(exactGlobalAssembler.matrix()); + + // ---------------------------------------------------------------- + // --- set up cg solver --- + // ---------------------------------------------------------------- + std::cout << std::endl; + std::cout << "---------------------------------------------" << std::endl; + std::cout << "Setting up the fine level CG solver:" << std::endl; + std::cout << "---------------------------------------------" << std::endl << std::endl; + + const LevelInterfaceNetwork<GV>& fineLevelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(fineLevelIdx); + + LevelInterfaceNetworkProblem<DGBasis, MatrixType, VectorType> fineGlobalAssembler(fineLevelInterfaceNetwork); + LocalOscAssembler fineLocalAssembler(oscDataHierarchy[fineLevelIdx]->data(), oscDataHierarchy[fineLevelIdx]->mapper()); + + std::vector<std::shared_ptr<LocalInterfaceAssembler>> fineLocalInterfaceAssemblers(fineLevelInterfaceNetwork.size()); + for (size_t i=0; i<fineLocalInterfaceAssemblers.size(); i++) { + const int k = minCantorLevel + fineLevelInterfaceNetwork.getIntersectionsLevels().at(i); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1); + fineLocalInterfaceAssemblers[i] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); + } + + fineGlobalAssembler.assemble(fineLocalAssembler, fineLocalInterfaceAssemblers, l2FunctionalAssembler); + fineGlobalAssembler.solve(); + + const VectorType& fineSol = fineGlobalAssembler.getSol(); + const VectorType& fineRhs = fineGlobalAssembler.rhs(); + std::shared_ptr<DGBasis> fineBasis = fineGlobalAssembler.basis(); + + EnergyNorm<MatrixType, VectorType> fineEnergyNorm(fineGlobalAssembler.matrix()); + + std::cout << "Setting up the preconditioner!" << std::endl; + + BitVector activeLevels(fineLevelIdx+1, true); + + std::vector<std::shared_ptr<LocalOscAssembler>> localAssemblers; + std::vector<std::vector<std::shared_ptr<LocalInterfaceAssembler>>> localIntersectionAssemblers; + + localAssemblers.resize(activeLevels.size()); + localIntersectionAssemblers.resize(activeLevels.size()); + + for (size_t i=0; i<activeLevels.size(); i++) { + localAssemblers[i] = std::make_shared<LocalOscAssembler>(oscDataHierarchy[i]->data(), oscDataHierarchy[i]->mapper()); + + const LevelInterfaceNetwork<GV>& levelInterfaceNetwork = interfaceNetwork.levelInterfaceNetwork(i); + std::vector<std::shared_ptr<LocalInterfaceAssembler>>& levelLocalIntersectionAssemblers = localIntersectionAssemblers[i]; + levelLocalIntersectionAssemblers.resize(levelInterfaceNetwork.size()); + for (size_t j=0; j<levelLocalIntersectionAssemblers.size(); j++) { + const int k = minCantorLevel + levelInterfaceNetwork.getIntersectionsLevels().at(j); + const double pen = (1+1.0/c) * std::pow(1.0+c, k) * std::pow(2.0, k-1); + levelLocalIntersectionAssemblers[j] = std::make_shared<LocalInterfaceAssembler>(B, 2, pen); + } + } + + typedef MultilevelPatchPreconditioner<DGBasis, LocalOscAssembler, LocalInterfaceAssembler, MatrixType, VectorType> MultilevelPatchPreconditioner; + typedef typename MultilevelPatchPreconditioner::LevelPatchPreconditionerType LevelPatchPreconditioner; + MultilevelPatchPreconditioner::Mode preconditionerMode = MultilevelPatchPreconditioner::Mode::multiplicative; + + MultilevelPatchPreconditioner preconditioner(interfaceNetwork, activeLevels, localAssemblers, localIntersectionAssemblers, preconditionerMode); + for (size_t i=0; i<preconditioner.size()-1; i++) { + LevelPatchPreconditioner& levelFaultPreconditioner = *preconditioner.levelPatchPreconditioner(i); + + levelFaultPreconditioner.setPatchDepth(patchDepth); + levelFaultPreconditioner.setBoundaryMode(LevelPatchPreconditioner::BoundaryMode::homogeneous); + } + preconditioner.build(); + + std::cout << "Setup complete, starting preconditioned cg iteration!" << std::endl; + std::cout << std::endl << std::endl; + + // set initial iterate + VectorType initialX; + DGMGTransfer<DGBasis> coarseFineTransfer(*coarseBasis, *fineBasis); + coarseFineTransfer.prolong(coarseSol, initialX); + + VectorType rhs(fineRhs); + + DGMGTransfer<DGBasis> discMGTransfer(*fineBasis, *exactBasis); + + // solve + OscCGSolver<MatrixType, VectorType, DGBasis> + solver(*fineBasis, &fineGlobalAssembler.matrix(), &initialX, &rhs, &preconditioner, + maxIterations, solverTolerance, &exactEnergyNorm, + Solver::FULL, &exactSol, &fineSol, discMGTransfer, 1.0, true); //((oscGridN+0.0)/fineGridN) + + + solver.historyBuffer_ = resultPath + "solutions/level_" + std::to_string(fineCantorLevel) + ".vec"; + + solver.check(); + solver.preprocess(); + solver.solve(); + + std::cout << std::endl; + std::cout << "Total time: " << totalTimer.elapsed() << " seconds." << std::endl; + + std::cout.rdbuf(coutbuf); //reset to standard output again + + std::cout << "Problem " << problemCount << " done!" << std::endl; + problemCount++; + } +} catch (Dune::Exception e) { + + std::cout << e << std::endl; + +} +} diff --git a/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset new file mode 100644 index 0000000..6ecfab3 --- /dev/null +++ b/src/cantorfaultnetworks/sparsecantorfaultnetwork.parset @@ -0,0 +1,23 @@ +path = ../data/ +resultPath = ../cantorfaultnetworks/results/sparse/ + + +########################################### + +[problem0] +oscDataFile = oscDataLaplace32.mat + +# level resolution in 2^(-...) +coarseCantorLevel = 1 +fineCantorLevel = 3 +maxCantorLevel = 4 + +penaltyFactor = 1 +patchDepth = 1 + +# cg solver +nestedIteration = 0 +tol = 1e-12 +maxIterations = 8 + +########################################### -- GitLab