Commit 7c8d7cad authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'feature/pnm-noncreepingflow' into 'master'

Feature/pnm noncreepingflow

See merge request !2596
parents fcea1723 f9e31fee
......@@ -6,6 +6,7 @@ Differences Between DuMu<sup>x</sup> 3.4 and DuMu<sup>x</sup> 3.3
- A quasi-static 2p PNM for the creation of pc-Sw curves has been added.
- A new namespace `Dumux::PoreNetwork` has been introduced, containing all relevant classes and functions for pore-network modeling.
- An example introduces the use of the 1p PNM for the estimation of the upscaled Darcy permeability.
- Advection type including inertia effect for simulations in the non-creeping flow regime
- Note that this is still considered a rather _experimental_ feature. Everything within namespace `Dumux::PoreNetwork` might undergo (backward-compatibility breaking) changes _without prior notice_.
- __Several scripts have been translated to Python__:
- `getusedversions.sh` to extract the used Dumux/Dune versions of a module (new script: `bin/util/getusedversions.py`)
......
......@@ -150,6 +150,204 @@ public:
}
};
/*!
* \file
* \ingroup PoreNetworkFlux
* \brief Non-creeping flow flux law to describe the advective flux for pore-network models based on
* El-Zehairy et al.(2019).
* https://doi.org/10.1016/j.advwatres.2019.103378
*/
template <class ScalarT, class... TransmissibilityLawTypes>
class NonCreepingFlow
{
public:
//! Export the Scalar type
using Scalar = ScalarT;
//! Export the creeping flow transmissibility law
using TransmissibilityCreepingFlow = Detail::Transmissibility<TransmissibilityLawTypes...>;
//! Inherit transmissibility from creeping flow transmissibility law to cache non-creeping flow-related parameters
class Transmissibility: public TransmissibilityCreepingFlow
{
public:
class SinglePhaseCache : public TransmissibilityCreepingFlow::SinglePhaseCache
{
using ParentType = typename TransmissibilityCreepingFlow::SinglePhaseCache;
public:
using Scalar = ScalarT;
template<class Problem, class Element, class FVElementGeometry, class ElementVolumeVariables, class FluxVariablesCache>
void fill(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const typename FVElementGeometry::SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
const FluxVariablesCache& fluxVarsCache,
const int phaseIdx)
{
ParentType::fill(problem, element, fvGeometry,scvf, elemVolVars, fluxVarsCache, phaseIdx);
const auto elemSol = elementSolution(element, elemVolVars, fvGeometry);
for (const auto& scv : scvs(fvGeometry))
{
const auto localIdx = scv.indexInElement();
const Scalar throatToPoreAreaRatio = fluxVarsCache.throatCrossSectionalArea() / problem.spatialParams().poreCrossSectionalArea(element, scv, elemSol);
// dimensionless momentum coefficient
const Scalar momentumCoefficient = problem.spatialParams().momentumCoefficient(element, scv, elemSol);
// dimensionless kinetik-energy coefficient
const Scalar kineticEnergyCoefficient = problem.spatialParams().kineticEnergyCoefficient(element, scv, elemSol);
expansionCoefficient_[localIdx] = getExpansionCoefficient_(throatToPoreAreaRatio, momentumCoefficient, kineticEnergyCoefficient);
contractionCoefficient_[localIdx] = getContractionCoefficient_(throatToPoreAreaRatio, momentumCoefficient, kineticEnergyCoefficient);
}
}
Scalar expansionCoefficient(int downstreamIdx) const
{ return expansionCoefficient_[downstreamIdx]; }
Scalar contractionCoefficient(int upstreamIdx) const
{ return contractionCoefficient_[upstreamIdx]; }
private:
Scalar getExpansionCoefficient_(const Scalar throatToPoreAreaRatio, const Scalar momentumCoefficient, const Scalar kineticEnergyCoefficient) const
{
Scalar expansionCoefficient = throatToPoreAreaRatio * throatToPoreAreaRatio * (2 * momentumCoefficient - kineticEnergyCoefficient)
+ kineticEnergyCoefficient - 2 * momentumCoefficient * throatToPoreAreaRatio
- (1 - throatToPoreAreaRatio * throatToPoreAreaRatio);
return expansionCoefficient;
}
Scalar getContractionCoefficient_(const Scalar throatToPoreAreaRatio, const Scalar momentumCoefficient, const Scalar kineticEnergyCoefficient) const
{
const Scalar contractionAreaRatio = getContractionAreaRatio_(throatToPoreAreaRatio);
Scalar contractionCoefficient = (1 - (throatToPoreAreaRatio * throatToPoreAreaRatio * kineticEnergyCoefficient - 2 * momentumCoefficient /*+1-throatToPoreAreaRatio*throatToPoreAreaRatio*/)
* contractionAreaRatio * contractionAreaRatio - 2 * contractionAreaRatio) / (contractionAreaRatio * contractionAreaRatio);
return contractionCoefficient;
}
Scalar getContractionAreaRatio_(const Scalar throatToPoreAreaRatio) const
{
return 1.0-(1.0-throatToPoreAreaRatio)/(2.08*(1.0-throatToPoreAreaRatio)+0.5371);
}
std::array<Scalar, 2> expansionCoefficient_;
std::array<Scalar, 2> contractionCoefficient_;
};
};
/*!
* \brief Returns the advective flux of a fluid phase
* across the given sub-control volume face (corresponding to a pore throat).
* \note The flux is given in N*m, and can be converted
* into a volume flux (m^3/s) or mass flux (kg/s) by applying an upwind scheme
* for the mobility (1/viscosity) or the product of density and mobility, respectively.
*/
template<class Problem, class Element, class FVElementGeometry,
class ElementVolumeVariables, class SubControlVolumeFace, class ElemFluxVarsCache>
static Scalar flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const int phaseIdx,
const ElemFluxVarsCache& elemFluxVarsCache)
{
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
// Get the inside and outside volume variables
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
const auto& insideVolVars = elemVolVars[insideScv];
const auto& outsideVolVars = elemVolVars[outsideScv];
// calculate the pressure difference
Scalar deltaP = insideVolVars.pressure(phaseIdx) - outsideVolVars.pressure(phaseIdx);
// add gravity term
static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
if (enableGravity)
{
const Scalar rho = 0.5*insideVolVars.density(phaseIdx) + 0.5*outsideVolVars.density(phaseIdx);
// projected distance between pores in gravity field direction
const auto poreToPoreVector = fluxVarsCache.poreToPoreDistance() * scvf.unitOuterNormal();
const auto& gravityVector = problem.spatialParams().gravity(scvf.center());
deltaP += rho * (gravityVector * poreToPoreVector);
}
const Scalar creepingFlowTransmissibility = fluxVarsCache.transmissibility(phaseIdx);
const Scalar throatCrossSectionalArea = fluxVarsCache.throatCrossSectionalArea();
assert(scvf.insideScvIdx() == 0);
assert(scvf.outsideScvIdx() == 1);
// determine the flow direction to predict contraction and expansion
const auto [upstreamIdx, downstreamIdx] = deltaP > 0 ? std::pair(0, 1) : std::pair(1, 0);
const Scalar contractionCoefficient = fluxVarsCache.singlePhaseFlowVariables().contractionCoefficient(upstreamIdx);
const Scalar expansionCoefficient = fluxVarsCache.singlePhaseFlowVariables().expansionCoefficient(downstreamIdx);
const Scalar mu = elemVolVars[upstreamIdx].viscosity();
//! El-Zehairy et al.(2019): Eq.(14):
//! A0Q^2 + B0Q + C0 = 0 where Q is volumetric flowrate, A0 = [Kc + Ke] * rho /[2aij^2], Kc and Ke are contraction and expansion coefficient respectively,
//! aij is cross sectional area of the throat, B0 = 1/K0 (K0 is trnasmissibility used in paper) and C0 = -deltaP
//! In our implementation viscosity of fluid will be included using upwinding later. Thus, creepingFlowTransmissibility = mu * K0 and
//! the volumetric flowrate calculated here q = mu * Q. Substitution of new variables into El-Zehairy et al.(2019), Eq.(14) gives:
//! Aq^2 + Bq + C = 0 where
//! A = A0 / mu^2, B = B0/mu = 1/creepingFlowTransmissibility and C = C0
//! attention: the q, volumetric flowrate, calculated here is always positive and its sign needs to be determined based on flow direction
//! this approach is taken to prevent the term under the square root (discriminant) becoming negative
const Scalar A = (contractionCoefficient * elemVolVars[upstreamIdx].density() + expansionCoefficient * elemVolVars[downstreamIdx].density())
/ (2.0 * mu * mu * throatCrossSectionalArea * throatCrossSectionalArea);
const Scalar B = 1.0/creepingFlowTransmissibility;
using std::abs;
const Scalar C = -abs(deltaP);
using std::sqrt;
const auto discriminant = B*B - 4*A*C;
const auto q = (-B + sqrt(discriminant)) / (2*A);
//! give the volume flowrate proper sign based on flow direction.
if (upstreamIdx == 0)
return q;
else
return -q;
}
template<class Problem, class Element, class FVElementGeometry, class ElementVolumeVariables, class FluxVariablesCache>
static Scalar calculateTransmissibility(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const typename FVElementGeometry::SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
const FluxVariablesCache& fluxVarsCache,
const int phaseIdx)
{
static_assert(ElementVolumeVariables::VolumeVariables::numFluidPhases() == 1);
return Transmissibility::singlePhaseTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, phaseIdx);
}
template<class Problem, class Element, class FVElementGeometry, class ElementVolumeVariables, class FluxVariablesCache>
static std::array<Scalar, 2> calculateTransmissibilities(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const typename FVElementGeometry::SubControlVolumeFace& scvf,
const FluxVariablesCache& fluxVarsCache)
{
static_assert(ElementVolumeVariables::VolumeVariables::numFluidPhases() == 1);
const Scalar t = calculateTransmissibility(problem, element, fvGeometry, scvf, elemVolVars, fluxVarsCache, 0);
return {t, -t};
}
};
} // end namespace Dumux
......
add_subdirectory(nonisothermal)
add_subdirectory(noncreepingflow)
add_input_file_links()
dune_symlink_to_source_files(FILES grids)
......
add_subdirectory(nonisothermal)
add_input_file_links()
dune_add_test(NAME test_pnm_1p_noncreeping_flow
SOURCES main.cc
LABELS porenetworkflow
COMPILE_DEFINITIONS ISOTHERMAL=1
CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )"
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_pnm_1p_noncreepingflow-reference.vtp
${CMAKE_CURRENT_BINARY_DIR}/test_pnm_1p_noncreepingflow-00000.vtp
--command "${CMAKE_CURRENT_BINARY_DIR}/test_pnm_1p_noncreeping_flow")
// -*- 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
*
* \brief test for the pore-network model with non-creeping flow
*/
#include <config.h>
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/porenetwork/common/pnmvtkoutputmodule.hh>
#include <dumux/porenetwork/common/boundaryflux.hh>
#include <dumux/io/grid/porenetwork/gridmanager.hh>
#include "properties.hh"
int main(int argc, char** argv)
{
using namespace Dumux;
using TypeTag = Properties::TTag::PNMOnePNonCreepingProblem;
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
// print dumux start message
if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/true);
////////////////////////////////////////////////////////////
// parse the command line arguments and input file
////////////////////////////////////////////////////////////
// parse command line arguments
Parameters::init(argc, argv);
//////////////////////////////////////////////////////////////////////
// try to create a grid (from the given grid file or the input file)
/////////////////////////////////////////////////////////////////////
using GridManager = Dumux::PoreNetwork::GridManager<3>;
GridManager gridManager;
gridManager.init();
// we compute on the leaf grid view
const auto& leafGridView = gridManager.grid().leafGridView();
auto gridData = gridManager.getGridData();
// create the finite volume grid geometry
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
gridGeometry->update(*gridData);
// the spatial parameters
using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
auto spatialParams = std::make_shared<SpatialParams>(gridGeometry);
// the problem (boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry, spatialParams);
// the solution vector
using GridView = typename GridGeometry::GridView;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(leafGridView.size(GridView::dimension));
problem->applyInitialSolution(x);
auto xOld = x;
// the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
// the assembler
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
// helper class to calculate the boundary flux
const auto boundaryFlux = Dumux::PoreNetwork::BoundaryFlux(*gridVariables, assembler->localResidual(), x);
// the linear solver
using LinearSolver = ILU0RestartedGMResBackend;
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
nonLinearSolver.solve(x);
// output result to vtk
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
Dumux::PoreNetwork::VtkOutputModule<GridVariables, GetPropType<TypeTag, Properties::FluxVariables>, SolutionVector>
vtkWriter(*gridVariables, x, problem->name());
IOFields::initOutputModule(vtkWriter);
vtkWriter.write(0.0);
if (mpiHelper.rank() == 0)
{
std::cout << "cumulative outflux is: "
<< boundaryFlux.getFlux("max", 0, true) << std::endl;
DumuxMessage::print();
Parameters::print();
}
return 0;
}
add_input_file_links()
dune_add_test(NAME test_pnm_1pni_noncreeping_flow
SOURCES main.cc
LABELS porenetworkflow
COMPILE_DEFINITIONS ISOTHERMAL=0
CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )"
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_pnm_1pni_noncreepingflow-reference.vtp
${CMAKE_CURRENT_BINARY_DIR}/test_pnm_1pni_noncreepingflow-00014.vtp
--command "${CMAKE_CURRENT_BINARY_DIR}/test_pnm_1pni_noncreeping_flow")
// -*- 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
*
* \brief test for the non-isothermal pore-network model with non-creeping flow
*/
#include <config.h>
#include <ctime>
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/grid/io/file/dgfparser/dgfexception.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/porenetwork/common/pnmvtkoutputmodule.hh>
#include <dumux/porenetwork/common/boundaryflux.hh>
#include <dumux/io/grid/porenetwork/gridmanager.hh>
#include "../properties.hh"
int main(int argc, char** argv)
{
using namespace Dumux;
using TypeTag = Properties::TTag::PNMOnePNonCreepingProblem;
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
// print dumux start message
if (mpiHelper.rank() == 0)
DumuxMessage::print(/*firstCall=*/true);
////////////////////////////////////////////////////////////
// parse the command line arguments and input file
////////////////////////////////////////////////////////////
// parse command line arguments
Parameters::init(argc, argv);
//////////////////////////////////////////////////////////////////////
// try to create a grid (from the given grid file or the input file)
/////////////////////////////////////////////////////////////////////
using GridManager = Dumux::PoreNetwork::GridManager<3>;
GridManager gridManager;
gridManager.init();
// we compute on the leaf grid view
const auto& leafGridView = gridManager.grid().leafGridView();
auto gridData = gridManager.getGridData();
// create the finite volume grid geometry
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
gridGeometry->update(*gridData);
// the spatial parameters
using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
auto spatialParams = std::make_shared<SpatialParams>(gridGeometry);
// the problem (boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry, spatialParams);
// the solution vector
using GridView = typename GridGeometry::GridView;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(leafGridView.size(GridView::dimension));
problem->applyInitialSolution(x);
auto xOld = x;
// the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
// get some time loop parameters
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
// intialize the vtk output module
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
Dumux::PoreNetwork::VtkOutputModule<GridVariables, GetPropType<TypeTag, Properties::FluxVariables>, SolutionVector>
vtkWriter(*gridVariables, x, problem->name());
IOFields::initOutputModule(vtkWriter); //! Add model specific output fields
vtkWriter.write(0.0);
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
// the assembler with time loop for instationary problem
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
// the linear solver
using LinearSolver = ILU0BiCGSTABBackend;
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
// time loop
timeLoop->start(); do
{
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
// solve the non-linear system
nonLinearSolver.solve(x, *timeLoop);
// make the new solution the old solution
xOld = x;
gridVariables->advanceTimeStep();
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
// write vtk output
vtkWriter.write(timeLoop->time());
// report statistics of this time step
timeLoop->reportTimeStep();
// set new dt as suggested by newton solver
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
} while (!timeLoop->finished());
timeLoop->finalize(leafGridView.comm());
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
// print dumux end message
if (mpiHelper.rank() == 0)
{
Parameters::print();
DumuxMessage::print(/*firstCall=*/false);
}
return 0;
}
[TimeLoop]
DtInitial = 1e-3 # [s]
TEnd = 1 # [s]
[Grid]
File = ./../../grids/pnm_converted.dgf
PoreGeometry = Cube
ThroatCrossSectionShape = Square
[Problem]
Name = test_pnm_1pni_noncreepingflow # name passed to the output routines
Verbose = true
InletPressure = 1e5
OutletPressure = 0.9e5
EnableGravity = false
[Vtk]
AddVelocity = true
[Component]
SolidHeatCapacity = 1.0 # for compatibility
SolidDensity = 1.0 # for compatibility
SolidThermalConductivity = 1.0 # for compatibility
[Grid]
File = ./../grids/pnm_converted.dgf
PoreGeometry = Cube
ThroatCrossSectionShape = Square
[Problem]
Name = test_pnm_1p_noncreepingflow # name passed to the output routines
InletPressure = 1e5