diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh index b388c0710757ca90419df6e355388a8526b92c57..1aa128345d0b64d995f89faae3dcbf6d54f35cc1 100644 --- a/dumux/common/properties.hh +++ b/dumux/common/properties.hh @@ -147,6 +147,8 @@ NEW_PROP_TAG(NumSPhases); NEW_PROP_TAG(NonMineralizationVtkOutputFields); NEW_PROP_TAG(NonMineralizationVolumeVariables); +NEW_PROP_TAG(UseConstraintSolver); //! Determines whether the constraint solver should be used# + ///////////////////////////////////////////////////////////// // non-isothermal porous medium flow models ///////////////////////////////////////////////////////////// diff --git a/dumux/common/timeloop.hh b/dumux/common/timeloop.hh index 759cfec43285331e354ce6e008b943175ac24a98..7416df37e3e73790f8dbc1e6133bf8cfc4859c6c 100644 --- a/dumux/common/timeloop.hh +++ b/dumux/common/timeloop.hh @@ -365,14 +365,14 @@ public: //! Check point management, TimeLoop::isCheckPoint() has to be called after this! // if we reached a periodic check point - if (periodicCheckPoints_ && Dune::FloatCmp::eq(this->time(), lastPeriodicCheckPoint_ + deltaPeriodicCheckPoint_, 1e-3*this->timeStepSize())) + if (periodicCheckPoints_ && Dune::FloatCmp::eq(this->time(), lastPeriodicCheckPoint_ + deltaPeriodicCheckPoint_, 1e-8*this->timeStepSize())) { lastPeriodicCheckPoint_ += deltaPeriodicCheckPoint_; isCheckPoint_ = true; } // or a manually set check point - else if (!checkPoints_.empty() && Dune::FloatCmp::eq(this->time(), checkPoints_.front(), 1e-3*this->timeStepSize())) + else if (!checkPoints_.empty() && Dune::FloatCmp::eq(this->time(), checkPoints_.front(), 1e-8*this->timeStepSize())) { checkPoints_.pop(); isCheckPoint_ = true; diff --git a/dumux/porousmediumflow/3p3c/implicit/CMakeLists.txt b/dumux/porousmediumflow/3p3c/implicit/CMakeLists.txt index 8df41cae3571ff308f6d9f199c2198a74744618b..296992e3c515c9f96e2cd8e20b48d28ebddff17c 100644 --- a/dumux/porousmediumflow/3p3c/implicit/CMakeLists.txt +++ b/dumux/porousmediumflow/3p3c/implicit/CMakeLists.txt @@ -1,12 +1,10 @@ #install headers install(FILES -fluxvariables.hh indices.hh localresidual.hh model.hh -newtoncontroller.hh -properties.hh -propertydefaults.hh +primaryvariableswitch.hh volumevariables.hh +vtkoutputfields.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/3p3c/implicit) diff --git a/dumux/porousmediumflow/3p3c/implicit/indices.hh b/dumux/porousmediumflow/3p3c/implicit/indices.hh index cd832f5ec67605326b1a5d2be03546130b3c90c8..ec7525d8ddf39023f3b7c0dc64fc71ed3e23ac26 100644 --- a/dumux/porousmediumflow/3p3c/implicit/indices.hh +++ b/dumux/porousmediumflow/3p3c/implicit/indices.hh @@ -24,7 +24,7 @@ #ifndef DUMUX_3P3C_INDICES_HH #define DUMUX_3P3C_INDICES_HH -#include "properties.hh" +#include namespace Dumux { @@ -85,6 +85,6 @@ public: static const int contiGEqIdx = conti0EqIdx + gCompIdx; //!< index of the mass conservation equation for the air component }; -} +} // end namespace Dumux #endif diff --git a/dumux/porousmediumflow/3p3c/implicit/localresidual.hh b/dumux/porousmediumflow/3p3c/implicit/localresidual.hh index 5b4282d8b89740783ab9d549368b816683b06eb4..488868f327fb27aef3ca79339254bb286f6391eb 100644 --- a/dumux/porousmediumflow/3p3c/implicit/localresidual.hh +++ b/dumux/porousmediumflow/3p3c/implicit/localresidual.hh @@ -25,7 +25,7 @@ #ifndef DUMUX_3P3C_LOCAL_RESIDUAL_HH #define DUMUX_3P3C_LOCAL_RESIDUAL_HH -#include "properties.hh" +#include namespace Dumux { @@ -42,10 +42,12 @@ class ThreePThreeCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual { using ParentType = typename GET_PROP_TYPE(TypeTag, BaseLocalResidual); using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual); + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); + using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector); using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables); using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache); using Indices = typename GET_PROP_TYPE(TypeTag, Indices); @@ -76,6 +78,8 @@ class ThreePThreeCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual public: + using ParentType::ParentType; + /*! * \brief Evaluate the amount of all conservation quantities * (e.g. phase mass) within a sub-control volume. @@ -87,10 +91,11 @@ public: * \param scvIdx The SCV (sub-control-volume) index * \param usePrevSol Evaluate function with solution of current or previous time step */ - PrimaryVariables computeStorage(const SubControlVolume& scv, - const VolumeVariables& volVars) const + ResidualVector computeStorage(const Problem& problem, + const SubControlVolume& scv, + const VolumeVariables& volVars) const { - PrimaryVariables storage(0.0); + ResidualVector storage(0.0); // compute storage term of all components within all phases for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) @@ -123,17 +128,18 @@ public: * \param onBoundary A boolean variable to specify whether the flux variables * are calculated for interior SCV faces or boundary faces, default=false */ - PrimaryVariables computeFlux(const Element& element, - const FVElementGeometry& fvGeometry, - const ElementVolumeVariables& elemVolVars, - const SubControlVolumeFace& scvf, - const ElementFluxVariablesCache& elemFluxVarsCache) + ResidualVector computeFlux(const Problem& problem, + const Element& element, + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf, + const ElementFluxVariablesCache& elemFluxVarsCache) const { FluxVariables fluxVars; - fluxVars.init(this->problem(), element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache); + fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache); // get upwind weights into local scope - PrimaryVariables flux(0.0); + ResidualVector flux(0.0); // advective fluxes for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) @@ -155,23 +161,23 @@ public: //! Add diffusive energy fluxes. For isothermal model the contribution is zero. EnergyLocalResidual::heatConductionFlux(flux, fluxVars); + const auto diffusionFluxesWPhase = fluxVars.molecularDiffusionFlux(wPhaseIdx); + // diffusive fluxes - Scalar jGW = fluxVars.molecularDiffusionFlux(wPhaseIdx, gCompIdx); - Scalar jNW = fluxVars.molecularDiffusionFlux(wPhaseIdx, nCompIdx); + Scalar jGW = diffusionFluxesWPhase[gCompIdx]; + Scalar jNW = diffusionFluxesWPhase[nCompIdx]; Scalar jWW = -(jGW+jNW); - Scalar jWG = fluxVars.molecularDiffusionFlux(gPhaseIdx, wCompIdx); - Scalar jNG = fluxVars.molecularDiffusionFlux(gPhaseIdx, nCompIdx); - Scalar jGG = -(jWG+jNG); + const auto diffusionFluxesGPhase = fluxVars.molecularDiffusionFlux(gPhaseIdx); - // Scalar jWN = fluxVars.molecularDiffusion().flux(nPhaseIdx, wCompIdx); - // Scalar jGN = fluxVars.molecularDiffusion().flux(nPhaseIdx, gCompIdx); - // Scalar jNN = -(jGN+jWN); + Scalar jWG = diffusionFluxesGPhase[wCompIdx]; + Scalar jNG = diffusionFluxesGPhase[nCompIdx]; + Scalar jGG = -(jWG+jNG); // At the moment we do not consider diffusion in the NAPL phase - Scalar jWN = 0; - Scalar jGN = 0; - Scalar jNN = 0; + Scalar jWN = 0.0; + Scalar jGN = 0.0; + Scalar jNN = 0.0; flux[conti0EqIdx] += jWW+jWG+jWN; flux[conti1EqIdx] += jNW+jNG+jNN; diff --git a/dumux/porousmediumflow/3p3c/implicit/model.hh b/dumux/porousmediumflow/3p3c/implicit/model.hh index a6553ff5d9f4e7b7e20f4a31cfb13bf06ed7811c..4c9ff253f546a6815b055f93963bb958c3fb5d0b 100644 --- a/dumux/porousmediumflow/3p3c/implicit/model.hh +++ b/dumux/porousmediumflow/3p3c/implicit/model.hh @@ -18,23 +18,6 @@ *****************************************************************************/ /*! * \file - * - * \brief Adaption of the fully implicit scheme to the three-phase three-component - * flow model. - * - * The model is designed for simulating three fluid phases with water, gas, and - * a liquid contaminant (NAPL - non-aqueous phase liquid) - */ -#ifndef DUMUX_3P3C_MODEL_HH -#define DUMUX_3P3C_MODEL_HH - -#include -#include "properties.hh" -#include "primaryvariableswitch.hh" - -namespace Dumux -{ -/*! * \ingroup ThreePThreeCModel * \brief Adaption of the fully implicit scheme to the three-phase three-component * flow model. @@ -91,309 +74,159 @@ namespace Dumux *
  • Water and gas phases are present: Primary variables \f$(S_w\f$, \f$x_w^g\f$, \f$p_g)\f$.
    • * */ -template -class ThreePThreeCModel: public GET_PROP_TYPE(TypeTag, BaseModel) + +#ifndef DUMUX_3P3C_MODEL_HH +#define DUMUX_3P3C_MODEL_HH + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "indices.hh" +#include "volumevariables.hh" +#include "vtkoutputfields.hh" +#include "primaryvariableswitch.hh" +#include "localresidual.hh" + +namespace Dumux { +namespace Properties { + +//! The type tags for the implicit three-phase three-component problems +NEW_TYPE_TAG(ThreePThreeC, INHERITS_FROM(PorousMediumFlow)); + +//! The type tags for the corresponding non-isothermal problems +NEW_TYPE_TAG(ThreePThreeCNI, INHERITS_FROM(ThreePThreeC, NonIsothermal)); + +////////////////////////////////////////////////////////////////// +// Property values +////////////////////////////////////////////////////////////////// + +/*! + * \brief Set the property for the number of components. + * + * We just forward the number from the fluid system and use an static + * assert to make sure it is 3. + */ +SET_PROP(ThreePThreeC, NumComponents) { - // the parent class needs to access the variable switch - friend typename GET_PROP_TYPE(TypeTag, BaseModel); - - using ParentType = typename GET_PROP_TYPE(TypeTag, BaseModel); - using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); - using GridView = typename GET_PROP_TYPE(TypeTag, GridView); - using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); - using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables); - using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); - using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); - using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector); - using Indices = typename GET_PROP_TYPE(TypeTag, Indices); - - enum { - dim = GridView::dimension, - dimWorld = GridView::dimensionworld, - - numPhases = GET_PROP_VALUE(TypeTag, NumPhases), - numComponents = GET_PROP_VALUE(TypeTag, NumComponents), - - switch1Idx = Indices::switch1Idx, - switch2Idx = Indices::switch2Idx, - - wPhaseIdx = Indices::wPhaseIdx, - nPhaseIdx = Indices::nPhaseIdx, - gPhaseIdx = Indices::gPhaseIdx, - - wCompIdx = Indices::wCompIdx, - nCompIdx = Indices::nCompIdx, - gCompIdx = Indices::gCompIdx, - - threePhases = Indices::threePhases, - wPhaseOnly = Indices::wPhaseOnly, - gnPhaseOnly = Indices::gnPhaseOnly, - wnPhaseOnly = Indices::wnPhaseOnly, - gPhaseOnly = Indices::gPhaseOnly, - wgPhaseOnly = Indices::wgPhaseOnly - - }; - - using GlobalPosition = Dune::FieldVector; - enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; - enum { dofCodim = isBox ? dim : 0 }; + static const int value = 3; + static_assert(value == GET_PROP_TYPE(TypeTag, FluidSystem)::numComponents, + "Only fluid systems with 3 components are supported by the 3p3c model!"); +}; -public: +/*! + * \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(ThreePThreeC, NumPhases) +{ + static const int value = 3; + static_assert(value == GET_PROP_TYPE(TypeTag, FluidSystem)::numPhases, + "Only fluid systems with 3 phases are supported by the 3p3c model!"); +}; + +//! Set as default that no component mass balance is replaced by the total mass balance +SET_INT_PROP(ThreePThreeC, ReplaceCompEqIdx, GET_PROP_VALUE(TypeTag, NumComponents)); +/*! + * \brief The fluid state which is used by the volume variables to + * store the thermodynamic state. This should be chosen + * appropriately for the model ((non-)isothermal, equilibrium, ...). + * This can be done in the problem. + */ +SET_PROP(ThreePThreeC, FluidState){ + private: + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + public: + typedef CompositionalFluidState type; +}; + +SET_INT_PROP(ThreePThreeC, NumEq, 3); //!< set the number of equations to 3 + +//! The local residual function of the conservation equations +SET_TYPE_PROP(ThreePThreeC, LocalResidual, ThreePThreeCLocalResidual); + +//! Enable advection +SET_BOOL_PROP(ThreePThreeC, EnableAdvection, true); + +//! Enable molecular diffusion +SET_BOOL_PROP(ThreePThreeC, EnableMolecularDiffusion, true); + +//! Isothermal model by default +SET_BOOL_PROP(ThreePThreeC, EnableEnergyBalance, false); + +//! The primary variable switch for the 3p3c model +SET_TYPE_PROP(ThreePThreeC, PrimaryVariableSwitch, ThreePThreeCPrimaryVariableSwitch); + +//! The primary variables vector for the 3p3c model +SET_TYPE_PROP(ThreePThreeC, PrimaryVariables, SwitchablePrimaryVariables); + +//! the VolumeVariables property +SET_TYPE_PROP(ThreePThreeC, VolumeVariables, ThreePThreeCVolumeVariables); - /*! - * \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); }); - - - for (int i = 0; i < numPhases; ++i) - for (int j = 0; j < numComponents; ++j) - vtkOutputModule.addSecondaryVariable("x^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(i), - [i,j](const VolumeVariables& v){ return v.moleFraction(i,j); }); - - 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(); }); - } - - /*! - * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write. - */ - template - void addVtkOutputFields(VtkOutputModule& outputModule) const - { - auto& phasePresence = outputModule.createScalarField("phase presence", dofCodim); - for (std::size_t i = 0; i < phasePresence.size(); ++i) - phasePresence[i] = this->curSol()[i].state(); - } - - /*! - * \brief One Newton iteration was finished. - * \param uCurrent The solution after the current Newton iteration - */ - template - typename std::enable_if::type - newtonEndStep() - { - // \todo resize volvars vector if grid was adapted - - // update the variable switch - switchFlag_ = priVarSwitch_().update(this->problem_(), this->curSol()); - - // update the secondary variables if global caching is enabled - // \note we only updated if phase presence changed as the volume variables - // are already updated once by the switch - for (const auto& element : elements(this->problem_().gridView())) - { - // make sure FVElementGeometry & vol vars are bound to the element - auto fvGeometry = localView(this->fvGridGeometry()); - fvGeometry.bindElement(element); - - if (switchFlag_) - { - for (auto&& scv : scvs(fvGeometry)) - { - auto dofIdxGlobal = scv.dofIndex(); - if (priVarSwitch_().wasSwitched(dofIdxGlobal)) - { - const auto eIdx = this->problem_().elementMapper().index(element); - const auto elemSol = this->elementSolution(element, this->curSol()); - this->nonConstCurGlobalVolVars().volVars(eIdx, scv.indexInElement()).update(elemSol, - this->problem_(), - element, - scv); - } - } - } - - // handle the boundary volume variables for cell-centered models - if(!isBox) - { - for (auto&& scvf : scvfs(fvGeometry)) - { - // if we are not on a boundary, skip the rest - if (!scvf.boundary()) - continue; - - // check if boundary is a pure dirichlet boundary - const auto bcTypes = this->problem_().boundaryTypes(element, scvf); - if (bcTypes.hasOnlyDirichlet()) - { - const auto insideScvIdx = scvf.insideScvIdx(); - const auto& insideScv = fvGeometry.scv(insideScvIdx); - const auto elemSol = ElementSolutionVector{this->problem_().dirichlet(element, scvf)}; - - this->nonConstCurGlobalVolVars().volVars(scvf.outsideScvIdx(), 0/*indexInElement*/).update(elemSol, this->problem_(), element, insideScv); - } - } - } - } - } - - /*! - * \brief Compute the total storage inside one phase of all - * conservation quantities. - * - * \param storage Contains the storage of each component for one phase - * \param phaseIdx The phase index - */ - void globalPhaseStorage(PrimaryVariables &storage, const int phaseIdx) - { - storage = 0; - - for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior)) - { - this->localResidual().evalPhaseStorage(element, phaseIdx); - - for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i) - storage += this->localResidual().storageTerm()[i]; - } - if (this->gridView_().comm().size() > 1) - storage = this->gridView_().comm().sum(storage); - } - - - /*! - * \brief Called by the update() method if applying the newton - * \brief One Newton iteration was finished. - * \param uCurrent The solution after the current Newton iteration - */ - template - typename std::enable_if::type - newtonEndStep() - { - // update the variable switch - switchFlag_ = priVarSwitch_().update(this->problem_(), this->curSol()); - } - - /*! - * \brief Called by the update() method if applying the Newton - * method was unsuccessful. - */ - void updateFailed() - { - ParentType::updateFailed(); - // reset privar switch flag - switchFlag_ = false; - } - - /*! - * \brief Called by the problem if a time integration was - * successful, post processing of the solution is done and the - * result has been written to disk. - * - * This should prepare the model for the next time integration. - */ - void advanceTimeLevel() - { - ParentType::advanceTimeLevel(); - // reset privar switch flag - switchFlag_ = false; - } - - /*! - * \brief Returns true if the primary variables were switched for - * at least one dof after the last timestep. - */ - bool switched() const - { - return switchFlag_; - } - - /*! - * \brief Write the current solution to a restart file. - * - * \param outStream The output stream of one entity for the restart file - * \param entity The entity, either a vertex or an element - */ - template - void serializeEntity(std::ostream &outStream, const Entity &entity) - { - // write primary variables - ParentType::serializeEntity(outStream, entity); - - int dofIdxGlobal = this->dofMapper().index(entity); - - if (!outStream.good()) - DUNE_THROW(Dune::IOError, "Could not serialize entity " << dofIdxGlobal); - - outStream << this->curSol()[dofIdxGlobal].state() << " "; - } - - /*! - * \brief Reads the current solution from a restart file. - * - * \param inStream The input stream of one entity from the restart file - * \param entity The entity, either a vertex or an element - */ - template - void deserializeEntity(std::istream &inStream, const Entity &entity) - { - // read primary variables - ParentType::deserializeEntity(inStream, entity); - - // read phase presence - int dofIdxGlobal = this->dofMapper().index(entity); - - if (!inStream.good()) - DUNE_THROW(Dune::IOError, "Could not deserialize entity " << dofIdxGlobal); - - int phasePresence; - inStream >> phasePresence; - - this->curSol()[dofIdxGlobal].setState(phasePresence); - this->prevSol()[dofIdxGlobal].setState(phasePresence); - } - - const Dumux::ThreePThreeCPrimaryVariableSwitch& priVarSwitch() const - { return switch_; } - -protected: - - Dumux::ThreePThreeCPrimaryVariableSwitch& priVarSwitch_() - { return switch_; } - - /*! - * \brief Applies the initial solution for all vertices of the grid. - * - * \todo the initial condition needs to be unique for - * each vertex. we should think about the API... - */ - void applyInitialSolution_() - { - ParentType::applyInitialSolution_(); - - // initialize the primary variable switch - priVarSwitch_().init(this->problem_()); - } - - //! the class handling the primary variable switch - Dumux::ThreePThreeCPrimaryVariableSwitch switch_; - bool switchFlag_; +//! Determines whether a constraint solver should be used explicitly +SET_BOOL_PROP(ThreePThreeC, UseConstraintSolver, false); + +//! The indices required by the isothermal 3p3c model +SET_TYPE_PROP(ThreePThreeC, Indices, ThreePThreeCIndices); + +//! The spatial parameters to be employed. +//! Use ImplicitSpatialParams by default. +SET_TYPE_PROP(ThreePThreeC, SpatialParams, ImplicitSpatialParams); + +//! The model after Millington (1961) is used for the effective diffusivity +SET_PROP(ThreePThreeC, EffectiveDiffusivityModel) +{ private : + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + public: + typedef DiffusivityMillingtonQuirk type; }; -} +//! Set the vtk output fields specific to the ThreeP model +SET_TYPE_PROP(ThreePThreeC, VtkOutputFields, ThreePThreeCVtkOutputFields); + +//! Use mole fractions in the balance equations by default +SET_BOOL_PROP(ThreePThreeC, UseMoles, true); + +//! Somerton is used as default model to compute the effective thermal heat conductivity +SET_PROP(ThreePThreeCNI, ThermalConductivityModel) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; +public: + typedef ThermalConductivitySomerton type; +}; + +////////////////////////////////////////////////////////////////// +// Property values for isothermal model required for the general non-isothermal model +////////////////////////////////////////////////////////////////// + +//set isothermal VolumeVariables +SET_TYPE_PROP(ThreePThreeCNI, IsothermalVolumeVariables, ThreePThreeCVolumeVariables); + +//set isothermal LocalResidual +SET_TYPE_PROP(ThreePThreeCNI, IsothermalLocalResidual, ThreePThreeCLocalResidual); + +//set isothermal Indices +SET_TYPE_PROP(ThreePThreeCNI, IsothermalIndices, ThreePThreeCIndices); + +//set isothermal NumEq +SET_INT_PROP(ThreePThreeCNI, IsothermalNumEq, 3); + +//! Set the vtk output fields specific to the ThreeP model +SET_TYPE_PROP(ThreePThreeCNI, IsothermalVtkOutputFields, ThreePThreeCVtkOutputFields); + +} // end namespace Properties -#include "propertydefaults.hh" +} // end namespace Dumux #endif diff --git a/dumux/porousmediumflow/3p3c/implicit/newtoncontroller.hh b/dumux/porousmediumflow/3p3c/implicit/newtoncontroller.hh deleted file mode 100644 index 46483749a6d9358d4e2a65095a9ac16fd685b335..0000000000000000000000000000000000000000 --- a/dumux/porousmediumflow/3p3c/implicit/newtoncontroller.hh +++ /dev/null @@ -1,71 +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 . * - *****************************************************************************/ -/*! - * \file - * \brief A three-phase three-component specific controller for the newton solver. - * - * This controller 'knows' what a 'physically meaningful' solution is - * which allows the newton method to abort quicker if the solution is - * way out of bounds. - */ -#ifndef DUMUX_3P3C_NEWTON_CONTROLLER_HH -#define DUMUX_3P3C_NEWTON_CONTROLLER_HH - -#include "properties.hh" - -#include - -namespace Dumux { -/*! - * \ingroup Newton - * \ingroup ThreePThreeCModel - * \brief A three-phase three-component specific controller for the newton solver. - * - * This controller 'knows' what a 'physically meaningful' solution is - * which allows the newton method to abort quicker if the solution is - * way out of bounds. - */ -template -class ThreePThreeCNewtonController : public NewtonController -{ - typedef NewtonController ParentType; - typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector; - -public: - ThreePThreeCNewtonController(const Problem &problem) - : ParentType(problem) {}; - - - /*! - * \brief Returns true if the current solution can be considered to - * be accurate enough - */ - bool newtonConverged() - { - if (this->method().model().switched()) - return false; - - return ParentType::newtonConverged(); - }; -}; - -} // end namespace Dumux - -#endif diff --git a/dumux/porousmediumflow/3p3c/implicit/primaryvariableswitch.hh b/dumux/porousmediumflow/3p3c/implicit/primaryvariableswitch.hh index 9df5ff90eaf40eb7bc538d987f552476215ad738..cde9052b608674505f6abd4c51c82d8e1fad3fca 100644 --- a/dumux/porousmediumflow/3p3c/implicit/primaryvariableswitch.hh +++ b/dumux/porousmediumflow/3p3c/implicit/primaryvariableswitch.hh @@ -67,6 +67,9 @@ class ThreePThreeCPrimaryVariableSwitch : public Dumux::PrimaryVariableSwitch. * - *****************************************************************************/ -/*! - * \ingroup Properties - * \ingroup ImplicitProperties - * \ingroup ThreePThreeCModel - */ -/*! - * \file - * - * \brief Defines the properties required for the three-phase three-component - * fully implicit model. - */ -#ifndef DUMUX_3P3C_PROPERTIES_HH -#define DUMUX_3P3C_PROPERTIES_HH - -#include -#include -#include - -namespace Dumux -{ - -namespace Properties -{ -//! The type tags for the implicit three-phase three-component problems -NEW_TYPE_TAG(ThreePThreeC, INHERITS_FROM(PorousMediumFlow)); - -//! The type tags for the corresponding non-isothermal problems -NEW_TYPE_TAG(ThreePThreeCNI, INHERITS_FROM(ThreePThreeC, NonIsothermal)); - -} // end namespace Properties -} // end namespace Dumux - -#endif diff --git a/dumux/porousmediumflow/3p3c/implicit/propertydefaults.hh b/dumux/porousmediumflow/3p3c/implicit/propertydefaults.hh deleted file mode 100644 index 6c96495c46dd58dbcf8e134a996a5ba271f0b6bb..0000000000000000000000000000000000000000 --- a/dumux/porousmediumflow/3p3c/implicit/propertydefaults.hh +++ /dev/null @@ -1,214 +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 Properties - * \ingroup ImplicitProperties - * \ingroup ThreePThreeCModel - * \file - * - * \brief Defines default values for most properties required by the - * three-phase three-component fully implicit model. - */ -#ifndef DUMUX_3P3C_PROPERTY_DEFAULTS_HH -#define DUMUX_3P3C_PROPERTY_DEFAULTS_HH - -#include "indices.hh" - -#include "model.hh" -#include "volumevariables.hh" -#include "properties.hh" -#include "newtoncontroller.hh" -#include "primaryvariableswitch.hh" - -#include -#include -#include -#include -#include -#include -#include - -namespace Dumux -{ - -namespace Properties { -////////////////////////////////////////////////////////////////// -// Property values -////////////////////////////////////////////////////////////////// - -/*! - * \brief Set the property for the number of components. - * - * We just forward the number from the fluid system and use an static - * assert to make sure it is 3. - */ -SET_PROP(ThreePThreeC, 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 3p3c 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(ThreePThreeC, 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 3p3c model!"); -}; - -//! Set as default that no component mass balance is replaced by the total mass balance -SET_INT_PROP(ThreePThreeC, ReplaceCompEqIdx, 100); -/*! - * \brief The fluid state which is used by the volume variables to - * store the thermodynamic state. This should be chosen - * appropriately for the model ((non-)isothermal, equilibrium, ...). - * This can be done in the problem. - */ -SET_PROP(ThreePThreeC, FluidState){ - private: - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; - public: - typedef CompositionalFluidState type; -}; - -SET_INT_PROP(ThreePThreeC, NumEq, 3); //!< set the number of equations to 2 - -//! The local residual function of the conservation equations -SET_TYPE_PROP(ThreePThreeC, LocalResidual, CompositionalLocalResidual); - -//! Enable advection -SET_BOOL_PROP(ThreePThreeC, EnableAdvection, true); - -//! Enable molecular diffusion -SET_BOOL_PROP(ThreePThreeC, EnableMolecularDiffusion, true); - -//! Isothermal model by default -SET_BOOL_PROP(ThreePThreeC, EnableEnergyBalance, false); - -//! Use the 3p3c specific newton controller for the 3p3c model -SET_TYPE_PROP(ThreePThreeC, NewtonController, ThreePThreeCNewtonController); - -//! the Model property -SET_TYPE_PROP(ThreePThreeC, Model, ThreePThreeCModel); - -//! The primary variable switch for the 3p3c model -SET_TYPE_PROP(ThreePThreeC, PrimaryVariableSwitch, ThreePThreeCPrimaryVariableSwitch); - -//! The primary variables vector for the 2p2c model -SET_TYPE_PROP(ThreePThreeC, PrimaryVariables, SwitchablePrimaryVariables); - -//! the VolumeVariables property -SET_TYPE_PROP(ThreePThreeC, VolumeVariables, ThreePThreeCVolumeVariables); - -//! Determines whether a constraint solver should be used explicitly -SET_BOOL_PROP(ThreePThreeC, UseConstraintSolver, false); - -//! The indices required by the isothermal 3p3c model -SET_TYPE_PROP(ThreePThreeC, Indices, ThreePThreeCIndices); - -//! The spatial parameters to be employed. -//! Use ImplicitSpatialParams by default. -SET_TYPE_PROP(ThreePThreeC, SpatialParams, ImplicitSpatialParams); - -//! The model after Millington (1961) is used for the effective diffusivity -SET_PROP(ThreePThreeC, EffectiveDiffusivityModel) -{ private : - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - public: - typedef DiffusivityMillingtonQuirk type; -}; - -// disable velocity output by default - -//! 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 -// porous medium. Taking it as a constant is only a first approximation -// (Nield, Bejan, Convection in porous media, 2006, p. 10) -SET_SCALAR_PROP(BoxModel, SpatialParamsForchCoeff, 0.55); - -/*! - * \brief default value for tortuosity value (tau) used in macroscopic diffusion - * - * Value is 0.5 according to Carman 1937: Fluid flow through granular beds - * \cite carman1937 - */ -SET_SCALAR_PROP(ThreePThreeC, TauTortuosity, 0.5); - -//! Use mole fractions in the balance equations by default -SET_BOOL_PROP(ThreePThreeC, UseMoles, true); - -//! Somerton is used as default model to compute the effective thermal heat conductivity -SET_PROP(ThreePThreeCNI, 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(ThreePThreeCNI, NiOutputLevel, 0); - -////////////////////////////////////////////////////////////////// -// Property values for isothermal model required for the general non-isothermal model -////////////////////////////////////////////////////////////////// - -// set isothermal Model -SET_TYPE_PROP(ThreePThreeCNI, IsothermalModel, ThreePThreeCModel); - -//set isothermal VolumeVariables -SET_TYPE_PROP(ThreePThreeCNI, IsothermalVolumeVariables, ThreePThreeCVolumeVariables); - -//set isothermal LocalResidual -SET_TYPE_PROP(ThreePThreeCNI, IsothermalLocalResidual, CompositionalLocalResidual); - -//set isothermal Indices -SET_PROP(ThreePThreeCNI, IsothermalIndices) -{ - -public: - typedef ThreePThreeCIndices type; -}; - -//set isothermal NumEq -SET_INT_PROP(ThreePThreeCNI, IsothermalNumEq, 3); - -} - -} - -#endif diff --git a/dumux/porousmediumflow/3p3c/implicit/volumevariables.hh b/dumux/porousmediumflow/3p3c/implicit/volumevariables.hh index 7b158db51db3f8c804b7ea26456bea6fab5696e2..8850570fce4a86caa2fede675a4b93567f1278c3 100644 --- a/dumux/porousmediumflow/3p3c/implicit/volumevariables.hh +++ b/dumux/porousmediumflow/3p3c/implicit/volumevariables.hh @@ -25,12 +25,13 @@ #ifndef DUMUX_3P3C_VOLUME_VARIABLES_HH #define DUMUX_3P3C_VOLUME_VARIABLES_HH -#include +#include #include #include #include #include -#include "properties.hh" +#include +#include namespace Dumux { @@ -45,18 +46,17 @@ template class ThreePThreeCVolumeVariables : public ImplicitVolumeVariables { using ParentType = ImplicitVolumeVariables; + using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables); using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - using Model = typename GET_PROP_TYPE(TypeTag, Model); using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); - using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector); + using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector); using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw); - using MaterialLawParams = typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params; - using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector); + using MaterialLawParams = typename MaterialLaw::Params; // constraint solvers using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams); @@ -98,9 +98,6 @@ class ThreePThreeCVolumeVariables : public ImplicitVolumeVariables // universial gas constant static constexpr Scalar R = Dumux::Constants::R; - enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; - enum { dofCodim = isBox ? dim : 0 }; - public: using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState); @@ -108,7 +105,7 @@ public: /*! * \copydoc ImplicitVolumeVariables::update */ - void update(const ElementSolution &elemSol, + void update(const ElementSolutionVector &elemSol, const Problem &problem, const Element &element, const SubControlVolume& scv) diff --git a/dumux/porousmediumflow/3p3c/implicit/vtkoutputfields.hh b/dumux/porousmediumflow/3p3c/implicit/vtkoutputfields.hh new file mode 100644 index 0000000000000000000000000000000000000000..4c64a04e806752e8e42c6088216684e9668a85da --- /dev/null +++ b/dumux/porousmediumflow/3p3c/implicit/vtkoutputfields.hh @@ -0,0 +1,73 @@ +// -*- 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_THREEPTHREEC_VTK_OUTPUT_FIELDS_HH +#define DUMUX_THREEPTHREEC_VTK_OUTPUT_FIELDS_HH + +#include + +namespace Dumux +{ + +/*! + * \ingroup ThreePThreeC, InputOutput + * \brief Adds vtk output fields specific to the three-phase three-component model + */ +template +class ThreePThreeCVtkOutputFields +{ + using Indices = typename GET_PROP_TYPE(TypeTag, Indices); + using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); + using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); + + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases), + numComponents = GET_PROP_VALUE(TypeTag, NumComponents) + }; +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"); + + for (int i = 0; i < numPhases; ++i) + for (int j = 0; j < numComponents; ++j) + vtk.addVolumeVariable([i,j](const VolumeVariables& v){ return v.moleFraction(i,j); }, + "x^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(i)); + + vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.porosity(); },"porosity"); + vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.priVars().state(); },"phase presence"); + vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.permeability(); },"permeability"); + } +}; + +} // end namespace Dumux + +#endif diff --git a/test/porousmediumflow/3p3c/implicit/CMakeLists.txt b/test/porousmediumflow/3p3c/implicit/CMakeLists.txt index 1bc428487005c7fc8c635c4714b9dc19c26f97c0..d9a03cda4fcb9558b0bd73db2e79f2aa8df26cf9 100644 --- a/test/porousmediumflow/3p3c/implicit/CMakeLists.txt +++ b/test/porousmediumflow/3p3c/implicit/CMakeLists.txt @@ -1,48 +1,60 @@ add_input_file_links() # isothermal tests -add_dumux_test(test_box3p3c test_box3p3c test_box3p3c.cc - python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/infiltrationbox-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/box3p3c-00005.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_box3p3c -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/test_box3p3c_reference.input") +dune_add_test(NAME test_3p3c_box + SOURCES test_3p3c_fv.cc + COMPILE_DEFINITIONS TYPETAG=InfiltrationThreePThreeCBoxTypeTag + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/infiltrationbox-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/infiltration3p3cbox-00005.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_3p3c_box test_3p3c_fv.input -Problem.Name infiltration3p3cbox") -add_dumux_test(test_cc3p3c test_cc3p3c test_cc3p3c.cc - python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/infiltrationcc-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/cc3p3c-00005.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_cc3p3c -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/test_cc3p3c_reference.input") +dune_add_test(NAME test_3p3c_tpfa + SOURCES test_3p3c_fv.cc + COMPILE_DEFINITIONS TYPETAG=InfiltrationThreePThreeCCCTpfaTypeTag + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/infiltrationcc-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/infiltration3p3ccc-00005.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_3p3c_tpfa test_3p3c_fv.input -Problem.Name infiltration3p3ccc") # non-isothermal tests -add_dumux_test(test_box3p3cnicolumnxylol test_box3p3cnicolumnxylol test_box3p3cnicolumnxylol.cc - python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/column3p3cnibox-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/columnxylolbox-00062.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_box3p3cnicolumnxylol") +dune_add_test(NAME test_3p3cni_columnxylol_box + SOURCES test_3p3c_fv.cc + COMPILE_DEFINITIONS TYPETAG=ColumnBoxTypeTag + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/column3p3cnibox-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/columnxylol_box-00062.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_3p3cni_columnxylol_box test_columnxylol_fv.input -Problem.Name columnxylol_box") -add_dumux_test(test_cc3p3cnicolumnxylol test_cc3p3cnicolumnxylol test_cc3p3cnicolumnxylol.cc - python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/column3p3cnicc-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/columnxylolcc-00053.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_cc3p3cnicolumnxylol") +dune_add_test(NAME test_3p3cni_columnxylol_tpfa + SOURCES test_3p3c_fv.cc + COMPILE_DEFINITIONS TYPETAG=ColumnCCTpfaTypeTag + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/column3p3cnicc-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/columnxylol_tpfa-00053.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_3p3cni_columnxylol_tpfa test_columnxylol_fv.input -Problem.Name columnxylol_tpfa") -add_dumux_test(test_box3p3cnikuevette test_box3p3cnikuevette test_box3p3cnikuevette.cc - python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/kuevette3p3cnibox-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/kuevettebox-00004.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_box3p3cnikuevette") +dune_add_test(NAME test_3p3cni_kuevette_box + SOURCES test_3p3c_fv.cc + COMPILE_DEFINITIONS TYPETAG=KuevetteBoxTypeTag + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/kuevette3p3cnibox-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/kuevette_box-00004.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_3p3cni_kuevette_box test_kuvette_fv.input -Problem.Name kuevette_box") -add_dumux_test(test_cc3p3cnikuevette test_cc3p3cnikuevette test_cc3p3cnikuevette.cc - python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py - --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/kuevette3p3cnicc-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/kuevettecc-00004.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_cc3p3cnikuevette") +dune_add_test(NAME test_3p3cni_kuevette_tpfa + SOURCES test_3p3c_fv.cc + COMPILE_DEFINITIONS TYPETAG=KuevetteCCTpfaTypeTag + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/kuevette3p3cnicc-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/kuvette_tpfa-00004.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_3p3cni_kuevette_tpfa test_kuvette_fv.input -Problem.Name kuvette_tpfa") #install sources install(FILES @@ -52,10 +64,8 @@ infiltrationproblem.hh infiltrationspatialparameters.hh kuevetteproblem.hh kuevettespatialparams.hh -test_box3p3c.cc -test_box3p3cnicolumnxylol.cc -test_box3p3cnikuevette.cc -test_cc3p3c.cc -test_cc3p3cnicolumnxylol.cc -test_cc3p3cnikuevette.cc +test_3p3c_fv.cc +test_3p3c_fv.input +test_columnxylol_fv.input +test_kuvette_fv.input DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/implicit/3p3c) diff --git a/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh b/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh index c9df17677848c5e44a1fbef78a54ce2d082076dd..04d8d72321145acda774dde58e1e78f276cc0e72 100644 --- a/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh +++ b/test/porousmediumflow/3p3c/implicit/columnxylolproblem.hh @@ -26,9 +26,10 @@ #define DUMUX_COLUMNXYLOLPROBLEM_HH #include -#include +#include +#include #include -#include +#include #include "columnxylolspatialparams.hh" @@ -41,18 +42,18 @@ class ColumnProblem; namespace Properties { -NEW_TYPE_TAG(ColumnProblem, INHERITS_FROM(ThreePThreeCNI, ColumnSpatialParams)); -NEW_TYPE_TAG(ColumnBoxProblem, INHERITS_FROM(BoxModel, ColumnProblem)); -NEW_TYPE_TAG(ColumnCCProblem, INHERITS_FROM(CCTpfaModel, ColumnProblem)); +NEW_TYPE_TAG(ColumnTypeTag, INHERITS_FROM(ThreePThreeCNI, ColumnSpatialParams)); +NEW_TYPE_TAG(ColumnBoxTypeTag, INHERITS_FROM(BoxModel, ColumnTypeTag)); +NEW_TYPE_TAG(ColumnCCTpfaTypeTag, INHERITS_FROM(CCTpfaModel, ColumnTypeTag)); // Set the grid type -SET_TYPE_PROP(ColumnProblem, Grid, Dune::YaspGrid<2>); +SET_TYPE_PROP(ColumnTypeTag, Grid, Dune::YaspGrid<2>); // Set the problem property -SET_TYPE_PROP(ColumnProblem, Problem, ColumnProblem); +SET_TYPE_PROP(ColumnTypeTag, Problem, ColumnProblem); // Set the fluid system -SET_TYPE_PROP(ColumnProblem, +SET_TYPE_PROP(ColumnTypeTag, FluidSystem, FluidSystems::H2OAirXylene); } @@ -60,7 +61,6 @@ SET_TYPE_PROP(ColumnProblem, /*! * \ingroup ThreePThreeCModel - * \ingroup ImplicitTestProblems * \brief Non-isothermal injection problem where a water is injected into a * sand column with a NAPL contamination. * @@ -81,20 +81,16 @@ SET_TYPE_PROP(ColumnProblem, * To adjust the simulation time it is necessary to edit the input file. * * To run the simulation execute the following line in shell: - * ./test_box3p3cnicolumnxylol test_box3p3cnicolumnxylol.input or - * ./test_cc3p3cnicolumnxylol test_cc3p3cnicolumnxylol.input + * ./test_box3p3cnicolumnxylol test_columnxylol_fv.input or + * ./test_cc3p3cnicolumnxylol test_columnxylol_fv.input */ template -class ColumnProblem : public ImplicitPorousMediaProblem +class ColumnProblem : public PorousMediumFlowProblem { - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GridView::Grid Grid; - - typedef ImplicitPorousMediaProblem ParentType; - - // copy some indices for convenience - typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using ParentType = PorousMediumFlowProblem; + using Indices = typename GET_PROP_TYPE(TypeTag, Indices); enum { pressureIdx = Indices::pressureIdx, @@ -111,26 +107,16 @@ class ColumnProblem : public ImplicitPorousMediaProblem dimWorld = GridView::dimensionworld }; - - typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; - typedef typename GET_PROP_TYPE(TypeTag, NumEqVector) NeumannFluxes; - 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, FluidSystem) FluidSystem; + 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 Element = typename GridView::template Codim<0>::Entity; + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); - using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); - - typedef Dune::FieldVector GlobalPosition; - - enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; - enum { dofCodim = isBox ? dim : 0 }; + using GlobalPosition = Dune::FieldVector; public: /*! @@ -139,12 +125,11 @@ public: * \param timeManager The time manager * \param gridView The grid view */ - ColumnProblem(TimeManager &timeManager, const GridView &gridView) - : ParentType(timeManager, gridView) + ColumnProblem(std::shared_ptr fvGridGeometry) + : ParentType(fvGridGeometry) { FluidSystem::init(); - - name_ = GET_RUNTIME_PARAM(TypeTag, std::string, Problem.Name); + name_ = getParam("Problem.Name"); } /*! @@ -160,18 +145,6 @@ public: const std::string name() const { return name_; } - /*! - * \brief Returns the source term - * - * \param values Stores the source values for the conservation equations in - * \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$ - * \param globalPos The global position - */ - PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const - { - return PrimaryVariables(0.0); - } - // \} /*! @@ -208,7 +181,6 @@ public: PrimaryVariables dirichletAtPos( const GlobalPosition &globalPos) const { return initial_(globalPos); - } /*! @@ -225,15 +197,15 @@ public: * in normal direction of each phase. Negative values mean influx. */ NeumannFluxes neumann(const Element& element, - const FVElementGeometry& fvGeometry, - const ElementVolumeVariables& elemVolVars, - const SubControlVolumeFace& scvf) const + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf) const { PrimaryVariables values(0.0); const auto& globalPos = scvf.ipGlobal(); // negative values for injection - if (globalPos[1] > 1.2 - eps_) + if (globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_) { values[Indices::contiWEqIdx] = -0.395710; values[Indices::contiGEqIdx] = -0.000001; @@ -261,34 +233,26 @@ public: } - /*! + /*! * \brief Append all quantities of interest which can be derived * from the solution of the current time step to the VTK * writer. Adjust this in case of anisotropic permeabilities. */ - void addOutputVtkFields() + template + void addVtkFields(VTKWriter& vtk) { - // get the number of degrees of freedom - unsigned numDofs = this->model().numDofs(); - - // create the scalar field required for the permeabilities - typedef Dune::BlockVector > ScalarField; - ScalarField *Kxx = this->resultWriter().allocateManagedBuffer(numDofs); - - FVElementGeometry fvGeometry; + const auto& gg = this->fvGridGeometry(); + Kxx_.resize(gg.numDofs()); + vtk.addField(Kxx_, "permeability"); for (const auto& element : elements(this->gridView())) { - fvGeometry.update(this->gridView(), element); + auto fvGeometry = localView(gg); + fvGeometry.bindElement(element); - for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx) - { - int dofIdxGlobal = this->model().dofMapper().subIndex(element, scvIdx, dofCodim); - (*Kxx)[dofIdxGlobal] = this->spatialParams().intrinsicPermeability(element, fvGeometry, scvIdx); - } + for (const auto& scv : scvs(fvGeometry)) + Kxx_[scv.dofIndex()] = this->spatialParams().intrinsicPermeabilityAtPos(scv.dofPosition()); } - - this->resultWriter().attachDofData(*Kxx, "permeability", isBox); } private: @@ -298,53 +262,54 @@ private: { PrimaryVariables values; values.setState(threePhases); - Scalar y = globalPos[1]; + const auto y = globalPos[1]; + const auto yMax = this->fvGridGeometry().bBoxMax()[1]; values[temperatureIdx] = 296.15; values[pressureIdx] = 1.e5; values[switch1Idx] = 0.005; - if (y > 1.2 - eps_) + if (y > yMax - eps_) values[switch2Idx] = 0.112; - else if (y > 1.2 - 0.0148 - eps_) - values[switch2Idx] = 0 + ((1.2 - y)/0.0148)*0.112; - else if (y > 1.2 - 0.0296 - eps_) - values[switch2Idx] = 0.112 + (((1.2 - y) - 0.0148)/0.0148)*(0.120 - 0.112); - else if (y > 1.2 - 0.0444 - eps_) - values[switch2Idx] = 0.120 + (((1.2 - y) - 0.0296)/0.0148)*(0.125 - 0.120); - else if (y > 1.2 - 0.0592 - eps_) - values[switch2Idx] = 0.125 + (((1.2 - y) - 0.0444)/0.0148)*(0.137 - 0.125); - else if (y > 1.2 - 0.0740 - eps_) - values[switch2Idx] = 0.137 + (((1.2 - y) - 0.0592)/0.0148)*(0.150 - 0.137); - else if (y > 1.2 - 0.0888 - eps_) - values[switch2Idx] = 0.150 + (((1.2 - y) - 0.0740)/0.0148)*(0.165 - 0.150); - else if (y > 1.2 - 0.1036 - eps_) - values[switch2Idx] = 0.165 + (((1.2 - y) - 0.0888)/0.0148)*(0.182 - 0.165); - else if (y > 1.2 - 0.1184 - eps_) - values[switch2Idx] = 0.182 + (((1.2 - y) - 0.1036)/0.0148)*(0.202 - 0.182); - else if (y > 1.2 - 0.1332 - eps_) - values[switch2Idx] = 0.202 + (((1.2 - y) - 0.1184)/0.0148)*(0.226 - 0.202); - else if (y > 1.2 - 0.1480 - eps_) - values[switch2Idx] = 0.226 + (((1.2 - y) - 0.1332)/0.0148)*(0.257 - 0.226); - else if (y > 1.2 - 0.1628 - eps_) - values[switch2Idx] = 0.257 + (((1.2 - y) - 0.1480)/0.0148)*(0.297 - 0.257); - else if (y > 1.2 - 0.1776 - eps_) - values[switch2Idx] = 0.297 + (((1.2 - y) - 0.1628)/0.0148)*(0.352 - 0.297); - else if (y > 1.2 - 0.1924 - eps_) - values[switch2Idx] = 0.352 + (((1.2 - y) - 0.1776)/0.0148)*(0.426 - 0.352); - else if (y > 1.2 - 0.2072 - eps_) - values[switch2Idx] = 0.426 + (((1.2 - y) - 0.1924)/0.0148)*(0.522 - 0.426); - else if (y > 1.2 - 0.2220 - eps_) - values[switch2Idx] = 0.522 + (((1.2 - y) - 0.2072)/0.0148)*(0.640 - 0.522); - else if (y > 1.2 - 0.2368 - eps_) - values[switch2Idx] = 0.640 + (((1.2 - y) - 0.2220)/0.0148)*(0.767 - 0.640); - else if (y > 1.2 - 0.2516 - eps_) - values[switch2Idx] = 0.767 + (((1.2 - y) - 0.2368)/0.0148)*(0.878 - 0.767); - else if (y > 1.2 - 0.2664 - eps_) - values[switch2Idx] = 0.878 + (((1.2 - y) - 0.2516)/0.0148)*(0.953 - 0.878); - else if (y > 1.2 - 0.2812 - eps_) - values[switch2Idx] = 0.953 + (((1.2 - y) - 0.2664)/0.0148)*(0.988 - 0.953); - else if (y > 1.2 - 0.3000 - eps_) + else if (y > yMax - 0.0148 - eps_) + values[switch2Idx] = 0 + ((yMax - y)/0.0148)*0.112; + else if (y > yMax - 0.0296 - eps_) + values[switch2Idx] = 0.112 + (((yMax - y) - 0.0148)/0.0148)*(0.120 - 0.112); + else if (y > yMax - 0.0444 - eps_) + values[switch2Idx] = 0.120 + (((yMax - y) - 0.0296)/0.0148)*(0.125 - 0.120); + else if (y > yMax - 0.0592 - eps_) + values[switch2Idx] = 0.125 + (((yMax - y) - 0.0444)/0.0148)*(0.137 - 0.125); + else if (y > yMax - 0.0740 - eps_) + values[switch2Idx] = 0.137 + (((yMax - y) - 0.0592)/0.0148)*(0.150 - 0.137); + else if (y > yMax - 0.0888 - eps_) + values[switch2Idx] = 0.150 + (((yMax - y) - 0.0740)/0.0148)*(0.165 - 0.150); + else if (y > yMax - 0.1036 - eps_) + values[switch2Idx] = 0.165 + (((yMax - y) - 0.0888)/0.0148)*(0.182 - 0.165); + else if (y > yMax - 0.1184 - eps_) + values[switch2Idx] = 0.182 + (((yMax - y) - 0.1036)/0.0148)*(0.202 - 0.182); + else if (y > yMax - 0.1332 - eps_) + values[switch2Idx] = 0.202 + (((yMax - y) - 0.1184)/0.0148)*(0.226 - 0.202); + else if (y > yMax - 0.1480 - eps_) + values[switch2Idx] = 0.226 + (((yMax - y) - 0.1332)/0.0148)*(0.257 - 0.226); + else if (y > yMax - 0.1628 - eps_) + values[switch2Idx] = 0.257 + (((yMax - y) - 0.1480)/0.0148)*(0.297 - 0.257); + else if (y > yMax - 0.1776 - eps_) + values[switch2Idx] = 0.297 + (((yMax - y) - 0.1628)/0.0148)*(0.352 - 0.297); + else if (y > yMax - 0.1924 - eps_) + values[switch2Idx] = 0.352 + (((yMax - y) - 0.1776)/0.0148)*(0.426 - 0.352); + else if (y > yMax - 0.2072 - eps_) + values[switch2Idx] = 0.426 + (((yMax - y) - 0.1924)/0.0148)*(0.522 - 0.426); + else if (y > yMax - 0.2220 - eps_) + values[switch2Idx] = 0.522 + (((yMax - y) - 0.2072)/0.0148)*(0.640 - 0.522); + else if (y > yMax - 0.2368 - eps_) + values[switch2Idx] = 0.640 + (((yMax - y) - 0.2220)/0.0148)*(0.767 - 0.640); + else if (y > yMax - 0.2516 - eps_) + values[switch2Idx] = 0.767 + (((yMax - y) - 0.2368)/0.0148)*(0.878 - 0.767); + else if (y > yMax - 0.2664 - eps_) + values[switch2Idx] = 0.878 + (((yMax - y) - 0.2516)/0.0148)*(0.953 - 0.878); + else if (y > yMax - 0.2812 - eps_) + values[switch2Idx] = 0.953 + (((yMax - y) - 0.2664)/0.0148)*(0.988 - 0.953); + else if (y > yMax - 0.3000 - eps_) values[switch2Idx] = 0.988; else values[switch2Idx] = 1.e-4; @@ -353,7 +318,9 @@ private: static constexpr Scalar eps_ = 1e-6; std::string name_; + std::vector Kxx_; }; -} //end namespace + +} //end namespace Dumux #endif diff --git a/test/porousmediumflow/3p3c/implicit/columnxylolspatialparams.hh b/test/porousmediumflow/3p3c/implicit/columnxylolspatialparams.hh index d494c61fe812610566078d779b8f5c9bd4873ed2..dec0e465dae86df01015043d7ac982d1941a5321 100644 --- a/test/porousmediumflow/3p3c/implicit/columnxylolspatialparams.hh +++ b/test/porousmediumflow/3p3c/implicit/columnxylolspatialparams.hh @@ -24,7 +24,6 @@ #ifndef DUMUX_COLUMNXYLOL_SPATIAL_PARAMS_HH #define DUMUX_COLUMNXYLOL_SPATIAL_PARAMS_HH -#include #include #include #include @@ -51,57 +50,36 @@ SET_PROP(ColumnSpatialParams, 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); public: // define the material law parameterized by absolute saturations - typedef EffToAbsLaw type; + using type = EffToAbsLaw>; }; -} +} // end namespace Properties /*! * \ingroup ThreePThreeCModel - * \ingroup ImplicitTestProblems * \brief Definition of the spatial parameters for the column problem */ template class ColumnSpatialParams : public ImplicitSpatialParams { - typedef ImplicitSpatialParams ParentType; - - typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; - 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 GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); enum { dimWorld=GridView::dimensionworld, dim=GridView::dimension }; - - typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; - enum { - wPhaseIdx = Indices::wPhaseIdx, - nPhaseIdx = Indices::nPhaseIdx - }; - - typedef Dune::FieldVector GlobalPosition; - typedef Dune::FieldVector DimVector; - - - typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector; - - typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables; - typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables; - - typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; - typedef typename GridView::template Codim<0>::Entity Element; + using Element = typename GridView::template Codim<0>::Entity; using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector); using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GlobalPosition = Dune::FieldVector; public: - 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; using PermeabilityType = Scalar; /*! @@ -109,8 +87,8 @@ public: * * \param gridView The grid view */ - ColumnSpatialParams(const Problem& problem, const GridView &gridView) - : ParentType(problem, gridView) + ColumnSpatialParams(const Problem& problem) + : ParentType(problem) { // intrinsic permeabilities fineK_ = 1.4e-11; @@ -151,28 +129,20 @@ public: fineMaterialParams_.setRhoBulk(1500.); } - ~ColumnSpatialParams() - {} - - /*! - * \brief Update the spatial parameters with the flow solution - * after a timestep. + * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$ + * \note It is possibly solution dependent. * - * \param globalSolution The global solution vector + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \return permeability */ - void update(const SolutionVector &globalSolution) - { - } - - /*! - * \brief Apply the intrinsic permeability tensor to a pressure - * potential gradient. - * - * \param globalPos The global position - */ - Scalar permeabilityAtPos(const GlobalPosition& globalPos) const + PermeabilityType permeability(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const { + const auto& globalPos = scv.dofPosition(); if (isFineMaterial_(globalPos)) return fineK_; return coarseK_; @@ -181,14 +151,15 @@ public: /*! * \brief Define the porosity \f$[-]\f$ of the spatial parameters * - * \param element The finite element - * \param fvGeometry The finite volume geometry - * \param scvIdx The local index of the sub-control volume where - * the porosity needs to be defined + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. */ - Scalar porosityAtPos(const GlobalPosition& globalPos) const + Scalar porosity(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const { - + const auto& globalPos = scv.dofPosition(); if (isFineMaterial_(globalPos)) return finePorosity_; else @@ -197,13 +168,18 @@ public: /*! - * \brief return the parameter object for the Brooks-Corey material law which depends on the position + * \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.). * - * \param globalPos The global position + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \return the material parameters object */ - const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const + const MaterialLawParams& materialLawParams(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const { - + const auto& globalPos = scv.dofPosition(); if (isFineMaterial_(globalPos)) return fineMaterialParams_; else @@ -215,12 +191,15 @@ public: * * This is only required for non-isothermal models. * - * \param element The finite element - * \param fvGeometry The finite volume geometry - * \param scvIdx The local index of the sub-control volume + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. */ - Scalar solidHeatCapacityAtPos(const GlobalPosition& globalPos) const + Scalar solidHeatCapacity(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const { + const auto& globalPos = scv.dofPosition(); if (isFineMaterial_(globalPos)) return fineHeatCap_; else @@ -261,9 +240,7 @@ public: private: bool isFineMaterial_(const GlobalPosition &globalPos) const { - if (0.90 - eps_ <= globalPos[1]) - return true; - else return false; + return (0.90 - eps_ <= globalPos[1]); } Scalar fineK_; @@ -283,6 +260,6 @@ private: static constexpr Scalar eps_ = 1e-6; }; -} +} // end namespace Dumux #endif diff --git a/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh b/test/porousmediumflow/3p3c/implicit/infiltration3p3cproblem.hh similarity index 72% rename from test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh rename to test/porousmediumflow/3p3c/implicit/infiltration3p3cproblem.hh index aa550116967433284814e5b49c7daa7e3967bcb7..c0000696090a6aea9764a3aad5a1f9dbc9015f73 100644 --- a/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh +++ b/test/porousmediumflow/3p3c/implicit/infiltration3p3cproblem.hh @@ -22,44 +22,38 @@ * \brief Isothermal NAPL infiltration problem: LNAPL contaminates * the unsaturated and the saturated groundwater zone. */ -#ifndef DUMUX_INFILTRATIONPROBLEM_HH -#define DUMUX_INFILTRATIONPROBLEM_HH - -#include +#ifndef DUMUX_INFILTRATION_THREEPTHREEC_PROBLEM_HH +#define DUMUX_INFILTRATION_THREEPTHREEC_PROBLEM_HH +#include +#include +#include #include -#include +#include -#include "infiltrationspatialparameters.hh" +#include "infiltration3p3cspatialparams.hh" namespace Dumux { template -class InfiltrationProblem; +class InfiltrationThreePThreeCProblem; namespace Properties { -NEW_TYPE_TAG(InfiltrationProblem, INHERITS_FROM(ThreePThreeC, InfiltrationSpatialParams)); -NEW_TYPE_TAG(InfiltrationBoxProblem, INHERITS_FROM(BoxModel, InfiltrationProblem)); -NEW_TYPE_TAG(InfiltrationCCProblem, INHERITS_FROM(CCTpfaModel, InfiltrationProblem)); +NEW_TYPE_TAG(InfiltrationThreePThreeCTypeTag, INHERITS_FROM(ThreePThreeC, InfiltrationThreePThreeCSpatialParamsTypeTag)); +NEW_TYPE_TAG(InfiltrationThreePThreeCBoxTypeTag, INHERITS_FROM(BoxModel, InfiltrationThreePThreeCTypeTag)); +NEW_TYPE_TAG(InfiltrationThreePThreeCCCTpfaTypeTag, INHERITS_FROM(CCTpfaModel, InfiltrationThreePThreeCTypeTag)); // Set the grid type -SET_TYPE_PROP(InfiltrationProblem, Grid, Dune::YaspGrid<2>); +SET_TYPE_PROP(InfiltrationThreePThreeCTypeTag, Grid, Dune::YaspGrid<2>); // Set the problem property -SET_TYPE_PROP(InfiltrationProblem, Problem, InfiltrationProblem); +SET_TYPE_PROP(InfiltrationThreePThreeCTypeTag, Problem, InfiltrationThreePThreeCProblem); // Set the fluid system -SET_TYPE_PROP(InfiltrationProblem, +SET_TYPE_PROP(InfiltrationThreePThreeCTypeTag, FluidSystem, FluidSystems::H2OAirMesitylene); - - -// Maximum tolerated relative error in the Newton method -SET_SCALAR_PROP(InfiltrationProblem, NewtonMaxRelativeShift, 1e-8); - -// -1 backward differences, 0: central differences, +1: forward differences -SET_INT_PROP(InfiltrationProblem, ImplicitNumericDifferenceMethod, 0); } /*! @@ -93,21 +87,17 @@ SET_INT_PROP(InfiltrationProblem, ImplicitNumericDifferenceMethod, 0); * To run the simulation execute the following line in shell: * ./test_box3p3c test_box3p3c.input or * ./test_cc3p3c test_cc3p3c.input - */ + * */ template -class InfiltrationProblem : public ImplicitPorousMediaProblem +class InfiltrationThreePThreeCProblem : public PorousMediumFlowProblem { - typedef ImplicitPorousMediaProblem ParentType; - - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; + using ParentType = PorousMediumFlowProblem; - typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; - typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params 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, switch1Idx = Indices::switch1Idx, @@ -125,23 +115,16 @@ class InfiltrationProblem : 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, NumEqVector) NeumannFluxes; - 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, 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) }; - enum { dofCodim = isBox ? dim : 0 }; + using GlobalPosition = Dune::FieldVector; public: /*! @@ -150,8 +133,8 @@ public: * \param timeManager The time manager * \param gridView The grid view */ - InfiltrationProblem(TimeManager &timeManager, const GridView &gridView) - : ParentType(timeManager, gridView) + InfiltrationThreePThreeCProblem(std::shared_ptr fvGridGeometry) + : ParentType(fvGridGeometry) { temperature_ = 273.15 + 10.0; // -> 10 degrees Celsius FluidSystem::init(/*tempMin=*/temperature_ - 1, @@ -161,14 +144,9 @@ public: /*pressMax=*/3*1e5, /*nPress=*/200); - name_ = GET_RUNTIME_PARAM(TypeTag, std::string, Problem.Name); + name_ = getParam("Problem.Name"); } - /*! - * \name Problem parameters - */ - // \{ - /*! * \brief The problem name. * @@ -206,9 +184,9 @@ public: BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const { BoundaryTypes values; - if(globalPos[0] > 500. - eps_) + if(globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_) values.setAllDirichlet(); - else if(globalPos[0] < eps_) + else if(globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_) values.setAllDirichlet(); else values.setAllNeumann(); @@ -225,34 +203,7 @@ public: * For this method, the \a values parameter stores primary variables. */ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const - { - PrimaryVariables values; - values.setState(wgPhaseOnly); - - 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[switch1Idx] = sw; - values[switch2Idx] = 1.e-6; - }else { - values[pressureIdx] = 1e5 + 9.81 * 1000.0 * ((-5E-4*x+5) - y); - values[switch1Idx] = 1.-sgr; - values[switch2Idx] = 1.e-6; - } - - return values; - } + { return initial_(globalPos); } /*! * \brief Evaluate the boundary conditions for a neumann @@ -273,7 +224,7 @@ public: NeumannFluxes values(0.0); // negative values for injection - if ((globalPos[0] <= 75.0 + eps_) && (globalPos[0] >= 50.0 - eps_) && (globalPos[1] >= 10.0 - eps_)) + if ((globalPos[0] < 80.0 + eps_) && (globalPos[0] > 55.0 - eps_) && (globalPos[1] > 10.0 - eps_)) { values[contiWEqIdx] = -0.0; //mole flow conversion to mass flow with molar mass M(Mesit.)=0,120 kg/mol --> 1.2e-4 kg/(sm) @@ -366,6 +317,6 @@ private: static constexpr Scalar eps_ = 1e-6; std::string name_; }; -} //end namespace +} //end namespace Dumux #endif diff --git a/test/porousmediumflow/3p3c/implicit/infiltrationspatialparameters.hh b/test/porousmediumflow/3p3c/implicit/infiltration3p3cspatialparams.hh similarity index 72% rename from test/porousmediumflow/3p3c/implicit/infiltrationspatialparameters.hh rename to test/porousmediumflow/3p3c/implicit/infiltration3p3cspatialparams.hh index 9853ffc28ffbc73ff0477b0c4260b74bc791df04..4f09df6ead43bebbf4ee4c1727f27b047b30f188 100644 --- a/test/porousmediumflow/3p3c/implicit/infiltrationspatialparameters.hh +++ b/test/porousmediumflow/3p3c/implicit/infiltration3p3cspatialparams.hh @@ -23,43 +23,42 @@ * which uses the isothermal two-phase two component * fully implicit model. */ -#ifndef DUMUX_INFILTRATION_SPATIAL_PARAMETERS_HH -#define DUMUX_INFILTRATION_SPATIAL_PARAMETERS_HH +#ifndef DUMUX_INFILTRATION_THREEPTHREEC_SPATIAL_PARAMETERS_HH +#define DUMUX_INFILTRATION_THREEPTHREEC_SPATIAL_PARAMETERS_HH -#include #include #include #include #include +#include namespace Dumux { //forward declaration template -class InfiltrationSpatialParams; +class InfiltrationThreePThreeCSpatialParams; namespace Properties { // The spatial parameters TypeTag -NEW_TYPE_TAG(InfiltrationSpatialParams); +NEW_TYPE_TAG(InfiltrationThreePThreeCSpatialParamsTypeTag); // Set the spatial parameters -SET_TYPE_PROP(InfiltrationSpatialParams, SpatialParams, InfiltrationSpatialParams); +SET_TYPE_PROP(InfiltrationThreePThreeCSpatialParamsTypeTag, SpatialParams, InfiltrationThreePThreeCSpatialParams); // Set the material Law -SET_PROP(InfiltrationSpatialParams, MaterialLaw) +SET_PROP(InfiltrationThreePThreeCSpatialParamsTypeTag, 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); public: // define the material law parameterized by absolute saturations - typedef EffToAbsLaw type; + using type = EffToAbsLaw>; }; -} +} // end namespace Properties /*! * \ingroup ThreePThreeCModel @@ -67,41 +66,25 @@ SET_PROP(InfiltrationSpatialParams, MaterialLaw) * \brief Definition of the spatial parameters for the infiltration problem */ template -class InfiltrationSpatialParams : public ImplicitSpatialParams +class InfiltrationThreePThreeCSpatialParams : public ImplicitSpatialParams { using ParentType = ImplicitSpatialParams; - using Grid = typename GET_PROP_TYPE(TypeTag, Grid); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); - using CoordScalar = typename Grid::ctype; enum { dim = GridView::dimension, dimWorld=GridView::dimensionworld }; - using Indices = typename GET_PROP_TYPE(TypeTag, Indices); - enum { - wPhaseIdx = Indices::wPhaseIdx, - nPhaseIdx = Indices::nPhaseIdx - }; - - using GlobalPosition = Dune::FieldVector; - using DimVector = Dune::FieldVector; - using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); - using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables); - using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables); - using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + using GlobalPosition = Dune::FieldVector; using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); + using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector); using Element = typename GridView::template Codim<0>::Entity; - //get the material law from the property system - using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw); - using MaterialLawParams = typename MaterialLaw::Params; - - public: - + using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw); + using MaterialLawParams = typename MaterialLaw::Params; using PermeabilityType = Scalar; /*! @@ -109,8 +92,8 @@ public: * * \param gridView The grid view */ - InfiltrationSpatialParams(const Problem& problem, const GridView &gridView) - : ParentType(problem, gridView) + InfiltrationThreePThreeCSpatialParams(const Problem& problem) + : ParentType(problem) { // intrinsic permeabilities fineK_ = 1.e-11; @@ -135,13 +118,19 @@ public: } /*! - * \brief Returns the intrinsic permeability tensor \f$[m^2]\f$ + * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$ + * \note It is possibly solution dependent. * - * \param globalPos The global position + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \return permeability */ - - PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const + PermeabilityType permeability(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const { + const auto& globalPos = scv.dofPosition(); if (isFineMaterial_(globalPos)) return fineK_; return coarseK_; @@ -185,6 +174,6 @@ private: MaterialLawParams materialParams_; }; -} +} // end namespace Dumux #endif diff --git a/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh b/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh index 1a4458334d3a6b0a524db3dbd563ba53e4ae6755..4d066889879515ad939e08689d2e3ff360f3780b 100644 --- a/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh +++ b/test/porousmediumflow/3p3c/implicit/kuevetteproblem.hh @@ -29,9 +29,10 @@ #include #include -#include +#include +#include #include -#include +#include #include "kuevettespatialparams.hh" @@ -44,23 +45,20 @@ class KuevetteProblem; namespace Properties { -NEW_TYPE_TAG(KuevetteProblem, INHERITS_FROM(ThreePThreeCNI, KuevetteSpatialParams)); -NEW_TYPE_TAG(KuevetteBoxProblem, INHERITS_FROM(BoxModel, KuevetteProblem)); -NEW_TYPE_TAG(KuevetteCCProblem, INHERITS_FROM(CCTpfaModel, KuevetteProblem)); +NEW_TYPE_TAG(KuevetteTypeTag, INHERITS_FROM(ThreePThreeCNI, KuevetteSpatialParams)); +NEW_TYPE_TAG(KuevetteBoxTypeTag, INHERITS_FROM(BoxModel, KuevetteTypeTag)); +NEW_TYPE_TAG(KuevetteCCTpfaTypeTag, INHERITS_FROM(CCTpfaModel, KuevetteTypeTag)); // Set the grid type -SET_TYPE_PROP(KuevetteProblem, Grid, Dune::YaspGrid<2>); +SET_TYPE_PROP(KuevetteTypeTag, Grid, Dune::YaspGrid<2>); // Set the problem property -SET_TYPE_PROP(KuevetteProblem, Problem, KuevetteProblem); +SET_TYPE_PROP(KuevetteTypeTag, Problem, KuevetteProblem); // Set the fluid system -SET_TYPE_PROP(KuevetteProblem, +SET_TYPE_PROP(KuevetteTypeTag, FluidSystem, FluidSystems::H2OAirMesitylene); - -// set newton relative tolerance -SET_SCALAR_PROP(KuevetteProblem, NewtonMaxRelativeShift, 1e-6); } @@ -98,16 +96,14 @@ SET_SCALAR_PROP(KuevetteProblem, NewtonMaxRelativeShift, 1e-6); * ./test_cc3p3cnikuevette test_cc3p3cnikuevette.input */ template -class KuevetteProblem : public ImplicitPorousMediaProblem +class KuevetteProblem : public PorousMediumFlowProblem { - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GridView::Grid Grid; - - typedef ImplicitPorousMediaProblem ParentType; + using ParentType = PorousMediumFlowProblem; + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); // copy some indices for convenience - typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + using Indices = typename GET_PROP_TYPE(TypeTag, Indices); enum { pressureIdx = Indices::pressureIdx, @@ -125,28 +121,18 @@ class KuevetteProblem : public ImplicitPorousMediaProblem dimWorld = GridView::dimensionworld }; - - typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; - typedef typename GET_PROP_TYPE(TypeTag, NumEqVector) NeumannFluxes; - 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, FluidSystem) FluidSystem; + 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 Element = typename GridView::template Codim<0>::Entity; + using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); + using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry); + using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem); using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables); - using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace); - typedef Dune::FieldVector GlobalPosition; - - enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) }; - enum { dofCodim = isBox ? dim : 0 }; - + using GlobalPosition = Dune::FieldVector; public: /*! * \brief The constructor. @@ -154,14 +140,11 @@ public: * \param timeManager The time manager * \param gridView The grid view */ - KuevetteProblem(TimeManager &timeManager, const GridView &gridView) - : ParentType(timeManager, gridView) + KuevetteProblem(std::shared_ptr fvGridGeometry) + : ParentType(fvGridGeometry) { FluidSystem::init(); - - name_ = GET_RUNTIME_PARAM(TypeTag, std::string, Problem.Name); - episodeLength_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, TimeManager, EpisodeLength); - this->timeManager().startNextEpisode(episodeLength_); + name_ = getParam("Problem.Name"); } /*! @@ -177,19 +160,6 @@ public: const std::string name() const { return name_; } - /*! - * \brief Returns the source term at specific position in the domain. - * - * \param values The source values for the primary variables - * \param globalPos The position of the center of the finite volume - */ - PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const - { - return PrimaryVariables(0.0); - } - - - // \} /*! @@ -207,7 +177,7 @@ public: BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const { BoundaryTypes bcTypes; - if(globalPos[0] > this->bBoxMax()[0] - eps_) + if(globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_) bcTypes.setAllDirichlet(); else bcTypes.setAllNeumann(); @@ -242,9 +212,9 @@ public: * in normal direction of each phase. Negative values mean influx. */ NeumannFluxes neumann(const Element& element, - const FVElementGeometry& fvGeometry, - const ElementVolumeVariables& elemVolVars, - const SubControlVolumeFace& scvf) const + const FVElementGeometry& fvGeometry, + const ElementVolumeVariables& elemVolVars, + const SubControlVolumeFace& scvf) const { PrimaryVariables values(0.0); const auto& globalPos = scvf.ipGlobal(); @@ -286,56 +256,31 @@ public: * from the solution of the current time step to the VTK * writer. Adjust this in case of anisotropic permeabilities. */ - void addOutputVtkFields() + template + void addVtkFields(VTKWriter& vtk) { - // get the number of degrees of freedom - unsigned numDofs = this->model().numDofs(); - - // create the scalar field required for the permeabilities - typedef Dune::BlockVector > ScalarField; - ScalarField *Kxx = this->resultWriter().allocateManagedBuffer(numDofs); - - FVElementGeometry fvGeometry; + const auto& gg = this->fvGridGeometry(); + Kxx_.resize(gg.numDofs()); + vtk.addField(Kxx_, "permeability"); for (const auto& element : elements(this->gridView())) { - fvGeometry.update(this->gridView(), element); + auto fvGeometry = localView(gg); + fvGeometry.bindElement(element); - for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx) - { - int dofIdxGlobal = this->model().dofMapper().subIndex(element, scvIdx, dofCodim); - (*Kxx)[dofIdxGlobal] = this->spatialParams().intrinsicPermeability(element, fvGeometry, scvIdx); - } + for (const auto& scv : scvs(fvGeometry)) + Kxx_[scv.dofIndex()] = this->spatialParams().intrinsicPermeabilityAtPos(scv.dofPosition()); } - - this->resultWriter().attachDofData(*Kxx, "permeability", isBox); - } - - bool shouldWriteOutput() const - { - return this->timeManager().timeStepIndex() == 0 || - this->timeManager().episodeWillBeFinished() || - this->timeManager().willBeFinished(); - } - - void episodeEnd() - { - this->timeManager().startNextEpisode(episodeLength_); - } - - bool shouldWriteRestartFile() const - { - return false; } private: // checks, whether a point is located inside the contamination zone bool isInContaminationZone(const GlobalPosition &globalPos) const { - return (Dune::FloatCmp::ge(globalPos[0],0.2) - && Dune::FloatCmp::le(globalPos[0],0.8) - && Dune::FloatCmp::ge(globalPos[1],0.4) - && Dune::FloatCmp::le(globalPos[1],0.65)); + return (Dune::FloatCmp::ge(globalPos[0], 0.2) + && Dune::FloatCmp::le(globalPos[0], 0.8) + && Dune::FloatCmp::ge(globalPos[1], 0.4) + && Dune::FloatCmp::le(globalPos[1], 0.65)); } // internal method for the initial condition (reused for the @@ -346,7 +291,7 @@ private: if (isInContaminationZone(globalPos)) values.setState(threePhases); else - values.setState(wgPhaseOnly); + values.setState(wgPhaseOnly); values[pressureIdx] = 1e5 ; values[switch1Idx] = 0.12; @@ -357,13 +302,14 @@ private: { values[switch2Idx] = 0.07; } - return values; + return values; } static constexpr Scalar eps_ = 1e-6; std::string name_; - Scalar episodeLength_; + std::vector Kxx_; }; -} //end namespace + +} //end namespace Dumux #endif diff --git a/test/porousmediumflow/3p3c/implicit/kuevettespatialparams.hh b/test/porousmediumflow/3p3c/implicit/kuevettespatialparams.hh index 946ee1a21958ea0349e9ef982cdd8c78fe849070..316810e456ae7ed4f75bef42c4f62763cf2a724a 100644 --- a/test/porousmediumflow/3p3c/implicit/kuevettespatialparams.hh +++ b/test/porousmediumflow/3p3c/implicit/kuevettespatialparams.hh @@ -53,13 +53,12 @@ SET_PROP(KuevetteSpatialParams, 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); public: // define the material law parameterized by absolute saturations - typedef EffToAbsLaw type; + using type = EffToAbsLaw>; }; -} +} // end namespace Dumux /*! * \ingroup ThreePThreeCModel @@ -69,42 +68,31 @@ SET_PROP(KuevetteSpatialParams, MaterialLaw) template class KuevetteSpatialParams : public ImplicitSpatialParams { - typedef ImplicitSpatialParams ParentType; + using ParentType = ImplicitSpatialParams; + + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); + using Indices = typename GET_PROP_TYPE(TypeTag, Indices); - typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - typedef typename Grid::ctype CoordScalar; enum { dim=GridView::dimension, dimWorld=GridView::dimensionworld }; - typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; enum { wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx }; - typedef Dune::FieldVector GlobalPosition; - typedef Dune::FieldVector DimVector; - - - typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector; - - typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables; - typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables; - - typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry; - typedef typename GridView::template Codim<0>::Entity Element; + using Element = typename GridView::template Codim<0>::Entity; using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume); using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector); using Problem = typename GET_PROP_TYPE(TypeTag, Problem); - + using GlobalPosition = Dune::FieldVector; public: - 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; using PermeabilityType = Scalar; /*! @@ -112,8 +100,8 @@ public: * * \param gridView The grid view */ - KuevetteSpatialParams(const Problem& problem, const GridView &gridView) - : ParentType(problem, gridView) + KuevetteSpatialParams(const Problem& problem) + : ParentType(problem) { // intrinsic permeabilities fineK_ = 6.28e-12; @@ -150,28 +138,20 @@ public: fineMaterialParams_.setRhoBulk(1500.); } - ~KuevetteSpatialParams() - {} - - /*! - * \brief Update the spatial parameters with the flow solution - * after a timestep. + * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$ + * \note It is possibly solution dependent. * - * \param globalSolution The global solution vector + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \return permeability */ - void update(const SolutionVector &globalSolution) - { - } - - /*! - * \brief Apply the intrinsic permeability tensor to a pressure - * potential gradient. - * - * \param globalPos The global position - */ - Scalar permeabilityAtPos(const GlobalPosition& globalPos) const + PermeabilityType permeability(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const { + const auto& globalPos = scv.dofPosition(); if (isFineMaterial_(globalPos)) return fineK_; return coarseK_; @@ -180,13 +160,15 @@ public: /*! * \brief Define the porosity \f$[-]\f$ of the spatial parameters * - * \param element The finite element - * \param fvGeometry The finite volume geometry - * \param scvIdx The local index of the sub-control volume where - * the porosity needs to be defined + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. */ - Scalar porosityAtPos(const GlobalPosition& globalPos) const + Scalar porosity(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const { + const auto& globalPos = scv.dofPosition(); if (isFineMaterial_(globalPos)) return finePorosity_; else @@ -195,14 +177,18 @@ public: /*! - * \brief return the parameter object for the Brooks-Corey material law which depends on the position + * \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.). * - * \param element The current finite element - * \param fvGeometry The current finite volume geometry of the element - * \param scvIdx The index of the sub-control volume + * \param element The current element + * \param scv The sub-control volume inside the element. + * \param elemSol The solution at the dofs connected to the element. + * \return the material parameters object */ - const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const + const MaterialLawParams& materialLawParams(const Element& element, + const SubControlVolume& scv, + const ElementSolutionVector& elemSol) const { + const auto& globalPos = scv.dofPosition(); if (isFineMaterial_(globalPos)) return fineMaterialParams_; else @@ -257,12 +243,12 @@ public: private: bool isFineMaterial_(const GlobalPosition &globalPos) const { - return ((Dune::FloatCmp::ge(globalPos[0],0.13) - && Dune::FloatCmp::le(globalPos[0],1.20) - && Dune::FloatCmp::ge(globalPos[1],0.32) - && Dune::FloatCmp::le(globalPos[1],0.57)) - || (Dune::FloatCmp::ge(globalPos[0],1.20) - && Dune::FloatCmp::le(globalPos[1],0.15))); + return ((Dune::FloatCmp::ge(globalPos[0], 0.13) + && Dune::FloatCmp::le(globalPos[0], 1.24) + && Dune::FloatCmp::ge(globalPos[1], 0.32) + && Dune::FloatCmp::le(globalPos[1], 0.60)) + || (Dune::FloatCmp::ge(globalPos[0], 1.20) + && Dune::FloatCmp::le(globalPos[1], 0.15))); } Scalar fineK_; @@ -277,6 +263,6 @@ private: Scalar lambdaSolid_; }; -} +} // end namespace Dumux #endif diff --git a/test/porousmediumflow/3p3c/implicit/test_3p3c_fv.cc b/test/porousmediumflow/3p3c/implicit/test_3p3c_fv.cc new file mode 100644 index 0000000000000000000000000000000000000000..00f2252c5961121bdbc1cc229ee0f5ac984c4e7f --- /dev/null +++ b/test/porousmediumflow/3p3c/implicit/test_3p3c_fv.cc @@ -0,0 +1,254 @@ +// -*- 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 Test for the three-phase three-component box model + */ +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include + +#include "infiltration3p3cproblem.hh" +#include "kuevetteproblem.hh" +#include "columnxylolproblem.hh" + +/*! + * \brief Provides an interface for customizing error messages associated with + * reading in parameters. + * + * \param progName The name of the program, that was tried to be started. + * \param errorMsg The error message that was issued by the start function. + * Comprises the thing that went wrong and a general help message. + */ +void usage(const char *progName, const std::string &errorMsg) +{ + if (errorMsg.size() > 0) { + std::string errorMessageOut = "\nUsage: "; + errorMessageOut += progName; + errorMessageOut += " [options]\n"; + errorMessageOut += errorMsg; + errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" + "\t-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) try +{ + + using namespace Dumux; + + // define the type tag for this problem + using TypeTag = TTAG(TYPETAG); + + //////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////// + + // initialize MPI, finalize is done automatically on exit + const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv); + + // print dumux start message + if (mpiHelper.rank() == 0) + DumuxMessage::print(/*firstCall=*/true); + + // parse command line arguments and input file + Parameters::init(argc, argv, usage); + + // try to create a grid (from the given grid file or the input file) + 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(fvGridGeometry->numDofs()); + 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); + timeLoop->setPeriodicCheckPoint(getParam("TimeLoop.EpisodeLength", std::numeric_limits::max())); + + // 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->dofMapper()); + + // the non-linear solver + using NewtonController = PriVarSwitchNewtonController; + using NewtonMethod = Dumux::NewtonMethod; + auto newtonController = std::make_shared(leafGridView.comm(), timeLoop); + NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver); + + // time loop + timeLoop->start(); + while (!timeLoop->finished()) + { + // 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 + // if episode length was specificied output only at the end of episodes + if (!haveParam("TimeLoop.EpisodeLength") || timeLoop->isCheckPoint() || timeLoop->finished() || timeLoop->timeStepIndex() == 1) + 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())); + } + + 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/3p3c/implicit/test_box3p3c_reference.input b/test/porousmediumflow/3p3c/implicit/test_3p3c_fv.input similarity index 76% rename from test/porousmediumflow/3p3c/implicit/test_box3p3c_reference.input rename to test/porousmediumflow/3p3c/implicit/test_3p3c_fv.input index 26470f2cb7214b99a460b52c31b97c136ee57aff..ef445b1a852b728898f0ae79c34506e0d9d8c446 100644 --- a/test/porousmediumflow/3p3c/implicit/test_box3p3c_reference.input +++ b/test/porousmediumflow/3p3c/implicit/test_3p3c_fv.input @@ -1,4 +1,4 @@ -[TimeManager] +[TimeLoop] DtInitial = 86400 # [s] TEnd = 864000 # [s] @@ -7,8 +7,10 @@ UpperRight = 500 10 Cells = 40 4 [Problem] -Name = box3p3c +Name = infiltration [Implicit] NumericDifferenceMethod = 0 # -1 backward differences, 0: central differences, +1: forward differences +[Newton] +MaxRelativeShift = 1e-6 diff --git a/test/porousmediumflow/3p3c/implicit/test_box3p3c.cc b/test/porousmediumflow/3p3c/implicit/test_box3p3c.cc deleted file mode 100644 index ba46d56229dbd6b939c0c072c857f658ec2fe5c1..0000000000000000000000000000000000000000 --- a/test/porousmediumflow/3p3c/implicit/test_box3p3c.cc +++ /dev/null @@ -1,59 +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 . * - *****************************************************************************/ -/*! - * \file - * - * \brief Test for the three-phase three-component box model - */ -#include -#include "infiltrationproblem.hh" -#include - -/*! - * \brief Provides an interface for customizing error messages associated with - * reading in parameters. - * - * \param progName The name of the program, that was tried to be started. - * \param errorMsg The error message that was issued by the start function. - * Comprises the thing that went wrong and a general help message. - */ -void usage(const char *progName, const std::string &errorMsg) -{ - if (errorMsg.size() > 0) { - std::string errorMessageOut = "\nUsage: "; - errorMessageOut += progName; - errorMessageOut += " [options]\n"; - errorMessageOut += errorMsg; - errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" - "\t-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) -{ - typedef TTAG(InfiltrationBoxProblem) ProblemTypeTag; - return Dumux::start(argc, argv, usage); -} - diff --git a/test/porousmediumflow/3p3c/implicit/test_box3p3c.input b/test/porousmediumflow/3p3c/implicit/test_box3p3c.input deleted file mode 100644 index 80639b17c803e26b67c7a4061dbfe05ed3b7ca1d..0000000000000000000000000000000000000000 --- a/test/porousmediumflow/3p3c/implicit/test_box3p3c.input +++ /dev/null @@ -1,14 +0,0 @@ -[TimeManager] -DtInitial = 60 # [s] -TEnd = 864000 # [s] - -[Grid] -UpperRight = 500 10 -Cells = 250 10 - -[Problem] -Name = infiltrationbox - -[Implicit] -NumericDifferenceMethod = 0 # -1 backward differences, 0: central differences, +1: forward differences - diff --git a/test/porousmediumflow/3p3c/implicit/test_box3p3cnicolumnxylol.cc b/test/porousmediumflow/3p3c/implicit/test_box3p3cnicolumnxylol.cc deleted file mode 100644 index bda25900291e159d97fbdf8726d01c52fda6842d..0000000000000000000000000000000000000000 --- a/test/porousmediumflow/3p3c/implicit/test_box3p3cnicolumnxylol.cc +++ /dev/null @@ -1,59 +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 . * - *****************************************************************************/ -/*! - * \file - * - * \brief Test for the non-isothermal three-phase three-component box model. - */ -#include -#include "columnxylolproblem.hh" -#include - - -/*! - * \brief Provides an interface for customizing error messages associated with - * reading in parameters. - * - * \param progName The name of the program, that was tried to be started. - * \param errorMsg The error message that was issued by the start function. - * Comprises the thing that went wrong and a general help message. - */ -void usage(const char *progName, const std::string &errorMsg) -{ - if (errorMsg.size() > 0) { - std::string errorMessageOut = "\nUsage: "; - errorMessageOut += progName; - errorMessageOut += " [options]\n"; - errorMessageOut += errorMsg; - errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" - "\t-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) -{ - typedef TTAG(ColumnBoxProblem) ProblemTypeTag; - return Dumux::start(argc, argv, usage); -} diff --git a/test/porousmediumflow/3p3c/implicit/test_box3p3cnicolumnxylol.input b/test/porousmediumflow/3p3c/implicit/test_box3p3cnicolumnxylol.input deleted file mode 100644 index bed5d51223df232d4315240d586d8e35b04d9bac..0000000000000000000000000000000000000000 --- a/test/porousmediumflow/3p3c/implicit/test_box3p3cnicolumnxylol.input +++ /dev/null @@ -1,16 +0,0 @@ -[TimeManager] -DtInitial = 1 # [s] -TEnd = 200 # [s] - -[Grid] -UpperRight = 0.1 1.20 -Cells = 1 120 - -[Problem] -Name = columnxylolbox # name passed to the output routines - -[Implicit] -NumericDifferenceMethod = 0 # -1 backward differences, 0: central differences, +1: forward differences - -[TimeManager] -MaxTimeStepSize= 4.0 # Set the maximum time step diff --git a/test/porousmediumflow/3p3c/implicit/test_box3p3cnikuevette.cc b/test/porousmediumflow/3p3c/implicit/test_box3p3cnikuevette.cc deleted file mode 100644 index 0d7fe2aea06d64da35561f337a61c233da1b2a61..0000000000000000000000000000000000000000 --- a/test/porousmediumflow/3p3c/implicit/test_box3p3cnikuevette.cc +++ /dev/null @@ -1,59 +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 . * - *****************************************************************************/ -/*! - * \file - * - * \brief Test for the non-isothermal three-phase three-component box model. - */ -#include -#include "kuevetteproblem.hh" -#include - - -/*! - * \brief Provides an interface for customizing error messages associated with - * reading in parameters. - * - * \param progName The name of the program, that was tried to be started. - * \param errorMsg The error message that was issued by the start function. - * Comprises the thing that went wrong and a general help message. - */ -void usage(const char *progName, const std::string &errorMsg) -{ - if (errorMsg.size() > 0) { - std::string errorMessageOut = "\nUsage: "; - errorMessageOut += progName; - errorMessageOut += " [options]\n"; - errorMessageOut += errorMsg; - errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" - "\t-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) -{ - typedef TTAG(KuevetteBoxProblem) ProblemTypeTag; - return Dumux::start(argc, argv, usage); -} diff --git a/test/porousmediumflow/3p3c/implicit/test_box3p3cnikuevette.input b/test/porousmediumflow/3p3c/implicit/test_box3p3cnikuevette.input deleted file mode 100644 index 3e844993aa30d88519714b2b14097eb095edb1d5..0000000000000000000000000000000000000000 --- a/test/porousmediumflow/3p3c/implicit/test_box3p3cnikuevette.input +++ /dev/null @@ -1,15 +0,0 @@ -[TimeManager] -DtInitial = 1 # [s] -TEnd = 2700 # [s] -MaxTimeStepSize = 60 # [s] -EpisodeLength = 1000 # [s] - -[Grid] -UpperRight = 1.5 0.74 -Cells = 15 8 - -[Problem] -Name = kuevettebox # name passed to the output routines - -[Implicit] -NumericDifferenceMethod = 0 # use central differences (backward -1, forward +1) diff --git a/test/porousmediumflow/3p3c/implicit/test_cc3p3c.cc b/test/porousmediumflow/3p3c/implicit/test_cc3p3c.cc deleted file mode 100644 index 6178da9013a10eb18b37f7034ac44685b6fa0f03..0000000000000000000000000000000000000000 --- a/test/porousmediumflow/3p3c/implicit/test_cc3p3c.cc +++ /dev/null @@ -1,59 +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 . * - *****************************************************************************/ -/*! - * \file - * - * \brief Test for the three-phase three-component CC model - */ -#include -#include "infiltrationproblem.hh" -#include - -/*! - * \brief Provides an interface for customizing error messages associated with - * reading in parameters. - * - * \param progName The name of the program, that was tried to be started. - * \param errorMsg The error message that was issued by the start function. - * Comprises the thing that went wrong and a general help message. - */ -void usage(const char *progName, const std::string &errorMsg) -{ - if (errorMsg.size() > 0) { - std::string errorMessageOut = "\nUsage: "; - errorMessageOut += progName; - errorMessageOut += " [options]\n"; - errorMessageOut += errorMsg; - errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" - "\t-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) -{ - typedef TTAG(InfiltrationCCProblem) ProblemTypeTag; - return Dumux::start(argc, argv, usage); -} - diff --git a/test/porousmediumflow/3p3c/implicit/test_cc3p3c.input b/test/porousmediumflow/3p3c/implicit/test_cc3p3c.input deleted file mode 100644 index 9db6dcf902a4eac522f5dbdfc3280f0bc94af0bd..0000000000000000000000000000000000000000 --- a/test/porousmediumflow/3p3c/implicit/test_cc3p3c.input +++ /dev/null @@ -1,17 +0,0 @@ -[TimeManager] -DtInitial = 60 # [s] -TEnd = 864000 # [s] - -[Grid] -UpperRight = 500 10 -Cells = 250 10 - -[Problem] -Name = infiltrationcc - -[Implicit] -NumericDifferenceMethod = 0 # -1 backward differences, 0: central differences, +1: forward differences - -[Vtk] -AddPermeability = true -AddPorosity = true diff --git a/test/porousmediumflow/3p3c/implicit/test_cc3p3c_reference.input b/test/porousmediumflow/3p3c/implicit/test_cc3p3c_reference.input deleted file mode 100644 index 60cdfce4700b8c6a22ea88569cfddc5bb4dc45fa..0000000000000000000000000000000000000000 --- a/test/porousmediumflow/3p3c/implicit/test_cc3p3c_reference.input +++ /dev/null @@ -1,17 +0,0 @@ -[TimeManager] -DtInitial = 86400 # [s] -TEnd = 864000 # [s] - -[Grid] -UpperRight = 500 10 -Cells = 40 4 - -[Problem] -Name = cc3p3c - -[Implicit] -NumericDifferenceMethod = 0 # -1 backward differences, 0: central differences, +1: forward differences - -[Vtk] -AddPermeability = true -AddPorosity = true diff --git a/test/porousmediumflow/3p3c/implicit/test_cc3p3cnicolumnxylol.cc b/test/porousmediumflow/3p3c/implicit/test_cc3p3cnicolumnxylol.cc deleted file mode 100644 index 96bcd687b17e7ab05b6b7a8734f6f46e99f74663..0000000000000000000000000000000000000000 --- a/test/porousmediumflow/3p3c/implicit/test_cc3p3cnicolumnxylol.cc +++ /dev/null @@ -1,59 +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 . * - *****************************************************************************/ -/*! - * \file - * - * \brief Test for the non-isothermal three-phase three-component CC model. - */ -#include -#include "columnxylolproblem.hh" -#include - - -/*! - * \brief Provides an interface for customizing error messages associated with - * reading in parameters. - * - * \param progName The name of the program, that was tried to be started. - * \param errorMsg The error message that was issued by the start function. - * Comprises the thing that went wrong and a general help message. - */ -void usage(const char *progName, const std::string &errorMsg) -{ - if (errorMsg.size() > 0) { - std::string errorMessageOut = "\nUsage: "; - errorMessageOut += progName; - errorMessageOut += " [options]\n"; - errorMessageOut += errorMsg; - errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" - "\t-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) -{ - typedef TTAG(ColumnCCProblem) ProblemTypeTag; - return Dumux::start(argc, argv, usage); -} diff --git a/test/porousmediumflow/3p3c/implicit/test_cc3p3cnikuevette.cc b/test/porousmediumflow/3p3c/implicit/test_cc3p3cnikuevette.cc deleted file mode 100644 index 8a0abb015d4c61fdb07fd89d1d3507c75a6a441e..0000000000000000000000000000000000000000 --- a/test/porousmediumflow/3p3c/implicit/test_cc3p3cnikuevette.cc +++ /dev/null @@ -1,59 +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 . * - *****************************************************************************/ -/*! - * \file - * - * \brief Test for the non-isothermal three-phase three-component CC model. - */ -#include -#include "kuevetteproblem.hh" -#include - - -/*! - * \brief Provides an interface for customizing error messages associated with - * reading in parameters. - * - * \param progName The name of the program, that was tried to be started. - * \param errorMsg The error message that was issued by the start function. - * Comprises the thing that went wrong and a general help message. - */ -void usage(const char *progName, const std::string &errorMsg) -{ - if (errorMsg.size() > 0) { - std::string errorMessageOut = "\nUsage: "; - errorMessageOut += progName; - errorMessageOut += " [options]\n"; - errorMessageOut += errorMsg; - errorMessageOut += "\n\nThe list of mandatory options for this program is:\n" - "\t-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) -{ - typedef TTAG(KuevetteCCProblem) ProblemTypeTag; - return Dumux::start(argc, argv, usage); -} diff --git a/test/porousmediumflow/3p3c/implicit/test_cc3p3cnicolumnxylol.input b/test/porousmediumflow/3p3c/implicit/test_columnxylol_fv.input similarity index 60% rename from test/porousmediumflow/3p3c/implicit/test_cc3p3cnicolumnxylol.input rename to test/porousmediumflow/3p3c/implicit/test_columnxylol_fv.input index 73b95afdfb9c2dee5e284868e95f1171714f15b0..c1848acdb19e54f7bc697e71ff5daf79569d8ea3 100644 --- a/test/porousmediumflow/3p3c/implicit/test_cc3p3cnicolumnxylol.input +++ b/test/porousmediumflow/3p3c/implicit/test_columnxylol_fv.input @@ -1,16 +1,14 @@ -[TimeManager] +[TimeLoop] DtInitial = 1 # [s] TEnd = 200 # [s] +MaxTimeStepSize= 4 # Set the maximum time step size [Grid] UpperRight = 0.1 1.20 Cells = 1 120 [Problem] -Name = columnxylolcc # name passed to the output routines +Name = columnxylol # name passed to the output routines [Implicit] NumericDifferenceMethod = 0 # -1 backward differences, 0: central differences, +1: forward differences - -[TimeManager] -MaxTimeStepSize= 4. # Set the maximum time step diff --git a/test/porousmediumflow/3p3c/implicit/test_cc3p3cnikuevette.input b/test/porousmediumflow/3p3c/implicit/test_kuvette_fv.input similarity index 70% rename from test/porousmediumflow/3p3c/implicit/test_cc3p3cnikuevette.input rename to test/porousmediumflow/3p3c/implicit/test_kuvette_fv.input index ffe5648570eb4fab0f6feab5dab88bd63a6946cc..1738b0ae022d312f8d7a7c9cc4c9e6077951bd03 100644 --- a/test/porousmediumflow/3p3c/implicit/test_cc3p3cnikuevette.input +++ b/test/porousmediumflow/3p3c/implicit/test_kuvette_fv.input @@ -1,4 +1,4 @@ -[TimeManager] +[TimeLoop] DtInitial = 1 # [s] TEnd = 2700 # [s] MaxTimeStepSize = 60 # [s] @@ -9,7 +9,10 @@ UpperRight = 1.5 0.74 Cells = 15 8 [Problem] -Name = kuevettecc # name passed to the output routines +Name = kuevette # name passed to the output routines [Implicit] NumericDifferenceMethod = 0 # use central differences (backward -1, forward +1) + +[Newton] +MaxRelativeShift = 1e-6 diff --git a/test/references/infiltrationbox-reference.vtu b/test/references/infiltrationbox-reference.vtu index 6960151100cb4d38f1ec173d1ec662ca27eeea90..04c14f1c15cfd77e598b1f3967f117aa2d2603ba 100644 --- a/test/references/infiltrationbox-reference.vtu +++ b/test/references/infiltrationbox-reference.vtu @@ -383,26 +383,6 @@ 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 - - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 - 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 diff --git a/test/references/infiltrationcc-reference.vtu b/test/references/infiltrationcc-reference.vtu index b14a02e2cf10dc84a8213f07713584a61413fad1..977da799a25e7db5fcc152a3bc86d08ddcf3d9d7 100644 --- a/test/references/infiltrationcc-reference.vtu +++ b/test/references/infiltrationcc-reference.vtu @@ -307,22 +307,6 @@ 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 - - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 - 283.15 283.15 283.15 283.15 - 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6