Commit 8f6b7b3d authored by Katharina Heck's avatar Katharina Heck Committed by Kilian Weishaupt
Browse files

[test][multidomain] add 1p3c_1p3c test for maxwellstefan coupling

parent 969732f0
add_input_file_links()
dune_add_test(NAME test_md_darcy1p3c_stokes1p3c_horizontal
SOURCES main.cc
CMAKE_GUARD HAVE_UMFPACK
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_md_darcy1p3c_stokes1p3c_horizontal_stokes-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_md_darcy1p3c_stokes1p3c_horizontal_stokes-00015.vtu
${CMAKE_SOURCE_DIR}/test/references/test_md_darcy1p3c_stokes1p3c_horizontal_darcy-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_md_darcy1p3c_stokes1p3c_horizontal_darcy-00015.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_md_darcy1p3c_stokes1p3c_horizontal 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 fluid system for one phase with the components h2, n2 and co2
*/
#ifndef DUMUX_THREE_GAS_COMPONENT_FLUID_SYSTEM_HH
#define DUMUX_THREE_GAS_COMPONENT_FLUID_SYSTEM_HH
#include <dumux/material/fluidsystems/base.hh>
namespace Dumux {
namespace FluidSystems {
/*!
* \ingroup FluidSystems
* \brief A simple fluid system with one MaxwellStefan component.
* \todo doc me!
*/
template<class Scalar>
class H2N2CO2FluidSystem: public Base<Scalar, H2N2CO2FluidSystem<Scalar>>
{
using ThisType = H2N2CO2FluidSystem<Scalar>;
using Base = FluidSystems::Base<Scalar, ThisType>;
public:
//! The number of phases
static constexpr int numPhases = 1;
static constexpr int numComponents = 3;
static constexpr int H2Idx = 0;//first major component
static constexpr int N2Idx = 1;//second major component
static constexpr int CO2Idx = 2;//secondary component
//! Human readable component name (index compIdx) (for vtk output)
static std::string componentName(int compIdx)
{
switch (compIdx)
{
case H2Idx: return "H2";
case N2Idx: return "N2";
case CO2Idx: return "CO2";
}
DUNE_THROW(Dune::InvalidStateException, "Invalid compIdx index " << compIdx);
}
//! Human readable phase name (index phaseIdx) (for velocity vtk output)
static std::string phaseName(int phaseIdx = 0)
{ return "Gas"; }
//! Molar mass in kg/mol of the component with index compIdx
static Scalar molarMass(unsigned int compIdx)
{ return 0.02896; }
using Base::binaryDiffusionCoefficient;
/*!
* \brief Given a phase's composition, temperature and pressure,
* return the binary diffusion coefficient \f$\mathrm{[m^2/s]}\f$ for components
* \f$i\f$ and \f$j\f$ in this phase.
*
* \param fluidState An arbitrary fluid state
* \param phaseIdx The index of the fluid phase to consider
* \param compIIdx The index of the first component to consider
* \param compJIdx The index of the second component to consider
*/
template <class FluidState>
static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
int phaseIdx,
int compIIdx,
int compJIdx)
{
if (compIIdx > compJIdx)
{
using std::swap;
swap(compIIdx, compJIdx);
}
if (compIIdx == H2Idx && compJIdx == N2Idx)
return 83.3e-6;
if (compIIdx == H2Idx && compJIdx == CO2Idx)
return 68.0e-6;
if (compIIdx == N2Idx && compJIdx == CO2Idx)
return 16.8e-6;
DUNE_THROW(Dune::InvalidStateException,
"Binary diffusion coefficient of components "
<< compIIdx << " and " << compJIdx << " is undefined!\n");
}
using Base::density;
/*!
* \brief Given a phase's composition, temperature, pressure, and
* the partial pressures of all components, return its
* density \f$\mathrm{[kg/m^3]}\f$.
* \param phaseIdx index of the phase
* \param fluidState the fluid state
*
*/
template <class FluidState>
static Scalar density(const FluidState &fluidState,
const int phaseIdx)
{
return 1;
}
using Base::viscosity;
/*!
* \brief Calculate the dynamic viscosity of a fluid phase \f$\mathrm{[Pa*s]}\f$
* \param fluidState An arbitrary fluid state
* \param phaseIdx The index of the fluid phase to consider
*/
template <class FluidState>
static Scalar viscosity(const FluidState &fluidState,
int phaseIdx)
{
return 1e-6;
}
using Base::molarDensity;
/*!
* \brief The molar density \f$\rho_{mol,\alpha}\f$
* of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
*
* The molar density for the simple relation is defined by the
* mass density \f$\rho_\alpha\f$ and the molar mass of the main component \f$M_\kappa\f$:
*
* \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{M_\kappa} \;.\f]
*/
template <class FluidState>
static Scalar molarDensity(const FluidState &fluidState, int phaseIdx)
{
return density(fluidState, phaseIdx)/molarMass(0);
}
};
} // end namespace FluidSystems
} // end namespace Dumux
#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 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 (1p)
*/
#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::StokesOnePThreeCTypeTag>
{
using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOnePThreeCTypeTag>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::DarcyOnePThreeCTypeTag>
{
using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesOnePThreeCTypeTag, Properties::TTag::StokesOnePThreeCTypeTag, 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::StokesOnePThreeCTypeTag;
using DarcyTypeTag = Properties::TTag::DarcyOnePThreeCTypeTag;
// 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);
// 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<TimeLoop<Scalar>>(restartTime, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
// 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);
// time loop
timeLoop->start(); do
{
// set previous solution for storage evaluations
assembler->setPreviousSolution(solOld);
// 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 = 1e6 # s
DtInitial = 1000 # 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
[Darcy.Problem]
Name = darcy
Pressure = 1.0e5
InitialMoleFractionN2 = 0.2
InitialMoleFractionCO2 = 0.3
[SpatialParams]
AlphaBeaversJoseph = 1.0
Permeability = 1e-10 # m^2
Porosity = 0.3
Tortuosity = 0.5
[Problem]
Name = test_md_darcy1p3c_stokes1p3c_horizontal
EnableGravity = false
[Vtk]
AddVelocity = 1
[Newton]
MaxRelativeShift = 1e-10
[Assembly]
NumericDifference.BaseEpsilon = 1e-8
// -*- 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 Darcy test problem using Maxwell-Stefan diffusion.
*/
#ifndef DUMUX_DARCY_SUBPROBLEM_ONEPTHREEC_HH
#define DUMUX_DARCY_SUBPROBLEM_ONEPTHREEC_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 "../1p2c_1p2c/spatialparams.hh"
#include "h2n2co2fluidsystem.hh"
#include <dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh>
namespace Dumux
{
template <class TypeTag>
class DarcySubProblem;
namespace Properties
{
// Create new type tags
namespace TTag {
struct DarcyOnePThreeCTypeTag { using InheritsFrom = std::tuple<OnePNC, CCTpfaModel>; };
} // end namespace TTag
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::DarcyOnePThreeCTypeTag> { using type = Dumux::DarcySubProblem<TypeTag>; };
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::DarcyOnePThreeCTypeTag> { using type = FluidSystems::H2N2CO2FluidSystem<GetPropType<TypeTag, Properties::Scalar>>; };
// Use moles
template<class TypeTag>
struct UseMoles<TypeTag, TTag::DarcyOnePThreeCTypeTag> { static constexpr bool value = true; };
// Do not replace one equation with a total mass balance
template<class TypeTag>
struct ReplaceCompEqIdx<TypeTag, TTag::DarcyOnePThreeCTypeTag> { static constexpr int value = 3; };
//! Use a model with constant tortuosity for the effective diffusivity
SET_TYPE_PROP(DarcyOnePThreeCTypeTag, EffectiveDiffusivityModel,
DiffusivityConstantTortuosity<GetPropType<TypeTag, Properties::Scalar>>);
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::DarcyOnePThreeCTypeTag> { using type = Dune::YaspGrid<2>; };
//Set the diffusion type
template<class TypeTag>
struct MolecularDiffusionType<TypeTag, TTag::DarcyOnePThreeCTypeTag> { using type = MaxwellStefansLaw<TypeTag>; };
// Set the spatial paramaters type
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::DarcyOnePThreeCTypeTag>
{
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>;