Commit 969732f0 authored by Katharina Heck's avatar Katharina Heck Committed by Kilian Weishaupt
Browse files

[test][multidomain] add 1p2c_1p2c test to compare diffusion laws....

[test][multidomain] add 1p2c_1p2c test to compare diffusion laws. MaxwellStefan and Fick's Law produce the same solution for 2 components
parent ebbc432c
add_subdirectory("verticalflow")
add_subdirectory("horizontalflow")
add_subdirectory("diffusionlawcomparison")
add_input_file_links()
dune_add_test(NAME test_md_darcy1p2c_stokes1p2c_maxwellstefan
SOURCES main.cc
CMAKE_GUARD HAVE_UMFPACK
COMPILE_DEFINITIONS DIFFUSIONTYPE=MaxwellStefansLaw<TypeTag>
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_md_darcy1p2c_stokes1p2c_comparison_stokes-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_md_darcy1p2c_stokes1p2c_comparison_stokes-00020.vtu
${CMAKE_SOURCE_DIR}/test/references/test_md_darcy1p2c_stokes1p2c_comparison_darcy-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_md_darcy1p2c_stokes1p2c_comparison_darcy-00020.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_md_darcy1p2c_stokes1p2c_maxwellstefan params.input")
dune_add_test(NAME test_md_darcy1p2c_stokes1p2c_fickslaw
SOURCES main.cc
CMAKE_GUARD HAVE_UMFPACK
COMPILE_DEFINITIONS DIFFUSIONTYPE=FicksLaw<TypeTag>
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_md_darcy1p2c_stokes1p2c_comparison_stokes-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_md_darcy1p2c_stokes1p2c_comparison_stokes-00020.vtu
${CMAKE_SOURCE_DIR}/test/references/test_md_darcy1p2c_stokes1p2c_comparison_darcy-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_md_darcy1p2c_stokes1p2c_comparison_darcy-00020.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_md_darcy1p2c_stokes1p2c_fickslaw params.input")
// -*- 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 coupled Stokes/Darcy problem (1p2c) comparing different diffusion laws
*/
#include <config.h>
#include <ctime>
#include <iostream>
#include <fstream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/istl/io.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/geometry/diameter.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/discretization/method.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/multidomain/staggeredtraits.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/newtonsolver.hh>
#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
#include "problem_darcy.hh"
#include "problem_stokes.hh"
namespace Dumux {
namespace Properties {
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::StokesOnePTwoCTypeTag>
{
using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOnePTwoCTypeTag>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::DarcyOnePTwoCTypeTag>
{
using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesOnePTwoCTypeTag, Properties::TTag::StokesOnePTwoCTypeTag, TypeTag>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
} // end namespace Properties
} // end namespace Dumux
int main(int argc, char** argv) try
{
using namespace Dumux;
// 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 command line arguments and input file
Parameters::init(argc, argv);
// Define the sub problem type tags
using StokesTypeTag = Properties::TTag::StokesOnePTwoCTypeTag;
using DarcyTypeTag = Properties::TTag::DarcyOnePTwoCTypeTag;
// try to create a grid (from the given grid file or the input file)
// for both sub-domains
using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
DarcyGridManager darcyGridManager;
darcyGridManager.init("Darcy"); // pass parameter group
using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>;
StokesGridManager stokesGridManager;
stokesGridManager.init("Stokes"); // pass parameter group
// we compute on the leaf grid view
const auto& darcyGridView = darcyGridManager.grid().leafGridView();
const auto& stokesGridView = stokesGridManager.grid().leafGridView();
// create the finite volume grid geometry
using StokesFVGridGeometry = GetPropType<StokesTypeTag, Properties::FVGridGeometry>;
auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
stokesFvGridGeometry->update();
using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::FVGridGeometry>;
auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
darcyFvGridGeometry->update();
using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
// the coupling manager
using CouplingManager = StokesDarcyCouplingManager<Traits>;
auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
// the indices
constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
constexpr auto darcyIdx = CouplingManager::darcyIdx;
// the problems (initial and boundary conditions)
using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>;
auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
// initialize the fluidsystem (tabulation)
GetPropType<StokesTypeTag, Properties::FluidSystem>::init();
// get some time loop parameters
using Scalar = GetPropType<StokesTypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
// check if we are about to restart a previously interrupted simulation
Scalar restartTime = 0;
if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
restartTime = getParam<Scalar>("TimeLoop.Restart");
// instantiate time loop
auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
// control the injection period
const Scalar injectionBegin = getParam<Scalar>("Stokes.Problem.InjectionBegin");
const Scalar injectionEnd = getParam<Scalar>("Stokes.Problem.InjectionEnd");
if(injectionBegin > 0.0)
timeLoop->setCheckPoint({injectionBegin, injectionEnd});
else
timeLoop->setCheckPoint({injectionEnd});
// the solution vector
Traits::SolutionVector sol;
sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
const auto& cellCenterSol = sol[stokesCellCenterIdx];
const auto& faceSol = sol[stokesFaceIdx];
// apply initial solution for instationary problems
GetPropType<StokesTypeTag, Properties::SolutionVector> stokesSol;
std::get<0>(stokesSol) = cellCenterSol;
std::get<1>(stokesSol) = faceSol;
stokesProblem->applyInitialSolution(stokesSol);
sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
darcyProblem->applyInitialSolution(sol[darcyIdx]);
auto solDarcyOld = sol[darcyIdx];
auto solOld = sol;
couplingManager->init(stokesProblem, darcyProblem, sol);
// the grid variables
using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
stokesGridVariables->init(stokesSol);
using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
darcyGridVariables->init(sol[darcyIdx]);
// intialize the vtk output module
const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
StaggeredVtkOutputModule<StokesGridVariables, GetPropType<StokesTypeTag, Properties::SolutionVector>> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
stokesVtkWriter.write(0.0);
VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
darcyVtkWriter.write(0.0);
// the assembler with time loop for instationary problem
using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
stokesFvGridGeometry->faceFVGridGeometryPtr(),
darcyFvGridGeometry),
std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
stokesGridVariables->faceGridVariablesPtr(),
darcyGridVariables),
couplingManager,
timeLoop);
// the linear solver
using LinearSolver = UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
constexpr auto eps = 1e-6;
// time loop
timeLoop->start(); do
{
// set previous solution for storage evaluations
assembler->setPreviousSolution(solOld);
if(timeLoop->time() > injectionBegin - eps && timeLoop->time() < injectionEnd + eps)
stokesProblem->setInjectionState(true);
else
stokesProblem->setInjectionState(false);
// solve the non-linear system with time step control
nonLinearSolver.solve(sol, *timeLoop);
// make the new solution the old solution
solOld = sol;
stokesGridVariables->advanceTimeStep();
darcyGridVariables->advanceTimeStep();
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
// write vtk output
stokesVtkWriter.write(timeLoop->time());
darcyVtkWriter.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(stokesGridView.comm());
timeLoop->finalize(darcyGridView.comm());
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
////////////////////////////////////////////////////////////
// print dumux end message
if (mpiHelper.rank() == 0)
{
Parameters::print();
DumuxMessage::print(/*firstCall=*/false);
}
return 0;
} // end main
catch (Dumux::ParameterException &e)
{
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
catch (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 (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]
TEnd = 3e6 # s
DtInitial = 100 # s
[Darcy.Grid]
UpperRight = 1.0 1.0
Cells = 20 20
[Stokes.Grid]
LowerLeft = 0.0 1.0
UpperRight = 1.0 2.0
Cells = 20 20
[Stokes.Problem]
Name = stokes
Velocity = 1e-6
Pressure = 1.0e5
InletMoleFraction = 1e-3
InjectionBegin = 0
InjectionEnd = 3e6
[Darcy.Problem]
Name = darcy
Pressure = 1.0e5
InitialMoleFraction = 0.0
[SpatialParams]
AlphaBeaversJoseph = 1.0
Permeability = 1e-10 # m^2
Porosity = 0.3
Tortuosity = 0.5
[Problem]
Name = test_md_darcy1p2c_stokes1p2c_comparison
EnableGravity = false
EnableInertiaTerms = true
[Vtk]
AddVelocity = 1
[Assembly.NumericDifference]
BaseEpsilon = 1e-5
// -*- 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 simple Darcy test problem (cell-centered finite volume method) for the comparison of different diffusion laws.
*/
#ifndef DUMUX_DARCY_SUBPROBLEM_DIFFUSION_COMPARISON_HH
#define DUMUX_DARCY_SUBPROBLEM_DIFFUSION_COMPARISON_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/cellcentered/tpfa/properties.hh>
#include <dumux/flux/maxwellstefanslaw.hh>
#include <dumux/porousmediumflow/1pnc/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include "./../spatialparams.hh"
#include <dumux/material/fluidsystems/1padapter.hh>
#include <dumux/material/fluidsystems/h2oair.hh>
#include <dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh>
namespace Dumux
{
template <class TypeTag>
class DarcySubProblem;
namespace Properties
{
// Create new type tags
namespace TTag {
struct DarcyOnePTwoCTypeTag { using InheritsFrom = std::tuple<OnePNC, CCTpfaModel>; };
} // end namespace TTag
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::DarcyOnePTwoCTypeTag> { using type = Dumux::DarcySubProblem<TypeTag>; };
// The fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::DarcyOnePTwoCTypeTag>
{
using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
static constexpr auto phaseIdx = H2OAir::liquidPhaseIdx; // simulate the water phase
using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
};
// Use moles
template<class TypeTag>
struct UseMoles<TypeTag, TTag::DarcyOnePTwoCTypeTag> { static constexpr bool value = true; };
// Do not replace one equation with a total mass balance
template<class TypeTag>
struct ReplaceCompEqIdx<TypeTag, TTag::DarcyOnePTwoCTypeTag> { static constexpr int value = 3; };
//! Use a model with constant tortuosity for the effective diffusivity
SET_TYPE_PROP(DarcyOnePTwoCTypeTag, EffectiveDiffusivityModel,
DiffusivityConstantTortuosity<GetPropType<TypeTag, Properties::Scalar>>);
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::DarcyOnePTwoCTypeTag> { using type = Dune::YaspGrid<2>; };
// Set the diffusion type
template<class TypeTag>
struct MolecularDiffusionType<TypeTag, TTag::DarcyOnePTwoCTypeTag> { using type = DIFFUSIONTYPE; };
// Set the spatial paramaters type
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::DarcyOnePTwoCTypeTag>
{
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = OnePSpatialParams<FVGridGeometry, Scalar>;
};
}
template <class TypeTag>
class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
// copy some indices for convenience
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
enum {
// grid and world dimension
dim = GridView::dimension,
dimworld = GridView::dimensionworld,
// primary variable indices
conti0EqIdx = Indices::conti0EqIdx,
pressureIdx = Indices::pressureIdx,
};
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
public:
DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
std::shared_ptr<CouplingManager> couplingManager)
: ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
{
pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
initialMoleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.InitialMoleFraction");
}
/*!
* \name Simulation steering
*/
// \{
/*!
* \brief Returns true if a restart file should be written to
* disk.
*/
bool shouldWriteRestartFile() const
{ return false; }
/*!
* \name Problem parameters
*/
// \{
bool shouldWriteOutput() const // define output
{ return true; }
/*!
* \brief Return the temperature within the domain in [K].
*
*/
Scalar temperature() const
{ return 273.15 + 10; } // 10°C
// \}
/*!
* \name Boundary conditions
*/
// \{
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary control volume.
*
* \param element The element
* \param scvf The boundary sub control volume face
*/
BoundaryTypes boundaryTypes(const Element& element, const SubControlVolumeFace& scvf) const
{
BoundaryTypes values;
values.setAllNeumann();
if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
values.setAllCouplingNeumann();
return values;
}
/*!
* \brief Evaluate the boundary conditions for a Neumann control volume.
*
* \param element The element for which the Neumann boundary condition is set
* \param fvGeomentry The fvGeometry
* \param elemVolVars The element volume variables