Commit ef97a4d5 authored by Timo Koch's avatar Timo Koch Committed by Ned Coltman
Browse files

[test] Remove Maxwell-Stefan tracer test

The idea of the tracer model is really tied to Fickian diffusion and components
with small concentration diffusing in a bulk phase. Hence also it's name.
The Maxwell-Stefan test actual models 3 gaseous components where none of
the components is really the main phase but most importantly there is no
bulk phase with a main component that is not considered in the equations.
Thus the 1pnc model is the right choice for such a test.
This changes the test to derive from the 1pnc model. As for the given
setting with only Neumann boundaries the 1pnc wouldn't convergence for
a compressible fluid, the Ideal gas law is used to model a
compressible gas.

The same test actually exists for the 1pnc model already in 1pnc/implicit/1p3c.
Thus this commit simply removes the tracer test.
parent 29f9f790
add_subdirectory(2ptracer)
add_subdirectory(constvel)
add_subdirectory(multicomp)
add_input_file_links()
# tracer tests
dumux_add_test(NAME test_tracer_maxwellstefan_tpfa
LABELS porousmediumflow tracer
SOURCES main.cc
COMPILE_DEFINITIONS TYPETAG=MaxwellStefanTestCC
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_tracer_maxwellstefan_tpfa-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_tracer_maxwellstefan_tpfa-00026.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_tracer_maxwellstefan_tpfa params.input -Problem.Name test_tracer_maxwellstefan_tpfa")
dumux_add_test(NAME test_tracer_maxwellstefan_box
LABELS porousmediumflow tracer
SOURCES main.cc
COMPILE_DEFINITIONS TYPETAG=MaxwellStefanTestBox
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_tracer_maxwellstefan_box-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_tracer_maxwellstefan_box-00026.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_tracer_maxwellstefan_box params.input -Problem.Name test_tracer_maxwellstefan_box")
// -*- 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 multi-component tracer porous medium flow 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.hh>
#include <dune/istl/io.hh>
#include "problem.hh"
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/valgrind.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/linear/amgbackend.hh>
#include <dumux/linear/linearsolvertraits.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/discretization/method.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
/*!
* \brief Provides an interface for customizing error messages associated with
* reading in parameters.
*
* \param progName The name of the program, that was tried to be started.
* \param errorMsg The error message that was issued by the start function.
* Comprises the thing that went wrong and a general help message.
*/
void usage(const char *progName, const std::string &errorMsg)
{
if (errorMsg.size() > 0) {
std::string errorMessageOut = "\nUsage: ";
errorMessageOut += progName;
errorMessageOut += " [options]\n";
errorMessageOut += errorMsg;
errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n"
"\t-TimeManager.TEnd End of the simulation [s] \n"
"\t-TimeManager.DtInitial Initial timestep size [s] \n"
"\t-Grid.LowerLeft Lower left corner coordinates\n"
"\t-Grid.UpperRight Upper right corner coordinates\n"
"\t-Grid.Cells Number of cells in respective coordinate directions\n"
"\t definition in DGF format\n"
"\t-SpatialParams.LensLowerLeft coordinates of the lower left corner of the lens [m] \n"
"\t-SpatialParams.LensUpperRight coordinates of the upper right corner of the lens [m] \n"
"\t-SpatialParams.Permeability Permeability of the domain [m^2] \n"
"\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n";
std::cout << errorMessageOut
<< "\n";
}
}
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 command line arguments and input file
Parameters::init(argc, argv, usage);
// try to create a grid (from the given grid file or the input file)
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
gridManager.init();
////////////////////////////////////////////////////////////
// run 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);
// 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>;
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
IOFields::initOutputModule(vtkWriter); //! Add model specific output fields
vtkWriter.write(0.0);
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(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 = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
// the non-linear solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
// time loop
timeLoop->start(); do
{
// solve the non-linear system with time step control
nonLinearSolver.solve(x, *timeLoop);
// make the new solution the old solution
xOld = x;
gridVariables->advanceTimeStep();
problem->plotComponentsOverTime(x, timeLoop->time() + timeLoop->timeStepSize());
// 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 the 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;
} // 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]
DtInitial = 1 # [s]
TEnd = 70000 # [s]
[Grid]
UpperRight = 1 1
Cells = 30 30
[Problem]
Name = maxwellstefan # name passed to the output routines
[Vtk]
AddVelocity = true
[Newton]
MaxRelativeShift = 1e-12
// -*- 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 TracerTest
* \brief Definition of a problem for the MaxwellStefan problem:
* A rotating velocity field mixes a MaxwellStefan band in a porous groundwater reservoir.
*/
#ifndef DUMUX_MAXWELL_STEFAN_TEST_PROBLEM_HH
#define DUMUX_MAXWELL_STEFAN_TEST_PROBLEM_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/elementsolution.hh>
#include <dumux/discretization/box.hh>
#include <dumux/discretization/cctpfa.hh>
#include <dumux/porousmediumflow/tracer/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include "spatialparams.hh"
#include <dumux/flux/maxwellstefanslaw.hh>
#include <dumux/io/gnuplotinterface.hh>
#include <dumux/material/fluidsystems/base.hh>
namespace Dumux {
template <class TypeTag>
class MaxwellStefanTestProblem;
namespace Properties {
// Create new type tags
namespace TTag {
struct MaxwellStefanTest { using InheritsFrom = std::tuple<Tracer>; };
struct MaxwellStefanTestCC { using InheritsFrom = std::tuple<MaxwellStefanTest, CCTpfaModel>; };
struct MaxwellStefanTestBox { using InheritsFrom = std::tuple<MaxwellStefanTest, BoxModel>; };
} // end namespace TTag
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::MaxwellStefanTest> { using type = Dune::YaspGrid<2>; };
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::MaxwellStefanTest> { using type = MaxwellStefanTestProblem<TypeTag>; };
// Set the spatial parameters
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::MaxwellStefanTest>
{
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = MaxwellStefanTestSpatialParams<GridGeometry, Scalar>;
};
// Define whether mole(true) or mass (false) fractions are used
template<class TypeTag>
struct UseMoles<TypeTag, TTag::MaxwellStefanTest> { static constexpr bool value = true; };
//! Here we set FicksLaw or MaxwellStefansLaw
template<class TypeTag>
struct MolecularDiffusionType<TypeTag, TTag::MaxwellStefanTest> { using type = MaxwellStefansLaw<TypeTag>; };
//! A simple fluid system with one MaxwellStefan component
template<class TypeTag>
class MaxwellStefanTracerFluidSystem
: public FluidSystems::Base<GetPropType<TypeTag, Properties::Scalar>, MaxwellStefanTracerFluidSystem<TypeTag>>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Element = typename GridView::template Codim<0>::Entity;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
public:
static constexpr bool isTracerFluidSystem()
{ return true; }
//! The number of components
static constexpr int numComponents = 3;
static constexpr int compOneIdx = 0;//first major component
static constexpr int compTwoIdx = 1;//second major component
static constexpr int compThreeIdx = 2;//secondary component
//! Human readable component name (index compIdx) (for vtk output)
static std::string componentName(int compIdx)
{
switch (compIdx)
{
case compOneIdx: return "CompOne";
case compTwoIdx: return "CompTwo";
case compThreeIdx:return "CompThree";
}
DUNE_THROW(Dune::InvalidStateException, "Invalid component 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; /*air*/ }
//! Binary diffusion coefficient
//! (might depend on spatial parameters like pressure / temperature)
static Scalar binaryDiffusionCoefficient(unsigned int compIdx,
const Problem& problem,
const Element& element,
const SubControlVolume& scv)
{
if (compIdx == compOneIdx)
return 0;
if (compIdx == compTwoIdx)
return 83.3e-6;
if (compIdx == compThreeIdx)
return 68.0e-6;
DUNE_THROW(Dune::InvalidStateException,
"Binary diffusion coefficient of component "
<< compIdx <<" is undefined!\n");
}
//! Binary diffusion coefficient
//! (might depend on spatial parameters like pressure / temperature)
static Scalar binaryDiffusionCoefficient(unsigned int compIIdx,
unsigned int compJIdx,
const Problem& problem,
const Element& element,
const SubControlVolume& scv)
{
if (compIIdx > compJIdx)
{
using std::swap;
swap(compIIdx, compJIdx);
}
if (compIIdx == compOneIdx && compJIdx == compTwoIdx)
return 83.3e-6;
if (compIIdx == compOneIdx && compJIdx == compThreeIdx)
return 68.0e-6;
if (compIIdx == compTwoIdx && compJIdx == compThreeIdx)
return 16.8e-6;
DUNE_THROW(Dune::InvalidStateException,
"Binary diffusion coefficient of components "
<< compIIdx << " and " << compJIdx << " is undefined!\n");
}
/*!
* \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);
}
};
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::MaxwellStefanTest> { using type = MaxwellStefanTracerFluidSystem<TypeTag>; };
} // end namespace Properties
/*!
* \ingroup TracerTest
* \brief Definition of a problem for the MaxwellStefan problem.
*
* This problem uses the MaxwellStefan equations.
*
* To run the simulation execute the following line in shell:
* <tt>./test_boxmaxwellstefan -ParameterFile ./test_boxmaxwellstefan.input</tt> or
* <tt>./test_ccmaxwellstefan -ParameterFile ./test_ccMaxwellstefan.input</tt>
*/
template <class TypeTag>
class MaxwellStefanTestProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
//! property that defines whether mole or mass fractions are used
static constexpr bool useMoles = getPropValue<TypeTag, Properties::UseMoles>();
using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
MaxwellStefanTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
name_ = getParam<std::string>("Problem.Name");
// 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';
plotOutput_ = false;
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \brief The problem name.
* This is used as a prefix for files generated by the simulation.
*/
const std::string& name() const
{ return name_; }
//! Called after every time step
void plotComponentsOverTime(const SolutionVector& curSol, Scalar time)
{
if (plotOutput_)
{
Scalar x_CompThree_left = 0.0;
Scalar x_CompTwo_left = 0.0;
Scalar x_CompThree_right = 0.0;
Scalar x_CompTwo_right = 0.0;
Scalar x_CompOne_left = 0.0;
Scalar x_CompOne_right = 0.0;
Scalar i = 0.0;
Scalar j = 0.0;
if (!(time < 0.0))
{
for (const auto& element : elements(this->gridGeometry().gridView()))
{
auto fvGeometry = localView(this->gridGeometry());
fvGeometry.bindElement(element);
const auto elemSol = elementSolution(element, curSol, this->gridGeometry());
for (auto&& scv : scvs(fvGeometry))
{
const auto& globalPos = scv.dofPosition();
VolumeVariables volVars;
volVars.update(elemSol, *this, element, scv);
if (globalPos[0] < 0.5)
{
x_CompThree_left += volVars.moleFraction(0,2);
x_CompTwo_left += volVars.moleFraction(0,1);
x_CompOne_left += volVars.moleFraction(0,0);
i +=1;
}
else
{
x_CompThree_right += volVars.moleFraction(0,2);
x_CompTwo_right += volVars.moleFraction(0,1);
x_CompOne_right += volVars.moleFraction(0,0);
j +=1;
}
}
}
x_CompThree_left /= i;
x_CompTwo_left /= i;
x_CompOne_left /= i;
x_CompThree_right /= j;
x_CompTwo_right /= j;
x_CompOne_right /= j;
//do a gnuplot
x_.push_back(time); // in seconds
y_.push_back(x_CompTwo_left);
y2_.push_back(x_CompTwo_right);
y3_.push_back(x_CompThree_left);
y4_.push_back(x_CompThree_right);
y5_.push_back(x_CompOne_left);
y6_.push_back(x_CompOne_right);
gnuplot_.resetPlot();
gnuplot_.setXRange(0, std::min(time, 72000.0));
gnuplot_.setYRange(0.4, 0.6);
gnuplot_.setXlabel("time [s]");
gnuplot_.setYlabel("mole fraction mol/mol");
gnuplot_.addDataSetToPlot(x_, y_, "CompTwo_left.dat", "w l t 'CompTwo left'");
gnuplot_.addDataSetToPlot(x_, y2_, "CompTwo_right.dat", "w l t 'CompTwo right'");
gnuplot_.plot("mole_fraction_CompTwo");
gnuplot2_.resetPlot();
gnuplot2_.setXRange(0, std::min(time, 72000.0));
gnuplot2_.setYRange(0.0, 0.6);