[embedded] Prototype for removed singularity

test/multidomain/embedded/1d3d/1p_1p/CMakeLists.txt

-dune_add_test(SOURCES test_1p_1p.cc)
-dune_symlink_to_source_files(FILES "grids" "test_1p_1p.input")
-set(CMAKE_BUILD_TYPE Release)
+add_subdirectory("analytical")
+add_subdirectory("analytical2")

test/multidomain/embedded/1d3d/1p_1p/analytical/CMakeLists.txt

+dune_add_test(SOURCES test_1p_1p.cc)
+dune_symlink_to_source_files(FILES "grids" "test_1p_1p.input")
+set(CMAKE_BUILD_TYPE Release)

test/multidomain/embedded/1d3d/1p_1p/analytical2/CMakeLists.txt

+dune_add_test(NAME test_1d3d_analytical2
+              SOURCES test_1p_1p.cc)
+dune_symlink_to_source_files(FILES "grids" "test_1p_1p.input")
+set(CMAKE_BUILD_TYPE Release)

test/multidomain/embedded/1d3d/1p_1p/analytical2/bloodflowproblem.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief A test problem for the one-phase blood flow model:
+ * Blood is flowing through a 1d network grid.
+ */
+#ifndef DUMUX_BLOOD_FLOW_PROBLEM_HH
+#define DUMUX_BLOOD_FLOW_PROBLEM_HH
+
+#include <dumux/common/parameters.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+
+#include <dumux/porousmediumflow/1p/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/porousmediumflow/1p/incompressiblelocalresidual.hh>
+
+#include <dumux/material/components/constant.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+
+#include "bloodflowspatialparams.hh"
+
+namespace Dumux {
+// forward declaration
+template <class TypeTag> class BloodFlowProblem;
+
+namespace Properties {
+
+NEW_TYPE_TAG(BloodFlowTypeTag, INHERITS_FROM(OneP));
+NEW_TYPE_TAG(BloodFlowCCTypeTag, INHERITS_FROM(CCTpfaModel, BloodFlowTypeTag));
+
+// Set the grid type
+SET_TYPE_PROP(BloodFlowTypeTag, Grid, Dune::FoamGrid<1, 3>);
+
+SET_BOOL_PROP(BloodFlowTypeTag, EnableFVGridGeometryCache, true);
+SET_BOOL_PROP(BloodFlowTypeTag, EnableGridVolumeVariablesCache, true);
+SET_BOOL_PROP(BloodFlowTypeTag, EnableGridFluxVariablesCache, true);
+SET_BOOL_PROP(BloodFlowTypeTag, SolutionDependentAdvection, false);
+SET_BOOL_PROP(BloodFlowTypeTag, SolutionDependentMolecularDiffusion, false);
+SET_BOOL_PROP(BloodFlowTypeTag, SolutionDependentHeatConduction, false);
+
+// Set the problem property
+SET_TYPE_PROP(BloodFlowTypeTag, Problem, BloodFlowProblem<TypeTag>);
+
+// the fluid system
+SET_PROP(BloodFlowTypeTag, FluidSystem)
+{
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using type = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
+};
+
+// Set the problem property
+SET_TYPE_PROP(BloodFlowTypeTag, LocalResidual, OnePIncompressibleLocalResidual<TypeTag>);
+
+// Set the spatial parameters
+SET_TYPE_PROP(BloodFlowTypeTag, SpatialParams,
+              BloodFlowSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
+                                     typename GET_PROP_TYPE(TypeTag, Scalar)>);
+} // end namespace Properties
+
+/*!
+ * \ingroup OnePTests
+ * \brief Exact solution 1D-3D
+ */
+template <class TypeTag>
+class BloodFlowProblem : public PorousMediumFlowProblem<TypeTag>
+{
+    using ParentType = PorousMediumFlowProblem<TypeTag>;
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using PointSource = typename GET_PROP_TYPE(TypeTag, PointSource);
+    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+    using GridView = typename FVGridGeometry::GridView;
+    using FVElementGeometry = typename FVGridGeometry::LocalView;
+    using SubControlVolume = typename FVGridGeometry::SubControlVolume;
+    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+    using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using GlobalPosition = typename FVGridGeometry::GlobalCoordinate;
+    using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+
+public:
+    BloodFlowProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+                     std::shared_ptr<CouplingManager> couplingManager)
+    : ParentType(fvGridGeometry)
+    , couplingManager_(couplingManager)
+    {
+        //read parameters from input file
+        name_ = getParam<std::string>("Problem.Name") + "_1d";
+        p_in_ = getParam<Scalar>("BoundaryConditions1D.PressureInput");
+        exactPressure_.resize(this->fvGridGeometry().numDofs());
+
+        for (const auto& element : elements(this->fvGridGeometry().gridView()))
+        {
+            auto fvGeometry = localView(this->fvGridGeometry());
+            fvGeometry.bindElement(element);
+
+            for (auto&& scv : scvs(fvGeometry))
+                exactPressure_[scv.dofIndex()] = exactSolution(scv.dofPosition());
+        }
+    }
+
+    /*!
+     * \brief Return how much the domain is extruded at a given sub-control volume.
+     *
+     * The extrusion factor here makes extrudes the 1d line to a circular tube with
+     * cross-section area pi*r^2.
+     */
+    template<class ElementSolution>
+    Scalar extrusionFactor(const Element &element,
+                           const SubControlVolume &scv,
+                           const ElementSolution& elemSol) const
+    {
+        const auto eIdx = this->fvGridGeometry().elementMapper().index(element);
+        const auto radius = this->spatialParams().radius(eIdx);
+        return M_PI*radius*radius;
+    }
+
+    /*!
+     * \name Problem parameters
+     */
+    // \{
+
+    /*!
+     * \brief The problem name.
+     *
+     * This is used as a prefix for files generated by the simulation.
+     */
+    const std::string& name() const
+    { return name_; }
+
+    /*!
+     * \brief Return the temperature within the domain in [K].
+     *
+     */
+    Scalar temperature() const
+    { return 273.15 + 37.0; } // Body temperature
+
+    // \}
+    /*!
+     * \name Boundary conditions
+     */
+    // \{
+
+   /*!
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary segment.
+     *
+     * \param globalPos The global position
+     */
+    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+    {
+        BoundaryTypes values;
+        values.setAllDirichlet();
+        return values;
+    }
+
+    /*!
+     * \brief Evaluate the boundary conditions for a dirichlet
+     *        control volume.
+     *
+     * \param globalPos The global position
+     *
+     * For this method, the \a values parameter stores primary variables.
+     */
+    PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+    {
+        return PrimaryVariables(p_in_);
+    }
+
+
+    // \}
+
+    /*!
+     * \name Volume terms
+     */
+    // \{
+
+     /*!
+     * \brief Applies a vector of point sources. The point sources
+     *        are possibly solution dependent.
+     *
+     * \param pointSources A vector of PointSource s that contain
+              source values for all phases and space positions.
+     *
+     * For this method, the \a values method of the point source
+     * has to return the absolute mass rate in kg/s. Positive values mean
+     * that mass is created, negative ones mean that it vanishes.
+     */
+    void addPointSources(std::vector<PointSource>& pointSources) const
+    { pointSources = this->couplingManager().lowDimPointSources(); }
+
+    /*!
+     * \brief Evaluate the point sources (added by addPointSources)
+     *        for all phases within a given sub-control-volume.
+     *
+     * This is the method for the case where the point source is
+     * solution dependent and requires some quantities that
+     * are specific to the fully-implicit method.
+     *
+     * \param pointSource A single point source
+     * \param element The finite element
+     * \param fvGeometry The finite-volume geometry
+     * \param elemVolVars All volume variables for the element
+     * \param scv The sub-control volume within the element
+     *
+     * For this method, the \a values() method of the point sources returns
+     * the absolute rate mass generated or annihilate in kg/s. Positive values mean
+     * that mass is created, negative ones mean that it vanishes.
+     */
+    template<class ElementVolumeVariables>
+    void pointSource(PointSource& source,
+                     const Element &element,
+                     const FVElementGeometry& fvGeometry,
+                     const ElementVolumeVariables& elemVolVars,
+                     const SubControlVolume &scv) const
+    {
+        source = 0.0;
+    }
+
+    /*!
+     * \brief Evaluate the initial value for a control volume.
+     *
+     * For this method, the \a priVars parameter stores primary
+     * variables.
+     */
+    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+    { return PrimaryVariables(p_in_); }
+
+    // \}
+
+    //! The exact pressure solution
+    Scalar exactSolution(const GlobalPosition &globalPos) const
+    { return 1 + globalPos[2]; }
+
+    //! Called after every time step
+    //! Output the total global exchange term
+    void computeSourceIntegral(const SolutionVector& sol, const GridVariables& gridVars)
+    {
+        PrimaryVariables source(0.0);
+        for (const auto& element : elements(this->fvGridGeometry().gridView()))
+        {
+            auto fvGeometry = localView(this->fvGridGeometry());
+            fvGeometry.bindElement(element);
+
+            auto elemVolVars = localView(gridVars.curGridVolVars());
+            elemVolVars.bindElement(element, fvGeometry, sol);
+
+            for (auto&& scv : scvs(fvGeometry))
+            {
+                auto pointSources = this->scvPointSources(element, fvGeometry, elemVolVars, scv);
+                pointSources *= scv.volume()*elemVolVars[scv].extrusionFactor();
+                source += pointSources;
+            }
+        }
+
+        std::cout << "Global integrated source (1D): " << source << '\n';
+    }
+
+    /*!
+     * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
+     */
+    template<class VtkOutputModule>
+    void addVtkOutputFields(VtkOutputModule& vtk) const
+    {
+        vtk.addField(exactPressure_, "exact pressure");
+    }
+
+    //! Set the coupling manager
+    void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+    { couplingManager_ = cm; }
+
+    //! Get the coupling manager
+    const CouplingManager& couplingManager() const
+    { return *couplingManager_; }
+
+private:
+    std::vector<Scalar> exactPressure_;
+
+    static constexpr Scalar eps_ = 1.5e-7;
+    std::string name_;
+
+    Scalar p_in_;
+
+    std::shared_ptr<CouplingManager> couplingManager_;
+};
+
+} //end namespace Dumux
+
+#endif

