stokes-darcy.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 <cstddef>
#include <cmath>
#include <iostream>
#include <string>
 
#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/indices.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/stringutility.hh>
 
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/common/gridfactory.hh>
 
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/umfpack.hh>
 
#include <dune/functions/common/subdomain.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/restrictedbasis.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
 
#include <dune/assembler/backends/istlmatrixbackend.hh>
#include <dune/assembler/backends/istlvectorbackend.hh>
#include <dune/assembler/defaultglobalassembler.hh>
 
#include <dune/vtk/datacollectors/discontinuouslagrangedatacollector.hh>
#include <dune/vtk/datacollectors/lagrangedatacollector.hh>
#include <dune/vtk/vtkwriter.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/forms/forms.hh>
#include <dune/fufem/grid/boundarydata.hh>
#include <dune/fufem/grid/elementdata.hh>
#include <dune/fufem/logger.hh>
 
int main (int argc, char *argv[]) try
{
 
  auto log = Dune::Fufem::makeLogger(std::cout, "[% 8.3f] [% 8.3f] ");
 
  Dune::MPIHelper::instance(argc, argv);
 
  log("MPI initialized");
  auto parameters = Dune::ParameterTree();
  Dune::ParameterTreeParser::readOptions(argc, argv, parameters);
 
  auto nu           = parameters.get<double>("viscosity", 1.0);
  auto K            = parameters.get<double>("permeability", 1.0e-2);
  auto alpha_BJS    = parameters.get<double>("alpha_BJS", 1.0);
  auto g            = parameters.get<double>("gravity", 1.0);
  auto meshFile     = parameters.get<std::string>("meshFile", "stokes-darcy.msh");
  auto StokesID     = parameters.get<int>("StokesID", 1);
  auto DarcyID      = parameters.get<int>("DarcyID", 2);
  auto refinements  = parameters.get<std::size_t>("refinements", 2);
 
  log("Command line processed");
  constexpr std::size_t dim = 2;
 
  auto stokesVelocityDirichletValues = [](auto x)
  {
    return Dune::FieldVector<double,dim>{1.0, 0.0};
  };
 
  auto darcyDirichletValues = [](auto x)
  {
    return 0;
  };
  using Grid = Dune::UGGrid<dim>;
 
  auto factory = Dune::GridFactory<Grid>();
  auto reader = Dune::GmshReader<Grid>(meshFile, factory);
  auto gridPtr = factory.createGrid();
  auto& grid = *gridPtr;
  auto elementData = Dune::Fufem::ElementData(reader, factory, grid);
  auto boundaryData = Dune::Fufem::BoundaryData(reader, factory, grid);
  grid.globalRefine(refinements);
 
  auto gridView = grid.leafGridView();
 
  elementData.update(gridView);
 
  log("Grid refined");
  using Dune::Functions::Experimental::SubDomain;
  using Dune::Functions::Experimental::SubDomainInterface;
 
  auto stokesDomain = SubDomain(gridView);
  for(auto&& element : elements(gridView))
    if (elementData(element) == StokesID)
      stokesDomain.insertElement(element);
 
  log(Dune::formatString("Stokes domain contains %d elements", stokesDomain.indexSet().size(0)));
 
  auto darcyDomain = SubDomain(gridView);
  for(auto&& element : elements(gridView))
    if (elementData(element) == DarcyID)
      darcyDomain.insertElement(element);
 
  log(Dune::formatString("Darcy domain contains %d elements", darcyDomain.indexSet().size(0)));
  auto interface = SubDomainInterface(stokesDomain, darcyDomain);
 
  log("Interface created");
  auto stokesDirichletBoundary = Dune::Fufem::BoundaryPatch(gridView);
  for (const auto& intersection : Dune::Fufem::Boundary(gridView))
    if (boundaryData(intersection) == StokesID)
      stokesDirichletBoundary.insertFace(intersection);
 
  auto darcyDirichletBoundary = Dune::Fufem::BoundaryPatch(gridView);
  for (const auto& intersection : Dune::Fufem::Boundary(gridView))
    if (boundaryData(intersection) == DarcyID)
      darcyDirichletBoundary.insertFace(intersection);
 
  log("Dirichlet boundaries created");
  using namespace Dune::Functions::BasisFactory;
  using namespace Dune::Functions::Experimental::BasisFactory;
 
  constexpr std::size_t maxOrder = 2;
  auto basis = makeBasis(gridView,
      composite(
        restrict(
          composite(
            power<dim>(lagrange<maxOrder>(), flatInterleaved()),
            lagrange<maxOrder-1>(),
            flatLexicographic()
          ), stokesDomain),
        restrict(lagrange<maxOrder+1>(), darcyDomain),
        flatLexicographic()
      )
    );
 
  log(Dune::formatString("Basis created with dimension %d", basis.dimension()));
  using namespace Dune::Indices;
 
  auto stokesVelocityBasis = Dune::Functions::subspaceBasis(basis, _0, _0);
  auto stokesPressureBasis = Dune::Functions::subspaceBasis(basis, _0, _1);
  auto darcyBasis = Dune::Functions::subspaceBasis(basis, _1);
  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);
  computeBoundaryConstraints(constraints, stokesVelocityBasis, stokesVelocityDirichletValues, stokesDirichletBoundary);
  computeBoundaryConstraints(constraints, darcyBasis, darcyDirichletValues, darcyDirichletBoundary);
 
  log("Boundary constraints computed");
  using namespace Dune::Fufem::Forms;
 
  auto u    = trialFunction(stokesVelocityBasis);
  auto p    = trialFunction(stokesPressureBasis);
  auto phi  = trialFunction(darcyBasis);
 
  auto v    = testFunction(stokesVelocityBasis);
  auto q    = testFunction(stokesPressureBasis);
  auto psi  = testFunction(darcyBasis);
  auto E = [&](const auto& v) {
    return symmetrize(grad(v));
  };
 
  auto normal = faceNormal(gridView);
  auto P_tangent = compose([&](auto n) {
    Dune::FieldMatrix<double, dim, dim> T;
    for(auto i : Dune::range(dim))
      for(auto j : Dune::range(dim))
        T[i][j] = (i==j) - n[i]*n[j];
    return T;
  }, normal);
  double alpha = alpha_BJS/std::sqrt(K);
 
  auto A =  integrate( 2*nu*dot(E(u), E(v)) - div(u)*q - div(v)*p, stokesDomain )
          + integrate( dot(K*grad(phi), grad(psi)), darcyDomain )
          + integrate( alpha*dot(P_tangent*u, v), interface )
          + integrate( g*outside(phi)*dot(v, normal), interface )
          + integrate( - outside(psi)*dot(u, normal), interface );
  auto assembler = Dune::Assembler::Assembler(basis, basis);
 
  log("Assembler set up");
  auto matrixBackend = Dune::Assembler::ISTLMatrixBackend(stiffnessMatrix);
  auto rhsBackend = Dune::Assembler::ISTLVectorBackend(rhs);
  auto patternBuilder = matrixBackend.patternBuilder();
  patternBuilder.resize(basis, basis);
  assembler.assembleMatrixPattern(A, patternBuilder);
 
  log("Matrix pattern assembled");
 
  patternBuilder.setupMatrix();
 
  log("Matrix set up");
  assembler.assembleMatrixEntries(A, matrixBackend);
 
  log("Matrix assembled");
  rhsBackend.resize(basis);
  rhsBackend.setZero();
 
  log("RHS assembled");
  constraints.constrainLinearSystem(stiffnessMatrix, rhs);
 
  log("Boundary condition incorporated into system");
  auto solBackend = Dune::Assembler::ISTLVectorBackend(sol);
  solBackend.resize(basis);
  solBackend.setZero();
