poisson-periodic.cc

Commented version of this example This example on gitlab

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#include <config.h>
 
#include <vector>
#include <cmath>
#include <numbers>
#include <thread>
 
#include <dune/common/exceptions.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/stringutility.hh>
#include <dune/common/rangeutilities.hh>
 
#include <dune/geometry/type.hh>
 
#include <dune/grid/uggrid.hh>
 
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
 
#include <dune/functions/backends/istlvectorbackend.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
 
#include <dune/vtk/vtkwriter.hh>
#include <dune/vtk/datacollectors/lagrangedatacollector.hh>
 
#include <dune/assembler/backends/istlvectorbackend.hh>
#include <dune/assembler/backends/istlmatrixbackend.hh>
#include <dune/assembler/defaultglobalassembler.hh>
 
#include <dune/assembler/parallel/gridcoloring.hh>
#include <dune/assembler/parallel/coloredrangeexecutor.hh>
 
#include <dune/fufem/logger.hh>
#include <dune/fufem/parallel/elementcoloring.hh>
#include <dune/fufem/forms/forms.hh>
#include <dune/fufem/grid/boundarydata.hh>
#include <dune/fufem/domains/boundary.hh>
#include <dune/fufem/domains/boundarypatch.hh>
#include <dune/fufem/constraints/affineconstraints.hh>
#include <dune/fufem/constraints/boundaryconstraints.hh>
#include <dune/fufem/constraints/periodicconstraints.hh>
 
#include "utilities.hh"
 