test/multidomain/embedded/1d3d/1p_1p/analytical2/bloodflowspatialparams.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup OneTests
+ * \brief The spatial parameters class blood flow problem
+ */
+#ifndef DUMUX_BlOOD_FLOW_SPATIALPARAMS_HH
+#define DUMUX_BlOOD_FLOW_SPATIALPARAMS_HH
+
+#include <dumux/material/spatialparams/fv1p.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup OneTests
+ * \brief Definition of the spatial parameters for the blood flow problem
+ */
+template<class FVGridGeometry, class Scalar>
+class BloodFlowSpatialParams
+: public FVSpatialParamsOneP<FVGridGeometry, Scalar, BloodFlowSpatialParams<FVGridGeometry, Scalar>>
+{
+    using ThisType = BloodFlowSpatialParams<FVGridGeometry, Scalar>;
+    using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, ThisType>;
+    using GlobalPosition = typename FVGridGeometry::GlobalCoordinate;
+
+public:
+    // export permeability type
+    using PermeabilityType = Scalar;
+
+    BloodFlowSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry)
+    {
+        radius_ = getParam<Scalar>("SpatialParams.Radius");
+        permeability_ = getParam<Scalar>("SpatialParams.PermeabilityVessel", 1);
+    }
+
+    /*!
+     * \brief Return the intrinsic permeability for the current sub-control volume in [m^2].
+     *
+     * \param ipGlobal The integration point
+     */
+    PermeabilityType permeabilityAtPos(const GlobalPosition& ipGlobal) const
+    {
+        return permeability_;
+    }
+
+    //! we evaluate the permeability directly at the scvf since we have an analytical expression for it
+    static constexpr bool evaluatePermeabilityAtScvfIP()
+    { return true; }
+
+    /*!
+     * \brief Return the radius of the circular pipe for the current sub-control volume in [m].
+     *
+     * \param the index of the element
+     */
+    Scalar radius(unsigned int eIdxGlobal) const
+    {
+        return radius_;
+    }
+
+    /*!
+     * \brief Returns the porosity \f$[-]\f$
+     *
+     * \param globalPos the scv center
+     */
+    Scalar porosityAtPos(const GlobalPosition& globalPos) const
+    {
+        return 1.0;
+    }
+
+private:
+    Scalar radius_;
+    Scalar permeability_;
+};
+
+} // end namespace Dumux
+
+#endif

test/multidomain/embedded/1d3d/1p_1p/analytical2/darcyslaw.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup CCTpfaDiscretization
+ * \brief Darcy's law for cell-centered finite volume schemes with two-point flux approximation
+ */
+#ifndef MY_DUMUX_DISCRETIZATION_CC_TPFA_DARCYS_LAW_HH
+#define MY_DUMUX_DISCRETIZATION_CC_TPFA_DARCYS_LAW_HH
+
+#include <dumux/common/math.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/properties.hh>
+
+#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh>
+#include <dumux/discretization/cellcentered/tpfa/darcyslaw.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup CCTpfaDiscretization
+ * \brief Darcy's law for cell-centered finite volume schemes with two-point flux approximation
+ * \note Darcy's law is speialized for network and surface grids (i.e. if grid dim < dimWorld)
+ * \tparam Scalar the scalar type for scalar physical quantities
+ * \tparam FVGridGeometry the grid geometry
+ * \tparam isNetwork whether we are computing on a network grid embedded in a higher world dimension
+ */
+template<class Scalar, class FVGridGeometry, bool isNetwork = false>
+class MyCCTpfaDarcysLaw;
+
+/*!
+ * \ingroup CCTpfaDiscretization
+ * \brief Specialization of the CCTpfaDarcysLaw grids where dim=dimWorld
+ */
+template<class ScalarType, class FVGridGeometry>
+class MyCCTpfaDarcysLaw<ScalarType, FVGridGeometry, /*isNetwork*/ false>
+{
+    using ThisType = CCTpfaDarcysLaw<ScalarType, FVGridGeometry, /*isNetwork*/ false>;
+    using FVElementGeometry = typename FVGridGeometry::LocalView;
+    using SubControlVolume = typename FVGridGeometry::SubControlVolume;
+    using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
+    using GridView = typename FVGridGeometry::GridView;
+    using Element = typename GridView::template Codim<0>::Entity;
+
+    static constexpr int dim = GridView::dimension;
+    static constexpr int dimWorld = GridView::dimensionworld;
+
+    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+  public:
+    //! state the scalar type of the law
+    using Scalar = ScalarType;
+
+    //! state the discretization method this implementation belongs to
+    static const DiscretizationMethod discMethod = DiscretizationMethod::cctpfa;
+
+    //! state the type for the corresponding cache
+    using Cache = TpfaDarcysLawCache<ThisType, FVGridGeometry>;
+
+    //! Compute the advective flux
+    template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache>
+    static Scalar flux(const Problem& problem,
+                       const Element& element,
+                       const FVElementGeometry& fvGeometry,
+                       const ElementVolumeVariables& elemVolVars,
+                       const SubControlVolumeFace& scvf,
+                       int phaseIdx,
+                       const ElementFluxVarsCache& elemFluxVarsCache)
+    {
+        static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
+
+        const auto& fluxVarsCache = elemFluxVarsCache[scvf];
+
+        // Get the inside and outside volume variables
+        const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+        const auto& insideVolVars = elemVolVars[insideScv];
+        const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
+
+        if (enableGravity)
+        {
+            // do averaging for the density over all neighboring elements
+            const auto rho = scvf.boundary() ? outsideVolVars.density(phaseIdx)
+                                             : (insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx))*0.5;
+
+            // Obtain inside and outside pressures
+            const auto pInside = insideVolVars.pressure(phaseIdx);
+            const auto pOutside = outsideVolVars.pressure(phaseIdx);
+
+            const auto& tij = fluxVarsCache.advectionTij();
+            const auto& g = problem.gravityAtPos(scvf.ipGlobal());
+
+            //! compute alpha := n^T*K*g
+            const auto alpha_inside = vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g)*insideVolVars.extrusionFactor();
+
+            Scalar flux = tij*(pInside - pOutside) + rho*scvf.area()*alpha_inside;
+
+            //! On interior faces we have to add K-weighted gravitational contributions
+            if (!scvf.boundary())
+            {
+                const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
+                const auto outsideK = outsideVolVars.permeability();
+                const auto outsideTi = computeTpfaTransmissibility(scvf, outsideScv, outsideK, outsideVolVars.extrusionFactor());
+                const auto alpha_outside = vtmv(scvf.unitOuterNormal(), outsideK, g)*outsideVolVars.extrusionFactor();
+
+                flux += rho*tij/outsideTi*(alpha_inside - alpha_outside);
+            }
+
+            return flux;
+        }
+        else
+        {
+            // Obtain inside and outside pressures
+            const auto pInside = insideVolVars.pressure(phaseIdx);
+            const auto pOutside = outsideVolVars.pressure(phaseIdx);
+
+            // return flux
+            auto flux = fluxVarsCache.advectionTij()*(pInside - pOutside);
+
+            if (std::abs(scvf.unitOuterNormal()*GlobalPosition({0.0, 0.0, 1.0})) < 1e-10)
+            {
+                static const Scalar k = getParam<Scalar>("SpatialParams.PermeabilityTissue");
+
+
+                // additional fluxes due to point sources
+                auto key = std::make_pair(problem.fvGridGeometry().elementMapper().index(element), 0);
+                auto keyOutside = std::make_pair(problem.fvGridGeometry().elementMapper().index(problem.fvGridGeometry().element(scvf.outsideScvIdx())), 0);
+                if (problem.pointSourceMap().count(key) || problem.pointSourceMap().count(keyOutside))
+                {