diff --git a/CMakeLists.txt b/CMakeLists.txt index f88ea3925..febe89342 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -210,6 +210,7 @@ foreach(tapp co2injection_flash_ni_vcfv co2injection_ncp_ni_ecfv co2injection_pvs_ni_ecfv co2_ptflash_ecfv + co2_ptflash_ecfv_validation powerinjection_forchheimer_fd powerinjection_forchheimer_ad powerinjection_darcy_fd diff --git a/tests/co2_ptflash_ecfv_validation.cc b/tests/co2_ptflash_ecfv_validation.cc new file mode 100644 index 000000000..4dac736bc --- /dev/null +++ b/tests/co2_ptflash_ecfv_validation.cc @@ -0,0 +1,62 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM 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. + + OPM 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 OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief Box problem with two phases and multiple components. + * Solved with a PTFlash two phase solver. + */ +#include "config.h" + +#include +#include "problems/co2ptflashproblemvalidation.hh" + + +namespace Opm::Properties { + +namespace TTag { + struct CO2PTEcfvProblemValidation { + using InheritsFrom = std::tuple; +}; +} + +template +struct SpatialDiscretizationSplice +{ + using type = TTag::EcfvDiscretization; +}; + +template +struct LocalLinearizerSplice +{ + using type = TTag::AutoDiffLocalLinearizer; +}; + + +} // namespace Opm::Properties + +int main(int argc, char **argv) +{ + using EcfvProblemTypeTag = Opm::Properties::TTag::CO2PTEcfvProblemValidation; + return Opm::start(argc, argv); +} diff --git a/tests/problems/co2ptflashproblem.hh b/tests/problems/co2ptflashproblem.hh index 204115c35..7cffe0d6a 100644 --- a/tests/problems/co2ptflashproblem.hh +++ b/tests/problems/co2ptflashproblem.hh @@ -218,6 +218,26 @@ struct VtkWriteFilterVelocities { static constexpr bool value = true; }; +template +struct VtkWriteViscosities { + static constexpr bool value = true; +}; + +template +struct VtkWritePorosity { + static constexpr bool value = true; +}; + +template +struct VtkWriteIntrinsicPermeabilities { + static constexpr bool value = true; +}; + +template +struct VtkWriteMobilities { + static constexpr bool value = true; +}; + template struct VtkWritePotentialGradients { static constexpr bool value = true; diff --git a/tests/problems/co2ptflashproblemvalidation.hh b/tests/problems/co2ptflashproblemvalidation.hh new file mode 100644 index 000000000..f2b5baf65 --- /dev/null +++ b/tests/problems/co2ptflashproblemvalidation.hh @@ -0,0 +1,614 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM 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. + + OPM 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 OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \copydoc Opm::co2ptflashproblem + */ +#ifndef OPM_CO2PTFLASH_PROBLEM_HH +#define OPM_CO2PTFLASH_PROBLEM_HH + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace Opm { +template +class CO2PTProblem; +} // namespace Opm */ + +namespace Opm::Properties { + +namespace TTag { +struct CO2PTBaseProblem {}; +} // end namespace TTag + +template +struct Temperature { using type = UndefinedProperty; }; +template +struct SimulationName { using type = UndefinedProperty; }; +template +struct EpisodeLength { using type = UndefinedProperty;}; + +template +struct Initialpressure { using type = UndefinedProperty;}; + +// Set the grid type: --->2D +template +struct Grid { using type = Dune::YaspGrid; }; + +// Set the problem property +template +struct Problem +{ using type = Opm::CO2PTProblem; }; + +// Set flash solver +template +struct FlashSolver { +private: + using Scalar = GetPropType; + using FluidSystem = GetPropType; + using Evaluation = GetPropType; + +public: + using type = Opm::PTFlash; +}; + +// Set fluid configuration +template +struct FluidSystem +{ +private: + using Scalar = GetPropType; + +public: + using type = Opm::ThreeComponentFluidSystem; +}; + +// Set the material Law +template +struct MaterialLaw { +private: + using FluidSystem = GetPropType; + enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; + enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; + + using Scalar = GetPropType; + using Traits = Opm::TwoPhaseMaterialTraits; + + // define the material law which is parameterized by effective saturation + using EffMaterialLaw = Opm::NullMaterial; + //using EffMaterialLaw = Opm::BrooksCorey; + +public: + using type = EffMaterialLaw; +}; + +// Write the Newton convergence behavior to disk? +template +struct NewtonWriteConvergence { +static constexpr bool value = false; }; + +// Enable gravity false +template +struct EnableGravity { static constexpr bool value = false; +}; + +// set the defaults for the problem specific properties + template + struct Temperature { + using type = GetPropType; + static constexpr type value = 423.25;//TODO + }; + +template +struct Initialpressure { + using type = GetPropType; + static constexpr type value = 75.e5; +}; + +template +struct SimulationName { + static constexpr auto value = "co2_ptflash"; +}; + +// The default for the end time of the simulation +template +struct EndTime { + using type = GetPropType; + static constexpr type value = 1.296e8/2; +}; + +// this is kinds of telling the report step length +template +struct EpisodeLength { + using type = GetPropType; + static constexpr type value = 5*21600.0; +}; + +// convergence control +template +struct InitialTimeStepSize { + using type = GetPropType; + static constexpr type value = 21600.0; +}; + +template +struct LinearSolverTolerance { + using type = GetPropType; + static constexpr type value = 1e-3; +}; + +template +struct LinearSolverAbsTolerance { + using type = GetPropType; + static constexpr type value = 0.; +}; + +template +struct NewtonTolerance { + using type = GetPropType; + static constexpr type value = 1e-3; +}; + +template +struct NewtonMaxIterations { + using type = GetPropType; + static constexpr type value = 30; +}; + +template +struct NewtonTargetIterations { + using type = GetPropType; + static constexpr type value = 6; +}; + +// output +template +struct VtkWriteFilterVelocities { + static constexpr bool value = true; +}; + +template +struct VtkWritePotentialGradients { + static constexpr bool value = true; +}; + +template +struct VtkWriteTotalMassFractions { + static constexpr bool value = true; +}; + +template +struct VtkWriteTotalMoleFractions { + static constexpr bool value = true; +}; + +template +struct VtkWriteFugacityCoeffs { + static constexpr bool value = true; +}; + +template +struct VtkWriteLiquidMoleFractions { + static constexpr bool value = true; +}; + +template +struct VtkWriteEquilibriumConstants { + static constexpr bool value = true; +}; + +// mesh grid +template +struct Vanguard { + using type = Opm::StructuredGridVanguard; +}; + +template +struct DomainSizeX { + using type = GetPropType; + static constexpr type value = 1000; // meter +}; + +template +struct DomainSizeY { + using type = GetPropType; + static constexpr type value = 1.0; +}; + +// DomainSizeZ is not needed, while to keep structuredgridvanguard.hh compile +template +struct DomainSizeZ { + using type = GetPropType; + static constexpr type value = 0.1; +}; + +template +struct CellsX { static constexpr int value = 1000; }; +template +struct CellsY { static constexpr int value = 1; }; +// CellsZ is not needed, while to keep structuredgridvanguard.hh compile +template +struct CellsZ { static constexpr int value = 1; }; + +template +struct EnableEnergy { + static constexpr bool value = false; +}; + +} // namespace Opm::Properties + + + + +namespace Opm { +/*! + * \ingroup TestProblems + * + * \brief 3 component simple testproblem with ["CO2", "C1", "C10"] + * + */ +template +class CO2PTProblem : public GetPropType +{ + using ParentType = GetPropType; + + using Scalar = GetPropType; + using Evaluation = GetPropType; + using GridView = GetPropType; + using FluidSystem = GetPropType; + enum { dim = GridView::dimension }; + enum { dimWorld = GridView::dimensionworld }; + using Indices = GetPropType; + using PrimaryVariables = GetPropType; + using RateVector = GetPropType; + using BoundaryRateVector = GetPropType; + using MaterialLaw = GetPropType; + using Simulator = GetPropType; + using Model = GetPropType; + using MaterialLawParams = GetPropType; + + using Toolbox = Opm::MathToolbox; + using CoordScalar = typename GridView::ctype; + + enum { numPhases = FluidSystem::numPhases }; + enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; + enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; + enum { Comp2Idx = FluidSystem::Comp2Idx }; + enum { Comp1Idx = FluidSystem::Comp1Idx }; + enum { Comp0Idx = FluidSystem::Comp0Idx }; + enum { conti0EqIdx = Indices::conti0EqIdx }; + enum { contiCO2EqIdx = conti0EqIdx + Comp1Idx }; + enum { numComponents = getPropValue() }; + enum { enableEnergy = getPropValue() }; + enum { enableDiffusion = getPropValue() }; + enum { enableGravity = getPropValue() }; + + using GlobalPosition = Dune::FieldVector; + using DimMatrix = Dune::FieldMatrix; + using DimVector = Dune::FieldVector; + using ComponentVector = Dune::FieldVector; + using FlashSolver = GetPropType; + +public: + using FluidState = Opm::CompositionalFluidState; + /*! + * \copydoc Doxygen::defaultProblemConstructor + */ + explicit CO2PTProblem(Simulator& simulator) + : ParentType(simulator) + { + const Scalar epi_len = EWOMS_GET_PARAM(TypeTag, Scalar, EpisodeLength); + simulator.setEpisodeLength(epi_len); + } + + void initPetrophysics() + { + temperature_ = EWOMS_GET_PARAM(TypeTag, Scalar, Temperature); + K_ = this->toDimMatrix_(9.869232667160131e-14); + + porosity_ = 0.25; + } + + template + const DimVector& + gravity([[maybe_unused]]const Context& context, [[maybe_unused]] unsigned spaceIdx, [[maybe_unused]] unsigned timeIdx) const + { + return gravity(); + } + + const DimVector& gravity() const + { + return gravity_; + } + + /*! + * \copydoc FvBaseProblem::finishInit + */ + void finishInit() + { + ParentType::finishInit(); + // initialize fixed parameters; temperature, permeability, porosity + initPetrophysics(); + } + + /*! + * \copydoc co2ptflashproblem::registerParameters + */ + static void registerParameters() + { + ParentType::registerParameters(); + + EWOMS_REGISTER_PARAM(TypeTag, Scalar, Temperature, + "The temperature [K] in the reservoir"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, Initialpressure, + "The initial pressure [Pa s] in the reservoir"); + EWOMS_REGISTER_PARAM(TypeTag, + std::string, + SimulationName, + "The name of the simulation used for the output " + "files"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, EpisodeLength, + "Time interval [s] for episode length"); + } + + /*! + * \copydoc FvBaseProblem::name + */ + std::string name() const + { + std::ostringstream oss; + oss << EWOMS_GET_PARAM(TypeTag, std::string, SimulationName); + return oss.str(); + } + + // This method must be overridden for the simulator to continue with + // a new episode. We just start a new episode with the same length as + // the old one. + void endEpisode() + { + Scalar epi_len = EWOMS_GET_PARAM(TypeTag, Scalar, EpisodeLength); + this->simulator().startNextEpisode(epi_len); + } + + // only write output when episodes change, aka. report steps, and + // include the initial timestep too + bool shouldWriteOutput() + { + return this->simulator().episodeWillBeOver() || (this->simulator().timeStepIndex() == -1); + } + + // we don't care about doing restarts from every fifth timestep, it + // will just slow us down + bool shouldWriteRestartFile() + { + return false; + } + + /*! + * \copydoc FvBaseProblem::endTimeStep + */ + void endTimeStep() + { + Scalar tol = this->model().newtonMethod().tolerance() * 1e5; + this->model().checkConservativeness(tol); + + // Calculate storage terms + PrimaryVariables storageO, storageW; + this->model().globalPhaseStorage(storageO, oilPhaseIdx); + + // Calculate storage terms + PrimaryVariables storageL, storageG; + this->model().globalPhaseStorage(storageL, /*phaseIdx=*/0); + this->model().globalPhaseStorage(storageG, /*phaseIdx=*/1); + + // Write mass balance information for rank 0 + // if (this->gridView().comm().rank() == 0) { + // std::cout << "Storage: liquid=[" << storageL << "]" + // << " gas=[" << storageG << "]\n" << std::flush; + // } + } + + /*! + * \copydoc FvBaseProblem::initial + */ + template + void initial(PrimaryVariables& values, const Context& context, unsigned spaceIdx, unsigned timeIdx) const + { + Opm::CompositionalFluidState fs; + initialFluidState(fs, context, spaceIdx, timeIdx); + values.assignNaive(fs); + } + + // Constant temperature + template + Scalar temperature([[maybe_unused]] const Context& context, [[maybe_unused]] unsigned spaceIdx, [[maybe_unused]] unsigned timeIdx) const + { + return temperature_; + } + + + // Constant permeability + template + const DimMatrix& intrinsicPermeability([[maybe_unused]] const Context& context, + [[maybe_unused]] unsigned spaceIdx, + [[maybe_unused]] unsigned timeIdx) const + { + return K_; + } + + // Constant porosity + template + Scalar porosity([[maybe_unused]] const Context& context, [[maybe_unused]] unsigned spaceIdx, [[maybe_unused]] unsigned timeIdx) const + { + int spatialIdx = context.globalSpaceIndex(spaceIdx, timeIdx); + int inj = 0; + int prod = EWOMS_GET_PARAM(TypeTag, unsigned, CellsX) - 1; + if (spatialIdx == inj || spatialIdx == prod) { + return 1000.0; + } else { + return porosity_; + } + } + + /*! + * \copydoc FvBaseMultiPhaseProblem::materialLawParams + */ + template + const MaterialLawParams& materialLawParams([[maybe_unused]] const Context& context, + [[maybe_unused]] unsigned spaceIdx, + [[maybe_unused]] unsigned timeIdx) const + { + return this->mat_; + } + + + // No flow (introduce fake wells instead) + template + void boundary(BoundaryRateVector& values, + const Context& /*context*/, + unsigned /*spaceIdx*/, + unsigned /*timeIdx*/) const + { values.setNoFlow(); } + + // No source terms + template + void source(RateVector& source_rate, + [[maybe_unused]] const Context& context, + [[maybe_unused]] unsigned spaceIdx, + [[maybe_unused]] unsigned timeIdx) const + { + source_rate = Scalar(0.0); + } + +private: + + /*! + * \copydoc FvBaseProblem::initial + */ + template + void initialFluidState(FluidState& fs, const Context& context, unsigned spaceIdx, unsigned timeIdx) const + { + // input file order + // 1000*0.6 -- DECANE + // 1000*0.1 -- CO2 + // 1000*0.3 -- METHANE + // opm order: + // CO2, C1, C10 + int inj = 0; + int prod = EWOMS_GET_PARAM(TypeTag, unsigned, CellsX) - 1; + int spatialIdx = context.globalSpaceIndex(spaceIdx, timeIdx); + ComponentVector comp; + comp[0] = Evaluation::createVariable(0.1, 1); + comp[1] = Evaluation::createVariable(0.3, 2); + comp[2] = 1. - comp[0] - comp[1]; + if (spatialIdx == inj) { + comp[0] = Evaluation::createVariable(0.99, 1); + comp[1] = Evaluation::createVariable(0.01 - 1e-3, 2); + comp[2] = 1. - comp[0] - comp[1]; + } + ComponentVector sat; + sat[0] = 1.0; + sat[1] = 1.0 - sat[0]; + + Scalar p0 = EWOMS_GET_PARAM(TypeTag, Scalar, Initialpressure); + + if (spatialIdx == inj) { + p0 = 100.0e5; + } + if (spatialIdx == prod) { + p0 = 50e5; + } + Evaluation p_init = Evaluation::createVariable(p0, 0); + + fs.setPressure(FluidSystem::oilPhaseIdx, p_init); + fs.setPressure(FluidSystem::gasPhaseIdx, p_init); + + fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp0Idx, comp[0]); + fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp1Idx, comp[1]); + fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp2Idx, comp[2]); + + fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp0Idx, comp[0]); + fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp1Idx, comp[1]); + fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp2Idx, comp[2]); + + // It is used here only for calculate the z + fs.setSaturation(FluidSystem::oilPhaseIdx, sat[0]); + fs.setSaturation(FluidSystem::gasPhaseIdx, sat[1]); + + fs.setTemperature(temperature_); + + // ParameterCache paramCache; + { + typename FluidSystem::template ParameterCache paramCache; + paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx); + paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx); + fs.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx)); + fs.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx)); + fs.setViscosity(FluidSystem::oilPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::oilPhaseIdx)); + fs.setViscosity(FluidSystem::gasPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::gasPhaseIdx)); + } + + // Set initial K and L + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + const Evaluation Ktmp = fs.wilsonK_(compIdx); + fs.setKvalue(compIdx, Ktmp); + } + + const Evaluation& Ltmp = -1.0; + fs.setLvalue(Ltmp); + } + + DimMatrix K_; + Scalar porosity_; + Scalar temperature_; + MaterialLawParams mat_; + DimVector gravity_; +}; +} // namespace Opm + +#endif