int main (int argc, char *argv[]) try
{
 
  auto log = Dune::Fufem::makeLogger(std::cout, "[% 8.3f] [% 8.3f] ");
 
  auto parameters = Dune::ParameterTree();
  Dune::ParameterTreeParser::readOptions(argc, argv, parameters);
 
  auto refinements  = parameters.get<std::size_t>("refinements", 3);
  auto threadCount  = parameters.get<std::size_t>("threads", std::thread::hardware_concurrency());
  auto tol          = parameters.get<double>("tol", 1e-4);
  auto maxIter      = parameters.get<std::size_t>("maxIter", 1000);
  auto verbosity    = parameters.get<std::size_t>("verbosity", 2);
 
  log("Command line processed");
 
  Dune::MPIHelper::instance(argc, argv);
 
  log("MPI initialized");
  auto dirichletIndicator = [&](auto x) {
    return (x[1] <= 0.25) or (x[1] > 0.75);
  };
 
  auto dirichletValues = [](const auto& x){ return 0; };
 
  auto rightHandSide = [] (const auto& x) { return 20*(x[0]>=.5);};
  using Grid = Dune::UGGrid<2>;
 
  auto factory = Dune::GridFactory<Grid>();
  for(unsigned int k : Dune::range(9))
    factory.insertVertex({0.5*(k%3), 0.5*(k/3)});
  factory.insertElement(Dune::GeometryTypes::cube(2), {0, 1, 3, 4});
  factory.insertElement(Dune::GeometryTypes::cube(2), {1, 2, 4, 5});
  factory.insertElement(Dune::GeometryTypes::simplex(2), {3, 4, 6});
  factory.insertElement(Dune::GeometryTypes::simplex(2), {4, 7, 6});
  factory.insertElement(Dune::GeometryTypes::simplex(2), {4, 5, 7});
  factory.insertElement(Dune::GeometryTypes::simplex(2), {5, 8, 7});
  factory.insertBoundarySegment({0, 3});
  factory.insertBoundarySegment({3, 6});
  factory.insertBoundarySegment({2, 5});
  factory.insertBoundarySegment({5, 8});
  constexpr int leftID = 0;
  constexpr int rightID = 1;
  auto boundaryIDs = std::vector({leftID, leftID, rightID, rightID});
  auto periodicRightToLeft = [](auto x) {
    x[0] -= 1.0;
    return x;
  };
  auto gridPtr = std::unique_ptr(factory.createGrid());
  auto& grid = *gridPtr;
 
  log("Grid created");
  auto boundaryData = Dune::Fufem::BoundaryData(boundaryIDs, factory, grid);
  grid.globalRefine(refinements);
 
  log(Dune::formatString("Grid refined %d times", refinements));
 
  auto gridView = grid.leafGridView();
  auto dirichletPatch = Dune::Fufem::BoundaryPatch(gridView);
  auto periodicLeftPatch = Dune::Fufem::BoundaryPatch(gridView);
  auto periodicRightPatch = Dune::Fufem::BoundaryPatch(gridView);
  for (const auto& intersection : Dune::Fufem::Boundary(gridView))
  {
    if (dirichletIndicator(intersection.geometry().center()))
      dirichletPatch.insertFace(intersection);
    if (boundaryData(intersection) == leftID)
      periodicLeftPatch.insertFace(intersection);
    if (boundaryData(intersection) == rightID)
      periodicRightPatch.insertFace(intersection);
  }
 
  log("Boundary patches set up");
  auto coloredRange = Dune::Assembler::Experimental::coloredElementRange(gridView);
 
  log("Grid coloring computed");
 
  auto executor = Dune::Assembler::Experimental::ColoredRangeExecutor(coloredRange, threadCount);
 
  log(Dune::formatString("Using parallelism with %d threads", threadCount));
  using namespace Dune::Functions::BasisFactory;
  auto basis = makeBasis(gridView, lagrange<2>());
 
  log(Dune::formatString("Basis created with dimension %d", basis.dimension()));
  using Vector = Dune::BlockVector<double>;
  using Matrix = Dune::BCRSMatrix<double>;
 
  auto sol = Vector();
  auto rhs = Vector();
  auto stiffnessMatrix = Matrix();
  auto constraints = Dune::Fufem::makeAffineConstraints(basis);
 
  Dune::Fufem::computeBoundaryConstraints(constraints, basis, dirichletValues, dirichletPatch);
 
  log("Boundary constraints computed");
 
  Dune::Fufem::computePeriodicConstraints(constraints, basis, periodicLeftPatch, periodicRightPatch, periodicRightToLeft);
 
  log("Periodic constraints computed");
  using namespace Dune::Fufem::Forms;
 
  auto u = trialFunction(basis);
  auto v = testFunction(basis);
  auto f = Coefficient(Dune::Functions::makeGridViewFunction(rightHandSide, gridView));
 
  auto a = integrate(dot(grad(u), grad(v)));
  auto l = integrate(f*v);
  auto matrixBackend = Dune::Assembler::ISTLMatrixBackend(stiffnessMatrix);
  auto rhsBackend = Dune::Assembler::ISTLVectorBackend(rhs);
  auto assembler = Dune::Assembler::Assembler(basis, basis, executor);
 
  log("Assembler set up");
  auto patternBuilder = matrixBackend.patternBuilder();
 
  patternBuilder.resize(basis, basis);
 
  assembler.assembleMatrixPattern(a, patternBuilder);
 
  constraints.extendMatrixPattern(patternBuilder, basis);
 
  log("Matrix pattern assembled");
  patternBuilder.setupMatrix();
 
  log("Matrix set up");
 
  assembler.assembleMatrixEntries(a, matrixBackend);
 
  log("Matrix assembled");
  rhsBackend.resize(basis);
  rhsBackend.setZero();
 
  assembler.assembleVectorEntries(l, rhsBackend);
 
  log("RHS assembled");
  constraints.constrainLinearSystem(stiffnessMatrix, rhs);
 
  log("Boundary condition incorporated into system");
  Dune::Assembler::ISTLVectorBackend(sol).resize(basis);
  Dune::Assembler::ISTLVectorBackend(sol).setZero();
  auto matrixOperator = Dune::MatrixAdapter<Matrix,Vector,Vector>(stiffnessMatrix);
 
  auto preconditioner = makeAMGPreconditioner<Vector>(stiffnessMatrix);
 
  log("Preconditioner created");
 
  auto solver = Dune::CGSolver<Vector>(
                  matrixOperator,
                  preconditioner,
                  tol,
                  maxIter,
                  verbosity
                );
 
  log("Solver created");
 
  auto solverStatistics = Dune::InverseOperatorResult();
 
  solver.apply(sol, rhs, solverStatistics);
 
  log("Linear system solved");
  constraints.interpolate(sol);
 
  log("Constrained solution interpolated");
  auto vtkWriter = Dune::Vtk::UnstructuredGridWriter(Dune::Vtk::LagrangeDataCollector(basis.gridView(), 2));
  vtkWriter.addPointData(u|sol, Dune::VTK::FieldInfo("sol", Dune::VTK::FieldInfo::Type::scalar, 1));
  vtkWriter.write("poisson-periodic");
 
  log("Solution written to vtk file");
}
catch (Dune::Exception& e)
{
  std::cout << e.what() << std::endl;
  return 1;
}