#pragma once #include #include #include #include #include #include #include #include #include #include #include #include #include namespace AMDiS { template void ProblemStat::initialize( Flag initFlag, Self* adoptProblem, Flag adoptFlag) { // create grids if (!grid_) { if (initFlag.isSet(CREATE_MESH) || (!adoptFlag.isSet(INIT_MESH) && (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) { createGrid(); } if (adoptProblem && (adoptFlag.isSet(INIT_MESH) || adoptFlag.isSet(INIT_SYSTEM) || adoptFlag.isSet(INIT_FE_SPACE))) { adoptGrid(adoptProblem->grid_, adoptProblem->boundaryManager_); } } if (!grid_) warning("no grid created"); // create fespace if (!globalBasis_) { if (initFlag.isSet(INIT_FE_SPACE) || (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) { createGlobalBasis(); } if (adoptProblem && (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) { adoptGlobalBasis(adoptProblem->globalBasis_); } } if (!globalBasis_) warning("no globalBasis created\n"); // create system if (initFlag.isSet(INIT_SYSTEM)) createMatricesAndVectors(); if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) { systemMatrix_ = adoptProblem->systemMatrix_; solution_ = adoptProblem->solution_; rhs_ = adoptProblem->rhs_; estimates_ = adoptProblem->estimates_; } // create solver if (!linearSolver_) { if (initFlag.isSet(INIT_SOLVER)) createSolver(); if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) { test_exit(!linearSolver_, "solver already created\n"); linearSolver_ = adoptProblem->linearSolver_; } } if (!linearSolver_) { warning("no solver created\n"); } // create marker if (initFlag.isSet(INIT_MARKER)) createMarker(); if (adoptProblem && adoptFlag.isSet(INIT_MARKER)) marker_ = adoptProblem->marker_; // create file writer if (initFlag.isSet(INIT_FILEWRITER)) createFileWriter(); solution_->resizeZero(); } template void ProblemStat:: restore(Flag initFlag) { std::string grid_filename = Parameters::get(name_ + "->restore->grid").value(); std::string solution_filename = Parameters::get(name_ + "->restore->solution").value(); test_exit(filesystem::exists(grid_filename), "Restore file '{}' not found.", grid_filename); test_exit(filesystem::exists(solution_filename), "Restore file '{}' not found.", solution_filename); // TODO(SP): implement BAckupRestore independent of wrapped grid using HostGrid = typename Grid::HostGrid; // restore grid from file if (Dune::Capabilities::hasBackupRestoreFacilities::v) adoptGrid(std::shared_ptr(Dune::BackupRestoreFacility::restore(grid_filename))); else adoptGrid(std::shared_ptr(BackupRestoreByGridFactory::restore(grid_filename))); // create fespace if (initFlag.isSet(INIT_FE_SPACE) || initFlag.isSet(INIT_SYSTEM)) createGlobalBasis(); // create system if (initFlag.isSet(INIT_SYSTEM)) createMatricesAndVectors(); // create solver if (!linearSolver_ && initFlag.isSet(INIT_SOLVER)) createSolver(); // create marker if (initFlag.isSet(INIT_MARKER)) createMarker(); // create file writer if (initFlag.isSet(INIT_FILEWRITER)) createFileWriter(); solution_->resize(sizeInfo(*globalBasis_)); solution_->restore(solution_filename); } template void ProblemStat::createGrid() { Parameters::get(name_ + "->mesh", gridName_); MeshCreator creator(gridName_); grid_ = creator.create(); boundaryManager_ = std::make_shared>(grid_); if (!creator.boundaryIds().empty()) boundaryManager_->setBoundaryIds(creator.boundaryIds()); info(3,"Create grid:"); info(3,"#elements = {}" , grid_->size(0)); info(3,"#faces/edges = {}", grid_->size(1)); info(3,"#vertices = {}" , grid_->size(dim)); info(3,"overlap-size = {}", grid_->leafGridView().overlapSize(0)); info(3,"ghost-size = {}" , grid_->leafGridView().ghostSize(0)); info(3,""); } template using HasCreate = decltype(T::create(std::declval())); template void ProblemStat::createGlobalBasis() { createGlobalBasisImpl(Dune::Std::is_detected{}); initGlobalBasis(); } template void ProblemStat::createGlobalBasisImpl(std::true_type) { assert( bool(grid_) ); static_assert(std::is_same_v, ""); auto basis = Traits::create(name_, grid_->leafGridView()); globalBasis_ = std::make_shared(std::move(basis)); } template void ProblemStat::createGlobalBasisImpl(std::false_type) { error_exit("Cannot create GlobalBasis from type. Pass a BasisCreator instead!"); } template void ProblemStat::initGlobalBasis() {} template void ProblemStat::createMatricesAndVectors() { systemMatrix_ = std::make_shared(globalBasis_, globalBasis_); std::string symmetryStr = "unknown"; Parameters::get(name_ + "->symmetry", symmetryStr); systemMatrix_->setSymmetryStructure(symmetryStr); solution_ = std::make_shared(globalBasis_); rhs_ = std::make_shared(globalBasis_); auto localView = globalBasis_->localView(); Traversal::forEachNode(localView.tree(), [&,this](auto&&, auto treePath) -> void { std::string i = to_string(treePath); estimates_[i].resize(globalBasis_->gridView().indexSet().size(0)); for (std::size_t j = 0; j < estimates_[i].size(); j++) estimates_[i][j] = 0.0; // TODO: Remove when estimate() is implemented }); } template void ProblemStat::createSolver() { std::string solverName = "default"; Parameters::get(name_ + "->solver", solverName); linearSolver_ = std::make_shared(solverName, name_ + "->solver"); } template void ProblemStat::createMarker() { marker_.clear(); auto localView = globalBasis_->localView(); Traversal::forEachNode(localView.tree(), [&,this](auto&&, auto treePath) -> void { std::string componentName = name_ + "->marker[" + to_string(treePath) + "]"; if (!Parameters::get(componentName + "->strategy")) return; std::string tp = to_string(treePath); auto newMarker = EstimatorMarker::createMarker(componentName, tp, estimates_[tp], grid_); assert(bool(newMarker)); this->addMarker(std::move(newMarker)); }); } template void ProblemStat::createFileWriter() { FileWriterCreator creator(solution_, boundaryManager_); filewriter_.clear(); auto localView = globalBasis_->localView(); Traversal::forEachNode(localView.tree(), [&](auto const& /*node*/, auto treePath) -> void { std::string componentName = name_ + "->output[" + to_string(treePath) + "]"; auto format = Parameters::get>(componentName + "->format"); if (!format && to_string(treePath).empty()) { // alternative for root treepath componentName = name_ + "->output"; format = Parameters::get>(componentName + "->format"); } if (!format) return; for (std::string const& type : format.value()) { auto writer = creator.create(type, componentName, treePath); if (writer) filewriter_.push_back(std::move(writer)); } }); } // Adds a Dirichlet boundary condition template template void ProblemStat:: addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values) { static constexpr bool isValidPredicate = Concepts::Functor; static_assert( Concepts::Functor, "Function passed to addDirichletBC for `predicate` does not model the Functor concept"); static constexpr bool isValidTreePath = Concepts::ValidTreePath && Concepts::ValidTreePath; static_assert(isValidTreePath, "Invalid row and/or col treepath passed to addDirichletBC!"); if constexpr (isValidPredicate && isValidTreePath) { auto basis = Dune::Functions::subspaceBasis(*globalBasis_, makeTreePath(col)); auto valueGridFct = makeGridFunction(values, this->gridView()); constraints_.push_back(DirichletBC{ std::move(basis), {predicate}, valueGridFct}); } } // Adds a Dirichlet boundary condition template template void ProblemStat:: addDirichletBC(BoundaryType id, RowTreePath row, ColTreePath col, Values const& values) { static constexpr bool isValidTreePath = Concepts::ValidTreePath && Concepts::ValidTreePath; static_assert(isValidTreePath, "Invalid row and/or col treepath passed to addDirichletBC!"); if constexpr (isValidTreePath) { auto basis = Dune::Functions::subspaceBasis(*globalBasis_, makeTreePath(col)); auto valueGridFct = makeGridFunction(values, this->gridView()); constraints_.push_back(DirichletBC{ std::move(basis), {*boundaryManager_, id}, valueGridFct}); } } template void ProblemStat:: addPeriodicBC(BoundaryType id, WorldMatrix const& matrix, WorldVector const& vector) { auto localView = globalBasis_->localView(); auto basis = Dune::Functions::subspaceBasis(*globalBasis_, makeTreePath()); constraints_.push_back(makePeriodicBC( std::move(basis), {*boundaryManager_, id}, {matrix, vector})); } template void ProblemStat:: solve(AdaptInfo& /*adaptInfo*/, bool createMatrixData, bool storeMatrixData) { Dune::InverseOperatorResult stat; solution_->resize(); if (createMatrixData) linearSolver_->init(systemMatrix_->impl()); // solve the linear system linearSolver_->apply(solution_->impl(), rhs_->impl(), stat); if (!storeMatrixData) linearSolver_->finish(); info(2, "solution of discrete system needed {} seconds", stat.elapsed); info(1, "Residual norm: ||b-Ax|| = {}", residuum(systemMatrix_->impl(),solution_->impl(), rhs_->impl())); if (stat.reduction >= 0.0) info(2, "Relative residual norm: ||b-Ax||/||b|| = {}", stat.reduction); else info(2, "Relative residual norm: ||b-Ax||/||b|| = {}", relResiduum(systemMatrix_->impl(),solution_->impl(), rhs_->impl())); bool ignoreConverged = false; Parameters::get(name_ + "->solver->ignore converged", ignoreConverged); test_exit(stat.converged || ignoreConverged, "Could not solver the linear system!"); } template Flag ProblemStat:: markElements(AdaptInfo& adaptInfo) { Dune::Timer t; Flag markFlag = 0; for (auto& currentMarker : marker_) markFlag |= currentMarker.second->markGrid(adaptInfo); msg("markElements needed {} seconds", t.elapsed()); return markFlag; } template Flag ProblemStat:: globalCoarsen(int n) { Dune::Timer t; bool adapted = false; // TODO(FM): Find a less expensive alternative to the loop adaption for (int i = 0; i < n; ++i) { // mark all entities for coarsening for (const auto& element : elements(grid_->leafGridView())) grid_->mark(-1, element); bool adaptedInLoop = grid_->preAdapt(); adaptedInLoop |= grid_->adapt(); grid_->postAdapt(); if (!adaptedInLoop) break; else adapted = true; } msg("globalCoarsen needed {} seconds", t.elapsed()); return adapted ? MESH_ADAPTED : Flag(0); } // grid has globalRefine(int, AdaptDataHandleInterface&) template using HasGlobalRefineADHI = decltype( std::declval().globalRefine(1,std::declval())); template Flag ProblemStat:: globalRefine(int n) { Dune::Timer t; if constexpr (Dune::Std::is_detected::value) grid_->globalRefine(n, globalBasis_->globalRefineCallback()); else grid_->globalRefine(n); msg("globalRefine needed {} seconds", t.elapsed()); return n > 0 ? MESH_ADAPTED : Flag(0); } template Flag ProblemStat:: adaptGrid(AdaptInfo& /*adaptInfo*/) { Dune::Timer t; bool adapted = grid_->preAdapt(); adapted |= grid_->adapt(); grid_->postAdapt(); msg("adaptGrid needed {} seconds", t.elapsed()); return adapted ? MESH_ADAPTED : Flag(0); } template void ProblemStat:: buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector) { Dune::Timer t; Dune::Timer t2; // 0. initialize boundary condition and other constraints for (auto& bc : constraints_) bc.init(); t2.reset(); // 1. init matrix and rhs vector and initialize dirichlet boundary conditions systemMatrix_->init(); rhs_->init(sizeInfo(*globalBasis_), asmVector); // statistic about system size if (Environment::mpiSize() > 1) msg("{} local DOFs, {} global DOFs", rhs_->localSize(), rhs_->globalSize()); else msg("{} local DOFs", rhs_->localSize()); auto localMatOperators = localOperators(systemMatrix_->operators()); auto localVecOperators = localOperators(rhs_->operators()); // 2. traverse grid and assemble operators on the elements auto localView = globalBasis_->localView(); for (auto const& element : elements(gridView(), PartitionSet{})) { localView.bind(element); if (asmMatrix) systemMatrix_->assemble(localView, localView, localMatOperators); if (asmVector) rhs_->assemble(localView, localVecOperators); localView.unbind(); } // 3. finish matrix insertion and apply dirichlet boundary conditions systemMatrix_->finish(); rhs_->finish(); info(2," assemble operators needed {} seconds", t2.elapsed()); t2.reset(); solution_->resize(sizeInfo(*globalBasis_)); // 4. apply boundary condition and constraints to matrix, solution, and rhs for (auto& bc : constraints_) bc.apply(*systemMatrix_, *solution_, *rhs_); info(2," assemble boundary conditions needed {} seconds", t2.elapsed()); msg("fill-in of assembled matrix: {}", systemMatrix_->nnz()); msg("assemble needed {} seconds", t.elapsed()); } template void ProblemStat:: writeFiles(AdaptInfo& adaptInfo, bool force) { Dune::Timer t; for (auto writer : filewriter_) writer->write(adaptInfo, force); msg("writeFiles needed {} seconds", t.elapsed()); } } // end namespace AMDiS