diff --git a/dumux/porousmediumflow/3p/implicit/model.hh b/dumux/porousmediumflow/3p/implicit/model.hh index 2178999bcc702912f9dacbf7b1e70b1737b0acdc..b4920149a98f712a3e994641826a1195dbbd625e 100644 --- a/dumux/porousmediumflow/3p/implicit/model.hh +++ b/dumux/porousmediumflow/3p/implicit/model.hh @@ -27,12 +27,8 @@ #ifndef DUMUX_3P_MODEL_HH #define DUMUX_3P_MODEL_HH -#include -#include "properties.hh" - -namespace Dumux -{ /*! + * \file * \ingroup ThreePModel * \brief Adaption of the fully implicit scheme to the three-phase flow model. * @@ -57,56 +53,7 @@ namespace Dumux * The used primary variables are gas phase pressure \f$p_g\f$, * water saturation \f$S_w\f$ and NAPL saturation \f$S_n\f$. */ -template -class ThreePModel: public GET_PROP_TYPE(TypeTag, BaseModel) -{ - using ParentType = typename GET_PROP_TYPE(TypeTag, BaseModel); - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); - using Indices = typename GET_PROP_TYPE(TypeTag, Indices);; - - enum { - wPhaseIdx = Indices::wPhaseIdx, - nPhaseIdx = Indices::nPhaseIdx, - gPhaseIdx = Indices::gPhaseIdx, - }; - -public: - /*! - * \brief Apply the initial conditions to the model. - * - * \param problem The object representing the problem which needs to - * be simulated. - */ - void init(Problem& problem) - { - ParentType::init(problem); - - // register standardized vtk output fields - auto& vtkOutputModule = problem.vtkOutputModule(); - vtkOutputModule.addSecondaryVariable("sw", [](const VolumeVariables& v){ return v.saturation(wPhaseIdx); }); - vtkOutputModule.addSecondaryVariable("sn", [](const VolumeVariables& v){ return v.saturation(nPhaseIdx); }); - vtkOutputModule.addSecondaryVariable("sg", [](const VolumeVariables& v){ return v.saturation(gPhaseIdx); }); - vtkOutputModule.addSecondaryVariable("pw", [](const VolumeVariables& v){ return v.pressure(wPhaseIdx); }); - vtkOutputModule.addSecondaryVariable("pn", [](const VolumeVariables& v){ return v.pressure(nPhaseIdx); }); - vtkOutputModule.addSecondaryVariable("pg", [](const VolumeVariables& v){ return v.pressure(gPhaseIdx); }); - vtkOutputModule.addSecondaryVariable("rhow", [](const VolumeVariables& v){ return v.density(wPhaseIdx); }); - vtkOutputModule.addSecondaryVariable("rhon", [](const VolumeVariables& v){ return v.density(nPhaseIdx); }); - vtkOutputModule.addSecondaryVariable("rhog", [](const VolumeVariables& v){ return v.density(gPhaseIdx); }); - vtkOutputModule.addSecondaryVariable("porosity", [](const VolumeVariables& v){ return v.porosity(); }); - vtkOutputModule.addSecondaryVariable("permeability", [](const VolumeVariables& v){ return v.permeability(); }); - vtkOutputModule.addSecondaryVariable("temperature", [](const VolumeVariables& v){ return v.temperature(); }); -// vtkOutputModule.addSecondaryVariable("mobW", [](const VolumeVariables& v){ return v.mobility(wPhaseIdx); }); -// vtkOutputModule.addSecondaryVariable("mobN", [](const VolumeVariables& v){ return v.mobility(nPhaseIdx); }); -// vtkOutputModule.addSecondaryVariable("mobG", [](const VolumeVariables& v){ return v.mobility(gPhaseIdx); }); - // TODO: these lines are commented-out in order to comply with the "old" reference solution; - // can be changed some time, as well as the parameter names - } -}; - -} - -#include "propertydefaults.hh" +#include "properties.hh" #endif diff --git a/dumux/porousmediumflow/3p/implicit/properties.hh b/dumux/porousmediumflow/3p/implicit/properties.hh index fec41bc369f896592305ca3187785c1689fdcd46..faf1cd79f5cb07c580e5edb60ac30cd385d30efa 100644 --- a/dumux/porousmediumflow/3p/implicit/properties.hh +++ b/dumux/porousmediumflow/3p/implicit/properties.hh @@ -27,8 +27,23 @@ #ifndef DUMUX_3P_PROPERTIES_HH #define DUMUX_3P_PROPERTIES_HH +#include +#include + +#include +#include +#include + +#include +#include +#include #include +#include "indices.hh" +#include "volumevariables.hh" +#include "properties.hh" +#include "vtkoutputfields.hh" + namespace Dumux { @@ -37,34 +52,117 @@ namespace Properties ////////////////////////////////////////////////////////////////// // Type tags ////////////////////////////////////////////////////////////////// -//! The type tags for the implicit isothermal one-phase two-component problems -NEW_TYPE_TAG(ThreeP); -NEW_TYPE_TAG(BoxThreeP, INHERITS_FROM(BoxModel, ThreeP)); -NEW_TYPE_TAG(CCThreeP, INHERITS_FROM(CCTpfaModel, ThreeP)); - -//! The type tags for the corresponding non-isothermal problems +NEW_TYPE_TAG(ThreeP, INHERITS_FROM(PorousMediumFlow, NumericModel, LinearSolverTypeTag)); NEW_TYPE_TAG(ThreePNI, INHERITS_FROM(ThreeP, NonIsothermal)); -NEW_TYPE_TAG(BoxThreePNI, INHERITS_FROM(BoxModel, ThreePNI)); -NEW_TYPE_TAG(CCThreePNI, INHERITS_FROM(CCTpfaModel, ThreePNI)); +////////////////////////////////////////////////////////////////// +// Properties for the isothermal 3p model +////////////////////////////////////////////////////////////////// +/*! + * \brief Set the property for the number of fluid phases. + * + * We just forward the number from the fluid system and use an static + * assert to make sure it is 3. + */ +SET_PROP(ThreeP, NumPhases) +{ + private: + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + public: + static const int value = FluidSystem::numPhases; + static_assert(value == 3, + "Only fluid systems with 3 phases are supported by the 3p model!"); +}; + +SET_PROP(ThreeP, NumComponents) +{ + private: + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + public: + static const int value = FluidSystem::numComponents; + static_assert(value == 3, + "Only fluid systems with 3 components are supported by the 3p model!"); +}; + +SET_INT_PROP(ThreeP, NumEq, 3); //!< set the number of equations to 3 + +/*! + * \brief Set the property for the material parameters by extracting + * it from the material law. + */ +SET_TYPE_PROP(ThreeP, MaterialLawParams, typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params); + +//! The local residual function of the conservation equations +SET_TYPE_PROP(ThreeP, LocalResidual, ImmiscibleLocalResidual); + +//! Enable advection +SET_BOOL_PROP(ThreeP, EnableAdvection, true); + +//! disable molecular diffusion for the 3p model +SET_BOOL_PROP(ThreeP, EnableMolecularDiffusion, false); + +//! Isothermal model by default +SET_BOOL_PROP(ThreeP, EnableEnergyBalance, false); + +//! the VolumeVariables property +SET_TYPE_PROP(ThreeP, VolumeVariables, ThreePVolumeVariables); + +//! The indices required by the isothermal 3p model +SET_TYPE_PROP(ThreeP, Indices, ThreePIndices); + +//! The spatial parameters to be employed. +//! Use ImplicitSpatialParams by default. +SET_TYPE_PROP(ThreeP, SpatialParams, ImplicitSpatialParams); + +SET_PROP(ThreeP, FluidState) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; +public: + typedef ImmiscibleFluidState type; +}; + +//! Set the vtk output fields specific to the ThreeP model +SET_TYPE_PROP(ThreeP, VtkOutputFields, ThreePVtkOutputFields); + + +///////////////////////////////////////////////// +// Properties for the non-isothermal 3p model +///////////////////////////////////////////////// + +//! Somerton is used as default model to compute the effective thermal heat conductivity +SET_PROP(ThreePNI, ThermalConductivityModel) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; +public: + typedef ThermalConductivitySomerton type; +}; + +//! temperature is already written by the isothermal model +SET_BOOL_PROP(ThreePNI, NiOutputLevel, 0); ////////////////////////////////////////////////////////////////// -// Property tags +// Property values for isothermal model required for the general non-isothermal model ////////////////////////////////////////////////////////////////// -NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system -NEW_PROP_TAG(NumComponents); //!< Number of components in the system -NEW_PROP_TAG(Indices); //!< Enumerations for the model -NEW_PROP_TAG(SpatialParams); //!< The type of the spatial parameters -NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations -NEW_PROP_TAG(FluidState); //!< the phases' state +//set isothermal VolumeVariables +SET_TYPE_PROP(ThreePNI, IsothermalVolumeVariables, ThreePVolumeVariables); + +//set isothermal LocalResidual +SET_TYPE_PROP(ThreePNI, IsothermalLocalResidual, ImmiscibleLocalResidual); + +//set isothermal Indices +SET_TYPE_PROP(ThreePNI, IsothermalIndices, ThreePIndices); -NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (extracted from the spatial parameters) -NEW_PROP_TAG(MaterialLawParams); //!< The parameters of the material law (extracted from the spatial parameters) +//set isothermal NumEq +SET_INT_PROP(ThreePNI, IsothermalNumEq, 3); +} // end namespace Properties -NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem -NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient -} -} +} // end namespace Dumux #endif diff --git a/dumux/porousmediumflow/3p/implicit/propertydefaults.hh b/dumux/porousmediumflow/3p/implicit/propertydefaults.hh deleted file mode 100644 index e524fe4142e307452b0cecfa2e327b1c2c7ae007..0000000000000000000000000000000000000000 --- a/dumux/porousmediumflow/3p/implicit/propertydefaults.hh +++ /dev/null @@ -1,166 +0,0 @@ -// -*- 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 . * - *****************************************************************************/ -/*! - * \ingroup ThreePModel - */ -/*! - * \file - * - * \brief Defines default values for most properties required by the - * three-phase fully implicit model. - */ -#ifndef DUMUX_3P_PROPERTY_DEFAULTS_HH -#define DUMUX_3P_PROPERTY_DEFAULTS_HH - -#include "indices.hh" - -#include "model.hh" -#include "indices.hh" -#include "volumevariables.hh" -#include "properties.hh" - -#include -#include -#include -#include -#include -#include -#include - - -namespace Dumux -{ - -namespace Properties { -////////////////////////////////////////////////////////////////// -// Property values -////////////////////////////////////////////////////////////////// - -/*! - * \brief Set the property for the number of fluid phases. - * - * We just forward the number from the fluid system and use an static - * assert to make sure it is 2. - */ -SET_PROP(ThreeP, NumPhases) -{ - private: - typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - - public: - static const int value = FluidSystem::numPhases; - static_assert(value == 3, - "Only fluid systems with 3 phases are supported by the 3p model!"); -}; - -SET_PROP(ThreeP, NumComponents) -{ - private: - typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - - public: - static const int value = FluidSystem::numComponents; - static_assert(value == 3, - "Only fluid systems with 3 components are supported by the 3p model!"); -}; - -SET_INT_PROP(ThreeP, NumEq, 3); //!< set the number of equations to 2 - -/*! - * \brief Set the property for the material parameters by extracting - * it from the material law. - */ -SET_TYPE_PROP(ThreeP, MaterialLawParams, typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params); - -//! The local residual function of the conservation equations -SET_TYPE_PROP(ThreeP, LocalResidual, ImmiscibleLocalResidual); - -//! Enable advection -SET_BOOL_PROP(ThreeP, EnableAdvection, true); - -//! disable molecular diffusion for the 3p model -SET_BOOL_PROP(ThreeP, EnableMolecularDiffusion, false); - -//! Isothermal model by default -SET_BOOL_PROP(ThreeP, EnableEnergyBalance, false); - -//! the Model property -SET_TYPE_PROP(ThreeP, Model, ThreePModel); - -//! the VolumeVariables property -SET_TYPE_PROP(ThreeP, VolumeVariables, ThreePVolumeVariables); - -//! The indices required by the isothermal 3p model -SET_TYPE_PROP(ThreeP, Indices, ThreePIndices); - -//! The spatial parameters to be employed. -//! Use ImplicitSpatialParams by default. -SET_TYPE_PROP(ThreeP, SpatialParams, ImplicitSpatialParams); - -SET_PROP(ThreeP, FluidState) -{ -private: - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; -public: - typedef ImmiscibleFluidState type; -}; - -// disable velocity output by default - -// enable gravity by default -SET_BOOL_PROP(ThreeP, ProblemEnableGravity, true); - -//! Somerton is used as default model to compute the effective thermal heat conductivity -SET_PROP(ThreePNI, ThermalConductivityModel) -{ -private: - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; -public: - typedef ThermalConductivitySomerton type; -}; - -//! temperature is already written by the isothermal model -SET_BOOL_PROP(ThreePNI, NiOutputLevel, 0); - -////////////////////////////////////////////////////////////////// -// Property values for isothermal model required for the general non-isothermal model -////////////////////////////////////////////////////////////////// - -// set isothermal Model -SET_TYPE_PROP(ThreePNI, IsothermalModel, ThreePModel); - -//set isothermal VolumeVariables -SET_TYPE_PROP(ThreePNI, IsothermalVolumeVariables, ThreePVolumeVariables); - -//set isothermal LocalResidual -SET_TYPE_PROP(ThreePNI, IsothermalLocalResidual, ImmiscibleLocalResidual); - -//set isothermal Indices -SET_TYPE_PROP(ThreePNI, IsothermalIndices, ThreePIndices); - -//set isothermal NumEq -SET_INT_PROP(ThreePNI, IsothermalNumEq, 3); - -} - -} - -#endif diff --git a/dumux/porousmediumflow/3p/implicit/volumevariables.hh b/dumux/porousmediumflow/3p/implicit/volumevariables.hh index 758b74e1ddb9d08d52a9ccd0ad32dc6581731643..0d65c7e7e34213a94901a43e827eddb1cc3daa27 100644 --- a/dumux/porousmediumflow/3p/implicit/volumevariables.hh +++ b/dumux/porousmediumflow/3p/implicit/volumevariables.hh @@ -25,10 +25,10 @@ #ifndef DUMUX_3P_VOLUME_VARIABLES_HH #define DUMUX_3P_VOLUME_VARIABLES_HH -#include #include #include #include +#include #include "properties.hh" namespace Dumux @@ -75,7 +75,7 @@ class ThreePVolumeVariables : public ImplicitVolumeVariables static const Scalar R; // universal gas constant - enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; + enum { isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box }; enum { dofCodim = isBox ? dim : 0 }; public: diff --git a/dumux/porousmediumflow/3p/implicit/vtkoutputfields.hh b/dumux/porousmediumflow/3p/implicit/vtkoutputfields.hh new file mode 100644 index 0000000000000000000000000000000000000000..e6d9077063cce6e373d70d4783e748024c7e32fc --- /dev/null +++ b/dumux/porousmediumflow/3p/implicit/vtkoutputfields.hh @@ -0,0 +1,62 @@ +// -*- 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 . * + *****************************************************************************/ +/*! + * \file + * \brief Adds vtk output fields specific to the twop model + */ +#ifndef DUMUX_THREEP_VTK_OUTPUT_FIELDS_HH +#define DUMUX_THREEP_VTK_OUTPUT_FIELDS_HH + +#include + +namespace Dumux +{ + +/*! + * \ingroup ThreeP, InputOutput + * \brief Adds vtk output fields specific to the twop model + */ +template +class ThreePVtkOutputFields +{ + using Indices = typename GET_PROP_TYPE(TypeTag, Indices); + using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); +public: + template + static void init(VtkOutputModule& vtk) + { + // register standardized vtk output fields + vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.saturation(Indices::wPhaseIdx); }, "sw"); + vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.saturation(Indices::nPhaseIdx); },"sn"); + vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.saturation(Indices::gPhaseIdx); },"sg"); + vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.pressure(Indices::wPhaseIdx); },"pw"); + vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.pressure(Indices::nPhaseIdx); },"pn"); + vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.pressure(Indices::gPhaseIdx); },"pg"); + vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.density(Indices::wPhaseIdx); },"rhow"); + vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.density(Indices::nPhaseIdx); },"rhon"); + vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.density(Indices::gPhaseIdx); },"rhog"); + vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.porosity(); },"porosity"); + vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.permeability(); },"permeability"); + vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.temperature(); },"temperature"); + } +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/porousmediumflow/3p3c/implicit/properties.hh b/dumux/porousmediumflow/3p3c/implicit/properties.hh index 972e69f1b2fb2bdbbed44e62bcc107cbeaf22bdc..3ffcc2e23636f500ad991613ee37cb88b9a3b906 100644 --- a/dumux/porousmediumflow/3p3c/implicit/properties.hh +++ b/dumux/porousmediumflow/3p3c/implicit/properties.hh @@ -30,6 +30,9 @@ #ifndef DUMUX_3P3C_PROPERTIES_HH #define DUMUX_3P3C_PROPERTIES_HH +#include +#include +#include #include namespace Dumux @@ -37,42 +40,13 @@ namespace Dumux namespace Properties { -////////////////////////////////////////////////////////////////// -// Type tags -////////////////////////////////////////////////////////////////// - //! The type tags for the implicit three-phase three-component problems -NEW_TYPE_TAG(ThreePThreeC); -NEW_TYPE_TAG(BoxThreePThreeC, INHERITS_FROM(BoxModel, ThreePThreeC)); -NEW_TYPE_TAG(CCThreePThreeC, INHERITS_FROM(CCModel, ThreePThreeC)); +NEW_TYPE_TAG(ThreePThreeC, INHERITS_FROM(PorousMediumFlow, NumericModel, LinearSolverTypeTag)); //! The type tags for the corresponding non-isothermal problems NEW_TYPE_TAG(ThreePThreeCNI, INHERITS_FROM(ThreePThreeC, NonIsothermal)); -NEW_TYPE_TAG(BoxThreePThreeCNI, INHERITS_FROM(BoxModel, ThreePThreeCNI)); -NEW_TYPE_TAG(CCThreePThreeCNI, INHERITS_FROM(CCModel, ThreePThreeCNI)); - -////////////////////////////////////////////////////////////////// -// Property tags -////////////////////////////////////////////////////////////////// - -NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system -NEW_PROP_TAG(NumComponents); //!< Number of fluid components in the system -NEW_PROP_TAG(Indices); //!< Enumerations for the model -NEW_PROP_TAG(SpatialParams); //!< The type of the spatial parameters -NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations -NEW_PROP_TAG(FluidState); //!< Type of the fluid state to be used - -NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (extracted from the spatial parameters) -NEW_PROP_TAG(MaterialLawParams); //!< The parameters of the material law (extracted from the spatial parameters) -NEW_PROP_TAG(EffectiveDiffusivityModel); //!< The employed model for the computation of the effective diffusivity -NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem -NEW_PROP_TAG(UseConstraintSolver); //!< Determines whether a constraint solver should be used explicitly -NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient -NEW_PROP_TAG(TauTortuosity); //!< Tortuosity value (tau) used in macroscopic diffusion -NEW_PROP_TAG(UseMoles);//!< Defines whether mole (true) or mass (false) fractions are used -NEW_PROP_TAG(ReplaceCompEqIdx); //!< The index of the total mass balance equation, if one component balance is replaced (ReplaceCompEqIdx < NumComponents) -} -} +} // end namespace Properties +} // end namespace Dumux #endif diff --git a/dumux/porousmediumflow/3p3c/implicit/propertydefaults.hh b/dumux/porousmediumflow/3p3c/implicit/propertydefaults.hh index a1ff05ff0ea6efff7926c68a741c350253748bb4..0c8137cb8622de630b7355196703bdda459d3a86 100644 --- a/dumux/porousmediumflow/3p3c/implicit/propertydefaults.hh +++ b/dumux/porousmediumflow/3p3c/implicit/propertydefaults.hh @@ -158,9 +158,6 @@ SET_PROP(ThreePThreeC, EffectiveDiffusivityModel) // disable velocity output by default -// enable gravity by default -SET_BOOL_PROP(ThreePThreeC, ProblemEnableGravity, true); - //! default value for the forchheimer coefficient // Source: Ward, J.C. 1964 Turbulent flow in porous media. ASCE J. Hydraul. Div 90. // Actually the Forchheimer coefficient is also a function of the dimensions of the diff --git a/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh b/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh index 2ebb373e23d5834c65bb62e03375e6463a54fa44..c9c3fbf61883782256b8e805bc1e2a60d2fd75fb 100644 --- a/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh +++ b/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh @@ -26,14 +26,15 @@ #include +#include +#include +#include +#include #include -#include -#include -#include - #include #include #include + #include "3pnispatialparams.hh" @@ -45,10 +46,13 @@ class ThreePNIConductionProblem; namespace Properties { -NEW_TYPE_TAG(ThreePNIConductionProblem, INHERITS_FROM(ThreePNI, ThreePNISpatialParams)); -NEW_TYPE_TAG(ThreePNIConductionBoxProblem, INHERITS_FROM(BoxModel, ThreePNIConductionProblem)); -NEW_TYPE_TAG(ThreePNIConductionCCProblem, INHERITS_FROM(CCTpfaModel, ThreePNIConductionProblem)); -NEW_TYPE_TAG(ThreePNIConductionCCMpfaProblem, INHERITS_FROM(CCMpfaModel, ThreePNIConductionProblem)); + +NEW_PROP_TAG(FVGridGeometry); + +NEW_TYPE_TAG(ThreePNIConductionProblem, INHERITS_FROM(ThreePNI)); +NEW_TYPE_TAG(ThreePNIConductionBoxProblem, INHERITS_FROM(BoxModel, ThreePNIConductionProblem, ThreePNISpatialParams)); +NEW_TYPE_TAG(ThreePNIConductionCCProblem, INHERITS_FROM(CCTpfaModel, ThreePNIConductionProblem, ThreePNISpatialParams)); +NEW_TYPE_TAG(ThreePNIConductionCCMpfaProblem, INHERITS_FROM(CCMpfaModel, ThreePNIConductionProblem, ThreePNISpatialParams)); // Set the grid type SET_TYPE_PROP(ThreePNIConductionProblem, Grid, Dune::YaspGrid<2>); @@ -67,7 +71,7 @@ SET_TYPE_PROP(ThreePNIConductionProblem, SpatialParams, ThreePNISpatialParams); -} +}// end namespace Properties /*! @@ -94,34 +98,34 @@ SET_TYPE_PROP(ThreePNIConductionProblem, * ./test_cc3pniconduction -ParameterFile ./test_cc3pniconduction.input */ template -class ThreePNIConductionProblem : public ImplicitPorousMediaProblem +class ThreePNIConductionProblem : public PorousMediumFlowProblem { - using ParentType = ImplicitPorousMediaProblem; + using ParentType = PorousMediumFlowProblem; using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes); - using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager); using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel); using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); + using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector); + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using IapwsH2O = H2O; - using VtkOutputModule = typename GET_PROP_TYPE(TypeTag, VtkOutputModule); + using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector); // copy some indices for convenience using Indices = typename GET_PROP_TYPE(TypeTag, Indices); enum { // world dimension - dimWorld = GridView::dimensionworld + dimWorld = GridView::dimensionworld, + dim=GridView::dimension }; - enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; - enum { dofCodim = isBox ? dimWorld : 0 }; - enum { // index of the primary variables pressureIdx = Indices::pressureIdx, @@ -136,48 +140,34 @@ class ThreePNIConductionProblem : public ImplicitPorousMediaProblem using GlobalPosition = Dune::FieldVector; public: - ThreePNIConductionProblem(TimeManager &timeManager, const GridView &gridView) - : ParentType(timeManager, gridView) + ThreePNIConductionProblem(std::shared_ptr fvGridGeometry) + : ParentType(fvGridGeometry) { //initialize fluid system FluidSystem::init(); - name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, - std::string, - Problem, - Name); - outputInterval_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, - int, Problem, OutputInterval); - + name_ = getParam("Problem.Name"); temperatureHigh_ = 300.0; + temperatureExact_.resize(fvGridGeometry->numDofs()); + } - } - - - bool shouldWriteOutput() const + //! get the analytical temperature + const std::vector& getExactTemperature() { - return - this->timeManager().timeStepIndex() == 0 || - this->timeManager().timeStepIndex() % outputInterval_ == 0 || - this->timeManager().episodeWillBeFinished() || - this->timeManager().willBeFinished(); + return temperatureExact_; } - - /*! - * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write. - */ - void addVtkOutputFields(VtkOutputModule& outputModule) const + //! udpate the analytical temperature + void updateExactTemperature(const SolutionVector& curSol, Scalar time) { - auto& temperatureExact = outputModule.createScalarField("temperatureExact", dofCodim); + const auto someElement = *(elements(this->fvGridGeometry().gridView()).begin()); - const auto someElement = *(elements(this->gridView()).begin()); - const auto someElemSol = this->model().elementSolution(someElement, this->model().curSol()); + ElementSolutionVector someElemSol(someElement, curSol, this->fvGridGeometry()); const auto someInitSol = initialAtPos(someElement.geometry().center()); - auto someFvGeometry = localView(this->model().fvGridGeometry()); - someFvGeometry.bindElement(someElement); - const auto someScv = *(scvs(someFvGeometry).begin()); + auto fvGeometry = localView(this->fvGridGeometry()); + fvGeometry.bindElement(someElement); + const auto someScv = *(scvs(fvGeometry).begin()); VolumeVariables volVars; volVars.update(someElemSol, *this, someElement, someScv); @@ -189,26 +179,27 @@ public: const auto heatCapacityS = this->spatialParams().solidHeatCapacity(someElement, someScv, someElemSol); const auto storage = densityW*heatCapacityW*porosity + densityS*heatCapacityS*(1 - porosity); const auto effectiveThermalConductivity = ThermalConductivityModel::effectiveThermalConductivity(volVars, this->spatialParams(), - someElement, someFvGeometry, someScv); + someElement, fvGeometry, someScv); using std::max; - Scalar time = max(this->timeManager().time() + this->timeManager().timeStepSize(), 1e-10); - - for (const auto& element : elements(this->gridView())) + time = max(time, 1e-10); + for (const auto& element : elements(this->fvGridGeometry().gridView())) { - auto fvGeometry = localView(this->model().fvGridGeometry()); + auto fvGeometry = localView(this->fvGridGeometry()); fvGeometry.bindElement(element); for (auto&& scv : scvs(fvGeometry)) { - auto globalIdx = scv.dofIndex(); - const auto& globalPos = scv.dofPosition(); - using std::erf; - using std::sqrt; - temperatureExact[globalIdx] = temperatureHigh_ + (someInitSol[temperatureIdx] - temperatureHigh_) + auto globalIdx = scv.dofIndex(); + const auto& globalPos = scv.dofPosition(); + using std::erf; + using std::sqrt; + temperatureExact_[globalIdx] = temperatureHigh_ + (someInitSol[temperatureIdx] - temperatureHigh_) *erf(0.5*sqrt(globalPos[0]*globalPos[0]*storage/time/effectiveThermalConductivity)); + } } } + /*! * \name Problem parameters */ @@ -239,7 +230,7 @@ public: BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const { BoundaryTypes values; - if(globalPos[0] < eps_ || globalPos[0] > this->bBoxMax()[0] - eps_) + if(globalPos[0] < eps_ || globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_) { values.setAllDirichlet(); } @@ -265,22 +256,20 @@ public: return values; } - /*! + /*! * \brief Evaluate the boundary conditions for a neumann * boundary segment. * - * \param element The finite element - * \param fvGeometry The finite-volume geometry in the box scheme - * \param elemVolVars The element volume variables - * \param scvf The subcontrolvolume face - * Negative values mean influx. + * \param values Stores the Neumann values for the conservation equations in + * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$ + * \param globalPos The position of the integration point of the boundary segment. + * + * For this method, the \a values parameter stores the mass flux + * in normal direction of each phase. Negative values mean influx. */ - PrimaryVariables neumann(const Element &element, - const FVElementGeometry& fvGeometry, - const ElementVolumeVariables& elemVolVars, - const SubControlVolumeFace& scvf) const + NeumannFluxes neumannAtPos(const GlobalPosition &globalPos) const { - return PrimaryVariables(0.0); + return NeumannFluxes(0.0); } // \} @@ -317,20 +306,19 @@ public: { PrimaryVariables values; values[pressureIdx] = 1e5; // initial condition for the pressure - values[swIdx] = 1.; // initial condition for the wetting phase saturation + values[swIdx] = 1.0; // initial condition for the wetting phase saturation values[snIdx] = 1e-5; // initial condition for the non-wetting phase saturation - values[temperatureIdx] = 290.; + values[temperatureIdx] = 290; return values; } // \} private: - Scalar temperatureHigh_; static constexpr Scalar eps_ = 1e-6; std::string name_; - int outputInterval_; + std::vector temperatureExact_; }; } //end namespace diff --git a/test/porousmediumflow/3p/implicit/3pniconvectionproblem.hh b/test/porousmediumflow/3p/implicit/3pniconvectionproblem.hh index a88a062336ec26032fd1099d331964d341191be9..1178f4d46e104db1257cbad3c9e36f0a8789a46a 100644 --- a/test/porousmediumflow/3p/implicit/3pniconvectionproblem.hh +++ b/test/porousmediumflow/3p/implicit/3pniconvectionproblem.hh @@ -26,16 +26,16 @@ #include +#include +#include +#include +#include #include -#include -#include -#include - #include #include #include -#include "3pnispatialparams.hh" +#include "3pnispatialparams.hh" namespace Dumux { @@ -94,34 +94,32 @@ SET_TYPE_PROP(ThreePNIConvectionProblem, * ./test_cc3pcniconvection -ParameterFile ./test_cc3pniconvection.input */ template -class ThreePNIConvectionProblem : public ImplicitPorousMediaProblem +class ThreePNIConvectionProblem : public PorousMediumFlowProblem { - using ParentType = ImplicitPorousMediaProblem; + using ParentType = PorousMediumFlowProblem; using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes); - using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager); using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel); using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); + using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector); + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using IapwsH2O = H2O; - using VtkOutputModule = typename GET_PROP_TYPE(TypeTag, VtkOutputModule); // copy some indices for convenience - typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + using Indices = typename GET_PROP_TYPE(TypeTag, Indices); enum { // world dimension dimWorld = GridView::dimensionworld }; - enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; - enum { dofCodim = isBox ? dimWorld : 0 }; - enum { // index of the primary variables pressureIdx = Indices::pressureIdx, @@ -133,61 +131,46 @@ class ThreePNIConvectionProblem : public ImplicitPorousMediaProblem energyEqIdx = Indices::energyEqIdx }; - - typedef typename GridView::template Codim<0>::Entity Element; - typedef typename GridView::Intersection Intersection; - - typedef Dune::FieldVector GlobalPosition; - + using Element = typename GridView::template Codim<0>::Entity; + using Intersection = typename GridView::Intersection; + using GlobalPosition = Dune::FieldVector; public: - ThreePNIConvectionProblem(TimeManager &timeManager, const GridView &gridView) - : ParentType(timeManager, gridView) + ThreePNIConvectionProblem(std::shared_ptr fvGridGeometry) + : ParentType(fvGridGeometry) { //initialize fluid system FluidSystem::init(); - name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, - std::string, - Problem, Name); - outputInterval_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, - int, Problem, OutputInterval); - darcyVelocity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, - Scalar, Problem, DarcyVelocity); + name_ = getParam("Problem.Name"); + outputInterval_ = getParam("Problem.OutputInterval"); + darcyVelocity_ = getParam("Problem.DarcyVelocity"); temperatureHigh_ = 291.; temperatureLow_ = 290.; pressureHigh_ = 2e5; pressureLow_ = 1e5; - + temperatureExact_.resize(this->fvGridGeometry().numDofs()); } - - bool shouldWriteOutput() const + //! Get exact temperature vector for output + const std::vector& getExactTemperature() { - return - this->timeManager().timeStepIndex() == 0 || - this->timeManager().timeStepIndex() % outputInterval_ == 0 || - this->timeManager().episodeWillBeFinished() || - this->timeManager().willBeFinished(); + return temperatureExact_; } - - /*! - * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write. - */ - void addVtkOutputFields(VtkOutputModule& outputModule) const + //! udpate the analytical temperature + void updateExactTemperature(const SolutionVector& curSol, Scalar time) { - auto& temperatureExact = outputModule.createScalarField("temperatureExact", dofCodim); + const auto someElement = *(elements(this->fvGridGeometry().gridView()).begin()); - const auto someElement = *(elements(this->gridView()).begin()); - const auto someElemSol = this->model().elementSolution(someElement, this->model().curSol()); + ElementSolutionVector someElemSol(someElement, curSol, this->fvGridGeometry()); const auto someInitSol = initialAtPos(someElement.geometry().center()); - auto someFvGeometry = localView(this->model().fvGridGeometry()); - someFvGeometry.bindElement(someElement); - const auto someScv = *(scvs(someFvGeometry).begin()); + auto fvGeometry = localView(this->fvGridGeometry()); + fvGeometry.bindElement(someElement); + const auto someScv = *(scvs(fvGeometry).begin()); VolumeVariables volVars; volVars.update(someElemSol, *this, someElement, someScv); @@ -202,24 +185,23 @@ public: std::cout << "storage: " << storageTotal << '\n'; using std::max; - const Scalar time = max(this->timeManager().time() + this->timeManager().timeStepSize(), 1e-10); + time = max(time, 1e-10); const Scalar retardedFrontVelocity = darcyVelocity_*storageW/storageTotal/porosity; std::cout << "retarded velocity: " << retardedFrontVelocity << '\n'; - for (const auto& element : elements(this->gridView())) + for (const auto& element : elements(this->fvGridGeometry().gridView())) { - auto fvGeometry = localView(this->model().fvGridGeometry()); + auto fvGeometry = localView(this->fvGridGeometry()); fvGeometry.bindElement(element); for (auto&& scv : scvs(fvGeometry)) { auto dofIdxGlobal = scv.dofIndex(); auto dofPosition = scv.dofPosition(); - temperatureExact[dofIdxGlobal] = (dofPosition[0] < retardedFrontVelocity*time) ? temperatureHigh_ : temperatureLow_; + temperatureExact_[dofIdxGlobal] = (dofPosition[0] < retardedFrontVelocity*time) ? temperatureHigh_ : temperatureLow_; } } } - /*! * \name Problem parameters */ @@ -251,7 +233,7 @@ public: BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const { BoundaryTypes values; - if(globalPos[0] > this->bBoxMax()[0] - eps_) + if(globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_) { values.setAllDirichlet(); } @@ -309,23 +291,6 @@ public: */ // \{ - /*! - * \brief Evaluate the source term for all phases within a given - * sub-control-volume. - * - * \param globalPos The position for which the source should be evaluated - * - * Returns the rate mass of a component is generated or annihilate - * per volume unit. Positive values mean that mass is created, - * negative ones mean that it vanishes. - * - * The units must be according to either using mole or mass fractions. (mole/(m^3*s) or kg/(m^3*s)) - */ - PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const - { - return PrimaryVariables(0.0); - } - /*! * \brief Evaluate the initial value for a control volume. * @@ -353,7 +318,9 @@ private: static constexpr Scalar eps_ = 1e-6; std::string name_; int outputInterval_; + std::vector temperatureExact_; }; -} //end namespace +} //end namespace Dumux + #endif diff --git a/test/porousmediumflow/3p/implicit/3pnispatialparams.hh b/test/porousmediumflow/3p/implicit/3pnispatialparams.hh index 6fdd140d43ea3268e4537f0db0943b48dc9e8486..9a75e54ef610368360045b057d7b9337b24be64d 100644 --- a/test/porousmediumflow/3p/implicit/3pnispatialparams.hh +++ b/test/porousmediumflow/3p/implicit/3pnispatialparams.hh @@ -58,11 +58,11 @@ SET_PROP(ThreePNISpatialParams, MaterialLaw) private: // define the material law which is parameterized by effective // saturations - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef RegularizedParkerVanGen3P EffectiveLaw; + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using EffectiveLaw = RegularizedParkerVanGen3P; public: // define the material law parameterized by absolute saturations - typedef EffToAbsLaw type; + using type = EffToAbsLaw; }; } @@ -94,8 +94,8 @@ public: using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw); using MaterialLawParams = typename MaterialLaw::Params; - ThreePNISpatialParams(const Problem& problem, const GridView &gridView) - : ParentType(problem, gridView) + ThreePNISpatialParams(const Problem& problem) + : ParentType(problem) { permeability_ = 1e-10; porosity_ = 0.4; diff --git a/test/porousmediumflow/3p/implicit/CMakeLists.txt b/test/porousmediumflow/3p/implicit/CMakeLists.txt index 4cf0fd0b92e999191d0bc269090ca5c400be1593..fd9c80347c044963ad9c2f83f9396e9c42bb8b53 100644 --- a/test/porousmediumflow/3p/implicit/CMakeLists.txt +++ b/test/porousmediumflow/3p/implicit/CMakeLists.txt @@ -1,52 +1,58 @@ add_input_file_links() # isothermal tests -add_dumux_test(test_box3p test_box3p test_box3p.cc - python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/infiltration3pbox-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/infiltration3pbox-00008.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_box3p") +dune_add_test(NAME test_box3p + SOURCES test_box3p.cc + COMMAND python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/infiltration3pbox-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/infiltration3pbox-00008.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_box3p") -add_dumux_test(test_cc3p test_cc3p test_cc3p.cc - python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/infiltration3pcc-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/infiltration3pcc-00008.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_cc3p") +dune_add_test(NAME test_cc3p + SOURCES test_cc3p.cc + COMMAND python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/infiltration3pcc-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/infiltration3pcc-00008.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_cc3p") # non-isothermal tests -add_dumux_test(test_box3pniconvection test_box3pniconvection test_box3pniconvection.cc - python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/3pniconvectionbox-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_box3pniconvection-00010.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_box3pniconvection" - --zeroThreshold {"velocity_w \(m/s\)_1":1e-8}) +dune_add_test(NAME test_box3pniconvection + SOURCES test_box3pniconvection.cc + COMMAND python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/3pniconvectionbox-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_box3pniconvection-00010.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_box3pniconvection" + --zeroThreshold {"velocity_w \(m/s\)_1":1e-8}) -add_dumux_test(test_cc3pniconvection test_cc3pniconvection test_cc3pniconvection.cc - python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/3pniconvectioncc-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_cc3pniconvection-00010.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_cc3pniconvection" - --zeroThreshold {"velocity_w \(m/s\)_1":1e-8}) +dune_add_test(NAME test_cc3pniconvection + SOURCES test_cc3pniconvection.cc + COMMAND python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/3pniconvectioncc-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_cc3pniconvection-00010.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_cc3pniconvection" + --zeroThreshold {"velocity_w \(m/s\)_1":1e-8}) -add_dumux_test(test_box3pniconduction test_box3pniconduction test_box3pniconduction.cc - python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/3pniconductionbox-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_box3pniconduction-00006.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_box3pniconduction" - --zeroThreshold {"velocity_w \(m/s\)_1":1e-8}) +dune_add_test(NAME test_box3pniconduction + SOURCES test_box3pniconduction.cc + COMMAND python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/3pniconductionbox-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_box3pniconduction-00006.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_box3pniconduction" + --zeroThreshold {"velocity_w \(m/s\)_1":1e-8}) -add_dumux_test(test_cc3pniconduction test_cc3pniconduction test_cc3pniconduction.cc - python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/3pniconductioncc-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_cc3pniconduction-00006.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_cc3pniconduction" - --zeroThreshold {"velocity_w \(m/s\)_1":1e-8}) +dune_add_test(NAME test_cc3pniconduction + SOURCES test_cc3pniconduction.cc + COMMAND python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/3pniconductioncc-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_cc3pniconduction-00006.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_cc3pniconduction" + --zeroThreshold {"velocity_w \(m/s\)_1":1e-8}) #install sources install(FILES diff --git a/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh b/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh index 89a8f9c264f76c8907931a760f52080b2eaea98f..5948e370ba2fb5d5479c54ad9aa5e8361af48dd5 100644 --- a/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh +++ b/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh @@ -25,11 +25,13 @@ #ifndef DUMUX_INFILTRATION_THREEP_PROBLEM_HH #define DUMUX_INFILTRATION_THREEP_PROBLEM_HH +#include +#include +#include +#include #include -#include -#include - #include +#include #include "infiltration3pspatialparams.hh" @@ -55,10 +57,7 @@ SET_TYPE_PROP(InfiltrationThreePProblem, FluidSystem, FluidSystems::H2OAirMesitylene); - -// Maximum tolerated relative error in the Newton method -SET_SCALAR_PROP(InfiltrationThreePProblem, NewtonMaxRelativeShift, 1e-4); -} +}// end namespace Properties /*! * \ingroup ThreePModel @@ -91,18 +90,14 @@ SET_SCALAR_PROP(InfiltrationThreePProblem, NewtonMaxRelativeShift, 1e-4); * ./naplinfiltrationexercise -parameterFile naplinfiltrationexercise.input * */ template -class InfiltrationThreePProblem : public ImplicitPorousMediaProblem +class InfiltrationThreePProblem : public PorousMediumFlowProblem { - typedef ImplicitPorousMediaProblem ParentType; - - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + using ParentType = PorousMediumFlowProblem; - typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; - typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Indices = typename GET_PROP_TYPE(TypeTag, Indices); - // copy some indices for convenience - typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; enum { pressureIdx = Indices::pressureIdx, swIdx = Indices::swIdx, @@ -113,24 +108,16 @@ class InfiltrationThreePProblem : public ImplicitPorousMediaProblem dimWorld = GridView::dimensionworld }; + using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); + using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector); + using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes); + using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); - typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; - typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes; - typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager; - - typedef typename GridView::template Codim<0>::Entity Element; - typedef typename GridView::template Codim::Entity Vertex; - typedef typename GridView::Intersection Intersection; - - typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; - typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables; - typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; - typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace; - typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - - typedef Dune::FieldVector GlobalPosition; + using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw); + using MaterialLawParams = typename MaterialLaw::Params; - enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; + using GlobalPosition = Dune::FieldVector; public: /*! @@ -139,28 +126,31 @@ public: * \param timeManager The time manager * \param gridView The grid view */ - InfiltrationThreePProblem(TimeManager &timeManager, const GridView &gridView) - : ParentType(timeManager, gridView) + InfiltrationThreePProblem(std::shared_ptr fvGridGeometry) + : ParentType(fvGridGeometry) { temperature_ = 273.15 + 10.0; // -> 10 degrees Celsius FluidSystem::init(282.15, 284.15, 3, 8e4, 3e5, 200); - name_ = GET_RUNTIME_PARAM(TypeTag, std::string, Problem.Name); + name_ = getParam("Problem.Name"); this->spatialParams().plotMaterialLaw(); + time_ = 0.0; } /*! * \name Problem parameters */ // \{ + void setTime(Scalar time) + { time_ = time; } /*! * \brief The problem name. * * This is used as a prefix for files generated by the simulation. */ - const std::string name() const + const std::string& name() const { return name_; } /*! @@ -175,17 +165,6 @@ public: return temperature_; } - /*! - * \brief Returns the source term at specific position in the domain. - * - * \param values The source values for the primary variables - * \param globalPos The position - */ - PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const - { - return PrimaryVariables(0.0); - } - // \} /*! @@ -204,9 +183,8 @@ public: { BoundaryTypes values; - if(globalPos[0] > 500. - eps_) - values.setAllDirichlet(); - else if(globalPos[0] < eps_) + if(globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_ + || globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_) values.setAllDirichlet(); else values.setAllNeumann(); @@ -226,65 +204,32 @@ public: PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const { PrimaryVariables values(0.0); - - Scalar y = globalPos[1]; - Scalar x = globalPos[0]; - Scalar sw, swr=0.12, sgr=0.03; - - if(y > (-1.E-3*x+5) - eps_) - { - Scalar pc = 9.81 * 1000.0 * (y - (-5E-4*x+5)); - if (pc < 0.0) pc = 0.0; - - sw = invertPcgw_(pc, - this->spatialParams().materialLawParamsAtPos(globalPos)); - if (sw < swr) sw = swr; - if (sw > 1.-sgr) sw = 1.-sgr; - - values[pressureIdx] = 1e5 ; - values[swIdx] = sw; - values[snIdx] = 0.; - }else { - values[pressureIdx] = 1e5 + 9.81 * 1000.0 * ((-5E-4*x+5) - y); - values[swIdx] = 1.-sgr; - values[snIdx] = 0.; - } - + initial_(values, globalPos); return values; } - /*! * \brief Evaluate the boundary conditions for a neumann * boundary segment. * - * \param values The neumann values for the conservation equations - * \param element The finite element - * \param fvGeometry The finite-volume geometry in the box scheme - * \param intersection The intersection between element and boundary - * \param scvIdx The local vertex index - * \param boundaryFaceIdx The index of the boundary face + * \param values Stores the Neumann values for the conservation equations in + * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$ + * \param globalPos The position of the integration point of the boundary segment. * * For this method, the \a values parameter stores the mass flux * in normal direction of each phase. Negative values mean influx. */ - PrimaryVariables neumann(const Element &element, - const FVElementGeometry& fvGeometry, - const ElementVolumeVariables& elemVolVars, - const SubControlVolumeFace& scvf) const + NeumannFluxes neumannAtPos(const GlobalPosition &globalPos) const { - PrimaryVariables values(0.0); - - GlobalPosition globalPos = scvf.center(); - // TODO: BOX??? FACES SHOULD HAVE INTEGRATION POINT + NeumannFluxes values(0.0); // negative values for injection - if (this->timeManager().time()<2592000.) + if (time_ < 2592000.0 - eps_) { - if ((globalPos[0] <= 175.+eps_) && (globalPos[0] >= 155.-eps_) && (globalPos[1] >= 10.-eps_)) + if ((globalPos[0] < 175.0 + eps_) && (globalPos[0] > 155.0 - eps_) + && (globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_)) { - values[Indices::contiWEqIdx] = -0.0; - values[Indices::contiNEqIdx] = -0.001, // /*Molfluss, umr. über M(Mesit.)=0,120 kg/mol --> 1.2e-4 kg/(sm) - values[Indices::contiGEqIdx] = -0.0; + // mol fluxes, convert with M(Mesit.)=0,120 kg/mol --> 1.2e-4 kg/(sm) + values[Indices::contiNEqIdx] = -0.001; } } @@ -311,7 +256,6 @@ public: { PrimaryVariables values(0.0); initial_(values, globalPos); - return values; } @@ -323,6 +267,8 @@ public: Scalar temperature() const { return temperature_; } + + private: // internal method for the initial condition (reused for the // dirichlet conditions!) @@ -353,6 +299,7 @@ private: } } + // small solver inverting the pc curve static Scalar invertPcgw_(Scalar pcIn, const MaterialLawParams &pcParams) { Scalar lower,upper; @@ -383,7 +330,9 @@ private: Scalar temperature_; static constexpr Scalar eps_ = 1e-6; std::string name_; + Scalar time_; }; -} //end namespace + +} // end namespace Dumux #endif diff --git a/test/porousmediumflow/3p/implicit/infiltration3pspatialparams.hh b/test/porousmediumflow/3p/implicit/infiltration3pspatialparams.hh index 0e06b9ef59b69737e6409aea94255a0eb0aba4e9..b2f00a114c93523cdc57d1a0038fccddfeb7dcc3 100644 --- a/test/porousmediumflow/3p/implicit/infiltration3pspatialparams.hh +++ b/test/porousmediumflow/3p/implicit/infiltration3pspatialparams.hh @@ -52,11 +52,11 @@ SET_PROP(InfiltrationThreePSpatialParams, MaterialLaw) private: // define the material law which is parameterized by effective // saturations - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef RegularizedParkerVanGen3P EffectiveLaw; + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using EffectiveLaw = RegularizedParkerVanGen3P; public: // define the material law parameterized by absolute saturations - typedef EffToAbsLaw type; + using type = EffToAbsLaw; }; } @@ -69,64 +69,59 @@ SET_PROP(InfiltrationThreePSpatialParams, MaterialLaw) template class InfiltrationThreePSpatialParams : public ImplicitSpatialParams { - typedef ImplicitSpatialParams ParentType; - - typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; - typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename Grid::ctype CoordScalar; + using ParentType = ImplicitSpatialParams; + + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using Element = typename GridView::template Codim<0>::Entity; + using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); + using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag,ElementSolutionVector); enum { dimWorld=GridView::dimensionworld }; - typedef Dune::FieldVector GlobalPosition; - - typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables; - typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; - typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; - typedef typename GridView::template Codim<0>::Entity Element; - + using GlobalPosition = Dune::FieldVector; public: // export permeability type using PermeabilityType = Scalar; //get the material law from the property system - typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; - typedef typename MaterialLaw::Params MaterialLawParams; + using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw); + using MaterialLawParams = typename MaterialLaw::Params; /*! * \brief The constructor * * \param gridView The grid view */ - InfiltrationThreePSpatialParams(const Problem& problem, const GridView &gridView) - : ParentType(problem, gridView) + InfiltrationThreePSpatialParams(const Problem& problem) + : ParentType(problem) { // intrinsic permeabilities - fineK_ = GET_RUNTIME_PARAM(TypeTag, Scalar, permeability); - coarseK_ = GET_RUNTIME_PARAM(TypeTag, Scalar, permeability); + fineK_ = getParam("SpatialParams.permeability"); + coarseK_ = getParam("SpatialParams.permeability"); // porosities - porosity_ = GET_RUNTIME_PARAM(TypeTag, Scalar, porosity); - + porosity_ = getParam("SpatialParams.porosity"); + vGAlpha_ = getParam("SpatialParams.vanGenuchtenAlpha"); + vGN_ = getParam("SpatialParams.vanGenuchtenN"); // residual saturations materialParams_.setSwr(0.12); materialParams_.setSnr(0.07); materialParams_.setSgr(0.03); // parameters for the 3phase van Genuchten law - materialParams_.setVgAlpha(GET_RUNTIME_PARAM(TypeTag, Scalar, vanGenuchtenAlpha)); - materialParams_.setVgn(GET_RUNTIME_PARAM(TypeTag, Scalar, vanGenuchtenN)); + materialParams_.setVgAlpha(vGAlpha_); + materialParams_.setVgn(vGN_); materialParams_.setKrRegardsSnr(false); // parameters for adsorption materialParams_.setKdNAPL(0.); materialParams_.setRhoBulk(1500.); - plotFluidMatrixInteractions_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, bool, Output, - PlotFluidMatrixInteractions); + plotFluidMatrixInteractions_ = getParam("Output.PlotFluidMatrixInteractions"); } ~InfiltrationThreePSpatialParams() @@ -144,14 +139,19 @@ public: plotMaterialLaw.plotkr(materialParams_); } - /*! - * \brief Returns the scalar intrinsic permeability \f$[m^2]\f$ + /*! + * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$. * - * \param globalPos The global position + * \param element The element + * \param scv The sub control volume + * \param elemSol The element solution vector + * \return the intrinsic permeability */ - Scalar permeabilityAtPos(const GlobalPosition& globalPos) const + PermeabilityType permeability(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const { - if (isFineMaterial_(globalPos)) + if (isFineMaterial_(scv.dofPosition())) return fineK_; return coarseK_; } @@ -186,6 +186,8 @@ private: Scalar coarseK_; Scalar porosity_; + Scalar vGN_; + Scalar vGAlpha_; MaterialLawParams materialParams_; diff --git a/test/porousmediumflow/3p/implicit/test_box3p.cc b/test/porousmediumflow/3p/implicit/test_box3p.cc index 73718b474e2a89983275ddab5e3ddc05194ff41c..0af33e750d4ac7d5aa18a4e4c11bda0d7741c551 100644 --- a/test/porousmediumflow/3p/implicit/test_box3p.cc +++ b/test/porousmediumflow/3p/implicit/test_box3p.cc @@ -23,7 +23,32 @@ */ #include #include "infiltration3pproblem.hh" -#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include /*! * \brief Provides an interface for customizing error messages associated with @@ -51,9 +76,177 @@ void usage(const char *progName, const std::string &errorMsg) } } -int main(int argc, char** argv) +int main(int argc, char** argv) try { - typedef TTAG(InfiltrationThreePBoxProblem) ProblemTypeTag; - return Dumux::start(argc, argv, usage); -} + using namespace Dumux; + + // define the type tag for this problem + using TypeTag = TTAG(InfiltrationThreePBoxProblem); + + //////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////// + + // 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) + using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator); + GridCreator::makeGrid(); + GridCreator::loadBalance(); + + //////////////////////////////////////////////////////////// + // run instationary non-linear problem on this grid + //////////////////////////////////////////////////////////// + + // we compute on the leaf grid view + const auto& leafGridView = GridCreator::grid().leafGridView(); + + // create the finite volume grid geometry + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + auto fvGridGeometry = std::make_shared(leafGridView); + fvGridGeometry->update(); + + // the problem (initial and boundary conditions) + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + auto problem = std::make_shared(fvGridGeometry); + + // the solution vector + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + SolutionVector x(leafGridView.size(GridView::dimension)); + problem->applyInitialSolution(x); + auto xOld = x; + + // the grid variables + using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables); + auto gridVariables = std::make_shared(problem, fvGridGeometry); + gridVariables->init(x, xOld); + + // get some time loop parameters + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + const auto tEnd = getParam("TimeLoop.TEnd"); + const auto maxDivisions = getParam("TimeLoop.MaxTimeStepDivisions"); + const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); + auto dt = getParam("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("TimeLoop.Restart"); + + // intialize the vtk output module + using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields); + VtkOutputModule vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name()); + VtkOutputFields::init(vtkWriter); //! Add model specific output fields + vtkWriter.write(0.0); + + // instantiate time loop + auto timeLoop = std::make_shared>(restartTime, dt, tEnd); + timeLoop->setMaxTimeStepSize(maxDt); + + // the assembler with time loop for instationary problem + using Assembler = FVAssembler; + auto assembler = std::make_shared(problem, fvGridGeometry, gridVariables, timeLoop); + + // the linear solver + using LinearSolver = Dumux::AMGBackend; + auto linearSolver = std::make_shared(leafGridView, fvGridGeometry->elementMapper()); + + // the non-linear solver + using NewtonController = Dumux::NewtonController; + using NewtonMethod = Dumux::NewtonMethod; + auto newtonController = std::make_shared(leafGridView.comm(), timeLoop); + NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver); + + // time loop + timeLoop->start(); do + { + // set previous solution for storage evaluations + assembler->setPreviousSolution(xOld); + + // try solving the non-linear system + for (int i = 0; i < maxDivisions; ++i) + { + // linearize & solve + auto converged = nonLinearSolver.solve(x); + + if (converged) + break; + + if (!converged && i == maxDivisions-1) + DUNE_THROW(Dune::MathError, + "Newton solver didn't converge after " + << maxDivisions + << " time-step divisions. dt=" + << timeLoop->timeStepSize() + << ".\nThe solutions of the current and the previous time steps " + << "have been saved to restart files."); + } + + // make the new solution the old solution + xOld = x; + gridVariables->advanceTimeStep(); + + // 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 newton controller + timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize())); + problem->setTime(timeLoop->time()+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; + +} +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; +} diff --git a/test/porousmediumflow/3p/implicit/test_box3p.input b/test/porousmediumflow/3p/implicit/test_box3p.input index 6445dff4ae034abcd0ce6f15d05dc78e0c9e4349..b41123ccc30639f2208a2dbab7db7c3dc771c2b2 100644 --- a/test/porousmediumflow/3p/implicit/test_box3p.input +++ b/test/porousmediumflow/3p/implicit/test_box3p.input @@ -1,4 +1,4 @@ -[TimeManager] +[TimeLoop] DtInitial = 60 # [s] TEnd = 3600 # [s] @@ -9,7 +9,7 @@ Cells = 50 10 [Problem] Name = infiltration3pbox -[] +[SpatialParams] permeability = 1.e-11 # m^2 porosity = 0.40 vanGenuchtenAlpha = 0.0005 @@ -18,3 +18,5 @@ vanGenuchtenN = 4.0 [Output] PlotFluidMatrixInteractions = false +[Newton] +MaxRelativeShift = 1e-4 diff --git a/test/porousmediumflow/3p/implicit/test_box3pniconduction.cc b/test/porousmediumflow/3p/implicit/test_box3pniconduction.cc index b0bfeb26b0625c47a3ad5e34199774fb30a40e4a..ec1298f9a0b403a3b097b653111f3905ea1c2ffe 100644 --- a/test/porousmediumflow/3p/implicit/test_box3pniconduction.cc +++ b/test/porousmediumflow/3p/implicit/test_box3pniconduction.cc @@ -23,7 +23,32 @@ */ #include #include "3pniconductionproblem.hh" -#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include /*! * \brief Provides an interface for customizing error messages associated with @@ -41,17 +66,194 @@ void usage(const char *progName, const std::string &errorMsg) errorMessageOut += " [options]\n"; errorMessageOut += errorMsg; errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" - "\t-TimeManager.TEnd End of the simulation [s] \n" - "\t-TimeManager.DtInitial Initial timestep size [s] \n" - "\t-Grid.File Name of the file containing the grid \n" - "\t definition in DGF format\n"; + "\t-TimeManager.TEnd End of the simulation [s] \n" + "\t-TimeManager.DtInitial Initial timestep size [s] \n" + "\t-Grid.File Name of the file containing the grid \n" + "\t definition in DGF format\n"; + std::cout << errorMessageOut << "\n"; } } -int main(int argc, char** argv) +int main(int argc, char** argv) try +{ + + using namespace Dumux; + + // define the type tag for this problem + using TypeTag = TTAG(ThreePNIConductionBoxProblem); + + //////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////// + + // 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) + using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator); + GridCreator::makeGrid(); + GridCreator::loadBalance(); + + //////////////////////////////////////////////////////////// + // run instationary non-linear problem on this grid + //////////////////////////////////////////////////////////// + + // we compute on the leaf grid view + const auto& leafGridView = GridCreator::grid().leafGridView(); + + // create the finite volume grid geometry + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + auto fvGridGeometry = std::make_shared(leafGridView); + fvGridGeometry->update(); + + // the problem (initial and boundary conditions) + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + auto problem = std::make_shared(fvGridGeometry); + + // the solution vector + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + SolutionVector x(leafGridView.size(GridView::dimension)); + problem->applyInitialSolution(x); + auto xOld = x; + + // the grid variables + using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables); + auto gridVariables = std::make_shared(problem, fvGridGeometry); + gridVariables->init(x, xOld); + + // get some time loop parameters + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + const auto tEnd = getParam("TimeLoop.TEnd"); + const auto maxDivisions = getParam("TimeLoop.MaxTimeStepDivisions"); + const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); + auto dt = getParam("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("TimeLoop.Restart"); + + // intialize the vtk output module + using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields); + VtkOutputModule vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name()); + VtkOutputFields::init(vtkWriter); //! Add model specific output fields + vtkWriter.addField(problem->getExactTemperature(), "temperatureExact"); + vtkWriter.write(0.0); + // output every vtkOutputInterval time step + const auto vtkOutputInterval = getParam("Problem.OutputInterval"); + + // instantiate time loop + auto timeLoop = std::make_shared>(restartTime, dt, tEnd); + timeLoop->setMaxTimeStepSize(maxDt); + + // the assembler with time loop for instationary problem + using Assembler = FVAssembler; + auto assembler = std::make_shared(problem, fvGridGeometry, gridVariables, timeLoop); + + // the linear solver + using LinearSolver = Dumux::AMGBackend; + auto linearSolver = std::make_shared(leafGridView, fvGridGeometry->vertexMapper()); + + // the non-linear solver + using NewtonController = Dumux::NewtonController; + using NewtonMethod = Dumux::NewtonMethod; + auto newtonController = std::make_shared(leafGridView.comm(), timeLoop); + NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver); + + // time loop + timeLoop->start(); do + { + // set previous solution for storage evaluations + assembler->setPreviousSolution(xOld); + + // try solving the non-linear system + for (int i = 0; i < maxDivisions; ++i) + { + // linearize & solve + auto converged = nonLinearSolver.solve(x); + + if (converged) + break; + + if (!converged && i == maxDivisions-1) + DUNE_THROW(Dune::MathError, + "Newton solver didn't converge after " + << maxDivisions + << " time-step divisions. dt=" + << timeLoop->timeStepSize() + << ".\nThe solutions of the current and the previous time steps " + << "have been saved to restart files."); + } + + // compute the new analytical temperature field for the output + problem->updateExactTemperature(x, 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(); + + // report statistics of this time step + timeLoop->reportTimeStep(); + + // set new dt as suggested by newton controller + timeLoop->setTimeStepSize(newtonController->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; + +} +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 (...) { - typedef TTAG(ThreePNIConductionBoxProblem) ProblemTypeTag; - return Dumux::start(argc, argv, usage); + std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; + return 4; } diff --git a/test/porousmediumflow/3p/implicit/test_box3pniconduction.input b/test/porousmediumflow/3p/implicit/test_box3pniconduction.input index 2492431958113698456d3f5169e0ed5339560ccd..cddf9d588eb13142e985d6b7797882fd1d6a5f61 100644 --- a/test/porousmediumflow/3p/implicit/test_box3pniconduction.input +++ b/test/porousmediumflow/3p/implicit/test_box3pniconduction.input @@ -1,4 +1,4 @@ -[TimeManager] +[TimeLoop] DtInitial = 1 # [s] TEnd = 1e5 # [s] MaxTimeStepSize = 1e10 diff --git a/test/porousmediumflow/3p/implicit/test_box3pniconvection.cc b/test/porousmediumflow/3p/implicit/test_box3pniconvection.cc index 4b6f12df0276bd5eb334bbd33662c0ccd3a6d40c..f6ea0553f8f5e6aa8a3627c35ca1ba4d90f3b683 100644 --- a/test/porousmediumflow/3p/implicit/test_box3pniconvection.cc +++ b/test/porousmediumflow/3p/implicit/test_box3pniconvection.cc @@ -23,7 +23,32 @@ */ #include #include "3pniconvectionproblem.hh" -#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include /*! * \brief Provides an interface for customizing error messages associated with @@ -41,17 +66,192 @@ void usage(const char *progName, const std::string &errorMsg) errorMessageOut += " [options]\n"; errorMessageOut += errorMsg; errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" - "\t-TimeManager.TEnd End of the simulation [s] \n" - "\t-TimeManager.DtInitial Initial timestep size [s] \n" - "\t-Grid.File Name of the file containing the grid \n" - "\t definition in DGF format\n"; + "\t-TimeManager.TEnd End of the simulation [s] \n" + "\t-TimeManager.DtInitial Initial timestep size [s] \n" + "\t-Grid.File Name of the file containing the grid \n" + "\t definition in DGF format\n"; + std::cout << errorMessageOut << "\n"; } } -int main(int argc, char** argv) +int main(int argc, char** argv) try +{ + + using namespace Dumux; + + // define the type tag for this problem + using TypeTag = TTAG(ThreePNIConvectionBoxProblem); + + //////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////// + + // 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) + using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator); + GridCreator::makeGrid(); + GridCreator::loadBalance(); + + //////////////////////////////////////////////////////////// + // run instationary non-linear problem on this grid + //////////////////////////////////////////////////////////// + + // we compute on the leaf grid view + const auto& leafGridView = GridCreator::grid().leafGridView(); + + // create the finite volume grid geometry + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + auto fvGridGeometry = std::make_shared(leafGridView); + fvGridGeometry->update(); + + // the problem (initial and boundary conditions) + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + auto problem = std::make_shared(fvGridGeometry); + + // the solution vector + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + SolutionVector x(leafGridView.size(GridView::dimension)); + problem->applyInitialSolution(x); + auto xOld = x; + + // the grid variables + using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables); + auto gridVariables = std::make_shared(problem, fvGridGeometry); + gridVariables->init(x, xOld); + + // get some time loop parameters + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + const auto tEnd = getParam("TimeLoop.TEnd"); + const auto maxDivisions = getParam("TimeLoop.MaxTimeStepDivisions"); + const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); + auto dt = getParam("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("TimeLoop.Restart"); + + // intialize the vtk output module + using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields); + VtkOutputModule vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name()); + VtkOutputFields::init(vtkWriter); //! Add model specific output fields + vtkWriter.addField(problem->getExactTemperature(), "temperatureExact"); + vtkWriter.write(0.0); + // output every vtkOutputInterval time step + const int vtkOutputInterval = getParam("Problem.OutputInterval"); + + // instantiate time loop + auto timeLoop = std::make_shared>(restartTime, dt, tEnd); + timeLoop->setMaxTimeStepSize(maxDt); + + // the assembler with time loop for instationary problem + using Assembler = FVAssembler; + auto assembler = std::make_shared(problem, fvGridGeometry, gridVariables, timeLoop); + + // the linear solver + using LinearSolver = Dumux::AMGBackend; + auto linearSolver = std::make_shared(leafGridView, fvGridGeometry->vertexMapper()); + + // the non-linear solver + using NewtonController = Dumux::NewtonController; + using NewtonMethod = Dumux::NewtonMethod; + auto newtonController = std::make_shared(leafGridView.comm(), timeLoop); + NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver); + + // time loop + timeLoop->start(); do + { + // set previous solution for storage evaluations + assembler->setPreviousSolution(xOld); + + // try solving the non-linear system + for (int i = 0; i < maxDivisions; ++i) + { + // linearize & solve + auto converged = nonLinearSolver.solve(x); + + if (converged) + break; + + if (!converged && i == maxDivisions-1) + DUNE_THROW(Dune::MathError, + "Newton solver didn't converge after " + << maxDivisions + << " time-step divisions. dt=" + << timeLoop->timeStepSize() + << ".\nThe solutions of the current and the previous time steps " + << "have been saved to restart files."); + } + + // compute the new analytical temperature field for the output + problem->updateExactTemperature(x, 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(); + + // report statistics of this time step + timeLoop->reportTimeStep(); + + // set new dt as suggested by newton controller + timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize())); + + 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; + +} +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 (...) { - typedef TTAG(ThreePNIConvectionBoxProblem) ProblemTypeTag; - return Dumux::start(argc, argv, usage); + std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; + return 4; } diff --git a/test/porousmediumflow/3p/implicit/test_box3pniconvection.input b/test/porousmediumflow/3p/implicit/test_box3pniconvection.input index 5dc7d07dbdc79866c820bf651548a6a26481f09d..74c95f92a70ec36847d50733bb961979ed598397 100644 --- a/test/porousmediumflow/3p/implicit/test_box3pniconvection.input +++ b/test/porousmediumflow/3p/implicit/test_box3pniconvection.input @@ -1,4 +1,4 @@ -[TimeManager] +[TimeLoop] DtInitial = 1 # [s] TEnd = 3e4 # [s] MaxTimeStepSize = 1e3 diff --git a/test/porousmediumflow/3p/implicit/test_cc3p.cc b/test/porousmediumflow/3p/implicit/test_cc3p.cc index 7ad8e04e371041dc68ffd728e3ded8845685b187..d7c1e3d44e90bf2fad48fef1adab26787b939e25 100644 --- a/test/porousmediumflow/3p/implicit/test_cc3p.cc +++ b/test/porousmediumflow/3p/implicit/test_cc3p.cc @@ -23,7 +23,31 @@ */ #include #include "infiltration3pproblem.hh" -#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include /*! * \brief Provides an interface for customizing error messages associated with @@ -40,20 +64,191 @@ void usage(const char *progName, const std::string &errorMsg) errorMessageOut += progName; errorMessageOut += " [options]\n"; errorMessageOut += errorMsg; - errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" - "\t-TimeManager.TEnd End of the simulation [s] \n" - "\t-TimeManager.DtInitial Initial timestep size [s] \n" - "\t-Grid.File Name of the file containing the grid \n" - "\t definition in DGF format\n"; + 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) +int main(int argc, char** argv) try { - typedef TTAG(InfiltrationThreePCCProblem) ProblemTypeTag; - return Dumux::start(argc, argv, usage); -} + using namespace Dumux; + + // define the type tag for this problem + using TypeTag = TTAG(InfiltrationThreePCCProblem); + + // 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) + using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator); + GridCreator::makeGrid(); + GridCreator::loadBalance(); + + //////////////////////////////////////////////////////////// + // run instationary non-linear problem on this grid + //////////////////////////////////////////////////////////// + + // we compute on the leaf grid view + const auto& leafGridView = GridCreator::grid().leafGridView(); + + // create the finite volume grid geometry + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + auto fvGridGeometry = std::make_shared(leafGridView); + fvGridGeometry->update(); + + // the problem (initial and boundary conditions) + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + auto problem = std::make_shared(fvGridGeometry); + + // the solution vector + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + SolutionVector x(leafGridView.size(0)); + problem->applyInitialSolution(x); + auto xOld = x; + + // the grid variables + using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables); + auto gridVariables = std::make_shared(problem, fvGridGeometry); + gridVariables->init(x, xOld); + + // get some time loop parameters + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + const auto tEnd = getParam("TimeLoop.TEnd"); + const auto maxDivisions = getParam("TimeLoop.MaxTimeStepDivisions"); + const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); + auto dt = getParam("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("TimeLoop.Restart"); + + // intialize the vtk output module + using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields); + VtkOutputModule vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name()); + VtkOutputFields::init(vtkWriter); //! Add model specific output fields + vtkWriter.write(0.0); + + // instantiate time loop + auto timeLoop = std::make_shared>(restartTime, dt, tEnd); + timeLoop->setMaxTimeStepSize(maxDt); + + // the assembler with time loop for instationary problem + using Assembler = FVAssembler; + auto assembler = std::make_shared(problem, fvGridGeometry, gridVariables, timeLoop); + + // the linear solver + using LinearSolver = UMFPackBackend; + auto linearSolver = std::make_shared(); + // the non-linear solver + using NewtonController = Dumux::NewtonController; + using NewtonMethod = Dumux::NewtonMethod; + auto newtonController = std::make_shared(leafGridView.comm(), timeLoop); + NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver); + + // time loop + timeLoop->start(); do + { + // set previous solution for storage evaluations + assembler->setPreviousSolution(xOld); + + // set the end of the next time step as time in the problem to control + // the boundary conditions for the implicit Euler scheme + problem->setTime(timeLoop->time()+timeLoop->timeStepSize()); + + // try solving the non-linear system + for (int i = 0; i < maxDivisions; ++i) + { + // linearize & solve + auto converged = nonLinearSolver.solve(x); + + if (converged) + break; + + if (!converged && i == maxDivisions-1) + DUNE_THROW(Dune::MathError, + "Newton solver didn't converge after " + << maxDivisions + << " time-step divisions. dt=" + << timeLoop->timeStepSize() + << ".\nThe solutions of the current and the previous time steps " + << "have been saved to restart files."); + } + + // 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(); + + // set new dt as suggested by newton controller + timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize())); + + // write vtk output + 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; +} diff --git a/test/porousmediumflow/3p/implicit/test_cc3p.input b/test/porousmediumflow/3p/implicit/test_cc3p.input index 764e1a3f0c4ab0036e27c25bc2f9177e9a33673e..2fc97e822c6883ea43f696ddb6d15cca2a256f48 100644 --- a/test/porousmediumflow/3p/implicit/test_cc3p.input +++ b/test/porousmediumflow/3p/implicit/test_cc3p.input @@ -1,4 +1,4 @@ -[TimeManager] +[TimeLoop] DtInitial = 60 # [s] TEnd = 3600 # [s] @@ -9,7 +9,7 @@ Cells = 50 10 [Problem] Name = infiltration3pcc -[] +[SpatialParams] permeability = 1.e-11 # m^2 porosity = 0.40 vanGenuchtenAlpha = 0.0005 @@ -18,3 +18,5 @@ vanGenuchtenN = 4.0 [Output] PlotFluidMatrixInteractions = false +[Newton] +MaxRelativeShift = 1e-4 diff --git a/test/porousmediumflow/3p/implicit/test_cc3pniconduction.cc b/test/porousmediumflow/3p/implicit/test_cc3pniconduction.cc index 6c1fc640e642a690be25416d747a657a7f54acc6..df082b77967ba9af093c3f0af5780434826f215e 100644 --- a/test/porousmediumflow/3p/implicit/test_cc3pniconduction.cc +++ b/test/porousmediumflow/3p/implicit/test_cc3pniconduction.cc @@ -23,7 +23,31 @@ */ #include #include "3pniconductionproblem.hh" -#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include /*! * \brief Provides an interface for customizing error messages associated with @@ -40,18 +64,193 @@ void usage(const char *progName, const std::string &errorMsg) errorMessageOut += progName; errorMessageOut += " [options]\n"; errorMessageOut += errorMsg; - errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" - "\t-TimeManager.TEnd End of the simulation [s] \n" - "\t-TimeManager.DtInitial Initial timestep size [s] \n" - "\t-Grid.File Name of the file containing the grid \n" - "\t definition in DGF format\n"; + 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) +int main(int argc, char** argv) try +{ + using namespace Dumux; + + // define the type tag for this problem + using TypeTag = TTAG(ThreePNIConductionCCProblem); + +// 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) + using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator); + GridCreator::makeGrid(); + GridCreator::loadBalance(); + + //////////////////////////////////////////////////////////// + // run instationary non-linear problem on this grid + //////////////////////////////////////////////////////////// + + // we compute on the leaf grid view + const auto& leafGridView = GridCreator::grid().leafGridView(); + + // create the finite volume grid geometry + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + auto fvGridGeometry = std::make_shared(leafGridView); + fvGridGeometry->update(); + + // the problem (initial and boundary conditions) + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + auto problem = std::make_shared(fvGridGeometry); + + // the solution vector + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + SolutionVector x(leafGridView.size(0)); + problem->applyInitialSolution(x); + auto xOld = x; + + // the grid variables + using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables); + auto gridVariables = std::make_shared(problem, fvGridGeometry); + gridVariables->init(x, xOld); + + // get some time loop parameters + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + const auto tEnd = getParam("TimeLoop.TEnd"); + const auto maxDivisions = getParam("TimeLoop.MaxTimeStepDivisions"); + const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); + auto dt = getParam("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("TimeLoop.Restart"); + + // intialize the vtk output module + using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields); + VtkOutputModule vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name()); + VtkOutputFields::init(vtkWriter); //! Add model specific output fields + vtkWriter.addField(problem->getExactTemperature(), "temperatureExact"); + vtkWriter.write(0.0); + // output every vtkOutputInterval time step + const auto vtkOutputInterval = getParam("Problem.OutputInterval"); + + // instantiate time loop + auto timeLoop = std::make_shared>(restartTime, dt, tEnd); + timeLoop->setMaxTimeStepSize(maxDt); + + // the assembler with time loop for instationary problem + using Assembler = FVAssembler; + auto assembler = std::make_shared(problem, fvGridGeometry, gridVariables, timeLoop); + + // the linear solver + using LinearSolver = Dumux::AMGBackend; + auto linearSolver = std::make_shared(leafGridView, fvGridGeometry->elementMapper()); + + // the non-linear solver + using NewtonController = Dumux::NewtonController; + using NewtonMethod = Dumux::NewtonMethod; + auto newtonController = std::make_shared(leafGridView.comm(), timeLoop); + NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver); + + // time loop + timeLoop->start(); do + { + // set previous solution for storage evaluations + assembler->setPreviousSolution(xOld); + + // try solving the non-linear system + for (int i = 0; i < maxDivisions; ++i) + { + // linearize & solve + auto converged = nonLinearSolver.solve(x); + + if (converged) + break; + + if (!converged && i == maxDivisions-1) + DUNE_THROW(Dune::MathError, + "Newton solver didn't converge after " + << maxDivisions + << " time-step divisions. dt=" + << timeLoop->timeStepSize() + << ".\nThe solutions of the current and the previous time steps " + << "have been saved to restart files."); + } + + // update the exact time temperature + problem->updateExactTemperature(x, 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(); + + // report statistics of this time step + timeLoop->reportTimeStep(); + + // set new dt as suggested by newton controller + timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize())); + + 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 (...) { - typedef TTAG(ThreePNIConductionCCProblem) ProblemTypeTag; - return Dumux::start(argc, argv, usage); + std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; + return 4; } diff --git a/test/porousmediumflow/3p/implicit/test_cc3pniconduction.input b/test/porousmediumflow/3p/implicit/test_cc3pniconduction.input index 42fb551f9b6a4a78d44272188c88321a03f1e661..11dff4a2c571408c3a9079969f8aa1c93e9b7b72 100644 --- a/test/porousmediumflow/3p/implicit/test_cc3pniconduction.input +++ b/test/porousmediumflow/3p/implicit/test_cc3pniconduction.input @@ -1,4 +1,4 @@ -[TimeManager] +[TimeLoop] DtInitial = 1 # [s] TEnd = 1e5 # [s] MaxTimeStepSize = 1e10 diff --git a/test/porousmediumflow/3p/implicit/test_cc3pniconvection.cc b/test/porousmediumflow/3p/implicit/test_cc3pniconvection.cc index dce33c08ccc8e3cbc7044a0ecb07dd64509e2c73..f00dacb3d6fb3cad808dd5aa0457967a211f7cba 100644 --- a/test/porousmediumflow/3p/implicit/test_cc3pniconvection.cc +++ b/test/porousmediumflow/3p/implicit/test_cc3pniconvection.cc @@ -23,7 +23,31 @@ */ #include #include "3pniconvectionproblem.hh" -#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include /*! * \brief Provides an interface for customizing error messages associated with @@ -40,18 +64,194 @@ void usage(const char *progName, const std::string &errorMsg) errorMessageOut += progName; errorMessageOut += " [options]\n"; errorMessageOut += errorMsg; - errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" - "\t-TimeManager.TEnd End of the simulation [s] \n" - "\t-TimeManager.DtInitial Initial timestep size [s] \n" - "\t-Grid.File Name of the file containing the grid \n" - "\t definition in DGF format\n"; + 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) +int main(int argc, char** argv) try +{ + using namespace Dumux; + + // define the type tag for this problem + using TypeTag = TTAG(ThreePNIConvectionCCProblem); + + // 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) + using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator); + GridCreator::makeGrid(); + GridCreator::loadBalance(); + + //////////////////////////////////////////////////////////// + // run instationary non-linear problem on this grid + //////////////////////////////////////////////////////////// + + // we compute on the leaf grid view + const auto& leafGridView = GridCreator::grid().leafGridView(); + + // create the finite volume grid geometry + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + auto fvGridGeometry = std::make_shared(leafGridView); + fvGridGeometry->update(); + + // the problem (initial and boundary conditions) + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + auto problem = std::make_shared(fvGridGeometry); + + // the solution vector + using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); + SolutionVector x(leafGridView.size(0)); + problem->applyInitialSolution(x); + auto xOld = x; + + // the grid variables + using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables); + auto gridVariables = std::make_shared(problem, fvGridGeometry); + gridVariables->init(x, xOld); + + // get some time loop parameters + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + const auto tEnd = getParam("TimeLoop.TEnd"); + const auto maxDivisions = getParam("TimeLoop.MaxTimeStepDivisions"); + const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); + auto dt = getParam("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("TimeLoop.Restart"); + + // intialize the vtk output module + using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields); + VtkOutputModule vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name()); + VtkOutputFields::init(vtkWriter); //! Add model specific output fields + vtkWriter.addField(problem->getExactTemperature(), "temperatureExact"); + vtkWriter.write(0.0); + // output every vtkOutputInterval time step + const auto vtkOutputInterval = getParam("Problem.OutputInterval"); + + // instantiate time loop + auto timeLoop = std::make_shared>(restartTime, dt, tEnd); + timeLoop->setMaxTimeStepSize(maxDt); + + // the assembler with time loop for instationary problem + using Assembler = FVAssembler; + auto assembler = std::make_shared(problem, fvGridGeometry, gridVariables, timeLoop); + + // the linear solver + using LinearSolver = Dumux::AMGBackend; + auto linearSolver = std::make_shared(leafGridView, fvGridGeometry->elementMapper()); + + // the non-linear solver + using NewtonController = Dumux::NewtonController; + using NewtonMethod = Dumux::NewtonMethod; + auto newtonController = std::make_shared(leafGridView.comm(), timeLoop); + NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver); + + // time loop + timeLoop->start(); do + { + // set previous solution for storage evaluations + assembler->setPreviousSolution(xOld); + + // try solving the non-linear system + for (int i = 0; i < maxDivisions; ++i) + { + // linearize & solve + auto converged = nonLinearSolver.solve(x); + + if (converged) + break; + + if (!converged && i == maxDivisions-1) + DUNE_THROW(Dune::MathError, + "Newton solver didn't converge after " + << maxDivisions + << " time-step divisions. dt=" + << timeLoop->timeStepSize() + << ".\nThe solutions of the current and the previous time steps " + << "have been saved to restart files."); + } + + // compute the new analytical temperature field for the output + problem->updateExactTemperature(x, 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(); + + // report statistics of this time step + timeLoop->reportTimeStep(); + + // set new dt as suggested by newton controller + timeLoop->setTimeStepSize(newtonController->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 (...) { - typedef TTAG(ThreePNIConvectionCCProblem) ProblemTypeTag; - return Dumux::start(argc, argv, usage); + std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; + return 4; } diff --git a/test/porousmediumflow/3p/implicit/test_cc3pniconvection.input b/test/porousmediumflow/3p/implicit/test_cc3pniconvection.input index 2cd99fed09f035001799b6a185290f3904d4e5b9..fc7f8894399b094963f19235f7dbfec8cdaabeab 100644 --- a/test/porousmediumflow/3p/implicit/test_cc3pniconvection.input +++ b/test/porousmediumflow/3p/implicit/test_cc3pniconvection.input @@ -1,4 +1,4 @@ -[TimeManager] +[TimeLoop] DtInitial = 1 # [s] TEnd = 3e4 # [s] MaxTimeStepSize = 1e3 @@ -15,3 +15,4 @@ EnableGravity = 0 # disable gravity [Vtk] AddVelocity = 1 #Enable velocity output +