[test][stokes] Add simple upscaling application

Co-authored-by: Martin Schneider <martin.schneider@iws.uni-stuttgart.de>

test/freeflow/navierstokes/CMakeLists.txt

@@ -6,6 +6,7 @@ add_subdirectory(channel)
 add_subdirectory(donea)
 add_subdirectory(kovasznay)
 add_subdirectory(periodic)
+add_subdirectory(permeabilityupscaling)
 add_subdirectory(sincos)
 add_subdirectory(unstructured)
 

test/freeflow/navierstokes/permeabilityupscaling/CMakeLists.txt

+# SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+dune_symlink_to_source_files(FILES "params.input" "sphere.raw")
+dumux_add_test(NAME test_ff_stokes_permeability_upscaling
+               SOURCES main.cc
+               LABELS freeflow navierstokes
+               CMAKE_GUARD "( HAVE_UMFPACK AND dune-subgrid_FOUND )")

test/freeflow/navierstokes/permeabilityupscaling/create_sphere_image.py

+#!/usr/bin/env python3
+# SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+import numpy as np
+
+
+def mask_sphere_3d(img, center, radius):
+    """create a mask in sphere form with center and radius"""
+    z, y, x = img.shape
+    X, Y, Z = np.ogrid[:x, :y, :z]
+    dist2 = (X - center[0])**2 + (Y-center[1])**2 + (Z-center[2])**2
+    return (dist2 <= radius**2).T
+
+
+img = np.zeros((15, 11, 13), dtype=np.uint8) # order is Z-Y-X
+mask = mask_sphere_3d(img, center=(6, 5, 7), radius=7)
+mask2 = mask_sphere_3d(img, center=(6, 5, 7), radius=3)
+img[~mask] = 1
+img[mask2] = 1 # non-connected sphere in the middle for testing
+img.flatten().tofile("sphere.raw")

