Commit 46bf2b34 authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

Merge branch 'feature/tpfa-facet-gravity-on-surface-grids' into 'master'

Feature/tpfa facet gravity on surface grids

See merge request !2259
parents 94b468f4 3a270614
......@@ -367,6 +367,9 @@ public:
//! Export transmissibility storage type
using AdvectionTransmissibilityContainer = std::array<Scalar, 2>;
//! Export the type used for the gravity coefficients
using GravityCoefficients = std::array<Scalar, 1>;
//! update subject to a given problem
template< class Problem, class ElementVolumeVariables >
void updateAdvection(const Problem& problem,
......@@ -375,7 +378,7 @@ public:
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace &scvf)
{
tij_ = AdvectionType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf);
tij_ = AdvectionType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, g_);
}
//! We use the same name as in the TpfaDarcysLawCache so
......@@ -393,8 +396,13 @@ public:
Scalar advectionTijFacet() const
{ return tij_[facetTijIdx]; }
//! return the coefficients for the computation of gravity at the scvf
const GravityCoefficients& gravityCoefficients() const
{ return g_; }
private:
AdvectionTransmissibilityContainer tij_;
GravityCoefficients g_;
};
/*!
......@@ -436,24 +444,49 @@ class CCTpfaFacetCouplingDarcysLawImpl<ScalarType, GridGeometry, /*isNetwork*/tr
int phaseIdx,
const ElementFluxVarsCache& elemFluxVarsCache)
{
static const Scalar gravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
if (gravity)
DUNE_THROW(Dune::NotImplemented, "Gravity for darcys law with facet coupling on surface grids");
if (!problem.couplingManager().isOnInteriorBoundary(element, scvf))
return TpfaDarcysLaw::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache);
// On surface grids only xi = 1.0 can be used, as the coupling condition
// for xi != 1.0 does not generalize for surface grids where there can be
// seveal neighbor meeting at a branching point.
static const Scalar xi = getParamFromGroup<Scalar>(problem.paramGroup(), "FacetCoupling.Xi", 1.0);
if (Dune::FloatCmp::ne(xi, 1.0, 1e-6))
DUNE_THROW(Dune::InvalidStateException, "Xi != 1.0 cannot be used on surface grids");
// Obtain inside and fracture pressures
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf);
const auto pInside = insideVolVars.pressure(phaseIdx);
const auto pFacet = problem.couplingManager().getLowDimVolVars(element, scvf).pressure(phaseIdx);
const auto pFacet = facetVolVars.pressure(phaseIdx);
// return flux
// compute and return flux
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
if (scvf.boundary())
return fluxVarsCache.advectionTijInside()*pInside + fluxVarsCache.advectionTijFacet()*pFacet;
else
return fluxVarsCache.advectionTijInside()*pInside + fluxVarsCache.advectionTijFacet()*pFacet;
Scalar flux = fluxVarsCache.advectionTijInside()*pInside + fluxVarsCache.advectionTijFacet()*pFacet;
static const Scalar gravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
if (gravity)
{
// compute alpha := n^T*K*g and add to flux (use arithmetic mean for density)
const auto& g = problem.spatialParams().gravity(scvf.ipGlobal());
const auto rho = 0.5*(insideVolVars.density(phaseIdx) + facetVolVars.density(phaseIdx));
const auto rhoTimesArea = rho*Extrusion::area(scvf);
const auto alpha_inside = rhoTimesArea*insideVolVars.extrusionFactor()
*vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g);
flux += alpha_inside;
// maybe add further gravitational contributions
if ( !scvf.boundary() && problem.interiorBoundaryTypes(element, scvf).hasOnlyNeumann() )
{
const auto alpha_facet = rhoTimesArea*insideVolVars.extrusionFactor()
*vtmv(scvf.unitOuterNormal(), facetVolVars.permeability(), g);
flux -= fluxVarsCache.gravityCoefficients()[0]*(alpha_inside - alpha_facet);
}
}
return flux;
}
// The flux variables cache has to be bound to an element prior to flux calculations
......@@ -464,6 +497,20 @@ class CCTpfaFacetCouplingDarcysLawImpl<ScalarType, GridGeometry, /*isNetwork*/tr
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf)
{
typename Cache::GravityCoefficients g;
return calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, g);
}
// This overload additionally receives a container in which the coefficients required
// for the computation of the gravitational acceleration ar the scvf are stored
template< class Problem, class ElementVolumeVariables >
static TijContainer calculateTransmissibility(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
typename Cache::GravityCoefficients& g)
{
TijContainer tij;
if (!problem.couplingManager().isCoupled(element, scvf))
......@@ -504,7 +551,8 @@ class CCTpfaFacetCouplingDarcysLawImpl<ScalarType, GridGeometry, /*isNetwork*/tr
*vtmv(scvf.unitOuterNormal(), facetVolVars.permeability(), scvf.unitOuterNormal());
// TODO: check for division by zero??
tij[Cache::insideTijIdx] = wFacet*wIn/(wIn+wFacet);
g[0] = wIn/(wIn+wFacet);
tij[Cache::insideTijIdx] = wFacet*g[0];
tij[Cache::facetTijIdx] = -tij[Cache::insideTijIdx];
}
else if (iBcTypes.hasOnlyDirichlet())
......
......@@ -70,3 +70,39 @@ dumux_add_test(NAME test_md_facet_1p1p_gravity_xi066_mpfa
--command "${CMAKE_CURRENT_BINARY_DIR}/test_md_facet_1p1p_gravity_mpfa params.input \
-FacetCoupling.Xi 0.66
-Vtk.OutputName test_md_facet_1p1p_gravity_xi066_mpfa")
dumux_add_test(NAME test_md_facet_1p1p_gravity_surface_tpfa
LABELS multidomain multidomain_facet 1p
CMAKE_GUARD "( dune-foamgrid_FOUND AND dune-alugrid_FOUND )"
SOURCES main.cc
COMPILE_DEFINITIONS BULKTYPETAG=OnePBulkTpfa
COMPILE_DEFINITIONS LOWDIMTYPETAG=OnePLowDimTpfa
COMPILE_DEFINITIONS BULKGRIDTYPE=Dune::ALUGrid<2,3,Dune::cube,Dune::nonconforming>
COMPILE_DEFINITIONS LOWDIMGRIDTYPE=Dune::FoamGrid<1,3>
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_CURRENT_BINARY_DIR}/test_md_facet_1p1p_gravity_surface_tpfa_bulk-00000.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_md_facet_1p1p_gravity_surface_tpfa_bulk-00001.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_md_facet_1p1p_gravity_surface_tpfa_lowdim-00000.vtp
${CMAKE_CURRENT_BINARY_DIR}/test_md_facet_1p1p_gravity_surface_tpfa_lowdim-00001.vtp
--command "${CMAKE_CURRENT_BINARY_DIR}/test_md_facet_1p1p_gravity_surface_tpfa params.input \
-Vtk.OutputName test_md_facet_1p1p_gravity_surface_tpfa \
-Grid.File grids/gravity_surface.msh")
dumux_add_test(NAME test_md_facet_1p1p_gravity_surface_mpfa
LABELS multidomain multidomain_facet 1p
CMAKE_GUARD "( dune-foamgrid_FOUND AND dune-alugrid_FOUND )"
SOURCES main.cc
COMPILE_DEFINITIONS BULKTYPETAG=OnePBulkMpfa
COMPILE_DEFINITIONS LOWDIMTYPETAG=OnePLowDimMpfa
COMPILE_DEFINITIONS BULKGRIDTYPE=Dune::ALUGrid<2,3,Dune::cube,Dune::nonconforming>
COMPILE_DEFINITIONS LOWDIMGRIDTYPE=Dune::FoamGrid<1,3>
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_CURRENT_BINARY_DIR}/test_md_facet_1p1p_gravity_surface_mpfa_bulk-00000.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_md_facet_1p1p_gravity_surface_mpfa_bulk-00001.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_md_facet_1p1p_gravity_surface_mpfa_lowdim-00000.vtp
${CMAKE_CURRENT_BINARY_DIR}/test_md_facet_1p1p_gravity_surface_mpfa_lowdim-00001.vtp
--command "${CMAKE_CURRENT_BINARY_DIR}/test_md_facet_1p1p_gravity_surface_mpfa params.input \
-Vtk.OutputName test_md_facet_1p1p_gravity_surface_mpfa \
-Grid.File grids/gravity_surface.msh")
This diff is collapsed.
......@@ -45,8 +45,8 @@
#include <dumux/io/vtkoutputmodule.hh>
#include "problem_bulk.hh"
#include "problem_lowdim.hh"
#include "properties_bulk.hh"
#include "properties_lowdim.hh"
namespace Dumux {
......
......@@ -24,64 +24,9 @@
#ifndef DUMUX_TEST_TPFAFACETCOUPLING_ONEP_BULKPROBLEM_HH
#define DUMUX_TEST_TPFAFACETCOUPLING_ONEP_BULKPROBLEM_HH
#include <dune/alugrid/grid.hh>
#include <dumux/common/boundarytypes.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/multidomain/facet/box/properties.hh>
#include <dumux/multidomain/facet/cellcentered/tpfa/properties.hh>
#include <dumux/multidomain/facet/cellcentered/mpfa/properties.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/porousmediumflow/1p/model.hh>
#include "spatialparams.hh"
// default for the bulk grid type
#ifndef BULKGRIDTYPE
#define BULKGRIDTYPE Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>
#endif
namespace Dumux {
// forward declarations
template<class TypeTag> class OnePBulkProblem;
namespace Properties {
// create the type tag nodes
namespace TTag {
struct OnePBulk { using InheritsFrom = std::tuple<OneP>; };
struct OnePBulkTpfa { using InheritsFrom = std::tuple<CCTpfaFacetCouplingModel, OnePBulk>; };
struct OnePBulkMpfa { using InheritsFrom = std::tuple<CCMpfaFacetCouplingModel, OnePBulk>; };
} // end namespace TTag
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::OnePBulk> { using type = BULKGRIDTYPE; };
// Set the problem type
template<class TypeTag>
struct Problem<TypeTag, TTag::OnePBulk> { using type = OnePBulkProblem<TypeTag>; };
// set the spatial params
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::OnePBulk>
{
using type = OnePSpatialParams< GetPropType<TypeTag, Properties::GridGeometry>,
GetPropType<TypeTag, Properties::Scalar> >;
};
// the fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::OnePBulk>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = FluidSystems::OnePLiquid< Scalar, Components::Constant<1, Scalar> >;
};
} // end namespace Properties
/*!
* \ingroup FacetTests
......
......@@ -24,65 +24,9 @@
#ifndef DUMUX_TEST_TPFAFACETCOUPLING_ONEP_LOWDIMPROBLEM_HH
#define DUMUX_TEST_TPFAFACETCOUPLING_ONEP_LOWDIMPROBLEM_HH
#include <dune/foamgrid/foamgrid.hh>
#include <dumux/common/boundarytypes.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/discretization/box.hh>
#include <dumux/discretization/cctpfa.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/porousmediumflow/1p/model.hh>
#include "spatialparams.hh"
// default for the grid type
#ifndef LOWDIMGRIDTYPE
#define LOWDIMGRIDTYPE Dune::FoamGrid<1, 2>
#endif
namespace Dumux {
// forward declarations
template<class TypeTag> class OnePLowDimProblem;
namespace Properties {
// create the type tag nodes
namespace TTag {
struct OnePLowDim { using InheritsFrom = std::tuple<OneP>; };
struct OnePLowDimTpfa { using InheritsFrom = std::tuple<OnePLowDim, CCTpfaModel>; };
// we need an additional type tag for the test using mpfa in the bulk domain
struct OnePLowDimMpfa { using InheritsFrom = std::tuple<OnePLowDim, CCTpfaModel>; };
} // end namespace TTag
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::OnePLowDim> { using type = LOWDIMGRIDTYPE; };
// Set the problem type
template<class TypeTag>
struct Problem<TypeTag, TTag::OnePLowDim> { using type = OnePLowDimProblem<TypeTag>; };
// set the spatial params
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::OnePLowDim>
{
using type = OnePSpatialParams< GetPropType<TypeTag, Properties::GridGeometry>,
GetPropType<TypeTag, Properties::Scalar> >;
};
// the fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::OnePLowDim>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = FluidSystems::OnePLiquid< Scalar, Components::Constant<1, Scalar> >;
};
} // end namespace Properties
/*!
* \ingroup FacetTests
......
// -*- 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 3 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 FacetTests
* \brief The properties for the bulk domain in the single-phase facet coupling test.
*/
#ifndef DUMUX_TEST_TPFAFACETCOUPLING_ONEP_BULK_PROPERTIES_HH
#define DUMUX_TEST_TPFAFACETCOUPLING_ONEP_BULK_PROPERTIES_HH
#include <dune/alugrid/grid.hh>
#include <dumux/common/boundarytypes.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/porousmediumflow/1p/model.hh>
#include <dumux/multidomain/facet/box/properties.hh>
#include <dumux/multidomain/facet/cellcentered/tpfa/properties.hh>
#include <dumux/multidomain/facet/cellcentered/mpfa/properties.hh>
#include "spatialparams.hh"
#include "problem_bulk.hh"
// default for the bulk grid type
#ifndef BULKGRIDTYPE
#define BULKGRIDTYPE Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>
#endif
namespace Dumux::Properties {
// create the type tag nodes
namespace TTag {
struct OnePBulk { using InheritsFrom = std::tuple<OneP>; };
struct OnePBulkTpfa { using InheritsFrom = std::tuple<CCTpfaFacetCouplingModel, OnePBulk>; };
struct OnePBulkMpfa { using InheritsFrom = std::tuple<CCMpfaFacetCouplingModel, OnePBulk>; };
} // end namespace TTag
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::OnePBulk> { using type = BULKGRIDTYPE; };
// Set the problem type
template<class TypeTag>
struct Problem<TypeTag, TTag::OnePBulk> { using type = OnePBulkProblem<TypeTag>; };
// set the spatial params
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::OnePBulk>
{
using type = OnePSpatialParams< GetPropType<TypeTag, Properties::GridGeometry>,
GetPropType<TypeTag, Properties::Scalar> >;
};
// the fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::OnePBulk>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = FluidSystems::OnePLiquid< Scalar, Components::Constant<1, Scalar> >;
};
} // end namespace Dumux::Properties
#endif
// -*- 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 3 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 FacetTests
* \brief The properties for the lower-dimensional domain in the single-phase facet coupling test.
*/
#ifndef DUMUX_TEST_TPFAFACETCOUPLING_ONEP_LOWDIM_PROPERTIES_HH
#define DUMUX_TEST_TPFAFACETCOUPLING_ONEP_LOWDIM_PROPERTIES_HH
#include <dune/foamgrid/foamgrid.hh>
#include <dumux/common/boundarytypes.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/porousmediumflow/1p/model.hh>
#include <dumux/discretization/box.hh>
#include <dumux/discretization/cctpfa.hh>
#include "spatialparams.hh"
#include "problem_lowdim.hh"
// default for the grid type
#ifndef LOWDIMGRIDTYPE
#define LOWDIMGRIDTYPE Dune::FoamGrid<1, 2>
#endif
namespace Dumux::Properties {
// create the type tag nodes
namespace TTag {
struct OnePLowDim { using InheritsFrom = std::tuple<OneP>; };
struct OnePLowDimTpfa { using InheritsFrom = std::tuple<OnePLowDim, CCTpfaModel>; };
// we need an additional type tag for the test using mpfa in the bulk domain
struct OnePLowDimMpfa { using InheritsFrom = std::tuple<OnePLowDim, CCTpfaModel>; };
} // end namespace TTag
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::OnePLowDim> { using type = LOWDIMGRIDTYPE; };
// Set the problem type
template<class TypeTag>
struct Problem<TypeTag, TTag::OnePLowDim> { using type = OnePLowDimProblem<TypeTag>; };
// set the spatial params
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::OnePLowDim>
{
using type = OnePSpatialParams< GetPropType<TypeTag, Properties::GridGeometry>,
GetPropType<TypeTag, Properties::Scalar> >;
};
// the fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::OnePLowDim>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = FluidSystems::OnePLiquid< Scalar, Components::Constant<1, Scalar> >;
};
} // end namespace Dumux::Properties
#endif
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment