Hello everybody,
I have a very strange unusual flow distribution after injecting a gas inside my domain (see attachment).
First of all I have a rectangular domain. My boundary conditions are on the left and right side: Dirichlet ( with temperature gradient and hydrostatic pressure). The bottom and top side are considered as Neumann, no flow. My injection is in the middle of horizontal side.
During the injection time everything seems okay, the distribution is symmetric. After the injection my gas starts to migrate to left side. It should actually migrate to the top side. I tried that as non isothermal as well, but I'm still getting the same results.
I will be very thankful, if you offer my any suggestions to stop the migration to left side.
facts:
model: 2pnc /2pncni
fluidsystem: brineco2 (from stabel with hardcoded salinity)
dumux tag: dumux-course-2018
greetings
Anwar
// -*- 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 * \ingroup TwoPNCTest * \brief Problem where hydrogen is injected */ #ifndef DUMUX_TWOPNCNI_CPG_PROBLEM_HH #define DUMUX_TWOPNCNI_CPG_PROBLEM_HH
#include <dumux/discretization/cellcentered/tpfa/properties.hh>
#include <dumux/discretization/box/properties.hh>
#include <dumux/porousmediumflow/2pnc/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/discretization/evalsolution.hh>
#include <dumux/discretization/evalgradients.hh>
//#include <appl/material/fluidsystems/brineco2h2experiment.hh>
#include <dumux/material/fluidsystems/brineco2.hh>
#include "cpg_ccsspatialparams.hh"
#include <appl/2pnc/ccs/co2tables.hh>
//#include <dune/grid/uggrid.hh>
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/maxwellstefanslaw.hh>
#ifndef DIFFUSIONTYPE // default to Fick's law if not set through CMake
#define DIFFUSIONTYPE FicksLaw<TypeTag>
#endif
namespace Dumux {
template <class TypeTag>
class TwoPNCNICPGProblem;
namespace Properties {
//Non-Isothermal
NEW_TYPE_TAG(TwoPNCNICPGTypeTag, INHERITS_FROM(TwoPNC));
NEW_TYPE_TAG(TwoPNCNICPGCCTypeTag, INHERITS_FROM(CCTpfaModel, TwoPNCNICPGTypeTag));
// Set the grid type
//SET_TYPE_PROP(TwoPNCNICPGTypeTag, Grid, Dune::UGGrid<2>);
SET_TYPE_PROP(TwoPNCNICPGTypeTag, Grid, Dune::YaspGrid<2>);
// Set the problem property
SET_TYPE_PROP(TwoPNCNICPGTypeTag, Problem, TwoPNCNICPGProblem<TypeTag>);
//Set fluid configuration
SET_TYPE_PROP(TwoPNCNICPGTypeTag, FluidSystem, FluidSystems::BrineCO2<typename GET_PROP_TYPE(TypeTag, Scalar),
CO2TAB::CO2Tables>);
// Set the spatial parameters
SET_TYPE_PROP(TwoPNCNICPGTypeTag, SpatialParams, TwoPNCNICPGSpatialParams<TypeTag>);
// Define whether mole(true) or mass (false) fractions are used
SET_BOOL_PROP(TwoPNCNICPGTypeTag, UseMoles, true);
//! Here we set FicksLaw or TwoPNCDiffusionsLaw
SET_TYPE_PROP(TwoPNCNICPGTypeTag, MolecularDiffusionType, DIFFUSIONTYPE);
//! Set the default formulation to pw-Sn: This can be over written in the problem.
SET_PROP(TwoPNCNICPGTypeTag, Formulation)
{ static constexpr auto value = TwoPFormulation::p0s1; };
} // end namespace Properties
/*!
* \ingroup TwoPNCTest
* \brief DOC_ME
*/
template <class TypeTag>
class TwoPNCNICPGProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Intersection = typename GET_PROP_TYPE(TypeTag, GridView)::Intersection;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState) ;
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
enum {
// Grid and world dimension
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
// primary variable indices
pressureIdx = Indices::pressureIdx,
switchIdx = Indices::switchIdx,
// FluidSystem CPG BRINECO2H2
CO2Idx = FluidSystem::CO2Idx,
//H2Idx = FluidSystem::H2Idx,
BrineIdx = FluidSystem::BrineIdx,
// FluidSystem H2ON2O2 Test
// CO2Idx = FluidSystem::O2Idx,
// BrineIdx = FluidSystem::H2OIdx,
// H2Idx = FluidSystem::N2Idx,
// temperatureIdx = Indices::temperatureIdx,
// energyEqIdx = Indices::energyEqIdx,
// equation indices
contiCO2EqIdx = Indices::conti0EqIdx + CO2Idx,
//contiH2EqIdx = Indices::conti0EqIdx + H2Idx,
contiBrineEqIdx = Indices::conti0EqIdx + BrineIdx,
// phase presence index
firstPhaseOnly = Indices::firstPhaseOnly,
secondPhaseOnly = Indices::secondPhaseOnly,
bothPhases = Indices::bothPhases,
// wetting/non wetting
liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
gasPhaseIdx = FluidSystem::gasPhaseIdx,
// Well
//isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox)
};
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
//! property that defines whether mole or mass fractions are used
static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
static const bool isBox = GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod == DiscretizationMethod::box;
// world dimension to access gravity vector
//static constexpr int dimWorld = GridView::dimensionworld;
public:
/*!
* \brief The constructor
*
* \param timeManager The time manager
* \param gridView The grid view
*/
TwoPNCNICPGProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{
// initialize the tables of the fluid system
FluidSystem::init();
name_ = getParam<std::string>("Problem.Name");
simulationtime_ = getParam<Scalar>("TimeLoop.TEnd");
sourceinjection_ = getParam<Scalar>("Source.Sourceinjection");
injectionDuration_ = getParam<Scalar>("Source.InjectionDuration");
co2Temperature_ = getParam<Scalar>("Source.CO2Temperature");
co2Pressure_ = getParam<Scalar>("Source.CO2Pressure");
h2Temperature_ = getParam<Scalar>("Source.H2Temperature");
h2Pressure_ = getParam<Scalar>("Source.H2Pressure");
reservoirBottom_ = getParam<Scalar>("Reservoir.ReservoirBottom");
surfacePressure_ = getParam<Scalar>("Surface.SurfacePressure");
surfaceTemperature_ = getParam<Scalar>("Surface.SurfaceTemperature");
h2X_ = getParam<Scalar>("Source.h2X");
permeability_ = getParam<Scalar>("Reservoir.PermeabilityReservoir");
unsigned int codim = GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod == DiscretizationMethod::box ? dim : 0;
temperature_.resize(fvGridGeometry->gridView().size(codim));
//stateing in the console whether mole or mass fractions are used
if(useMoles)
{
std::cout<<"problem uses mole-fractions"<<std::endl;
}
else
{
std::cout<<"problem uses mass-fractions"<<std::endl;
}
}
/*!
* \name Problem parameters
*/
// \{
/*!
* \brief Returns the problem name
*
* This is used as a prefix for files generated by the simulation.
*/
const std::string& name() const
{ return name_; }
/*!
* \brief Returns the temperature [K]
*/
Scalar temperatureAtPos(const GlobalPosition &globalPos) const
{ return initialTemperatureField_(globalPos); }
const std::vector<Scalar>& getTemperature()
{ return temperature_; }
/*!
* \brief Write the relevant secondary variables of the current
* solution into an VTK output file.
*/
void updateVtkOutput(const SolutionVector& curSol)
{
for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
const auto elemSol = elementSolution(element, curSol, this->fvGridGeometry());
auto fvGeometry = localView(this->fvGridGeometry());
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
VolumeVariables volVars;
volVars.update(elemSol, *this, element, scv);
const auto dofIdxGlobal = scv.dofIndex();
temperature_[dofIdxGlobal] = initialTemperatureField_(scv.dofPosition());
//[dofIdxGlobal] = initialTemperatureField_(scv.dofPosition());
}
}
}
/*!
* \brief Append all quantities of interest which can be derived
* from the solution of the current time step to the VTK
* writer.
*/
template<class VTKWriter>
void addFieldsToWriter(VTKWriter& vtk)
{
//const auto numElements = this->fvGridGeometry().gridView().size(0);
//const auto numDofs = this->fvGridGeometry().numDofs();
//moleFracCO2InBrine_.resize(numElements);
//vtk.addField(moleFracCO2InBrine_, "x_CO2w");
vtk.addVolumeVariable([](const VolumeVariables& v){ return v.enthalpy(BrineIdx); }, "enthalpyW");
vtk.addVolumeVariable([](const VolumeVariables& v){ return v.enthalpy(CO2Idx); }, "enthalpyN");
//vtk.addVolumeVariable([](const VolumeVariables& v){ return v.saturation(FS::phase1Idx)*; }, "Sn*")
// //x_g_Co2
//
//
//
//
//
//
// const auto& gridView = this->fvGridGeometry().gridView();
// for (const auto& element : elements(gridView))
// {
//
// auto fvGeometry = localView(this->fvGridGeometry());
// fvGeometry.bindElement(element);
//
// for (const auto& scv : scvs(fvGeometry))
// {
// const auto dofIdxGlobal = scv.dofIndex();
// vtkBoxVolume_[dofIdxGlobal] += scv.volume();
//
// }
// //
// //
// }
// const auto& gridView = this->fvGridGeometry().gridView();
// for (const auto& element : elements(gridView))
// {
// const auto eIdx = this->fvGridGeometry().elementMapper().index(element);
//
// moleFracCO2InBrine_[eIdx] = ;
// }
}
/*!
* \name Boundary conditions
*/
// \{
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
* \param values The boundary types for the conservation equations
* \param globalPos The position for which the bc type should be evaluated
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
//left side
if(globalPos[0] < eps_ )
{
values.setAllDirichlet();
}
//right side
else if (globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_)
{
values.setAllDirichlet();
}
//top & bottom side
else
{
values.setAllNeumann();
}
return values;
}
/*!
* \brief Evaluates the boundary conditions for a Dirichlet
* boundary segment
*
* \param values Stores the Dirichlet values for the conservation equations in
* \f$ [ \textnormal{unit of primary variable} ] \f$
* \param globalPos The global position
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
return initial_(globalPos);
}
/*!
* \brief Evaluate the boundary conditions for a Neumann
* boundary segment.
*
* For this method, the \a priVars parameter stores the mass flux
* in normal direction of each component. Negative values mean
* influx.
*
* The units must be according to either using mole or mass fractions. (mole/(m^2*s) or kg/(m^2*s))
*/
PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
{
// // initialize values to zero, i.e. no-flow Neumann boundary conditions
PrimaryVariables values(0.0);
//
//
// if (time_ < injectionDuration_ && globalPos[0] > 685 +eps_ && globalPos[0] <715 && globalPos[1] > eps_)
// {
// // Aus dumux course
// // inject nitrogen. negative values mean injection
// // units kg/(s*m^2)
// values[contiCO2EqIdx] = sourceinjection_;
// //values[contiH2EqIdx] = h2X_*sourceinjection_;
// //to be edited
// values[energyEqIdx] = 20;
// std::cout<< "Marker Injektion" << "\n";
// }
//
return values;
}
NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
{
NumEqVector values(0.0);
if (time_ < injectionDuration_ && globalPos[0] > 685 && globalPos[0] <715 && globalPos[1] < 200 && globalPos[1] > 150 )
{
values[contiCO2EqIdx]= 0.0015;
// values[contiH2EqIdx]= values[contiCO2EqIdx]*0.01;
//values[energyEqIdx]= 20;
//values[contiH2EqIdx] * FluidSystem::molarMass(H2Idx) * Components::H2<Scalar>::gasEnthalpy(h2Temperature_, h2Pressure_);
}
return values;
}
/*!
* \name Volume terms
*/
// \{
/*!
* \brief Evaluates the initial values for a control volume
*
* \param values Stores the initial values for the conservation equations in
* \f$ [ \textnormal{unit of primary variables} ] \f$
* \param globalPos The global position
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
return initial_(globalPos);
}
//! set the time for the time dependent boundary conditions (called from main)
void setTime(Scalar time)
{ time_ = time; }
private:
/*!
* \brief Evaluates the initial values for a control volume
*
* The internal method for the initial condition
*
* \param values Stores the initial values for the conservation equations in
* \f$ [ \textnormal{unit of primary variables} ] \f$
* \param globalPos The global position
*/
PrimaryVariables initial_(const GlobalPosition &globalPos) const
{
PrimaryVariables priVars(0.0);
priVars.setState(Indices::firstPhaseOnly);
const Scalar densityW = 1000;
const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(reservoirBottom_ - globalPos[dimWorld-1]);
priVars[pressureIdx] = pw;
//mole-fraction formulation
priVars[Indices::switchIdx] = 1e-8;
priVars[CO2Idx] = 1e-8;
//non-isothermal
//priVars[temperatureIdx] = initialTemperatureField_(globalPos);
return priVars;
}
Scalar initialTemperatureField_(const GlobalPosition globalPos) const
{
return surfaceTemperature_ + (reservoirBottom_ - globalPos[1])*0.03;
}
Scalar initialPressureField_(const GlobalPosition globalPos) const
{
const Scalar densityW = 1000;
return surfacePressure_ - densityW * 9.81 * (reservoirBottom_ - globalPos[1]);
}
static constexpr Scalar eps_ = 1e-6;
Scalar sourceinjection_;
Scalar h2X_;
Scalar time_;
Scalar injectionDuration_;
Scalar reservoirBottom_;
Scalar surfacePressure_;
Scalar surfaceTemperature_;
Scalar injectionPressure_ ;
Scalar co2Temperature_;
Scalar co2Pressure_;
Scalar h2Temperature_;
Scalar h2Pressure_;
Scalar simulationtime_;
Scalar permeability_;
Scalar moleFracCO2InBrine_;
Scalar freqOutput_;
std::string name_ ;
std::vector<Scalar> Kxx_;
std::vector<double> temperature_;
std::vector<Scalar> vtkKxx_, vtkPorosity_, vtkBoxVolume_;
};
} //end namespace Dumux
#endif
2pnc_cpg_ccs_problem_experimental.input
Description: Binary data
// -*- 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 * \ingroup TwoPNCTest * \brief Definition of the spatial parameters for the fuel cell * problem which uses the isothermal/non-insothermal 2pnc box model */ #ifndef DUMUX_TWOPNCNI_CPG_SPATIAL_PARAMS_HH #define DUMUX_TWOPNCNI_CPG_SPATIAL_PARAMS_HH #include <dumux/material/spatialparams/fv.hh> #include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh> #include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh> #include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh> #include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh> #include <dumux/io/gnuplotinterface.hh> #include <dumux/io/ploteffectivediffusivitymodel.hh> #include <dumux/io/plotmateriallaw.hh> #include <dumux/io/plotthermalconductivitymodel.hh> namespace Dumux { /*! * \ingroup TwoPNCTest * \brief Definition of the spatial parameters for the TwoPNCDiffusion * problem which uses the isothermal 2p2c box model */ template<class TypeTag> class TwoPNCNICPGSpatialParams : public FVSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry), typename GET_PROP_TYPE(TypeTag, Scalar), TwoPNCNICPGSpatialParams<TypeTag>> { using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); using GridView = typename FVGridGeometry::GridView; using ParentType = FVSpatialParams<FVGridGeometry, Scalar, TwoPNCNICPGSpatialParams<TypeTag>>; static constexpr int dimWorld = GridView::dimensionworld; using Element = typename GridView::template Codim<0>::Entity; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>; using EffectiveLaw = RegularizedBrooksCorey<Scalar>; public: using PermeabilityType = DimWorldMatrix; using MaterialLaw = EffToAbsLaw<EffectiveLaw>; using MaterialLawParams = typename MaterialLaw::Params; /*! * \brief The constructor * * \param gridView The grid view */ TwoPNCNICPGSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry) : ParentType(fvGridGeometry), K_(0) { // intrinsic permeabilities permeability_ = getParam<Scalar>("Reservoir.PermeabilityReservoir"); // porosities porosity_ = getParam<Scalar>("Reservoir.PorosityReservoir"); // Rock solidDensity_ = getParam<Scalar>("Component.SolidDensity"); solidThermalConductivity_ =getParam<Scalar>("Component.SolidThermalConductivity"); solidHeatCapacity_ =getParam<Scalar>("Component.SolidHeatCapacity"); // residual saturations materialParams_.setSwr(0.02); materialParams_.setSnr(0.0); //parameters for the vanGenuchten law materialParams_.setPe(1e2); // alpha = 1/pcb materialParams_.setLambda(1.3); } /*! * \brief Returns the hydraulic conductivity \f$[m^2]\f$ * * \param globalPos The global position */ DimWorldMatrix permeabilityAtPos(const GlobalPosition& globalPos) const { return permeability_; } /*! * \brief Define the porosity \f$[-]\f$ of the spatial parameters * * \param globalPos The global position */ Scalar porosityAtPos(const GlobalPosition& globalPos) const { return porosity_ ; } /*! * \brief return the parameter object for the Brooks-Corey material law which depends on the position * * \param globalPos The global position */ const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const { return materialParams_; } /*! * \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix. * * This is only required for non-isothermal models. * * \param globalPos The global position */ Scalar solidHeatCapacityAtPos(const GlobalPosition& globalPos) const { return solidHeatCapacity_; // specific heat capacity of granite [J / (kg K)] } /*! * \brief Returns the mass density \f$[kg / m^3]\f$ of the rock matrix. * * This is only required for non-isothermal models. * * \param globalPos The global position */ Scalar solidDensityAtPos(const GlobalPosition& globalPos) const { return solidDensity_; // density of granite [kg/m^3] -> Like M.Saar } /*! * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material. * * \param globalPos The global position */ Scalar solidThermalConductivityAtPos(const GlobalPosition& globalPos) const { return solidThermalConductivity_; } /*! * \brief Function for defining which phase is to be considered as the wetting phase. * * \return the wetting phase index * \param globalPos The position of the center of the element */ template<class FluidSystem> int wettingPhaseAtPos(const GlobalPosition& globalPos) const { return FluidSystem::BrineIdx; } // { return FluidSystem::H2OIdx; } private: DimWorldMatrix K_; Scalar permeability_ ; Scalar porosity_; Scalar solidDensity_ , solidThermalConductivity_ ,solidHeatCapacity_ ; static constexpr Scalar eps_ = 1e-6; MaterialLawParams materialParams_; }; }//end namespace #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 Test for the 2pnc cc model used for water management in PEM fuel cells. */ #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 "2pnc_cpg_ccs_problem_experimental.hh" #include <dumux/common/properties.hh> #include <dumux/common/parameters.hh> #include <dumux/common/valgrind.hh> #include <dumux/common/dumuxmessage.hh> #include <dumux/common/defaultusagemessage.hh> #include <dumux/linear/amgbackend.hh> #include <dumux/nonlinear/privarswitchnewtonsolver.hh> #include <dumux/assembly/fvassembler.hh> #include <dumux/assembly/diffmethod.hh> #include <dumux/discretization/methods.hh> #include <dumux/io/grid/gridmanager.hh> #include <dumux/io/vtkoutputmodule.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 options for this program is:\n" "\t-ParameterFile Parameter file (Input file) \n"; std::cout << errorMessageOut << "\n"; } } int main(int argc, char** argv) try { using namespace Dumux; // define the type tag for this problem using TypeTag = TTAG(TwoPNCNICPGCCTypeTag); // 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<typename GET_PROP_TYPE(TypeTag, 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 FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView); fvGridGeometry->update(); // the problem (initial and boundary conditions) using Problem = typename GET_PROP_TYPE(TypeTag, Problem); auto problem = std::make_shared<Problem>(fvGridGeometry); // the solution vector using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); SolutionVector x(fvGridGeometry->numDofs()); problem->applyInitialSolution(x); auto xOld = x; // the grid variables using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables); auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry); gridVariables->init(x, xOld); // get some time loop parameters using Scalar = typename GET_PROP_TYPE(TypeTag, 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"); // initialize the vtk output module using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields); VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name()); VtkOutputFields::init(vtkWriter); //!< Add model specific output fields //problem->addVtkFields(vtkWriter); //!< Add problem specific output fields vtkWriter.write(0.0); // output every vtkOutputInterval time step const auto vtkOutputInterval = getParam<int>("Problem.OutputInterval"); // instantiate time loop auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, 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, fvGridGeometry, gridVariables, timeLoop); // the linear solver using LinearSolver = AMGBackend<TypeTag>; auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper()); // the non-linear solver using NewtonSolver = PriVarSwitchNewtonSolver<Assembler, LinearSolver, typename GET_PROP_TYPE(TypeTag, PrimaryVariableSwitch)>; NewtonSolver nonLinearSolver(assembler, linearSolver); //the convergence writer using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using NewtonConvergenceWriter = Dumux::NewtonConvergenceWriter<GridView, SolutionVector>; auto convergenceWriter = std::make_shared<NewtonConvergenceWriter>(leafGridView, fvGridGeometry->numDofs()); // time loop timeLoop->start(); do { // set previous solution for storage evaluations assembler->setPreviousSolution(xOld); // solve the non-linear system with time step control nonLinearSolver.solve(x, *timeLoop); //set time in problem (is used in time-dependent Neumann boundary condition) problem->setTime(timeLoop->time()+timeLoop->timeStepSize()); // make the new solution the old solution xOld = x; gridVariables->advanceTimeStep(); // advance to the time loop to the next step timeLoop->advanceTimeStep(); // update the output fields before write problem->updateVtkOutput(x); // report statistics of this time step timeLoop->reportTimeStep(); // set new dt as suggested by the newton solver timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize())); // write vtk output if(timeLoop->timeStepIndex()==0 || timeLoop->timeStepIndex() % vtkOutputInterval == 0 || timeLoop->willBeFinished()) vtkWriter.write(timeLoop->time()); } 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; }
_______________________________________________ Dumux mailing list [email protected] https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
