diff --git a/examples/problems/co2injectiontestflash.hh b/examples/problems/co2injectiontestflash.hh new file mode 100644 index 000000000..89f53674e --- /dev/null +++ b/examples/problems/co2injectiontestflash.hh @@ -0,0 +1,66 @@ +// -*- 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::Co2InjectionTestFlash + */ +#ifndef EWOMS_CO2_INJECTION_TEST_FLASH_HH +#define EWOMS_CO2_INJECTION_TEST_FLASH_HH + +#include + +namespace Opm { +/*! + * \brief Flash solver used by the CO2 injection test problem. + * + * This class is just the NCP flash solver with the guessInitial() + * method that is adapted to the pressure regime of the CO2 injection + * problem. + */ +template +class Co2InjectionTestFlash : public Opm::NcpFlash +{ + using ParentType = Opm::NcpFlash; + + enum { numPhases = FluidSystem::numPhases }; + +public: + /*! + * \brief Guess initial values for all quantities. + */ + template + static void guessInitial(FluidState& fluidState, const ComponentVector& globalMolarities) + { + ParentType::guessInitial(fluidState, globalMolarities); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + // pressure. use something close to the reservoir pressure as initial guess + fluidState.setPressure(phaseIdx, 100e5); + } + } +}; + +} // namespace Opm + +#endif // EWOMS_CO2_INJECTION_TEST_FLASH_HH diff --git a/examples/problems/simpletestproblem.hh b/examples/problems/simpletestproblem.hh new file mode 100644 index 000000000..d55493d01 --- /dev/null +++ b/examples/problems/simpletestproblem.hh @@ -0,0 +1,667 @@ +// -*- 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::simpletestproblem + */ +#ifndef EWOMS_SIMPLETEST_PROBLEM_HH +#define EWOMS_SIMPLETEST_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 SimpleTest; +} // namespace Opm + +namespace Opm::Properties { + +namespace TTag { +struct SimpleTest {}; +} // end namespace TTag + +// declare the "simpletest" problem specify property tags +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: --->1D +template +struct Grid { using type = Dune::YaspGrid; }; + +// Set the problem property +template +struct Problem +{ using type = Opm::SimpleTest; }; + +// 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 saturations + + 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 = 1e5; +}; + +template +struct SimulationName { + static constexpr auto value = "simpletest"; +}; + +// The default for the end time of the simulation +template +struct EndTime { + using type = GetPropType; + static constexpr type value = 3 * 24. * 60. * 60.;//3 days +}; + +// convergence control +template +struct InitialTimeStepSize { + using type = GetPropType; + //static constexpr type value = 30; + static constexpr type value = 1 * 24. * 60. * 60.; +}; + +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 MaxTimeStepSize { + using type = GetPropType; + static constexpr type value = 60 * 60; +}; + +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; +}; + +// write restart for every hour +template +struct EpisodeLength { + using type = GetPropType; + static constexpr type value = 60. * 60.; +}; + +// mesh grid +template +struct Vanguard { + using type = Opm::StructuredGridVanguard; +}; + +//\Note: from the Julia code, the problem is a 1D problem with 3X1 cell. +//\Note: DomainSizeX is 3.0 meters +//\Note: DomainSizeY is 1.0 meters +template +struct DomainSizeX { + using type = GetPropType; + static constexpr type value = 3; // meter +}; + +template +struct DomainSizeY { + using type = GetPropType; + static constexpr type value = 1.0; +}; + +template +struct DomainSizeZ { + using type = GetPropType; + static constexpr type value = 1.0; +}; + +template +struct CellsX { static constexpr int value = 3; }; +template +struct CellsY { static constexpr int value = 1; }; +template +struct CellsZ { static constexpr int value = 1; }; + + +// compositional, with diffusion +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 SimpleTest : 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 + */ + SimpleTest(Simulator& simulator) + : ParentType(simulator) + { + 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.1; + + } + + 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 FvBaseMultiPhaseProblem::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); + std::cout << "primary variables for cell " << context.globalSpaceIndex(spaceIdx, timeIdx) << ": " << values << "\n"; + } + + // 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 = 2; + if (spatialIdx == inj || spatialIdx == prod) + return 1.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 + { + using Scalar = double; + using FluidSystem = Opm::ThreeComponentFluidSystem; + + constexpr auto numComponents = FluidSystem::numComponents; + using Evaluation = Opm::DenseAd::Evaluation; + typedef Dune::FieldVector ComponentVector; + + // input from Olav + //z0 = [0.5, 0.3, 0.2] + //zi = [0.99, 0.01-1e-3, 1e-3] + //p0 = 75e5 + //T0 = 423.25 + int inj = 0; + int prod = 2; + int spatialIdx = context.globalSpaceIndex(spaceIdx, timeIdx); + ComponentVector comp; + comp[0] = Evaluation::createVariable(0.5, 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]; + // TODO: should we put the derivative against the temperature here? + const Scalar temp = 423.25; + + // TODO: no capillary pressure for now + Scalar p0 = 75e5; //CONVERGENCE FAILURE WITH 75 + + //\Note, for an AD variable, if we multiply it with 2, the derivative will also be scalced with 2, + //\Note, so we should not do it. + if (spatialIdx == inj){ + p0 *= 2.0; + } + if (spatialIdx == prod) { + p0 *= 0.5; + } + 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(temp); + + // 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)); + } + + // ComponentVector zInit(0.); // TODO; zInit needs to be normalized. + // { + // Scalar sumMoles = 0.0; + // for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + // for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + // Scalar tmp = Opm::getValue(fs.molarity(phaseIdx, compIdx) * fs.saturation(phaseIdx)); + // zInit[compIdx] += Opm::max(tmp, 1e-8); + // sumMoles += tmp; + // } + // } + // zInit /= sumMoles; + // // initialize the derivatives + // // TODO: the derivative eventually should be from the reservoir flow equations + // Evaluation z_last = 1.; + // for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) { + // zInit[compIdx] = Evaluation::createVariable(Opm::getValue(zInit[compIdx]), compIdx + 1); + // z_last -= zInit[compIdx]; + // } + // zInit[numComponents - 1] = z_last; + // } + + // TODO: only, p, z need the derivatives. + const double flash_tolerance = 1.e-12; // just to test the setup in co2-compositional + //const int flash_verbosity = 1; + const std::string flash_twophase_method = "newton"; // "ssi" + //const std::string flash_twophase_method = "ssi"; + // const std::string flash_twophase_method = "ssi+newton"; + + // 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 \ No newline at end of file diff --git a/opm/models/common/darcyfluxmodule.hh b/opm/models/common/darcyfluxmodule.hh index 097908543..94a466ca0 100644 --- a/opm/models/common/darcyfluxmodule.hh +++ b/opm/models/common/darcyfluxmodule.hh @@ -244,7 +244,9 @@ protected: Evaluation pStatIn; if (std::is_same::value || - interiorDofIdx_ == static_cast(focusDofIdx)) + interiorDofIdx_ == static_cast(focusDofIdx) || + (phaseIdx == 0 && (intQuantsIn.fluidState().L() == 1 && intQuantsEx.fluidState().L() == 0)) || + (phaseIdx == 1 && (intQuantsIn.fluidState().L() == 0 && intQuantsEx.fluidState().L() == 1))) { const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx); pStatIn = - rhoIn*(gIn*distVecIn); @@ -259,7 +261,9 @@ protected: Evaluation pStatEx; if (std::is_same::value || - exteriorDofIdx_ == static_cast(focusDofIdx)) + exteriorDofIdx_ == static_cast(focusDofIdx) || + (phaseIdx == 0 && (intQuantsIn.fluidState().L() == 0 && intQuantsEx.fluidState().L() == 1)) || + (phaseIdx == 1 && (intQuantsIn.fluidState().L() == 1 && intQuantsEx.fluidState().L() == 0))) { const Evaluation& rhoEx = intQuantsEx.fluidState().density(phaseIdx); pStatEx = - rhoEx*(gEx*distVecEx); @@ -274,7 +278,22 @@ protected: // control volume centers and the length (pStaticExterior - // pStaticInterior)/distanceInteriorToExterior Dune::FieldVector f(distVecTotal); - f *= (pStatEx - pStatIn)/absDistTotalSquared; + if (phaseIdx == 0 && (intQuantsIn.fluidState().L() == 1 && intQuantsEx.fluidState().L() == 0)) + f *= -(pStatIn + pStatIn)/absDistTotalSquared; + else if (phaseIdx == 0 && (intQuantsIn.fluidState().L() == 0 && intQuantsEx.fluidState().L() == 1)) + f *= (pStatEx + pStatEx)/absDistTotalSquared; + else if (phaseIdx == 0 && (intQuantsIn.fluidState().L() == 0 && intQuantsEx.fluidState().L() == 0)) + f *= 0.0; + + else if (phaseIdx == 1 && (intQuantsIn.fluidState().L() == 0 && intQuantsEx.fluidState().L() == 1)) + f *= -(pStatIn + pStatIn)/absDistTotalSquared; + else if (phaseIdx == 1 && (intQuantsIn.fluidState().L() == 1 && intQuantsEx.fluidState().L() == 0)) + f *= (pStatEx + pStatEx)/absDistTotalSquared; + else if (phaseIdx == 1 && (intQuantsIn.fluidState().L() == 1 && intQuantsEx.fluidState().L() == 1)) + f *= 0.0; + + else + f *= (pStatEx - pStatIn)/absDistTotalSquared; // calculate the final potential gradient for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) @@ -559,4 +578,4 @@ protected: } // namespace Opm -#endif +#endif \ No newline at end of file diff --git a/opm/models/discretization/common/fvbaseadlocallinearizer.hh b/opm/models/discretization/common/fvbaseadlocallinearizer.hh index 69fb9dfd3..8c205245f 100644 --- a/opm/models/discretization/common/fvbaseadlocallinearizer.hh +++ b/opm/models/discretization/common/fvbaseadlocallinearizer.hh @@ -273,18 +273,26 @@ protected: { const auto& resid = localResidual_.residual(); - for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++) + for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++){ residual_[focusDofIdx][eqIdx] = resid[focusDofIdx][eqIdx].value(); - + std::cout << "residual = " << residual_[focusDofIdx][eqIdx] << std::endl; + } + + int spatialIdx = elemCtx.globalSpaceIndex(focusDofIdx, /*timeIdx=*/0); + std::cout << "Focus cell = " << spatialIdx << std::endl; size_t numDof = elemCtx.numDof(/*timeIdx=*/0); for (unsigned dofIdx = 0; dofIdx < numDof; dofIdx++) { for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++) { for (unsigned pvIdx = 0; pvIdx < numEq; pvIdx++) { + int spatialIdx2 = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); // A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of // the residual function 'eqIdx' for the degree of freedom 'dofIdx' // with regard to the focus variable 'pvIdx' of the degree of freedom // 'focusDofIdx' jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] = resid[dofIdx][eqIdx].derivative(pvIdx); + //std::cout << "J[" << spatialIdx2 << "][" << spatialIdx << "][" << eqIdx << "][" << pvIdx << "] = " << jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] << std::endl; + std::cout << "J[" << dofIdx+1 << "," << focusDofIdx+1 << "][" << eqIdx+1 << "," << pvIdx+1 << "] = " << jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] << std::endl; + Valgrind::CheckDefined(jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx]); } } diff --git a/opm/models/discretization/common/fvbasediscretization.hh b/opm/models/discretization/common/fvbasediscretization.hh index 7c85bedd3..51149086f 100644 --- a/opm/models/discretization/common/fvbasediscretization.hh +++ b/opm/models/discretization/common/fvbasediscretization.hh @@ -274,12 +274,12 @@ struct EnableConstraints { static constexpr // impact because of the intensive quantity cache will cause additional pressure on the // CPU caches... template -struct EnableIntensiveQuantityCache { static constexpr bool value = false; }; +struct EnableIntensiveQuantityCache { static constexpr bool value = true; }; // do not use thermodynamic hints by default. If you enable this, make sure to also // enable the intensive quantity cache above to avoid getting an exception... template -struct EnableThermodynamicHints { static constexpr bool value = false; }; +struct EnableThermodynamicHints { static constexpr bool value = true; }; // if the deflection of the newton method is large, we do not need to solve the linear // approximation accurately. Assuming that the value for the current solution is quite @@ -309,9 +309,9 @@ struct TimeDiscHistorySize { static constex template struct ExtensiveStorageTerm { static constexpr bool value = false; }; -// use volumetric residuals is default +// use volumetric residuals is not default template -struct UseVolumetricResidual { static constexpr bool value = true; }; +struct UseVolumetricResidual { static constexpr bool value = false; }; //! eWoms is mainly targeted at research, so experimental features are enabled by //! default. diff --git a/opm/models/discretization/common/fvbaseelementcontext.hh b/opm/models/discretization/common/fvbaseelementcontext.hh index 55e1cdc51..bfaa9cf78 100644 --- a/opm/models/discretization/common/fvbaseelementcontext.hh +++ b/opm/models/discretization/common/fvbaseelementcontext.hh @@ -554,8 +554,8 @@ protected: const PrimaryVariables& dofSol = globalSol[globalIdx]; dofVars_[dofIdx].priVars[timeIdx] = &dofSol; - dofVars_[dofIdx].thermodynamicHint[timeIdx] = - model().thermodynamicHint(globalIdx, timeIdx); + dofVars_[dofIdx].thermodynamicHint[timeIdx] = model().thermodynamicHint(globalIdx, timeIdx); + dofVars_[dofIdx].thermodynamicHint[1] = model().thermodynamicHint(globalIdx, 1); const auto *cachedIntQuants = model().cachedIntensiveQuantities(globalIdx, timeIdx); if (cachedIntQuants) { diff --git a/opm/models/discretization/common/fvbasefdlocallinearizer.hh b/opm/models/discretization/common/fvbasefdlocallinearizer.hh index 27082d7e7..7111bd44a 100644 --- a/opm/models/discretization/common/fvbasefdlocallinearizer.hh +++ b/opm/models/discretization/common/fvbasefdlocallinearizer.hh @@ -490,6 +490,8 @@ protected: unsigned pvIdx) { size_t numDof = elemCtx.numDof(/*timeIdx=*/0); + int spatialIdx = elemCtx.globalSpaceIndex(focusDofIdx, /*timeIdx=*/0); + std::cout << "Focus cell = " << spatialIdx << std::endl; for (unsigned dofIdx = 0; dofIdx < numDof; dofIdx++) { for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++) { // A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of the @@ -497,6 +499,7 @@ protected: // regard to the primary variable 'pvIdx' of the degree of freedom // 'focusDofIdx' jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] = derivResidual_[dofIdx][eqIdx]; + std::cout << "J[" << dofIdx << "][" << focusDofIdx << "][" << eqIdx << "][" << pvIdx << "] = " << jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] << std::endl; Valgrind::CheckDefined(jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx]); } } diff --git a/opm/models/discretization/common/fvbaselocalresidual.hh b/opm/models/discretization/common/fvbaselocalresidual.hh index e8c4d0716..f852210b0 100644 --- a/opm/models/discretization/common/fvbaselocalresidual.hh +++ b/opm/models/discretization/common/fvbaselocalresidual.hh @@ -161,6 +161,8 @@ public: assert(residual.size() == elemCtx.numDof(/*timeIdx=*/0)); residual = 0.0; + std::cout << "residual from elemtCtx" << std::endl; + // evaluate the flux terms asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0); @@ -324,8 +326,10 @@ public: assert(alpha > 0.0); assert(isfinite(alpha)); - for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) + for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx){ flux[eqIdx] *= alpha; + std::cout << " flux in eqIdx " << eqIdx << " is: " << flux[eqIdx] << std::endl; + } // The balance equation for a finite volume is given by // diff --git a/opm/models/io/vtkcompositionmodule.hh b/opm/models/io/vtkcompositionmodule.hh index d307b19f7..83cfb46c2 100644 --- a/opm/models/io/vtkcompositionmodule.hh +++ b/opm/models/io/vtkcompositionmodule.hh @@ -59,6 +59,10 @@ template struct VtkWriteFugacities { using type = UndefinedProperty; }; template struct VtkWriteFugacityCoeffs { using type = UndefinedProperty; }; +template +struct VtkWriteLiquidMoleFractions { using type = UndefinedProperty; }; +template +struct VtkWriteEquilibriumConstants { using type = UndefinedProperty; }; // set default values for what quantities to output template @@ -75,6 +79,10 @@ template struct VtkWriteFugacities { static constexpr bool value = false; }; template struct VtkWriteFugacityCoeffs { static constexpr bool value = false; }; +template +struct VtkWriteLiquidMoleFractions { static constexpr bool value = false; }; +template +struct VtkWriteEquilibriumConstants { static constexpr bool value = false; }; } // namespace Opm::Properties @@ -112,6 +120,7 @@ class VtkCompositionModule : public BaseOutputModule using ComponentBuffer = typename ParentType::ComponentBuffer; using PhaseComponentBuffer = typename ParentType::PhaseComponentBuffer; + using ScalarBuffer = typename ParentType::ScalarBuffer; public: VtkCompositionModule(const Simulator& simulator) @@ -137,6 +146,10 @@ public: "Include component fugacities in the VTK output files"); EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteFugacityCoeffs, "Include component fugacity coefficients in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteLiquidMoleFractions, + "Include liquid mole fractions (L) in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteEquilibriumConstants, + "Include equilibrium constants (K) in the VTK output files"); } /*! @@ -155,7 +168,10 @@ public: this->resizeComponentBuffer_(totalMoleFrac_); if (molarityOutput_()) this->resizePhaseComponentBuffer_(molarity_); - + if (LOutput_()) + this->resizeScalarBuffer_(L_); + if (equilConstOutput_()) + this->resizeComponentBuffer_(K_); if (fugacityOutput_()) this->resizeComponentBuffer_(fugacity_); if (fugacityCoeffOutput_()) @@ -178,6 +194,9 @@ public: const auto& intQuants = elemCtx.intensiveQuantities(i, /*timeIdx=*/0); const auto& fs = intQuants.fluidState(); + if (LOutput_()) + L_[I] = Toolbox::value(fs.L()); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { if (moleFracOutput_()) @@ -222,6 +241,10 @@ public: } if (fugacityOutput_()) fugacity_[compIdx][I] = Toolbox::value(intQuants.fluidState().fugacity(/*phaseIdx=*/0, compIdx)); + + if (equilConstOutput_()) + K_[compIdx][I] = Toolbox::value(fs.K(compIdx)); + } } } @@ -236,7 +259,7 @@ public: return; } - if (moleFracOutput_()) + if (moleFracOutput_()) this->commitPhaseComponentBuffer_(baseWriter, "moleFrac_%s^%s", moleFrac_); if (massFracOutput_()) this->commitPhaseComponentBuffer_(baseWriter, "massFrac_%s^%s", massFrac_); @@ -246,11 +269,15 @@ public: this->commitComponentBuffer_(baseWriter, "totalMassFrac^%s", totalMassFrac_); if (totalMoleFracOutput_()) this->commitComponentBuffer_(baseWriter, "totalMoleFrac^%s", totalMoleFrac_); + if (equilConstOutput_()) + this->commitComponentBuffer_(baseWriter, "K^%s", K_); if (fugacityOutput_()) this->commitComponentBuffer_(baseWriter, "fugacity^%s", fugacity_); if (fugacityCoeffOutput_()) this->commitPhaseComponentBuffer_(baseWriter, "fugacityCoeff_%s^%s", fugacityCoeff_); + if (LOutput_()) + this->commitScalarBuffer_(baseWriter, "L", L_); } private: @@ -295,15 +322,30 @@ private: static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteFugacityCoeffs); return val; } + + static bool LOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteLiquidMoleFractions); + return val; + } + + static bool equilConstOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteEquilibriumConstants); + return val; + } PhaseComponentBuffer moleFrac_; PhaseComponentBuffer massFrac_; PhaseComponentBuffer molarity_; ComponentBuffer totalMassFrac_; ComponentBuffer totalMoleFrac_; + ComponentBuffer K_; ComponentBuffer fugacity_; PhaseComponentBuffer fugacityCoeff_; + + ScalarBuffer L_; }; } // namespace Opm diff --git a/opm/models/nonlinear/newtonmethod.hh b/opm/models/nonlinear/newtonmethod.hh index df85b45e9..8e6d22ef1 100644 --- a/opm/models/nonlinear/newtonmethod.hh +++ b/opm/models/nonlinear/newtonmethod.hh @@ -44,6 +44,7 @@ #include #include #include +//#include #include #include @@ -128,7 +129,7 @@ class NewtonMethod using Communicator = typename Dune::MPIHelper::MPICommunicator; using CollectiveCommunication = typename Dune::Communication; - + // using PrintMatrix = typename Dune::printSparseMatrix; public: NewtonMethod(Simulator& simulator) : simulator_(simulator) @@ -319,6 +320,7 @@ public: linearSolver_.getResidual(residual); solveTimer_.stop(); + //Dune::storeMatrixMarket(linearSolver_.overlappingMatrix_, "mymatrix" + ".mm"); // The preSolve_() method usually computes the errors, but it can do // something else in addition. TODO: should its costs be counted to // the linearization or to the update? diff --git a/opm/models/ptflash/flashboundaryratevector.hh b/opm/models/ptflash/flashboundaryratevector.hh new file mode 100644 index 000000000..c932dc995 --- /dev/null +++ b/opm/models/ptflash/flashboundaryratevector.hh @@ -0,0 +1,205 @@ +// -*- 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::FlashBoundaryRateVector + */ +#ifndef EWOMS_FLASH_BOUNDARY_RATE_VECTOR_HH +#define EWOMS_FLASH_BOUNDARY_RATE_VECTOR_HH + +#include "flashproperties.hh" + +#include +#include + +namespace Opm { + +/*! + * \ingroup FlashModel + * + * \brief Implements a boundary vector for the fully implicit + * compositional multi-phase model which is based on flash + * calculations. + */ +template +class FlashBoundaryRateVector : public GetPropType +{ + using ParentType = GetPropType; + using ExtensiveQuantities = GetPropType; + using FluidSystem = GetPropType; + using Scalar = GetPropType; + using Evaluation = GetPropType; + using Indices = GetPropType; + + enum { numEq = getPropValue() }; + enum { numPhases = getPropValue() }; + enum { numComponents = getPropValue() }; + enum { conti0EqIdx = Indices::conti0EqIdx }; + enum { enableEnergy = getPropValue() }; + + using EnergyModule = Opm::EnergyModule; + using Toolbox = Opm::MathToolbox; + +public: + FlashBoundaryRateVector() : ParentType() + {} + + /*! + * \copydoc + * ImmiscibleBoundaryRateVector::ImmiscibleBoundaryRateVector(Scalar) + */ + FlashBoundaryRateVector(const Evaluation& value) : ParentType(value) + {} + + /*! + * \copydoc ImmiscibleBoundaryRateVector::ImmiscibleBoundaryRateVector(const + * ImmiscibleBoundaryRateVector& ) + */ + FlashBoundaryRateVector(const FlashBoundaryRateVector& value) = default; + FlashBoundaryRateVector& operator=(const FlashBoundaryRateVector& value) = default; + + /*! + * \copydoc ImmiscibleBoundaryRateVector::setFreeFlow + */ + template + void setFreeFlow(const Context& context, + unsigned bfIdx, + unsigned timeIdx, + const FluidState& fluidState) + { + ExtensiveQuantities extQuants; + extQuants.updateBoundary(context, bfIdx, timeIdx, fluidState); + const auto& insideIntQuants = context.intensiveQuantities(bfIdx, timeIdx); + unsigned focusDofIdx = context.focusDofIndex(); + unsigned interiorDofIdx = context.interiorScvIndex(bfIdx, timeIdx); + + //////// + // advective fluxes of all components in all phases + //////// + (*this) = Evaluation(0.0); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + Evaluation density; + if (fluidState.pressure(phaseIdx) > insideIntQuants.fluidState().pressure(phaseIdx)) { + if (focusDofIdx == interiorDofIdx) + density = fluidState.density(phaseIdx); + else + density = Opm::getValue(fluidState.density(phaseIdx)); + } + else if (focusDofIdx == interiorDofIdx) + density = insideIntQuants.fluidState().density(phaseIdx); + else + density = Opm::getValue(insideIntQuants.fluidState().density(phaseIdx)); + + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + Evaluation molarity; + if (fluidState.pressure(phaseIdx) > insideIntQuants.fluidState().pressure(phaseIdx)) { + if (focusDofIdx == interiorDofIdx) + molarity = fluidState.molarity(phaseIdx, compIdx); + else + molarity = Opm::getValue(fluidState.molarity(phaseIdx, compIdx)); + } + else if (focusDofIdx == interiorDofIdx) + molarity = insideIntQuants.fluidState().molarity(phaseIdx, compIdx); + else + molarity = Opm::getValue(insideIntQuants.fluidState().molarity(phaseIdx, compIdx)); + + // add advective flux of current component in current + // phase + (*this)[conti0EqIdx + compIdx] += extQuants.volumeFlux(phaseIdx)*molarity; + } + + if (enableEnergy) { + Evaluation specificEnthalpy; + if (fluidState.pressure(phaseIdx) > insideIntQuants.fluidState().pressure(phaseIdx)) { + if (focusDofIdx == interiorDofIdx) + specificEnthalpy = fluidState.enthalpy(phaseIdx); + else + specificEnthalpy = Opm::getValue(fluidState.enthalpy(phaseIdx)); + } + else if (focusDofIdx == interiorDofIdx) + specificEnthalpy = insideIntQuants.fluidState().enthalpy(phaseIdx); + else + specificEnthalpy = Opm::getValue(insideIntQuants.fluidState().enthalpy(phaseIdx)); + + Evaluation enthalpyRate = density*extQuants.volumeFlux(phaseIdx)*specificEnthalpy; + EnergyModule::addToEnthalpyRate(*this, enthalpyRate); + } + } + + // thermal conduction + EnergyModule::addToEnthalpyRate(*this, EnergyModule::thermalConductionRate(extQuants)); + +#ifndef NDEBUG + for (unsigned i = 0; i < numEq; ++i) { + Opm::Valgrind::CheckDefined((*this)[i]); + } +#endif + } + + /*! + * \copydoc ImmiscibleBoundaryRateVector::setInFlow + */ + template + void setInFlow(const Context& context, + unsigned bfIdx, + unsigned timeIdx, + const FluidState& fluidState) + { + this->setFreeFlow(context, bfIdx, timeIdx, fluidState); + + // we only allow fluxes in the direction opposite to the outer unit normal + for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { + Evaluation& val = this->operator[](eqIdx); + val = Toolbox::min(0.0, val); + } + } + + /*! + * \copydoc ImmiscibleBoundaryRateVector::setOutFlow + */ + template + void setOutFlow(const Context& context, + unsigned bfIdx, + unsigned timeIdx, + const FluidState& fluidState) + { + this->setFreeFlow(context, bfIdx, timeIdx, fluidState); + + // we only allow fluxes in the same direction as the outer unit normal + for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { + Evaluation& val = this->operator[](eqIdx); + val = Toolbox::max(0.0, val); + } + } + + /*! + * \copydoc ImmiscibleBoundaryRateVector::setNoFlow + */ + void setNoFlow() + { (*this) = Evaluation(0.0); } +}; + +} // namespace Opm + +#endif diff --git a/opm/models/ptflash/flashextensivequantities.hh b/opm/models/ptflash/flashextensivequantities.hh new file mode 100644 index 000000000..344e8c66d --- /dev/null +++ b/opm/models/ptflash/flashextensivequantities.hh @@ -0,0 +1,96 @@ +// -*- 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::FlashExtensiveQuantities + */ +#ifndef EWOMS_FLASH_EXTENSIVE_QUANTITIES_HH +#define EWOMS_FLASH_EXTENSIVE_QUANTITIES_HH + +#include "flashproperties.hh" + +#include +#include +#include + +namespace Opm { + +/*! + * \ingroup FlashModel + * \ingroup ExtensiveQuantities + * + * \brief This template class contains the data which is required to + * calculate all fluxes of components over a face of a finite + * volume for the compositional multi-phase model based on + * flash calculations. + * + * This means pressure and concentration gradients, phase densities at + * the integration point, etc. + */ +template +class FlashExtensiveQuantities + : public MultiPhaseBaseExtensiveQuantities + , public EnergyExtensiveQuantities()> + , public DiffusionExtensiveQuantities()> +{ + using ParentType = MultiPhaseBaseExtensiveQuantities; + + using ElementContext = GetPropType; + using FluidSystem = GetPropType; + + enum { enableDiffusion = getPropValue() }; + using DiffusionExtensiveQuantities = Opm::DiffusionExtensiveQuantities; + + enum { enableEnergy = getPropValue() }; + using EnergyExtensiveQuantities = Opm::EnergyExtensiveQuantities; + +public: + /*! + * \copydoc MultiPhaseBaseExtensiveQuantities::update + */ + void update(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx) + { + ParentType::update(elemCtx, scvfIdx, timeIdx); + DiffusionExtensiveQuantities::update_(elemCtx, scvfIdx, timeIdx); + EnergyExtensiveQuantities::update_(elemCtx, scvfIdx, timeIdx); + } + + /*! + * \copydoc MultiPhaseBaseExtensiveQuantities::updateBoundary + */ + template + void updateBoundary(const Context& context, + unsigned bfIdx, + unsigned timeIdx, + const FluidState& fluidState) + { + ParentType::updateBoundary(context, bfIdx, timeIdx, fluidState); + DiffusionExtensiveQuantities::updateBoundary_(context, bfIdx, timeIdx, fluidState); + EnergyExtensiveQuantities::updateBoundary_(context, bfIdx, timeIdx, fluidState); + } +}; + +} // namespace Opm + +#endif diff --git a/opm/models/ptflash/flashindices.hh b/opm/models/ptflash/flashindices.hh new file mode 100644 index 000000000..e2cdefc1a --- /dev/null +++ b/opm/models/ptflash/flashindices.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: +/* + 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::FlashIndices + */ +#ifndef EWOMS_FLASH_INDICES_HH +#define EWOMS_FLASH_INDICES_HH + +#include "flashproperties.hh" +#include + +namespace Opm { + +/*! + * \ingroup FlashModel + * + * \brief Defines the primary variable and equation indices for the + * compositional multi-phase model based on flash calculations. + * + * \tparam PVOffset The first index in a primary variable vector. + */ +template +class FlashIndices + : public EnergyIndices(), + getPropValue()> +{ + enum { numComponents = getPropValue() }; + enum { enableEnergy = getPropValue() }; + using EnergyIndices = Opm::EnergyIndices; + +public: + //! number of equations/primary variables + static const int numEq = numComponents + EnergyIndices::numEq_; + + // Primary variable indices + + //! Index of the total concentration of the first component in the + //! pore space. + static const int pressure0Idx = PVOffset; + static const int z0Idx = pressure0Idx+1; // numcomponents-1 variables follow + + // equation indices + + //! Index of the mass conservation equation for the first + //! component. + static const int conti0EqIdx = PVOffset; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/ptflash/flashintensivequantities.hh b/opm/models/ptflash/flashintensivequantities.hh new file mode 100644 index 000000000..7a4562001 --- /dev/null +++ b/opm/models/ptflash/flashintensivequantities.hh @@ -0,0 +1,331 @@ +// -*- 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::FlashIntensiveQuantities + */ +#ifndef EWOMS_FLASH_INTENSIVE_QUANTITIES_HH +#define EWOMS_FLASH_INTENSIVE_QUANTITIES_HH + +#include "flashproperties.hh" +#include "flashindices.hh" + +#include +#include + +#include +#include + +#include +#include + +namespace Opm { + +/*! + * \ingroup FlashModel + * \ingroup IntensiveQuantities + * + * \brief Contains the intensive quantities of the flash-based compositional multi-phase model + */ +template +class FlashIntensiveQuantities + : public GetPropType + , public DiffusionIntensiveQuantities() > + , public EnergyIntensiveQuantities() > + , public GetPropType::FluxIntensiveQuantities +{ + using ParentType = GetPropType; + + using ElementContext = GetPropType; + using MaterialLaw = GetPropType; + using MaterialLawParams = GetPropType; + using Indices = GetPropType; + using FluxModule = GetPropType; + using GridView = GetPropType; + using ThreadManager = GetPropType; + + // primary variable indices + enum { z0Idx = Indices::z0Idx }; + enum { numPhases = getPropValue() }; + enum { numComponents = getPropValue() }; + enum { enableDiffusion = getPropValue() }; + enum { enableEnergy = getPropValue() }; + enum { dimWorld = GridView::dimensionworld }; + enum { pressure0Idx = Indices::pressure0Idx }; + + using Scalar = GetPropType; + using Evaluation = GetPropType; + using FluidSystem = GetPropType; + using FlashSolver = GetPropType; + using Problem = GetPropType; + using ComponentVector = Dune::FieldVector; + using DimMatrix = Dune::FieldMatrix; + + using FluxIntensiveQuantities = typename FluxModule::FluxIntensiveQuantities; + using DiffusionIntensiveQuantities = Opm::DiffusionIntensiveQuantities; + using EnergyIntensiveQuantities = Opm::EnergyIntensiveQuantities; + +public: + //! The type of the object returned by the fluidState() method + using FluidState = Opm::CompositionalFluidState; + + FlashIntensiveQuantities() + { } + + FlashIntensiveQuantities(const FlashIntensiveQuantities& other) = default; + + FlashIntensiveQuantities& operator=(const FlashIntensiveQuantities& other) = default; + + /*! + * \copydoc IntensiveQuantities::update + */ + void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx) + { + ParentType::update(elemCtx, dofIdx, timeIdx); + EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx); + + const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx); + const auto& problem = elemCtx.problem(); + + Scalar flashTolerance = 1.e-12;//EWOMS_GET_PARAM(TypeTag, Scalar, FlashTolerance); + int flashVerbosity = 0;//EWOMS_GET_PARAM(TypeTag, int, FlashVerbosity); + std::string flashTwoPhaseMethod = EWOMS_GET_PARAM(TypeTag, std::string, FlashTwoPhaseMethod); + + // extract the total molar densities of the components + ComponentVector z(0.); + { + Evaluation lastZ = 1.0; + for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) { + z[compIdx] = priVars.makeEvaluation(z0Idx + compIdx, timeIdx); + lastZ -= z[compIdx]; + } + z[numComponents - 1] = lastZ; + + Evaluation sumz = 0.0; + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + z[compIdx] = Opm::max(z[compIdx], 1e-8); + sumz +=z[compIdx]; + } + z /= sumz; + + } + + Evaluation p = priVars.makeEvaluation(pressure0Idx, timeIdx); + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + fluidState_.setPressure(phaseIdx, p); + + // Get initial K and L from storage initially (if enabled) + const auto *hint = elemCtx.thermodynamicHint(dofIdx, timeIdx); + const auto *hint2 = elemCtx.thermodynamicHint(dofIdx, 1); + if (hint) { + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + const Evaluation& Ktmp = hint->fluidState().K(compIdx); + fluidState_.setKvalue(compIdx, Ktmp); + } + const Evaluation& Ltmp = hint->fluidState().L(); + fluidState_.setLvalue(Ltmp); + } + else if (hint2) { + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + const Evaluation& Ktmp = hint2->fluidState().K(compIdx); + fluidState_.setKvalue(compIdx, Ktmp); + } + const Evaluation& Ltmp = hint2->fluidState().L(); + fluidState_.setLvalue(Ltmp); + } + else { + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + const Evaluation Ktmp = fluidState_.wilsonK_(compIdx); + fluidState_.setKvalue(compIdx, Ktmp); + } + const Evaluation& Ltmp = -1.0; + fluidState_.setLvalue(Ltmp); + } + + ///////////// + // Compute the phase compositions and densities + ///////////// + int spatialIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx); + //FlashSolver::solve(fluidState_, z, spatialIdx, flashVerbosity, flashTwoPhaseMethod, flashTolerance); + //Flash::solve(fluidState_, z, spatialIdx, flashVerbosity, flashTwoPhaseMethod, flashTolerance); + + + using Flash = Opm::PTFlash; + FlashSolver::solve(fluidState_, z, spatialIdx, flashTwoPhaseMethod, flashTolerance, flashVerbosity); + +//printing of flashresult after solve + // std::cout << " \n After flashsolve for cell " << spatialIdx << std::endl; + // ComponentVector x, y; + // Evaluation L0 = fluidState_.L(); + // for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) { + // x[comp_idx] = fluidState_.moleFraction(FluidSystem::oilPhaseIdx, comp_idx); + // y[comp_idx] = fluidState_.moleFraction(FluidSystem::gasPhaseIdx, comp_idx); + // } + // for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) { + // std::cout << " x for component: " << comp_idx << " is " << x[comp_idx] << std::endl; + // for (int i = 0; i < 3; ++i) { + // std::cout << " x deriv " << i << " is: " << x[comp_idx].derivative(i) << std::endl; + // } + + // std::cout << " y for component: " << comp_idx << "is " << y[comp_idx] << std::endl; + // for (int i = 0; i < 3; ++i) { + // std::cout << " y deriv " << i << " is: " << y[comp_idx].derivative(i) << std::endl; + // } + // } + // std::cout << " L is " << L0 << std::endl; + // for (int i = 0; i < L0.size(); ++i) { + // std::cout << " L deriv " << i << " is: " << L0.derivative(i) << std::endl; + // } + // //end printinting 1 + + // Update phases + typename FluidSystem::template ParameterCache paramCache; + paramCache.updatePhase(fluidState_, FluidSystem::oilPhaseIdx); + + const Scalar R = Opm::Constants::R; + Evaluation Z_L = (paramCache.molarVolume(FluidSystem::oilPhaseIdx) * fluidState_.pressure(FluidSystem::oilPhaseIdx) )/ + (R * fluidState_.temperature(FluidSystem::oilPhaseIdx)); + paramCache.updatePhase(fluidState_, FluidSystem::gasPhaseIdx); + Evaluation Z_V = (paramCache.molarVolume(FluidSystem::gasPhaseIdx) * fluidState_.pressure(FluidSystem::gasPhaseIdx) )/ + (R * fluidState_.temperature(FluidSystem::gasPhaseIdx)); + + + // Update saturation + Evaluation L = fluidState_.L(); + Evaluation So = Opm::max((L*Z_L/(L*Z_L+(1-L)*Z_V)), 0.0); + Evaluation Sg = Opm::max(1-So, 0.0); + Scalar sumS = Opm::getValue(So) + Opm::getValue(Sg); + So /= sumS; + Sg /= sumS; + + fluidState_.setSaturation(0, So); + fluidState_.setSaturation(1, Sg); + + fluidState_.setCompressFactor(0, Z_L); + fluidState_.setCompressFactor(1, Z_V); + + // Print saturation + if (flashVerbosity == 5) { + std::cout << "So = " << So < +#include + +#include + +namespace Opm { +/*! + * \ingroup FlashModel + * + * \brief Calculates the local residual of the compositional multi-phase + * model based on flash calculations. + */ +template +class FlashLocalResidual: public GetPropType +{ + using Evaluation = GetPropType; + using EqVector = GetPropType; + using RateVector = GetPropType; + using Indices = GetPropType; + using IntensiveQuantities = GetPropType; + using ElementContext = GetPropType; + + enum { numEq = getPropValue() }; + enum { numPhases = getPropValue() }; + enum { numComponents = getPropValue() }; + enum { conti0EqIdx = Indices::conti0EqIdx }; + + enum { enableDiffusion = getPropValue() }; + using DiffusionModule = Opm::DiffusionModule; + + enum { enableEnergy = getPropValue() }; + using EnergyModule = Opm::EnergyModule; + + using Toolbox = Opm::MathToolbox; + +public: + /*! + * \copydoc ImmiscibleLocalResidual::addPhaseStorage + */ + template + void addPhaseStorage(Dune::FieldVector& storage, + const ElementContext& elemCtx, + unsigned dofIdx, + unsigned timeIdx, + unsigned phaseIdx) const + { + const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx); + const auto& fs = intQuants.fluidState(); + + // compute storage term of all components within all phases + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + unsigned eqIdx = conti0EqIdx + compIdx; + storage[eqIdx] += + Toolbox::template decay(fs.massFraction(phaseIdx, compIdx)) + * Toolbox::template decay(fs.density(phaseIdx)) + * Toolbox::template decay(fs.saturation(phaseIdx)) + * Toolbox::template decay(intQuants.porosity()); + } + + EnergyModule::addPhaseStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx), phaseIdx); + } + + /*! + * \copydoc FvBaseLocalResidual::computeStorage + */ + template + void computeStorage(Dune::FieldVector& storage, + const ElementContext& elemCtx, + unsigned dofIdx, + unsigned timeIdx) const + { + storage = 0; + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + addPhaseStorage(storage, elemCtx, dofIdx, timeIdx, phaseIdx); + + EnergyModule::addSolidEnergyStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx)); + } + + /*! + * \copydoc FvBaseLocalResidual::computeFlux + */ + void computeFlux(RateVector& flux, + const ElementContext& elemCtx, + unsigned scvfIdx, + unsigned timeIdx) const + { + flux = 0.0; + addAdvectiveFlux(flux, elemCtx, scvfIdx, timeIdx); + Opm::Valgrind::CheckDefined(flux); + + addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx); + Opm::Valgrind::CheckDefined(flux); + } + + /*! + * \copydoc ImmiscibleLocalResidual::addAdvectiveFlux + */ + void addAdvectiveFlux(RateVector& flux, + const ElementContext& elemCtx, + unsigned scvfIdx, + unsigned timeIdx) const + { + const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx); + + unsigned focusDofIdx = elemCtx.focusDofIndex(); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + // data attached to upstream and the finite volume of the current phase + unsigned upIdx = static_cast(extQuants.upstreamIndex(phaseIdx)); + const IntensiveQuantities& up = elemCtx.intensiveQuantities(upIdx, timeIdx); + + // this is a bit hacky because it is specific to the element-centered + // finite volume scheme. (N.B. that if finite differences are used to + // linearize the system of equations, it does not matter.) + if (upIdx == focusDofIdx) { + Evaluation tmp = + up.fluidState().density(phaseIdx) + * extQuants.volumeFlux(phaseIdx); + + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + flux[conti0EqIdx + compIdx] += + tmp*up.fluidState().massFraction(phaseIdx, compIdx); + } + } + else { + Evaluation tmp = + Toolbox::value(up.fluidState().density(phaseIdx) + * extQuants.volumeFlux(phaseIdx)); + + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + flux[conti0EqIdx + compIdx] += + tmp*Toolbox::value(up.fluidState().massFraction(phaseIdx, compIdx)); + } + } + } + + EnergyModule::addAdvectiveFlux(flux, elemCtx, scvfIdx, timeIdx); + } + + /*! + * \copydoc ImmiscibleLocalResidual::addDiffusiveFlux + */ + void addDiffusiveFlux(RateVector& flux, + const ElementContext& elemCtx, + unsigned scvfIdx, + unsigned timeIdx) const + { + DiffusionModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx); + EnergyModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx); + } + + /*! + * \copydoc ImmiscibleLocalResidual::computeSource + */ + void computeSource(RateVector& source, + const ElementContext& elemCtx, + unsigned dofIdx, + unsigned timeIdx) const + { + Opm::Valgrind::SetUndefined(source); + elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx); + Opm::Valgrind::CheckDefined(source); + } +}; + +} // namespace Opm + +#endif diff --git a/opm/models/ptflash/flashmodel.hh b/opm/models/ptflash/flashmodel.hh new file mode 100644 index 000000000..25b7a783d --- /dev/null +++ b/opm/models/ptflash/flashmodel.hh @@ -0,0 +1,355 @@ +// -*- 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::FlashModel + */ +#ifndef EWOMS_FLASH_MODEL_HH +#define EWOMS_FLASH_MODEL_HH + +#include + +#include "flashproperties.hh" +#include "flashprimaryvariables.hh" +#include "flashlocalresidual.hh" +#include "flashratevector.hh" +#include "flashboundaryratevector.hh" +#include "flashintensivequantities.hh" +#include "flashextensivequantities.hh" +#include "flashindices.hh" +#include "flashnewtonmethod.hh" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace Opm { +template +class FlashModel; +} + +namespace Opm::Properties { + +namespace TTag { +//! The type tag for the isothermal single phase problems +struct FlashModel { using InheritsFrom = std::tuple; }; +} // namespace TTag + +//! Use the FlashLocalResidual function for the flash model +template +struct LocalResidual { using type = Opm::FlashLocalResidual; }; + +//! Use the PT flash specific newton method for the flash model +template +struct NewtonMethod { using type = Opm::FlashNewtonMethod; }; + +//! Use the Pt flash solver by default +template +struct FlashSolver +{ using type = Opm::PTFlash, + GetPropType>; }; + +//! Let the flash solver choose its tolerance by default +template +struct FlashTolerance +{ + using type = GetPropType; + static constexpr type value = -1.0; +}; + +// Flash solver verbosity +template +struct FlashVerbosity { static constexpr int value = 0; }; + +// Flash two-phase method +template +struct FlashTwoPhaseMethod { static constexpr auto value = "ssi"; }; + +//! the Model property +template +struct Model { using type = Opm::FlashModel; }; + +//! the PrimaryVariables property +template +struct PrimaryVariables { using type = Opm::FlashPrimaryVariables; }; + +//! the RateVector property +template +struct RateVector { using type = Opm::FlashRateVector; }; + +//! the BoundaryRateVector property +template +struct BoundaryRateVector { using type = Opm::FlashBoundaryRateVector; }; + +//! the IntensiveQuantities property +template +struct IntensiveQuantities { using type = Opm::FlashIntensiveQuantities; }; + +//! the ExtensiveQuantities property +template +struct ExtensiveQuantities { using type = Opm::FlashExtensiveQuantities; }; + +//! The indices required by the flash-baseed isothermal compositional model +template +struct Indices { using type = Opm::FlashIndices; }; + +// The updates of intensive quantities tend to be _very_ expensive for this +// model, so let's try to minimize the number of required ones +template +struct EnableIntensiveQuantityCache { static constexpr bool value = true; }; + +// since thermodynamic hints are basically free if the cache for intensive quantities is +// enabled, and this model usually shows quite a performance improvment if they are +// enabled, let's enable them by default. +template +struct EnableThermodynamicHints { static constexpr bool value = true; }; + +// disable molecular diffusion by default +template +struct EnableDiffusion { static constexpr bool value = false; }; + +//! Disable the energy equation by default +template +struct EnableEnergy { static constexpr bool value = false; }; + +} // namespace Opm::Properties + +namespace Opm { + +/*! + * \ingroup FlashModel + * + * \brief A compositional multi-phase model based on flash-calculations + * + * This model assumes a flow of \f$M \geq 1\f$ fluid phases + * \f$\alpha\f$, each of which is assumed to be a mixture \f$N \geq + * M\f$ chemical species (denoted by the upper index \f$\kappa\f$). + * + * By default, the standard multi-phase Darcy approach is used to determine + * the velocity, i.e. + * \f[ + * \mathbf{v}_\alpha = + * - \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K} + * \left(\mathbf{grad}\, p_\alpha + * - \varrho_{\alpha} \mathbf{g} \right) \;, + * \f] + * although the actual approach which is used can be specified via the + * \c FluxModule property. For example, the velocity model can by + * changed to the Forchheimer approach by + * \code + * template +struct FluxModule { using type = Opm::ForchheimerFluxModule; }; + * \endcode + * + * The core of the model is the conservation mass of each component by + * means of the equation + * \f[ + * \sum_\alpha \frac{\partial\;\phi c_\alpha^\kappa S_\alpha }{\partial t} + * - \sum_\alpha \mathrm{div} \left\{ c_\alpha^\kappa \mathbf{v}_\alpha \right\} + * - q^\kappa = 0 \;. + * \f] + * + * To determine the quanties that occur in the equations above, this + * model uses flash calculations. A flash solver starts with + * the total mass or molar mass per volume for each component and, + * calculates the compositions, saturations and pressures of all + * phases at a given temperature. For this the flash solver has to use + * some model assumptions internally. (Often these are the same + * primary variable switching or NCP assumptions as used by the other + * fully implicit compositional multi-phase models provided by eWoms.) + * + * Using flash calculations for the flow model has some disadvantages: + * - The accuracy of the flash solver needs to be sufficient to + * calculate the parital derivatives using numerical differentiation + * which are required for the Newton scheme. + * - Flash calculations tend to be quite computationally expensive and + * are often numerically unstable. + * + * It is thus adviced to increase the target tolerance of the Newton + * scheme or a to use type for scalar values which exhibits higher + * precision than the standard \c double (e.g. \c quad) if this model + * ought to be used. + * + * The model uses the following primary variables: + * - The total molar concentration of each component: + * \f$c^\kappa = \sum_\alpha S_\alpha x_\alpha^\kappa \rho_{mol, \alpha}\f$ + * - The absolute temperature $T$ in Kelvins if the energy equation enabled. + */ +template +class FlashModel + : public MultiPhaseBaseModel +{ + using ParentType = MultiPhaseBaseModel; + + using Scalar = GetPropType; + using FluidSystem = GetPropType; + using Simulator = GetPropType; + + using Indices = GetPropType; + + enum { numComponents = getPropValue() }; + enum { enableDiffusion = getPropValue() }; + enum { enableEnergy = getPropValue() }; + + + using EnergyModule = Opm::EnergyModule; + +public: + FlashModel(Simulator& simulator) + : ParentType(simulator) + {} + + /*! + * \brief Register all run-time parameters for the immiscible model. + */ + static void registerParameters() + { + ParentType::registerParameters(); + + // register runtime parameters of the VTK output modules + Opm::VtkCompositionModule::registerParameters(); + + if (enableDiffusion) + Opm::VtkDiffusionModule::registerParameters(); + + if (enableEnergy) + Opm::VtkEnergyModule::registerParameters(); + + EWOMS_REGISTER_PARAM(TypeTag, Scalar, FlashTolerance, + "The maximum tolerance for the flash solver to " + "consider the solution converged"); + EWOMS_REGISTER_PARAM(TypeTag, int, FlashVerbosity, + "Flash solver verbosity level"); + EWOMS_REGISTER_PARAM(TypeTag, std::string, FlashTwoPhaseMethod, + "Method for solving vapor-liquid composition"); + } + + /*! + * \copydoc FvBaseDiscretization::name + */ + static std::string name() + { return "flash"; } + + /*! + * \copydoc FvBaseDiscretization::primaryVarName + */ + std::string primaryVarName(unsigned pvIdx) const + { + const std::string& tmp = EnergyModule::primaryVarName(pvIdx); + if (tmp != "") + return tmp; + + std::ostringstream oss; + if (Indices::z0Idx <= pvIdx && pvIdx < Indices::z0Idx + numComponents - 1) + oss << "z_," << FluidSystem::componentName(/*compIdx=*/pvIdx - Indices::z0Idx); + else if (pvIdx==Indices::pressure0Idx) + oss << "pressure_" << FluidSystem::phaseName(0); + else + assert(false); + + return oss.str(); + } + + /*! + * \copydoc FvBaseDiscretization::eqName + */ + std::string eqName(unsigned eqIdx) const + { + const std::string& tmp = EnergyModule::eqName(eqIdx); + if (tmp != "") + return tmp; + + std::ostringstream oss; + if (Indices::conti0EqIdx <= eqIdx && eqIdx < Indices::conti0EqIdx + + numComponents) { + unsigned compIdx = eqIdx - Indices::conti0EqIdx; + oss << "continuity^" << FluidSystem::componentName(compIdx); + } + else + assert(false); + + return oss.str(); + } + + /*! + * \copydoc FvBaseDiscretization::primaryVarWeight + */ + Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const + { + Scalar tmp = EnergyModule::primaryVarWeight(*this, globalDofIdx, pvIdx); + if (tmp > 0) + return tmp; + + unsigned compIdx = pvIdx - Indices::conti0EqIdx; + + // make all kg equal. also, divide the weight of all total + // compositions by 100 to make the relative errors more + // comparable to the ones of the other models (at 10% porosity + // the medium is fully saturated with water at atmospheric + // conditions if 100 kg/m^3 are present!) + return FluidSystem::molarMass(compIdx) / 100.0; + } + + /*! + * \copydoc FvBaseDiscretization::eqWeight + */ + Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const + { + Scalar tmp = EnergyModule::eqWeight(*this, globalDofIdx, eqIdx); + if (tmp > 0) + return tmp; + + unsigned compIdx = eqIdx - Indices::conti0EqIdx; + + // make all kg equal + return FluidSystem::molarMass(compIdx); + } + + void registerOutputModules_() + { + ParentType::registerOutputModules_(); + + // add the VTK output modules which are meaningful for the model + this->addOutputModule(new Opm::VtkCompositionModule(this->simulator_)); + if (enableDiffusion) + this->addOutputModule(new Opm::VtkDiffusionModule(this->simulator_)); + if (enableEnergy) + this->addOutputModule(new Opm::VtkEnergyModule(this->simulator_)); + } +}; + +} // namespace Opm + +#endif diff --git a/opm/models/ptflash/flashnewtonmethod.hh b/opm/models/ptflash/flashnewtonmethod.hh new file mode 100644 index 000000000..ffc4c0222 --- /dev/null +++ b/opm/models/ptflash/flashnewtonmethod.hh @@ -0,0 +1,144 @@ +// -*- 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::FlashNewtonMethod + */ +#ifndef EWOMS_FLASH_NEWTON_METHOD_HH +#define EWOMS_FLASH_NEWTON_METHOD_HH + +#include "flashproperties.hh" + +#include + +#include + +#include + +namespace Opm::Properties { + +template +struct DiscNewtonMethod; + +} // namespace Opm::Properties + +namespace Opm { +/*! + * \ingroup FlashModel + * + * \brief A Newton solver specific to the NCP model. + */ + +template +class FlashNewtonMethod : public GetPropType +{ + using ParentType = GetPropType; + + using PrimaryVariables = GetPropType; + using EqVector = GetPropType; + using Simulator = GetPropType; + using Scalar = GetPropType; + using Indices = GetPropType; + + enum { pressure0Idx = Indices::pressure0Idx }; + enum { z0Idx = Indices::z0Idx }; + enum { numComponents = getPropValue() }; + +public: + /*! + * \copydoc FvBaseNewtonMethod::FvBaseNewtonMethod(Problem& ) + */ + FlashNewtonMethod(Simulator& simulator) : ParentType(simulator) + {} + +protected: + friend ParentType; + friend NewtonMethod; + + /*! + * \copydoc FvBaseNewtonMethod::updatePrimaryVariables_ + */ + void updatePrimaryVariables_(unsigned globalDofIdx, + PrimaryVariables& nextValue, + const PrimaryVariables& currentValue, + const EqVector& update, + const EqVector& currentResidual) + { + // normal Newton-Raphson update + nextValue = currentValue; + nextValue -= update; + + //// + // Pressure updates + //// + // limit pressure reference change to 20% of the total value per iteration + clampValue_(nextValue[pressure0Idx], + currentValue[pressure0Idx]*0.8, + currentValue[pressure0Idx]*1.2); + + //// + // z updates + //// + // restrict update to at most 0.1 + Scalar maxDeltaZ = 0.0; // in update vector + Scalar sumDeltaZ = 0.0; // changes in last component (not in update vector) + for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) { + maxDeltaZ = std::max(std::abs(update[z0Idx + compIdx]), maxDeltaZ); + sumDeltaZ += update[z0Idx + compIdx]; + } + maxDeltaZ = std::max(std::abs(-sumDeltaZ), maxDeltaZ); + + // if max. update is above 0.1, restrict that one to 0.1 and adjust the rest accordingly (s.t. last comp. update is sum of the changes in update vector) + if (maxDeltaZ > 0.1) { + Scalar alpha = 0.2/maxDeltaZ; + for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) { + nextValue[z0Idx + compIdx] = currentValue[z0Idx + compIdx] - alpha*update[z0Idx + compIdx]; + } + } + + // ensure that z-values are less than tol or more than 1-tol + Scalar tol = 1e-8; + for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) { + clampValue_(nextValue[z0Idx + compIdx], tol, 1-tol); + } + + // normalize s.t. z sums to 1.0 + // Scalar lastZ = 1.0; + // Scalar sumZ = 0.0; + // for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) { + // lastZ -= nextValue[z0Idx + compIdx]; + // sumZ += nextValue[z0Idx + compIdx]; + // } + // sumZ += lastZ; + // for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) { + // nextValue[z0Idx + compIdx] /= sumZ; + // } + } +private: + void clampValue_(Scalar& val, Scalar minVal, Scalar maxVal) const + { val = std::max(minVal, std::min(val, maxVal)); } + +}; // class FlashNewtonMethod +} // namespace Opm +#endif \ No newline at end of file diff --git a/opm/models/ptflash/flashprimaryvariables.hh b/opm/models/ptflash/flashprimaryvariables.hh new file mode 100644 index 000000000..21b1a3a90 --- /dev/null +++ b/opm/models/ptflash/flashprimaryvariables.hh @@ -0,0 +1,160 @@ +// -*- 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::FlashPrimaryVariables + */ +#ifndef EWOMS_FLASH_PRIMARY_VARIABLES_HH +#define EWOMS_FLASH_PRIMARY_VARIABLES_HH + +#include "flashindices.hh" +#include "flashproperties.hh" + +#include +#include + +#include +#include +#include + +#include + +#include + +namespace Opm { + +/*! + * \ingroup FlashModel + * + * \brief Represents the primary variables used by the compositional + * flow model based on flash calculations. + * + * This class is basically a Dune::FieldVector which can retrieve its + * contents from an aribitatry fluid state. + */ +template +class FlashPrimaryVariables : public FvBasePrimaryVariables +{ + using ParentType = FvBasePrimaryVariables; + + using Scalar = GetPropType; + using Evaluation = GetPropType; + using MaterialLawParams = GetPropType; + using FluidSystem = GetPropType; + using Indices = GetPropType; + + // primary variable indices + enum { z0Idx = Indices::z0Idx }; + enum { pressure0Idx = Indices::pressure0Idx }; + + enum { numPhases = getPropValue() }; + enum { numComponents = getPropValue() }; + + using Toolbox = typename Opm::MathToolbox; + using EnergyModule = Opm::EnergyModule()>; + +public: + FlashPrimaryVariables() : ParentType() + { Opm::Valgrind::SetDefined(*this); } + + /*! + * \copydoc ImmisciblePrimaryVariables::ImmisciblePrimaryVariables(Scalar) + */ + FlashPrimaryVariables(Scalar value) : ParentType(value) + { + Opm::Valgrind::CheckDefined(value); + Opm::Valgrind::SetDefined(*this); + } + + /*! + * \copydoc ImmisciblePrimaryVariables::ImmisciblePrimaryVariables(const + * ImmisciblePrimaryVariables& ) + */ + FlashPrimaryVariables(const FlashPrimaryVariables& value) = default; + FlashPrimaryVariables& operator=(const FlashPrimaryVariables& value) = default; + + /*! + * \copydoc ImmisciblePrimaryVariables::assignMassConservative + */ + template + void assignMassConservative(const FluidState& fluidState, + const MaterialLawParams&, + bool = false) + { + // there is no difference between naive and mass conservative + // assignment in the flash model. (we only need the total + // concentrations.) + assignNaive(fluidState); + } + + /*! + * \copydoc ImmisciblePrimaryVariables::assignNaive + */ + template + void assignNaive(const FluidState& fluidState) + { + // reset everything + (*this) = 0.0; + + // assign the phase temperatures. this is out-sourced to + // the energy module + EnergyModule::setPriVarTemperatures(*this, fluidState); + + // determine the phase presence. + Dune::FieldVector z(0.0); + Scalar sumMoles = 0.0; + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + Scalar tmp = Opm::getValue(fluidState.molarity(phaseIdx, compIdx) * fluidState.saturation(phaseIdx)); + z[compIdx] += Opm::max(tmp, 1e-8); + sumMoles += tmp; + } + } + z /= sumMoles; + + for (int i = 0; i < numComponents - 1; ++i) + (*this)[z0Idx + i] = z[i]; + (*this)[pressure0Idx] = Opm::getValue(fluidState.pressure(0)); + } + + /*! + * \brief Prints the names of the primary variables and their values. + * + * \param os The \c std::ostream which should be used for the output. + */ + void print(std::ostream& os = std::cout) const + { + os << "(p_" << FluidSystem::phaseName(0) << " = " + << this->operator[](pressure0Idx); + for (unsigned compIdx = 0; compIdx < numComponents - 2; ++compIdx) { + os << ", z_" << FluidSystem::componentName(compIdx) << " = " + << this->operator[](z0Idx + compIdx); + } + os << ")" << std::flush; + } +}; + +} // namespace Opm + +#endif diff --git a/opm/models/ptflash/flashproperties.hh b/opm/models/ptflash/flashproperties.hh new file mode 100644 index 000000000..7f682a432 --- /dev/null +++ b/opm/models/ptflash/flashproperties.hh @@ -0,0 +1,55 @@ +// -*- 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 + * \ingroup FlashModel + * + * \brief Declares the properties required by the compositional + * multi-phase model based on flash calculations. + */ +#ifndef EWOMS_FLASH_PROPERTIES_HH +#define EWOMS_FLASH_PROPERTIES_HH + +#include +#include +#include +#include + +namespace Opm::Properties { + +//! The type of the flash constraint solver +template +struct FlashSolver { using type = UndefinedProperty; }; +//! The maximum accepted error of the flash solver +template +struct FlashTolerance { using type = UndefinedProperty; }; +//! The verbosity level of the flash solver +template +struct FlashVerbosity { using type = UndefinedProperty; }; +//! Two-phase flash method +template +struct FlashTwoPhaseMethod { using type = UndefinedProperty; }; + +} // namespace Opm::Properties + +#endif diff --git a/opm/models/ptflash/flashratevector.hh b/opm/models/ptflash/flashratevector.hh new file mode 100644 index 000000000..8fe1a83b4 --- /dev/null +++ b/opm/models/ptflash/flashratevector.hh @@ -0,0 +1,143 @@ +// -*- 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::FlashRateVector + */ +#ifndef EWOMS_FLASH_RATE_VECTOR_HH +#define EWOMS_FLASH_RATE_VECTOR_HH + +#include + +#include +#include + +#include "flashintensivequantities.hh" + +namespace Opm { + +/*! + * \ingroup FlashModel + * + * \copydoc Opm::ImmiscibleRateVector + */ +template +class FlashRateVector + : public Dune::FieldVector, + getPropValue()> +{ + using Scalar = GetPropType; + using Evaluation = GetPropType; + using FluidSystem = GetPropType; + using Indices = GetPropType; + + enum { conti0EqIdx = Indices::conti0EqIdx }; + enum { numComponents = getPropValue() }; + enum { numEq = getPropValue() }; + + using ParentType = Dune::FieldVector; + using EnergyModule = Opm::EnergyModule()>; + +public: + FlashRateVector() : ParentType() + { Opm::Valgrind::SetUndefined(*this); } + + /*! + * \copydoc ImmiscibleRateVector::ImmiscibleRateVector(Scalar) + */ + FlashRateVector(const Evaluation& value) : ParentType(value) + {} + + /*! + * \copydoc ImmiscibleRateVector::ImmiscibleRateVector(const + * ImmiscibleRateVector& ) + */ + FlashRateVector(const FlashRateVector& value) : ParentType(value) + {} + + /*! + * \copydoc ImmiscibleRateVector::setMassRate + */ + void setMassRate(const ParentType& value) + { + // convert to molar rates + ParentType molarRate(value); + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) + molarRate[conti0EqIdx + compIdx] /= FluidSystem::molarMass(compIdx); + + setMolarRate(molarRate); + } + + /*! + * \copydoc ImmiscibleRateVector::setMolarRate + */ + void setMolarRate(const ParentType& value) + { ParentType::operator=(value); } + + /*! + * \copydoc ImmiscibleRateVector::setEnthalpyRate + */ + void setEnthalpyRate(const Evaluation& rate) + { EnergyModule::setEnthalpyRate(*this, rate); } + + /*! + * \copydoc ImmiscibleRateVector::setVolumetricRate + */ + template + void setVolumetricRate(const FluidState& fluidState, unsigned phaseIdx, const RhsEval& volume) + { + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) + (*this)[conti0EqIdx + compIdx] = + fluidState.density(phaseIdx, compIdx) + * fluidState.moleFraction(phaseIdx, compIdx) + * volume; + + EnergyModule::setEnthalpyRate(*this, fluidState, phaseIdx, volume); + } + + /*! + * \brief Assignment operator from a scalar or a function evaluation + */ + template + FlashRateVector& operator=(const RhsEval& value) + { + for (unsigned i=0; i < this->size(); ++i) + (*this)[i] = value; + return *this; + } + + /*! + * \brief Assignment operator from another rate vector + */ + FlashRateVector& operator=(const FlashRateVector& other) + { + for (unsigned i=0; i < this->size(); ++i) + (*this)[i] = other[i]; + return *this; + } +}; + +} // namespace Opm + +#endif diff --git a/opm/simulators/linalg/parallelbasebackend.hh b/opm/simulators/linalg/parallelbasebackend.hh index 3319a3a18..2d248ac42 100644 --- a/opm/simulators/linalg/parallelbasebackend.hh +++ b/opm/simulators/linalg/parallelbasebackend.hh @@ -37,6 +37,8 @@ #include #include #include +//#include +//#include #include #include @@ -261,6 +263,9 @@ public: ParallelScalarProduct parScalarProduct(overlappingMatrix_->overlap()); ParallelOperator parOperator(*overlappingMatrix_); + //const auto& matrix = parOperator.getMatrix(); + //Dune::storeMatrixMarket(*overlappingMatrix_, std::string("mymatrix.mm")); + // retrieve the linear solver auto solver = asImp_().prepareSolver_(parOperator, parScalarProduct, @@ -278,6 +283,7 @@ public: // copy the result back to the non-overlapping vector overlappingx_->assignTo(x); + overlappingx_->print(); // return the result of the solver return result.first; diff --git a/opm/simulators/linalg/parallelbicgstabbackend.hh b/opm/simulators/linalg/parallelbicgstabbackend.hh index 0ea3885a4..5eda6600c 100644 --- a/opm/simulators/linalg/parallelbicgstabbackend.hh +++ b/opm/simulators/linalg/parallelbicgstabbackend.hh @@ -156,6 +156,9 @@ protected: bicgstabSolver->setMaxIterations(EWOMS_GET_PARAM(TypeTag, int, LinearSolverMaxIterations)); bicgstabSolver->setLinearOperator(&parOperator); bicgstabSolver->setRhs(this->overlappingb_); + // const auto& matrix = parOperator.getMatrix(); + // Dune::storeMatrixMarket(matrix, "mymatrix" + ".mm"); + return bicgstabSolver; } diff --git a/tests/simpletest_ptflash_ecfv.cc b/tests/simpletest_ptflash_ecfv.cc new file mode 100644 index 000000000..5eef271d1 --- /dev/null +++ b/tests/simpletest_ptflash_ecfv.cc @@ -0,0 +1,64 @@ +// -*- 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 + */ +#include "config.h" + +#include +//#include +//#include +#include "problems/simpletestproblem.hh" + + +namespace Opm::Properties { + +namespace TTag { + struct SimpleTestEcfvProblem +{ + using InheritsFrom = std::tuple; +}; +} + +template +struct SpatialDiscretizationSplice +{ + using type = TTag::EcfvDiscretization; +}; +//TESTAD +template +struct LocalLinearizerSplice +{ + using type = TTag::AutoDiffLocalLinearizer; +}; + + +} // namespace Opm::Properties + +int main(int argc, char **argv) +{ + using EcfvProblemTypeTag = Opm::Properties::TTag::SimpleTestEcfvProblem; + return Opm::start(argc, argv); +} \ No newline at end of file