From ccc452f4d403c4963f873a12a0663bbd8264e500 Mon Sep 17 00:00:00 2001 From: David Landa Marban Date: Mon, 24 Mar 2025 21:55:47 +0100 Subject: [PATCH 1/3] Updating MICP --- flow/flow_micp.cpp | 13 + ...BlackOilIntensiveQuantitiesGlobalIndex.hpp | 5 + opm/models/blackoil/blackoilbrinemodules.hh | 2 +- .../blackoil/blackoildiffusionmodule.hh | 270 +++++++----- .../blackoil/blackoildispersionmodule.hh | 184 ++++---- opm/models/blackoil/blackoilindices.hh | 32 +- .../blackoil/blackoilintensivequantities.hh | 19 +- .../blackoil/blackoillocalresidualtpfa.hh | 40 +- opm/models/blackoil/blackoilmicpmodules.hh | 396 ++++++++++-------- opm/models/blackoil/blackoilmicpparams.cpp | 96 +++-- opm/models/blackoil/blackoilmicpparams.hpp | 40 +- opm/models/blackoil/blackoilnewtonmethod.hpp | 37 +- .../blackoil/blackoilonephaseindices.hh | 32 +- .../blackoil/blackoiltwophaseindices.hh | 32 +- .../discretization/common/tpfalinearizer.hh | 13 +- opm/models/io/vtkblackoilmicpmodule.hpp | 3 +- opm/simulators/flow/FIPContainer.cpp | 85 ++++ opm/simulators/flow/FIPContainer.hpp | 20 + opm/simulators/flow/FlowGenericProblem.hpp | 12 +- opm/simulators/flow/FlowProblem.hpp | 17 +- opm/simulators/flow/FlowProblemBlackoil.hpp | 58 ++- opm/simulators/flow/FlowProblemComp.hpp | 15 + .../flow/GenericOutputBlackoilModule.cpp | 2 +- opm/simulators/flow/OutputBlackoilModule.hpp | 117 +++++- opm/simulators/utils/BlackoilPhases.cpp | 1 + opm/simulators/utils/BlackoilPhases.hpp | 5 +- opm/simulators/utils/phaseUsageFromDeck.cpp | 22 + opm/simulators/utils/phaseUsageFromDeck.hpp | 3 + .../wells/MultisegmentWell_impl.hpp | 4 + opm/simulators/wells/PerfData.cpp | 16 +- opm/simulators/wells/PerfData.hpp | 8 +- opm/simulators/wells/SingleWellState.cpp | 18 + opm/simulators/wells/SingleWellState.hpp | 3 + .../wells/StandardWellConnections.cpp | 12 +- .../wells/StandardWellConnections.hpp | 5 +- opm/simulators/wells/StandardWell_impl.hpp | 5 +- opm/simulators/wells/WellInterfaceGeneric.cpp | 7 +- opm/simulators/wells/WellInterface_impl.hpp | 6 - opm/simulators/wells/WellState.cpp | 14 + 39 files changed, 1119 insertions(+), 550 deletions(-) diff --git a/flow/flow_micp.cpp b/flow/flow_micp.cpp index c10fa4504d..aeef32e657 100644 --- a/flow/flow_micp.cpp +++ b/flow/flow_micp.cpp @@ -24,6 +24,8 @@ #include #include #include +#include +#include namespace Opm { namespace Properties { @@ -36,6 +38,10 @@ template struct EnableMICP { static constexpr bool value = true; }; +template +struct Linearizer { using type = TpfaLinearizer; }; +template +struct LocalResidual { using type = BlackOilLocalResidualTPFA; }; //! The indices required by the model template struct Indices @@ -58,6 +64,13 @@ public: /*enabledCompIdx=*/FluidSystem::waterCompIdx, 5>; //Five MICP components }; + +template +struct EnableDiffusion { static constexpr bool value = true; }; + +template +struct EnableDispersion { static constexpr bool value = true; }; + }} namespace Opm { diff --git a/flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp b/flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp index 20f5b0d308..d20e805a7b 100644 --- a/flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp +++ b/flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp @@ -512,6 +512,11 @@ public: Scalar referencePorosity() const { return referencePorosity_; } + const Evaluation& permFactor() const + { + throw std::logic_error("permFactor() is not yet implemented for compositional modeling"); + } + private: friend BlackOilSolventIntensiveQuantities; friend BlackOilExtboIntensiveQuantities; diff --git a/opm/models/blackoil/blackoilbrinemodules.hh b/opm/models/blackoil/blackoilbrinemodules.hh index 204a75bffe..2e4281de72 100644 --- a/opm/models/blackoil/blackoilbrinemodules.hh +++ b/opm/models/blackoil/blackoilbrinemodules.hh @@ -486,6 +486,6 @@ public: }; -} // namespace Ewoms +} // namespace Opm #endif diff --git a/opm/models/blackoil/blackoildiffusionmodule.hh b/opm/models/blackoil/blackoildiffusionmodule.hh index 4d77be5917..bb52c0bb31 100644 --- a/opm/models/blackoil/blackoildiffusionmodule.hh +++ b/opm/models/blackoil/blackoildiffusionmodule.hh @@ -28,6 +28,7 @@ #ifndef OPM_BLACKOIL_DIFFUSION_MODULE_HH #define OPM_BLACKOIL_DIFFUSION_MODULE_HH +#include #include #include @@ -100,11 +101,19 @@ class BlackOilDiffusionModule using RateVector = GetPropType; using FluidSystem = GetPropType; using Indices = GetPropType; + using IntensiveQuantities = GetPropType; enum { numPhases = FluidSystem::numPhases }; enum { numComponents = FluidSystem::numComponents }; enum { conti0EqIdx = Indices::conti0EqIdx }; + enum { enableMICP = getPropValue() }; + + static constexpr unsigned contiMicrobialEqIdx = Indices::contiMicrobialEqIdx; + static constexpr unsigned contiOxygenEqIdx = Indices::contiOxygenEqIdx; + static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx; + static constexpr unsigned contiUreaEqIdx = Indices::contiUreaEqIdx; + using Toolbox = MathToolbox; public: @@ -157,86 +166,113 @@ public: if(!FluidSystem::enableDiffusion()) return; const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx); - const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState(); - const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState(); + const auto& inIq = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx); + const auto& exIq = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx); const auto& diffusivity = extQuants.diffusivity(); const auto& effectiveDiffusionCoefficient = extQuants.effectiveDiffusionCoefficient(); - addDiffusiveFlux(flux, fluidStateI, fluidStateJ, diffusivity, effectiveDiffusionCoefficient); + addDiffusiveFlux(flux, inIq, exIq, diffusivity, effectiveDiffusionCoefficient); } - template + template static void addDiffusiveFlux(RateVector& flux, - const FluidState& fluidStateI, - const FluidState& fluidStateJ, + const IntensiveQuantities& inIq, + const IntensiveQuantities& exIq, const Evaluation& diffusivity, const EvaluationArray& effectiveDiffusionCoefficient) { - unsigned pvtRegionIndex = fluidStateI.pvtRegionIndex(); - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } + const auto& inFs = inIq.fluidState(); + const auto& exFs = exIq.fluidState(); + Evaluation diffR = 0.0; + if constexpr(enableMICP) { + // The diffusion coefficients are given for mass concentrations (urea scaled by a factor of 10) + Evaluation bAvg = (inFs.invB(waterPhaseIdx) + Toolbox::value(exFs.invB(waterPhaseIdx))) / 2; + diffR = inIq.microbialConcentration() - Toolbox::value(exIq.microbialConcentration()); + flux[contiMicrobialEqIdx] += + bAvg + * diffR + * diffusivity + * effectiveDiffusionCoefficient[waterPhaseIdx][0]; + diffR = inIq.oxygenConcentration() - Toolbox::value(exIq.oxygenConcentration()); + flux[contiOxygenEqIdx] += + bAvg + * diffR + * diffusivity + * effectiveDiffusionCoefficient[waterPhaseIdx][1]; + diffR = inIq.ureaConcentration() - Toolbox::value(exIq.ureaConcentration()); + // dividing by scaling factor 10 (see WellInterfaceGeneric.cpp) + flux[contiUreaEqIdx] += + bAvg + * diffR + * diffusivity + * effectiveDiffusionCoefficient[waterPhaseIdx][2] / 10; + } + else { + unsigned pvtRegionIndex = inFs.pvtRegionIndex(); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } - // no diffusion in water for blackoil models - if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { - continue; - } + // no diffusion in water for blackoil models + if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { + continue; + } - // no diffusion in gas phase in water + gas system. - if (FluidSystem::gasPhaseIdx == phaseIdx && !FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - continue; - } + // no diffusion in gas phase in water + gas system. + if (FluidSystem::gasPhaseIdx == phaseIdx && !FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + continue; + } - // arithmetic mean of the phase's b factor weighed by saturation - Evaluation bSAvg = fluidStateI.saturation(phaseIdx) * fluidStateI.invB(phaseIdx); - bSAvg += Toolbox::value(fluidStateJ.saturation(phaseIdx)) * Toolbox::value(fluidStateJ.invB(phaseIdx)); - bSAvg /= 2; + // arithmetic mean of the phase's b factor weighed by saturation + Evaluation bSAvg = inFs.saturation(phaseIdx) * inFs.invB(phaseIdx); + bSAvg += Toolbox::value(exFs.saturation(phaseIdx)) * Toolbox::value(exFs.invB(phaseIdx)); + bSAvg /= 2; - // phase not present, skip - if(bSAvg < 1.0e-6) - continue; - Evaluation convFactor = 1.0; - Evaluation diffR = 0.0; - if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) { - Evaluation rsAvg = (fluidStateI.Rs() + Toolbox::value(fluidStateJ.Rs())) / 2; - convFactor = 1.0 / (toFractionGasOil(pvtRegionIndex) + rsAvg); - diffR = fluidStateI.Rs() - Toolbox::value(fluidStateJ.Rs()); - } - if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) { - Evaluation rvAvg = (fluidStateI.Rv() + Toolbox::value(fluidStateJ.Rv())) / 2; - convFactor = toFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toFractionGasOil(pvtRegionIndex)); - diffR = fluidStateI.Rv() - Toolbox::value(fluidStateJ.Rv()); - } - if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) { - Evaluation rsAvg = (fluidStateI.Rsw() + Toolbox::value(fluidStateJ.Rsw())) / 2; - convFactor = 1.0 / (toFractionGasWater(pvtRegionIndex) + rsAvg); - diffR = fluidStateI.Rsw() - Toolbox::value(fluidStateJ.Rsw()); - } - if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) { - Evaluation rvAvg = (fluidStateI.Rvw() + Toolbox::value(fluidStateJ.Rvw())) / 2; - convFactor = toFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toFractionGasWater(pvtRegionIndex)); - diffR = fluidStateI.Rvw() - Toolbox::value(fluidStateJ.Rvw()); - } + // phase not present, skip + if(bSAvg < 1.0e-6) + continue; + Evaluation convFactor = 1.0; + if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) { + Evaluation rsAvg = (inFs.Rs() + Toolbox::value(exFs.Rs())) / 2; + convFactor = 1.0 / (toFractionGasOil(pvtRegionIndex) + rsAvg); + diffR = inFs.Rs() - Toolbox::value(exFs.Rs()); + } + if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) { + Evaluation rvAvg = (inFs.Rv() + Toolbox::value(exFs.Rv())) / 2; + convFactor = toFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toFractionGasOil(pvtRegionIndex)); + diffR = inFs.Rv() - Toolbox::value(exFs.Rv()); + } + if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) { + Evaluation rsAvg = (inFs.Rsw() + Toolbox::value(exFs.Rsw())) / 2; + convFactor = 1.0 / (toFractionGasWater(pvtRegionIndex) + rsAvg); + diffR = inFs.Rsw() - Toolbox::value(exFs.Rsw()); + } + if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) { + Evaluation rvAvg = (inFs.Rvw() + Toolbox::value(exFs.Rvw())) / 2; + convFactor = toFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toFractionGasWater(pvtRegionIndex)); + diffR = inFs.Rvw() - Toolbox::value(exFs.Rvw()); + } - // mass flux of solvent component (oil in oil or gas in gas) - unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx); - unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx); - flux[conti0EqIdx + activeSolventCompIdx] += - - bSAvg - * convFactor - * diffR - * diffusivity - * effectiveDiffusionCoefficient[phaseIdx][solventCompIdx]; + // mass flux of solvent component (oil in oil or gas in gas) + unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx); + unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx); + flux[conti0EqIdx + activeSolventCompIdx] += + - bSAvg + * convFactor + * diffR + * diffusivity + * effectiveDiffusionCoefficient[phaseIdx][solventCompIdx]; - // mass flux of solute component (gas in oil or oil in gas) - unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx); - unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx); - flux[conti0EqIdx + activeSoluteCompIdx] += - bSAvg - * diffR - * convFactor - * diffusivity - * effectiveDiffusionCoefficient[phaseIdx][soluteCompIdx]; + // mass flux of solute component (gas in oil or oil in gas) + unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx); + unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx); + flux[conti0EqIdx + activeSoluteCompIdx] += + bSAvg + * diffR + * convFactor + * diffusivity + * effectiveDiffusionCoefficient[phaseIdx][soluteCompIdx]; + } } } @@ -339,8 +375,11 @@ class BlackOilDiffusionIntensiveQuantities using ElementContext = GetPropType; using FluidSystem = GetPropType; using IntensiveQuantities = GetPropType; + using MICPModule = BlackOilMICPModule; enum { numPhases = FluidSystem::numPhases }; enum { numComponents = FluidSystem::numComponents }; + enum { enableMICP = getPropValue() }; + static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx; public: BlackOilDiffusionIntensiveQuantities() = default; @@ -418,34 +457,42 @@ protected: const IntensiveQuantities& intQuants) { using Toolbox = MathToolbox; - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } + if constexpr(enableMICP) { + unsigned pvtRegionIndex = intQuants.fluidState().pvtRegionIndex(); + diffusionCoefficient_[waterPhaseIdx][0] = MICPModule::microbialDiffusion(pvtRegionIndex); + diffusionCoefficient_[waterPhaseIdx][1] = MICPModule::oxygenDiffusion(pvtRegionIndex); + diffusionCoefficient_[waterPhaseIdx][2] = MICPModule::ureaDiffusion(pvtRegionIndex); + } + else { + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } - // no diffusion in water for blackoil models - if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { - continue; - } + // no diffusion in water for blackoil models + if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { + continue; + } - // Based on Millington, R. J., & Quirk, J. P. (1961). - // \Note: it is possible to use NumericalConstants later - // constexpr auto& numconst = GetPropValue; - constexpr double myeps = 0.0001; //numconst.blackoildiffusionmoduleeps; - const Evaluation& base = - Toolbox::max(myeps, //0.0001, - intQuants.porosity() - * intQuants.fluidState().saturation(phaseIdx)); - tortuosity_[phaseIdx] = - 1.0 / (intQuants.porosity() * intQuants.porosity()) - * Toolbox::pow(base, 10.0/3.0); + // Based on Millington, R. J., & Quirk, J. P. (1961). + // \Note: it is possible to use NumericalConstants later + // constexpr auto& numconst = GetPropValue; + constexpr double myeps = 0.0001; //numconst.blackoildiffusionmoduleeps; + const Evaluation& base = + Toolbox::max(myeps, //0.0001, + intQuants.porosity() + * intQuants.fluidState().saturation(phaseIdx)); + tortuosity_[phaseIdx] = + 1.0 / (intQuants.porosity() * intQuants.porosity()) + * Toolbox::pow(base, 10.0/3.0); - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - diffusionCoefficient_[phaseIdx][compIdx] = - FluidSystem::diffusionCoefficient(fluidState, - paramCache, - phaseIdx, - compIdx); + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + diffusionCoefficient_[phaseIdx][compIdx] = + FluidSystem::diffusionCoefficient(fluidState, + paramCache, + phaseIdx, + compIdx); + } } } } @@ -535,6 +582,8 @@ class BlackOilDiffusionExtensiveQuantities enum { dimWorld = GridView::dimensionworld }; enum { numPhases = getPropValue() }; enum { numComponents = getPropValue() }; + enum { enableMICP = getPropValue() }; + static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx; using DimVector = Dune::FieldVector; using DimEvalVector = Dune::FieldVector; @@ -569,21 +618,30 @@ public: const IntensiveQuantities& intQuantsInside, const IntensiveQuantities& intQuantsOutside) { // opm-models expects per area flux - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - // no diffusion in water for blackoil models - if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { - continue; - } + // use the arithmetic average for the effective + // diffusion coefficients. + if constexpr(enableMICP) { for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - // use the arithmetic average for the effective - // diffusion coefficients. - effectiveDiffusionCoefficient[phaseIdx][compIdx] = 0.5 * - ( intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx) + - intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx) ); - Valgrind::CheckDefined(effectiveDiffusionCoefficient[phaseIdx][compIdx]); + effectiveDiffusionCoefficient[waterPhaseIdx][compIdx] = 0.5 * + ( intQuantsInside.effectiveDiffusionCoefficient(waterPhaseIdx, compIdx) + + intQuantsOutside.effectiveDiffusionCoefficient(waterPhaseIdx, compIdx) ); + } + } + else { + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + // no diffusion in water for blackoil models + if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { + continue; + } + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + effectiveDiffusionCoefficient[phaseIdx][compIdx] = 0.5 * + ( intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx) + + intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx) ); + Valgrind::CheckDefined(effectiveDiffusionCoefficient[phaseIdx][compIdx]); + } } } } @@ -608,7 +666,7 @@ public: { return diffusivity_; } /*! - * \brief The effective diffusion coeffcient of a component in a + * \brief The effective diffusion coefficient of a component in a * fluid phase at the face's integration point * * \copydoc Doxygen::phaseIdxParam diff --git a/opm/models/blackoil/blackoildispersionmodule.hh b/opm/models/blackoil/blackoildispersionmodule.hh index c524b75356..70b5554c7a 100644 --- a/opm/models/blackoil/blackoildispersionmodule.hh +++ b/opm/models/blackoil/blackoildispersionmodule.hh @@ -62,6 +62,7 @@ class BlackOilDispersionModule using RateVector = GetPropType; using FluidSystem = GetPropType; using Evaluation = GetPropType; + using IntensiveQuantities = GetPropType; enum { numPhases = FluidSystem::numPhases }; @@ -85,10 +86,10 @@ public: unsigned) {} - template + template static void addDispersiveFlux(RateVector&, - const FluidState&, - const FluidState&, + const IntensiveQuantities&, + const IntensiveQuantities&, const Evaluation&, const Scalar&) {} @@ -116,6 +117,12 @@ class BlackOilDispersionModule enum { numComponents = FluidSystem::numComponents }; enum { conti0EqIdx = Indices::conti0EqIdx }; enum { enableDispersion = getPropValue() }; + enum { enableMICP = getPropValue() }; + + static constexpr unsigned contiMicrobialEqIdx = Indices::contiMicrobialEqIdx; + static constexpr unsigned contiOxygenEqIdx = Indices::contiOxygenEqIdx; + static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx; + static constexpr unsigned contiUreaEqIdx = Indices::contiUreaEqIdx; using Toolbox = MathToolbox; @@ -148,11 +155,11 @@ public: return; } const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx); - const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState(); - const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState(); + const auto& inIq = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx); + const auto& exIq = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx); const auto& dispersivity = extQuants.dispersivity(); const auto& normVelocityAvg = extQuants.normVelocityAvg(); - addDispersiveFlux(flux, fluidStateI, fluidStateJ, dispersivity, normVelocityAvg); + addDispersiveFlux(flux, inIq, exIq, dispersivity, normVelocityAvg); } /*! @@ -175,86 +182,105 @@ public: * \f$D_{ij}\f$ is computed in ecltransmissibility_impl.hh, using the dispersion coefficients \f$E_i\f$ * and \f$E_j\f$. */ - template + template static void addDispersiveFlux(RateVector& flux, - const FluidState& fluidStateI, - const FluidState& fluidStateJ, + const IntensiveQuantities& inIq, + const IntensiveQuantities& exIq, const Evaluation& dispersivity, const Scalar& normVelocityAvg) { - unsigned pvtRegionIndex = fluidStateI.pvtRegionIndex(); - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } + const auto& inFs = inIq.fluidState(); + const auto& exFs = exIq.fluidState(); + Evaluation diffR = 0.0; + if constexpr(enableMICP) { + // The dispersion coefficients are given for mass concentrations + Evaluation bAvg = (inFs.invB(waterPhaseIdx) + Toolbox::value(exFs.invB(waterPhaseIdx))) / 2; + diffR = inIq.microbialConcentration()- Toolbox::value(exIq.microbialConcentration()); + flux[contiMicrobialEqIdx] += + bAvg + * normVelocityAvg[waterPhaseIdx] + * dispersivity + * diffR; + diffR = inIq.oxygenConcentration()- Toolbox::value(exIq.oxygenConcentration()); + flux[contiOxygenEqIdx] += + bAvg + * normVelocityAvg[waterPhaseIdx] + * dispersivity + * diffR; + diffR = inIq.ureaConcentration()- Toolbox::value(exIq.ureaConcentration()); + // dividing by scaling factor 10 (see WellInterfaceGeneric.cpp) + flux[contiUreaEqIdx] += + bAvg + * normVelocityAvg[waterPhaseIdx] + * dispersivity + * diffR / 10; + } + else { + unsigned pvtRegionIndex = inFs.pvtRegionIndex(); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } - // no dispersion in water for blackoil models unless water can contain dissolved gas - if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { - continue; - } + // no dispersion in water for blackoil models unless water can contain dissolved gas + if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { + continue; + } - // Adding dispersion in the gas phase leads to - // convergence issues and unphysical results. - // We disable dispersion in the gas phase for now - // See comment below - if (FluidSystem::gasPhaseIdx == phaseIdx) { - continue; - } + // adding dispersion in the gas phase leads to + // convergence issues and unphysical results. + // we disable dispersion in the gas phase for now + if (FluidSystem::gasPhaseIdx == phaseIdx) { + continue; + } - // no dispersion in gas for blackoil models unless gas can contain evaporated water or oil - // phase check disabled due to if above, reenable when removing unconditional gas phase disablement - if ((!FluidSystem::enableVaporizedWater() && !FluidSystem::enableVaporizedOil()) - /*&& FluidSystem::gasPhaseIdx == phaseIdx*/) { - continue; - } + // arithmetic mean of the phase's b factor + Evaluation bAvg = inFs.invB(phaseIdx); + bAvg += Toolbox::value(exFs.invB(phaseIdx)); + bAvg /= 2; - // arithmetic mean of the phase's b factor - Evaluation bAvg = fluidStateI.invB(phaseIdx); - bAvg += Toolbox::value(fluidStateJ.invB(phaseIdx)); - bAvg /= 2; + Evaluation convFactor = 1.0; + if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) { + Evaluation rsAvg = (inFs.Rs() + Toolbox::value(exFs.Rs())) / 2; + convFactor = 1.0 / (toMassFractionGasOil(pvtRegionIndex) + rsAvg); + diffR = inFs.Rs() - Toolbox::value(exFs.Rs()); + } + if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) { + Evaluation rvAvg = (inFs.Rv() + Toolbox::value(exFs.Rv())) / 2; + convFactor = toMassFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toMassFractionGasOil(pvtRegionIndex)); + diffR = inFs.Rv() - Toolbox::value(exFs.Rv()); + } + if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) { + Evaluation rsAvg = (inFs.Rsw() + Toolbox::value(exFs.Rsw())) / 2; + convFactor = 1.0 / (toMassFractionGasWater(pvtRegionIndex) + rsAvg); + diffR = inFs.Rsw() - Toolbox::value(exFs.Rsw()); + } + if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) { + Evaluation rvAvg = (inFs.Rvw() + Toolbox::value(exFs.Rvw())) / 2; + convFactor = toMassFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toMassFractionGasWater(pvtRegionIndex)); + diffR = inFs.Rvw() - Toolbox::value(exFs.Rvw()); + } - Evaluation convFactor = 1.0; - Evaluation diffR = 0.0; - if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) { - Evaluation rsAvg = (fluidStateI.Rs() + Toolbox::value(fluidStateJ.Rs())) / 2; - convFactor = 1.0 / (toMassFractionGasOil(pvtRegionIndex) + rsAvg); - diffR = fluidStateI.Rs() - Toolbox::value(fluidStateJ.Rs()); - } - if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) { - Evaluation rvAvg = (fluidStateI.Rv() + Toolbox::value(fluidStateJ.Rv())) / 2; - convFactor = toMassFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toMassFractionGasOil(pvtRegionIndex)); - diffR = fluidStateI.Rv() - Toolbox::value(fluidStateJ.Rv()); - } - if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) { - Evaluation rsAvg = (fluidStateI.Rsw() + Toolbox::value(fluidStateJ.Rsw())) / 2; - convFactor = 1.0 / (toMassFractionGasWater(pvtRegionIndex) + rsAvg); - diffR = fluidStateI.Rsw() - Toolbox::value(fluidStateJ.Rsw()); - } - if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) { - Evaluation rvAvg = (fluidStateI.Rvw() + Toolbox::value(fluidStateJ.Rvw())) / 2; - convFactor = toMassFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toMassFractionGasWater(pvtRegionIndex)); - diffR = fluidStateI.Rvw() - Toolbox::value(fluidStateJ.Rvw()); - } + // mass flux of solvent component + unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx); + unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx); + flux[conti0EqIdx + activeSolventCompIdx] += + - bAvg + * normVelocityAvg[phaseIdx] + * convFactor + * dispersivity + * diffR; - // mass flux of solvent component - unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx); - unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx); - flux[conti0EqIdx + activeSolventCompIdx] += - - bAvg - * normVelocityAvg[phaseIdx] - * convFactor - * dispersivity - * diffR; - - // mass flux of solute component - unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx); - unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx); - flux[conti0EqIdx + activeSoluteCompIdx] += - bAvg - * normVelocityAvg[phaseIdx] - * convFactor - * dispersivity - * diffR; + // mass flux of solute component + unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx); + unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx); + flux[conti0EqIdx + activeSoluteCompIdx] += + bAvg + * normVelocityAvg[phaseIdx] + * convFactor + * dispersivity + * diffR; + } } } @@ -498,8 +524,10 @@ public: if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { continue; } - // no dispersion in gas for blackoil models unless gas can contain evaporated water or oil - if ((!FluidSystem::enableVaporizedWater() && !FluidSystem::enableVaporizedOil()) && FluidSystem::gasPhaseIdx == phaseIdx) { + // adding dispersion in the gas phase leads to + // convergence issues and unphysical results. + // we disable dispersion in the gas phase for now + if (FluidSystem::gasPhaseIdx == phaseIdx) { continue; } // use the arithmetic average for the effective diff --git a/opm/models/blackoil/blackoilindices.hh b/opm/models/blackoil/blackoilindices.hh index 53f4ed6bba..d0760ea60f 100644 --- a/opm/models/blackoil/blackoilindices.hh +++ b/opm/models/blackoil/blackoilindices.hh @@ -147,24 +147,24 @@ struct BlackOilIndices numPolymers > 1 ? polymerConcentrationIdx + 1 : -1000; //! Index of the primary variable for the first MICP component - static constexpr int microbialConcentrationIdx = + static constexpr int biofilmConcentrationIdx = enableMICP ? PVOffset + numPhases + numSolvents : -1000; //! Index of the primary variable for the second MICP component - static constexpr int oxygenConcentrationIdx = - numMICPs > 1 ? microbialConcentrationIdx + 1 : -1000; + static constexpr int microbialConcentrationIdx = + numMICPs > 1 ? biofilmConcentrationIdx + 1 : -1000; //! Index of the primary variable for the third MICP component - static constexpr int ureaConcentrationIdx = - numMICPs > 2 ? oxygenConcentrationIdx + 1 : -1000; + static constexpr int oxygenConcentrationIdx = + numMICPs > 2 ? microbialConcentrationIdx + 1 : -1000; //! Index of the primary variable for the fourth MICP component - static constexpr int biofilmConcentrationIdx = - numMICPs > 3 ? ureaConcentrationIdx + 1 : -1000; + static constexpr int ureaConcentrationIdx = + numMICPs > 3 ? oxygenConcentrationIdx + 1 : -1000; //! Index of the primary variable for the fifth MICP component static constexpr int calciteConcentrationIdx = - numMICPs > 4 ? biofilmConcentrationIdx + 1 : -1000; + numMICPs > 4 ? ureaConcentrationIdx + 1 : -1000; //! Index of the primary variable for the foam static constexpr int foamConcentrationIdx = @@ -204,24 +204,24 @@ struct BlackOilIndices numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000; //! Index of the continuity equation for the first MICP component - static constexpr int contiMicrobialEqIdx = + static constexpr int contiBiofilmEqIdx = enableMICP ? PVOffset + numPhases + numSolvents + numExtbos : -1000; //! Index of the continuity equation for the second MICP component - static constexpr int contiOxygenEqIdx = - numMICPs > 1 ? contiMicrobialEqIdx + 1 : -1000; + static constexpr int contiMicrobialEqIdx = + numMICPs > 1 ? contiBiofilmEqIdx + 1 : -1000; //! Index of the continuity equation for the third MICP component - static constexpr int contiUreaEqIdx = - numMICPs > 2 ? contiOxygenEqIdx + 1 : -1000; + static constexpr int contiOxygenEqIdx = + numMICPs > 2 ? contiMicrobialEqIdx + 1 : -1000; //! Index of the continuity equation for the fourth MICP component - static constexpr int contiBiofilmEqIdx = - numMICPs > 3 ? contiUreaEqIdx + 1 : -1000; + static constexpr int contiUreaEqIdx = + numMICPs > 3 ? contiOxygenEqIdx + 1 : -1000; //! Index of the continuity equation for the fifth MICP component static constexpr int contiCalciteEqIdx = - numMICPs > 4 ? contiBiofilmEqIdx + 1 : -1000; + numMICPs > 4 ? contiUreaEqIdx + 1 : -1000; //! Index of the continuity equation for the foam component static constexpr int contiFoamEqIdx = diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index a55d31d9fa..074a460736 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -129,6 +129,8 @@ class BlackOilIntensiveQuantities using DirectionalMobilityPtr = Opm::Utility::CopyablePtr>; using BrineModule = BlackOilBrineModule; + using BrineIntQua = BlackOilBrineIntensiveQuantities; + using MICPIntQua = BlackOilMICPIntensiveQuantities; public: @@ -506,11 +508,11 @@ public: // deal with water induced rock compaction porosity_ *= problem.template rockCompPoroMultiplier(*this, globalSpaceIdx); - // the MICP processes change the porosity + // deal with MICP if constexpr (enableMICP){ - Evaluation biofilm_ = priVars.makeEvaluation(Indices::biofilmConcentrationIdx, timeIdx, linearizationType); - Evaluation calcite_ = priVars.makeEvaluation(Indices::calciteConcentrationIdx, timeIdx, linearizationType); - porosity_ += - biofilm_ - calcite_; + Evaluation biofilm_ = priVars.makeEvaluation(Indices::biofilmConcentrationIdx, timeIdx, linearizationType); + Evaluation calcite_ = priVars.makeEvaluation(Indices::calciteConcentrationIdx, timeIdx, linearizationType); + porosity_ += max(- biofilm_ - calcite_, -referencePorosity_+1e-8); } // deal with salt-precipitation @@ -692,6 +694,15 @@ public: Scalar referencePorosity() const { return referencePorosity_; } + const Evaluation& permFactor() const + { + if constexpr (enableMICP) + return MICPIntQua::permFactor(); + if constexpr (enableSaltPrecipitation) + return BrineIntQua::permFactor(); + throw std::logic_error("permFactor() called but salt precipitation or MICP are disabled"); + } + private: friend BlackOilSolventIntensiveQuantities; friend BlackOilExtboIntensiveQuantities; diff --git a/opm/models/blackoil/blackoillocalresidualtpfa.hh b/opm/models/blackoil/blackoillocalresidualtpfa.hh index 553e1374e9..6dfbb9d263 100644 --- a/opm/models/blackoil/blackoillocalresidualtpfa.hh +++ b/opm/models/blackoil/blackoillocalresidualtpfa.hh @@ -372,7 +372,10 @@ public: const IntensiveQuantities& up = (upIdx == interiorDofIdx) ? intQuantsIn : intQuantsEx; unsigned globalUpIndex = (upIdx == interiorDofIdx) ? globalIndexIn : globalIndexEx; // Use arithmetic average (more accurate with harmonic, but that requires recomputing the transmissbility) - const Evaluation transMult = (intQuantsIn.rockCompTransMultiplier() + Toolbox::value(intQuantsEx.rockCompTransMultiplier()))/2; + Evaluation transMult = (intQuantsIn.rockCompTransMultiplier() + Toolbox::value(intQuantsEx.rockCompTransMultiplier()))/2; + if constexpr (enableMICP) { + transMult *= (intQuantsIn.permFactor() + Toolbox::value(intQuantsEx.permFactor()))/2; + } Evaluation darcyFlux; if (globalUpIndex == globalIndexIn) { darcyFlux = pressureDifference * up.mobility(phaseIdx, facedir) * transMult * (-trans / faceArea); @@ -396,6 +399,10 @@ public: EnergyModule::template addPhaseEnthalpyFluxes_( flux, phaseIdx, darcyFlux, up.fluidState()); } + if constexpr (enableMICP) { + MICPModule::template addMICPFluxes_( + flux, darcyFlux, intQuantsIn); + } } else { const auto& invB = getInvB_(up.fluidState(), phaseIdx, pvtRegionIdx); const auto& surfaceVolumeFlux = invB * darcyFlux; @@ -406,6 +413,11 @@ public: addPhaseEnthalpyFluxes_ (flux,phaseIdx,darcyFlux, up.fluidState()); } + if constexpr (enableMICP) { + MICPModule::template + addMICPFluxes_ + (flux, darcyFlux, intQuantsEx); + } } } @@ -478,8 +490,8 @@ public: const Scalar diffusivity = nbInfo.diffusivity; const Scalar tmpdiffusivity = diffusivity / faceArea; DiffusionModule::addDiffusiveFlux(flux, - intQuantsIn.fluidState(), - intQuantsEx.fluidState(), + intQuantsIn, + intQuantsEx, tmpdiffusivity, effectiveDiffusionCoefficient); @@ -491,16 +503,12 @@ public: const Scalar dispersivity = nbInfo.dispersivity; const Scalar tmpdispersivity = dispersivity / faceArea; DispersionModule::addDispersiveFlux(flux, - intQuantsIn.fluidState(), - intQuantsEx.fluidState(), + intQuantsIn, + intQuantsEx, tmpdispersivity, normVelocityAvg); } - // deal with micp (if present) - static_assert(!enableMICP, "Relevant computeFlux() method must be implemented for this module before enabling."); - // MICPModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx); - } template @@ -612,7 +620,7 @@ public: // conductive heat flux from boundary if constexpr(enableEnergy){ Evaluation heatFlux; - // avoid overload of functions with same numeber of elements in eclproblem + // avoid overload of functions with same number of elements in eclproblem Scalar alpha = problem.eclTransmissibilities().thermalHalfTransBoundary(globalSpaceIdx, bdyInfo.boundaryFaceIndex); unsigned inIdx = 0;//dummy // always calculated with derivatives of this cell @@ -627,7 +635,7 @@ public: static_assert(!enableSolvent, "Relevant treatment of boundary conditions must be implemented before enabling."); static_assert(!enablePolymer, "Relevant treatment of boundary conditions must be implemented before enabling."); - static_assert(!enableMICP, "Relevant treatment of boundary conditions must be implemented before enabling."); + //static_assert(!enableMICP, "Relevant treatment of boundary conditions must be implemented before enabling."); // make sure that the right mass conservation quantities are used adaptMassConservationQuantities_(bdyFlux, insideIntQuants.pvtRegionIndex()); @@ -677,6 +685,7 @@ public: static void computeSource(RateVector& source, const Problem& problem, + const IntensiveQuantities& insideIntQuants, unsigned globalSpaceIdex, unsigned timeIdx) { @@ -685,9 +694,7 @@ public: problem.source(source, globalSpaceIdex, timeIdx); // deal with MICP (if present) - // deal with micp (if present) - static_assert(!enableMICP, "Relevant addSource() method must be implemented for this module before enabling."); - // MICPModule::addSource(source, elemCtx, dofIdx, timeIdx); + MICPModule::addSource(source, problem, insideIntQuants, globalSpaceIdex); // scale the source term of the energy equation if (enableEnergy) @@ -696,6 +703,7 @@ public: static void computeSourceDense(RateVector& source, const Problem& problem, + const IntensiveQuantities& insideIntQuants, unsigned globalSpaceIdex, unsigned timeIdx) { @@ -703,9 +711,7 @@ public: problem.addToSourceDense(source, globalSpaceIdex, timeIdx); // deal with MICP (if present) - // deal with micp (if present) - static_assert(!enableMICP, "Relevant addSource() method must be implemented for this module before enabling."); - // MICPModule::addSource(source, elemCtx, dofIdx, timeIdx); + MICPModule::addSource(source, problem, insideIntQuants, globalSpaceIdex); // scale the source term of the energy equation if (enableEnergy) diff --git a/opm/models/blackoil/blackoilmicpmodules.hh b/opm/models/blackoil/blackoilmicpmodules.hh index 507eef3520..3ec8c036f8 100644 --- a/opm/models/blackoil/blackoilmicpmodules.hh +++ b/opm/models/blackoil/blackoilmicpmodules.hh @@ -25,8 +25,8 @@ * * \brief Contains the classes required to extend the black-oil model by MICP. */ -#ifndef EWOMS_BLACK_OIL_MICP_MODULE_HH -#define EWOMS_BLACK_OIL_MICP_MODULE_HH +#ifndef OPM_BLACK_OIL_MICP_MODULE_HH +#define OPM_BLACK_OIL_MICP_MODULE_HH #include @@ -59,9 +59,14 @@ class BlackOilMICPModule using EqVector = GetPropType; using RateVector = GetPropType; using Indices = GetPropType; + using Problem = GetPropType; using Toolbox = MathToolbox; + using TabulatedFunction = typename BlackOilMICPParams::TabulatedFunction; + + enum { waterCompIdx = FluidSystem::waterCompIdx }; + static constexpr unsigned microbialConcentrationIdx = Indices::microbialConcentrationIdx; static constexpr unsigned oxygenConcentrationIdx = Indices::oxygenConcentrationIdx; static constexpr unsigned ureaConcentrationIdx = Indices::ureaConcentrationIdx; @@ -85,30 +90,13 @@ public: params_ = params; } - /*! - * \brief The simulator stops if "clogging" has been (almost) reached in any of the cells. - * - * I.e., porosity - biofilm - calcite < tol_clgg, where tol_clgg is a given tolerance. In the - * implemented model a permebaility-porosity relatonship is used where a minimum - * permeability value is reached if porosity - biofilm - calcite < phi_crit. - */ - static void checkCloggingMICP(const Model& model, const Scalar phi, unsigned dofIdx) - { - const PrimaryVariables& priVars = model.solution(/*timeIdx=*/1)[dofIdx]; - if (phi - priVars[biofilmConcentrationIdx] - priVars[calciteConcentrationIdx] < toleranceBeforeClogging()) - throw std::logic_error("Clogging has been (almost) reached in at least one cell\n"); - } - /*! * \brief Register all run-time parameters for the black-oil MICP module. */ static void registerParameters() { - if (!enableMICP) - // MICP has been disabled at compile time - return; - - VtkBlackOilMICPModule::registerParameters(); + if constexpr (enableMICP) + VtkBlackOilMICPModule::registerParameters(); } /*! @@ -117,20 +105,17 @@ public: static void registerOutputModules(Model& model, Simulator& simulator) { - if (!enableMICP) - // MICP has been disabled at compile time - return; - - model.addOutputModule(new VtkBlackOilMICPModule(simulator)); + if constexpr (enableMICP) + model.addOutputModule(new VtkBlackOilMICPModule(simulator)); } static bool eqApplies(unsigned eqIdx) { - if (!enableMICP) + if constexpr (enableMICP) + return eqIdx == contiMicrobialEqIdx || eqIdx == contiOxygenEqIdx || eqIdx == contiUreaEqIdx + || eqIdx == contiBiofilmEqIdx || eqIdx == contiCalciteEqIdx; + else return false; - - // All MICP components are true here - return eqIdx == contiMicrobialEqIdx || eqIdx == contiOxygenEqIdx || eqIdx == contiUreaEqIdx || eqIdx == contiBiofilmEqIdx || eqIdx == contiCalciteEqIdx; } static Scalar eqWeight([[maybe_unused]] unsigned eqIdx) @@ -146,208 +131,240 @@ public: static void addStorage(Dune::FieldVector& storage, const IntensiveQuantities& intQuants) { - if (!enableMICP) - return; - - LhsEval surfaceVolumeWater = Toolbox::template decay(intQuants.porosity()); - // avoid singular matrix if no water is present. - surfaceVolumeWater = max(surfaceVolumeWater, 1e-10); - // Suspended microbes in water phase - const LhsEval massMicrobes = surfaceVolumeWater * Toolbox::template decay(intQuants.microbialConcentration()); - LhsEval accumulationMicrobes = massMicrobes; - storage[contiMicrobialEqIdx] += accumulationMicrobes; - // Oxygen in water phase - const LhsEval massOxygen = surfaceVolumeWater * Toolbox::template decay(intQuants.oxygenConcentration()); - LhsEval accumulationOxygen = massOxygen; - storage[contiOxygenEqIdx] += accumulationOxygen; - // Urea in water phase - const LhsEval massUrea = surfaceVolumeWater * Toolbox::template decay(intQuants.ureaConcentration()); - LhsEval accumulationUrea = massUrea; - storage[contiUreaEqIdx] += accumulationUrea; - // Biofilm - const LhsEval massBiofilm = Toolbox::template decay(intQuants.biofilmConcentration()); - LhsEval accumulationBiofilm = massBiofilm; - storage[contiBiofilmEqIdx] += accumulationBiofilm; - // Calcite - const LhsEval massCalcite = Toolbox::template decay(intQuants.calciteConcentration()); - LhsEval accumulationCalcite = massCalcite; - storage[contiCalciteEqIdx] += accumulationCalcite; + if constexpr (enableMICP) { + const auto& fs = intQuants.fluidState(); + LhsEval surfaceVolumeWater = Toolbox::template decay(fs.invB(waterPhaseIdx)) + * Toolbox::template decay(intQuants.porosity()); + // avoid singular matrix if no water is present. + surfaceVolumeWater = max(surfaceVolumeWater, 1e-10); + // suspended microbes in water phase + const LhsEval massMicrobes = surfaceVolumeWater * Toolbox::template decay(intQuants.microbialConcentration()); + LhsEval accumulationMicrobes = massMicrobes; + storage[contiMicrobialEqIdx] += accumulationMicrobes; + // oxygen in water phase + const LhsEval massOxygen = surfaceVolumeWater * Toolbox::template decay(intQuants.oxygenConcentration()); + LhsEval accumulationOxygen = massOxygen; + storage[contiOxygenEqIdx] += accumulationOxygen; + // urea in water phase + const LhsEval massUrea = surfaceVolumeWater * Toolbox::template decay(intQuants.ureaConcentration()); + LhsEval accumulationUrea = massUrea; + storage[contiUreaEqIdx] += accumulationUrea; + // biofilm + const LhsEval massBiofilm = Toolbox::template decay(intQuants.biofilmConcentration()); + LhsEval accumulationBiofilm = massBiofilm; + storage[contiBiofilmEqIdx] += accumulationBiofilm; + // calcite + const LhsEval massCalcite = Toolbox::template decay(intQuants.calciteConcentration()); + LhsEval accumulationCalcite = massCalcite; + storage[contiCalciteEqIdx] += accumulationCalcite; + } } - static void computeFlux(RateVector& flux, - const ElementContext& elemCtx, - unsigned scvfIdx, - unsigned timeIdx) + template + static void addMICPFluxes_(RateVector& flux, + const Eval& volumeFlux, + const IntensiveQuantities& upFs) + { + if constexpr (enableMICP) { + flux[contiMicrobialEqIdx] += + decay(upFs.microbialConcentration()) + * volumeFlux; + flux[contiOxygenEqIdx] += + decay(upFs.oxygenConcentration()) + * volumeFlux; + flux[contiUreaEqIdx] += + decay(upFs.ureaConcentration()) + * volumeFlux; + } + } + + static void computeFlux([[maybe_unused]] RateVector& flux, + [[maybe_unused]] const ElementContext& elemCtx, + [[maybe_unused]] unsigned scvfIdx, + [[maybe_unused]] unsigned timeIdx) { - if (!enableMICP) - return; + if constexpr (enableMICP) { + flux[contiMicrobialEqIdx] = 0.0; + flux[contiOxygenEqIdx] = 0.0; + flux[contiUreaEqIdx] = 0.0; + const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx); + unsigned focusIdx = elemCtx.focusDofIndex(); + unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx); + if (upIdx == focusIdx) + addMICPFluxes_(flux, elemCtx, scvfIdx, timeIdx); + else + addMICPFluxes_(flux, elemCtx, scvfIdx, timeIdx); + } + } + template + static void addMICPFluxes_(RateVector& flux, + const ElementContext& elemCtx, + unsigned scvfIdx, + unsigned timeIdx) + { const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx); - - const unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx); - const unsigned inIdx = extQuants.interiorIndex(); + unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx); const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx); - - - if (upIdx == inIdx) { - flux[contiMicrobialEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * up.microbialConcentration(); - flux[contiOxygenEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * up.oxygenConcentration(); - flux[contiUreaEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * up.ureaConcentration(); - } - else { - flux[contiMicrobialEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * decay(up.microbialConcentration()); - flux[contiOxygenEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * decay(up.oxygenConcentration()); - flux[contiUreaEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * decay(up.ureaConcentration()); - } + const auto& fs = up.fluidState(); + const auto& volFlux = extQuants.volumeFlux(waterPhaseIdx); + addMICPFluxes_(flux, volFlux, fs); } - // See https://doi.org/10.1016/j.ijggc.2021.103256 for the micp processes in the model. + // see https://doi.org/10.1016/j.ijggc.2021.103256 for the MICP processes in the model static void addSource(RateVector& source, - const ElementContext& elemCtx, - unsigned dofIdx, - unsigned timeIdx) - + const Problem& problem, + const IntensiveQuantities& intQuants, + unsigned globalSpaceIdex) { - if (!enableMICP) - return; + if constexpr (enableMICP) { + const auto& velocityInf = problem.model().linearizer().getVelocityInfo(); + auto velocityInfos = velocityInf[globalSpaceIdex]; + Scalar normVelocityCell = 0.0; + for (const auto& velocityInfo : velocityInfos) { + normVelocityCell = max(normVelocityCell, std::abs(velocityInfo.velocity[waterPhaseIdx])); + } - // compute dpW (max norm of the pressure gradient in the cell center) - const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx); - const auto& K = elemCtx.problem().intrinsicPermeability(elemCtx, dofIdx, 0); - size_t numInteriorFaces = elemCtx.numInteriorFaces(timeIdx); - Evaluation dpW = 0; - for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) { - const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx); - unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx); - const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx); - const Evaluation& mobWater = up.mobility(waterPhaseIdx); + // get the model parameters + const auto b = intQuants.fluidState().invB(waterPhaseIdx); + unsigned satnumIdx = problem.satnumRegionIndex(globalSpaceIdex); + Scalar k_a = microbialAttachmentRate(satnumIdx); + Scalar k_d = microbialDeathRate(satnumIdx); + Scalar rho_b = densityBiofilm(satnumIdx); + Scalar rho_c = densityCalcite(satnumIdx); + Scalar k_str = detachmentRate(satnumIdx); + Scalar eta = detachmentExponent(satnumIdx); + Scalar k_o = halfVelocityOxygen(satnumIdx); + Scalar k_u = halfVelocityUrea(satnumIdx) / 10;// dividing by scaling factor 10 (see WellInterfaceGeneric.cpp) + Scalar mu = maximumGrowthRate(satnumIdx); + Scalar mu_u = maximumUreaUtilization(satnumIdx) / 10;// dividing by scaling factor 10 + Scalar Y_sb = yieldGrowthCoefficient(satnumIdx); + Scalar F = oxygenConsumptionFactor(satnumIdx); + Scalar Y_uc = yieldUreaToCalciteCoefficient(satnumIdx) * 10;// multiplying by scaling factor 10 - // compute water velocity from flux - Evaluation waterVolumeVelocity = extQuants.volumeFlux(waterPhaseIdx) / (K[0][0] * mobWater); - dpW = std::max(dpW, abs(waterVolumeVelocity)); + // compute the processes + source[Indices::contiMicrobialEqIdx] += intQuants.microbialConcentration() * intQuants.porosity() * b * + (Y_sb * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()) - k_d - k_a) + + rho_b * intQuants.biofilmConcentration() * k_str * pow(normVelocityCell, eta); + + source[Indices::contiOxygenEqIdx] -= (intQuants.microbialConcentration() * intQuants.porosity() * b + rho_b * intQuants.biofilmConcentration()) * + F * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()); + + source[Indices::contiUreaEqIdx] -= rho_b * intQuants.biofilmConcentration() * mu_u * intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration()); + + source[Indices::contiBiofilmEqIdx] += intQuants.biofilmConcentration() * (Y_sb * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()) - k_d + - k_str * pow(normVelocityCell, eta) - Y_uc * (rho_b / rho_c) * intQuants.biofilmConcentration() * mu_u * + (intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration())) / (intQuants.porosity() + intQuants.biofilmConcentration())) + + k_a * intQuants.microbialConcentration() * intQuants.porosity() * b / rho_b; + + source[Indices::contiCalciteEqIdx] += (rho_b / rho_c) * intQuants.biofilmConcentration() * Y_uc * mu_u * intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration()); } - - // get the model parameters - Scalar k_a = microbialAttachmentRate(); - Scalar k_d = microbialDeathRate(); - Scalar rho_b = densityBiofilm(); - Scalar rho_c = densityCalcite(); - Scalar k_str = detachmentRate(); - Scalar k_o = halfVelocityOxygen(); - Scalar k_u = halfVelocityUrea() / 10.0;//Dividing by scaling factor 10 (see WellInterface_impl.hpp) - Scalar mu = maximumGrowthRate(); - Scalar mu_u = maximumUreaUtilization() / 10.0;//Dividing by scaling factor 10 (see WellInterface_impl.hpp) - Scalar Y_sb = yieldGrowthCoefficient(); - Scalar F = oxygenConsumptionFactor(); - Scalar Y_uc = 1.67 * 10; //Multiplying by scaling factor 10 (see WellInterface_impl.hpp) - - // compute the processes - source[Indices::contiMicrobialEqIdx] += intQuants.microbialConcentration() * intQuants.porosity() * - (Y_sb * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()) - k_d - k_a) - + rho_b * intQuants.biofilmConcentration() * k_str * pow(intQuants.porosity() * dpW, 0.58); - - source[Indices::contiOxygenEqIdx] -= (intQuants.microbialConcentration() * intQuants.porosity() + rho_b * intQuants.biofilmConcentration()) * - F * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()); - - source[Indices::contiUreaEqIdx] -= rho_b * intQuants.biofilmConcentration() * mu_u * intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration()); - - source[Indices::contiBiofilmEqIdx] += intQuants.biofilmConcentration() * (Y_sb * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()) - k_d - - k_str * pow(intQuants.porosity() * dpW, 0.58) - Y_uc * (rho_b / rho_c) * intQuants.biofilmConcentration() * mu_u * - (intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration())) / (intQuants.porosity() + intQuants.biofilmConcentration())) - + k_a * intQuants.microbialConcentration() * intQuants.porosity() / rho_b; - - source[Indices::contiCalciteEqIdx] += (rho_b / rho_c) * intQuants.biofilmConcentration() * Y_uc * mu_u * intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration()); } - static const Scalar densityBiofilm() + static void addSource([[maybe_unused]] RateVector& source, + [[maybe_unused]] const ElementContext& elemCtx, + [[maybe_unused]] unsigned dofIdx, + [[maybe_unused]] unsigned timeIdx) { - return params_.densityBiofilm_; + if constexpr (enableMICP) { + const auto& problem = elemCtx.problem(); + const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx); + addSource(source, problem, intQuants, dofIdx); + } } - static const Scalar densityCalcite() + static const Scalar densityBiofilm(unsigned satnumRegionIdx) { - return params_.densityCalcite_; + return params_.densityBiofilm_[satnumRegionIdx]; } - static const Scalar detachmentRate() + static const Scalar densityCalcite(unsigned satnumRegionIdx) { - return params_.detachmentRate_; + return params_.densityCalcite_[satnumRegionIdx]; } - static const Scalar criticalPorosity() + static const Scalar detachmentRate(unsigned satnumRegionIdx) { - return params_.criticalPorosity_; + return params_.detachmentRate_[satnumRegionIdx]; } - static const Scalar fittingFactor() + static const Scalar detachmentExponent(unsigned satnumRegionIdx) { - return params_.fittingFactor_; + return params_.detachmentExponent_[satnumRegionIdx]; } - static const Scalar halfVelocityOxygen() + static const Scalar halfVelocityOxygen(unsigned satnumRegionIdx) { - return params_.halfVelocityOxygen_; + return params_.halfVelocityOxygen_[satnumRegionIdx]; } - static const Scalar halfVelocityUrea() + static const Scalar halfVelocityUrea(unsigned satnumRegionIdx) { - return params_.halfVelocityUrea_; + return params_.halfVelocityUrea_[satnumRegionIdx]; } - static const Scalar maximumGrowthRate() + static const Scalar maximumGrowthRate(unsigned satnumRegionIdx) { - return params_.maximumGrowthRate_; + return params_.maximumGrowthRate_[satnumRegionIdx]; } - static const Scalar maximumOxygenConcentration() + static const Scalar maximumUreaUtilization(unsigned satnumRegionIdx) { - return params_.maximumOxygenConcentration_; + return params_.maximumUreaUtilization_[satnumRegionIdx]; } - static const Scalar maximumUreaConcentration() + static const Scalar microbialAttachmentRate(unsigned satnumRegionIdx) { - return params_.maximumUreaConcentration_ / 10.0;//Dividing by scaling factor 10 (see WellInterface_impl.hpp); + return params_.microbialAttachmentRate_[satnumRegionIdx]; } - static const Scalar maximumUreaUtilization() + static const Scalar microbialDeathRate(unsigned satnumRegionIdx) { - return params_.maximumUreaUtilization_; + return params_.microbialDeathRate_[satnumRegionIdx]; } - static const Scalar microbialAttachmentRate() + static const Scalar oxygenConsumptionFactor(unsigned satnumRegionIdx) { - return params_.microbialAttachmentRate_; + return params_.oxygenConsumptionFactor_[satnumRegionIdx]; } - static const Scalar microbialDeathRate() + static const Scalar yieldGrowthCoefficient(unsigned satnumRegionIdx) { - return params_.microbialDeathRate_; + return params_.yieldGrowthCoefficient_[satnumRegionIdx]; } - static const Scalar minimumPermeability() + static const Scalar yieldUreaToCalciteCoefficient(unsigned satnumRegionIdx) { - return params_.minimumPermeability_; + return params_.yieldUreaToCalciteCoefficient_[satnumRegionIdx]; } - static const Scalar oxygenConsumptionFactor() + static const Scalar microbialDiffusion(unsigned pvtRegionIdx) { - return params_.oxygenConsumptionFactor_; + return params_.microbialDiffusion_[pvtRegionIdx]; } - static const Scalar toleranceBeforeClogging() + static const Scalar oxygenDiffusion(unsigned pvtRegionIdx) { - return params_.toleranceBeforeClogging_; + return params_.oxygenDiffusion_[pvtRegionIdx]; } - static const Scalar yieldGrowthCoefficient() + static const Scalar ureaDiffusion(unsigned pvtRegionIdx) { - return params_.yieldGrowthCoefficient_; + return params_.ureaDiffusion_[pvtRegionIdx]; } - static const std::vector phi() + static const TabulatedFunction& permfactTable(const ElementContext& elemCtx, + unsigned scvIdx, + unsigned timeIdx) { - return params_.phi_; + unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx); + return params_.permfactTable_[satnumRegionIdx]; + } + + static const TabulatedFunction& permfactTable(unsigned satnumRegionIdx) + { + return params_.permfactTable_[satnumRegionIdx]; } private: @@ -396,17 +413,12 @@ public: * */ void MICPPropertiesUpdate_(const ElementContext& elemCtx, - unsigned dofIdx, - unsigned timeIdx) + unsigned dofIdx, + unsigned timeIdx) { const auto linearizationType = elemCtx.linearizationType(); const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx); - const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx); - const auto& K = elemCtx.problem().intrinsicPermeability(elemCtx, dofIdx, timeIdx); - Scalar referencePorosity_ = elemCtx.problem().porosity(elemCtx, dofIdx, timeIdx); - Scalar eta = MICPModule::fittingFactor(); - Scalar k_min = MICPModule::minimumPermeability(); - Scalar phi_crit = MICPModule::criticalPorosity(); + const Scalar referencePorosity_ = elemCtx.problem().referencePorosity(dofIdx, timeIdx); microbialConcentration_ = priVars.makeEvaluation(microbialConcentrationIdx, timeIdx, linearizationType); oxygenConcentration_ = priVars.makeEvaluation(oxygenConcentrationIdx, timeIdx, linearizationType); @@ -414,9 +426,13 @@ public: biofilmConcentration_ = priVars.makeEvaluation(biofilmConcentrationIdx, timeIdx, linearizationType); calciteConcentration_ = priVars.makeEvaluation(calciteConcentrationIdx, timeIdx, linearizationType); - // Permeability reduction due to MICP, by adjusting the water mobility - asImp_().mobility_[waterPhaseIdx] *= max((pow((intQuants.porosity() - phi_crit) / (referencePorosity_ - phi_crit), eta) + k_min / K[0][0])/(1. + k_min / K[0][0]), k_min / K[0][0]); + const Evaluation porosityFactor = min(1.0 - (biofilmConcentration_+calciteConcentration_)/(referencePorosity_+1e-8), 1.0); //phi/phi_0 + unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, dofIdx, timeIdx); + const auto& permfactTable = MICPModule::permfactTable(satnumRegionIdx); + permFactor_ = permfactTable.eval(porosityFactor, /*extrapolation=*/true); + biofilmMass_ = referencePorosity_*biofilmConcentration_*MICPModule::densityBiofilm(satnumRegionIdx); + calciteMass_ = referencePorosity_*calciteConcentration_*MICPModule::densityCalcite(satnumRegionIdx); } const Evaluation& microbialConcentration() const @@ -434,16 +450,24 @@ public: const Evaluation& calciteConcentration() const { return calciteConcentration_; } + const Evaluation biofilmMass() const + { return biofilmMass_; } + + const Evaluation calciteMass() const + { return calciteMass_; } + + const Evaluation& permFactor() const + { return permFactor_; } protected: - Implementation& asImp_() - { return *static_cast(this); } - Evaluation microbialConcentration_; Evaluation oxygenConcentration_; Evaluation ureaConcentration_; Evaluation biofilmConcentration_; Evaluation calciteConcentration_; + Evaluation biofilmMass_; + Evaluation calciteMass_; + Evaluation permFactor_; }; @@ -456,8 +480,8 @@ class BlackOilMICPIntensiveQuantities public: void MICPPropertiesUpdate_(const ElementContext&, - unsigned, - unsigned) + unsigned, + unsigned) { } const Evaluation& microbialConcentration() const @@ -474,6 +498,15 @@ public: const Evaluation& calciteConcentration() const { throw std::logic_error("calciteConcentration() called but MICP is disabled"); } + + const Evaluation& biofilmMass() const + { throw std::logic_error("biofilmMass() called but MICP is disabled"); } + + const Evaluation& calciteMass() const + { throw std::logic_error("calciteMass() called but MICP is disabled"); } + + const Evaluation& permFactor() const + { throw std::logic_error("permFactor() called but MICP is disabled"); } }; /*! @@ -487,11 +520,6 @@ template ; - -private: - Implementation& asImp_() - { return *static_cast(this); } - }; template diff --git a/opm/models/blackoil/blackoilmicpparams.cpp b/opm/models/blackoil/blackoilmicpparams.cpp index 9afaa55edc..1ac9e937e2 100644 --- a/opm/models/blackoil/blackoilmicpparams.cpp +++ b/opm/models/blackoil/blackoilmicpparams.cpp @@ -26,7 +26,11 @@ #if HAVE_ECL_INPUT #include -#include +#include +#include +#include +#include +#include #endif #include @@ -59,33 +63,71 @@ initFromState(const EclipseState& eclState) if (!eclState.runspec().micp()) return; // MICP treatment is supposed to be disabled*/ - // initialize the objects which deal with the MICPpara keyword - const auto& MICPpara = eclState.getMICPpara(); - densityBiofilm_ = MICPpara.getDensityBiofilm(); - densityCalcite_ = MICPpara.getDensityCalcite(); - detachmentRate_ = MICPpara.getDetachmentRate(); - criticalPorosity_ = MICPpara.getCriticalPorosity(); - fittingFactor_ = MICPpara.getFittingFactor(); - halfVelocityOxygen_ = MICPpara.getHalfVelocityOxygen(); - halfVelocityUrea_ = MICPpara.getHalfVelocityUrea(); - maximumGrowthRate_ = MICPpara.getMaximumGrowthRate(); - maximumUreaUtilization_ = MICPpara.getMaximumUreaUtilization(); - microbialAttachmentRate_ = MICPpara.getMicrobialAttachmentRate(); - microbialDeathRate_ = MICPpara.getMicrobialDeathRate(); - minimumPermeability_ = MICPpara.getMinimumPermeability(); - oxygenConsumptionFactor_ = MICPpara.getOxygenConsumptionFactor(); - yieldGrowthCoefficient_ = MICPpara.getYieldGrowthCoefficient(); - maximumOxygenConcentration_ = MICPpara.getMaximumOxygenConcentration(); - maximumUreaConcentration_ = MICPpara.getMaximumUreaConcentration(); - toleranceBeforeClogging_ = MICPpara.getToleranceBeforeClogging(); + const auto& tableManager = eclState.getTableManager(); + unsigned numSatRegions = tableManager.getTabdims().getNumSatTables(); + unsigned numPvtRegions = tableManager.getTabdims().getNumPVTTables(); - // obtain the porosity for the clamp in the blackoilnewtonmethod - if constexpr (std::is_same_v) { - const auto phi = eclState.fieldProps().get_double("PORO"); - phi_.resize(phi.size()); - std::copy(phi.begin(), phi.end(), phi_.begin()); - } else { - phi_ = eclState.fieldProps().get_double("PORO"); + // initialize the objects which deal with the MICP parameters + const TableContainer& biofilmTables = tableManager.getBiofilmTables(); + if (biofilmTables.empty()) { + throw std::runtime_error("MICP requires the BIOFPARA keyword"); + } + densityBiofilm_.resize(numSatRegions); + densityCalcite_.resize(numSatRegions); + detachmentRate_.resize(numSatRegions); + detachmentExponent_.resize(numSatRegions); + halfVelocityOxygen_.resize(numSatRegions); + halfVelocityUrea_.resize(numSatRegions); + maximumGrowthRate_.resize(numSatRegions); + maximumUreaUtilization_.resize(numSatRegions); + microbialAttachmentRate_.resize(numSatRegions); + microbialDeathRate_.resize(numSatRegions); + oxygenConsumptionFactor_.resize(numSatRegions); + yieldGrowthCoefficient_.resize(numSatRegions); + yieldUreaToCalciteCoefficient_.resize(numSatRegions); + for (unsigned stnRegionIdx = 0; stnRegionIdx < numSatRegions; ++stnRegionIdx) { + const BiofilmTable& biofilmTable = biofilmTables.getTable(stnRegionIdx); + densityBiofilm_[stnRegionIdx] = biofilmTable.getDensityBiofilm().front(); + densityCalcite_[stnRegionIdx] = biofilmTable.getDensityCalcite().front(); + detachmentRate_[stnRegionIdx] = biofilmTable.getDetachmentRate().front(); + detachmentExponent_[stnRegionIdx] = biofilmTable.getDetachmentExponent().front(); + halfVelocityOxygen_[stnRegionIdx] = biofilmTable.getHalfVelocityOxygen().front(); + halfVelocityUrea_[stnRegionIdx] = biofilmTable.getHalfVelocityUrea().front(); + maximumGrowthRate_[stnRegionIdx] = biofilmTable.getMaximumGrowthRate().front(); + maximumUreaUtilization_[stnRegionIdx] = biofilmTable.getMaximumUreaUtilization().front(); + microbialAttachmentRate_[stnRegionIdx] = biofilmTable.getMicrobialAttachmentRate().front(); + microbialDeathRate_[stnRegionIdx] = biofilmTable.getMicrobialDeathRate().front(); + oxygenConsumptionFactor_[stnRegionIdx] = biofilmTable.getOxygenConsumptionFactor().front(); + yieldGrowthCoefficient_[stnRegionIdx] = biofilmTable.getYieldGrowthCoefficient().front(); + yieldUreaToCalciteCoefficient_[stnRegionIdx] = biofilmTable.getYieldUreaToCalciteCoefficient().front(); + } + + const TableContainer& diffMICPTables = tableManager.getDiffMICPTables(); + microbialDiffusion_.resize(numPvtRegions); + oxygenDiffusion_.resize(numPvtRegions); + ureaDiffusion_.resize(numPvtRegions); + for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++pvtRegionIdx) { + if (!diffMICPTables.empty()) { + const DiffMICPTable& diffMICPTable = diffMICPTables.getTable(pvtRegionIdx); + microbialDiffusion_[pvtRegionIdx] = diffMICPTable.getMicrobialDiffusion().front(); + oxygenDiffusion_[pvtRegionIdx] = diffMICPTable.getOxygenDiffusion().front(); + ureaDiffusion_[pvtRegionIdx] = diffMICPTable.getUreaDiffusion().front(); + } + else { + microbialDiffusion_[pvtRegionIdx] = 0.0; + oxygenDiffusion_[pvtRegionIdx] = 0.0; + ureaDiffusion_[pvtRegionIdx] = 0.0; + } + } + + const TableContainer& permfactTables = tableManager.getPermfactTables(); + if (permfactTables.empty()) { + throw std::runtime_error("MICP requires the PERMFACT keyword"); + } + permfactTable_.resize(numSatRegions); + for (std::size_t i = 0; i < permfactTables.size(); ++i) { + const PermfactTable& permfactTable = permfactTables.getTable(i); + permfactTable_[i].setXYContainers(permfactTable.getPorosityChangeColumn(), permfactTable.getPermeabilityMultiplierColumn()); } } #endif diff --git a/opm/models/blackoil/blackoilmicpparams.hpp b/opm/models/blackoil/blackoilmicpparams.hpp index 13dcc67b5c..02a8c8643f 100644 --- a/opm/models/blackoil/blackoilmicpparams.hpp +++ b/opm/models/blackoil/blackoilmicpparams.hpp @@ -28,6 +28,8 @@ #ifndef OPM_BLACK_OIL_MICP_PARAMS_HPP #define OPM_BLACK_OIL_MICP_PARAMS_HPP +#include + #include namespace Opm { @@ -45,24 +47,26 @@ struct BlackOilMICPParams void initFromState(const EclipseState& eclState); #endif - Scalar densityBiofilm_; - Scalar densityCalcite_; - Scalar detachmentRate_; - Scalar criticalPorosity_; - Scalar fittingFactor_; - Scalar halfVelocityOxygen_; - Scalar halfVelocityUrea_; - Scalar maximumGrowthRate_; - Scalar maximumUreaUtilization_; - Scalar microbialAttachmentRate_; - Scalar microbialDeathRate_; - Scalar minimumPermeability_; - Scalar oxygenConsumptionFactor_; - Scalar yieldGrowthCoefficient_; - Scalar maximumOxygenConcentration_; - Scalar maximumUreaConcentration_; - Scalar toleranceBeforeClogging_; - std::vector phi_; + using TabulatedFunction = Tabulated1DFunction; + + std::vector permfactTable_; + + std::vector densityBiofilm_; + std::vector densityCalcite_; + std::vector detachmentRate_; + std::vector detachmentExponent_; + std::vector halfVelocityOxygen_; + std::vector halfVelocityUrea_; + std::vector maximumGrowthRate_; + std::vector maximumUreaUtilization_; + std::vector microbialAttachmentRate_; + std::vector microbialDeathRate_; + std::vector oxygenConsumptionFactor_; + std::vector yieldGrowthCoefficient_; + std::vector yieldUreaToCalciteCoefficient_; + std::vector microbialDiffusion_; + std::vector oxygenDiffusion_; + std::vector ureaDiffusion_; }; } // namespace Opm diff --git a/opm/models/blackoil/blackoilnewtonmethod.hpp b/opm/models/blackoil/blackoilnewtonmethod.hpp index cbbc1bdbd2..19bac6eb4e 100644 --- a/opm/models/blackoil/blackoilnewtonmethod.hpp +++ b/opm/models/blackoil/blackoilnewtonmethod.hpp @@ -373,25 +373,24 @@ protected: } - // Limit the variables to [0, cmax] values to improve the convergence. - // For the microorganisms we set this value equal to the biomass density value. - // For the oxygen and urea we set this value to the maximum injected - // concentration (the urea concentration has been scaled by 10). For - // the biofilm and calcite, we set this value equal to the porosity minus the clogging tolerance. - if (enableMICP && pvIdx == Indices::microbialConcentrationIdx) { - nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::densityBiofilm()); - } - if (enableMICP && pvIdx == Indices::oxygenConcentrationIdx) { - nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::maximumOxygenConcentration()); - } - if (enableMICP && pvIdx == Indices::ureaConcentrationIdx) { - nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::maximumUreaConcentration()); - } - if (enableMICP && pvIdx == Indices::biofilmConcentrationIdx) { - nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::phi()[globalDofIdx] - MICPModule::toleranceBeforeClogging()); - } - if (enableMICP && pvIdx == Indices::calciteConcentrationIdx) { - nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::phi()[globalDofIdx] - MICPModule::toleranceBeforeClogging()); + // keep the values above 0 + // for the biofilm and calcite, we set an upper limit + if constexpr (enableMICP) { + if (pvIdx == Indices::microbialConcentrationIdx) { + nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0}); + } + if (pvIdx == Indices::oxygenConcentrationIdx) { + nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0}); + } + if (pvIdx == Indices::ureaConcentrationIdx) { + nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0}); + } + if (pvIdx == Indices::biofilmConcentrationIdx) { + nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, this->problem().referencePorosity(globalDofIdx, 0) - 1.e-8); + } + if (pvIdx == Indices::calciteConcentrationIdx) { + nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, this->problem().referencePorosity(globalDofIdx, 0) - 1.e-8); + } } } diff --git a/opm/models/blackoil/blackoilonephaseindices.hh b/opm/models/blackoil/blackoilonephaseindices.hh index 863c8b00ea..9305fdc4b7 100644 --- a/opm/models/blackoil/blackoilonephaseindices.hh +++ b/opm/models/blackoil/blackoilonephaseindices.hh @@ -145,24 +145,24 @@ struct BlackOilOnePhaseIndices numPolymers > 1 ? polymerConcentrationIdx + 1 : -1000; //! Index of the primary variable for the first MICP component - static constexpr int microbialConcentrationIdx = + static constexpr int biofilmConcentrationIdx = enableMICP ? PVOffset + numPhases + numSolvents : -1000; //! Index of the primary variable for the second MICP component - static constexpr int oxygenConcentrationIdx = - numMICPs > 1 ? microbialConcentrationIdx + 1 : -1000; + static constexpr int microbialConcentrationIdx = + numMICPs > 1 ? biofilmConcentrationIdx + 1 : -1000; //! Index of the primary variable for the third MICP component - static constexpr int ureaConcentrationIdx = - numMICPs > 2 ? oxygenConcentrationIdx + 1 : -1000; + static constexpr int oxygenConcentrationIdx = + numMICPs > 2 ? microbialConcentrationIdx + 1 : -1000; //! Index of the primary variable for the fourth MICP component - static constexpr int biofilmConcentrationIdx = - numMICPs > 3 ? ureaConcentrationIdx + 1 : -1000; + static constexpr int ureaConcentrationIdx = + numMICPs > 3 ? oxygenConcentrationIdx + 1 : -1000; //! Index of the primary variable for the fifth MICP component static constexpr int calciteConcentrationIdx = - numMICPs > 4 ? biofilmConcentrationIdx + 1 : -1000; + numMICPs > 4 ? ureaConcentrationIdx + 1 : -1000; //! Index of the primary variable for the foam static constexpr int foamConcentrationIdx = @@ -221,24 +221,24 @@ struct BlackOilOnePhaseIndices numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000; //! Index of the continuity equation for the first MICP component - static constexpr int contiMicrobialEqIdx = + static constexpr int contiBiofilmEqIdx = enableMICP ? PVOffset + numPhases + numSolvents : -1000; //! Index of the continuity equation for the second MICP component - static constexpr int contiOxygenEqIdx = - numMICPs > 1 ? contiMicrobialEqIdx + 1 : -1000; + static constexpr int contiMicrobialEqIdx = + numMICPs > 1 ? contiBiofilmEqIdx + 1 : -1000; //! Index of the continuity equation for the third MICP component - static constexpr int contiUreaEqIdx = - numMICPs > 2 ? contiOxygenEqIdx + 1 : -1000; + static constexpr int contiOxygenEqIdx = + numMICPs > 2 ? contiMicrobialEqIdx + 1 : -1000; //! Index of the continuity equation for the fourth MICP component - static constexpr int contiBiofilmEqIdx = - numMICPs > 3 ? contiUreaEqIdx + 1 : -1000; + static constexpr int contiUreaEqIdx = + numMICPs > 3 ? contiOxygenEqIdx + 1 : -1000; //! Index of the continuity equation for the fifth MICP component static constexpr int contiCalciteEqIdx = - numMICPs > 4 ? contiBiofilmEqIdx + 1 : -1000; + numMICPs > 4 ? contiUreaEqIdx + 1 : -1000; //! Index of the continuity equation for the foam component static constexpr int contiFoamEqIdx = diff --git a/opm/models/blackoil/blackoiltwophaseindices.hh b/opm/models/blackoil/blackoiltwophaseindices.hh index 9d112a1488..27d1630bb6 100644 --- a/opm/models/blackoil/blackoiltwophaseindices.hh +++ b/opm/models/blackoil/blackoiltwophaseindices.hh @@ -145,24 +145,24 @@ struct BlackOilTwoPhaseIndices numPolymers > 1 ? polymerConcentrationIdx + 1 : -1000; //! Index of the primary variable for the first MICP component - static constexpr int microbialConcentrationIdx = + static constexpr int biofilmConcentrationIdx = enableMICP ? PVOffset + numPhases + numSolvents : -1000; //! Index of the primary variable for the second MICP component - static constexpr int oxygenConcentrationIdx = - numMICPs > 1 ? microbialConcentrationIdx + 1 : -1000; + static constexpr int microbialConcentrationIdx = + numMICPs > 1 ? biofilmConcentrationIdx + 1 : -1000; //! Index of the primary variable for the third MICP component - static constexpr int ureaConcentrationIdx = - numMICPs > 2 ? oxygenConcentrationIdx + 1 : -1000; + static constexpr int oxygenConcentrationIdx = + numMICPs > 2 ? microbialConcentrationIdx + 1 : -1000; //! Index of the primary variable for the fourth MICP component - static constexpr int biofilmConcentrationIdx = - numMICPs > 3 ? ureaConcentrationIdx + 1 : -1000; + static constexpr int ureaConcentrationIdx = + numMICPs > 3 ? oxygenConcentrationIdx + 1 : -1000; //! Index of the primary variable for the fifth MICP component static constexpr int calciteConcentrationIdx = - numMICPs > 4 ? biofilmConcentrationIdx + 1 : -1000; + numMICPs > 4 ? ureaConcentrationIdx + 1 : -1000; //! Index of the primary variable for the foam static constexpr int foamConcentrationIdx = @@ -240,24 +240,24 @@ struct BlackOilTwoPhaseIndices numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000; //! Index of the continuity equation for the first MICP component - static constexpr int contiMicrobialEqIdx = + static constexpr int contiBiofilmEqIdx = enableMICP ? PVOffset + numPhases + numSolvents : -1000; //! Index of the continuity equation for the second MICP component - static constexpr int contiOxygenEqIdx = - numMICPs > 1 ? contiMicrobialEqIdx + 1 : -1000; + static constexpr int contiMicrobialEqIdx = + numMICPs > 1 ? contiBiofilmEqIdx + 1 : -1000; //! Index of the continuity equation for the third MICP component - static constexpr int contiUreaEqIdx = - numMICPs > 2 ? contiOxygenEqIdx + 1 : -1000; + static constexpr int contiOxygenEqIdx = + numMICPs > 2 ? contiMicrobialEqIdx + 1 : -1000; //! Index of the continuity equation for the fourth MICP component - static constexpr int contiBiofilmEqIdx = - numMICPs > 3 ? contiUreaEqIdx + 1 : -1000; + static constexpr int contiUreaEqIdx = + numMICPs > 3 ? contiOxygenEqIdx + 1 : -1000; //! Index of the continuity equation for the fifth MICP component static constexpr int contiCalciteEqIdx = - numMICPs > 4 ? contiBiofilmEqIdx + 1 : -1000; + numMICPs > 4 ? contiUreaEqIdx + 1 : -1000; //! Index of the continuity equation for the foam component static constexpr int contiFoamEqIdx = diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index 54a58b83f5..9116df9e86 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -108,6 +108,7 @@ class TpfaLinearizer static const bool linearizeNonLocalElements = getPropValue(); static const bool enableEnergy = getPropValue(); static const bool enableDiffusion = getPropValue(); + static const bool enableMICP = getPropValue(); // copying the linearizer is not a good idea TpfaLinearizer(const TpfaLinearizer&) = delete; @@ -543,7 +544,7 @@ private: simulator_().problem().eclWriter().outputModule().getFlows().hasBlockFlows(); const bool anyFlores = simulator_().problem().eclWriter().outputModule().getFlows().anyFlores(); const bool enableDispersion = simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion(); - if (((!anyFlows || !flowsInfo_.empty()) && (!anyFlores || !floresInfo_.empty())) && !enableDispersion) { + if (((!anyFlows || !flowsInfo_.empty()) && (!anyFlores || !floresInfo_.empty())) && (!enableDispersion && !enableMICP)) { return; } const auto& model = model_(); @@ -569,7 +570,7 @@ private: if (anyFlores) { floresInfo_.reserve(numCells, 6 * numCells); } - if (enableDispersion) { + if (enableDispersion || enableMICP) { velocityInfo_.reserve(numCells, 6 * numCells); } @@ -616,7 +617,7 @@ private: if (anyFlores) { floresInfo_.appendRow(loc_flinfo.begin(), loc_flinfo.end()); } - if (enableDispersion) { + if (enableDispersion || enableMICP) { velocityInfo_.appendRow(loc_vlinfo.begin(), loc_vlinfo.end()); } } @@ -757,7 +758,7 @@ private: const IntensiveQuantities& intQuantsEx = model_().intensiveQuantities(globJ, /*timeIdx*/ 0); LocalResidual::computeFlux(adres,darcyFlux, globI, globJ, intQuantsIn, intQuantsEx, nbInfo.res_nbinfo, problem_().moduleParams()); adres *= nbInfo.res_nbinfo.faceArea; - if (enableDispersion) { + if (enableDispersion || enableMICP) { for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) { velocityInfo_[globI][loc].velocity[phaseIdx] = darcyFlux[phaseIdx].value() / nbInfo.res_nbinfo.faceArea; } @@ -831,9 +832,9 @@ private: bMat = 0.0; adres = 0.0; if (separateSparseSourceTerms_) { - LocalResidual::computeSourceDense(adres, problem_(), globI, 0); + LocalResidual::computeSourceDense(adres, problem_(), intQuantsIn, globI, 0); } else { - LocalResidual::computeSource(adres, problem_(), globI, 0); + LocalResidual::computeSource(adres, problem_(), intQuantsIn, globI, 0); } adres *= -volume; setResAndJacobi(res, bMat, adres); diff --git a/opm/models/io/vtkblackoilmicpmodule.hpp b/opm/models/io/vtkblackoilmicpmodule.hpp index 8f492c90de..4e791b8cb7 100644 --- a/opm/models/io/vtkblackoilmicpmodule.hpp +++ b/opm/models/io/vtkblackoilmicpmodule.hpp @@ -141,8 +141,9 @@ public: } if (params_.ureaConcentrationOutput_) { + // multypliging by scaling factor 10 (see WellInterfaceGeneric.cpp) ureaConcentration_[globalDofIdx] = - 10 * scalarValue(intQuants.ureaConcentration());//Multypliging by scaling factor 10 (see WellInterface_impl.hpp) + 10 * scalarValue(intQuants.ureaConcentration()); } if (params_.biofilmConcentrationOutput_) { diff --git a/opm/simulators/flow/FIPContainer.cpp b/opm/simulators/flow/FIPContainer.cpp index 95d6e37efa..78ee3aa6ec 100644 --- a/opm/simulators/flow/FIPContainer.cpp +++ b/opm/simulators/flow/FIPContainer.cpp @@ -337,6 +337,91 @@ assignCo2InWater(const unsigned globalDofIdx, } } +template +bool +FIPContainer:: +hasMicrobialMass() const +{ + return has(Inplace::Phase::MicrobialMass); +} + +template +void +FIPContainer:: +assignMicrobialMass(const unsigned globalDofIdx, + const Scalar microbialMass) +{ + this->fip_[Inplace::Phase::MicrobialMass][globalDofIdx] = microbialMass; +} + +template +bool +FIPContainer:: +hasOxygenMass() const +{ + return has(Inplace::Phase::OxygenMass); +} + +template +void +FIPContainer:: +assignOxygenMass(const unsigned globalDofIdx, + const Scalar oxygenMass) +{ + this->fip_[Inplace::Phase::OxygenMass][globalDofIdx] = oxygenMass; +} + +template +bool +FIPContainer:: +hasUreaMass() const +{ + return has(Inplace::Phase::UreaMass); +} + +template +void +FIPContainer:: +assignUreaMass(const unsigned globalDofIdx, + const Scalar ureaMass) +{ + this->fip_[Inplace::Phase::UreaMass][globalDofIdx] = ureaMass; +} + +template +bool +FIPContainer:: +hasBiofilmMass() const +{ + return has(Inplace::Phase::BiofilmMass); +} + +template +void +FIPContainer:: +assignBiofilmMass(const unsigned globalDofIdx, + const Scalar biofilmMass) +{ + this->fip_[Inplace::Phase::BiofilmMass][globalDofIdx] = biofilmMass; +} + +template +bool +FIPContainer:: +hasCalciteMass() const +{ + return has(Inplace::Phase::CalciteMass); +} + +template +void +FIPContainer:: +assignCalciteMass(const unsigned globalDofIdx, + const Scalar calciteMass) +{ + this->fip_[Inplace::Phase::CalciteMass][globalDofIdx] = calciteMass; +} + template void FIPContainer:: diff --git a/opm/simulators/flow/FIPContainer.hpp b/opm/simulators/flow/FIPContainer.hpp index a104e07129..53af3bba28 100644 --- a/opm/simulators/flow/FIPContainer.hpp +++ b/opm/simulators/flow/FIPContainer.hpp @@ -86,6 +86,26 @@ public: const Scalar co2InWater, const Scalar mM); + bool hasMicrobialMass() const; + void assignMicrobialMass(const unsigned globalDofIdx, + const Scalar microbialMass); + + bool hasOxygenMass() const; + void assignOxygenMass(const unsigned globalDofIdx, + const Scalar oxygenMass); + + bool hasUreaMass() const; + void assignUreaMass(const unsigned globalDofIdx, + const Scalar ureaMass); + + bool hasBiofilmMass() const; + void assignBiofilmMass(const unsigned globalDofIdx, + const Scalar biofilmMass); + + bool hasCalciteMass() const; + void assignCalciteMass(const unsigned globalDofIdx, + const Scalar calciteMass); + void assignGasWater(const unsigned globalDofIdx, const std::array& fip, const Scalar gasInPlaceWater, diff --git a/opm/simulators/flow/FlowGenericProblem.hpp b/opm/simulators/flow/FlowGenericProblem.hpp index cddb75992c..0fd2c82abc 100644 --- a/opm/simulators/flow/FlowGenericProblem.hpp +++ b/opm/simulators/flow/FlowGenericProblem.hpp @@ -173,7 +173,7 @@ public: /*! * \brief Returns the initial polymer concentration for a given a cell index */ - Scalar polymerConcentration(unsigned elemIdx) const; + Scalar polymerConcentration(unsigned elemIdx) const; /*! * \brief Returns the polymer molecule weight for a given cell index @@ -184,27 +184,27 @@ public: /*! * \brief Returns the initial microbial concentration for a given a cell index */ - Scalar microbialConcentration(unsigned elemIdx) const; + Scalar microbialConcentration(unsigned elemIdx) const; /*! * \brief Returns the initial oxygen concentration for a given a cell index */ - Scalar oxygenConcentration(unsigned elemIdx) const; + Scalar oxygenConcentration(unsigned elemIdx) const; /*! * \brief Returns the initial urea concentration for a given a cell index */ - Scalar ureaConcentration(unsigned elemIdx) const; + Scalar ureaConcentration(unsigned elemIdx) const; /*! * \brief Returns the initial biofilm concentration for a given a cell index */ - Scalar biofilmConcentration(unsigned elemIdx) const; + Scalar biofilmConcentration(unsigned elemIdx) const; /*! * \brief Returns the initial calcite concentration for a given a cell index */ - Scalar calciteConcentration(unsigned elemIdx) const; + Scalar calciteConcentration(unsigned elemIdx) const; /*! * \brief Returns the index the relevant PVT region given a cell index diff --git a/opm/simulators/flow/FlowProblem.hpp b/opm/simulators/flow/FlowProblem.hpp index bca9d13d39..25102a1451 100644 --- a/opm/simulators/flow/FlowProblem.hpp +++ b/opm/simulators/flow/FlowProblem.hpp @@ -1102,7 +1102,7 @@ public: /*! - * \brief Return the well transmissibility multiplier due to rock changues. + * \brief Return the well transmissibility multiplier due to rock changes. */ template LhsEval wellTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const @@ -1153,6 +1153,15 @@ public: case BCComponent::POLYMER: this->handlePolymerBC(bc, rate); break; + case BCComponent::MICR: + this->handleMicrBC(bc, rate); + break; + case BCComponent::OXYG: + this->handleOxygBC(bc, rate); + break; + case BCComponent::UREA: + this->handleUreaBC(bc, rate); + break; case BCComponent::NONE: throw std::logic_error("you need to specify the component when RATE type is set in BC"); break; @@ -1701,6 +1710,12 @@ protected: virtual void handlePolymerBC(const BCProp::BCFace&, RateVector&) const = 0; + virtual void handleMicrBC(const BCProp::BCFace&, RateVector&) const = 0; + + virtual void handleOxygBC(const BCProp::BCFace&, RateVector&) const = 0; + + virtual void handleUreaBC(const BCProp::BCFace&, RateVector&) const = 0; + BCData bcindex_; bool nonTrivialBoundaryConditions_ = false; bool first_step_ = true; diff --git a/opm/simulators/flow/FlowProblemBlackoil.hpp b/opm/simulators/flow/FlowProblemBlackoil.hpp index 0cb2d8e01f..e536c47e65 100644 --- a/opm/simulators/flow/FlowProblemBlackoil.hpp +++ b/opm/simulators/flow/FlowProblemBlackoil.hpp @@ -486,17 +486,6 @@ public: }); }); } - - // Deal with "clogging" for the MICP model - if constexpr (enableMICP) { - auto& model = this->model(); - const auto& residual = model.linearizer().residual(); - - for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); ++globalDofIdx) { - auto& phi = this->referencePorosity_[/*timeIdx=*/1][globalDofIdx]; - MICPModule::checkCloggingMICP(model, phi, globalDofIdx); - } - } } /*! @@ -629,6 +618,12 @@ public: if constexpr (enablePolymer) { rate[Indices::polymerConcentrationIdx] += source.rate({ijk, SourceComponent::POLYMER}) / this->model().dofTotalVolume(globalDofIdx); } + if constexpr (enableMICP) { + rate[Indices::microbialConcentrationIdx] += source.rate({ijk, SourceComponent::MICR}) / this->model().dofTotalVolume(globalDofIdx); + rate[Indices::oxygenConcentrationIdx] += source.rate({ijk, SourceComponent::OXYG}) / this->model().dofTotalVolume(globalDofIdx); + // dividing by scaling factor 10 (see WellInterfaceGeneric.cpp) + rate[Indices::ureaConcentrationIdx] += source.rate({ijk, SourceComponent::UREA}) / (this->model().dofTotalVolume(globalDofIdx) * 10); + } if constexpr (enableEnergy) { for (unsigned i = 0; i < phidx_map.size(); ++i) { const auto phaseIdx = phidx_map[i]; @@ -692,7 +687,7 @@ public: LhsEval permFactTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const { OPM_TIMEBLOCK_LOCAL(permFactTransMultiplier); - if (!enableSaltPrecipitation) + if (!enableSaltPrecipitation && !enableMICP) return 1.0; const auto& fs = intQuants.fluidState(); @@ -700,8 +695,15 @@ public: LhsEval porosityFactor = decay(1. - fs.saltSaturation()); porosityFactor = min(porosityFactor, 1.0); - const auto& permfactTable = BrineModule::permfactTable(tableIdx); - return permfactTable.eval(porosityFactor, /*extrapolation=*/true); + if constexpr (enableSaltPrecipitation) { + const auto& permfactTable = BrineModule::permfactTable(tableIdx); + return permfactTable.eval(porosityFactor, /*extrapolation=*/true); + } + if constexpr (enableMICP){ + const auto& permfactTable = MICPModule::permfactTable(tableIdx); + return permfactTable.eval(porosityFactor, /*extrapolation=*/true); + } + return 1.0; } // temporary solution to facilitate output of initial state from flow @@ -763,6 +765,9 @@ public: break; case BCComponent::SOLVENT: case BCComponent::POLYMER: + case BCComponent::MICR: + case BCComponent::OXYG: + case BCComponent::UREA: case BCComponent::NONE: throw std::logic_error("you need to specify a valid component (OIL, WATER or GAS) when DIRICHLET type is set in BC"); } @@ -1516,6 +1521,31 @@ protected: rate[Indices::polymerConcentrationIdx] = bc.rate; } + void handleMicrBC(const BCProp::BCFace& bc, RateVector& rate) const override + { + if constexpr (!enableMICP) + throw std::logic_error("MICP is disabled and you're trying to add microbes to BC"); + + rate[Indices::microbialConcentrationIdx] = bc.rate; + } + + void handleOxygBC(const BCProp::BCFace& bc, RateVector& rate) const override + { + if constexpr (!enableMICP) + throw std::logic_error("MICP is disabled and you're trying to add oxygen to BC"); + + rate[Indices::oxygenConcentrationIdx] = bc.rate; + } + + void handleUreaBC(const BCProp::BCFace& bc, RateVector& rate) const override + { + if constexpr (!enableMICP) + throw std::logic_error("MICP is disabled and you're trying to add urea to BC"); + + // dividing by scaling factor 10 (see WellInterfaceGeneric.cpp) + rate[Indices::ureaConcentrationIdx] = bc.rate / 10; + } + void updateExplicitQuantities_(const bool first_step_after_restart) { OPM_TIMEBLOCK(updateExplicitQuantities); diff --git a/opm/simulators/flow/FlowProblemComp.hpp b/opm/simulators/flow/FlowProblemComp.hpp index 212d2869ea..fc561849a9 100644 --- a/opm/simulators/flow/FlowProblemComp.hpp +++ b/opm/simulators/flow/FlowProblemComp.hpp @@ -607,6 +607,21 @@ private: throw std::logic_error("polymer is disabled for compositional modeling and you're trying to add polymer to BC"); } + void handleMicrBC(const BCProp::BCFace& /* bc */, RateVector& /* rate */) const override + { + throw std::logic_error("MICP is disabled for compositional modeling and you're trying to add microbes to BC"); + } + + void handleOxygBC(const BCProp::BCFace& /* bc */, RateVector& /* rate */) const override + { + throw std::logic_error("MICP is disabled for compositional modeling and you're trying to add oxygen to BC"); + } + + void handleUreaBC(const BCProp::BCFace& /* bc */, RateVector& /* rate */) const override + { + throw std::logic_error("MICP is disabled for compositional modeling and you're trying to add urea to BC"); + } + FlowThresholdPressure thresholdPressures_; std::vector initialFluidStates_; diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index 8b581fb113..59577e1574 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -708,7 +708,7 @@ doAllocBuffers(const unsigned bufferSize, Entry{&cFoam_, "", enableFoam_}, Entry{&cSalt_, "", enableBrine_}, Entry{&pSalt_, "", enableSaltPrecipitation_}, - Entry{&permFact_, "", enableSaltPrecipitation_}, + Entry{&permFact_, "", enableSaltPrecipitation_ || enableMICP_}, Entry{&soMax_, "", oilvap.getType() == OilVapP::VAPPARS}, Entry{&soMax_, "", hysteresisConfig && hysteresisConfig->enableNonWettingHysteresis() && diff --git a/opm/simulators/flow/OutputBlackoilModule.hpp b/opm/simulators/flow/OutputBlackoilModule.hpp index 74a27aa8e1..57bd5ec4c3 100644 --- a/opm/simulators/flow/OutputBlackoilModule.hpp +++ b/opm/simulators/flow/OutputBlackoilModule.hpp @@ -114,6 +114,7 @@ class OutputBlackOilModule : public GenericOutputBlackoilModule() }; + enum { enableMICP = getPropValue() }; template static Scalar value_or_zero(const VectorType& v) @@ -873,6 +874,25 @@ private: { this->updateCO2InWater(globalDofIdx, pv, fs); } + + if constexpr(enableMICP) { + const auto surfVolWat = pv * getValue(fs.invB(waterPhaseIdx)); + if (this->fipC_.hasMicrobialMass()) { + this->updateMicrobialMass(globalDofIdx, intQuants, surfVolWat); + } + if (this->fipC_.hasOxygenMass()) { + this->updateOxygenMass(globalDofIdx, intQuants, surfVolWat); + } + if (this->fipC_.hasUreaMass()) { + this->updateUreaMass(globalDofIdx, intQuants, surfVolWat); + } + if (this->fipC_.hasBiofilmMass()) { + this->updateBiofilmMass(globalDofIdx, intQuants, totVolume); + } + if (this->fipC_.hasCalciteMass()) { + this->updateCalciteMass(globalDofIdx, intQuants, totVolume); + } + } } template @@ -991,6 +1011,57 @@ private: return xoG * pv * rhoo * so / mM; } + template + void updateMicrobialMass(const unsigned globalDofIdx, + const IntensiveQuantities& intQuants, + const double surfVolWat) + { + const Scalar mass = surfVolWat * intQuants.microbialConcentration().value(); + + this->fipC_.assignMicrobialMass(globalDofIdx, mass); + } + + template + void updateOxygenMass(const unsigned globalDofIdx, + const IntensiveQuantities& intQuants, + const double surfVolWat) + { + const Scalar mass = surfVolWat * intQuants.oxygenConcentration().value(); + + this->fipC_.assignOxygenMass(globalDofIdx, mass); + } + + template + void updateUreaMass(const unsigned globalDofIdx, + const IntensiveQuantities& intQuants, + const double surfVolWat) + { + // rescaling back the urea concentration (see WellInterfaceGeneric.cpp) + const Scalar mass = surfVolWat * intQuants.ureaConcentration().value() * 10; + + this->fipC_.assignUreaMass(globalDofIdx, mass); + } + + template + void updateBiofilmMass(const unsigned globalDofIdx, + const IntensiveQuantities& intQuants, + const double totVolume) + { + const Scalar mass = totVolume * intQuants.biofilmMass().value(); + + this->fipC_.assignBiofilmMass(globalDofIdx, mass); + } + + template + void updateCalciteMass(const unsigned globalDofIdx, + const IntensiveQuantities& intQuants, + const double totVolume) + { + const Scalar mass = totVolume * intQuants.calciteMass().value(); + + this->fipC_.assignCalciteMass(globalDofIdx, mass); + } + //! \brief Setup extractors for element-level data. void setupElementExtractors_() { @@ -1441,7 +1512,7 @@ private: micpC.assign(ectx.globalDofIdx, ectx.intQuants.microbialConcentration().value(), ectx.intQuants.oxygenConcentration().value(), - // Rescaling back the urea concentration (see WellInterface_impl.hpp) + // rescaling back the urea concentration (see WellInterfaceGeneric.cpp) 10 * ectx.intQuants.ureaConcentration().value(), ectx.intQuants.biofilmConcentration().value(), ectx.intQuants.calciteConcentration().value()); @@ -1958,6 +2029,50 @@ private: } } }, + Entry{ScalarEntry{"BMMIP", + [&model = this->simulator_.model()](const Context& ectx) + { + return getValue(ectx.intQuants.microbialConcentration()) * + getValue(ectx.intQuants.porosity()) * + model.dofTotalVolume(ectx.globalDofIdx); + } + } + }, + Entry{ScalarEntry{"BMOIP", + [&model = this->simulator_.model()](const Context& ectx) + { + return getValue(ectx.intQuants.oxygenConcentration()) * + getValue(ectx.intQuants.porosity()) * + model.dofTotalVolume(ectx.globalDofIdx); + } + } + }, + Entry{ScalarEntry{"BMUIP", + [&model = this->simulator_.model()](const Context& ectx) + { + // rescaling back the urea concentration (see WellInterfaceGeneric.cpp) + return getValue(ectx.intQuants.ureaConcentration()) * + getValue(ectx.intQuants.porosity()) * + model.dofTotalVolume(ectx.globalDofIdx) * 10; + } + } + }, + Entry{ScalarEntry{"BMBIP", + [&model = this->simulator_.model()](const Context& ectx) + { + return model.dofTotalVolume(ectx.globalDofIdx) * + getValue(ectx.intQuants.biofilmMass()); + } + } + }, + Entry{ScalarEntry{"BMCIP", + [&model = this->simulator_.model()](const Context& ectx) + { + return model.dofTotalVolume(ectx.globalDofIdx) * + getValue(ectx.intQuants.calciteMass()); + } + } + }, }; this->blockExtractors_ = BlockExtractor::setupExecMap(this->blockData_, handlers); diff --git a/opm/simulators/utils/BlackoilPhases.cpp b/opm/simulators/utils/BlackoilPhases.cpp index be5971a1b5..70b58807fe 100644 --- a/opm/simulators/utils/BlackoilPhases.cpp +++ b/opm/simulators/utils/BlackoilPhases.cpp @@ -55,4 +55,5 @@ Opm::PhaseUsage::PhaseUsage(std::vector phases) this->has_foam = this->phase_used[BlackoilPhases::Foam]; this->has_brine = this->phase_used[BlackoilPhases::Brine]; this->has_zFraction = this->phase_used[BlackoilPhases::ZFraction]; + this->has_micp = this->phase_used[BlackoilPhases::MICP]; } diff --git a/opm/simulators/utils/BlackoilPhases.hpp b/opm/simulators/utils/BlackoilPhases.hpp index 39e0a62e82..8a99423a07 100644 --- a/opm/simulators/utils/BlackoilPhases.hpp +++ b/opm/simulators/utils/BlackoilPhases.hpp @@ -36,10 +36,10 @@ namespace Opm // to and from active ones. That said, they are not considered by num_phases or // MaxNumPhases. The crypto phases which are currently implemented are solvent, // polymer, energy, polymer molecular weight, foam and brine. - static const int NumCryptoPhases = 7; + static const int NumCryptoPhases = 8; // enum ComponentIndex { Water = 0, Oil = 1, Gas = 2 }; - enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2, Solvent = 3, Polymer = 4, Energy = 5, PolymerMW = 6, Foam = 7, Brine = 8, ZFraction = 9 }; + enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2, Solvent = 3, Polymer = 4, Energy = 5, PolymerMW = 6, Foam = 7, Brine = 8, ZFraction = 9, MICP = 10 }; }; struct PhaseUsage : public BlackoilPhases @@ -60,6 +60,7 @@ namespace Opm bool has_foam{}; bool has_brine{}; bool has_zFraction{}; + bool has_micp{}; }; diff --git a/opm/simulators/utils/phaseUsageFromDeck.cpp b/opm/simulators/utils/phaseUsageFromDeck.cpp index c9bd30e763..8cdfb21db2 100644 --- a/opm/simulators/utils/phaseUsageFromDeck.cpp +++ b/opm/simulators/utils/phaseUsageFromDeck.cpp @@ -134,10 +134,32 @@ PhaseUsage phaseUsage(const Phases& phases) return pu; } +PhaseUsage phaseUsageMICP() +{ + PhaseUsage pu; + pu.phase_used.fill(0); + for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases + BlackoilPhases::NumCryptoPhases; ++phaseIdx) { + pu.phase_used[phaseIdx] = 0; + pu.phase_pos[phaseIdx] = -1; + } + pu.phase_used[BlackoilPhases::Aqua] = 1; + pu.phase_used[BlackoilPhases::MICP] = 1; + pu.phase_pos[BlackoilPhases::Aqua] = 0; + pu.phase_pos[BlackoilPhases::MICP] = 1; + pu.num_phases = 1; + pu.has_micp = true; + + return pu; +} + PhaseUsage phaseUsageFromDeck(const EclipseState& eclipseState) { const auto& phases = eclipseState.runspec().phases(); + if (eclipseState.runspec().micp()) { + return phaseUsageMICP(); + } + return phaseUsage(phases); } diff --git a/opm/simulators/utils/phaseUsageFromDeck.hpp b/opm/simulators/utils/phaseUsageFromDeck.hpp index ae7403c600..473982a5f7 100644 --- a/opm/simulators/utils/phaseUsageFromDeck.hpp +++ b/opm/simulators/utils/phaseUsageFromDeck.hpp @@ -39,6 +39,9 @@ namespace Opm /// Looks at presence of WATER, OIL and GAS keywords in deck /// to determine active phases. PhaseUsage phaseUsageFromDeck(const Deck& deck); + + /// Handels MICP + PhaseUsage phaseUsageMICP(const Phases& phases); } #endif // OPM_PHASEUSAGEFROMDECK_HEADER_INCLUDED diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index b8f62dee1a..80fd2b8593 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -98,6 +98,10 @@ namespace Opm OPM_THROW(std::runtime_error, "water evaporation is not supported by multisegment well yet"); } + if constexpr (Base::has_micp) { + OPM_THROW(std::runtime_error, "MICP is not supported by multisegment well yet"); + } + if(this->rsRvInj() > 0) { OPM_THROW(std::runtime_error, "dissolved gas/ vapporized oil in injected oil/gas not supported by multisegment well yet." diff --git a/opm/simulators/wells/PerfData.cpp b/opm/simulators/wells/PerfData.cpp index 97096a93ff..34a05acffa 100644 --- a/opm/simulators/wells/PerfData.cpp +++ b/opm/simulators/wells/PerfData.cpp @@ -42,7 +42,9 @@ PerfData::PerfData(const std::size_t num_perf, , polymer_rates(num_perf) , brine_rates(num_perf) , prod_index(num_perf * num_phases) - , micp_rates(num_perf) + , microbial_rates(num_perf) + , oxygen_rates(num_perf) + , urea_rates(num_perf) , cell_index(num_perf) , connection_transmissibility_factor(num_perf) , connection_d_factor(num_perf) @@ -79,7 +81,9 @@ PerfData PerfData::serializationTestObject() result.polymer_rates = {10.0, 11.0, 12.0}; result.brine_rates = {13.0}; result.prod_index = {14.0, 15.0}; - result.micp_rates = {16.0}; + result.microbial_rates = {16.0}; + result.oxygen_rates = {16.0}; + result.urea_rates = {16.0}; result.cell_index = {17, 18, 19, 20}; result.connection_transmissibility_factor = {21.0}; result.connection_d_factor = {21.5}; @@ -127,7 +131,9 @@ bool PerfData::try_assign(const PerfData& other) this->polymer_rates = other.polymer_rates; this->brine_rates = other.brine_rates; this->prod_index = other.prod_index; - this->micp_rates = other.micp_rates; + this->microbial_rates = other.microbial_rates; + this->oxygen_rates = other.oxygen_rates; + this->urea_rates = other.urea_rates; this->water_throughput = other.water_throughput; this->skin_pressure = other.skin_pressure; this->water_velocity = other.water_velocity; @@ -150,7 +156,9 @@ bool PerfData::operator==(const PerfData& rhs) const && (this->polymer_rates == rhs.polymer_rates) && (this->brine_rates == rhs.brine_rates) && (this->prod_index == rhs.prod_index) - && (this->micp_rates == rhs.micp_rates) + && (this->microbial_rates == rhs.microbial_rates) + && (this->oxygen_rates == rhs.oxygen_rates) + && (this->urea_rates == rhs.urea_rates) && (this->cell_index == rhs.cell_index) && (this->connection_transmissibility_factor == rhs.connection_transmissibility_factor) && (this->connection_d_factor == rhs.connection_d_factor) diff --git a/opm/simulators/wells/PerfData.hpp b/opm/simulators/wells/PerfData.hpp index c3b99d3938..eefcd8da83 100644 --- a/opm/simulators/wells/PerfData.hpp +++ b/opm/simulators/wells/PerfData.hpp @@ -66,7 +66,9 @@ public: serializer(polymer_rates); serializer(brine_rates); serializer(prod_index); - serializer(micp_rates); + serializer(microbial_rates); + serializer(oxygen_rates); + serializer(urea_rates); serializer(cell_index); serializer(connection_transmissibility_factor); serializer(connection_d_factor); @@ -97,7 +99,9 @@ public: std::vector polymer_rates{}; std::vector brine_rates{}; std::vector prod_index{}; - std::vector micp_rates{}; + std::vector microbial_rates{}; + std::vector oxygen_rates{}; + std::vector urea_rates{}; std::vector cell_index{}; std::vector connection_transmissibility_factor{}; std::vector connection_d_factor{}; diff --git a/opm/simulators/wells/SingleWellState.cpp b/opm/simulators/wells/SingleWellState.cpp index 5696f489e6..6aafc518a3 100644 --- a/opm/simulators/wells/SingleWellState.cpp +++ b/opm/simulators/wells/SingleWellState.cpp @@ -193,6 +193,24 @@ Scalar SingleWellState::sum_solvent_rates() const return this->sum_connection_rates(this->perf_data.solvent_rates); } +template +Scalar SingleWellState::sum_microbial_rates() const +{ + return this->sum_connection_rates(this->perf_data.microbial_rates); +} + +template +Scalar SingleWellState::sum_oxygen_rates() const +{ + return this->sum_connection_rates(this->perf_data.oxygen_rates); +} + +template +Scalar SingleWellState::sum_urea_rates() const +{ + return this->sum_connection_rates(this->perf_data.urea_rates); +} + template Scalar SingleWellState::sum_filtrate_rate() const { diff --git a/opm/simulators/wells/SingleWellState.hpp b/opm/simulators/wells/SingleWellState.hpp index f96e3522a6..93446174be 100644 --- a/opm/simulators/wells/SingleWellState.hpp +++ b/opm/simulators/wells/SingleWellState.hpp @@ -144,6 +144,9 @@ public: Scalar sum_solvent_rates() const; Scalar sum_polymer_rates() const; Scalar sum_brine_rates() const; + Scalar sum_microbial_rates() const; + Scalar sum_oxygen_rates() const; + Scalar sum_urea_rates() const; Scalar sum_filtrate_rate() const; Scalar sum_filtrate_total() const; diff --git a/opm/simulators/wells/StandardWellConnections.cpp b/opm/simulators/wells/StandardWellConnections.cpp index d333ee0fce..1de062ba09 100644 --- a/opm/simulators/wells/StandardWellConnections.cpp +++ b/opm/simulators/wells/StandardWellConnections.cpp @@ -756,7 +756,10 @@ std::tuple::Eval, typename StandardWellConnections::Eval, typename StandardWellConnections::Eval> StandardWellConnections:: -connectionRatesMICP(const std::vector& cq_s, +connectionRatesMICP(Scalar& rate_m, + Scalar& rate_o, + Scalar& rate_u, + const std::vector& cq_s, const std::variant& microbialConcentration, const std::variant& oxygenConcentration, const std::variant& ureaConcentration) const @@ -769,6 +772,8 @@ connectionRatesMICP(const std::vector& cq_s, cq_s_microbe *= std::get(microbialConcentration); } + rate_m = cq_s_microbe.value(); + EvalWell cq_s_oxygen = cq_s[waterCompIdx]; if (well_.isInjector()) { cq_s_oxygen *= std::get(oxygenConcentration); @@ -776,6 +781,8 @@ connectionRatesMICP(const std::vector& cq_s, cq_s_oxygen *= std::get(oxygenConcentration); } + rate_o = cq_s_oxygen.value(); + EvalWell cq_s_urea = cq_s[waterCompIdx]; if (well_.isInjector()) { cq_s_urea *= std::get(ureaConcentration); @@ -783,6 +790,9 @@ connectionRatesMICP(const std::vector& cq_s, cq_s_urea *= std::get(ureaConcentration); } + // rescaling back the urea concentration (see WellInterfaceGeneric.cpp) + rate_u = cq_s_urea.value() * 10; + return {well_.restrictEval(cq_s_microbe), well_.restrictEval(cq_s_oxygen), well_.restrictEval(cq_s_urea)}; diff --git a/opm/simulators/wells/StandardWellConnections.hpp b/opm/simulators/wells/StandardWellConnections.hpp index 56016732f7..86cefe4435 100644 --- a/opm/simulators/wells/StandardWellConnections.hpp +++ b/opm/simulators/wells/StandardWellConnections.hpp @@ -124,7 +124,10 @@ public: const std::variant& polymerConcentration) const; std::tuple - connectionRatesMICP(const std::vector& cq_s, + connectionRatesMICP(Scalar& rate_m, + Scalar& rate_o, + Scalar& rate_u, + const std::vector& cq_s, const std::variant& microbialConcentration, const std::variant& oxygenConcentration, const std::variant& ureaConcentration) const; diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 70199a5e15..eecf2aa270 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -623,7 +623,10 @@ namespace Opm std::tie(connectionRates[perf][Indices::contiMicrobialEqIdx], connectionRates[perf][Indices::contiOxygenEqIdx], connectionRates[perf][Indices::contiUreaEqIdx]) = - this->connections_.connectionRatesMICP(cq_s, + this->connections_.connectionRatesMICP(perf_data.microbial_rates[perf], + perf_data.oxygen_rates[perf], + perf_data.urea_rates[perf], + cq_s, microbialConcentration, oxygenConcentration, ureaConcentration); diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index c4f6f6c0d9..302d73de08 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -785,6 +785,11 @@ Scalar WellInterfaceGeneric::wpolymer_() const } } +// The urea injection concentration is scaled down by a factor of 10, since its value +// can be much bigger than 1 (not doing this slows the simulations). The corresponding +// values are scaled accordingly when computing the reactions and also when writing +// the output files. + template Scalar WellInterfaceGeneric::wurea_() const { @@ -792,7 +797,7 @@ Scalar WellInterfaceGeneric::wurea_() const if (injectorType == InjectorType::WATER) { WellMICPProperties urea = this->well_ecl_.getMICPProperties(); - const Scalar urea_injection_concentration = urea.m_ureaConcentration / 10.; //Dividing by scaling factor 10 + const Scalar urea_injection_concentration = urea.m_ureaConcentration / 10;//Dividing by scaling factor 10 return urea_injection_concentration; } else { // Not a water injection well => no urea. diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 65a43bffc2..7465d339d5 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -173,12 +173,6 @@ namespace Opm return 0.0; } - // The urea injection concentration is scaled down by a factor of 10, since its value - // can be much bigger than 1 (not doing this slows the simulations). The - // corresponding values are scaled accordingly in blackoilmicpmodules.hh when computing - // the reactions and also when writing the output files (vtk and eclipse format, i.e., - // vtkblackoilmicpmodule.hh and ecloutputblackoilmodel.hh respectively). - template typename WellInterface::Scalar WellInterface:: diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp index f3d94f7fd1..efa93802ca 100644 --- a/opm/simulators/wells/WellState.cpp +++ b/opm/simulators/wells/WellState.cpp @@ -577,6 +577,12 @@ WellState::report(const int* globalCellIdxMap, well.rates.set(rt::brine, ws.sum_brine_rates()); } + if (pu.has_micp) { + well.rates.set(rt::microbial, ws.sum_microbial_rates()); + well.rates.set(rt::oxygen, ws.sum_oxygen_rates()); + well.rates.set(rt::urea, ws.sum_urea_rates()); + } + if (ws.producer) { well.rates.set(rt::alq, getALQ(wname)); } @@ -1123,6 +1129,14 @@ reportConnectionPressuresAndRates(const std::size_t well_index, connections[i].rates.set(rt::solvent, perf_data.solvent_rates[i]); } } + + if (pu.has_micp) { + for (auto i = 0*num_perf_well; i < num_perf_well; ++i) { + connections[i].rates.set(rt::microbial, perf_data.microbial_rates[i]); + connections[i].rates.set(rt::oxygen, perf_data.oxygen_rates[i]); + connections[i].rates.set(rt::urea, perf_data.urea_rates[i]); + } + } } template From 1541e84f7df6ffaa6002d2fdf89f61e18899367d Mon Sep 17 00:00:00 2001 From: David Landa Marban Date: Wed, 26 Mar 2025 14:28:36 +0100 Subject: [PATCH 2/3] Adressing akva2 comments --- flow/flow_micp.cpp | 9 +++--- ...BlackOilIntensiveQuantitiesGlobalIndex.hpp | 1 + opm/models/blackoil/blackoilindices.hh | 32 +++++++++---------- .../blackoil/blackoilintensivequantities.hh | 10 ++++-- opm/models/blackoil/blackoilmicpparams.cpp | 2 +- .../blackoil/blackoilonephaseindices.hh | 32 +++++++++---------- .../blackoil/blackoiltwophaseindices.hh | 32 +++++++++---------- opm/models/io/vtkblackoilmicpmodule.hpp | 2 +- .../flow/FlowGenericProblem_impl.hpp | 7 +++- opm/simulators/utils/phaseUsageFromDeck.hpp | 2 +- opm/simulators/wells/WellInterfaceGeneric.cpp | 2 +- 11 files changed, 71 insertions(+), 60 deletions(-) diff --git a/flow/flow_micp.cpp b/flow/flow_micp.cpp index aeef32e657..48edcd0bb6 100644 --- a/flow/flow_micp.cpp +++ b/flow/flow_micp.cpp @@ -18,15 +18,16 @@ #include +#include + #include #include - -#include -#include -#include #include #include +#include +#include + namespace Opm { namespace Properties { namespace TTag { diff --git a/flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp b/flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp index d20e805a7b..a080953d0f 100644 --- a/flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp +++ b/flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp @@ -53,6 +53,7 @@ #include +#include #include #include diff --git a/opm/models/blackoil/blackoilindices.hh b/opm/models/blackoil/blackoilindices.hh index d0760ea60f..53f4ed6bba 100644 --- a/opm/models/blackoil/blackoilindices.hh +++ b/opm/models/blackoil/blackoilindices.hh @@ -147,24 +147,24 @@ struct BlackOilIndices numPolymers > 1 ? polymerConcentrationIdx + 1 : -1000; //! Index of the primary variable for the first MICP component - static constexpr int biofilmConcentrationIdx = + static constexpr int microbialConcentrationIdx = enableMICP ? PVOffset + numPhases + numSolvents : -1000; //! Index of the primary variable for the second MICP component - static constexpr int microbialConcentrationIdx = - numMICPs > 1 ? biofilmConcentrationIdx + 1 : -1000; + static constexpr int oxygenConcentrationIdx = + numMICPs > 1 ? microbialConcentrationIdx + 1 : -1000; //! Index of the primary variable for the third MICP component - static constexpr int oxygenConcentrationIdx = - numMICPs > 2 ? microbialConcentrationIdx + 1 : -1000; + static constexpr int ureaConcentrationIdx = + numMICPs > 2 ? oxygenConcentrationIdx + 1 : -1000; //! Index of the primary variable for the fourth MICP component - static constexpr int ureaConcentrationIdx = - numMICPs > 3 ? oxygenConcentrationIdx + 1 : -1000; + static constexpr int biofilmConcentrationIdx = + numMICPs > 3 ? ureaConcentrationIdx + 1 : -1000; //! Index of the primary variable for the fifth MICP component static constexpr int calciteConcentrationIdx = - numMICPs > 4 ? ureaConcentrationIdx + 1 : -1000; + numMICPs > 4 ? biofilmConcentrationIdx + 1 : -1000; //! Index of the primary variable for the foam static constexpr int foamConcentrationIdx = @@ -204,24 +204,24 @@ struct BlackOilIndices numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000; //! Index of the continuity equation for the first MICP component - static constexpr int contiBiofilmEqIdx = + static constexpr int contiMicrobialEqIdx = enableMICP ? PVOffset + numPhases + numSolvents + numExtbos : -1000; //! Index of the continuity equation for the second MICP component - static constexpr int contiMicrobialEqIdx = - numMICPs > 1 ? contiBiofilmEqIdx + 1 : -1000; + static constexpr int contiOxygenEqIdx = + numMICPs > 1 ? contiMicrobialEqIdx + 1 : -1000; //! Index of the continuity equation for the third MICP component - static constexpr int contiOxygenEqIdx = - numMICPs > 2 ? contiMicrobialEqIdx + 1 : -1000; + static constexpr int contiUreaEqIdx = + numMICPs > 2 ? contiOxygenEqIdx + 1 : -1000; //! Index of the continuity equation for the fourth MICP component - static constexpr int contiUreaEqIdx = - numMICPs > 3 ? contiOxygenEqIdx + 1 : -1000; + static constexpr int contiBiofilmEqIdx = + numMICPs > 3 ? contiUreaEqIdx + 1 : -1000; //! Index of the continuity equation for the fifth MICP component static constexpr int contiCalciteEqIdx = - numMICPs > 4 ? contiUreaEqIdx + 1 : -1000; + numMICPs > 4 ? contiBiofilmEqIdx + 1 : -1000; //! Index of the continuity equation for the foam component static constexpr int contiFoamEqIdx = diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 074a460736..de87a64747 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -696,11 +696,15 @@ public: const Evaluation& permFactor() const { - if constexpr (enableMICP) + if constexpr (enableMICP) { return MICPIntQua::permFactor(); - if constexpr (enableSaltPrecipitation) + } + else if constexpr (enableSaltPrecipitation) { return BrineIntQua::permFactor(); - throw std::logic_error("permFactor() called but salt precipitation or MICP are disabled"); + } + else { + throw std::logic_error("permFactor() called but salt precipitation or MICP are disabled"); + } } private: diff --git a/opm/models/blackoil/blackoilmicpparams.cpp b/opm/models/blackoil/blackoilmicpparams.cpp index 1ac9e937e2..fd6c2da7de 100644 --- a/opm/models/blackoil/blackoilmicpparams.cpp +++ b/opm/models/blackoil/blackoilmicpparams.cpp @@ -29,8 +29,8 @@ #include #include #include -#include #include +#include #endif #include diff --git a/opm/models/blackoil/blackoilonephaseindices.hh b/opm/models/blackoil/blackoilonephaseindices.hh index 9305fdc4b7..863c8b00ea 100644 --- a/opm/models/blackoil/blackoilonephaseindices.hh +++ b/opm/models/blackoil/blackoilonephaseindices.hh @@ -145,24 +145,24 @@ struct BlackOilOnePhaseIndices numPolymers > 1 ? polymerConcentrationIdx + 1 : -1000; //! Index of the primary variable for the first MICP component - static constexpr int biofilmConcentrationIdx = + static constexpr int microbialConcentrationIdx = enableMICP ? PVOffset + numPhases + numSolvents : -1000; //! Index of the primary variable for the second MICP component - static constexpr int microbialConcentrationIdx = - numMICPs > 1 ? biofilmConcentrationIdx + 1 : -1000; + static constexpr int oxygenConcentrationIdx = + numMICPs > 1 ? microbialConcentrationIdx + 1 : -1000; //! Index of the primary variable for the third MICP component - static constexpr int oxygenConcentrationIdx = - numMICPs > 2 ? microbialConcentrationIdx + 1 : -1000; + static constexpr int ureaConcentrationIdx = + numMICPs > 2 ? oxygenConcentrationIdx + 1 : -1000; //! Index of the primary variable for the fourth MICP component - static constexpr int ureaConcentrationIdx = - numMICPs > 3 ? oxygenConcentrationIdx + 1 : -1000; + static constexpr int biofilmConcentrationIdx = + numMICPs > 3 ? ureaConcentrationIdx + 1 : -1000; //! Index of the primary variable for the fifth MICP component static constexpr int calciteConcentrationIdx = - numMICPs > 4 ? ureaConcentrationIdx + 1 : -1000; + numMICPs > 4 ? biofilmConcentrationIdx + 1 : -1000; //! Index of the primary variable for the foam static constexpr int foamConcentrationIdx = @@ -221,24 +221,24 @@ struct BlackOilOnePhaseIndices numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000; //! Index of the continuity equation for the first MICP component - static constexpr int contiBiofilmEqIdx = + static constexpr int contiMicrobialEqIdx = enableMICP ? PVOffset + numPhases + numSolvents : -1000; //! Index of the continuity equation for the second MICP component - static constexpr int contiMicrobialEqIdx = - numMICPs > 1 ? contiBiofilmEqIdx + 1 : -1000; + static constexpr int contiOxygenEqIdx = + numMICPs > 1 ? contiMicrobialEqIdx + 1 : -1000; //! Index of the continuity equation for the third MICP component - static constexpr int contiOxygenEqIdx = - numMICPs > 2 ? contiMicrobialEqIdx + 1 : -1000; + static constexpr int contiUreaEqIdx = + numMICPs > 2 ? contiOxygenEqIdx + 1 : -1000; //! Index of the continuity equation for the fourth MICP component - static constexpr int contiUreaEqIdx = - numMICPs > 3 ? contiOxygenEqIdx + 1 : -1000; + static constexpr int contiBiofilmEqIdx = + numMICPs > 3 ? contiUreaEqIdx + 1 : -1000; //! Index of the continuity equation for the fifth MICP component static constexpr int contiCalciteEqIdx = - numMICPs > 4 ? contiUreaEqIdx + 1 : -1000; + numMICPs > 4 ? contiBiofilmEqIdx + 1 : -1000; //! Index of the continuity equation for the foam component static constexpr int contiFoamEqIdx = diff --git a/opm/models/blackoil/blackoiltwophaseindices.hh b/opm/models/blackoil/blackoiltwophaseindices.hh index 27d1630bb6..9d112a1488 100644 --- a/opm/models/blackoil/blackoiltwophaseindices.hh +++ b/opm/models/blackoil/blackoiltwophaseindices.hh @@ -145,24 +145,24 @@ struct BlackOilTwoPhaseIndices numPolymers > 1 ? polymerConcentrationIdx + 1 : -1000; //! Index of the primary variable for the first MICP component - static constexpr int biofilmConcentrationIdx = + static constexpr int microbialConcentrationIdx = enableMICP ? PVOffset + numPhases + numSolvents : -1000; //! Index of the primary variable for the second MICP component - static constexpr int microbialConcentrationIdx = - numMICPs > 1 ? biofilmConcentrationIdx + 1 : -1000; + static constexpr int oxygenConcentrationIdx = + numMICPs > 1 ? microbialConcentrationIdx + 1 : -1000; //! Index of the primary variable for the third MICP component - static constexpr int oxygenConcentrationIdx = - numMICPs > 2 ? microbialConcentrationIdx + 1 : -1000; + static constexpr int ureaConcentrationIdx = + numMICPs > 2 ? oxygenConcentrationIdx + 1 : -1000; //! Index of the primary variable for the fourth MICP component - static constexpr int ureaConcentrationIdx = - numMICPs > 3 ? oxygenConcentrationIdx + 1 : -1000; + static constexpr int biofilmConcentrationIdx = + numMICPs > 3 ? ureaConcentrationIdx + 1 : -1000; //! Index of the primary variable for the fifth MICP component static constexpr int calciteConcentrationIdx = - numMICPs > 4 ? ureaConcentrationIdx + 1 : -1000; + numMICPs > 4 ? biofilmConcentrationIdx + 1 : -1000; //! Index of the primary variable for the foam static constexpr int foamConcentrationIdx = @@ -240,24 +240,24 @@ struct BlackOilTwoPhaseIndices numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000; //! Index of the continuity equation for the first MICP component - static constexpr int contiBiofilmEqIdx = + static constexpr int contiMicrobialEqIdx = enableMICP ? PVOffset + numPhases + numSolvents : -1000; //! Index of the continuity equation for the second MICP component - static constexpr int contiMicrobialEqIdx = - numMICPs > 1 ? contiBiofilmEqIdx + 1 : -1000; + static constexpr int contiOxygenEqIdx = + numMICPs > 1 ? contiMicrobialEqIdx + 1 : -1000; //! Index of the continuity equation for the third MICP component - static constexpr int contiOxygenEqIdx = - numMICPs > 2 ? contiMicrobialEqIdx + 1 : -1000; + static constexpr int contiUreaEqIdx = + numMICPs > 2 ? contiOxygenEqIdx + 1 : -1000; //! Index of the continuity equation for the fourth MICP component - static constexpr int contiUreaEqIdx = - numMICPs > 3 ? contiOxygenEqIdx + 1 : -1000; + static constexpr int contiBiofilmEqIdx = + numMICPs > 3 ? contiUreaEqIdx + 1 : -1000; //! Index of the continuity equation for the fifth MICP component static constexpr int contiCalciteEqIdx = - numMICPs > 4 ? contiUreaEqIdx + 1 : -1000; + numMICPs > 4 ? contiBiofilmEqIdx + 1 : -1000; //! Index of the continuity equation for the foam component static constexpr int contiFoamEqIdx = diff --git a/opm/models/io/vtkblackoilmicpmodule.hpp b/opm/models/io/vtkblackoilmicpmodule.hpp index 4e791b8cb7..713f294b1c 100644 --- a/opm/models/io/vtkblackoilmicpmodule.hpp +++ b/opm/models/io/vtkblackoilmicpmodule.hpp @@ -141,7 +141,7 @@ public: } if (params_.ureaConcentrationOutput_) { - // multypliging by scaling factor 10 (see WellInterfaceGeneric.cpp) + // multiplying by scaling factor 10 (see WellInterfaceGeneric.cpp) ureaConcentration_[globalDofIdx] = 10 * scalarValue(intQuants.ureaConcentration()); } diff --git a/opm/simulators/flow/FlowGenericProblem_impl.hpp b/opm/simulators/flow/FlowGenericProblem_impl.hpp index ac2850f3ed..d61acbde92 100644 --- a/opm/simulators/flow/FlowGenericProblem_impl.hpp +++ b/opm/simulators/flow/FlowGenericProblem_impl.hpp @@ -527,7 +527,12 @@ readBlackoilExtentionsInitialConditions_(std::size_t numDof, micp_.oxygenConcentration.resize(numDof, 0.0); } if (eclState_.fieldProps().has_double("SUREA")) { - micp_.ureaConcentration = getArray(eclState_.fieldProps().get_double("SUREA")); + // dividing by scaling factor 10 (see WellInterfaceGeneric.cpp) + micp_.ureaConcentration.resize(numDof, 0.1); + const auto& urea = getArray(eclState_.fieldProps().get_double("SUREA")); + for (std::size_t xIdx = 0; xIdx < numDof; ++xIdx) { + micp_.ureaConcentration[xIdx] *= urea[xIdx]; + } } else { micp_.ureaConcentration.resize(numDof, 0.0); } diff --git a/opm/simulators/utils/phaseUsageFromDeck.hpp b/opm/simulators/utils/phaseUsageFromDeck.hpp index 473982a5f7..245cba3e66 100644 --- a/opm/simulators/utils/phaseUsageFromDeck.hpp +++ b/opm/simulators/utils/phaseUsageFromDeck.hpp @@ -40,7 +40,7 @@ namespace Opm /// to determine active phases. PhaseUsage phaseUsageFromDeck(const Deck& deck); - /// Handels MICP + /// Handles MICP PhaseUsage phaseUsageMICP(const Phases& phases); } diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index 302d73de08..a16e73f30c 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -797,7 +797,7 @@ Scalar WellInterfaceGeneric::wurea_() const if (injectorType == InjectorType::WATER) { WellMICPProperties urea = this->well_ecl_.getMICPProperties(); - const Scalar urea_injection_concentration = urea.m_ureaConcentration / 10;//Dividing by scaling factor 10 + const Scalar urea_injection_concentration = urea.m_ureaConcentration / 10; return urea_injection_concentration; } else { // Not a water injection well => no urea. From 1031e993b46ad75639e0666f8eb0113ebdc0d109 Mon Sep 17 00:00:00 2001 From: David Landa Marban Date: Wed, 2 Apr 2025 15:16:45 +0200 Subject: [PATCH 3/3] Addressing totto82 comments --- .../blackoil/blackoilboundaryratevector.hh | 2 + .../blackoil/blackoildiffusionmodule.hh | 226 +++++++++--------- .../blackoil/blackoildispersionmodule.hh | 131 +++++----- .../blackoil/blackoilintensivequantities.hh | 7 +- opm/models/blackoil/blackoillocalresidual.hh | 2 +- .../blackoil/blackoillocalresidualtpfa.hh | 15 +- opm/models/blackoil/blackoilmicpmodules.hh | 54 +++-- opm/models/blackoil/blackoilmodel.hh | 13 + opm/models/blackoil/blackoilnewtonmethod.hpp | 12 +- opm/models/blackoil/blackoilproperties.hh | 3 + opm/models/io/vtkblackoilmicpmodule.hpp | 3 +- .../flow/FlowGenericProblem_impl.hpp | 7 +- opm/simulators/flow/FlowProblemBlackoil.hpp | 8 +- opm/simulators/flow/OutputBlackoilModule.hpp | 9 +- opm/simulators/utils/BlackoilPhases.cpp | 2 +- opm/simulators/utils/BlackoilPhases.hpp | 6 +- opm/simulators/utils/phaseUsageFromDeck.cpp | 27 +-- opm/simulators/utils/phaseUsageFromDeck.hpp | 3 - .../wells/StandardWellConnections.cpp | 3 +- opm/simulators/wells/WellInterfaceGeneric.cpp | 7 +- 20 files changed, 278 insertions(+), 262 deletions(-) diff --git a/opm/models/blackoil/blackoilboundaryratevector.hh b/opm/models/blackoil/blackoilboundaryratevector.hh index 60c586801c..97731ff365 100644 --- a/opm/models/blackoil/blackoilboundaryratevector.hh +++ b/opm/models/blackoil/blackoilboundaryratevector.hh @@ -182,6 +182,8 @@ public: (*this)[Indices::contiMicrobialEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.microbialConcentration(); (*this)[Indices::contiOxygenEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.oxygenConcentration(); (*this)[Indices::contiUreaEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.ureaConcentration(); + // since the urea concentration can be much larger than 1, then we apply a scaling factor + (*this)[Indices::contiUreaEqIdx] *= getPropValue(); } // make sure that the right mass conservation quantities are used diff --git a/opm/models/blackoil/blackoildiffusionmodule.hh b/opm/models/blackoil/blackoildiffusionmodule.hh index bb52c0bb31..48308d0e27 100644 --- a/opm/models/blackoil/blackoildiffusionmodule.hh +++ b/opm/models/blackoil/blackoildiffusionmodule.hh @@ -184,95 +184,94 @@ public: const auto& exFs = exIq.fluidState(); Evaluation diffR = 0.0; if constexpr(enableMICP) { - // The diffusion coefficients are given for mass concentrations (urea scaled by a factor of 10) + // The diffusion coefficients are given for mass concentrations Evaluation bAvg = (inFs.invB(waterPhaseIdx) + Toolbox::value(exFs.invB(waterPhaseIdx))) / 2; diffR = inIq.microbialConcentration() - Toolbox::value(exIq.microbialConcentration()); flux[contiMicrobialEqIdx] += bAvg * diffR * diffusivity - * effectiveDiffusionCoefficient[waterPhaseIdx][0]; + * effectiveDiffusionCoefficient[waterPhaseIdx][contiMicrobialEqIdx - 1]; diffR = inIq.oxygenConcentration() - Toolbox::value(exIq.oxygenConcentration()); flux[contiOxygenEqIdx] += bAvg * diffR * diffusivity - * effectiveDiffusionCoefficient[waterPhaseIdx][1]; + * effectiveDiffusionCoefficient[waterPhaseIdx][contiOxygenEqIdx - 1]; diffR = inIq.ureaConcentration() - Toolbox::value(exIq.ureaConcentration()); - // dividing by scaling factor 10 (see WellInterfaceGeneric.cpp) flux[contiUreaEqIdx] += bAvg * diffR * diffusivity - * effectiveDiffusionCoefficient[waterPhaseIdx][2] / 10; + * effectiveDiffusionCoefficient[waterPhaseIdx][contiUreaEqIdx - 1]; + return; } - else { - unsigned pvtRegionIndex = inFs.pvtRegionIndex(); - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - // no diffusion in water for blackoil models - if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { - continue; - } - - // no diffusion in gas phase in water + gas system. - if (FluidSystem::gasPhaseIdx == phaseIdx && !FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - continue; - } - - // arithmetic mean of the phase's b factor weighed by saturation - Evaluation bSAvg = inFs.saturation(phaseIdx) * inFs.invB(phaseIdx); - bSAvg += Toolbox::value(exFs.saturation(phaseIdx)) * Toolbox::value(exFs.invB(phaseIdx)); - bSAvg /= 2; - - // phase not present, skip - if(bSAvg < 1.0e-6) - continue; - Evaluation convFactor = 1.0; - if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) { - Evaluation rsAvg = (inFs.Rs() + Toolbox::value(exFs.Rs())) / 2; - convFactor = 1.0 / (toFractionGasOil(pvtRegionIndex) + rsAvg); - diffR = inFs.Rs() - Toolbox::value(exFs.Rs()); - } - if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) { - Evaluation rvAvg = (inFs.Rv() + Toolbox::value(exFs.Rv())) / 2; - convFactor = toFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toFractionGasOil(pvtRegionIndex)); - diffR = inFs.Rv() - Toolbox::value(exFs.Rv()); - } - if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) { - Evaluation rsAvg = (inFs.Rsw() + Toolbox::value(exFs.Rsw())) / 2; - convFactor = 1.0 / (toFractionGasWater(pvtRegionIndex) + rsAvg); - diffR = inFs.Rsw() - Toolbox::value(exFs.Rsw()); - } - if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) { - Evaluation rvAvg = (inFs.Rvw() + Toolbox::value(exFs.Rvw())) / 2; - convFactor = toFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toFractionGasWater(pvtRegionIndex)); - diffR = inFs.Rvw() - Toolbox::value(exFs.Rvw()); - } - - // mass flux of solvent component (oil in oil or gas in gas) - unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx); - unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx); - flux[conti0EqIdx + activeSolventCompIdx] += - - bSAvg - * convFactor - * diffR - * diffusivity - * effectiveDiffusionCoefficient[phaseIdx][solventCompIdx]; - - // mass flux of solute component (gas in oil or oil in gas) - unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx); - unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx); - flux[conti0EqIdx + activeSoluteCompIdx] += - bSAvg - * diffR - * convFactor - * diffusivity - * effectiveDiffusionCoefficient[phaseIdx][soluteCompIdx]; + unsigned pvtRegionIndex = inFs.pvtRegionIndex(); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; } + + // no diffusion in water for blackoil models + if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { + continue; + } + + // no diffusion in gas phase in water + gas system. + if (FluidSystem::gasPhaseIdx == phaseIdx && !FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + continue; + } + + // arithmetic mean of the phase's b factor weighed by saturation + Evaluation bSAvg = inFs.saturation(phaseIdx) * inFs.invB(phaseIdx); + bSAvg += Toolbox::value(exFs.saturation(phaseIdx)) * Toolbox::value(exFs.invB(phaseIdx)); + bSAvg /= 2; + + // phase not present, skip + if(bSAvg < 1.0e-6) + continue; + Evaluation convFactor = 1.0; + if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) { + Evaluation rsAvg = (inFs.Rs() + Toolbox::value(exFs.Rs())) / 2; + convFactor = 1.0 / (toFractionGasOil(pvtRegionIndex) + rsAvg); + diffR = inFs.Rs() - Toolbox::value(exFs.Rs()); + } + if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) { + Evaluation rvAvg = (inFs.Rv() + Toolbox::value(exFs.Rv())) / 2; + convFactor = toFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toFractionGasOil(pvtRegionIndex)); + diffR = inFs.Rv() - Toolbox::value(exFs.Rv()); + } + if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) { + Evaluation rsAvg = (inFs.Rsw() + Toolbox::value(exFs.Rsw())) / 2; + convFactor = 1.0 / (toFractionGasWater(pvtRegionIndex) + rsAvg); + diffR = inFs.Rsw() - Toolbox::value(exFs.Rsw()); + } + if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) { + Evaluation rvAvg = (inFs.Rvw() + Toolbox::value(exFs.Rvw())) / 2; + convFactor = toFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toFractionGasWater(pvtRegionIndex)); + diffR = inFs.Rvw() - Toolbox::value(exFs.Rvw()); + } + + // mass flux of solvent component (oil in oil or gas in gas) + unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx); + unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx); + flux[conti0EqIdx + activeSolventCompIdx] += + - bSAvg + * convFactor + * diffR + * diffusivity + * effectiveDiffusionCoefficient[phaseIdx][solventCompIdx]; + + // mass flux of solute component (gas in oil or oil in gas) + unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx); + unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx); + flux[conti0EqIdx + activeSoluteCompIdx] += + bSAvg + * diffR + * convFactor + * diffusivity + * effectiveDiffusionCoefficient[phaseIdx][soluteCompIdx]; } } @@ -376,9 +375,11 @@ class BlackOilDiffusionIntensiveQuantities using FluidSystem = GetPropType; using IntensiveQuantities = GetPropType; using MICPModule = BlackOilMICPModule; + enum { numPhases = FluidSystem::numPhases }; enum { numComponents = FluidSystem::numComponents }; enum { enableMICP = getPropValue() }; + static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx; public: @@ -462,37 +463,37 @@ protected: diffusionCoefficient_[waterPhaseIdx][0] = MICPModule::microbialDiffusion(pvtRegionIndex); diffusionCoefficient_[waterPhaseIdx][1] = MICPModule::oxygenDiffusion(pvtRegionIndex); diffusionCoefficient_[waterPhaseIdx][2] = MICPModule::ureaDiffusion(pvtRegionIndex); + return; } - else { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - // no diffusion in water for blackoil models - if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { - continue; - } + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } - // Based on Millington, R. J., & Quirk, J. P. (1961). - // \Note: it is possible to use NumericalConstants later - // constexpr auto& numconst = GetPropValue; - constexpr double myeps = 0.0001; //numconst.blackoildiffusionmoduleeps; - const Evaluation& base = - Toolbox::max(myeps, //0.0001, - intQuants.porosity() - * intQuants.fluidState().saturation(phaseIdx)); - tortuosity_[phaseIdx] = - 1.0 / (intQuants.porosity() * intQuants.porosity()) - * Toolbox::pow(base, 10.0/3.0); + // no diffusion in water for blackoil models + if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { + continue; + } - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - diffusionCoefficient_[phaseIdx][compIdx] = - FluidSystem::diffusionCoefficient(fluidState, - paramCache, - phaseIdx, - compIdx); - } + // Based on Millington, R. J., & Quirk, J. P. (1961). + // \Note: it is possible to use NumericalConstants later + // constexpr auto& numconst = GetPropValue; + constexpr double myeps = 0.0001; //numconst.blackoildiffusionmoduleeps; + const Evaluation& base = + Toolbox::max(myeps, //0.0001, + intQuants.porosity() + * intQuants.fluidState().saturation(phaseIdx)); + tortuosity_[phaseIdx] = + 1.0 / (intQuants.porosity() * intQuants.porosity()) + * Toolbox::pow(base, 10.0/3.0); + + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + diffusionCoefficient_[phaseIdx][compIdx] = + FluidSystem::diffusionCoefficient(fluidState, + paramCache, + phaseIdx, + compIdx); } } } @@ -583,6 +584,7 @@ class BlackOilDiffusionExtensiveQuantities enum { numPhases = getPropValue() }; enum { numComponents = getPropValue() }; enum { enableMICP = getPropValue() }; + static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx; using DimVector = Dune::FieldVector; @@ -626,22 +628,22 @@ public: ( intQuantsInside.effectiveDiffusionCoefficient(waterPhaseIdx, compIdx) + intQuantsOutside.effectiveDiffusionCoefficient(waterPhaseIdx, compIdx) ); } + return; } - else { - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - // no diffusion in water for blackoil models - if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { - continue; - } - for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - effectiveDiffusionCoefficient[phaseIdx][compIdx] = 0.5 * - ( intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx) + - intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx) ); - Valgrind::CheckDefined(effectiveDiffusionCoefficient[phaseIdx][compIdx]); - } + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + // no diffusion in water for blackoil models + if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { + continue; + } + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + effectiveDiffusionCoefficient[phaseIdx][compIdx] = 0.5 * + ( intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx) + + intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx) ); + Valgrind::CheckDefined(effectiveDiffusionCoefficient[phaseIdx][compIdx]); } } } diff --git a/opm/models/blackoil/blackoildispersionmodule.hh b/opm/models/blackoil/blackoildispersionmodule.hh index 70b5554c7a..e16329b5aa 100644 --- a/opm/models/blackoil/blackoildispersionmodule.hh +++ b/opm/models/blackoil/blackoildispersionmodule.hh @@ -208,79 +208,78 @@ public: * dispersivity * diffR; diffR = inIq.ureaConcentration()- Toolbox::value(exIq.ureaConcentration()); - // dividing by scaling factor 10 (see WellInterfaceGeneric.cpp) flux[contiUreaEqIdx] += bAvg * normVelocityAvg[waterPhaseIdx] * dispersivity - * diffR / 10; + * diffR; + return; } - else { - unsigned pvtRegionIndex = inFs.pvtRegionIndex(); - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - // no dispersion in water for blackoil models unless water can contain dissolved gas - if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { - continue; - } - - // adding dispersion in the gas phase leads to - // convergence issues and unphysical results. - // we disable dispersion in the gas phase for now - if (FluidSystem::gasPhaseIdx == phaseIdx) { - continue; - } - - // arithmetic mean of the phase's b factor - Evaluation bAvg = inFs.invB(phaseIdx); - bAvg += Toolbox::value(exFs.invB(phaseIdx)); - bAvg /= 2; - - Evaluation convFactor = 1.0; - if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) { - Evaluation rsAvg = (inFs.Rs() + Toolbox::value(exFs.Rs())) / 2; - convFactor = 1.0 / (toMassFractionGasOil(pvtRegionIndex) + rsAvg); - diffR = inFs.Rs() - Toolbox::value(exFs.Rs()); - } - if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) { - Evaluation rvAvg = (inFs.Rv() + Toolbox::value(exFs.Rv())) / 2; - convFactor = toMassFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toMassFractionGasOil(pvtRegionIndex)); - diffR = inFs.Rv() - Toolbox::value(exFs.Rv()); - } - if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) { - Evaluation rsAvg = (inFs.Rsw() + Toolbox::value(exFs.Rsw())) / 2; - convFactor = 1.0 / (toMassFractionGasWater(pvtRegionIndex) + rsAvg); - diffR = inFs.Rsw() - Toolbox::value(exFs.Rsw()); - } - if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) { - Evaluation rvAvg = (inFs.Rvw() + Toolbox::value(exFs.Rvw())) / 2; - convFactor = toMassFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toMassFractionGasWater(pvtRegionIndex)); - diffR = inFs.Rvw() - Toolbox::value(exFs.Rvw()); - } - - // mass flux of solvent component - unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx); - unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx); - flux[conti0EqIdx + activeSolventCompIdx] += - - bAvg - * normVelocityAvg[phaseIdx] - * convFactor - * dispersivity - * diffR; - - // mass flux of solute component - unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx); - unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx); - flux[conti0EqIdx + activeSoluteCompIdx] += - bAvg - * normVelocityAvg[phaseIdx] - * convFactor - * dispersivity - * diffR; + unsigned pvtRegionIndex = inFs.pvtRegionIndex(); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; } + + // no dispersion in water for blackoil models unless water can contain dissolved gas + if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) { + continue; + } + + // adding dispersion in the gas phase leads to + // convergence issues and unphysical results. + // we disable dispersion in the gas phase for now + if (FluidSystem::gasPhaseIdx == phaseIdx) { + continue; + } + + // arithmetic mean of the phase's b factor + Evaluation bAvg = inFs.invB(phaseIdx); + bAvg += Toolbox::value(exFs.invB(phaseIdx)); + bAvg /= 2; + + Evaluation convFactor = 1.0; + if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) { + Evaluation rsAvg = (inFs.Rs() + Toolbox::value(exFs.Rs())) / 2; + convFactor = 1.0 / (toMassFractionGasOil(pvtRegionIndex) + rsAvg); + diffR = inFs.Rs() - Toolbox::value(exFs.Rs()); + } + if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) { + Evaluation rvAvg = (inFs.Rv() + Toolbox::value(exFs.Rv())) / 2; + convFactor = toMassFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toMassFractionGasOil(pvtRegionIndex)); + diffR = inFs.Rv() - Toolbox::value(exFs.Rv()); + } + if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) { + Evaluation rsAvg = (inFs.Rsw() + Toolbox::value(exFs.Rsw())) / 2; + convFactor = 1.0 / (toMassFractionGasWater(pvtRegionIndex) + rsAvg); + diffR = inFs.Rsw() - Toolbox::value(exFs.Rsw()); + } + if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) { + Evaluation rvAvg = (inFs.Rvw() + Toolbox::value(exFs.Rvw())) / 2; + convFactor = toMassFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toMassFractionGasWater(pvtRegionIndex)); + diffR = inFs.Rvw() - Toolbox::value(exFs.Rvw()); + } + + // mass flux of solvent component + unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx); + unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx); + flux[conti0EqIdx + activeSolventCompIdx] += + - bAvg + * normVelocityAvg[phaseIdx] + * convFactor + * dispersivity + * diffR; + + // mass flux of solute component + unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx); + unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx); + flux[conti0EqIdx + activeSoluteCompIdx] += + bAvg + * normVelocityAvg[phaseIdx] + * convFactor + * dispersivity + * diffR; } } diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index de87a64747..bb6824d886 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -510,9 +510,10 @@ public: // deal with MICP if constexpr (enableMICP){ - Evaluation biofilm_ = priVars.makeEvaluation(Indices::biofilmConcentrationIdx, timeIdx, linearizationType); - Evaluation calcite_ = priVars.makeEvaluation(Indices::calciteConcentrationIdx, timeIdx, linearizationType); - porosity_ += max(- biofilm_ - calcite_, -referencePorosity_+1e-8); + Evaluation biofilm_ = priVars.makeEvaluation(Indices::biofilmConcentrationIdx, timeIdx, linearizationType); + Evaluation calcite_ = priVars.makeEvaluation(Indices::calciteConcentrationIdx, timeIdx, linearizationType); + // minimum porosity of 1e-8 to prevent numerical issues + porosity_ -= min(biofilm_ + calcite_, referencePorosity_ - 1e-8); } // deal with salt-precipitation diff --git a/opm/models/blackoil/blackoillocalresidual.hh b/opm/models/blackoil/blackoillocalresidual.hh index 44110bc86a..4ff2a7664c 100644 --- a/opm/models/blackoil/blackoillocalresidual.hh +++ b/opm/models/blackoil/blackoillocalresidual.hh @@ -255,7 +255,7 @@ public: MICPModule::addSource(source, elemCtx, dofIdx, timeIdx); // scale the source term of the energy equation - if (enableEnergy) + if constexpr(enableEnergy) source[Indices::contiEnergyEqIdx] *= getPropValue(); } diff --git a/opm/models/blackoil/blackoillocalresidualtpfa.hh b/opm/models/blackoil/blackoillocalresidualtpfa.hh index 6dfbb9d263..238c13e6c2 100644 --- a/opm/models/blackoil/blackoillocalresidualtpfa.hh +++ b/opm/models/blackoil/blackoillocalresidualtpfa.hh @@ -481,8 +481,6 @@ public: static_assert(!enableBrine, "Relevant computeFlux() method must be implemented for this module before enabling."); // BrineModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx); - - // deal with diffusion (if present). opm-models expects per area flux (added in the tmpdiffusivity). if constexpr(enableDiffusion){ typename DiffusionModule::ExtensiveQuantities::EvaluationArray effectiveDiffusionCoefficient; @@ -496,6 +494,7 @@ public: effectiveDiffusionCoefficient); } + // deal with dispersion (if present). opm-models expects per area flux (added in the tmpdispersivity). if constexpr(enableDispersion){ typename DispersionModule::ExtensiveQuantities::ScalarArray normVelocityAvg; @@ -509,6 +508,11 @@ public: normVelocityAvg); } + + // apply the scaling for the urea equation in MICP + if constexpr (enableMICP) { + MICPModule::applyScaling(flux); + } } template @@ -635,7 +639,6 @@ public: static_assert(!enableSolvent, "Relevant treatment of boundary conditions must be implemented before enabling."); static_assert(!enablePolymer, "Relevant treatment of boundary conditions must be implemented before enabling."); - //static_assert(!enableMICP, "Relevant treatment of boundary conditions must be implemented before enabling."); // make sure that the right mass conservation quantities are used adaptMassConservationQuantities_(bdyFlux, insideIntQuants.pvtRegionIndex()); @@ -697,7 +700,7 @@ public: MICPModule::addSource(source, problem, insideIntQuants, globalSpaceIdex); // scale the source term of the energy equation - if (enableEnergy) + if constexpr(enableEnergy) source[Indices::contiEnergyEqIdx] *= getPropValue(); } @@ -714,7 +717,7 @@ public: MICPModule::addSource(source, problem, insideIntQuants, globalSpaceIdex); // scale the source term of the energy equation - if (enableEnergy) + if constexpr(enableEnergy) source[Indices::contiEnergyEqIdx] *= getPropValue(); } @@ -735,7 +738,7 @@ public: // scale the source term of the energy equation if constexpr(enableEnergy) - source[Indices::contiEnergyEqIdx] *= getPropValue(); + source[Indices::contiEnergyEqIdx] *= getPropValue(); } template diff --git a/opm/models/blackoil/blackoilmicpmodules.hh b/opm/models/blackoil/blackoilmicpmodules.hh index 3ec8c036f8..630191b095 100644 --- a/opm/models/blackoil/blackoilmicpmodules.hh +++ b/opm/models/blackoil/blackoilmicpmodules.hh @@ -135,7 +135,7 @@ public: const auto& fs = intQuants.fluidState(); LhsEval surfaceVolumeWater = Toolbox::template decay(fs.invB(waterPhaseIdx)) * Toolbox::template decay(intQuants.porosity()); - // avoid singular matrix if no water is present. + // avoid singular matrix if no water is present surfaceVolumeWater = max(surfaceVolumeWater, 1e-10); // suspended microbes in water phase const LhsEval massMicrobes = surfaceVolumeWater * Toolbox::template decay(intQuants.microbialConcentration()); @@ -145,10 +145,11 @@ public: const LhsEval massOxygen = surfaceVolumeWater * Toolbox::template decay(intQuants.oxygenConcentration()); LhsEval accumulationOxygen = massOxygen; storage[contiOxygenEqIdx] += accumulationOxygen; - // urea in water phase + // urea in water phase (applying the scaling factor for the urea equation) const LhsEval massUrea = surfaceVolumeWater * Toolbox::template decay(intQuants.ureaConcentration()); LhsEval accumulationUrea = massUrea; storage[contiUreaEqIdx] += accumulationUrea; + storage[contiUreaEqIdx] *= getPropValue(); // biofilm const LhsEval massBiofilm = Toolbox::template decay(intQuants.biofilmConcentration()); LhsEval accumulationBiofilm = massBiofilm; @@ -178,6 +179,14 @@ public: } } + // since the urea concentration can be much larger than 1, then we apply a scaling factor + static void applyScaling(RateVector& flux) + { + if constexpr (enableMICP) { + flux[contiUreaEqIdx] *= getPropValue(); + } + } + static void computeFlux([[maybe_unused]] RateVector& flux, [[maybe_unused]] const ElementContext& elemCtx, [[maybe_unused]] unsigned scvfIdx, @@ -236,29 +245,42 @@ public: Scalar k_str = detachmentRate(satnumIdx); Scalar eta = detachmentExponent(satnumIdx); Scalar k_o = halfVelocityOxygen(satnumIdx); - Scalar k_u = halfVelocityUrea(satnumIdx) / 10;// dividing by scaling factor 10 (see WellInterfaceGeneric.cpp) + Scalar k_u = halfVelocityUrea(satnumIdx); Scalar mu = maximumGrowthRate(satnumIdx); - Scalar mu_u = maximumUreaUtilization(satnumIdx) / 10;// dividing by scaling factor 10 + Scalar mu_u = maximumUreaUtilization(satnumIdx); Scalar Y_sb = yieldGrowthCoefficient(satnumIdx); Scalar F = oxygenConsumptionFactor(satnumIdx); - Scalar Y_uc = yieldUreaToCalciteCoefficient(satnumIdx) * 10;// multiplying by scaling factor 10 + Scalar Y_uc = yieldUreaToCalciteCoefficient(satnumIdx); + + // compute Monod terms (the negative region is replaced by a straight line) + // Schäfer et al (1998) https://doi.org/10.1016/S0169-7722(97)00060-0 + Evaluation k_g = mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()); + Evaluation k_c = mu_u * intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration()); + if (intQuants.oxygenConcentration() < 0) + k_g = mu * intQuants.oxygenConcentration() / k_o; + if (intQuants.ureaConcentration() < 0) + k_c = mu_u * intQuants.ureaConcentration() / k_u; // compute the processes - source[Indices::contiMicrobialEqIdx] += intQuants.microbialConcentration() * intQuants.porosity() * b * - (Y_sb * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()) - k_d - k_a) - + rho_b * intQuants.biofilmConcentration() * k_str * pow(normVelocityCell, eta); + source[Indices::contiMicrobialEqIdx] += intQuants.microbialConcentration() * intQuants.porosity() * + b * (Y_sb * k_g - k_d - k_a) + + rho_b * intQuants.biofilmConcentration() * k_str * pow(normVelocityCell, eta); - source[Indices::contiOxygenEqIdx] -= (intQuants.microbialConcentration() * intQuants.porosity() * b + rho_b * intQuants.biofilmConcentration()) * - F * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()); + source[Indices::contiOxygenEqIdx] -= (intQuants.microbialConcentration() * intQuants.porosity() * + b + rho_b * intQuants.biofilmConcentration()) * F * k_g; - source[Indices::contiUreaEqIdx] -= rho_b * intQuants.biofilmConcentration() * mu_u * intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration()); + source[Indices::contiUreaEqIdx] -= rho_b * intQuants.biofilmConcentration() * k_c; - source[Indices::contiBiofilmEqIdx] += intQuants.biofilmConcentration() * (Y_sb * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()) - k_d - - k_str * pow(normVelocityCell, eta) - Y_uc * (rho_b / rho_c) * intQuants.biofilmConcentration() * mu_u * - (intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration())) / (intQuants.porosity() + intQuants.biofilmConcentration())) - + k_a * intQuants.microbialConcentration() * intQuants.porosity() * b / rho_b; + source[Indices::contiBiofilmEqIdx] += intQuants.biofilmConcentration() * (Y_sb * k_g - k_d - + k_str * pow(normVelocityCell, eta) - Y_uc * (rho_b / rho_c) * + intQuants.biofilmConcentration() * k_c / (intQuants.porosity() + + intQuants.biofilmConcentration())) + k_a * intQuants.microbialConcentration() * + intQuants.porosity() * b / rho_b; - source[Indices::contiCalciteEqIdx] += (rho_b / rho_c) * intQuants.biofilmConcentration() * Y_uc * mu_u * intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration()); + source[Indices::contiCalciteEqIdx] += (rho_b / rho_c) * intQuants.biofilmConcentration() * Y_uc * k_c; + + // since the urea concentration can be much larger than 1, then we apply a scaling factor + source[Indices::contiUreaEqIdx] *= getPropValue(); } } diff --git a/opm/models/blackoil/blackoilmodel.hh b/opm/models/blackoil/blackoilmodel.hh index 0f875f7aab..4a562c0aa2 100644 --- a/opm/models/blackoil/blackoilmodel.hh +++ b/opm/models/blackoil/blackoilmodel.hh @@ -197,6 +197,19 @@ public: static constexpr Scalar value = 1.0/(30.0*4184.0*alpha); }; +//! similarly to the energy equation, the urea equation in MICP is scaled by a factor of 0.1 +template +struct BlackOilUreaScalingFactor +{ +private: + using Scalar = GetPropType; + static constexpr Scalar alpha = getPropValue() ? 1000.0 : 1.0; + +public: + using type = Scalar; + static constexpr Scalar value = 1.0/(10.0*alpha); +}; + // by default, ebos formulates the conservation equations in terms of mass not surface // volumes template diff --git a/opm/models/blackoil/blackoilnewtonmethod.hpp b/opm/models/blackoil/blackoilnewtonmethod.hpp index 19bac6eb4e..716d466cad 100644 --- a/opm/models/blackoil/blackoilnewtonmethod.hpp +++ b/opm/models/blackoil/blackoilnewtonmethod.hpp @@ -372,9 +372,11 @@ protected: nextValue[pvIdx] = std::clamp(nextValue[pvIdx], bparams_.pressMin_, bparams_.pressMax_); } - // keep the values above 0 - // for the biofilm and calcite, we set an upper limit + // for the biofilm and calcite, we set an upper limit equal to the initial porosity + // minus 1e-8. This prevents singularities (e.g., one of the calcite source term is + // evaluated at 1/(iniPoro - calcite)). The value 1e-8 is taken from the salt precipitation + // clapping above. if constexpr (enableMICP) { if (pvIdx == Indices::microbialConcentrationIdx) { nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0}); @@ -386,10 +388,12 @@ protected: nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0}); } if (pvIdx == Indices::biofilmConcentrationIdx) { - nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, this->problem().referencePorosity(globalDofIdx, 0) - 1.e-8); + nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, + this->problem().referencePorosity(globalDofIdx, 0) - 1e-8); } if (pvIdx == Indices::calciteConcentrationIdx) { - nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, this->problem().referencePorosity(globalDofIdx, 0) - 1.e-8); + nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, + this->problem().referencePorosity(globalDofIdx, 0) - 1e-8); } } } diff --git a/opm/models/blackoil/blackoilproperties.hh b/opm/models/blackoil/blackoilproperties.hh index 3094741656..1a68cb80f4 100644 --- a/opm/models/blackoil/blackoilproperties.hh +++ b/opm/models/blackoil/blackoilproperties.hh @@ -90,6 +90,9 @@ struct EnableMech { using type = UndefinedProperty; }; template struct BlackOilEnergyScalingFactor { using type = UndefinedProperty; }; +//! Similarly to the energy equation, a scaling is applied to the urea equation in MICP +template +struct BlackOilUreaScalingFactor { using type = UndefinedProperty; }; } // namespace Opm::Properties diff --git a/opm/models/io/vtkblackoilmicpmodule.hpp b/opm/models/io/vtkblackoilmicpmodule.hpp index 713f294b1c..6e74b99f7b 100644 --- a/opm/models/io/vtkblackoilmicpmodule.hpp +++ b/opm/models/io/vtkblackoilmicpmodule.hpp @@ -141,9 +141,8 @@ public: } if (params_.ureaConcentrationOutput_) { - // multiplying by scaling factor 10 (see WellInterfaceGeneric.cpp) ureaConcentration_[globalDofIdx] = - 10 * scalarValue(intQuants.ureaConcentration()); + scalarValue(intQuants.ureaConcentration()); } if (params_.biofilmConcentrationOutput_) { diff --git a/opm/simulators/flow/FlowGenericProblem_impl.hpp b/opm/simulators/flow/FlowGenericProblem_impl.hpp index d61acbde92..ac2850f3ed 100644 --- a/opm/simulators/flow/FlowGenericProblem_impl.hpp +++ b/opm/simulators/flow/FlowGenericProblem_impl.hpp @@ -527,12 +527,7 @@ readBlackoilExtentionsInitialConditions_(std::size_t numDof, micp_.oxygenConcentration.resize(numDof, 0.0); } if (eclState_.fieldProps().has_double("SUREA")) { - // dividing by scaling factor 10 (see WellInterfaceGeneric.cpp) - micp_.ureaConcentration.resize(numDof, 0.1); - const auto& urea = getArray(eclState_.fieldProps().get_double("SUREA")); - for (std::size_t xIdx = 0; xIdx < numDof; ++xIdx) { - micp_.ureaConcentration[xIdx] *= urea[xIdx]; - } + micp_.ureaConcentration = getArray(eclState_.fieldProps().get_double("SUREA")); } else { micp_.ureaConcentration.resize(numDof, 0.0); } diff --git a/opm/simulators/flow/FlowProblemBlackoil.hpp b/opm/simulators/flow/FlowProblemBlackoil.hpp index e536c47e65..98319b37cf 100644 --- a/opm/simulators/flow/FlowProblemBlackoil.hpp +++ b/opm/simulators/flow/FlowProblemBlackoil.hpp @@ -621,8 +621,7 @@ public: if constexpr (enableMICP) { rate[Indices::microbialConcentrationIdx] += source.rate({ijk, SourceComponent::MICR}) / this->model().dofTotalVolume(globalDofIdx); rate[Indices::oxygenConcentrationIdx] += source.rate({ijk, SourceComponent::OXYG}) / this->model().dofTotalVolume(globalDofIdx); - // dividing by scaling factor 10 (see WellInterfaceGeneric.cpp) - rate[Indices::ureaConcentrationIdx] += source.rate({ijk, SourceComponent::UREA}) / (this->model().dofTotalVolume(globalDofIdx) * 10); + rate[Indices::ureaConcentrationIdx] += source.rate({ijk, SourceComponent::UREA}) / (this->model().dofTotalVolume(globalDofIdx)); } if constexpr (enableEnergy) { for (unsigned i = 0; i < phidx_map.size(); ++i) { @@ -1542,8 +1541,9 @@ protected: if constexpr (!enableMICP) throw std::logic_error("MICP is disabled and you're trying to add urea to BC"); - // dividing by scaling factor 10 (see WellInterfaceGeneric.cpp) - rate[Indices::ureaConcentrationIdx] = bc.rate / 10; + rate[Indices::ureaConcentrationIdx] = bc.rate; + // since the urea concentration can be much larger than 1, then we apply a scaling factor + rate[Indices::ureaConcentrationIdx] *= getPropValue(); } void updateExplicitQuantities_(const bool first_step_after_restart) diff --git a/opm/simulators/flow/OutputBlackoilModule.hpp b/opm/simulators/flow/OutputBlackoilModule.hpp index 57bd5ec4c3..d33e1e48f9 100644 --- a/opm/simulators/flow/OutputBlackoilModule.hpp +++ b/opm/simulators/flow/OutputBlackoilModule.hpp @@ -1036,8 +1036,7 @@ private: const IntensiveQuantities& intQuants, const double surfVolWat) { - // rescaling back the urea concentration (see WellInterfaceGeneric.cpp) - const Scalar mass = surfVolWat * intQuants.ureaConcentration().value() * 10; + const Scalar mass = surfVolWat * intQuants.ureaConcentration().value(); this->fipC_.assignUreaMass(globalDofIdx, mass); } @@ -1512,8 +1511,7 @@ private: micpC.assign(ectx.globalDofIdx, ectx.intQuants.microbialConcentration().value(), ectx.intQuants.oxygenConcentration().value(), - // rescaling back the urea concentration (see WellInterfaceGeneric.cpp) - 10 * ectx.intQuants.ureaConcentration().value(), + ectx.intQuants.ureaConcentration().value(), ectx.intQuants.biofilmConcentration().value(), ectx.intQuants.calciteConcentration().value()); }, this->micpC_.allocated() @@ -2050,10 +2048,9 @@ private: Entry{ScalarEntry{"BMUIP", [&model = this->simulator_.model()](const Context& ectx) { - // rescaling back the urea concentration (see WellInterfaceGeneric.cpp) return getValue(ectx.intQuants.ureaConcentration()) * getValue(ectx.intQuants.porosity()) * - model.dofTotalVolume(ectx.globalDofIdx) * 10; + model.dofTotalVolume(ectx.globalDofIdx) * 1; } } }, diff --git a/opm/simulators/utils/BlackoilPhases.cpp b/opm/simulators/utils/BlackoilPhases.cpp index 70b58807fe..6ca942c7c2 100644 --- a/opm/simulators/utils/BlackoilPhases.cpp +++ b/opm/simulators/utils/BlackoilPhases.cpp @@ -55,5 +55,5 @@ Opm::PhaseUsage::PhaseUsage(std::vector phases) this->has_foam = this->phase_used[BlackoilPhases::Foam]; this->has_brine = this->phase_used[BlackoilPhases::Brine]; this->has_zFraction = this->phase_used[BlackoilPhases::ZFraction]; - this->has_micp = this->phase_used[BlackoilPhases::MICP]; + this->has_micp = false; } diff --git a/opm/simulators/utils/BlackoilPhases.hpp b/opm/simulators/utils/BlackoilPhases.hpp index 8a99423a07..848484013f 100644 --- a/opm/simulators/utils/BlackoilPhases.hpp +++ b/opm/simulators/utils/BlackoilPhases.hpp @@ -36,10 +36,10 @@ namespace Opm // to and from active ones. That said, they are not considered by num_phases or // MaxNumPhases. The crypto phases which are currently implemented are solvent, // polymer, energy, polymer molecular weight, foam and brine. - static const int NumCryptoPhases = 8; + static const int NumCryptoPhases = 7; // enum ComponentIndex { Water = 0, Oil = 1, Gas = 2 }; - enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2, Solvent = 3, Polymer = 4, Energy = 5, PolymerMW = 6, Foam = 7, Brine = 8, ZFraction = 9, MICP = 10 }; + enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2, Solvent = 3, Polymer = 4, Energy = 5, PolymerMW = 6, Foam = 7, Brine = 8, ZFraction = 9 }; }; struct PhaseUsage : public BlackoilPhases @@ -66,4 +66,4 @@ namespace Opm } // namespace Opm -#endif // OPM_BLACKOILPHASES_HEADER_INCLUDED +#endif // OPM_BLACKOILPHASES_HEADER_INCLUDED \ No newline at end of file diff --git a/opm/simulators/utils/phaseUsageFromDeck.cpp b/opm/simulators/utils/phaseUsageFromDeck.cpp index 8cdfb21db2..732360695e 100644 --- a/opm/simulators/utils/phaseUsageFromDeck.cpp +++ b/opm/simulators/utils/phaseUsageFromDeck.cpp @@ -134,33 +134,18 @@ PhaseUsage phaseUsage(const Phases& phases) return pu; } -PhaseUsage phaseUsageMICP() -{ - PhaseUsage pu; - pu.phase_used.fill(0); - for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases + BlackoilPhases::NumCryptoPhases; ++phaseIdx) { - pu.phase_used[phaseIdx] = 0; - pu.phase_pos[phaseIdx] = -1; - } - pu.phase_used[BlackoilPhases::Aqua] = 1; - pu.phase_used[BlackoilPhases::MICP] = 1; - pu.phase_pos[BlackoilPhases::Aqua] = 0; - pu.phase_pos[BlackoilPhases::MICP] = 1; - pu.num_phases = 1; - pu.has_micp = true; - - return pu; -} - PhaseUsage phaseUsageFromDeck(const EclipseState& eclipseState) { const auto& phases = eclipseState.runspec().phases(); if (eclipseState.runspec().micp()) { - return phaseUsageMICP(); + PhaseUsage pu = phaseUsage(phases); + pu.has_micp = true; + return pu; + } + else { + return phaseUsage(phases); } - - return phaseUsage(phases); } /// Looks at presence of WATER, OIL and GAS keywords in deck diff --git a/opm/simulators/utils/phaseUsageFromDeck.hpp b/opm/simulators/utils/phaseUsageFromDeck.hpp index 245cba3e66..ae7403c600 100644 --- a/opm/simulators/utils/phaseUsageFromDeck.hpp +++ b/opm/simulators/utils/phaseUsageFromDeck.hpp @@ -39,9 +39,6 @@ namespace Opm /// Looks at presence of WATER, OIL and GAS keywords in deck /// to determine active phases. PhaseUsage phaseUsageFromDeck(const Deck& deck); - - /// Handles MICP - PhaseUsage phaseUsageMICP(const Phases& phases); } #endif // OPM_PHASEUSAGEFROMDECK_HEADER_INCLUDED diff --git a/opm/simulators/wells/StandardWellConnections.cpp b/opm/simulators/wells/StandardWellConnections.cpp index 1de062ba09..c3c51b750f 100644 --- a/opm/simulators/wells/StandardWellConnections.cpp +++ b/opm/simulators/wells/StandardWellConnections.cpp @@ -790,8 +790,7 @@ connectionRatesMICP(Scalar& rate_m, cq_s_urea *= std::get(ureaConcentration); } - // rescaling back the urea concentration (see WellInterfaceGeneric.cpp) - rate_u = cq_s_urea.value() * 10; + rate_u = cq_s_urea.value(); return {well_.restrictEval(cq_s_microbe), well_.restrictEval(cq_s_oxygen), diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index a16e73f30c..c1f03c61e0 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -785,11 +785,6 @@ Scalar WellInterfaceGeneric::wpolymer_() const } } -// The urea injection concentration is scaled down by a factor of 10, since its value -// can be much bigger than 1 (not doing this slows the simulations). The corresponding -// values are scaled accordingly when computing the reactions and also when writing -// the output files. - template Scalar WellInterfaceGeneric::wurea_() const { @@ -797,7 +792,7 @@ Scalar WellInterfaceGeneric::wurea_() const if (injectorType == InjectorType::WATER) { WellMICPProperties urea = this->well_ecl_.getMICPProperties(); - const Scalar urea_injection_concentration = urea.m_ureaConcentration / 10; + const Scalar urea_injection_concentration = urea.m_ureaConcentration; return urea_injection_concentration; } else { // Not a water injection well => no urea.