Commit 9e251a1e authored by Dennis Gläser's avatar Dennis Gläser
Browse files

Merge branch 'feature/tracer-multiphase-test-zero-saturation' into 'master'

Feature/tracer multiphase test zero saturation

Closes #704

See merge request !1610
parents 16aeb73a f0f325d4
......@@ -344,7 +344,7 @@ public:
// assemble the undeflected residual
const auto residual = this->evalLocalResidual()[0];
const auto storageResidual = this->evalLocalStorageResidual();
const auto storageResidual = this->evalLocalStorageResidual()[0];
//////////////////////////////////////////////////////////////////////////////////////////////////
// Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
......@@ -386,7 +386,7 @@ public:
// the residual with the deflected primary variables
elemSol[0][pvIdx] = priVar;
curVolVars.update(elemSol, this->problem(), element, scv);
return this->evalLocalStorageResidual();
return this->evalLocalStorageResidual()[0];
};
// for non-ghosts compute the derivative numerically
......
......@@ -25,6 +25,8 @@
#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_LOCAL_ASSEMBLER_HH
#define DUMUX_DISCRETIZATION_CC_MPFA_O_LOCAL_ASSEMBLER_HH
#include <algorithm>
#include <dumux/discretization/cellcentered/mpfa/localassemblerbase.hh>
#include <dumux/discretization/cellcentered/mpfa/localassemblerhelper.hh>
#include <dumux/discretization/cellcentered/mpfa/computetransmissibility.hh>
......@@ -104,6 +106,14 @@ public:
A[zeroRowIndices.first] = 0.0;
handle.CA()[faceIdx] = 0.0;
handle.T()[faceIdx] = 0.0;
// reset outside transmissibilities on surface grids
static constexpr int dim = IV::Traits::GridView::dimension;
static constexpr int dimWorld = IV::Traits::GridView::dimensionworld;
if constexpr (dim < dimWorld)
std::for_each( handle.tijOutside()[faceIdx].begin(),
handle.tijOutside()[faceIdx].end(),
[] (auto& outsideTij) { outsideTij = 0.0; } );
}
}
}
......@@ -250,9 +260,11 @@ private:
{
if (abs(wijk[faceIdx][0][localDir]) < wijZeroThresh)
{
insideZeroWij = !curGlobalScvf.boundary();
if (!curIsDirichlet)
{
insideZeroWij = true;
faceMarkers.emplace_back( std::make_pair(curLocalDofIdx, faceIdx) );
}
}
}
......@@ -308,7 +320,8 @@ private:
const auto& otherLocalScvf = iv.localScvf(otherLocalScvfIdx);
const auto otherLocalDofIdx = otherLocalScvf.localDofIndex();
if (otherLocalDofIdx == curLocalDofIdx && insideZeroWij)
// check for zero transmissibilities (skip if inside has been zero already)
if (otherLocalDofIdx == curLocalDofIdx && !insideZeroWij)
if (abs(wijk[faceIdx][idxOnScvf][localDir]) < wijZeroThresh)
faceMarkers.emplace_back( std::make_pair(curLocalDofIdx, faceIdx) );
......
......@@ -628,22 +628,30 @@ private:
// Effective diffusion coefficients might get zero if saturation = 0.
// Compute epsilon to detect obsolete rows in the iv-local matrices during assembly
if constexpr (ModelTraits::numFluidPhases() > 1)
const auto& scv = *scvs(fvGeometry()).begin();
const auto& scvf = *scvfs(fvGeometry()).begin();
const auto& vv = elemVolVars()[scv];
const auto D = [&] ()
{
const auto& scv = *scvs(fvGeometry()).begin();
const auto& scvf = *scvfs(fvGeometry()).begin();
const auto& vv = elemVolVars()[scv];
const auto& D = vv.diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
const auto tij = computeTpfaTransmissibility(scvf, scv, D, vv.extrusionFactor())*scvf.area();
// use transmissibility with molecular coefficient for epsilon estimate
localAssembler.assembleMatrices(handle.diffusionHandle(), iv, getD, tij*1e-7);
}
else
localAssembler.assembleMatrices(handle.diffusionHandle(), iv, getD);
// diffusion coefficients below 1e-20 are treated as zeroes!!
using std::max;
if constexpr (!FluidSystem::isTracerFluidSystem())
return max(1e-20, vv.diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx));
else
return max(1e-20, vv.diffusionCoefficient(0, 0, compIdx));
} ();
auto eps = 1e-7*computeTpfaTransmissibility(scvf, scv, D, vv.extrusionFactor())*scvf.area();
localAssembler.assembleMatrices(handle.diffusionHandle(), iv, getD, eps);
}
// assemble vector of mole fractions
auto getMassOrMoleFraction = [phaseIdx, compIdx] (const auto& volVars) { return (DiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged) ?volVars.massFraction(phaseIdx, compIdx) : volVars.moleFraction(phaseIdx, compIdx); };
auto getMassOrMoleFraction = [phaseIdx, compIdx] (const auto& volVars)
{
return (DiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged) ? volVars.massFraction(phaseIdx, compIdx) :
volVars.moleFraction(phaseIdx, compIdx);
};
localAssembler.assembleU(handle.diffusionHandle(), iv, getMassOrMoleFraction);
}
}
......
......@@ -26,6 +26,8 @@
#ifndef DUMUX_TRACER_LOCAL_RESIDUAL_HH
#define DUMUX_TRACER_LOCAL_RESIDUAL_HH
#include <dune/common/exceptions.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/discretization/method.hh>
......@@ -80,6 +82,12 @@ public:
{
NumEqVector storage(0.0);
// regularize saturation so we don't get singular matrices when the saturation is zero
// note that the fluxes will still be zero (zero effective diffusion coefficient),
// and we still solve the equation storage = 0 yielding the correct result
using std::max;
const Scalar saturation = max(1e-8, volVars.saturation(phaseIdx));
// formulation with mole balances
if (useMoles)
{
......@@ -87,7 +95,7 @@ public:
storage[compIdx] += volVars.porosity()
* volVars.molarDensity(phaseIdx)
* volVars.moleFraction(phaseIdx, compIdx)
* volVars.saturation(phaseIdx);
* saturation;
}
// formulation with mass balances
else
......@@ -96,8 +104,7 @@ public:
storage[compIdx] += volVars.porosity()
* volVars.density(phaseIdx)
* volVars.massFraction(phaseIdx, compIdx)
* volVars.saturation(phaseIdx);
* saturation;
}
return storage;
......@@ -190,9 +197,15 @@ public:
const VolumeVariables& curVolVars,
const SubControlVolume& scv) const
{
// regularize saturation so we don't get singular matrices when the saturation is zero
// note that the fluxes will still be zero (zero effective diffusion coefficient),
// and we still solve the equation storage = 0 yielding the correct result
using std::max;
const auto saturation = max(1e-8, curVolVars.saturation(phaseIdx));
const auto porosity = curVolVars.porosity();
const auto rho = useMoles ? curVolVars.molarDensity() : curVolVars.density();
const auto d_storage = scv.volume()*porosity*rho/this->timeLoop().timeStepSize();
const auto d_storage = scv.volume()*porosity*rho*saturation/this->timeLoop().timeStepSize();
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
partialDerivatives[compIdx][compIdx] += d_storage;
......@@ -229,6 +242,9 @@ public:
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
if constexpr (FVElementGeometry::GridGeometry::discMethod != DiscretizationMethod::cctpfa)
DUNE_THROW(Dune::NotImplemented, "Analytic flux differentiation only implemented for tpfa");
// advective term: we do the same for all tracer components
auto rho = [](const VolumeVariables& volVars)
{ return useMoles ? volVars.molarDensity() : volVars.density(); };
......@@ -273,7 +289,8 @@ public:
DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
derivativeMatrices[scvf.insideScvIdx()][compIdx][compIdx] += (advDerivII + diffDeriv);
derivativeMatrices[scvf.outsideScvIdx()][compIdx][compIdx] += (advDerivIJ - diffDeriv);
if (!scvf.boundary())
derivativeMatrices[scvf.outsideScvIdx()][compIdx][compIdx] += (advDerivIJ - diffDeriv);
}
}
......
add_subdirectory(2ptracer)
add_subdirectory(constvel)
add_subdirectory(multiphase)
dune_symlink_to_source_files(FILES "params.input")
dune_add_test(NAME test_tracer_multiphase_tpfa
LABELS porousmediumflow tracer
SOURCES main.cc
COMPILE_DEFINITIONS TYPETAG=TracerTestTpfa
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_tracer_multiphase_tpfa-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_tracer_multiphase_tpfa-00010.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_tracer_multiphase_tpfa params.input -Problem.Name test_tracer_multiphase_tpfa")
dune_add_test(NAME test_tracer_multiphase_mpfa
LABELS porousmediumflow tracer
SOURCES main.cc
COMPILE_DEFINITIONS TYPETAG=TracerTestMpfa
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_tracer_multiphase_tpfa-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_tracer_multiphase_mpfa-00010.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_tracer_multiphase_mpfa params.input -Problem.Name test_tracer_multiphase_mpfa")
// -*- 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 TracerTests
* \brief The test fluidsystem
*/
#ifndef DUMUX_TRACER_MULTIPHASE_TEST_FLUIDSYSTEM_HH
#define DUMUX_TRACER_MULTIPHASE_TEST_FLUIDSYSTEM_HH
#include <dumux/material/fluidsystems/base.hh>
namespace Dumux::FluidSystems {
//! A simple fluid system with two independent tracer components
template<class Scalar>
class TracerTest : public Base<Scalar, TracerTest<Scalar>>
{
public:
static constexpr bool isTracerFluidSystem()
{ return true; }
//! The number of components
static constexpr int numComponents = 2;
//! Human readable component name (index compIdx) (for vtk output)
static std::string componentName(int compIdx)
{ return "tracer_" + std::to_string(compIdx); }
//! Human readable phase name (index phaseIdx) (for velocity vtk output)
static std::string phaseName(int phaseIdx = 0)
{ return "aq"; }
//! Molar mass in kg/mol of the component with index compIdx
static Scalar molarMass(unsigned int compIdx)
{ return 0.300; }
//! Binary diffusion coefficient
//! (might depend on spatial parameters like pressure / temperature)
template<class Problem, class Element, class SubControlVolume>
static Scalar binaryDiffusionCoefficient(unsigned int compIdx,
const Problem& problem,
const Element& element,
const SubControlVolume& scv)
{
static const Scalar D = getParam<Scalar>("Problem.D");
static const Scalar D2 = getParam<Scalar>("Problem.D2");
if (compIdx == 0)
return D;
else
return D2;
}
/*!
* \copydoc Dumux::FluidSystems::Base::isCompressible
*/
static constexpr bool isCompressible(int phaseIdx)
{ return false; }
/*!
* \copydoc Dumux::FluidSystems::Base::viscosityIsConstant
*/
static constexpr bool viscosityIsConstant(int phaseIdx)
{ return true; }
};
} // end namespace Dumux::FluidSystems
#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 TracerTests
* \brief Test for the tracer CC model.
*/
#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/vtksequencewriter.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/linear/pdesolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
#include "properties.hh"
int main(int argc, char** argv) try
{
using namespace Dumux;
//! define the type tag for this problem
using TypeTag = Properties::TTag::TYPETAG;
//! 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
Parameters::init(argc, argv);
// try to create a grid (from the given grid file or the input file)
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
// setup instationary non-linear problem on this grid
////////////////////////////////////////////////////////////
//! we compute on the leaf grid view
const auto& leafGridView = gridManager.grid().leafGridView();
//! create the finite volume grid geometry
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
gridGeometry->update();
//! the problem (initial and boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
//! the solution vector
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(gridGeometry->numDofs());
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);
//! intialize the vtk output module
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
vtkWriter.addVolumeVariable([](const auto& v){ return v.saturation(); }, "S_aq");
vtkWriter.write(0.0);
//! get some time loop parameters
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
auto dt = getParam<Scalar>("TimeLoop.Dt");
//! instantiate time loop
auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, dt, tEnd);
//! the assembler with time loop for instationary problem
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric, /*implicit=*/false>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop);
assembler->setPreviousSolution(xOld);
//! the linear solver
using LinearSolver = ExplicitDiagonalSolver;
auto linearSolver = std::make_shared<LinearSolver>();
//! the system solver
LinearPDESolver<Assembler, LinearSolver> solver(assembler, linearSolver);
/////////////////////////////////////////////////////////////////////////////////////////////////
// run instationary non-linear simulation
/////////////////////////////////////////////////////////////////////////////////////////////////
//! set some check points for the time loop
timeLoop->setPeriodicCheckPoint(tEnd/10.0);
//! start the time loop
timeLoop->start();
while (!timeLoop->finished())
{
// assemble & solve
solver.solve(x);
// make the new solution the old solution
xOld = x;
gridVariables->advanceTimeStep();
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
// report statistics of this time step
timeLoop->reportTimeStep();
// write vtk output on check points
if (timeLoop->isCheckPoint() || timeLoop->finished())
vtkWriter.write(timeLoop->time());
// set new dt TODO implement CFL-criterion
timeLoop->setTimeStepSize(dt);
}
timeLoop->finalize(leafGridView.comm());
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
//! print dumux end message
if (mpiHelper.rank() == 0)
Parameters::print();
return 0;
}
catch (const Dumux::ParameterException &e)
{
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
catch (const Dune::DGFException & e)
{
std::cerr << "DGF exception thrown (" << e <<
"). Most likely, the DGF file name is wrong "
"or the DGF file is corrupted, "
"e.g. missing hash at end of file or wrong number (dimensions) of entries."
<< " ---> Abort!" << std::endl;
return 2;
}
catch (const Dune::Exception &e)
{
std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
return 3;
}
catch (...)
{
std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
return 4;
}
[TimeLoop]
Dt = 10000 # [s]
TEnd = 1e7 # [s]
[Grid]
UpperRight = 1 1
Cells = 50 50
[Problem]
Name = tracer
D = 0
D2 = 1e-8
FrontVelocity = 1e-7
Porosity = 0.5
[Vtk]
AddVelocity = true
[Assembly]
NumericDifference.BaseEpsilon = 1
// -*- 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 TracerTests
* \brief Definition of a problem for the tracer problem:
* A rotating velocity field mixes a tracer band in a porous groundwater reservoir.
*/
#ifndef DUMUX_TRACER_MULTIPHASE_TEST_PROBLEM_HH
#define DUMUX_TRACER_MULTIPHASE_TEST_PROBLEM_HH
#include <dumux/porousmediumflow/problem.hh>
namespace Dumux {
/*!
* \ingroup TracerTests
* \brief A problem for the multiphase tracer model
*/
template <class TypeTag>
class TracerTest : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using GridView = typename GridGeometry::GridView;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
//! property that defines whether mole or mass fractions are used
static constexpr bool useMoles = getPropValue<TypeTag, Properties::UseMoles>();
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
TracerTest(std::shared_ptr<const GridGeometry> gg)
: ParentType(gg)
{
// stating in the console whether mole or mass fractions are used
if(useMoles)
std::cout<<"problem uses mole fractions" << '\n';
else
std::cout<<"problem uses mass fractions" << '\n';
}
/*!
* \brief The boundary types
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const