test/freeflow/navierstokes/permeabilityupscaling/main.cc

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup NavierStokesTests
+ * \brief Pore flow test for permeability upscaling
+ */
+
+#include <config.h>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/common/float_cmp.hh>
+
+#include <dumux/common/initialize.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/indextraits.hh>
+#include <dumux/common/properties.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager_sub.hh>
+#include <dumux/io/grid/gridmanager_yasp.hh>
+
+#include <dumux/linear/istlsolvers.hh>
+#include <dumux/linear/linearsolvertraits.hh>
+#include <dumux/linear/linearalgebratraits.hh>
+
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/traits.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+
+#include <dumux/freeflow/navierstokes/velocityoutput.hh>
+#include <dumux/freeflow/navierstokes/fluxoveraxisalignedsurface.hh>
+
+#include "properties.hh"
+
+int main(int argc, char** argv)
+{
+    using namespace Dumux;
+
+    // define the type tag for this problem
+    using MomentumTypeTag = Properties::TTag::PoreFlowTestMomentum;
+    using MassTypeTag = Properties::TTag::PoreFlowTestMass;
+
+    // maybe initialize MPI and/or multithreading backend
+    initialize(argc, argv);
+    const auto& mpiHelper = Dune::MPIHelper::instance();
+
+    // print dumux start message
+    if (mpiHelper.rank() == 0)
+        DumuxMessage::print(/*firstCall=*/true);
+
+    // parse command line arguments and input file
+    Parameters::init(argc, argv);
+
+    // create a grid
+    using Scalar = GetPropType<MassTypeTag, Properties::Scalar>;
+    using Grid = GetPropType<MassTypeTag, Properties::Grid>;
+    Dumux::GridManager<Grid> gridManager;
+    gridManager.init();
+
+    // we compute on the leaf grid view
+    const auto& leafGridView = gridManager.grid().leafGridView();
+
+    // create the finite volume grid geometry
+    using MomentumGridGeometry = GetPropType<MomentumTypeTag, Properties::GridGeometry>;
+    auto momentumGridGeometry = std::make_shared<MomentumGridGeometry>(leafGridView);
+    using MassGridGeometry = GetPropType<MassTypeTag, Properties::GridGeometry>;
+    auto massGridGeometry = std::make_shared<MassGridGeometry>(leafGridView);
+
+    // mass-momentum coupling manager
+    using CouplingManager = GetPropType<MomentumTypeTag, Properties::CouplingManager>;
+    auto couplingManager = std::make_shared<CouplingManager>();
+
+    // the problems (boundary conditions)
+    using MomentumProblem = GetPropType<MomentumTypeTag, Properties::Problem>;
+    auto momentumProblem = std::make_shared<MomentumProblem>(momentumGridGeometry, couplingManager);
+    using MassProblem = GetPropType<MassTypeTag, Properties::Problem>;
+    auto massProblem = std::make_shared<MassProblem>(massGridGeometry, couplingManager);
+
+    // the solution vector
+    constexpr auto momentumIdx = CouplingManager::freeFlowMomentumIndex;
+    constexpr auto massIdx = CouplingManager::freeFlowMassIndex;
+    using Traits = MultiDomainTraits<MomentumTypeTag, MassTypeTag>;
+    using SolutionVector = typename Traits::SolutionVector;
+    SolutionVector x;
+    x[momentumIdx].resize(momentumGridGeometry->numDofs());
+    x[massIdx].resize(massGridGeometry->numDofs());
+    x = 0.0;
+
+    // the grid variables
+    using MomentumGridVariables = GetPropType<MomentumTypeTag, Properties::GridVariables>;
+    auto momentumGridVariables = std::make_shared<MomentumGridVariables>(momentumProblem, momentumGridGeometry);
+    using MassGridVariables = GetPropType<MassTypeTag, Properties::GridVariables>;
+    auto massGridVariables = std::make_shared<MassGridVariables>(massProblem, massGridGeometry);
+
+    // compute coupling stencil and afterwards initialize grid variables (need coupling information)
+    couplingManager->init(momentumProblem, massProblem, std::make_tuple(momentumGridVariables, massGridVariables), x);
+    massGridVariables->init(x[massIdx]);
+    momentumGridVariables->init(x[momentumIdx]);
+
+    using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+    auto assembler = std::make_shared<Assembler>(std::make_tuple(momentumProblem, massProblem),
+                                                 std::make_tuple(momentumGridGeometry, massGridGeometry),
+                                                 std::make_tuple(momentumGridVariables, massGridVariables),
+                                                 couplingManager);
+
+    // initialize the vtk output module
+    using IOFields = GetPropType<MassTypeTag, Properties::IOFields>;
+    VtkOutputModule vtkWriter(*massGridVariables, x[massIdx], massProblem->name());
+    IOFields::initOutputModule(vtkWriter); // Add model specific output fields
+    vtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<MassGridVariables>>());
+    vtkWriter.write(0.0);
+
+    // the linear solver
+    using LinearSolver = UMFPackIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
+    auto linearSolver = std::make_shared<LinearSolver>();
+
+    // solve the non-linear system
+    using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+    NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+    nonLinearSolver.solve(x);
+
+    // post-processing and output
+    vtkWriter.write(1.0);
+
+    // set up two planes over which fluxes are calculated
+    FluxOverAxisAlignedSurface flux(*massGridVariables, x[massIdx], assembler->localResidual(massIdx));
+
+    using GridView = typename GetPropType<MassTypeTag, Properties::GridGeometry>::GridView;
+    using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
+
+    const Scalar xMin = massGridGeometry->bBoxMin()[0];
+    const Scalar xMax = massGridGeometry->bBoxMax()[0];
+    const Scalar yMin = massGridGeometry->bBoxMin()[1];
+    const Scalar yMax = massGridGeometry->bBoxMax()[1];
+    const Scalar zMin = massGridGeometry->bBoxMin()[2];
+    const Scalar zMax = massGridGeometry->bBoxMax()[2];
+
+    // the first plane is at the inlet
+    const auto inletLowerLeft = GlobalPosition{xMin, yMin, zMin};
+    const auto inletUpperRight = GlobalPosition{xMin, yMax, zMax};
+    flux.addAxisAlignedSurface("inlet", inletLowerLeft, inletUpperRight);
+
+    // The last plane is placed at the outlet of the channel.
+    const auto outletLowerLeft = GlobalPosition{xMax, yMin, zMin};
+    const auto outletUpperRight = GlobalPosition{xMax, yMax, zMax};
+    flux.addAxisAlignedSurface("outlet", outletLowerLeft, outletUpperRight);
+
+    // calculate and print mass fluxes over the planes
+    flux.calculateAllFluxes();
+    std::cout << "mass flux at inlet is: " << flux.flux("inlet") << " µm^3/s" <<  std::endl;
+    std::cout << "mass flux at outlet is: " << flux.flux("outlet") << " µm^3/s" << std::endl;
+
+    // make sure solver is converged enough
+    if (Dune::FloatCmp::ne(flux.flux("inlet"), -flux.flux("outlet"), 1e-5))
+        DUNE_THROW(Dune::Exception, "Inlet and outlet fluxes do not match by " << flux.flux("inlet")+flux.flux("outlet"));
+
+    const auto cells = getParam<std::array<double, Grid::dimension>>("Grid.Cells");
+    const auto dims = getParam<std::array<double, Grid::dimension>>("Grid.PixelDimensions");
+    const auto totalAreaOutlet = dims[1]*cells[1]*dims[2]*cells[2];
+    // viscosity, density, and pressure gradient are 1.0
+    // With the pressure difference across the domain defined as 1/m,
+    // the permeability can be determined by dividing the flux at the outlet by the area of the outlet.
+    const auto K = std::abs(flux.flux("outlet")/totalAreaOutlet)*1e-12; // domain is in µm
+    // The porosity can be calculated for uniform sized grid cells by dividing
+    // the number of cells in the leafGridView by the number of cells used in the original hostgrid.
+    const auto phi = double(leafGridView.size(0))/double(cells[0]*cells[1]*cells[2]);
+    std::cout << "Permeability is: " << K << " m^2" << std::endl;
+    std::cout << "Porosity is: " << phi << std::endl;
+
+    // test against reference solution (hard-coded here)
+    // note that for the test the geometry is very coarse to minimize runtime
+    // therefore these reference values are just regression test references
+    // computed with this test at the time it was set up
+    const auto KRef = 3.43761e-13;
+    const auto phiRef = 0.439627;
+    if (Dune::FloatCmp::ne(K, KRef, 1e-5))
+        DUNE_THROW(Dune::Exception, "Permeability " << K << " doesn't match reference " << KRef << " by " << K-KRef);
+    if (Dune::FloatCmp::ne(phi, phiRef, 1e-5))
+        DUNE_THROW(Dune::Exception, "Porosity " << phi << " doesn't match reference " << phiRef << " by " << phi-phiRef);
+
+    ////////////////////////////////////////////////////////////
+    // finalize, print dumux message to say goodbye
+    ////////////////////////////////////////////////////////////
+
+    // print dumux end message
+    if (mpiHelper.rank() == 0)
+    {
+        Parameters::print();
+        DumuxMessage::print(/*firstCall=*/false);
+    }
+
+    return 0;
+}

test/freeflow/navierstokes/permeabilityupscaling/params.input

+[Grid]
+BinaryMask = sphere.raw
+Cells = 13 11 15
+PixelDimensions = 2 2 2
+
+[Problem]
+Name = test_stokes_upscaling
+EnableGravity = false
+EnableInertiaTerms = false
+
+[Component]
+# this test requires a unit fluid
+LiquidDensity = 1
+LiquidKinematicViscosity = 1
+
+[ Newton ]
+MaxSteps = 20
+MaxRelativeShift = 1e-8
+
+[Vtk]
+WriteFaceData = false
+
+[FluxOverPlane]
+Verbose = true

test/freeflow/navierstokes/permeabilityupscaling/problem.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup NavierStokesTests
+ * \brief Pore flow test for the staggered grid (Navier-)Stokes model.
+ */
+
+#ifndef DUMUX_TEST_STOKES_PERMEABILITY_UPSCALING_PROBLEM_HH
+#define DUMUX_TEST_STOKES_PERMEABILITY_UPSCALING_PROBLEM_HH
+
+#include <dune/common/float_cmp.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+
+#include <dumux/freeflow/navierstokes/momentum/fluxhelper.hh>
+#include <dumux/freeflow/navierstokes/scalarfluxhelper.hh>
+#include <dumux/freeflow/navierstokes/mass/1p/advectiveflux.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup NavierStokesTests
+ * \brief  Test problem of the pore flow test for permeability upscaling
+ */
+template <class TypeTag, class BaseProblem>
+class PoreFlowTestProblem : public BaseProblem
+{
+    using ParentType = BaseProblem;
+
+    using BoundaryTypes = typename ParentType::BoundaryTypes;
+    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+    using FVElementGeometry = typename GridGeometry::LocalView;
+    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+    using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+    using InitialValues = typename ParentType::InitialValues;
+    using Sources = typename ParentType::Sources;
+    using DirichletValues = typename ParentType::DirichletValues;
+    using BoundaryFluxes = typename ParentType::BoundaryFluxes;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+
+    static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
+    using Element = typename FVElementGeometry::Element;
+    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+    using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+
+public:
+    PoreFlowTestProblem(std::shared_ptr<const GridGeometry> gridGeometry,
+                        std::shared_ptr<CouplingManager> couplingManager)
+    : ParentType(gridGeometry, couplingManager)
+    {
+        // gradP is 1/m
+        deltaP_ = (this->gridGeometry().bBoxMax()[0]-this->gridGeometry().bBoxMin()[0]);
+    }
+
+    /*!
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary control volume.
+     *
+     * \param globalPos The position of the center of the finite volume
+     */
+    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+    {
+        BoundaryTypes values;
+
+        if constexpr (ParentType::isMomentumProblem())
+        {
+            values.setAllDirichlet();
+            if (isOutlet_(globalPos) || isInlet_(globalPos))
+                values.setAllNeumann();
+        }
+        else
+            values.setNeumann(Indices::conti0EqIdx);
+
+        return values;
+    }
+
+    /*!
+     * \brief Evaluates the boundary conditions for a Dirichlet control volume.
+     * \param globalPos The center of the finite volume which ought to be set.
+     */
+    DirichletValues dirichletAtPos(const GlobalPosition &globalPos) const
+    {
+        // no-flow/no-slip
+        return DirichletValues(0.0);
+    }
+
+    /*!
+     * \brief Evaluates the boundary conditions for a Neumann control volume
+     */
+    template<class ElementVolumeVariables, class ElementFluxVariablesCache>
+    BoundaryFluxes neumann(const Element& element,
+                           const FVElementGeometry& fvGeometry,
+                           const ElementVolumeVariables& elemVolVars,
+                           const ElementFluxVariablesCache& elemFluxVarsCache,
+                           const SubControlVolumeFace& scvf) const
+    {
+        BoundaryFluxes values(0.0);
+        const auto& globalPos = scvf.ipGlobal();
+
+        if constexpr (ParentType::isMomentumProblem())
+        {
+            const auto p = isInlet_(globalPos) ? deltaP_ : 0.0;
+            using FluxHelper = NavierStokesMomentumBoundaryFlux<typename GridGeometry::DiscretizationMethod>;
+            values = FluxHelper::fixedPressureMomentumFlux(
+                *this, fvGeometry, scvf, elemVolVars, elemFluxVarsCache, p
+            );
+        }
+        else
+        {
+            values = NavierStokesScalarBoundaryFluxHelper<AdvectiveFlux<ModelTraits>>::scalarOutflowFlux(
+                *this, element, fvGeometry, scvf, elemVolVars
+            );
+        }
+
+        return values;
+    }
+
+private:
+
+    bool isInlet_(const GlobalPosition& globalPos) const
+    { return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
+
+    bool isOutlet_(const GlobalPosition& globalPos) const
+    { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
+
+    static constexpr Scalar eps_ = 1e-8;
+    Scalar deltaP_;
+};
+
+} // end namespace Dumux
+
+#endif

test/freeflow/navierstokes/permeabilityupscaling/properties.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup NavierStokesTests
+ * \brief The properties of the pore flow test for permeability upscaling
+ */
+#ifndef DUMUX_UPSCALING_TEST_PROPERTIES_HH
+#define DUMUX_UPSCALING_TEST_PROPERTIES_HH
+
+#ifndef ENABLECACHING
+#define ENABLECACHING 1
+#endif
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/subgrid/subgrid.hh>
+
+#include <dumux/discretization/fcstaggered.hh>
+#include <dumux/discretization/cctpfa.hh>
+
+#include <dumux/freeflow/navierstokes/momentum/model.hh>
+#include <dumux/freeflow/navierstokes/mass/1p/model.hh>
+#include <dumux/freeflow/navierstokes/momentum/problem.hh>
+#include <dumux/freeflow/navierstokes/mass/problem.hh>
+#include <dumux/multidomain/traits.hh>
+#include <dumux/multidomain/freeflow/couplingmanager.hh>
+#include <dumux/material/components/constant.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+
+#include "problem.hh"
+namespace Dumux::Properties {
+
+// Create new type tags
+namespace TTag {
+struct PoreFlowTest {};
+struct PoreFlowTestMomentum { using InheritsFrom = std::tuple<PoreFlowTest, NavierStokesMomentum, FaceCenteredStaggeredModel>; };
+struct PoreFlowTestMass { using InheritsFrom = std::tuple<PoreFlowTest, NavierStokesMassOneP, CCTpfaModel>; };
+} // end namespace TTag
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::PoreFlowTestMomentum>
+{ using type = PoreFlowTestProblem<TypeTag, Dumux::NavierStokesMomentumProblem<TypeTag>>; };
+
+template<class TypeTag>
+struct Problem<TypeTag, TTag::PoreFlowTestMass>
+{ using type = PoreFlowTestProblem<TypeTag, Dumux::NavierStokesMassProblem<TypeTag>>; };
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::PoreFlowTest>
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using type = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::PoreFlowTest>
+{
+    static constexpr int dim = 3;
+    using HostGrid = Dune::YaspGrid<dim, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, dim> >;
+    using type = Dune::SubGrid<HostGrid::dimension, HostGrid>;
+};
+
+template<class TypeTag>
+struct EnableGridGeometryCache<TypeTag, TTag::PoreFlowTest> { static constexpr bool value = ENABLECACHING; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::PoreFlowTest> { static constexpr bool value = ENABLECACHING; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::PoreFlowTest> { static constexpr bool value = ENABLECACHING; };
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::PoreFlowTest>
+{
+    using Traits = MultiDomainTraits<TTag::PoreFlowTestMomentum, TTag::PoreFlowTestMass>;
+    using type = FreeFlowCouplingManager<Traits>;
+};
+
+} // end namespace Dumux::Properties
+
+#endif