#if HAVE_UMFPACK
  auto solver = Dune::UMFPack<Matrix>();
  solver.setMatrix(stiffnessMatrix);
#else
  auto matrixOperator = Dune::MatrixAdapter<Matrix,Vector,Vector>(stiffnessMatrix);
  auto preconditioner = Dune::SeqILU<Matrix,Vector,Vector>(stiffnessMatrix,1.0);
  auto solver = Dune::RestartedGMResSolver<Vector>(
                  matrixOperator,
                  preconditioner,
                  1e-6,
                  1000,
                  1000,
                  2);
#endif
 
  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::DiscontinuousLagrangeDataCollector(basis.gridView(), maxOrder));
  vtkWriter.addPointData(u|sol, Dune::VTK::FieldInfo("stokes_velocity", Dune::VTK::FieldInfo::Type::vector, dim));
  vtkWriter.addPointData(p|sol, Dune::VTK::FieldInfo("stokes_pressure", Dune::VTK::FieldInfo::Type::scalar, 1));
  vtkWriter.addPointData(phi|sol, Dune::VTK::FieldInfo("darcy_pressure", Dune::VTK::FieldInfo::Type::scalar, 1));
  vtkWriter.addPointData(-K*jacobian(phi)|sol, Dune::VTK::FieldInfo("darcy_velocity", Dune::VTK::FieldInfo::Type::vector, dim));
  vtkWriter.addPointData((p+phi)|sol, Dune::VTK::FieldInfo("pressure", Dune::VTK::FieldInfo::Type::scalar, 1));
  vtkWriter.addPointData((u-K*jacobian(phi))|sol, Dune::VTK::FieldInfo("velocity", Dune::VTK::FieldInfo::Type::vector, dim));
  vtkWriter.write("stokes-darcy");
 
  log("Solution written to vtk file");
}
catch (Dune::Exception& e)
{
  std::cout << e.what() << std::endl;
  return 1;
}