diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index ad147010f22..bd3e4c672f4 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -494,6 +494,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_glift1.cpp tests/test_graphcoloring.cpp tests/test_GroupState.cpp + tests/test_impesweights.cpp tests/test_injection_topup_phase_validation.cpp tests/test_interregflows.cpp tests/test_invert.cpp diff --git a/doc/man1/flow.1 b/doc/man1/flow.1 index f2acf13bb4c..b1a4051fdd2 100644 --- a/doc/man1/flow.1 +++ b/doc/man1/flow.1 @@ -176,7 +176,7 @@ Injection multiplier oscillation threshold (used for multiplier dampening). Defa Set compatibility mode for the SKIP100/SKIP300 keywords. Options are 100 (skip SKIP100..ENDSKIP, keep SKIP300..ENDSKIP) [default], 300 (skip SKIP300..ENDSKIP, keep SKIP100..ENDSKIP) and all (skip both SKIP100..ENDSKIP and SKIP300..ENDSKIP) . Default: "100" .TP \fB\-\-linear\-solver\fR=\fI\,STRING\/\fR -Configuration of solver. Valid options are: cprw (default), ilu0, dilu, cpr (an alias for cprw), cpr_quasiimpes, cpr_trueimpes, cpr_trueimpesanalytic, amg or hybrid (experimental). Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, ending with '.json.'. Default: "cprw" +Configuration of solver. Valid options are: cprw (default), ilu0, dilu, cpr (an alias for cprw), cpr_quasiimpes, cpr_trueimpes, cpr_trueimpesanalytic, cpr_coatsblackoil, cpr_coatscompositional, amg or hybrid (experimental). Alternatively, you can request a configuration to be read from a JSON file by giving the filename here, ending with '.json.'. Default: "cprw" .TP \fB\-\-linear\-solver\-ignore\-convergence\-failure\fR=\fI\,BOOLEAN\/\fR Continue with the simulation like nothing happened after the linear solver did not converge. Default: false diff --git a/opm/simulators/linalg/FlowLinearSolverParameters.cpp b/opm/simulators/linalg/FlowLinearSolverParameters.cpp index c5d9c1d12b1..fe1ededb027 100644 --- a/opm/simulators/linalg/FlowLinearSolverParameters.cpp +++ b/opm/simulators/linalg/FlowLinearSolverParameters.cpp @@ -135,7 +135,7 @@ void FlowLinearSolverParameters::registerParameters() Parameters::Register ("Configuration of solver. Valid options are: cprw (default), " "ilu0, dilu, cpr (an alias for cprw), cpr_quasiimpes, " - "cpr_trueimpes, cpr_trueimpesanalytic, amg or hybrid (experimental). " + "cpr_trueimpes, cpr_trueimpesanalytic, cpr_coatsblackoil, cpr_coatscompositional, amg or hybrid (experimental). " "Alternatively, you can request a configuration to be read from a " "JSON file by giving the filename here, ending with '.json'"); Parameters::Register diff --git a/opm/simulators/linalg/ISTLSolver.hpp b/opm/simulators/linalg/ISTLSolver.hpp index aeafcb202ba..9e2dc18edbf 100644 --- a/opm/simulators/linalg/ISTLSolver.hpp +++ b/opm/simulators/linalg/ISTLSolver.hpp @@ -632,11 +632,49 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, ); return weights; }; + } else if (weightsType == "coatsblackoil") { + // Coats (1980) pressure-equation weights for the blackoil model. + // Extends the analytic formula to include vaporized water (Rvw) + // and dissolved gas in water (Rsw) cross-terms. + weightsCalculator = + [this, pressIndex, enableThreadParallel] + { + Vector weights(rhs_->size()); + ElementContext elemCtx(simulator_); + Amg::getCoatsWeightsBlackoil(pressIndex, + weights, + elemCtx, + simulator_.model(), + *element_chunks_, + enableThreadParallel + ); + return weights; + }; + } else if (weightsType == "coatscompositional") { + // Coats (1980) pressure-equation weights for the compositional + // (ptflash) model. Uses saturation-weighted specific volumes as + // weights so that the weighted sum of component mass balances + // yields the total pore-volume (volumetric) balance. + weightsCalculator = + [this, pressIndex, enableThreadParallel] + { + Vector weights(rhs_->size()); + ElementContext elemCtx(simulator_); + Amg::getCoatsWeightsCompositional(pressIndex, + weights, + elemCtx, + simulator_.model(), + *element_chunks_, + enableThreadParallel + ); + return weights; + }; } else { OPM_THROW(std::invalid_argument, "Weights type " + weightsType + "not implemented for cpr." - " Please use quasiimpes, trueimpes or trueimpesanalytic."); + " Please use quasiimpes, trueimpes, trueimpesanalytic," + " coatsblackoil or coatscompositional."); } } return weightsCalculator; diff --git a/opm/simulators/linalg/getQuasiImpesWeights.hpp b/opm/simulators/linalg/getQuasiImpesWeights.hpp index 786faf1d507..d79e8b13a80 100644 --- a/opm/simulators/linalg/getQuasiImpesWeights.hpp +++ b/opm/simulators/linalg/getQuasiImpesWeights.hpp @@ -157,6 +157,31 @@ namespace Amg } #endif + /// \brief Compute IMPES pressure-equation weights via AD storage Jacobian (derivative-based). + /// + /// For each cell, the storage Jacobian dM/dx (where M_i is the accumulation term + /// for equation i and x_j is primary variable j) is assembled from the automatic + /// differentiation (AD) derivatives of the storage residual. The transposed system + /// + /// (dM/dx)^T w = e_pressure + /// + /// is then solved for the weight vector w. This is the "numerical" or + /// derivative-based IMPES approach: it works for any blackoil extension + /// (dissolved gas in water, vaporised water, thermal, solvents, etc.) + /// without requiring a closed-form analytic formula. + /// + /// This method is more expensive than the analytic alternatives + /// (getTrueImpesWeightsAnalytic, getCoatsWeightsBlackoil) because it + /// requires an element-context update and a small dense solve per cell. + /// + /// References: + /// - Wallis, J.R. (1983). "Incomplete Gaussian Elimination as a + /// Preconditioning for Generalized Conjugate Gradient Acceleration." + /// SPE-12265-MS. https://doi.org/10.2118/12265-MS + /// - Cao, H., Tchelepi, H.A., Wallis, J.R., Yardumian, H. (2005). + /// "Parallel Scalable Unstructured CPR-Type Linear Solver for + /// Reservoir Simulation." SPE-96809-MS. + /// https://doi.org/10.2118/96809-MS template void getTrueImpesWeights(int pressureVarIndex, Vector& weights, const ElementContext& elemCtx, @@ -224,6 +249,40 @@ namespace Amg OPM_END_PARALLEL_TRY_CATCH("getTrueImpesWeights() failed: ", elemCtx.simulator().vanguard().grid().comm()); } + /// \brief Compute analytic IMPES/Coats pressure-equation weights for the blackoil model. + /// + /// Equivalent to getCoatsWeightsBlackoil() for the standard blackoil model + /// (no dissolved gas in water, no vaporised water). The weights are: + /// + /// w_water = B_w + /// w_oil = (B_o - R_s B_g) / (1 - R_s R_v) + /// w_gas = (B_g - R_v B_o) / (1 - R_s R_v) + /// + /// where B_alpha = 1/invB_alpha is the formation-volume factor. + /// These weights satisfy the volumetric-balance linear system + /// + /// | invBw 0 0 | | w_w | | 1 | + /// | 0 invBo Rs*invBo| | w_o | = | 1 | + /// | 0 Rv*invBg invBg | | w_g | | 1 | + /// + /// so that the weighted sum of the component mass balances equals the + /// total reservoir pore-volume balance phi*(S_w + S_o + S_g). + /// + /// Undersaturation is handled by switching off the appropriate ratio: + /// - gas dissolved only (GasMeaning::Rs): set R_v = 0 + /// - oil vaporised only (GasMeaning::Rv): set R_s = 0 + /// + /// For the extended blackoil model with dissolved gas in water (R_sw) + /// and vaporised water in gas (R_vw), use getCoatsWeightsBlackoil() + /// instead, which adds the full cross-terms. + /// + /// References: + /// - Wallis, J.R. (1983). "Incomplete Gaussian Elimination as a + /// Preconditioning for Generalized Conjugate Gradient Acceleration." + /// SPE-12265-MS. https://doi.org/10.2118/12265-MS + /// - Coats, K.H. (1980). "An Equation of State Compositional Model." + /// SPE Journal, 20(5), 363-376. SPE-8284-PA. + /// https://doi.org/10.2118/8284-PA template void getTrueImpesWeightsAnalytic(int /*pressureVarIndex*/, Vector& weights, @@ -234,10 +293,10 @@ namespace Amg { // The sequential residual is a linear combination of the // mass balance residuals, with coefficients equal to (for - // water, oil, gas): - // 1/bw, - // (1/bo - rs/bg)/(1-rs*rv) - // (1/bg - rv/bo)/(1-rs*rv) + // water, oil, gas) -- here bw means B_w (formation-volume factor): + // B_w = 1/invBw + // (B_o - Rs*B_g)/(1-Rs*Rv) = (1/invBo - Rs/invBg)/(1-Rs*Rv) + // (B_g - Rv*B_o)/(1-Rs*Rv) = (1/invBg - Rv/invBo)/(1-Rs*Rv) // These coefficients must be applied for both the residual and // Jacobian. using FluidSystem = typename Model::FluidSystem; @@ -282,10 +341,10 @@ namespace Amg double rv = Toolbox::template decay(fs.Rv()); const auto& priVars = solution[index]; if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) { - rs = 0.0; + rs = 0.0;// then oil saturation is zero an rs undetermined, in simulator it is saturated } if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) { - rv = 0.0; + rv = 0.0;// then gas saturation is zero but rv is undetermined, in simulator it is saturated } if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { @@ -306,12 +365,328 @@ namespace Amg (1 / fs.invB(FluidSystem::gasPhaseIdx) - rv / fs.invB(FluidSystem::oilPhaseIdx)) / denominator); } + // in the understaturated case it is possible to modify the weights between saturated and undersaturated equation. + // a natural since bouth only depend on one saturation. + weights[index] = bweights; } } OPM_END_PARALLEL_TRY_CATCH("getTrueImpesAnalyticWeights() failed: ", elemCtx.simulator().vanguard().grid().comm()); } + + /// \brief Compute Coats pressure-equation weights for the blackoil model, + /// including the extended Rsw/Rvw cross terms. + /// + /// Based on Coats (1980), the pressure equation is formed by taking a linear + /// combination of the component mass-balance equations with coefficients that + /// convert them into a total reservoir-volume balance: + /// + /// sum_i w_i * (component i mass balance) = d(total pore volume)/dt + flux terms + /// + /// The weights ensure that the weighted combination of storage terms satisfies + /// + /// w_w * invBw * S_w + w_o * invBo * S_o + w_g * invBg * S_g = S_w + S_o + S_g + /// + /// For a 3-phase blackoil system the weights are derived by solving: + /// + /// | invBw 0 Rsw*invBw | | ww | | 1 | + /// | 0 invBo Rs*invBo | | wo | = | 1 | + /// | Rvw*invBg Rv*invBg invBg | | wg | | 1 | + /// + /// giving (with D = 1 - Rs*Rv - Rvw*Rsw): + /// ww = [ Bw*(1 - Rs*Rv) - Rsw*(Bg - Rv*Bo) ] / D + /// wo = [ Bo*(1 - Rvw*Rsw) - Rs*(Bg - Rvw*Bw) ] / D + /// wg = ( Bg - Rvw*Bw - Rv*Bo ) / D + /// + /// where B_alpha = 1/invB_alpha. The function handles all active-phase + /// combinations and the optional dissolved-gas-in-water (Rsw) and + /// vaporized-water-in-gas (Rvw) extensions. + /// + /// Without Rsw and Rvw this formula reduces to getTrueImpesWeightsAnalytic(). + /// + /// Undersaturation is handled by switching off the appropriate ratio: + /// - oil undersaturated (GasMeaning::Rv encoding): set R_s = 0 + /// - gas undersaturated (GasMeaning::Rs encoding): set R_v = 0 + /// + /// References: + /// - Coats, K.H. (1980). "An Equation of State Compositional Model." + /// SPE Journal, 20(5), 363-376. SPE-8284-PA. + /// https://doi.org/10.2118/8284-PA + /// - Coats, K.H. (1999). "A Note on IMPES and Some IMPES-Based + /// Simulation Models." SPE Journal, 5(3), 245-251. SPE-65998-PA. + /// https://doi.org/10.2118/65998-PA + template + void getCoatsWeightsBlackoil(int /*pressureVarIndex*/, + Vector& weights, + const ElementContext& elemCtx, + const Model& model, + const ElementChunksType& element_chunks, + [[maybe_unused]] bool enable_thread_parallel) + { + using FluidSystem = typename Model::FluidSystem; + using LhsEval = double; + using PrimaryVariables = typename Model::PrimaryVariables; + using VectorBlockType = typename Vector::block_type; + using Evaluation = + typename std::decay_t::block_type; + using Toolbox = MathToolbox; + + const auto& solution = model.solution(/*timeIdx=*/0); + VectorBlockType bweights; + + OPM_BEGIN_PARALLEL_TRY_CATCH(); +#ifdef _OPENMP +#pragma omp parallel for private(bweights) if(enable_thread_parallel) +#endif + for (const auto& chunk : element_chunks) { + ElementContext localElemCtx(elemCtx.simulator()); + + for (const auto& elem : chunk) { + localElemCtx.updatePrimaryStencil(elem); + localElemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + + const auto index = localElemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& intQuants = localElemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + const auto& priVars = solution[index]; + + bweights = 0.0; + + const bool waterActive = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); + const bool oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); + const bool gasActive = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + + // Formation volume factors B_alpha = 1 / invB_alpha + const double Bw = waterActive + ? Toolbox::template decay(1.0 / fs.invB(FluidSystem::waterPhaseIdx)) + : 0.0; + const double Bo = oilActive + ? Toolbox::template decay(1.0 / fs.invB(FluidSystem::oilPhaseIdx)) + : 0.0; + const double Bg = gasActive + ? Toolbox::template decay(1.0 / fs.invB(FluidSystem::gasPhaseIdx)) + : 0.0; + + // Solution-gas/vaporized-oil ratios. Zero out if the primary + // variable encodes the saturated value instead (switching variables). + double rs = 0.0, rv = 0.0, rsw = 0.0, rvw = 0.0; + + if (oilActive && gasActive) { + if (FluidSystem::enableDissolvedGas()) { + rs = Toolbox::template decay(fs.Rs()); + // When the gas primary variable represents Rv (oil vaporized in gas), + // the oil phase is undersaturated, so Rs should be treated as zero. + if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) { + rs = 0.0; + } + } + if (FluidSystem::enableVaporizedOil()) { + rv = Toolbox::template decay(fs.Rv()); + // When the gas primary variable represents Rs (gas dissolved in oil), + // the gas phase is undersaturated, so Rv should be treated as zero. + if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) { + rv = 0.0; + } + } + } + if (waterActive && gasActive) { + if (FluidSystem::enableDissolvedGasInWater()) { + rsw = Toolbox::template decay(fs.Rsw()); + } + if (FluidSystem::enableVaporizedWater()) { + rvw = Toolbox::template decay(fs.Rvw()); + } + } + + // Denominator D = 1 - Rs*Rv - Rvw*Rsw + const double D = 1.0 - rs * rv - rvw * rsw; + + // Compute weights by solving the linear system derived from the + // requirement that sum_i w_i * M_i = total reservoir saturation. + // + // 3-phase (w, o, g) with full cross-terms: + // ww = [ Bw*(1 - Rs*Rv) - Rsw*(Bg - Rv*Bo) ] / D + // wo = [ Bo*(1 - Rvw*Rsw) - Rs*(Bg - Rvw*Bw) ] / D + // wg = ( Bg - Rvw*Bw - Rv*Bo ) / D + // + // Degenerate (fewer active phases) cases collapse correctly: + // water-only: ww = Bw + // gas+oil: wo = (Bo - Rs*Bg)/(1-Rs*Rv), wg = (Bg - Rv*Bo)/(1-Rs*Rv) + // water+gas: ww = (Bw - Rsw*Bg)/(1-Rsw*Rvw), wg = (Bg - Rvw*Bw)/(1-Rsw*Rvw) + // water+oil: ww = Bw, wo = Bo (no cross-dissolution) + + if (waterActive) { + const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx( + FluidSystem::solventComponentIndex(FluidSystem::waterPhaseIdx)); + bweights[activeCompIdx] = static_cast( + (Bw * (1.0 - rs * rv) - rsw * (Bg - rv * Bo)) / D); + } + if (oilActive) { + const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx( + FluidSystem::solventComponentIndex(FluidSystem::oilPhaseIdx)); + bweights[activeCompIdx] = static_cast( + (Bo * (1.0 - rvw * rsw) - rs * (Bg - rvw * Bw)) / D); + } + if (gasActive) { + const unsigned activeCompIdx = FluidSystem::canonicalToActiveCompIdx( + FluidSystem::solventComponentIndex(FluidSystem::gasPhaseIdx)); + bweights[activeCompIdx] = static_cast( + (Bg - rvw * Bw - rv * Bo) / D); + } + + // Normalize so that the maximum absolute weight equals one. + const double abs_max = + *std::ranges::max_element(bweights, + [](double a, double b) + { return std::fabs(a) < std::fabs(b); }); + if (std::fabs(abs_max) > 0.0) { + bweights /= std::fabs(abs_max); + } + + weights[index] = bweights; + } + } + OPM_END_PARALLEL_TRY_CATCH("getCoatsWeightsBlackoil() failed: ", + elemCtx.simulator().vanguard().grid().comm()); + } + + /// \brief Compute Coats pressure-equation weights for the compositional (ptflash) model. + /// + /// Based on Coats (1980), the pressure equation is obtained by forming a + /// volumetric balance from the component mass balances. Unlike the blackoil + /// formulation (getCoatsWeightsBlackoil), the compositional model uses + /// saturation-weighted specific volumes derived from the phase compositions + /// and densities provided by the ptflash equilibrium calculation. + /// + /// References: + /// - Coats, K.H. (1980). "An Equation of State Compositional Model." + /// SPE Journal, 20(5), 363-376. SPE-8284-PA. + /// https://doi.org/10.2118/8284-PA + /// + /// For each hydrocarbon + /// component i with mass-based storage term + /// + /// M_i = phi * sum_{alpha in HC} omega_{alpha,i} * rho_alpha * S_alpha + /// + /// the Coats weight is chosen so that + /// + /// sum_i w_i * M_i = phi * (S_o + S_g) (total HC pore volume) + /// + /// which is satisfied by taking + /// + /// w_i = (sum_{alpha in HC} omega_{alpha,i} * S_alpha) + /// / (sum_{alpha in HC} omega_{alpha,i} * rho_alpha * S_alpha) + /// = saturation-weighted mass fraction / saturation-weighted mass density + /// + /// For the water equation (if present): w_water = 1 / rho_w + /// + /// The energy equation (if present) is left with zero weight because + /// temperature preconditioning is handled separately. + template + void getCoatsWeightsCompositional(int /*pressureVarIndex*/, + Vector& weights, + const ElementContext& elemCtx, + const Model& model, + const ElementChunksType& element_chunks, + [[maybe_unused]] bool enable_thread_parallel) + { + using FluidSystem = typename Model::FluidSystem; + using LhsEval = double; + using VectorBlockType = typename Vector::block_type; + using Evaluation = + typename std::decay_t::block_type; + using Toolbox = MathToolbox; + using Indices = typename Model::Indices; + + constexpr int numComponents = FluidSystem::numComponents; + constexpr int numPhases = FluidSystem::numPhases; + constexpr bool waterEnabled = Indices::waterEnabled; + + VectorBlockType bweights; + + OPM_BEGIN_PARALLEL_TRY_CATCH(); +#ifdef _OPENMP +#pragma omp parallel for private(bweights) if(enable_thread_parallel) +#endif + for (const auto& chunk : element_chunks) { + ElementContext localElemCtx(elemCtx.simulator()); + + for (const auto& elem : chunk) { + localElemCtx.updatePrimaryStencil(elem); + localElemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + + const auto index = localElemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& intQuants = localElemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + + bweights = 0.0; + + // Iterate over HC components (equations conti0EqIdx + 0 .. numComponents-1). + // For each component i compute: + // num_i = sum_{alpha in HC} massFraction(alpha,i) * saturation(alpha) + // den_i = sum_{alpha in HC} massFraction(alpha,i) * density(alpha) * saturation(alpha) + // w_i = num_i / den_i (= saturation-weighted specific volume) + // + // When den_i is zero (component absent from all HC phases), set w_i = 0. + + for (int compIdx = 0; compIdx < numComponents; ++compIdx) { + double num_i = 0.0; + double den_i = 0.0; + + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + // Skip the immiscible water phase for HC component weights + if (waterEnabled && + phaseIdx == static_cast(FluidSystem::waterPhaseIdx)) + { + continue; + } + const double S_alpha = + Toolbox::template decay(fs.saturation(phaseIdx)); + const double omega_ai = + Toolbox::template decay(fs.massFraction(phaseIdx, compIdx)); + const double rho_alpha = + Toolbox::template decay(fs.density(phaseIdx)); + + num_i += omega_ai * S_alpha; + den_i += omega_ai * rho_alpha * S_alpha; + } + + const int eqIdx = Indices::conti0EqIdx + compIdx; + bweights[eqIdx] = (den_i > 0.0) ? static_cast(num_i / den_i) + : static_cast(0.0); + } + + // Water equation weight: w_water = 1 / rho_w (ensures phi*S_w contribution) + if constexpr (waterEnabled) { + const int waterEqIdx = Indices::conti0EqIdx + numComponents; + const double rho_w = + Toolbox::template decay( + fs.density(FluidSystem::waterPhaseIdx)); + bweights[waterEqIdx] = (rho_w > 0.0) + ? static_cast(1.0 / rho_w) + : static_cast(0.0); + } + + // Energy equation (if present, numEq > numComponents + waterEnabled): + // left at zero -- temperature preconditioning is a separate task. + + // Normalize so that the maximum absolute weight equals one. + const double abs_max = + *std::ranges::max_element(bweights, + [](double a, double b) + { return std::fabs(a) < std::fabs(b); }); + if (std::fabs(abs_max) > 0.0) { + bweights /= std::fabs(abs_max); + } + + weights[index] = bweights; + } + } + OPM_END_PARALLEL_TRY_CATCH("getCoatsWeightsCompositional() failed: ", + elemCtx.simulator().vanguard().grid().comm()); + } + } // namespace Amg } // namespace Opm diff --git a/opm/simulators/linalg/gpuistl/ISTLSolverGPUISTL.hpp b/opm/simulators/linalg/gpuistl/ISTLSolverGPUISTL.hpp index d4f4d012cbf..21cb177c3ef 100644 --- a/opm/simulators/linalg/gpuistl/ISTLSolverGPUISTL.hpp +++ b/opm/simulators/linalg/gpuistl/ISTLSolverGPUISTL.hpp @@ -408,6 +408,46 @@ class ISTLSolverGPUISTL : public AbstractISTLSolvercopyFromHostAsync(m_cpuWeights); + return *m_weights; + }; + } else if (weightsType == "coatsblackoil") { + // Create CPU vector for the weights and initialize GPU vector + m_cpuWeights.resize(m_matrix->N()); + m_pinnedWeightsMemory = std::make_unique>( + const_cast(&m_cpuWeights[0][0]), m_cpuWeights.dim()); + m_weights.emplace(m_cpuWeights); + const bool enableThreadParallel = m_parameters.cpr_weights_thread_parallel_; + // CPU implementation wrapped for GPU + weightsCalculator = [this, enableThreadParallel]() -> GPUVector& { + ElementContext elemCtx(m_simulator); + Amg::getCoatsWeightsBlackoil(pressureIndex, + m_cpuWeights, + elemCtx, m_simulator.model(), + m_element_chunks, + enableThreadParallel); + + // Copy CPU vector to GPU vector using main stream and asynchronous transfer + m_weights->copyFromHostAsync(m_cpuWeights); + return *m_weights; + }; + } else if (weightsType == "coatscompositional") { + // Create CPU vector for the weights and initialize GPU vector + m_cpuWeights.resize(m_matrix->N()); + m_pinnedWeightsMemory = std::make_unique>( + const_cast(&m_cpuWeights[0][0]), m_cpuWeights.dim()); + m_weights.emplace(m_cpuWeights); + const bool enableThreadParallel = m_parameters.cpr_weights_thread_parallel_; + // CPU implementation wrapped for GPU + weightsCalculator = [this, enableThreadParallel]() -> GPUVector& { + ElementContext elemCtx(m_simulator); + Amg::getCoatsWeightsCompositional(pressureIndex, + m_cpuWeights, + elemCtx, m_simulator.model(), + m_element_chunks, + enableThreadParallel); + // Copy CPU vector to GPU vector using main stream and asynchronous transfer m_weights->copyFromHostAsync(m_cpuWeights); return *m_weights; @@ -416,7 +456,8 @@ class ISTLSolverGPUISTL : public AbstractISTLSolver. +*/ + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include +#include +#include + +#define BOOST_TEST_MODULE ImpesWeightsTest +#include + +namespace { + +constexpr double tolerance = 1.0e-12; + +struct Weights { + double water = 0.0; + double oil = 0.0; + double gas = 0.0; +}; + +Weights analyticBlackoilWeights(const double Bw, + const double Bo, + const double Bg, + const double rs, + const double rv, + const bool waterActive = true, + const bool oilActive = true, + const bool gasActive = true) +{ + Weights weights; + const double denominator = (oilActive && gasActive) ? (1.0 - rs * rv) : 1.0; + + if (waterActive) { + weights.water = Bw; + } + if (oilActive) { + weights.oil = (Bo - rs * Bg) / denominator; + } + if (gasActive) { + weights.gas = (Bg - rv * Bo) / denominator; + } + + return weights; +} + +Weights coatsBlackoilWeights(const double Bw, + const double Bo, + const double Bg, + const double rs, + const double rv, + const double rsw, + const double rvw, + const bool waterActive = true, + const bool oilActive = true, + const bool gasActive = true) +{ + Weights weights; + const double denominator = 1.0 - rs * rv - rvw * rsw; + + if (waterActive) { + weights.water = (Bw * (1.0 - rs * rv) - rsw * (Bg - rv * Bo)) / denominator; + } + if (oilActive) { + weights.oil = (Bo * (1.0 - rvw * rsw) - rs * (Bg - rvw * Bw)) / denominator; + } + if (gasActive) { + weights.gas = (Bg - rvw * Bw - rv * Bo) / denominator; + } + + return weights; +} + +double absMax(const Weights& weights) +{ + return std::max({std::fabs(weights.water), std::fabs(weights.oil), std::fabs(weights.gas)}); +} + +Weights normalizeByAbsMax(Weights weights) +{ + const double scale = absMax(weights); + if (scale > 0.0) { + weights.water /= scale; + weights.oil /= scale; + weights.gas /= scale; + } + + return weights; +} + +void checkStandardBlackoilBalance(const Weights& weights, + const double Bw, + const double Bo, + const double Bg, + const double rs, + const double rv, + const bool waterActive = true, + const bool oilActive = true, + const bool gasActive = true) +{ + if (waterActive) { + BOOST_CHECK_CLOSE_FRACTION(weights.water / Bw, 1.0, tolerance); + } + if (oilActive) { + BOOST_CHECK_CLOSE_FRACTION((weights.oil + rs * weights.gas) / Bo, 1.0, tolerance); + } + if (gasActive) { + BOOST_CHECK_CLOSE_FRACTION((rv * weights.oil + weights.gas) / Bg, 1.0, tolerance); + } +} + +void checkExtendedBlackoilBalance(const Weights& weights, + const double Bw, + const double Bo, + const double Bg, + const double rs, + const double rv, + const double rsw, + const double rvw, + const bool waterActive = true, + const bool oilActive = true, + const bool gasActive = true) +{ + if (waterActive) { + BOOST_CHECK_CLOSE_FRACTION((weights.water + rsw * weights.gas) / Bw, 1.0, tolerance); + } + if (oilActive) { + BOOST_CHECK_CLOSE_FRACTION((weights.oil + rs * weights.gas) / Bo, 1.0, tolerance); + } + if (gasActive) { + BOOST_CHECK_CLOSE_FRACTION((rvw * weights.water + rv * weights.oil + weights.gas) / Bg, + 1.0, + tolerance); + } +} + +} // anonymous namespace + +BOOST_AUTO_TEST_CASE(StandardBlackoilWeightsSatisfyVolumeBalance) +{ + const double Bw = 1.05; + const double Bo = 1.20; + const double Bg = 0.85; + const double rs = 0.25; + const double rv = 0.15; + + const auto weights = analyticBlackoilWeights(Bw, Bo, Bg, rs, rv); + + checkStandardBlackoilBalance(weights, Bw, Bo, Bg, rs, rv); +} + +BOOST_AUTO_TEST_CASE(UndersaturatedLimitsReduceToExpectedFormulas) +{ + const double Bw = 1.02; + const double Bo = 1.18; + const double Bg = 0.81; + + { + const double rs = 0.0; + const double rv = 0.12; + const auto weights = analyticBlackoilWeights(Bw, Bo, Bg, rs, rv); + + BOOST_CHECK_CLOSE_FRACTION(weights.water, Bw, tolerance); + BOOST_CHECK_CLOSE_FRACTION(weights.oil, Bo, tolerance); + BOOST_CHECK_CLOSE_FRACTION(weights.gas, Bg - rv * Bo, tolerance); + checkStandardBlackoilBalance(weights, Bw, Bo, Bg, rs, rv); + } + + { + const double rs = 0.18; + const double rv = 0.0; + const auto weights = analyticBlackoilWeights(Bw, Bo, Bg, rs, rv); + + BOOST_CHECK_CLOSE_FRACTION(weights.water, Bw, tolerance); + BOOST_CHECK_CLOSE_FRACTION(weights.oil, Bo - rs * Bg, tolerance); + BOOST_CHECK_CLOSE_FRACTION(weights.gas, Bg, tolerance); + checkStandardBlackoilBalance(weights, Bw, Bo, Bg, rs, rv); + } +} + +BOOST_AUTO_TEST_CASE(AnalyticAndCoatsMatchWithoutWaterGasCrossTerms) +{ + const double Bw = 1.08; + const double Bo = 1.24; + const double Bg = 0.79; + const double rs = 0.22; + const double rv = 0.11; + + const auto analytic = analyticBlackoilWeights(Bw, Bo, Bg, rs, rv); + const auto coats = coatsBlackoilWeights(Bw, Bo, Bg, rs, rv, 0.0, 0.0); + const auto normalizedCoats = normalizeByAbsMax(coats); + const auto normalizedAnalytic = normalizeByAbsMax(analytic); + + BOOST_CHECK_CLOSE_FRACTION(coats.water, analytic.water, tolerance); + BOOST_CHECK_CLOSE_FRACTION(coats.oil, analytic.oil, tolerance); + BOOST_CHECK_CLOSE_FRACTION(coats.gas, analytic.gas, tolerance); + + BOOST_CHECK_CLOSE_FRACTION(normalizedCoats.water, normalizedAnalytic.water, tolerance); + BOOST_CHECK_CLOSE_FRACTION(normalizedCoats.oil, normalizedAnalytic.oil, tolerance); + BOOST_CHECK_CLOSE_FRACTION(normalizedCoats.gas, normalizedAnalytic.gas, tolerance); +} + +BOOST_AUTO_TEST_CASE(ExtendedCoatsWeightsSatisfyVolumeBalance) +{ + const double Bw = 1.04; + const double Bo = 1.19; + const double Bg = 0.83; + const double rs = 0.21; + const double rv = 0.09; + const double rsw = 0.04; + const double rvw = 0.03; + + const auto weights = coatsBlackoilWeights(Bw, Bo, Bg, rs, rv, rsw, rvw); + + checkExtendedBlackoilBalance(weights, Bw, Bo, Bg, rs, rv, rsw, rvw); +} + +BOOST_AUTO_TEST_CASE(SingleAndTwoPhaseLimitsRemainConsistent) +{ + { + const double Bw = 1.07; + const auto weights = coatsBlackoilWeights(Bw, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + true, false, false); + + BOOST_CHECK_CLOSE_FRACTION(weights.water, Bw, tolerance); + BOOST_CHECK_SMALL(weights.oil, tolerance); + BOOST_CHECK_SMALL(weights.gas, tolerance); + checkExtendedBlackoilBalance(weights, Bw, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, + true, false, false); + } + + { + const double Bo = 1.23; + const double Bg = 0.82; + const double rs = 0.19; + const double rv = 0.08; + const auto weights = analyticBlackoilWeights(0.0, Bo, Bg, rs, rv, + false, true, true); + + BOOST_CHECK_SMALL(weights.water, tolerance); + checkStandardBlackoilBalance(weights, 1.0, Bo, Bg, rs, rv, + false, true, true); + } +} \ No newline at end of file