diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index b07e2d5..36c7dd9 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -12,7 +12,7 @@ jobs:
strategy:
fail-fast: false
matrix:
- python-version: ['3.8']
+ python-version: ['3.10']
os: [ubuntu-latest]
runs-on: ${{ matrix.os }}
diff --git a/examples/cemracs2023/ml/recovery.py b/examples/cemracs2023/ml/recovery.py
index c980d36..e5a2772 100644
--- a/examples/cemracs2023/ml/recovery.py
+++ b/examples/cemracs2023/ml/recovery.py
@@ -6,10 +6,11 @@
"""
import os
+
+import matplotlib.pyplot as plt
import numpy as np
from ecl.summary import EclSum
from mako.template import Template
-import matplotlib.pyplot as plt
npoints, npruns = 20, 5
tmin, tmax = 0, 30
@@ -18,7 +19,7 @@
FLOW = "/Users/dmar/Github/opm/build/opm-simulators/bin/flow"
FLAGS = (
" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0"
- + " --ecl-enable-drift-compensation=0 --newton-max-iterations=50"
+ + " --enable-drift-compensation=0 --newton-max-iterations=50"
+ " --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5"
+ " --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true"
+ " --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false"
diff --git a/examples/cemracs2023/ml_cyclic_h2/h2.mako b/examples/cemracs2023/ml_cyclic_h2/h2.mako
index f6d0868..012baa6 100644
--- a/examples/cemracs2023/ml_cyclic_h2/h2.mako
+++ b/examples/cemracs2023/ml_cyclic_h2/h2.mako
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-${flow} --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+${flow} --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
h2store okoroafor2023 #Model (co2store/h2store)
diff --git a/examples/cemracs2023/ml_example_co2/co2.mako b/examples/cemracs2023/ml_example_co2/co2.mako
index 7c2fc1d..8b7454e 100644
--- a/examples/cemracs2023/ml_example_co2/co2.mako
+++ b/examples/cemracs2023/ml_example_co2/co2.mako
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-${flow} --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+${flow} --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
co2store base #Model (co2store/h2store)
diff --git a/examples/cemracs2023/ml_example_h2/h2.mako b/examples/cemracs2023/ml_example_h2/h2.mako
index 84d76e6..de832be 100644
--- a/examples/cemracs2023/ml_example_h2/h2.mako
+++ b/examples/cemracs2023/ml_example_h2/h2.mako
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-${flow} --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+${flow} --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
h2store okoroafor2023 #Model (co2store/h2store)
diff --git a/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_finescale_wellmodel.mako b/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_finescale_wellmodel.mako
index cbe0f81..97eff29 100644
--- a/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_finescale_wellmodel.mako
+++ b/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_finescale_wellmodel.mako
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-${flow} --ml-wi-filename="" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+${flow} --ml-wi-filename="" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
co2store base #Model (co2store/h2store)
diff --git a/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_finescale_wellmodel_large_reservoir.mako b/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_finescale_wellmodel_large_reservoir.mako
index d7f80aa..51555d4 100644
--- a/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_finescale_wellmodel_large_reservoir.mako
+++ b/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_finescale_wellmodel_large_reservoir.mako
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-${flow} --ml-wi-filename="" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+${flow} --ml-wi-filename="" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
co2store base #Model (co2store/h2store)
diff --git a/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_flow_wellmodel.mako b/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_flow_wellmodel.mako
index 4169ad3..59db0ba 100644
--- a/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_flow_wellmodel.mako
+++ b/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_flow_wellmodel.mako
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-${flow} --ml-wi-filename="" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+${flow} --ml-wi-filename="" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
co2store base #Model (co2store/h2store)
diff --git a/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_flow_wellmodel_large_reservoir.mako b/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_flow_wellmodel_large_reservoir.mako
index 13c456d..f6680f6 100644
--- a/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_flow_wellmodel_large_reservoir.mako
+++ b/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_flow_wellmodel_large_reservoir.mako
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-${flow} --ml-wi-filename="" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+${flow} --ml-wi-filename="" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
co2store base #Model (co2store/h2store)
diff --git a/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_ml_wellmodel.mako b/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_ml_wellmodel.mako
index f6891df..738aac7 100644
--- a/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_ml_wellmodel.mako
+++ b/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_ml_wellmodel.mako
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-${flow} --ml-wi-filename="${pwd}/model_pressure_radius_WI/WI.model" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+${flow} --ml-wi-filename="${pwd}/model_pressure_radius_WI/WI.model" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
co2store base #Model (co2store/h2store)
diff --git a/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_ml_wellmodel_large_reservoir.mako b/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_ml_wellmodel_large_reservoir.mako
index 21b440b..6afcbf1 100644
--- a/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_ml_wellmodel_large_reservoir.mako
+++ b/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/co2_3d_ml_wellmodel_large_reservoir.mako
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-${flow} --ml-wi-filename="${pwd}/model_pressure_radius_WI/WI.model" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+${flow} --ml-wi-filename="${pwd}/model_pressure_radius_WI/WI.model" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
co2store base #Model (co2store/h2store)
diff --git a/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/run_ensemble.py b/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/run_ensemble.py
index 8601436..8ef6be6 100644
--- a/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/run_ensemble.py
+++ b/examples/cemracs2023/peaceman_well-model_nn/2023-08-10_ml_well_model_CO2_realistic_critical_point/run_ensemble.py
@@ -10,6 +10,7 @@
1. WI [m*s]
"""
+
from __future__ import annotations
import logging
@@ -54,13 +55,15 @@
X: float = 2.500000e-01 # Outer coordinates of first cell.
Y: float = -1.443376e-01
WELL_RADIUS: float = math.sqrt(X**2 + Y**2) # unit: [m]; Fixed during training.
-DENSITY: float = 12.9788 # unit: kg/m^3; for 72 bar, 30.9780 °C. Is this at surface conditions or not?
+DENSITY: float = (
+ 12.9788 # unit: kg/m^3; for 72 bar, 30.9780 °C. Is this at surface conditions or not?
+)
VISCOSITY: float = 1.52786e-05 # unit: Pa*s; for 72 bar, 30.9780 °C
FLOW = "flow"
FLAGS = (
" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0"
- + " --ecl-enable-drift-compensation=0 --newton-max-iterations=50"
+ + " --enable-drift-compensation=0 --newton-max-iterations=50"
+ " --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5"
+ " --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true"
+ " --enable-opm-rst-file=true --linear-solver=cprw"
diff --git a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/co2_3d_flow_wellmodel.mako b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/co2_3d_flow_wellmodel.mako
index 77a7d2a..a046624 100644
--- a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/co2_3d_flow_wellmodel.mako
+++ b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/co2_3d_flow_wellmodel.mako
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-${flow} --ml-wi-filename="" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+${flow} --ml-wi-filename="" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
co2store base #Model (co2store/h2store)
diff --git a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/co2_3d_ml_wellmodel.mako b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/co2_3d_ml_wellmodel.mako
index 281cabc..6848622 100644
--- a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/co2_3d_ml_wellmodel.mako
+++ b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/co2_3d_ml_wellmodel.mako
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-${flow} --ml-wi-filename="${pwd}/water_re.modelPeaceman" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+${flow} --ml-wi-filename="${pwd}/water_re.modelPeaceman" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
co2store base #Model (co2store/h2store)
diff --git a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/co2_nearwell.mako b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/co2_nearwell.mako
index 890b2c9..7cb5938 100644
--- a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/co2_nearwell.mako
+++ b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/co2_nearwell.mako
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-${flow} --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+${flow} --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
co2store base #Model (co2store/h2store)
diff --git a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/example.py b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/example.py
index 2d2bab5..df740c5 100644
--- a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/example.py
+++ b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_co2/example.py
@@ -24,12 +24,15 @@ def compute_peaceman(k_h: float, r_e: float, r_w: float) -> float:
from the Peaceman well model.
.. math::
WI\cdot\frac{\mu}{\rho} = \frac{2\pi hk}{\ln (r_e/r_w)}
- Parameters:
- k_h: Permeability times the cell thickness (thickness fix to 1 m).
- r_e: Equivalent well-block radius.
- r_w: Wellbore radius.
+
+ Args:
+ k_h (float): Permeability times the cell thickness (thickness fix to 1 m).
+ r_e (float): Equivalent well-block radius.
+ r_w (float): Wellbore radius.
+
Returns:
- :math:`WI\cdot\frac{\mu}{\rho}`
+ float: :math:`WI\cdot\frac{\mu}{\rho}`
+
"""
w_i = (2 * math.pi * k_h) / (math.log(r_e / r_w))
return w_i
diff --git a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/example.py b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/example.py
index fb09a67..7f6453c 100644
--- a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/example.py
+++ b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/example.py
@@ -7,13 +7,14 @@
branches. This can be achieve by running the build_dune_and_opm-flow_ml-peaceman-branch.bash
"""
-import os
import math
+import os
+
+import matplotlib
+import matplotlib.pyplot as plt
import numpy as np
from ecl.eclfile import EclFile
from mako.template import Template
-import matplotlib
-import matplotlib.pyplot as plt
np.random.seed(7)
@@ -23,12 +24,14 @@ def compute_peaceman(k_h: float, r_e: float, r_w: float) -> float:
from the Peaceman well model.
.. math::
WI\cdot\frac{\mu}{\rho} = \frac{2\pi hk}{\ln (r_e/r_w)}
- Parameters:
- k_h: Permeability times the cell thickness (thickness fix to 1 m).
- r_e: Equivalent well-block radius.
- r_w: Wellbore radius.
+
+ Args:
+ k_h (float): Permeability times the cell thickness (thickness fix to 1 m).
+ r_e (float): Equivalent well-block radius.
+ r_w (float): Wellbore radius.
+
Returns:
- :math:`WI\cdot\frac{\mu}{\rho}`
+ float: :math:`WI\cdot\frac{\mu}{\rho}`
"""
w_i = (2 * math.pi * k_h) / (math.log(r_e / r_w))
return w_i
@@ -37,19 +40,21 @@ def compute_peaceman(k_h: float, r_e: float, r_w: float) -> float:
# Give the path to the OPM Flow ML repos
PYOPM = "/Users/dmar/Github/pyopmnearwell/examples/cemracs2023/peaceman_well-model_nn"
FLOW = f"{PYOPM}/build/opm-simulators/bin/flow_gaswater_dissolution_diffuse"
-PERMEABILITY = 1e-12 # K between 1e-13 to 1e-12 m2
-WELLRADI = 0.1 # Between 0.025 to 0.125 m
-RATE = 1.668e-05 # [sm3/s] Fix to 1.665e-2 kg/s, ref_dens = 998.108 kg/sm3
+PERMEABILITY = 1e-12 # K between 1e-13 to 1e-12 m2
+WELLRADI = 0.1 # Between 0.025 to 0.125 m
+RATE = 1.668e-05 # [sm3/s] Fix to 1.665e-2 kg/s, ref_dens = 998.108 kg/sm3
BOFAC = 1.0
VISCOSCITY = 0.6532 * 0.001
-CLIP = 2 # Remove the distances with value less than this [m]
+CLIP = 2 # Remove the distances with value less than this [m]
-#Run the main routine
+# Run the main routine
NPOINTS, NPRUNS = 20, 5
-PERMEABILITYHS = np.linspace(1e-13, 1e-11, NPOINTS) # K between 1e-13 to 1e-12 m2 and H betweem 1 to 10, i,e, [1e-13,1e-11]
-#WELLRADIS = [0.025, .05, .1, .125]
-#WELLRADIS = np.linspace(.05, .125, NPOINTS)
-WELLRADIS = [0.1, .125]
+PERMEABILITYHS = np.linspace(
+ 1e-13, 1e-11, NPOINTS
+) # K between 1e-13 to 1e-12 m2 and H betweem 1 to 10, i,e, [1e-13,1e-11]
+# WELLRADIS = [0.025, .05, .1, .125]
+# WELLRADIS = np.linspace(.05, .125, NPOINTS)
+WELLRADIS = [0.1, 0.125]
nradis = len(WELLRADIS)
nperm = len(PERMEABILITYHS)
wi_simulated = []
@@ -61,7 +66,13 @@ def compute_peaceman(k_h: float, r_e: float, r_w: float) -> float:
for k, WELLRADI in enumerate(WELLRADIS):
r_w.append(WELLRADI)
for i, PERMEABILITYH in enumerate(PERMEABILITYHS):
- var = {"flow": FLOW, "perm": PERMEABILITYH, "radius": WELLRADI, "rate": RATE, "pwd": os.getcwd()}
+ var = {
+ "flow": FLOW,
+ "perm": PERMEABILITYH,
+ "radius": WELLRADI,
+ "rate": RATE,
+ "pwd": os.getcwd(),
+ }
filledtemplate = mytemplate.render(**var)
with open(
f"h2o_{i+k*nperm}.txt",
@@ -89,17 +100,21 @@ def compute_peaceman(k_h: float, r_e: float, r_w: float) -> float:
wi_analytical.append([])
for r in r_e[-1]:
wi_analytical[-1].append(
- compute_peaceman(PERMEABILITYH, r,r_wi[NPRUNS*i+j]) * BOFAC / VISCOSCITY
+ compute_peaceman(PERMEABILITYH, r, r_wi[NPRUNS * i + j])
+ * BOFAC
+ / VISCOSCITY
)
rst = EclFile(f"./h2o_{NPRUNS*i+j}/output/H2O_{NPRUNS*i+j}.UNRST")
pressure = np.array(rst.iget_kw("PRESSURE")[-1])
pw = pressure[0]
- cell_pressures = pressure[len(pressure)-len(r_e[-1]):]
- wi_simulated.append(RATE / ((pw - cell_pressures) * 1e5)) # 1e5 to connvert from bar to Pascals
+ cell_pressures = pressure[len(pressure) - len(r_e[-1]) :]
+ wi_simulated.append(
+ RATE / ((pw - cell_pressures) * 1e5)
+ ) # 1e5 to connvert from bar to Pascals
axis.plot(
r_e[-1],
wi_simulated[-1],
- color=matplotlib.colormaps["tab20"].colors[(NPRUNS*i+j)%20],
+ color=matplotlib.colormaps["tab20"].colors[(NPRUNS * i + j) % 20],
linestyle="",
marker="*",
markersize=5,
@@ -108,7 +123,7 @@ def compute_peaceman(k_h: float, r_e: float, r_w: float) -> float:
axis.plot(
r_e[-1],
wi_analytical[-1],
- color=matplotlib.colormaps["tab20"].colors[(NPRUNS*i+j)%20],
+ color=matplotlib.colormaps["tab20"].colors[(NPRUNS * i + j) % 20],
linestyle="",
marker=".",
markersize=5,
@@ -117,7 +132,13 @@ def compute_peaceman(k_h: float, r_e: float, r_w: float) -> float:
os.system(f"rm -rf h2o_{NPRUNS*i+j} h2o_{NPRUNS*i+j}.txt")
# Write the configuration files for the comparison in the 3D reservoir
-var = {"flow": FLOW, "perm": PERMEABILITY, "radius": .1, "rate": RATE, "pwd": os.getcwd()}
+var = {
+ "flow": FLOW,
+ "perm": PERMEABILITY,
+ "radius": 0.1,
+ "rate": RATE,
+ "pwd": os.getcwd(),
+}
for name in ["3d_flow_wellmodel", "3d_ml_wellmodel"]:
mytemplate = Template(filename=f"h2o_{name}.mako")
filledtemplate = mytemplate.render(**var)
diff --git a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/h2o_3d_flow_wellmodel.mako b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/h2o_3d_flow_wellmodel.mako
index 00bb6f4..166be55 100644
--- a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/h2o_3d_flow_wellmodel.mako
+++ b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/h2o_3d_flow_wellmodel.mako
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-${flow} --ml-wi-filename="" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+${flow} --ml-wi-filename="" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
co2store base #Model (co2store/h2store)
diff --git a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/h2o_3d_ml_wellmodel.mako b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/h2o_3d_ml_wellmodel.mako
index c51a7ce..562827b 100644
--- a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/h2o_3d_ml_wellmodel.mako
+++ b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/h2o_3d_ml_wellmodel.mako
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-${flow} --ml-wi-filename="${pwd}/water_re.modelPeaceman" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+${flow} --ml-wi-filename="${pwd}/water_re.modelPeaceman" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
co2store base #Model (co2store/h2store)
diff --git a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/h2o_nearwell.mako b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/h2o_nearwell.mako
index 6ba7d3c..e00d045 100644
--- a/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/h2o_nearwell.mako
+++ b/examples/cemracs2023/peaceman_well-model_nn/example_framework_peaceman_water/h2o_nearwell.mako
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-${flow} --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+${flow} --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
co2store base #Model (co2store/h2store)
diff --git a/examples/h2.txt b/examples/h2.txt
index 99084c9..5df95b7 100644
--- a/examples/h2.txt
+++ b/examples/h2.txt
@@ -1,5 +1,5 @@
"""Set the full path to the flow executable and flags"""
-flow --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --ecl-enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
+flow --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0 --enable-drift-compensation=0 --newton-max-iterations=50 --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-5 --relaxed-well-flow-tol=1e-5 --use-multisegment-well=false --enable-tuning=true --enable-opm-rst-file=true --linear-solver=cprw --enable-well-operability-check=false --min-time-step-before-shutting-problematic-wells-in-days=1e-99
"""Set the model parameters"""
h2store base #Model (co2store/h2store)
diff --git a/requirements.txt b/requirements.txt
index 4db8f13..ce7d854 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,10 +1,9 @@
-ecl
resdata
opm
mako
matplotlib
pandas
scipy
-tensorflow
+tensorflow==2.13.0
keras_tuner
scikit-learn
\ No newline at end of file
diff --git a/src/pyopmnearwell/core/pyopmnearwell.py b/src/pyopmnearwell/core/pyopmnearwell.py
index 027e199..41c156b 100644
--- a/src/pyopmnearwell/core/pyopmnearwell.py
+++ b/src/pyopmnearwell/core/pyopmnearwell.py
@@ -5,6 +5,7 @@
import argparse
import os
import pathlib
+from typing import Any
from pyopmnearwell.utils.inputvalues import process_input
from pyopmnearwell.utils.runs import plotting, simulations
@@ -72,7 +73,7 @@ def pyopmnearwell():
help="Scale for the x axis in the figures: 'normal' or 'log' ('normal' by default)",
)
cmdargs = vars(parser.parse_known_args()[0])
- dic = {
+ dic: dict[str, Any] = {
"pat": os.path.split(os.path.dirname(__file__))[0]
} # Path to the pyopmnearwell folder
dic["exe"] = os.getcwd() # Path to the folder of the input.txt file
diff --git a/src/pyopmnearwell/ml/analysis.py b/src/pyopmnearwell/ml/analysis.py
new file mode 100644
index 0000000..fbe3287
--- /dev/null
+++ b/src/pyopmnearwell/ml/analysis.py
@@ -0,0 +1,226 @@
+# pylint: disable=fixme
+"""Analyze the sensitivity of a neural network to its inputs.
+
+Inspiration taken from
+https://f0nzie.github.io/machine_learning_compilation/sensitivity-analysis-for-a-neural-network.html
+
+"""
+
+from __future__ import annotations
+
+import math
+import pathlib
+from typing import Literal, Optional
+
+import numpy as np
+from matplotlib import pyplot as plt
+from tensorflow import keras
+
+from pyopmnearwell.utils import plotting
+
+
+def sensitivity_analysis(
+ model: keras.Model,
+ resolution_1: int = 20,
+ resolution_2: int = 20,
+ mode: (
+ Literal["homogeneous", "random_uniform", "random_normal"] | float
+ ) = "homogeneous",
+) -> tuple[np.ndarray, np.ndarray]:
+ """Perform a sensitivity analysis of a neural network.
+
+ For each input variable, vary from a min to a max value and measure how the output
+ of the network changes. The other input variables are kept constant meanwhile.
+
+ Note: It is assumed that the network has a single output.
+
+ Args:
+ model (keras.Model): Neural network.
+ resolution_1 (int): Number of different values that the fixed inputs take.
+ Equals the number of plotted lines in ``plot_analysis``.
+ resolution_2 (int): Number of different values the variable input takes. Equals
+ the number of datapoints along the x-axis ``plot_analysis``.
+ mode (Literal["homogeneous", "random_uniform", "random_normal"] | float): Sets
+ the sampling mode for the fixed inputs.
+ - "homogeneous":
+ - "random_uniform":
+ - "random_normal"
+ - float: All fixed inputs are set to this value.
+ Default is "homogeneous".
+
+ Returns:
+ tuple[np.ndarray, np.ndarray]: Output and inputs from the sensitivity analysis.
+ First axis is the input variable that is varying, second axis is the variation,
+ third axis is the variation for the fixed variables. The input array contains an
+ additional axis in case of input dimension > 1.
+
+ """
+ # Get the number of input variables.
+ num_inputs = model.input_shape[1]
+
+ # Get the minimum and maximum values for each input variable
+ # min_values = np.min(model.input, axis=0)
+ # max_values = np.max(model.input, axis=0)
+
+ # We assume the model is normalized such that each input variable has min value -1
+ # and max value 1.
+ # TODO: Add options for when this is not the case.
+ min_values: np.ndarray = np.full((num_inputs,), -1)
+ max_values: np.ndarray = np.full((num_inputs,), 1)
+
+ # Initialize the inputs and outputs array.
+ inputs: np.ndarray = np.zeros((num_inputs, resolution_1, resolution_2, num_inputs))
+ outputs: np.ndarray = np.zeros((num_inputs, resolution_1, resolution_2))
+
+ # Initialize random generators.
+ if isinstance(mode, str) and mode.startswith("random"):
+ rng: np.random.Generator = np.random.default_rng()
+
+ # Iterate over all input variables. For each variable and constant input vector, we
+ # create a batch of input arrays and evaluate outside the innermost for loop,
+ # instead of evaluating on single-point batches.
+ for i in range(num_inputs):
+ # Create input arrays for the fixed variables.
+ if mode == "homogeneous":
+ # fixed_inputs: np.ndarray = np.repeat(
+ # np.linspace(min_values, max_values, resolution_1)[..., None],
+ # num_inputs,
+ # axis=-1,
+ # )
+ fixed_inputs: np.ndarray = np.linspace(min_values, max_values, resolution_1)
+ elif mode == "random_uniform":
+ # Get uniform distribution on [0,1) and scale to (min_value, max_value).
+ # pylint: disable-next=possibly-used-before-assignment
+ fixed_inputs = rng.uniform(
+ min_values, max_values, size=(resolution_1, num_inputs)
+ )
+ elif mode == "random_normal":
+ # Set standard deviation s.t. ~95% of the inputs are inside the [min_values,
+ # max_values] interval.
+ fixed_inputs = rng.normal(
+ max_values - min_values,
+ (max_values - min_values) / 4,
+ size=(resolution_1, num_inputs),
+ )
+ elif isinstance(mode, float):
+ fixed_inputs = np.full((resolution_1, num_inputs), mode)
+ else:
+ raise ValueError("'mode' has invalid value")
+ for j, fixed_input in enumerate(fixed_inputs):
+ assert fixed_input.shape == (num_inputs,)
+ fixed_input = np.tile(fixed_input, (resolution_2, 1))
+
+ # Replace i-th input with varying value.
+ fixed_input[..., i] = np.linspace(
+ min_values[i], max_values[i], resolution_2
+ )
+
+ predictions = model.predict(fixed_input)
+
+ assert fixed_input.shape == (resolution_2, num_inputs)
+ assert predictions.shape == (resolution_2, 1)
+
+ outputs[i, j] = predictions.flatten()
+ inputs[i, j] = fixed_input
+
+ return outputs, inputs
+
+
+def plot_analysis(
+ outputs: np.ndarray,
+ inputs: np.ndarray,
+ savepath: str | pathlib.Path,
+ feature_names: Optional[list[str]] = None,
+ main_plot: Optional[tuple[np.ndarray, np.ndarray]] = None,
+ **kwargs,
+) -> None:
+ r"""
+ Plot the analysis of the model outputs against inputs.
+
+ Args:
+ outputs (np.ndarray): The model outputs.
+ inputs (np.ndarray): The model inputs. Shape is ``(num_inputs, resolution_1,
+ resolution_2)``.
+ savepath (str | pathlib.Path): The path to save the plot.
+ feature_names (Optional[list[str]]): The names of the input features. If None,
+ default names (:math:`x_1,x_2,\dots`) will be used. Default is None.
+ main_plot (Optional[tuple[np.ndarray, np.ndarray]]): Add output and input for a
+ main plot that is specifically highlighted. Default is None.
+ **kwargs:
+ - legend (bool): Whether to plot a legend. Default
+ is True.
+
+ Returns:
+ None
+
+ """
+ if feature_names is None:
+ feature_names = [rf"$x_{{{i}}}$" for i in range(inputs.shape[-1])]
+ elif len(feature_names) != inputs.shape[-1]:
+ raise ValueError
+
+ # Max ``max_columns`` (default=3) columns and increase figure size.
+ num_rows: int = math.ceil(outputs.shape[0] / kwargs.get("max_columns", 3))
+ num_columns: int = min(outputs.shape[0], kwargs.get("max_columns", 3))
+
+ fig, axes = plt.subplots(
+ num_rows,
+ num_columns,
+ sharex=True,
+ sharey=True,
+ # NOTE: Setting the figure size here does not work.
+ # figheight=max(num_rows * 3, 5),
+ # figwidth=max(num_columns * 3, 5)
+ # size_inches=(max(num_columns * 3, 5), max(num_rows * 3, 5)),
+ )
+ fig.set_size_inches(w=max(num_columns * 3, 5), h=max(num_rows * 3, 5))
+
+ # mypy complains that ``axes`` might not have an attribute ``.flatten()``, but this
+ # shouldn't be a problem.
+ for i, ax in enumerate(
+ axes.flatten() # type: ignore # pylint: disable=invalid-name
+ ):
+ # If not all rows are full, there are more axes than features. Skip the last few
+ # axes.
+ if i >= len(feature_names):
+ break
+ for j, color in enumerate(
+ # Ignore mypy complaining about ``cm`` not being a module.
+ # pylint: disable-next=no-member
+ plt.cm.Blues(np.linspace(0.3, 0.7, outputs.shape[1])), # type: ignore
+ ):
+ ax.plot(
+ inputs[i, j, :, i],
+ outputs[i, j],
+ color=color,
+ linestyle=":",
+ label=f"{inputs[i, j, 0, i - 1]}",
+ )
+ if main_plot is not None:
+ ax.plot(
+ main_plot[1][i, 0, :, i],
+ main_plot[0][i, 0],
+ color="blue",
+ label="0",
+ )
+ ax.set_title(feature_names[i])
+
+ # Create common axes for all subplots. Add a big axis, hide frame.
+ fig.add_subplot(111, frameon=False)
+ # Set common x and y labels.
+ plt.tick_params(
+ labelcolor="none",
+ which="both",
+ top=False,
+ bottom=False,
+ left=False,
+ right=False,
+ )
+ plt.xlabel("Explanatory")
+ plt.ylabel("Response")
+
+ if kwargs.get("legend", True):
+ # pylint: disable-next=undefined-loop-variable
+ ax.legend(title="Splits", loc="center left", bbox_to_anchor=(1, 0.5))
+
+ plotting.save_fig_and_data(fig, savepath)
diff --git a/src/pyopmnearwell/ml/ensemble.py b/src/pyopmnearwell/ml/ensemble.py
index 67f4196..20412bc 100644
--- a/src/pyopmnearwell/ml/ensemble.py
+++ b/src/pyopmnearwell/ml/ensemble.py
@@ -17,12 +17,12 @@
import numpy as np
import tensorflow as tf
-from ecl.eclfile.ecl_file import open_ecl_file
-from ecl.summary.ecl_sum import EclSum
from mako import exceptions
from mako.template import Template
+from resdata import FileMode
+from resdata.resfile import ResdataFile
+from resdata.summary import Summary
-from pyopmnearwell.utils import units
from pyopmnearwell.utils.formulas import area_squaredcircle, pyopmnearwell_correction
from pyopmnearwell.utils.inputvalues import readthefirstpart, readthesecondpart
from pyopmnearwell.utils.writefile import reservoir_files
@@ -34,7 +34,7 @@
FLAGS = (
" --linear-solver-reduction=1e-5 --relaxed-max-pv-fraction=0"
- + " --ecl-enable-drift-compensation=0 --newton-max-iterations=50"
+ + " --enable-drift-compensation=0 --newton-max-iterations=50"
+ " --newton-min-iterations=5 --tolerance-mb=1e-7 --tolerance-wells=1e-1"
+ " --relaxed-well-flow-tol=1e-3 --use-multisegment-well=false --enable-tuning=true"
+ " --enable-opm-rst-file=true --linear-solver=cprw"
@@ -44,12 +44,20 @@
def create_ensemble(
- runspecs: dict[str, Any], efficient_sampling: Optional[list[str]] = None
+ runspecs: dict[str, Any],
+ efficient_sampling: Optional[list[str]] = None,
+ seed: Optional[int] = None,
) -> list[dict[str, Any]]:
"""Create an ensemble.
- Note: It is assumed that the user provides the variables in the correct units for
- pyopmnearwell.
+ Note:
+ - It is assumed that the user provides the variables in the correct units for
+ pyopmnearwell.
+ - If the variable name starts with ``"PERM"`` or ``"LOG"``, the distribution for
+ random sampling is log uniform.
+ - If the variable name starts with ``"INT"``, the distribution for random
+ sampling is uniform on integers.
+ - Else, the distribution for random sampling is uniformly distributed.
Args:
runspecs (dict[str, Any]): Dictionary with at least the following keys:
@@ -70,6 +78,8 @@ def create_ensemble(
vertical permeabilities. Only 10 layers with 10 samples generate a grid of
10^10 values. By sampling directly instead of generating the grid first, it
is possible to deal with the complexity.
+ seed: (Optional[int]): Seed for the ``np.random.Generator``. Is passed to
+ ``memory_efficient_sample`` as well. Default is ``None``.
Note: The ensemble is generated as the cartesian product of all variable ranges. The
total number of ensemble members is thus the product of all individual
@@ -90,17 +100,23 @@ def create_ensemble(
variables: dict[str, np.ndarray] = OrderedDict()
+ rng: np.random.Generator = np.random.default_rng(seed=seed)
+
# Generate random value ranges for all variables.
logger.info("Generating value ranges for all variables")
for variable, (min_val, max_val, npoints) in runspecs["variables"].items():
- if variable.startswith("PERM"):
- # Generate a log uniform distribution for permeabilities.
+ if variable.startswith("PERM") or variable.startswith("LOG"):
+ # Generate a log uniform distribution for permeabilities and other values
+ # that should be log uniform sampled.
variables[variable] = np.exp(
- np.random.uniform(math.log(min_val), math.log(max_val), npoints)
+ rng.uniform(math.log(min_val), math.log(max_val), npoints)
)
+ elif variable.startswith("INT"):
+ # Generate a uniform distribution on integers.
+ variables[variable] = rng.integers(min_val, max_val, npoints)
else:
# Generate a uniform distribution for all other variables.
- variables[variable] = np.random.uniform(min_val, max_val, npoints)
+ variables[variable] = rng.uniform(min_val, max_val, npoints)
constants: dict[str, float] = runspecs["constants"]
# Differentiate between variables whose data ranges are sampled memory efficiently
@@ -138,6 +154,7 @@ def create_ensemble(
sampled_variables: np.ndarray | list = memory_efficient_sample(
np.array(list(variables_to_sample.values())),
num_members=runspecs["npoints"],
+ seed=seed,
)
else:
sampled_variables = []
@@ -156,15 +173,15 @@ def create_ensemble(
ensemble.append(member)
# Sample a subset of the ensemble if the ensemble size is larger than wanted.
- ensemble_size: int = np.prod(
- [npoints for _, __, npoints in runspecs["variables"].values()]
+ ensemble_size: int = int(
+ np.prod([npoints for _, __, npoints in runspecs["variables"].values()])
)
if runspecs["npoints"] < ensemble_size:
logger.info(
"Sample size is larger than ensemble size. Randomly selecting a subset"
)
- idx = np.random.randint(
+ idx = rng.integers(
len(ensemble),
size=runspecs["npoints"],
)
@@ -180,7 +197,9 @@ def create_ensemble(
return ensemble
-def memory_efficient_sample(variables: np.ndarray, num_members: int) -> np.ndarray:
+def memory_efficient_sample(
+ variables: np.ndarray, num_members: int, seed: Optional[int] = None
+) -> np.ndarray:
"""Sample all variables individually.
Note: Requires that all variables arrays have the same length.
@@ -188,40 +207,43 @@ def memory_efficient_sample(variables: np.ndarray, num_members: int) -> np.ndarr
Args:
variables (np.ndarray), (``shape=(num_variables, len_variables)``): _description_
num_members (int): _description_
+ seed: (Optional[int]): Seed for the ``np.random.Generator``. Default is
+ ``None``.
Returns:
np.ndarray (``shape=()``):
"""
- indices: np.ndarray = np.random.randint(
+ rng: np.random.Generator = np.random.default_rng(seed=seed)
+ indices: np.ndarray = rng.integers(
0, variables.shape[-1], size=(variables.shape[0], num_members)
)
return variables[np.arange(variables.shape[0])[..., None], indices]
def setup_ensemble(
- ensemble_path: pathlib.Path,
+ ensemble_path: str | pathlib.Path,
ensemble: list[dict[str, Any]],
makofile: str | pathlib.Path,
- recalc_grid: bool = False,
- recalc_tables: bool = False,
- recalc_sections: bool = False,
+ **kwargs,
) -> None:
"""Create a deck file for each ensemble member.
Args:
- ensemble_path (pathlib.Path): The path to the ensemble directory.
+ ensemble_path (str | pathlib.Path): The path to the ensemble directory.
ensemble (list[dict[str, Any]]): A list of dictionaries containing the
parameters for each ensemble member. Usually generated by
``create_ensemble``.
makofile (str | pathlib.Path): The path to the Mako template file for the
pyopmnearwell deck for ensemble members.
- recalc_grid (bool, optional): Whether to recalculate ``GRID.INC`` for each
- ensemble member. Defaults to False.
- recalc_tables (bool, optional): Whether to recalculate ``TABLES.INC`` for each
- ensemble member. Defaults to False.
- recalc_sections (bool, optional): Whether to recalculate ``GEOLOGY.INC`` and
- ``REGIONS.INC`` for each ensemble member. Defaults to False.
+ **kwargs: kwargs are passed to ``reservoir_files``. Possible kwargs are:
+ - recalc_grid (bool, optional): Whether to recalculate ``GRID.INC`` for each
+ ensemble member. Defaults to False.
+ - recalc_tables (bool, optional): Whether to recalculate ``TABLES.INC`` for
+ each ensemble member. Defaults to False.
+ - recalc_sections (bool, optional): Whether to recalculate ``GEOLOGY.INC``
+ and ``REGIONS.INC`` for each ensemble member. Defaults to False.
+
Raises:
Exception: If there is an error rendering the Mako template.
@@ -232,6 +254,12 @@ def setup_ensemble(
# Ensure ``ensemble_path`` is a ``Path`` object.
ensemble_path = pathlib.Path(ensemble_path)
+ # Update kwargs with the (future) relative path to the first first ensemble member
+ # from any other ensemble member.
+ kwargs.update(
+ {"inc_folder": pathlib.Path("..") / ".." / "runfiles_0" / "preprocessing"}
+ )
+
logger.info(f"Filling templates for {len(ensemble)} members")
mytemplate: Template = Template(filename=str(makofile))
for i, member in enumerate(ensemble):
@@ -246,7 +274,8 @@ def setup_ensemble(
(ensemble_path / f"runfiles_{i}" / "jobs").mkdir(exist_ok=True)
lol = []
- for row in csv.reader(filledtemplate.split("\n"), delimiter="#"):
+ # Ignore Pylance complaining that the argument might be ``list[bool]``
+ for row in csv.reader(filledtemplate.split("\n"), delimiter="#"): # type: ignore
lol.append(row)
dic, index = readthefirstpart(
lol,
@@ -264,10 +293,11 @@ def setup_ensemble(
else:
reservoir_files(
dic,
- recalc_grid=recalc_grid,
- recalc_tables=recalc_tables,
- recalc_sections=recalc_sections,
- inc_folder=pathlib.Path("..") / ".." / "runfiles_0" / "preprocessing",
+ **kwargs,
+ # recalc_grid=recalc_grid,
+ # recalc_tables=recalc_tables,
+ # recalc_sections=recalc_sections,
+ # inc_folder=pathlib.Path("..") / ".." / "runfiles_0" / "preprocessing",
)
# pyopmnearwell creates these unneeded folders, so we remove them.
try:
@@ -278,26 +308,63 @@ def setup_ensemble(
logger.info(f"Filled templates for {len(ensemble)} members")
+def get_flags(
+ makofile: str | pathlib.Path,
+) -> str:
+ """Extract OPM Flow run flags from a makofile.
+
+
+ Args:
+ makofile (str | pathlib.Path): Path to the makofile.
+
+ Returns:
+ str: All flags that are passed to OPM Flow.
+
+ """
+ # Ensure ``makofile`` is a ``Path`` object.
+ makofile = pathlib.Path(makofile)
+ with makofile.open("r") as f:
+ # Second line has command line arguments.
+ next(iter(f))
+ line: str = next(iter(f))
+ # Strip first argument, which is just the ``flow`` command.
+ command_line_args: list[str] = line.split(" ")[1:]
+ # Strip the newline at the end.
+ return (" ").join(command_line_args).rstrip()
+
+
def run_ensemble(
flow_path: str | pathlib.Path,
- ensemble_path: pathlib.Path,
+ ensemble_path: str | pathlib.Path,
runspecs: dict[str, Any],
ecl_keywords: list[str],
init_keywords: list[str],
summary_keywords: list[str],
num_report_steps: Optional[int] = None,
+ keep_result_files: bool = False,
+ **kwargs,
) -> dict[str, Any]:
- """Run OPM flow for each ensemble member and store data.
+ """Run OPM Flow for each ensemble member and store data.
+
+ Note: The initial time step (i.e., t=0) is always disregarded.
Args:
- flow_path (_type_): _description_
- ensemble_path (pathlib.Path): _description_
+ flow_path (str | pathlib.Path): _description_
+ ensemble_path (str | pathlib.Path): _description_
runspecs (dict[str, Any]): _description_
ecl_keywords (list[str]): _description_
init_keywords (list[str]): _description_
summary_keywords (list[str]): _description_
- num_report_steps (Optional[int], optional): Disregard an ensemble simulation, if
+ num_report_steps (Optional[int], optional): Disregard an ensemble simulation if
it did not run to the last report step. Defaults to None.
+ keep_result_files (bool): Keep result files of all ensemble members, not
+ only the first one. Defaults to False.
+ **kwargs: Possible parameters are:
+ - step_size_time (int): Save data only for every ``step_size_time`` report
+ step. Default is 1.
+ - step_size_cell (int): Save data only for every ``step_size_cell`` grid
+ cell. Default is 1.
+ - flags (str): Flags to run OPM Flow with.
Returns:
dict[str, Any]: _description_
@@ -311,84 +378,105 @@ def run_ensemble(
}
num_disregarded_runs: int = 0
+ # Get **kwargs that determine how many report steps and cells shall be skipped when
+ # extracting the data.
+ step_size_time: int = kwargs.get("step_size_time", 1)
+ step_size_cell: int = kwargs.get("step_size_cell", 1)
+
for i in range(round(runspecs["npoints"] / runspecs["npruns"])):
command = " ".join(
[
f"{flow_path}"
+ f" {ensemble_path / f'runfiles_{j}' / 'preprocessing' / f'RUN_{j}.DATA'}"
+ f" --output-dir={ensemble_path / f'results_{j}'}"
- + f" {FLAGS} & "
+ + f" {kwargs.get('flags', '')} & "
for j in range(runspecs["npruns"] * i, runspecs["npruns"] * (i + 1))
]
)
+ # TODO: Possibly better to use subprocess?
os.system(command + "wait")
for j in range(runspecs["npruns"] * i, runspecs["npruns"] * (i + 1)):
simulation_finished: bool = True
- with open_ecl_file(
- str(ensemble_path / f"results_{j}" / f"RUN_{j}.UNRST")
- ) as ecl_file:
- # Skip result, if the simulation did not run to the last time step.
+ resdata_file: ResdataFile = ResdataFile(
+ str(ensemble_path / f"results_{j}" / f"RUN_{j}.UNRST"),
+ flags=FileMode.CLOSE_STREAM,
+ )
+ # Skip result, if the simulation did not run to the last time step.
+ if (
+ num_report_steps is not None
+ and resdata_file.num_report_steps() < num_report_steps
+ ):
+ simulation_finished = False
+
+ # Check again if the simulation data is available for all time steps.
+ # It seems that sometimes the keyword array has zero report steps, even
+ # though `resdata_file.num_report_steps()` is nonzero.
+ member_data: dict[str, np.ndarray] = {}
+ for keyword in ecl_keywords:
+ # Append the data corresponding to the keyword for all chosen report
+ # steps and cells. Disregard the zeroth time step.
+ member_data[keyword] = np.array(resdata_file.iget_kw(keyword))[
+ 1::step_size_time, ::step_size_cell
+ ]
if (
num_report_steps is not None
- and ecl_file.num_report_steps() < num_report_steps
+ and member_data[keyword].shape[0]
+ < num_report_steps // step_size_time
):
simulation_finished = False
- # Check again if the simulation data is available for all time steps.
- # It seems that sometimes the keyword array has zero report steps, even
- # though `ecl_file.num_report_steps()` is nonzero.
- member_data: dict[str, np.ndarray] = {}
- for keyword in ecl_keywords:
- # Append the data corresponding to the keyword for all report steps
- # and cells. Disregard the zeroth time step.
- member_data[keyword] = np.array(ecl_file.iget_kw(keyword))[1:, ...]
- if (
- num_report_steps is not None
- and member_data[keyword].shape[0] < num_report_steps
- ):
- simulation_finished = False
-
- # Disregard the result if an `inf` value is returned.
- elif np.any(np.isinf(member_data[keyword])):
- simulation_finished = False
+ # Disregard the result if an `inf` value is returned.
+ elif np.any(np.isinf(member_data[keyword])):
+ simulation_finished = False
# Only append data if the simulation finished.
if simulation_finished:
for keyword in ecl_keywords:
data[keyword].append(member_data[keyword])
- with open_ecl_file(
- ensemble_path / f"results_{j}" / f"RUN_{j}.INIT"
- ) as init_file:
+ # Get additional data from init and summary file.
+ if len(init_keywords) > 0:
+ init_file: ResdataFile = ResdataFile(
+ str(ensemble_path / f"results_{j}" / f"RUN_{j}.INIT"),
+ flags=FileMode.CLOSE_STREAM,
+ )
for keyword in init_keywords:
- # Append the data corresponding to the keyword for all cells.
- data[keyword].append(np.array(init_file.iget_kw(keyword)))
-
- summary_file: EclSum = EclSum(
- ensemble_path / f"results_{j}" / f"RUN_{j}.SMSPEC"
- )
- for keyword in summary_keywords:
- # Append the data corresponding to the keyword for all report steps
- # (not for all time steps). The ``*.SMSPEC`` file does not include
- # the zeroth report step. Add a dimension to make it
- # broadcastable to data from the ``*.UNRST`` and ``*.INIT`` files.
- data[keyword].append(
- np.array(
- summary_file.get_values(keyword, report_only=True)[
- ..., None
- ]
+ # Append the data corresponding to the keyword for all chosen
+ # cells.
+ # NOTE: The array has shape ``[1, num_cells]``, hence no axis
+ # needs to be added.
+ data[keyword].append(
+ np.array(init_file.iget_kw(keyword))[::step_size_cell]
)
+
+ if len(summary_keywords):
+ # TODO: Check if lazyload option for ``Summary`` is faster or
+ # slower.
+ summary_file: Summary = Summary(
+ str(ensemble_path / f"results_{j}" / f"RUN_{j}.SMSPEC")
)
- # NOTE: There does not seem to be a way to close an ``EclSum`` object.
- # Also, there is no context manager.
+ for keyword in summary_keywords:
+ # Append the data corresponding to the keyword for all chosen report
+ # steps (not for all time steps). The ``*.SMSPEC`` file does not
+ # include the zeroth report step. Add a dimension to make the array
+ # broadcastable to data from the ``*.UNRST`` and ``*.INIT`` files.
+ data[keyword].append(
+ np.array(
+ summary_file.get_values(keyword, report_only=True)
+ )[::step_size_time, None]
+ )
+ # NOTE: There does not seem to be a way to specify that a
+ # ``Summary`` object shall be closed after use. Also, there is no
+ # context manager for ``Summary`` and ``ResdataFile`` objects.
else:
num_disregarded_runs += 1
logger.info(f"Disregarded ensemble run {j}")
+
# Remove the run files and result folder (except for the first one that
# remains to check if everything went right).
- if j > 0:
+ if not keep_result_files and j > 0:
shutil.rmtree(ensemble_path / f"results_{j}")
shutil.rmtree(ensemble_path / f"runfiles_{j}")
logger.info(f"Disregarded {num_disregarded_runs} of {runspecs['npoints']} runs")
@@ -402,7 +490,7 @@ def calculate_radii(
# num_dims: int = 1,
return_outer_inner: bool = False,
triangle_grid: bool = False,
- theta: float = math.pi / 3,
+ angle: float = math.pi / 3,
) -> np.ndarray | tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Calculates the radii of the cells in a grid grom a given .
@@ -416,8 +504,8 @@ def calculate_radii(
triangle_grid (bool, optional): Whether the grid is a triangle grid. If
True, transform altitudes of the triangle grid to radii of a raidal grid
with equal solution. Defaults to False.
- theta (float, optional): Angle between the horizontal axis and the edge of the
- triangle grid. Defaults to ``math.pi/3``.
+ angle (float, optional): Angle between both sides of the triangle grid. Defaults
+ to ``math.pi/3``.
Returns:
@@ -447,7 +535,7 @@ def calculate_radii(
)
)
if triangle_grid:
- radii *= pyopmnearwell_correction(theta)
+ radii *= pyopmnearwell_correction(angle)
inner_radii: np.ndarray = radii[:-1]
outer_radii: np.ndarray = radii[1:]
radii_t: np.ndarray = (inner_radii + outer_radii) / 2 # unit: [m]
@@ -457,85 +545,69 @@ def calculate_radii(
def calculate_WI(
- data: dict[str, Any],
- runspecs: dict[str, Any],
- num_zcells: int,
+ pressures: np.ndarray,
+ injection_rates: float | np.ndarray,
) -> tuple[np.ndarray, list[int]]:
r"""Calculate the well index (WI) for a given dataset.
- The well index (WI) is calculated using the following formula:
+ The well index (WI) is calculated using the following formula:
.. math::
WI = \frac{q}{{p_w - p_{gb}}}
- Note: In 3D this will more probably than not fail.
+ Note:
+ - The unit of ``WI_array`` will depend on the units of ``pressures`` and
+ ``injection_rates``.
+ - In 3D this might fail. The user is responsible to fix the array shapes before
+ passing to this function.
Args:
- data (dict[str, Any]): Data generated by ``run_ensemble``.
- runspecs (dict[str, Any]): Must contain a key "INJECTION_RATE_PER_SECOND". Unit
- [m^3/s].
+ pressures (np.ndarray): First axis are the ensemble members. Last axis is
+ assumed to be the x-axis. Must contain the well cells (i.e., well pressures
+ values) at ``pressures[...,0]``.
+ injection_rates (float | np.ndarray): Injection rate. If an ``np.ndarray``, it
+ must have shape broadcastable to ``pressures.shape``.
+
Returns:
- WI_array (numpy.ndarray): ``shape=(,num_cells - 1)``
- An array of well index values for each data point in the dataset. In the
- correct unit for OPM. Unit [m^4*s/kg].
+ WI_array (numpy.ndarray): ``shape=(...,num_x_cells - 1)``
+ An array of well index values for each data point in the dataset.
failed_indices (list[int]): Indices for the ensemble members where WI could not
- be computed. E.g., if the simmulation went wrong and the pressure difference is
- zero.
+ be computed. E.g., if the simmulation went wrong and the pressure difference
+ is zero.
Raises:
ValueError: If no data is found for the 'pressure' keyword in the dataset.
+
"""
- # Transform pressure values from [bar] (unit in *.UNRST files) to [Pa] (unit to
- # calculate the WI and internally used by OPM).
- pressure_data = (
- np.array(data["PRESSURE"]) * units.BAR_TO_PASCAL
- ) # ``shape=(runspecs["npoints"], num_time_data - 1, num_grid_cell_data)``;
- # unit [Pa]
- if len(pressure_data) == 0:
- raise ValueError("No data found for the 'pressure' keyword.")
# Calculate WI for each ensemble member.
WI_values: list[np.ndarray] = []
failed_indices: list[int] = []
- for i, member_pressure_data in enumerate(pressure_data):
- p_w = member_pressure_data[
- ..., 0 :: int(member_pressure_data.shape[-1] / num_zcells)
- ] # Bottom hole pressure extracted at the well blocks.
+ if isinstance(injection_rates, float):
+ injection_rates = np.full((pressures.shape[0],), injection_rates)
+
+ for i, (member_pressure, member_injection_rate) in enumerate(
+ zip(pressures, injection_rates) # type: ignore
+ ):
+ p_w = member_pressure[
+ ..., 0
+ ] # Bottom hole/well pressure extracted at the well cells.
p_gb = np.delete(
- member_pressure_data,
- np.arange(
- 0,
- member_pressure_data.shape[-1],
- int(member_pressure_data.shape[-1] / num_zcells),
- ),
- axis=-1,
+ member_pressure, 0, axis=-1
) # Cell pressures of all but the well blocks.
try:
- if "INJECTION_RATE_PER_SECOND" in runspecs["constants"]:
- # Injection rate is constant for all ensemble members.
- WI = runspecs["constants"]["INJECTION_RATE_PER_SECOND"] / (
- np.tile(p_w, int(member_pressure_data.shape[-1] / num_zcells) - 1)
- - p_gb
- ) # ``shape=(runspecs["npoints"], num_time_data - 1, num_grid_cell_data)``;
- # unit [m^4*s/kg]
- elif "INJECTION_RATE_PER_SECOND" in data:
- # Get i-th injection rate in the ensemble.
- WI = data["INJECTION_RATE_PER_SECOND"][i] / (
- np.tile(p_w, int(member_pressure_data.shape[-1] / num_zcells) - 1)
- - p_gb
- ) # ``shape=(runspecs["npoints"], num_time_data - 1, num_grid_cell_data)``;
- # unit [m^4*s/kg]
- if np.any(np.isnan(WI)):
- raise ValueError("WI is nan.")
+ WI: np.ndarray = member_injection_rate / (p_w - p_gb)
+ if np.any(np.logical_or(np.isnan(WI), np.isinf(WI))):
+ raise ValueError("WI is nan/inf.")
WI_values.append(WI)
except ValueError:
failed_indices.append(i)
return (
- np.array(WI_values),
+ np.array(WI_values), # ``shape=(...,num_x_cells - 1)``
failed_indices,
- ) # ``shape=(,num_cells - 1)``; ; unit [m^4*s/kg]
+ )
def extract_features(
@@ -543,7 +615,7 @@ def extract_features(
keywords: list[str],
keyword_scalings: Optional[dict[str, float]] = None,
) -> np.ndarray:
- r"""Extract features into a numpy array.
+ r"""Extract features from a ``run_ensemble`` run into a numpy array.
Note: The features are in the units used in ``*.UNRST`` and ``*.SMSPEC`` files. The
user is responsible for transforming them to the units needed. This may depend on
@@ -562,11 +634,11 @@ def extract_features(
keyword_scalings (Optional[dict[str, float]]): Scalings for the features.
Returns:
- feature_array: (numpy.ndarray): ``shape=(,num_cells - 1)``
- An array of input features each data point in the dataset.
+ feature_array: (numpy.ndarray): ``shape=(ensemble_size, num_report_steps, num_cells, num_features)``
+ An array of input features for each data point in the dataset.
Raises:
- ValueError: If no data is found for the 'pressure' keyword in the dataset.
+ ValueError: If no data is found for one of the keywords.
"""
if keyword_scalings is None:
@@ -594,34 +666,58 @@ def extract_features(
def integrate_fine_scale_value(
radial_values: np.ndarray,
radii: np.ndarray,
- block_sidelength: float,
+ block_sidelengths: float | np.ndarray,
axis: int = -1,
) -> np.ndarray:
- """Integrate a fine scale value across all radial cells lying in a square grid
+ """Integrate a fine scale value across all radial cells covering a square grid
block.
+ This function correctly takes only fractions of the integrated values for radial
+ cells that are only partially inside the square block.
+
Args:
radial_values (np.ndarray): Cell values for the radial cells.
radii (np.ndarray): Array of radii for inner and outer radius of the radial
- cells.
- block_sidelength (float): The sidelength of the square grid block.
+ cells. Has to be ordered from low to high.
+ block_sidelengths (float | np.ndarray): The sidelengths of the square grid
+ blocks. The length of this array determines the new length of the integrated
+ axis.
+ # TODO: Update the tests for this new functionality, i.e., that
+ block_sidelengths determins the return shape.
axis (int): Axis to integrate along.
Returns:
float: The integrated value of fine-scale data.
+ Raise:
+ ValueError: If the radial cells do not cover the square grid block.
+
"""
+ if isinstance(block_sidelengths, float):
+ block_sidelengths = np.array([block_sidelengths])
+
+ # Return 0 for empty radial_values.
+ # TODO: Should this be changed to raise an error as well, if block_sidelength > 0?
if radial_values.shape[0] > 0:
assert radial_values.shape[axis] == radii.shape[0] - 1
else:
return np.zeros_like(radial_values)
+ if np.any(np.sqrt(2 * block_sidelengths**2) > 2 * radii[-1]):
+ raise ValueError(
+ "The disks defined by the radii do not cover all square blocks."
+ )
+
+ integrated_values_lst: list[np.ndarray] = []
# Ignore mypy complaining. ``area_squaredcircle`` returns an np.ndarray in this
# case. This can be removed, once the typing in ``formulas.py`` is more strict.
- cell_areas: np.ndarray = area_squaredcircle( # type: ignore
- radii[1::], block_sidelength
- ) - area_squaredcircle(radii[:-1:], block_sidelength)
- return np.sum(radial_values * cell_areas, axis=axis)
+ # For some reason Pylance thinks that ``block_sidelengths`` is ``int``. Ignore this.
+ for block_sidelength in block_sidelengths: # type: ignore
+ cell_areas: np.ndarray = area_squaredcircle( # type: ignore
+ radii[1::], block_sidelength
+ ) - area_squaredcircle(radii[:-1:], block_sidelength)
+ integrated_values_lst.append(np.sum(radial_values * cell_areas, axis=axis))
+ return np.stack(integrated_values_lst, axis=axis)
def store_dataset(
diff --git a/src/pyopmnearwell/ml/integration.py b/src/pyopmnearwell/ml/integration.py
index a81a529..621d791 100644
--- a/src/pyopmnearwell/ml/integration.py
+++ b/src/pyopmnearwell/ml/integration.py
@@ -2,7 +2,16 @@
# SPDX-FileCopyrightText: 2023 UiB
# SPDX-License-Identifier: GPL-3.0
# pylint: skip-file
-""""Run simulations with machine learned well models."""
+""""Run simulations with machine learned well models.
+
+Check the https://github.com/cssr-tools/ML_near_well/ repo for examples on how to use
+the ``recompile_flow`` and ``run_integration`` functions.
+
+``recompile_flow`` can only used for well models at the moment, however the
+functionality could easily be extended to replace other parts of the OPM simulator
+before recompiling.
+
+"""
from __future__ import annotations
@@ -21,16 +30,12 @@
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)
-dirname: pathlib.Path = pathlib.Path(__file__).parent
-
def recompile_flow(
scalingsfile: pathlib.Path,
opm_path: pathlib.Path,
- StandardWell_impl_template: Literal[
- "co2_3_inputs", "co2_5_inputs", "h2o_2_inputs", "co2_local_stencil"
- ],
- StandardWell_template: Literal["base", "local_stencil"] = "base",
+ StandardWell_impl_template: pathlib.Path,
+ StandardWell_template: pathlib.Path,
stencil_size: int = 3,
local_feature_names: Optional[list[str]] = None,
) -> None:
@@ -44,11 +49,9 @@ def recompile_flow(
scalingsfile (pathlib.Path): Path to the csv file containing the input and
output scalings for the model.
opm_path (pathlib.Path): Path to a OPM installation with ml functionality.
- StandardWell_impl_template (Literal["co2_3_inputs", "co2_5_inputs",
- "h2o_2_inputs", "co2_local_stencil"]): Template for
+ StandardWell_impl_template (pathlib.Path): Template for
``StandardWell_impl.hpp``. Decides the neural network architecture.
- StandardWell_template (Literal["base", "local_stencil"], optional): Template for
- ``StandardWell.hpp``. Defaults to "base".
+ StandardWell_template (pathlib.Path): Template for ``StandardWell.hpp``.
stencil_size (int, optional): The size of the vertical stencil of the model.
Defaults to 3.
local_feature_names (Optional[list[str]], optional): List of local feature names
@@ -71,7 +74,6 @@ def recompile_flow(
opm_well_path: pathlib.Path = (
opm_path / "opm-simulators" / "opm" / "simulators" / "wells"
)
- template_path: pathlib.Path = dirname / ".." / "templates"
# Get the scaling and write it to the C++ mako that integrates nn into OPM.
feature_min: list[float] = []
@@ -112,19 +114,13 @@ def recompile_flow(
"stencil_size": stencil_size,
"cell_feature_names": local_feature_names,
}
- filledtemplate: str = fill_template(
- var,
- filename=str(
- template_path / "standardwell_impl" / f"{StandardWell_impl_template}.mako"
- ),
- )
+
+ # Fill templates and copy into OPM installation.
+ filledtemplate: str = fill_template(var, filename=str(StandardWell_impl_template))
with (opm_well_path / "StandardWell_impl.hpp").open("w", encoding="utf-8") as file:
file.write(filledtemplate)
- shutil.copyfile(
- template_path / "standardwell" / f"{StandardWell_template}.hpp",
- opm_well_path / "StandardWell.hpp",
- )
+ shutil.copyfile(StandardWell_template, opm_well_path / "StandardWell.hpp")
# Recompile flow.
os.chdir(opm_path / "build" / "opm-simulators")
@@ -136,7 +132,7 @@ def run_integration(
) -> None:
"""Runs ``pyopmnearwell`` simulations for the specified runspecs.
- Note: All "variables" need to
+ Note: All "variables" in runspecs need to have the same number of values.
Args:
runspecs (dict[str, Any]): Contains at least two keys "variables" and
@@ -176,7 +172,8 @@ def run_integration(
print(exceptions.text_error_template().render())
raise error
with (savepath / f"run_{i}.txt").open("w", encoding="utf-8") as file:
- file.write(filledtemplate)
+ # We assume that filledtemplate is a string and ignore Pylance complaining.
+ file.write(filledtemplate) # type: ignore
# Use our pyopmnearwell friend to run the 3D simulations and compare the
# results.
diff --git a/src/pyopmnearwell/ml/nn.py b/src/pyopmnearwell/ml/nn.py
index f768df4..c8a0b23 100644
--- a/src/pyopmnearwell/ml/nn.py
+++ b/src/pyopmnearwell/ml/nn.py
@@ -7,20 +7,21 @@
import math
import pathlib
from functools import partial
-from typing import Any, Literal, Optional
+from typing import Any, Literal, Optional, TypeAlias
import keras_tuner
import numpy as np
import pandas as pd
import tensorflow as tf
+from pyopmnearwell.ml.kerasify import export_model
from sklearn.preprocessing import MinMaxScaler
from tensorflow import keras
-from pyopmnearwell.ml.kerasify import export_model
-
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)
+ArrayLike: TypeAlias = tf.Tensor | np.ndarray
+
def get_FCNN(
ninputs: int,
@@ -32,8 +33,7 @@ def get_FCNN(
kernel_initializer: Literal["glorot_normal", "glorot_uniform"] = "glorot_normal",
normalization: bool = False,
) -> keras.Model:
- """
- Returns a fully connected neural network with the specified architecture.
+ """Return a fully connected neural network with the specified architecture.
Args:
ninputs (int): Number of inputs to the model.
@@ -72,19 +72,164 @@ def get_FCNN(
return model
+def get_RNN(
+ ninputs: int,
+ noutputs: int,
+ units: int = 20,
+ saved_model: Optional[str] = None,
+ activation: Literal["sigmoid", "relu", "tanh"] = "tanh",
+ kernel_initializer: Literal["glorot_normal", "glorot_uniform"] = "glorot_uniform",
+) -> keras.Model:
+ """Return a recurrent neural network with the specified architecture.
+
+ Args:
+ ninputs (int): Number of inputs to the model.
+ noutputs (int): Number of outputs from the model.
+ units (int, optional): Size of internal model state. Defaults to 20.
+ hidden_dim (int, optional): Number of neurons in each hidden layer. Defaults to
+ 10.
+ saved_model (str, optional): Path to a saved model to load weights from.
+ Defaults to None.
+ activation (Literal["sigmoid", "relu", "tanh"], optional): Activation function
+ to use in the hidden layers. Defaults to "sigmoid".
+ kernel_initializer (Literal["glorot_normal", "glorot_uniform"], optional):
+ Weight initialization method to use in the hidden layers. Defaults to
+ "glorot_normal".
+
+ Returns:
+ keras.Model: A fully connected neural network.
+
+ """
+ model: keras.Model = keras.Sequential()
+ # RNN as model head.
+ model.add(
+ keras.layers.SimpleRNN(
+ units,
+ input_shape=(None, ninputs),
+ return_sequences=True,
+ activation=activation,
+ kernel_initializer=kernel_initializer,
+ )
+ )
+ # FCNN to as model tail.
+ model.add(keras.layers.TimeDistributed(keras.layers.Dense(10, activation="relu")))
+ model.add(keras.layers.TimeDistributed(keras.layers.Dense(10, activation="relu")))
+ model.add(keras.layers.TimeDistributed(keras.layers.Dense(10, activation="relu")))
+ model.add(keras.layers.Dense(noutputs))
+ if saved_model is not None:
+ model.load_weights(saved_model)
+ return model
+
+
+def get_GRU(
+ ninputs: int,
+ noutputs: int,
+ units: int = 20,
+ saved_model: Optional[str] = None,
+ activation: Literal["sigmoid", "relu", "tanh"] = "tanh",
+ kernel_initializer: Literal["glorot_normal", "glorot_uniform"] = "glorot_uniform",
+) -> keras.Model:
+ """Return a recurrent neural network with the specified architecture.
+
+ Args:
+ ninputs (int): Number of inputs to the model.
+ noutputs (int): Number of outputs from the model.
+ units (int, optional): Size of internal model state. Defaults to 20.
+ hidden_dim (int, optional): Number of neurons in each hidden layer. Defaults to
+ 10.
+ saved_model (str, optional): Path to a saved model to load weights from.
+ Defaults to None.
+ activation (Literal["sigmoid", "relu", "tanh"], optional): Activation function
+ to use in the hidden layers. Defaults to "sigmoid".
+ kernel_initializer (Literal["glorot_normal", "glorot_uniform"], optional):
+ Weight initialization method to use in the hidden layers. Defaults to
+ "glorot_normal".
+
+ Returns:
+ keras.Model: A fully connected neural network.
+
+ """
+ model: keras.Model = keras.Sequential()
+ # RNN as model head.
+ model.add(
+ keras.layers.GRU(
+ units,
+ input_shape=(None, ninputs),
+ return_sequences=True,
+ activation=activation,
+ kernel_initializer=kernel_initializer,
+ )
+ )
+ # FCNN to as model tail.
+ model.add(keras.layers.TimeDistributed(keras.layers.Dense(10, activation="relu")))
+ model.add(keras.layers.TimeDistributed(keras.layers.Dense(10, activation="relu")))
+ model.add(keras.layers.TimeDistributed(keras.layers.Dense(10, activation="relu")))
+ model.add(keras.layers.Dense(noutputs))
+ if saved_model is not None:
+ model.load_weights(saved_model)
+ return model
+
+
+def get_LSTM(
+ ninputs: int,
+ noutputs: int,
+ units: int = 20,
+ saved_model: Optional[str] = None,
+ activation: Literal["sigmoid", "relu", "tanh"] = "tanh",
+ kernel_initializer: Literal["glorot_normal", "glorot_uniform"] = "glorot_uniform",
+) -> keras.Model:
+ """Return a recurrent neural network with the specified architecture.
+
+ Args:
+ ninputs (int): Number of inputs to the model.
+ noutputs (int): Number of outputs from the model.
+ units (int, optional): Size of internal model state. Defaults to 20.
+ hidden_dim (int, optional): Number of neurons in each hidden layer. Defaults to
+ 10.
+ saved_model (str, optional): Path to a saved model to load weights from.
+ Defaults to None.
+ activation (Literal["sigmoid", "relu", "tanh"], optional): Activation function
+ to use in the hidden layers. Defaults to "sigmoid".
+ kernel_initializer (Literal["glorot_normal", "glorot_uniform"], optional):
+ Weight initialization method to use in the hidden layers. Defaults to
+ "glorot_normal".
+
+ Returns:
+ keras.Model: A fully connected neural network.
+
+ """
+ model: keras.Model = keras.Sequential()
+ # RNN as model head.
+ model.add(
+ keras.layers.LSTM(
+ units,
+ input_shape=(None, ninputs),
+ return_sequences=True,
+ activation=activation,
+ kernel_initializer=kernel_initializer,
+ )
+ )
+ # FCNN to as model tail.
+ model.add(keras.layers.TimeDistributed(keras.layers.Dense(10, activation="relu")))
+ model.add(keras.layers.TimeDistributed(keras.layers.Dense(10, activation="relu")))
+ model.add(keras.layers.TimeDistributed(keras.layers.Dense(10, activation="relu")))
+ model.add(keras.layers.Dense(noutputs))
+ if saved_model is not None:
+ model.load_weights(saved_model)
+ return model
+
+
def scale_and_prepare_dataset(
dsfile: str | pathlib.Path,
feature_names: list[str],
- savepath: pathlib.Path,
+ savepath: str | pathlib.Path,
train_split: float = 0.9,
val_split: Optional[float] = 0.1,
test_split: Optional[float] = None,
- conv_input: bool = False,
- conv_output: bool = False,
shuffle: Literal["first", "last", "false"] = "first",
feature_range: tuple[float, float] = (-1, 1),
target_range: tuple[float, float] = (-1, 1),
- scale: bool = False,
+ scale: bool = True,
**kwargs,
) -> (
tuple[tuple[np.ndarray, np.ndarray], tuple[np.ndarray, np.ndarray]]
@@ -103,10 +248,6 @@ def scale_and_prepare_dataset(
train_split (float, optional): Train split. Defaults to 0.9.
val_split (float, optional): Val split. Defaults to 0.1.
test_split (float, optional): Test split. Defaults to None.
- conv_input (bool, optional): Whether the input is for a convolutional neural
- network. Defaults to False.
- conv_output (bool, optional): Whether the output is for a convolutional neural
- network. Defaults to False.
shuffle (Literal["first", "last", "false"], optional): Options for shuffling the
dataset:
- "first": The dataset gets shuffled before the split.
@@ -117,10 +258,10 @@ def scale_and_prepare_dataset(
Defaults to (-1, 1).
target_range (tuple[float, float], optional): Target range of target scaling.
Defaults to (-1, 1)
- scale (bool, optional): Whether to scale the dataset. Defaults to False.
+ scale (bool, optional): Whether to scale the dataset. Defaults to True.
Returns:
- tuple[tuple[np.ndarray, np.ndarray], np.ndarray]
+ tuple[tuple[np.ndarray, np.ndarray], tuple[np.ndarray, np.ndarray]]
| tuple[
tuple[np.ndarray, np.ndarray],
tuple[np.ndarray, np.ndarray],
@@ -136,12 +277,12 @@ def scale_and_prepare_dataset(
ds: tf.data.Dataset = tf.data.Dataset.load(str(dsfile))
features, targets = next(iter(ds.batch(batch_size=len(ds)).as_numpy_iterator()))
- # Reshape features and targets for multidimensional input and output (e.g., for
- # CNNs).
- if conv_input:
- features = features.reshape((-1, features.shape[-1]))
- if conv_output:
- targets = targets.reshape((-1, targets.shape[-1]))
+ # Save feature and targets shape, e.g., for multidimensional data.
+ features_shape: tuple = features.shape
+ targets_shape: tuple = targets.shape
+ # Reshape to one-dimensional data for the scalers to work.
+ features = features.reshape(-1, features.shape[-1])
+ targets = targets.reshape(-1, targets.shape[-1])
if len(feature_names) > features.shape[-1]:
raise ValueError("Too many feature names.")
@@ -284,6 +425,12 @@ def scale_and_prepare_dataset(
(1, train_targets.shape[-1])
)
+ # Reshape to one-dimensional data.
+ train_features = train_features.reshape(-1, features_shape[-1])
+ train_targets = train_targets.reshape(-1, targets_shape[-1])
+ val_features = val_features.reshape(-1, features_shape[-1])
+ val_targets = val_targets.reshape(-1, targets_shape[-1])
+
# Scale the features and targets.
logger.info("Scaling data")
train_features = feature_scaler.transform(train_features)
@@ -291,14 +438,27 @@ def scale_and_prepare_dataset(
val_features = feature_scaler.transform(val_features)
val_targets = target_scaler.transform(val_targets)
+ # Reshape to original shape
+ train_features = train_features.reshape(-1, *features_shape[1:])
+ train_targets = train_targets.reshape(-1, *targets_shape[1:])
+ val_features = val_features.reshape(-1, *features_shape[1:])
+ val_targets = val_targets.reshape(-1, *targets_shape[1:])
+
# Only return test ds if ``test_split > 0``.
# Ignore mypy complaining that ``test_split`` can be None.
if test_split > 0: # type: ignore
test_features, test_targets = next(
iter(test_ds.batch(batch_size=len(test_ds)).as_numpy_iterator())
)
+ test_features = test_features.reshape(-1, features_shape[-1])
+ test_targets = test_targets.reshape(-1, targets_shape[-1])
+
test_features = feature_scaler.transform(test_features)
test_targets = target_scaler.transform(test_targets)
+
+ test_features = test_features.reshape(-1, *features_shape[1:])
+ test_targets = test_targets.reshape(-1, *targets_shape[1:])
+
return (
(train_features, train_targets),
(val_features, val_targets),
@@ -310,9 +470,9 @@ def scale_and_prepare_dataset(
def train(
model: keras.Model,
- train_data: tuple[tf.Tensor, tf.Tensor],
- val_data: tuple[tf.Tensor, tf.Tensor],
- savepath: pathlib.Path,
+ train_data: tuple[ArrayLike, ArrayLike],
+ val_data: tuple[ArrayLike, ArrayLike],
+ savepath: str | pathlib.Path,
lr: float = 0.1,
epochs: int = 500,
bs: int = 64,
@@ -323,14 +483,14 @@ def train(
"mse", "MeanAbsolutePercentageError", "MeanSquaredLogarithmicError"
] = "mse",
recompile_model: bool = True,
- **train_kwargs,
+ **kwargs,
) -> None:
"""Train a tensorflow model on the provided training data and save the best model.
Args:
model (tf.Module): Model to be trained.
- train_data (tuple[tf.Tensor, tf.Tensor]): Training features and targets.
- val_data (tuple[tf.Tensor, tf.Tensor]): Validation features and targets.
+ train_data (tuple[ArrayLike, ArrayLike]): Training features and targets.
+ val_data (tuple[ArrayLike, ArrayLike]): Validation features and targets.
savepath (pathlib.Path): Savepath for models and logging.
lr (float, optional): Initial learning rate. Defaults to 0.1.
epochs (_type_, optional): Training epochs. Defaults to 500.
@@ -346,8 +506,7 @@ def train(
recompile_model (bool, optional): Whether to recompile the model before
training. Can e.g., be set to false, if the model is built and compiled by a
different function. Defaults to True.
- **train_kwargs: Additional keyword arguments to be passed to the ``model.fit()``
- method.
+ **kwargs: Get passed to the ``model.fit()`` method.
Returns:
None
@@ -412,12 +571,13 @@ def train(
early_stopping_callback,
tensorboard_callback,
],
- **train_kwargs,
+ **kwargs,
)
model.save_weights(savepath / "finalmodel")
# Load the best model and save to OPM format.
model.load_weights(savepath / "bestmodel")
+ model.save(savepath / "bestmodel.keras")
if kerasify:
export_model(model, savepath / "WI.model")
@@ -426,6 +586,7 @@ def build_model(
hp: keras_tuner.HyperParameters,
ninputs: int,
noutputs: int,
+ lr_tune: float = 0.1,
) -> tf.Module:
"""Build and compile a FCNN with the given hyperparameters.
@@ -439,7 +600,8 @@ def build_model(
"""
# Get hyperparameters.
- depth: int = hp.Int("depth", min_value=3, max_value=20, step=2)
+ # depth: int = hp.Int("depth", min_value=3, max_value=20, step=2)
+ depth: int = hp.Int("depth", min_value=3, max_value=6, step=2)
hidden_size: int = hp.Int("hidden_size", min_value=5, max_value=50, step=5)
activation: Literal["sigmoid", "relu", "tanh"] = hp.Choice(
"activation", ["sigmoid", "relu", "tanh"]
@@ -462,7 +624,7 @@ def build_model(
)
model.compile(
loss=loss,
- optimizer=tf.keras.optimizers.Adam(learning_rate=0.1),
+ optimizer=tf.keras.optimizers.Adam(learning_rate=lr_tune),
)
return model
@@ -470,25 +632,32 @@ def build_model(
def tune(
ninputs: int,
noutputs: int,
- train_data: tuple[tf.Tensor, tf.Tensor],
- val_data: tuple[tf.Tensor, tf.Tensor],
+ train_data: tuple[ArrayLike, ArrayLike],
+ val_data: tuple[ArrayLike, ArrayLike],
+ savepath: str | pathlib.Path,
objective: Literal["loss", "val_loss"] = "val_loss",
- **tune_kwargs,
-) -> tuple[tf.Module, keras_tuner.Tuner]:
+ max_trials: int = 5,
+ executions_per_trial: int = 1,
+ sample_weight: ArrayLike = np.array([1.0]),
+ lr_tune: float = 0.1,
+ **kwargs,
+) -> tuple[keras.Model, keras_tuner.Tuner]:
"""
Tune the hyperparameters of a neural network model using random search.
Args:
ninputs (int): Number of input features to the model.
noutputs (int): Number of output features to the model.
- train_data (tuple[tf.Tensor, tf.Tensor]): Tuple of training input and target
+ train_data (tuple[ArrayLike, ArrayLike]): Tuple of training input and target
data.
- val_data (tuple[tf.Tensor, tf.Tensor],): Tuple of validation input and target
+ val_data (tuple[ArrayLike, ArrayLike],): Tuple of validation input and target
data.
objective (Literal["loss", "val_loss"], optional): Objective for search.
Defaults to ``"val_loss"``.
- **tune_kwargs: Additional keyword arguments to pass to the tuner's search
- method.
+ max_trials (int): Default is 5.
+ executions_per_trial (int): Default is 1.
+ sample_weight:(ArrayLike): Default is ``np.array([1.0])``.
+ **kwargs: Get passed to the tuner's search method.
Returns:
tf.Module: The model compiled with the best hyperparameters.
@@ -500,13 +669,19 @@ def tune(
"""
# Define the tuner and start a search.
tuner = keras_tuner.RandomSearch(
- hypermodel=partial(build_model, ninputs=ninputs, noutputs=noutputs),
+ hypermodel=partial(
+ build_model,
+ ninputs=ninputs,
+ noutputs=noutputs,
+ lr_tune=lr_tune,
+ ),
objective=objective,
- max_trials=20,
- executions_per_trial=2,
+ max_trials=max_trials,
+ executions_per_trial=executions_per_trial,
overwrite=True,
- directory="my_dir",
- project_name="helloworld",
+ directory=savepath,
+ project_name="tuner",
+ **kwargs,
)
tuner.search_space_summary()
@@ -515,8 +690,14 @@ def tune(
if not isinstance(val_data, tuple) or len(val_data) != 2:
raise ValueError("val_data must be a tuple of two tensors.")
+ # If kwargs contains epochs, this will fail.
tuner.search(
- train_data[0], train_data[1], epochs=20, validation_data=val_data, **tune_kwargs
+ train_data[0],
+ train_data[1],
+ epochs=20,
+ validation_data=val_data,
+ sample_weight=sample_weight,
+ **kwargs,
)
tuner.results_summary()
@@ -527,7 +708,7 @@ def tune(
return model, tuner
-def save_tune_results(tuner: keras_tuner.Tuner, savepath: pathlib.Path) -> None:
+def save_tune_results(tuner: keras_tuner.Tuner, savepath: str | pathlib.Path) -> None:
# Ensure ``savepath`` is a ``Path`` object.
savepath = pathlib.Path(savepath)
@@ -546,19 +727,19 @@ def save_tune_results(tuner: keras_tuner.Tuner, savepath: pathlib.Path) -> None:
def scale_and_evaluate(
model: keras.Model,
- model_input: np.ndarray,
- scalingsfile: pathlib.Path,
-) -> np.ndarray:
+ model_input: ArrayLike,
+ scalingsfile: str | pathlib.Path,
+) -> tf.Tensor:
"""Scale the input, evaluate with the model and scale the output.
Args:
model (tf.keras.Model): A Keras model to evaluate the input with.
- model_input (np.ndarray): Input tensor. Can be a batch.
- scalingsfile (pathlib.Path): The path to the CSV file containing the
+ model_input (ArrayLike): Input tensor. Can be a batch.
+ scalingsfile (str | pathlib.Path): The path to the CSV file containing the
scaling parameters for MinMaxScaling.
Returns:
- np.ndarray: The model's output, scaled back to the original range.
+ tf.Tensor: The model's output, scaled back to the original range.
Raises:
FileNotFoundError: If ``scalingsfile`` does not exist.
@@ -601,23 +782,53 @@ def scale_and_evaluate(
feature_scaler: MinMaxScaler = MinMaxScaler(feature_range)
feature_scaler.data_min_ = np.array(feature_min)
feature_scaler.data_max_ = np.array(feature_max)
- feature_scaler.scale_ = (feature_range[1] - feature_range[0]) / (
- feature_scaler.data_max_ - feature_scaler.data_min_
- )
+ feature_scaler.scale_ = (
+ feature_range[1] - feature_range[0]
+ ) / handle_zeros_in_scale(feature_scaler.data_max_ - feature_scaler.data_min_)
feature_scaler.min_ = (
feature_range[0] - feature_scaler.data_min_ * feature_scaler.scale_
)
target_scaler: MinMaxScaler = MinMaxScaler(target_range)
target_scaler.data_min_ = np.array(target_min)
target_scaler.data_max_ = np.array(target_max)
- target_scaler.scale_ = (target_range[1] - target_range[0]) / (
+ target_scaler.scale_ = (target_range[1] - target_range[0]) / handle_zeros_in_scale(
target_scaler.data_max_ - target_scaler.data_min_
)
target_scaler.min_ = (
target_range[0] - target_scaler.data_min_ * target_scaler.scale_
)
+ # Save feature and targets shape and reshape to one-dimensional data, e.g., for
+ # multidimensional data. The scaler accepts at most two dimensions.
+ input_shape: tuple = model_input.shape
+ model_input = model_input.reshape(-1, input_shape[-1])
+
+ scaled_input: np.ndarray = feature_scaler.transform(model_input)
+ scaled_input = scaled_input.reshape(input_shape)
+
# Run model.
- scaled_input: tf.Tensor = feature_scaler.transform(model_input)
output: tf.Tensor = model(scaled_input)
- return target_scaler.inverse_transform(output)
+
+ output_shape: tuple = output.shape
+ output = tf.reshape(output, (-1, output_shape[-1]))
+
+ unscaled_output: tf.Tensor = target_scaler.inverse_transform(output)
+ unscaled_output = tf.reshape(unscaled_output, output_shape)
+
+ return unscaled_output
+
+
+def handle_zeros_in_scale(scale: ArrayLike) -> np.ndarray:
+ """Set scales of near constant features to 1.
+
+ Note: This behavior is in line with `sklearn.preprocessing.MinMaxScaler`.
+
+ Args:
+ scale (ArrayLike): The scale array.
+
+ Returns:
+ np.ndarray: The modified scale array.
+
+ """
+ # ``atol`` must be very low s.t. this works for permeability in [m^2] for example.
+ return np.where(np.isclose(scale, 0.0, atol=1e-20), np.ones_like(scale), scale)
diff --git a/src/pyopmnearwell/ml/data.py b/src/pyopmnearwell/ml/resdata_dataset.py
similarity index 71%
rename from src/pyopmnearwell/ml/data.py
rename to src/pyopmnearwell/ml/resdata_dataset.py
index 0009201..e67b238 100644
--- a/src/pyopmnearwell/ml/data.py
+++ b/src/pyopmnearwell/ml/resdata_dataset.py
@@ -1,3 +1,4 @@
+# pylint: disable=fixme
"""This module provides functionality to parse ``*.UNRST`` files for given keywords and
transform the extracted values into a tensorflow dataset.
@@ -5,7 +6,8 @@
different from the default one. The lines that need to be changed are marked with
# MANUAL CHANGES.
-Deprecated: This module is deprecated in favor of the ``ensemble`` module.
+Deprecated: This module is deprecated in favor of the ``ensemble`` module which can run
+an ensemble of pyopmnearwell decks AND extract data afterwards.
"""
@@ -18,20 +20,21 @@
import numpy as np
import tensorflow as tf
-from ecl.eclfile.ecl_file import EclFile, open_ecl_file
+from resdata import FileMode
+from resdata.resfile import ResdataFile
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)
-class EclDataSet: # pylint: disable=R0902
- """Generate samples for a ``tf.data.Dataset`` from a folder of ``*.UNRST`` files.
+class ResDataSet: # pylint: disable=R0902
+ """Generate samples for a ``tf.data.Dataset`` from a folder of ``.UNRST`` files.
Example:
After instantiation of the class, it can be passed to
``tf.data.Dataset.from_generator()``, to create a dataset that ``tensorflow``
can work with.
- >>> data = EclDataSet(path, input_kws, target_kws)
+ >>> data = ResDataSet(path, input_kws, target_kws)
>>> data.read_data()
>>> ds = tf.data.Dataset.from_generator(
>>> data,
@@ -40,10 +43,11 @@ class EclDataSet: # pylint: disable=R0902
To save the dataset, use
>>> ds.save(path)
- Afterwards, the ``*.UNRST`` files used to generated the dataset can be deleted.
+ Afterwards, the ``.UNRST`` files used to generated the dataset can be deleted.
"""
+ # Typing for instance attributes.
features: tf.Tensor
"""Stores all inputs of the dataset.
@@ -62,7 +66,7 @@ def __init__( # pylint: disable=R0913
path: str,
input_kws: list[str],
target_kws: list[str],
- file_format: Literal["ecl", "opm"] = "ecl",
+ file_format: Literal["resdata", "opm"] = "resdata",
dtype=tf.float32,
shuffle_on_epoch_end: bool = False,
read_data_on_init: bool = True,
@@ -71,17 +75,17 @@ def __init__( # pylint: disable=R0913
Parameters:
path: _description_
- input_kws: Keywords for attributes of the ``EclFile`` that shall become
+ input_kws: Keywords for attributes of the ``.UNRST`` file that shall become
model input.
- target_kws: Keywords for attributes of the ``EclFile`` that shall become
+ target_kws: Keywords for attributes of the ``.UNRST`` file that shall become
targets for model training.
- type: _description_. Defaults to ``"ecl"``.
- read_data_on_init: Reads data from ``*.UNRST`` files in ``path`` on
+ type: _description_. Defaults to ``"resdata"``.
+ read_data_on_init: Reads data from ``.UNRST`` files in ``path`` on
instantiation. Disable for testing/debugging. Defaults to ``True``.
Warning:
- As of now, ``type`` is always assumed to be ``ecl``. ``opm`` is not
+ As of now, ``type`` is always assumed to be ``resdata``. ``opm`` is not
implemented yet.
Returns:
@@ -95,28 +99,30 @@ def __init__( # pylint: disable=R0913
if read_data_on_init:
self.read_data()
self.shuffle_on_epoch_end: bool = shuffle_on_epoch_end
- self.file_format: Literal["ecl", "opm"] = file_format
+ self.file_format: Literal["resdata", "opm"] = file_format
def read_data(self):
- """Create a ``tensorflow`` dataset from a folder of ``ecl`` or ``opm`` files."""
+ """Create a ``tensorflow`` dataset from a folder of ``resdata`` or ``opm`` files."""
logger.info("Generating datapoints...")
_features_lst: list[tf.Tensor] = []
_targets_lst: list[tf.Tensor] = []
for filename in os.listdir(self.path):
if filename.endswith("UNRST"):
- with open_ecl_file(os.path.join(self.path, filename)) as ecl_file:
- try:
- feature, target = self.EclFile_to_datapoint(ecl_file)
- _features_lst.append(feature)
- _targets_lst.append(target)
- logger.info( # pylint: disable=W1203
- f"Generated a datapoint from {filename}."
- )
- except KeyError as keyerror:
- logger.info( # pylint: disable=W1203
- f"{filename} has no keyword {keyerror}."
- )
- continue
+ resdata_file: ResdataFile = ResdataFile(
+ os.path.join(self.path, filename), flags=FileMode.CLOSE_STREAM
+ )
+ try:
+ feature, target = self.ResdataFile_to_datapoint(resdata_file)
+ _features_lst.append(feature)
+ _targets_lst.append(target)
+ logger.info( # pylint: disable=W1203
+ f"Generated a datapoint from {filename}."
+ )
+ except KeyError as keyerror:
+ logger.info( # pylint: disable=W1203
+ f"{filename} has no keyword {keyerror}."
+ )
+ continue
if len(_features_lst) > 0 and len(_targets_lst) > 0:
# Transform the lists into a tensor
@@ -138,25 +144,24 @@ def read_data(self):
Generated an empty dataset for now."""
)
- def EclFile_to_datapoint( # pylint: disable= C0103
- self, ecl_file: EclFile
+ def ResdataFile_to_datapoint( # pylint: disable= C0103
+ self, resdata_file: ResdataFile
) -> tuple[tf.Tensor, tf.Tensor]:
- """Extract values from an ``EclFile`` to form an (input, target) tuple of
- tensors.
+ """Extract values from an ``ResdataFile`` object to form an (input, target)
+ tuple of tensors.
-
- Parameters:
- ecl_file: _description_
+ Args:
+ ecl_file (EclFile): _description_
Raises:
- KeyError: If ``ecl_file`` does not have either of the keywords in
+ KeyError: If ``resdata_file`` does not have either of the keywords in
``self.input_kws`` or ``self.target_kws``
Returns:
- A tuple containing the input and target tensor. The former has shape
- ``(ecl_file.num_report_steps(), num_cells, len(input_kws))``, while the
- latter has shape
- ``(ecl_file.num_report_steps(), num_cells, len(target_kws))``.
+ tuple[tf.Tensor, tf.Tensor]: A tuple containing the input and target tensor.
+ The former has shape ``(resdata_file.num_report_steps(), num_cells,
+ len(input_kws))``, while the latter has shape
+ ``(resdata_file.num_report_steps(), num_cells, len(target_kws))``.
"""
# Only add the datapoint if all input features and targets are
@@ -167,13 +172,13 @@ def EclFile_to_datapoint( # pylint: disable= C0103
# corresponding to the keyword, with shape
# ``(ecl_file.num_report_steps(), num_cells)``.
for input_kw in self.input_kws:
- if ecl_file.has_kw(input_kw):
- feature[input_kw] = np.array(ecl_file.iget_kw(input_kw))
+ if resdata_file.has_kw(input_kw):
+ feature[input_kw] = np.array(resdata_file.iget_kw(input_kw))
else:
raise KeyError(input_kw)
for target_kw in self.target_kws:
- if ecl_file.has_kw(target_kw):
- target[target_kw] = np.array(ecl_file.iget_kw(target_kw))
+ if resdata_file.has_kw(target_kw):
+ target[target_kw] = np.array(resdata_file.iget_kw(target_kw))
else:
raise KeyError(target_kw)
# Stack the different input and target properties into one input
@@ -220,13 +225,13 @@ def on_epoch_end(self):
def main(args): # pylint: disable=W0621
"""Create a dataset from the given arguments and store it"""
- data = EclDataSet(args.path, args.input_kws, args.target_kws, args.file_format)
+ data = ResDataSet(args.path, args.input_kws, args.target_kws, args.file_format)
assert len(data) > 0
dataset = tf.data.Dataset.from_generator(
data,
output_signature=(
- tf.TensorSpec.from_tensor(data[0][0]),
- tf.TensorSpec.from_tensor(data[0][1]),
+ tf.TensorSpec.from_tensor(data[0][0]), # type: ignore
+ tf.TensorSpec.from_tensor(data[0][1]), # type: ignore
),
)
# Manually set the dataset cardinality.
diff --git a/src/pyopmnearwell/ml/upscale.py b/src/pyopmnearwell/ml/upscale.py
new file mode 100644
index 0000000..973e29e
--- /dev/null
+++ b/src/pyopmnearwell/ml/upscale.py
@@ -0,0 +1,475 @@
+# pylint: disable=fixme,missing-function-docstring,no-member, pointless-string-statement
+"""Functionality to upscale data from an ensemble run on a radial grid to a cartesian
+grid.
+
+Note: pylint: no-member is disabled, because it complains about the ``BaseUpscaler``
+missing instance attributes, which are taken care of by the ``Upscaler`` protocol.
+pylint: pointless-string-statement is disabled, as it complains an attribute docstring
+in the ``Upscaler`` protocol.
+
+"""
+
+import math
+import pathlib
+from abc import ABC, abstractmethod
+from typing import Protocol
+
+import numpy as np
+
+# ``tqdm`` is not a dependency. Up to the user to install it.
+try:
+ # Avoid some mypy trouble.
+ import tqdm.autonotebook as tqdm # type: ignore
+except ImportError:
+ _IS_TQDM_AVAILABLE: bool = False
+else:
+ _IS_TQDM_AVAILABLE = True
+
+from pyopmnearwell.ml import ensemble
+from pyopmnearwell.utils import formulas, units
+
+
+class Upscaler(Protocol):
+ """Protocol class for upscalers.
+
+ This class is used for typing of abstract attributes of ``BaseUpscaler``. MyPy will
+ check if subclasses of ``BaseUpscaler`` implement the following instance attributes:
+ - ``num_timesteps``
+ - ``num_layers``
+ - ``num_zcells``
+ - ``num_xcells``
+ - ``single_feature_shape``
+ However, in comparison to missing abstract functions, no runtime error will be
+ raised if they are missing.
+
+ See, e.g., https://stackoverflow.com/a/75253719 for an explanation.
+
+ Note: As of now (Python 3.11), each instance method of ``BaseUpscaler`` or a
+ subclass making use of one of the attributes, needs to have its self argumented
+ annotated with ``Upscaler``. In future python versions it should be possible to
+ use ``class BaseUpscaler[Upscaler](ABC):...`` instead and remove these
+ annotations.
+
+ """
+
+ @property
+ def num_timesteps(self) -> int: ...
+
+ @property
+ def num_layers(self) -> int: ...
+
+ @property
+ def num_zcells(self) -> int: ...
+
+ @property
+ def num_xcells(self) -> int: ...
+
+ @property
+ def single_feature_shape(self) -> tuple: ...
+
+ @property
+ def angle(self) -> float: ...
+
+ """Angle of the cake radial grid. Default is 60°."""
+
+
+class BaseUpscaler(ABC):
+ """Extract and upscale data from an array of ensemble data.
+
+ This base class provides several methods to extract features from fine-scale radial
+ simulations and upscale to coarse cartesian cells. Depending on the type of data,
+ this is done by averaging/summing/etc. values along all cells that correspond to a
+ coarse cell.
+
+ Additionally, the sparsity of the dataset can be increased by taking only some
+ timesteps/horizontal cells.
+
+ The upscaled data is usually provided in form of two ``np.ndarrays``, one for
+ features and one for targets.
+
+ Subclasses need to implement ``__init__`` and (if needed) ``create_ds`` methods.
+
+ The feature array will have shape ``(num_ensemble_runs, num_timesteps/step_size_t,
+ num_layers, num_xcells/step_size_x, num_features)``.
+ The target array will have shape ``(num_ensemble_runs, num_timesteps/step_size_t,
+ num_layers, num_xcells/step_size_x, 1)``
+
+ Note: All methods assume that all cells have the same height. if this is not the
+ case, the methods must be overridden.
+
+ """
+
+ @abstractmethod
+ def __init__(self: Upscaler): # pylint: disable=missing-function-docstring
+ return
+
+ @abstractmethod
+ def create_ds(self: Upscaler): # pylint: disable=missing-function-docstring
+ return
+
+ def reduce_data_size(
+ self: Upscaler,
+ feature: np.ndarray,
+ step_size_x: int = 1,
+ step_size_t: int = 1,
+ random: bool = False,
+ ) -> np.ndarray:
+ """Reduce the size of the input feature array by selecting elements with a
+ fixed step size.
+
+ Args:
+ feature (np.ndarray): The input feature array.
+ step_size_x (int, optional): The step size for the x-axis. Defaults to 1.
+ step_size_t (int, optional): The step size for the t-axis. Defaults to 1.
+ random (bool, optional): If True, select elements randomly instead of using
+ a fixed step size. Defaults to False. Not implemented yet.
+
+ Returns:
+ np.ndarray: The reduced feature array.
+
+ """
+ # TODO: Add option to have a random selection of elements for each member,
+ # instead of a fixed stepsize. This needs to be the same for each feature, hence
+ # the current idea does not work.
+ if random:
+ pass
+ return feature[:, ::step_size_t, ::, ::step_size_x]
+
+ def get_vertically_averaged_values(
+ self: Upscaler,
+ features: np.ndarray,
+ feature_index,
+ disregard_first_xcell: bool = True,
+ ) -> np.ndarray:
+ """Average features vertically inside each layer.
+
+ Args:
+ features (np.ndarray): _description_
+ feature_index (int): _description_.
+ disregard_first_xcell (bool): __description__. Default is True.
+
+ Returns:
+ np.ndarray (``shape = (num_ensemble_runs, num_timesteps, num_layers,
+ num_xcells)``):
+
+ """
+ # Innermost cells (well cells) get disregarded.
+ feature: np.ndarray = np.average(features[..., feature_index], axis=-2)
+
+ if disregard_first_xcell:
+ feature = feature[..., 1:]
+ assert feature.shape == self.single_feature_shape
+
+ else:
+ assert feature.shape[:-1] == self.single_feature_shape[:-1]
+ assert feature.shape[-1] == self.single_feature_shape[-1] + 1
+
+ return feature
+
+ def get_radii(
+ self: Upscaler, radii_file: pathlib.Path
+ ) -> tuple[np.ndarray, np.ndarray]:
+ """Get full list of cell radii."""
+ cell_center_radii, inner_radii, outer_radii = ensemble.calculate_radii(
+ radii_file,
+ num_cells=self.num_xcells + 1,
+ return_outer_inner=True,
+ triangle_grid=True,
+ angle=self.angle,
+ )
+ # Innermost cells (well cells) get disregarded. If the cell_boundary_radii are
+ # used for integrating, this is actually needed, as the content of the well does
+ # not count towards the well block, hence the well radius should not be
+ # integrated.
+ cell_center_radii = cell_center_radii[1:]
+ inner_radii = inner_radii[1:]
+ outer_radii = outer_radii[1:]
+
+ cell_boundary_radii: np.ndarray = np.append(inner_radii, outer_radii[-1])
+ assert cell_center_radii.shape == (self.num_xcells,)
+ assert cell_boundary_radii.shape == (self.num_xcells + 1,)
+ return cell_center_radii, cell_boundary_radii
+
+ def get_timesteps(self: Upscaler, simulation_length: float) -> np.ndarray:
+ """_summary_
+
+ Returns:
+ np.ndarray (``shape = (num_timesteps,) ``): Unit: [d].
+
+ """
+ timesteps: np.ndarray = np.linspace(0, simulation_length, self.num_timesteps)
+ assert timesteps.shape == self.num_timesteps
+ return timesteps
+
+ def get_horizontically_integrated_values( # pylint: disable=too-many-arguments
+ self: Upscaler,
+ features: np.ndarray,
+ cell_center_radii: np.ndarray,
+ cell_boundary_radii: np.ndarray,
+ feature_index: int,
+ disregard_first_xcell: bool = True,
+ ):
+ """Integrate feature horizontically along layers and divide by equivalent
+ cartesian block area.
+
+ Note:
+ - Before integrating, the feature is averaged vertically inside each layer.
+ - As the feature is averaged and not summed, the integration takes place in
+ 2D with vertically averaged values. Hence it suffices to divide by area
+ and not by volume.
+
+ Args:
+ features (np.ndarray): _description_
+ cell_center_radii (np.ndarray):
+ cell_boundary_radii (np.ndarray):
+ feature_index (int): _description_. Default is 1.
+ disregard_first_xcell (bool): __description__. Default is True.
+
+ Returns:
+ np.ndarray (``shape = (num_ensemble_runs, num_timesteps, num_layers,
+ num_xcells)``): Features values average for each cell.
+
+ """
+ # Average along vertical cells in a layer.
+ feature: np.ndarray = np.average(features[..., feature_index], axis=-2)
+
+ if disregard_first_xcell:
+ feature = feature[..., 1:]
+
+ # Integrate horizontically along layers and divide by equivalent cartesian block
+ # area.
+ block_sidelengths: np.ndarray = formulas.cell_size(cell_center_radii) # type: ignore
+ # feature_lst: list[np.ndarray] = []
+ # for i in range(feature.shape[-1]):
+ # feature_lst.append(
+ # ensemble.integrate_fine_scale_value(
+ # feature[..., : i + 1],
+ # cell_boundary_radii[: i + 2],
+ # block_sidelengths[: i + 1],
+ # )
+ # / (block_sidelengths[i] ** 2)
+ # )
+ # averaged_feature: np.ndarray = np.concatenate(feature_lst, axis=-1)
+ # The horizontal dimension is the last axis of each feature, hence we pass
+ # ``axis=-1``.
+ integrated_feature: np.ndarray = ensemble.integrate_fine_scale_value(
+ feature, cell_boundary_radii, block_sidelengths, axis=-1
+ ) / (block_sidelengths**2)
+
+ assert integrated_feature.shape == self.single_feature_shape
+ return integrated_feature
+
+ def get_homogeneous_values(
+ self: Upscaler, features, feature_index, disregard_first_xcell: bool = True
+ ):
+ """Get a feature that is homogeneous inside a layer.
+
+ Note: Since the feature is equal inside a layer, this method takes the first
+ value for each layer.
+
+ Args:
+ features (np.ndarray): _description_
+ feature_index (int): _description_.
+ disregard_first_xcell (bool): __description__. Default is True.
+
+ Returns:
+ np.ndarray (``shape = (num_ensemble_runs, num_timesteps, num_layers,
+ num_xcells)``):
+
+ """
+ # Innermost cells (well cells) get disregarded.
+ feature: np.ndarray = features[..., feature_index][..., 0, :]
+
+ if disregard_first_xcell:
+ feature = feature[..., 1:]
+
+ assert feature.shape == self.single_feature_shape
+ return feature
+
+ def get_analytical_PI( # pylint: disable=invalid-name
+ self: Upscaler,
+ permeabilities: np.ndarray,
+ cell_heights: np.ndarray,
+ radii: np.ndarray,
+ well_radius: float,
+ ) -> np.ndarray:
+ """_summary_
+
+ _extended_summary_
+
+ Args:
+ permeabilities (np.ndarray): Unit has to be [m^2]!
+ cell_heights (np.ndarray): Unit [m].
+ radii (np.ndarray): _description_
+ well_radius (float): _description_
+
+ Returns:
+ np.ndarray: _description_
+ """
+ analytical_PI: np.ndarray = formulas.peaceman_matrix_WI( # type: ignore
+ k_h=permeabilities * cell_heights,
+ r_e=radii,
+ r_w=well_radius,
+ )
+ assert analytical_PI.shape == self.single_feature_shape
+ return analytical_PI
+
+ # pylint: disable-next=invalid-name, too-many-arguments, too-many-locals
+ def get_analytical_WI(
+ self: Upscaler,
+ pressures: np.ndarray,
+ saturations: np.ndarray,
+ permeabilities: np.ndarray,
+ temperature: float,
+ surface_density: float,
+ radii: np.ndarray,
+ well_radius: float,
+ # pylint: disable-next=invalid-name
+ OPM: pathlib.Path,
+ ) -> np.ndarray:
+ """_summary_
+
+ _extended_summary_
+
+ Args:
+ pressures (np.ndarray): _description_
+ saturations (np.ndarray): _description_
+ permeabilities (np.ndarray): Unit has to be [mD]!
+ temperature (float): _description_
+ surface_density (float): _description_
+ radii (np.ndarray): _description_
+ well_radius (float): _description_
+ OPM (pathlib.Path): _description_
+
+ Returns:
+ np.ndarray: _description_
+
+ """
+ densities_lst: list[list[float]] = []
+ viscosities_lst: list[list[float]] = []
+ # ``formulas.co2brinepvt`` calls CO2BRINEPVT from OPM Flow for each pressure,
+ # which is quite inefficient. If tqdm is available, we show a progress bar.
+ if _IS_TQDM_AVAILABLE:
+ pressures_bar = tqdm.tqdm(pressures.flatten())
+ else:
+ pressures_bar = pressures.flatten()
+
+ for pressure in pressures_bar:
+ # Evaluate density and viscosity.
+ density_tuple: list[float] = []
+ viscosity_tuple: list[float] = []
+
+ for phase in ["water", "CO2"]:
+ density_tuple.append(
+ formulas.co2brinepvt(
+ pressure=pressure,
+ temperature=temperature + units.CELSIUS_TO_KELVIN,
+ phase_property="density",
+ phase=phase, # type: ignore
+ OPM=OPM,
+ )
+ )
+
+ viscosity_tuple.append(
+ formulas.co2brinepvt(
+ pressure=pressure,
+ temperature=temperature + units.CELSIUS_TO_KELVIN,
+ phase_property="viscosity",
+ phase=phase, # type: ignore
+ OPM=OPM,
+ )
+ )
+ densities_lst.append(density_tuple)
+ viscosities_lst.append(viscosity_tuple)
+
+ densities_shape = list(pressures.shape)
+ densities_shape.extend([2])
+ densities: np.ndarray = np.array(densities_lst).reshape(densities_shape)
+ viscosities: np.ndarray = np.array(viscosities_lst).reshape(densities_shape)
+
+ # Calculate the well index from Peaceman. The analytical well index is in [m*s],
+ # hence we need to devide by surface density to transform to [m^4*s/kg].
+ # pylint: disable-next=invalid-name
+ analytical_WI: np.ndarray = ( # type: ignore
+ # Ignore unsupported operand types for *. Fixing this would be quite
+ # complex.
+ formulas.two_phase_peaceman_WI( # type: ignore
+ k_h=permeabilities
+ * units.MILIDARCY_TO_M2
+ * (
+ self.num_zcells / self.num_layers
+ ), # TODO: This is FALSE, it must be height instead!
+ r_e=radii,
+ r_w=well_radius,
+ rho_1=densities[..., 0],
+ mu_1=viscosities[..., 0],
+ k_r1=(1 - saturations) ** 2,
+ rho_2=densities[..., 1],
+ mu_2=viscosities[..., 1],
+ k_r2=saturations**2,
+ )
+ / surface_density
+ )
+
+ # TODO: Add this assert again. It's turned off s.t. the analytical WI can be
+ # computed after size reduction of features.
+ # assert analytical_WI.shape == self.single_feature_shape
+ return analytical_WI
+
+ def get_data_WI( # pylint: disable=invalid-name
+ self: Upscaler,
+ features: np.ndarray,
+ pressure_index: int,
+ inj_rate_index: int,
+ angle: float = math.pi / 3,
+ ) -> np.ndarray:
+ """Calculate data-driven WI from pressure and flow rate.
+
+ Similar functionality to ``ensemble.calculate_WI``, but can additionally treat
+ multiple vertical cells in a layer correctly.
+
+ Note:
+ - Pressures get averaged over each layer, injection rates get summed over
+ each layer.
+ - The method automatically scales the near-well injection rate from the cake
+ grid with angle ``angle`` to a 360° well. Furthermore, the rate is scaled
+ from rate-per-day (which the results are in) to rate-per-second, which OPM
+ Flow uses internally for the WI.
+
+
+ Args:
+ features (np.ndarray): _description_
+ pressure_index (int): _description_
+ inj_rate_index (int): _description_
+
+ Returns:
+ np.ndarray: _description_
+
+ """
+ # Take the pressure values of the well blocks as bhp. Average along each layer.
+ bhps: np.ndarray = np.average(features[..., pressure_index], axis=-2)[..., 0][
+ ..., None
+ ] # ``shape = (num_completed_runs, num_timesteps, num_layers, 1)``
+
+ # Ge the pressure values of all other blocks. Average along each layer.
+ pressures: np.ndarray = np.average(features[..., pressure_index], axis=-2)[
+ ..., 1:
+ ] # ``shape = (num_completed_runs, num_timesteps, num_layers, num_xcells)``
+
+ # Get the individual injection rates per second for each layer. Sum across a
+ # layer to get the rates for a full layer. Multiply by
+ # ``(math.pi * 2 / angle)`` to transform from a cake of given ``angle`` to a
+ # full radial model and convert to rate per second.
+ injection_rate_per_second_per_cell: np.ndarray = (
+ np.sum(features[..., inj_rate_index], axis=-2)[..., 0]
+ * (math.pi * 2 / angle)
+ * units.Q_per_day_to_Q_per_seconds
+ )[
+ ..., None
+ ] # ``shape = (num_completed_runs, num_timesteps, num_layers, 1)``
+
+ # Check that we do not divide by zero.
+ assert np.all(bhps - pressures)
+ WI_data: np.ndarray = injection_rate_per_second_per_cell / (bhps - pressures)
+ assert WI_data.shape == self.single_feature_shape
+ return WI_data
diff --git a/src/pyopmnearwell/ml/utils.py b/src/pyopmnearwell/ml/utils.py
new file mode 100644
index 0000000..2a24979
--- /dev/null
+++ b/src/pyopmnearwell/ml/utils.py
@@ -0,0 +1,24 @@
+"""Utilility functions for the ensemble and ML capabilities.
+
+Note: ``ml.ensemble`` makes use of ``np.random.default_rng``, which ignores the global
+seed of ``numpy``. Make sure to set them locally for full determinism.
+
+"""
+
+from typing import Optional
+
+import tensorflow as tf
+
+
+def enable_determinism(seed: Optional[int] = None):
+ """Set a seed for python, numpy, and tensorflow and enable deterministic behavior.
+
+ Args:
+ seed: (Optional[int]): Seed for the ``np.random.Generator``. Default is
+ ``None``.
+
+ """
+ # ``tf.keras.utils.set_random_seed`` sets the python, numpy, and tensorflow seed
+ # simultaneously.
+ tf.keras.utils.set_random_seed(seed=seed)
+ tf.config.experimental.enable_op_determinism()
diff --git a/src/pyopmnearwell/templates/co2store/base.mako b/src/pyopmnearwell/templates/co2store/base.mako
index 2b6e7e1..5db7da0 100644
--- a/src/pyopmnearwell/templates/co2store/base.mako
+++ b/src/pyopmnearwell/templates/co2store/base.mako
@@ -242,6 +242,14 @@ WCONPROD
%endif
/
%endif
+-- Close the specified connections
+% if len(dic['inj'][j]) >= 6:
+WELOPEN
+% for i in range(1, dic['noCells'][2] + 1):
+'INJ0' ${'SHUT' if i in dic['inj'][j][5:] else 'OPEN'} 0 0 ${i} /
+% endfor
+/
+% endif
TSTEP
${mt.floor(dic['inj'][j][0]/dic['inj'][j][1])}*${dic['inj'][j][1]}
/
diff --git a/src/pyopmnearwell/templates/co2store/no_disgas_no_diffusion.mako b/src/pyopmnearwell/templates/co2store/no_disgas_no_diffusion.mako
index d43c8b0..32409e0 100644
--- a/src/pyopmnearwell/templates/co2store/no_disgas_no_diffusion.mako
+++ b/src/pyopmnearwell/templates/co2store/no_disgas_no_diffusion.mako
@@ -192,8 +192,7 @@ COMPDAT
% if dic['grid'] == 'core':
'INJ0' 1 ${1+mt.floor(dic['noCells'][2]/2)} ${1+mt.floor(dic['noCells'][2]/2)} ${1+mt.floor(dic['noCells'][2]/2)} 'OPEN' 1* ${dic["jfactor"]} /
% else:
---'INJ0' ${max(1, 1+mt.floor(dic['noCells'][1]/2))} ${max(1, 1+mt.floor(dic['noCells'][1]/2))} 1 ${dic['noCells'][2]} 'OPEN' 1* ${dic["jfactor"]*2*mt.pi*dic['rock'][0][0]*dic['dims'][2]/dic['noCells'][2]} /
-'INJ0' ${max(1, 1+mt.floor(dic['noCells'][1]/2))} ${max(1, 1+mt.floor(dic['noCells'][1]/2))} 1 ${dic['noCells'][2]} 'OPEN' 1* ${dic["jfactor"]} /
+'INJ0' ${max(1, 1+mt.floor(dic['noCells'][1]/2))} ${max(1, 1+mt.floor(dic['noCells'][1]/2))} 1 ${dic['noCells'][2]} 'OPEN' 1* ${dic["jfactor"]} ${dic['diameter']}/
% endif
% endif
@@ -222,7 +221,7 @@ WCONINJE
% else:
'INJ0' 'WATER' ${'OPEN' if dic['inj'][j][4] > 0 else 'SHUT'}
'RATE' ${f"{dic['inj'][j][4] / 998.108 : E}"} 1* 400/
-%endif
+% endif
/
% if dic["pvMult"] == 0 or dic['grid']== 'core':
WCONPROD
@@ -233,9 +232,17 @@ WCONPROD
'PRO3' ${'OPEN' if dic['inj'][j][4] > 0 else 'SHUT'} 'BHP' 5* ${dic['pressure']}/
% else:
'PRO0' ${'OPEN' if dic['inj'][j][4] > 0 else 'SHUT'} 'BHP' 5* ${dic['pressure']}/
-%endif
+% endif
/
-%endif
+% endif
+-- Close the specified connections
+% if len(dic['inj'][j]) >= 6:
+WELOPEN
+% for i in range(1, dic['noCells'][2] + 1):
+'INJ0' ${'SHUT' if i in dic['inj'][j][5:] else 'OPEN'} 0 0 ${i} /
+% endfor
+/
+% endif
TSTEP
${mt.floor(dic['inj'][j][0]/dic['inj'][j][1])}*${dic['inj'][j][1]}
/
diff --git a/src/pyopmnearwell/templates/standardwell/base.hpp b/src/pyopmnearwell/templates/standardwell/base.hpp
deleted file mode 100644
index c15e536..0000000
--- a/src/pyopmnearwell/templates/standardwell/base.hpp
+++ /dev/null
@@ -1,513 +0,0 @@
-/*
- Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
- Copyright 2017 Statoil ASA.
- Copyright 2016 - 2017 IRIS AS.
-
- This file is part of the Open Porous Media project (OPM).
-
- OPM is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OPM is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with OPM. If not, see .
-*/
-
-
-#ifndef OPM_STANDARDWELL_HEADER_INCLUDED
-#define OPM_STANDARDWELL_HEADER_INCLUDED
-
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-
-#include
-#include
-#include
-#include
-#include
-#include
-
-#include
-#include
-
-#include
-
-#include
-#include
-
-#include
-#include
-
-namespace Opm
-{
-
- template
- class StandardWell : public WellInterface
- , public StandardWellEval,
- GetPropType,
- GetPropType>
- {
-
- public:
- using Base = WellInterface;
- using StdWellEval = StandardWellEval,
- GetPropType,
- GetPropType>;
-
- // TODO: some functions working with AD variables handles only with values (double) without
- // dealing with derivatives. It can be beneficial to make functions can work with either AD or scalar value.
- // And also, it can also be beneficial to make these functions hanle different types of AD variables.
- using typename Base::Simulator;
- using typename Base::IntensiveQuantities;
- using typename Base::FluidSystem;
- using typename Base::MaterialLaw;
- using typename Base::ModelParameters;
- using typename Base::Indices;
- using typename Base::RateConverterType;
- using typename Base::SparseMatrixAdapter;
- using typename Base::FluidState;
- using typename Base::RateVector;
-
- using Base::has_solvent;
- using Base::has_zFraction;
- using Base::has_polymer;
- using Base::has_polymermw;
- using Base::has_foam;
- using Base::has_brine;
- using Base::has_energy;
- using Base::has_micp;
-
- using PolymerModule = BlackOilPolymerModule;
- using FoamModule = BlackOilFoamModule;
- using BrineModule = BlackOilBrineModule;
- using typename Base::PressureMatrix;
-
- // number of the conservation equations
- static constexpr int numWellConservationEq = Indices::numPhases + Indices::numSolvents;
- // number of the well control equations
- static constexpr int numWellControlEq = 1;
- // number of the well equations that will always be used
- // based on the solution strategy, there might be other well equations be introduced
- static constexpr int numStaticWellEq = numWellConservationEq + numWellControlEq;
-
- // the index for Bhp in primary variables and also the index of well control equation
- // they both will be the last one in their respective system.
- // TODO: we should have indices for the well equations and well primary variables separately
- static constexpr int Bhp = numStaticWellEq - numWellControlEq;
-
- using StdWellEval::WQTotal;
-
- using typename Base::Scalar;
-
-
- using Base::name;
- using Base::Water;
- using Base::Oil;
- using Base::Gas;
-
- using typename Base::BVector;
-
- using Eval = typename StdWellEval::Eval;
- using EvalWell = typename StdWellEval::EvalWell;
- using BVectorWell = typename StdWellEval::BVectorWell;
-
- StandardWell(const Well& well,
- const ParallelWellInfo& pw_info,
- const int time_step,
- const ModelParameters& param,
- const RateConverterType& rate_converter,
- const int pvtRegionIdx,
- const int num_components,
- const int num_phases,
- const int index_of_well,
- const std::vector& perf_data);
-
- virtual void init(const PhaseUsage* phase_usage_arg,
- const std::vector& depth_arg,
- const double gravity_arg,
- const int num_cells,
- const std::vector< Scalar >& B_avg,
- const bool changed_to_open_this_step) override;
-
-
- void initPrimaryVariablesEvaluation() override;
-
- /// check whether the well equations get converged for this well
- virtual ConvergenceReport getWellConvergence(const SummaryState& summary_state,
- const WellState& well_state,
- const std::vector& B_avg,
- DeferredLogger& deferred_logger,
- const bool relax_tolerance) const override;
-
- /// Ax = Ax - C D^-1 B x
- virtual void apply(const BVector& x, BVector& Ax) const override;
- /// r = r - C D^-1 Rw
- virtual void apply(BVector& r) const override;
-
- /// using the solution x to recover the solution xw for wells and applying
- /// xw to update Well State
- void recoverWellSolutionAndUpdateWellState(const SummaryState& summary_state,
- const BVector& x,
- WellState& well_state,
- DeferredLogger& deferred_logger) override;
-
- /// computing the well potentials for group control
- virtual void computeWellPotentials(const Simulator& ebosSimulator,
- const WellState& well_state,
- std::vector& well_potentials,
- DeferredLogger& deferred_logger) /* const */ override;
-
- void updatePrimaryVariables(const SummaryState& summary_state,
- const WellState& well_state,
- DeferredLogger& deferred_logger) override;
-
- virtual void solveEqAndUpdateWellState(const SummaryState& summary_state,
- WellState& well_state,
- DeferredLogger& deferred_logger) override;
-
- virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
- const WellState& well_state,
- DeferredLogger& deferred_logger) override; // should be const?
-
- virtual void updateProductivityIndex(const Simulator& ebosSimulator,
- const WellProdIndexCalculator& wellPICalc,
- WellState& well_state,
- DeferredLogger& deferred_logger) const override;
-
- virtual double connectionDensity(const int globalConnIdx,
- const int openConnIdx) const override;
-
- virtual void addWellContributions(SparseMatrixAdapter& mat) const override;
-
- virtual void addWellPressureEquations(PressureMatrix& mat,
- const BVector& x,
- const int pressureVarIndex,
- const bool use_well_weights,
- const WellState& well_state) const override;
-
- // iterate well equations with the specified control until converged
- bool iterateWellEqWithControl(const Simulator& ebosSimulator,
- const double dt,
- const Well::InjectionControls& inj_controls,
- const Well::ProductionControls& prod_controls,
- WellState& well_state,
- const GroupState& group_state,
- DeferredLogger& deferred_logger) override;
-
- /// \brief Wether the Jacobian will also have well contributions in it.
- virtual bool jacobianContainsWellContributions() const override
- {
- return this->param_.matrix_add_well_contributions_;
- }
-
- /* returns BHP */
- double computeWellRatesAndBhpWithThpAlqProd(const Simulator &ebos_simulator,
- const SummaryState &summary_state,
- DeferredLogger &deferred_logger,
- std::vector &potentials,
- double alq) const;
-
- void computeWellRatesWithThpAlqProd(
- const Simulator &ebos_simulator,
- const SummaryState &summary_state,
- DeferredLogger &deferred_logger,
- std::vector &potentials,
- double alq) const;
-
- std::optional computeBhpAtThpLimitProdWithAlq(
- const Simulator& ebos_simulator,
- const SummaryState& summary_state,
- const double alq_value,
- DeferredLogger& deferred_logger) const override;
-
- virtual void computeWellRatesWithBhp(
- const Simulator& ebosSimulator,
- const double& bhp,
- std::vector& well_flux,
- DeferredLogger& deferred_logger) const override;
-
- // NOTE: These cannot be protected since they are used by GasLiftRuntime
- using Base::phaseUsage;
- using Base::vfp_properties_;
-
- virtual std::vector computeCurrentWellRates(const Simulator& ebosSimulator,
- DeferredLogger& deferred_logger) const override;
-
- std::vector getPrimaryVars() const override;
-
- int setPrimaryVars(std::vector::const_iterator it) override;
-
- protected:
- bool regularize_;
-
- // updating the well_state based on well solution dwells
- void updateWellState(const SummaryState& summary_state,
- const BVectorWell& dwells,
- WellState& well_state,
- DeferredLogger& deferred_logger);
-
- // calculate the properties for the well connections
- // to calulate the pressure difference between well connections.
- using WellConnectionProps = typename StdWellEval::StdWellConnections::Properties;
- void computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator,
- const WellState& well_state,
- WellConnectionProps& props) const;
-
- void computeWellConnectionDensitesPressures(const Simulator& ebosSimulator,
- const WellState& well_state,
- const WellConnectionProps& props,
- DeferredLogger& deferred_logger);
-
- void computeWellConnectionPressures(const Simulator& ebosSimulator,
- const WellState& well_state,
- DeferredLogger& deferred_logger);
-
- template
- void computePerfRate(const IntensiveQuantities& intQuants,
- const std::vector& mob,
- const Value& bhp,
- const double Tw,
- const int perf,
- const bool allow_cf,
- const double elapsed_time,
- std::vector& cq_s,
- PerforationRates& perf_rates,
- DeferredLogger& deferred_logger) const;
-
- template
- void computePerfRate(const std::vector& mob,
- const Value& pressure,
- const Value& bhp,
- const Value& rs,
- const Value& rv,
- const Value& rvw,
- const Value& rsw,
- std::vector& b_perfcells_dense,
- const double Tw,
- const int perf,
- const bool allow_cf,
- const Value& skin_pressure,
- const std::vector& cmix_s,
- const double elapsed_time,
- std::vector& cq_s,
- PerforationRates& perf_rates,
- DeferredLogger& deferred_logger) const;
-
- void computeWellRatesWithBhpIterations(const Simulator& ebosSimulator,
- const double& bhp,
- std::vector& well_flux,
- DeferredLogger& deferred_logger) const override;
-
- std::vector computeWellPotentialWithTHP(
- const Simulator& ebosSimulator,
- DeferredLogger& deferred_logger,
- const WellState &well_state) const;
-
- virtual double getRefDensity() const override;
-
- // get the mobility for specific perforation
- template
- void getMobility(const Simulator& ebosSimulator,
- const int perf,
- std::vector& mob,
- DeferredLogger& deferred_logger) const;
-
- void updateWaterMobilityWithPolymer(const Simulator& ebos_simulator,
- const int perf,
- std::vector& mob_water,
- DeferredLogger& deferred_logger) const;
-
- void updatePrimaryVariablesNewton(const BVectorWell& dwells,
- const bool stop_or_zero_rate_target,
- DeferredLogger& deferred_logger);
-
- void updateWellStateFromPrimaryVariables(const bool stop_or_zero_rate_target,
- WellState& well_state,
- const SummaryState& summary_state,
- DeferredLogger& deferred_logger) const;
-
- virtual void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
- const double dt,
- const Well::InjectionControls& inj_controls,
- const Well::ProductionControls& prod_controls,
- WellState& well_state,
- const GroupState& group_state,
- DeferredLogger& deferred_logger) override;
-
- void assembleWellEqWithoutIterationImpl(const Simulator& ebosSimulator,
- const double dt,
- const Well::InjectionControls& inj_controls,
- const Well::ProductionControls& prod_controls,
- WellState& well_state,
- const GroupState& group_state,
- DeferredLogger& deferred_logger);
-
- void calculateSinglePerf(const Simulator& ebosSimulator,
- const int perf,
- WellState& well_state,
- std::vector& connectionRates,
- std::vector& cq_s,
- EvalWell& water_flux_s,
- EvalWell& cq_s_zfrac_effective,
- DeferredLogger& deferred_logger) const;
-
- // check whether the well is operable under BHP limit with current reservoir condition
- virtual void checkOperabilityUnderBHPLimit(const WellState& well_state, const Simulator& ebos_simulator, DeferredLogger& deferred_logger) override;
-
- // check whether the well is operable under THP limit with current reservoir condition
- virtual void checkOperabilityUnderTHPLimit(const Simulator& ebos_simulator, const WellState& well_state, DeferredLogger& deferred_logger) override;
-
- // updating the inflow based on the current reservoir condition
- virtual void updateIPR(const Simulator& ebos_simulator, DeferredLogger& deferred_logger) const override;
-
- // for a well, when all drawdown are in the wrong direction, then this well will not
- // be able to produce/inject .
- bool allDrawDownWrongDirection(const Simulator& ebos_simulator) const;
-
- // whether the well can produce / inject based on the current well state (bhp)
- bool canProduceInjectWithCurrentBhp(const Simulator& ebos_simulator,
- const WellState& well_state,
- DeferredLogger& deferred_logger);
-
- // turn on crossflow to avoid singular well equations
- // when the well is banned from cross-flow and the BHP is not properly initialized,
- // we turn on crossflow to avoid singular well equations. It can result in wrong-signed
- // well rates, it can cause problem for THP calculation
- // TODO: looking for better alternative to avoid wrong-signed well rates
- bool openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const;
-
- // calculate the skin pressure based on water velocity, throughput and polymer concentration.
- // throughput is used to describe the formation damage during water/polymer injection.
- // calculated skin pressure will be applied to the drawdown during perforation rate calculation
- // to handle the effect from formation damage.
- EvalWell pskin(const double throuhgput,
- const EvalWell& water_velocity,
- const EvalWell& poly_inj_conc,
- DeferredLogger& deferred_logger) const;
-
- // calculate the skin pressure based on water velocity, throughput during water injection.
- EvalWell pskinwater(const double throughput,
- const EvalWell& water_velocity,
- DeferredLogger& deferred_logger) const;
-
- // calculate the injecting polymer molecular weight based on the througput and water velocity
- EvalWell wpolymermw(const double throughput,
- const EvalWell& water_velocity,
- DeferredLogger& deferred_logger) const;
-
- // modify the water rate for polymer injectivity study
- void handleInjectivityRate(const Simulator& ebosSimulator,
- const int perf,
- std::vector& cq_s) const;
-
- // handle the extra equations for polymer injectivity study
- void handleInjectivityEquations(const Simulator& ebosSimulator,
- const WellState& well_state,
- const int perf,
- const EvalWell& water_flux_s,
- DeferredLogger& deferred_logger);
-
- virtual void updateWaterThroughput(const double dt, WellState& well_state) const override;
-
- // checking convergence of extra equations, if there are any
- void checkConvergenceExtraEqs(const std::vector& res,
- ConvergenceReport& report) const;
-
- // updating the connectionRates_ related polymer molecular weight
- void updateConnectionRatePolyMW(const EvalWell& cq_s_poly,
- const IntensiveQuantities& int_quants,
- const WellState& well_state,
- const int perf,
- std::vector& connectionRates,
- DeferredLogger& deferred_logger) const;
-
-
- std::optional computeBhpAtThpLimitProd(const WellState& well_state,
- const Simulator& ebos_simulator,
- const SummaryState& summary_state,
- DeferredLogger& deferred_logger) const;
-
- std::optional computeBhpAtThpLimitInj(const Simulator& ebos_simulator,
- const SummaryState& summary_state,
- DeferredLogger& deferred_logger) const;
- template
- Value wellIndexEval(const int perf, const double elapsed_time, const Value& pressure) const;
-
- private:
-
- template
- Value scaleFunction(Value X, double min, double max) const;
-
- template
- Value unscaleFunction(Value X, double min, double max) const;
-
- Eval connectionRateEnergy(const double maxOilSaturation,
- const std::vector& cq_s,
- const IntensiveQuantities& intQuants,
- DeferredLogger& deferred_logger) const;
-
- template
- void gasOilPerfRateInj(const std::vector& cq_s,
- PerforationRates& perf_rates,
- const Value& rv,
- const Value& rs,
- const Value& pressure,
- const Value& rvw,
- DeferredLogger& deferred_logger) const;
-
- template
- void gasOilPerfRateProd(std::vector& cq_s,
- PerforationRates& perf_rates,
- const Value& rv,
- const Value& rs,
- const Value& rvw) const;
-
- template
- void gasWaterPerfRateProd(std::vector& cq_s,
- PerforationRates& perf_rates,
- const Value& rvw,
- const Value& rsw) const;
-
- template
- void gasWaterPerfRateInj(const std::vector& cq_s,
- PerforationRates& perf_rates,
- const Value& rvw,
- const Value& rsw,
- const Value& pressure,
- DeferredLogger& deferred_logger) const;
-
- template
- void disOilVapWatVolumeRatio(Value& volumeRatio,
- const Value& rvw,
- const Value& rsw,
- const Value& pressure,
- const std::vector& cmix_s,
- const std::vector& b_perfcells_dense,
- DeferredLogger& deferred_logger) const;
-
- template
- void gasOilVolumeRatio(Value& volumeRatio,
- const Value& rv,
- const Value& rs,
- const Value& pressure,
- const std::vector& cmix_s,
- const std::vector& b_perfcells_dense,
- DeferredLogger& deferred_logger) const;
- };
-
-}
-
-#include "StandardWell_impl.hpp"
-
-#endif // OPM_STANDARDWELL_HEADER_INCLUDED
diff --git a/src/pyopmnearwell/templates/standardwell/local_stencil.hpp b/src/pyopmnearwell/templates/standardwell/local_stencil.hpp
deleted file mode 100644
index 4b8bbfc..0000000
--- a/src/pyopmnearwell/templates/standardwell/local_stencil.hpp
+++ /dev/null
@@ -1,515 +0,0 @@
-/*
- Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
- Copyright 2017 Statoil ASA.
- Copyright 2016 - 2017 IRIS AS.
-
- This file is part of the Open Porous Media project (OPM).
-
- OPM is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OPM is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with OPM. If not, see .
-*/
-
-
-#ifndef OPM_STANDARDWELL_HEADER_INCLUDED
-#define OPM_STANDARDWELL_HEADER_INCLUDED
-
-#include
-#include
-#include
-#include
-#include
-#include
-#include
-
-#include
-#include
-#include
-#include
-#include
-#include
-
-#include
-#include
-
-#include
-
-#include
-#include
-
-#include
-#include
-
-namespace Opm
-{
-
- template
- class StandardWell : public WellInterface
- , public StandardWellEval,
- GetPropType,
- GetPropType>
- {
-
- public:
- using Base = WellInterface;
- using StdWellEval = StandardWellEval,
- GetPropType,
- GetPropType>;
-
- // TODO: some functions working with AD variables handles only with values (double) without
- // dealing with derivatives. It can be beneficial to make functions can work with either AD or scalar value.
- // And also, it can also be beneficial to make these functions hanle different types of AD variables.
- using typename Base::Simulator;
- using typename Base::IntensiveQuantities;
- using typename Base::FluidSystem;
- using typename Base::MaterialLaw;
- using typename Base::ModelParameters;
- using typename Base::Indices;
- using typename Base::RateConverterType;
- using typename Base::SparseMatrixAdapter;
- using typename Base::FluidState;
- using typename Base::RateVector;
-
- using Base::has_solvent;
- using Base::has_zFraction;
- using Base::has_polymer;
- using Base::has_polymermw;
- using Base::has_foam;
- using Base::has_brine;
- using Base::has_energy;
- using Base::has_micp;
-
- using PolymerModule = BlackOilPolymerModule;
- using FoamModule = BlackOilFoamModule;
- using BrineModule = BlackOilBrineModule;
- using typename Base::PressureMatrix;
-
- // number of the conservation equations
- static constexpr int numWellConservationEq = Indices::numPhases + Indices::numSolvents;
- // number of the well control equations
- static constexpr int numWellControlEq = 1;
- // number of the well equations that will always be used
- // based on the solution strategy, there might be other well equations be introduced
- static constexpr int numStaticWellEq = numWellConservationEq + numWellControlEq;
-
- // the index for Bhp in primary variables and also the index of well control equation
- // they both will be the last one in their respective system.
- // TODO: we should have indices for the well equations and well primary variables separately
- static constexpr int Bhp = numStaticWellEq - numWellControlEq;
-
- using StdWellEval::WQTotal;
-
- using typename Base::Scalar;
-
-
- using Base::name;
- using Base::Water;
- using Base::Oil;
- using Base::Gas;
-
- using typename Base::BVector;
-
- using Eval = typename StdWellEval::Eval;
- using EvalWell = typename StdWellEval::EvalWell;
- using BVectorWell = typename StdWellEval::BVectorWell;
-
- StandardWell(const Well& well,
- const ParallelWellInfo& pw_info,
- const int time_step,
- const ModelParameters& param,
- const RateConverterType& rate_converter,
- const int pvtRegionIdx,
- const int num_components,
- const int num_phases,
- const int index_of_well,
- const std::vector& perf_data);
-
- virtual void init(const PhaseUsage* phase_usage_arg,
- const std::vector& depth_arg,
- const double gravity_arg,
- const int num_cells,
- const std::vector< Scalar >& B_avg,
- const bool changed_to_open_this_step) override;
-
-
- void initPrimaryVariablesEvaluation() override;
-
- /// check whether the well equations get converged for this well
- virtual ConvergenceReport getWellConvergence(const SummaryState& summary_state,
- const WellState& well_state,
- const std::vector& B_avg,
- DeferredLogger& deferred_logger,
- const bool relax_tolerance) const override;
-
- /// Ax = Ax - C D^-1 B x
- virtual void apply(const BVector& x, BVector& Ax) const override;
- /// r = r - C D^-1 Rw
- virtual void apply(BVector& r) const override;
-
- /// using the solution x to recover the solution xw for wells and applying
- /// xw to update Well State
- void recoverWellSolutionAndUpdateWellState(const SummaryState& summary_state,
- const BVector& x,
- WellState& well_state,
- DeferredLogger& deferred_logger) override;
-
- /// computing the well potentials for group control
- virtual void computeWellPotentials(const Simulator& ebosSimulator,
- const WellState& well_state,
- std::vector& well_potentials,
- DeferredLogger& deferred_logger) /* const */ override;
-
- void updatePrimaryVariables(const SummaryState& summary_state,
- const WellState& well_state,
- DeferredLogger& deferred_logger) override;
-
- virtual void solveEqAndUpdateWellState(const SummaryState& summary_state,
- WellState& well_state,
- DeferredLogger& deferred_logger) override;
-
- virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
- const WellState& well_state,
- DeferredLogger& deferred_logger) override; // should be const?
-
- virtual void updateProductivityIndex(const Simulator& ebosSimulator,
- const WellProdIndexCalculator& wellPICalc,
- WellState& well_state,
- DeferredLogger& deferred_logger) const override;
-
- virtual double connectionDensity(const int globalConnIdx,
- const int openConnIdx) const override;
-
- virtual void addWellContributions(SparseMatrixAdapter& mat) const override;
-
- virtual void addWellPressureEquations(PressureMatrix& mat,
- const BVector& x,
- const int pressureVarIndex,
- const bool use_well_weights,
- const WellState& well_state) const override;
-
- // iterate well equations with the specified control until converged
- bool iterateWellEqWithControl(const Simulator& ebosSimulator,
- const double dt,
- const Well::InjectionControls& inj_controls,
- const Well::ProductionControls& prod_controls,
- WellState& well_state,
- const GroupState& group_state,
- DeferredLogger& deferred_logger) override;
-
- /// \brief Wether the Jacobian will also have well contributions in it.
- virtual bool jacobianContainsWellContributions() const override
- {
- return this->param_.matrix_add_well_contributions_;
- }
-
- /* returns BHP */
- double computeWellRatesAndBhpWithThpAlqProd(const Simulator &ebos_simulator,
- const SummaryState &summary_state,
- DeferredLogger &deferred_logger,
- std::vector &potentials,
- double alq) const;
-
- void computeWellRatesWithThpAlqProd(
- const Simulator &ebos_simulator,
- const SummaryState &summary_state,
- DeferredLogger &deferred_logger,
- std::vector &potentials,
- double alq) const;
-
- std::optional computeBhpAtThpLimitProdWithAlq(
- const Simulator& ebos_simulator,
- const SummaryState& summary_state,
- const double alq_value,
- DeferredLogger& deferred_logger) const override;
-
- virtual void computeWellRatesWithBhp(
- const Simulator& ebosSimulator,
- const double& bhp,
- std::vector& well_flux,
- DeferredLogger& deferred_logger) const override;
-
- // NOTE: These cannot be protected since they are used by GasLiftRuntime
- using Base::phaseUsage;
- using Base::vfp_properties_;
-
- virtual std::vector computeCurrentWellRates(const Simulator& ebosSimulator,
- DeferredLogger& deferred_logger) const override;
-
- std::vector getPrimaryVars() const override;
-
- int setPrimaryVars(std::vector::const_iterator it) override;
-
- protected:
- bool regularize_;
-
- // updating the well_state based on well solution dwells
- void updateWellState(const SummaryState& summary_state,
- const BVectorWell& dwells,
- WellState& well_state,
- DeferredLogger& deferred_logger);
-
- // calculate the properties for the well connections
- // to calulate the pressure difference between well connections.
- using WellConnectionProps = typename StdWellEval::StdWellConnections::Properties;
- void computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator,
- const WellState& well_state,
- WellConnectionProps& props) const;
-
- void computeWellConnectionDensitesPressures(const Simulator& ebosSimulator,
- const WellState& well_state,
- const WellConnectionProps& props,
- DeferredLogger& deferred_logger);
-
- void computeWellConnectionPressures(const Simulator& ebosSimulator,
- const WellState& well_state,
- DeferredLogger& deferred_logger);
-
- template
- void computePerfRate(const IntensiveQuantities& intQuants,
- const std::vector& mob,
- const Value& bhp,
- const double Tw,
- const int perf,
- const bool allow_cf,
- const double elapsed_time,
- std::vector& cq_s,
- PerforationRates& perf_rates,
- DeferredLogger& deferred_logger,
- const Simulator& ebosSimulator) const;
-
- template
- void computePerfRate(const std::vector& mob,
- const Value& pressure,
- const Value& bhp,
- const Value& rs,
- const Value& rv,
- const Value& rvw,
- const Value& rsw,
- std::vector& b_perfcells_dense,
- const double Tw,
- const int perf,
- const bool allow_cf,
- const Value& skin_pressure,
- const std::vector& cmix_s,
- const double elapsed_time,
- std::vector& cq_s,
- PerforationRates& perf_rates,
- DeferredLogger& deferred_logger,
- const Simulator& ebosSimulator) const;
-
- void computeWellRatesWithBhpIterations(const Simulator& ebosSimulator,
- const double& bhp,
- std::vector& well_flux,
- DeferredLogger& deferred_logger) const override;
-
- std::vector computeWellPotentialWithTHP(
- const Simulator& ebosSimulator,
- DeferredLogger& deferred_logger,
- const WellState &well_state) const;
-
- virtual double getRefDensity() const override;
-
- // get the mobility for specific perforation
- template
- void getMobility(const Simulator& ebosSimulator,
- const int perf,
- std::vector& mob,
- DeferredLogger& deferred_logger) const;
-
- void updateWaterMobilityWithPolymer(const Simulator& ebos_simulator,
- const int perf,
- std::vector& mob_water,
- DeferredLogger& deferred_logger) const;
-
- void updatePrimaryVariablesNewton(const BVectorWell& dwells,
- const bool stop_or_zero_rate_target,
- DeferredLogger& deferred_logger);
-
- void updateWellStateFromPrimaryVariables(const bool stop_or_zero_rate_target,
- WellState& well_state,
- const SummaryState& summary_state,
- DeferredLogger& deferred_logger) const;
-
- virtual void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
- const double dt,
- const Well::InjectionControls& inj_controls,
- const Well::ProductionControls& prod_controls,
- WellState& well_state,
- const GroupState& group_state,
- DeferredLogger& deferred_logger) override;
-
- void assembleWellEqWithoutIterationImpl(const Simulator& ebosSimulator,
- const double dt,
- const Well::InjectionControls& inj_controls,
- const Well::ProductionControls& prod_controls,
- WellState& well_state,
- const GroupState& group_state,
- DeferredLogger& deferred_logger);
-
- void calculateSinglePerf(const Simulator& ebosSimulator,
- const int perf,
- WellState& well_state,
- std::vector& connectionRates,
- std::vector& cq_s,
- EvalWell& water_flux_s,
- EvalWell& cq_s_zfrac_effective,
- DeferredLogger& deferred_logger) const;
-
- // check whether the well is operable under BHP limit with current reservoir condition
- virtual void checkOperabilityUnderBHPLimit(const WellState& well_state, const Simulator& ebos_simulator, DeferredLogger& deferred_logger) override;
-
- // check whether the well is operable under THP limit with current reservoir condition
- virtual void checkOperabilityUnderTHPLimit(const Simulator& ebos_simulator, const WellState& well_state, DeferredLogger& deferred_logger) override;
-
- // updating the inflow based on the current reservoir condition
- virtual void updateIPR(const Simulator& ebos_simulator, DeferredLogger& deferred_logger) const override;
-
- // for a well, when all drawdown are in the wrong direction, then this well will not
- // be able to produce/inject .
- bool allDrawDownWrongDirection(const Simulator& ebos_simulator) const;
-
- // whether the well can produce / inject based on the current well state (bhp)
- bool canProduceInjectWithCurrentBhp(const Simulator& ebos_simulator,
- const WellState& well_state,
- DeferredLogger& deferred_logger);
-
- // turn on crossflow to avoid singular well equations
- // when the well is banned from cross-flow and the BHP is not properly initialized,
- // we turn on crossflow to avoid singular well equations. It can result in wrong-signed
- // well rates, it can cause problem for THP calculation
- // TODO: looking for better alternative to avoid wrong-signed well rates
- bool openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const;
-
- // calculate the skin pressure based on water velocity, throughput and polymer concentration.
- // throughput is used to describe the formation damage during water/polymer injection.
- // calculated skin pressure will be applied to the drawdown during perforation rate calculation
- // to handle the effect from formation damage.
- EvalWell pskin(const double throuhgput,
- const EvalWell& water_velocity,
- const EvalWell& poly_inj_conc,
- DeferredLogger& deferred_logger) const;
-
- // calculate the skin pressure based on water velocity, throughput during water injection.
- EvalWell pskinwater(const double throughput,
- const EvalWell& water_velocity,
- DeferredLogger& deferred_logger) const;
-
- // calculate the injecting polymer molecular weight based on the througput and water velocity
- EvalWell wpolymermw(const double throughput,
- const EvalWell& water_velocity,
- DeferredLogger& deferred_logger) const;
-
- // modify the water rate for polymer injectivity study
- void handleInjectivityRate(const Simulator& ebosSimulator,
- const int perf,
- std::vector& cq_s) const;
-
- // handle the extra equations for polymer injectivity study
- void handleInjectivityEquations(const Simulator& ebosSimulator,
- const WellState& well_state,
- const int perf,
- const EvalWell& water_flux_s,
- DeferredLogger& deferred_logger);
-
- virtual void updateWaterThroughput(const double dt, WellState& well_state) const override;
-
- // checking convergence of extra equations, if there are any
- void checkConvergenceExtraEqs(const std::vector& res,
- ConvergenceReport& report) const;
-
- // updating the connectionRates_ related polymer molecular weight
- void updateConnectionRatePolyMW(const EvalWell& cq_s_poly,
- const IntensiveQuantities& int_quants,
- const WellState& well_state,
- const int perf,
- std::vector& connectionRates,
- DeferredLogger& deferred_logger) const;
-
-
- std::optional computeBhpAtThpLimitProd(const WellState& well_state,
- const Simulator& ebos_simulator,
- const SummaryState& summary_state,
- DeferredLogger& deferred_logger) const;
-
- std::optional computeBhpAtThpLimitInj(const Simulator& ebos_simulator,
- const SummaryState& summary_state,
- DeferredLogger& deferred_logger) const;
- template
- Value wellIndexEval(const int perf, const Simulator& ebosSimulator) const;
-
- private:
-
- template
- Value scaleFunction(Value X, double min, double max, double range_min, double range_max) const;
-
- template
- Value unscaleFunction(Value X, double min, double max, double range_min, double range_max) const;
-
- Eval connectionRateEnergy(const double maxOilSaturation,
- const std::vector& cq_s,
- const IntensiveQuantities& intQuants,
- DeferredLogger& deferred_logger) const;
-
- template
- void gasOilPerfRateInj(const std::vector& cq_s,
- PerforationRates& perf_rates,
- const Value& rv,
- const Value& rs,
- const Value& pressure,
- const Value& rvw,
- DeferredLogger& deferred_logger) const;
-
- template
- void gasOilPerfRateProd(std::vector& cq_s,
- PerforationRates& perf_rates,
- const Value& rv,
- const Value& rs,
- const Value& rvw) const;
-
- template
- void gasWaterPerfRateProd(std::vector& cq_s,
- PerforationRates& perf_rates,
- const Value& rvw,
- const Value& rsw) const;
-
- template
- void gasWaterPerfRateInj(const std::vector& cq_s,
- PerforationRates& perf_rates,
- const Value& rvw,
- const Value& rsw,
- const Value& pressure,
- DeferredLogger& deferred_logger) const;
-
- template
- void disOilVapWatVolumeRatio(Value& volumeRatio,
- const Value& rvw,
- const Value& rsw,
- const Value& pressure,
- const std::vector& cmix_s,
- const std::vector& b_perfcells_dense,
- DeferredLogger& deferred_logger) const;
-
- template
- void gasOilVolumeRatio(Value& volumeRatio,
- const Value& rv,
- const Value& rs,
- const Value& pressure,
- const std::vector& cmix_s,
- const std::vector& b_perfcells_dense,
- DeferredLogger& deferred_logger) const;
- };
-
-}
-
-#include "StandardWell_impl.hpp"
-
-#endif // OPM_STANDARDWELL_HEADER_INCLUDED
diff --git a/src/pyopmnearwell/templates/standardwell_impl/co2_3_inputs.mako b/src/pyopmnearwell/templates/standardwell_impl/co2_3_inputs.mako
deleted file mode 100644
index 8bd3b0a..0000000
--- a/src/pyopmnearwell/templates/standardwell_impl/co2_3_inputs.mako
+++ /dev/null
@@ -1,2583 +0,0 @@
-/*
- Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
- Copyright 2017 Statoil ASA.
- Copyright 2016 - 2017 IRIS AS.
-
- This file is part of the Open Porous Media project (OPM).
-
- OPM is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- OPM is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with OPM. If not, see .
-*/
-
-#include
-
-#include
-
-#include
-
-#include
-#include
-#include
-#include
-#include
-
-#include
-
-#include
-
-#include
-#include
-#include
-
-
-namespace {
-
-template
-auto dValueError(const dValue& d,
- const std::string& name,
- const std::string& methodName,
- const Value& Rs,
- const Value& Rv,
- const Value& pressure)
-{
- return fmt::format("Problematic d value {} obtained for well {}"
- " during {} calculations with rs {}"
- ", rv {} and pressure {}."
- " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
- " for this connection.", d, name, methodName, Rs, Rv, pressure);
-}
-
-}
-
-namespace Opm
-{
-
- template
- StandardWell::
- StandardWell(const Well& well,
- const ParallelWellInfo& pw_info,
- const int time_step,
- const ModelParameters& param,
- const RateConverterType& rate_converter,
- const int pvtRegionIdx,
- const int num_components,
- const int num_phases,
- const int index_of_well,
- const std::vector& perf_data)
- : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
- , StdWellEval(static_cast&>(*this))
- , regularize_(false)
- {
- assert(this->num_components_ == numWellConservationEq);
- }
-
-
-
-
-
- template
- void
- StandardWell::
- init(const PhaseUsage* phase_usage_arg,
- const std::vector& depth_arg,
- const double gravity_arg,
- const int num_cells,
- const std::vector< Scalar >& B_avg,
- const bool changed_to_open_this_step)
- {
- Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg, changed_to_open_this_step);
- this->StdWellEval::init(this->perf_depth_, depth_arg, num_cells, Base::has_polymermw);
- }
-
-
-
-
-
- template
- void StandardWell::
- initPrimaryVariablesEvaluation()
- {
- this->primary_variables_.init();
- }
-
-
-
-
-
- template
- template
- void
- StandardWell::
- computePerfRate(const IntensiveQuantities& intQuants,
- const std::vector& mob,
- const Value& bhp,
- const double Tw,
- const int perf,
- const bool allow_cf,
- const double elapsed_time,
- std::vector& cq_s,
- PerforationRates& perf_rates,
- DeferredLogger& deferred_logger) const
- {
- auto obtain = [this](const Eval& value)
- {
- if constexpr (std::is_same_v) {
- static_cast(this); // suppress clang warning
- return getValue(value);
- } else {
- return this->extendEval(value);
- }
- };
- auto obtainN = [](const auto& value)
- {
- if constexpr (std::is_same_v) {
- return getValue(value);
- } else {
- return value;
- }
- };
- auto zeroElem = [this]()
- {
- if constexpr (std::is_same_v) {
- static_cast(this); // suppress clang warning
- return 0.0;
- } else {
- return Value{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
- }
- };
-
- const auto& fs = intQuants.fluidState();
- const Value pressure = obtain(this->getPerfCellPressure(fs));
- const Value rs = obtain(fs.Rs());
- const Value rv = obtain(fs.Rv());
- const Value rvw = obtain(fs.Rvw());
- const Value rsw = obtain(fs.Rsw());
-
- std::vector b_perfcells_dense(this->numComponents(), zeroElem());
- for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
- if (!FluidSystem::phaseIsActive(phaseIdx)) {
- continue;
- }
- const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
- b_perfcells_dense[compIdx] = obtain(fs.invB(phaseIdx));
- }
- if constexpr (has_solvent) {
- b_perfcells_dense[Indices::contiSolventEqIdx] = obtain(intQuants.solventInverseFormationVolumeFactor());
- }
-
- if constexpr (has_zFraction) {
- if (this->isInjector()) {
- const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
- b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
- b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
- }
- }
-
- Value skin_pressure = zeroElem();
- if (has_polymermw) {
- if (this->isInjector()) {
- const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
- skin_pressure = obtainN(this->primary_variables_.eval(pskin_index));
- }
- }
-
- // surface volume fraction of fluids within wellbore
- std::vector cmix_s(this->numComponents(), zeroElem());
- for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
- cmix_s[componentIdx] = obtainN(this->primary_variables_.surfaceVolumeFraction(componentIdx));
- }
-
- computePerfRate(mob,
- pressure,
- bhp,
- rs,
- rv,
- rvw,
- rsw,
- b_perfcells_dense,
- Tw,
- perf,
- allow_cf,
- skin_pressure,
- cmix_s,
- elapsed_time,
- cq_s,
- perf_rates,
- deferred_logger);
- }
-
- template
- template
- void
- StandardWell::
- computePerfRate(const std::vector& mob,
- const Value& pressure,
- const Value& bhp,
- const Value& rs,
- const Value& rv,
- const Value& rvw,
- const Value& rsw,
- std::vector& b_perfcells_dense,
- const double Tw,
- const int perf,
- const bool allow_cf,
- const Value& skin_pressure,
- const std::vector& cmix_s,
- const double elapsed_time,
- std::vector& cq_s,
- PerforationRates& perf_rates,
- DeferredLogger& deferred_logger) const
- {
- // Pressure drawdown (also used to determine direction of flow)
- const Value well_pressure = bhp + this->connections_.pressure_diff(perf);
- Value drawdown = pressure - well_pressure;
- if (this->isInjector()) {
- drawdown += skin_pressure;
- }
-
- // producing perforations
- if (drawdown > 0) {
- // Do nothing if crossflow is not allowed
- if (!allow_cf && this->isInjector()) {
- return;
- }
-
- // compute component volumetric rates at standard conditions
- for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
- const Value cq_p = - Tw * (mob[componentIdx] * drawdown);
- cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
- }
-
- if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
- gasOilPerfRateProd(cq_s, perf_rates, rv, rs, rvw);
- } else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
- gasWaterPerfRateProd(cq_s, perf_rates, rvw, rsw);
- }
- } else {
- // Do nothing if crossflow is not allowed
- if (!allow_cf && this->isProducer()) {
- return;
- }
-
- // Using total mobilities
- Value total_mob_dense = mob[0];
- for (int componentIdx = 1; componentIdx < this->numComponents(); ++componentIdx) {
- total_mob_dense += mob[componentIdx];
- }
-
- // Use well index from ML
- if (!Base::ml_wi_filename_.empty()) {
- Value WI;
- if (allow_cf) {
- OPM_DEFLOG_THROW(std::runtime_error,
- fmt::format("ML well index only implemented with no cross flow. Well {}", name()),
- deferred_logger);
- }
- if constexpr (std::is_same_v) {
- WI = this->extendEval(wellIndexEval(perf, elapsed_time, Base::restrictEval(pressure)));
- } else {
- WI = wellIndexEval(perf, elapsed_time, pressure);
- }
- auto injectorType = this->well_ecl_.injectorType();
- if (injectorType == InjectorType::WATER) {
- const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
- //cq_s[waterCompIdx] = - WI * (total_mob_dense * drawdown) * b_perfcells_dense[waterCompIdx];
- cq_s[waterCompIdx] = - WI * drawdown;
- std::cout << cq_s[waterCompIdx] << " " << - Tw * (total_mob_dense * drawdown)*b_perfcells_dense[waterCompIdx] << " " << cmix_s[waterCompIdx] << " " << b_perfcells_dense[waterCompIdx] << " " << total_mob_dense << std::endl;
- std::cout << WI << " " << (Tw*total_mob_dense*b_perfcells_dense[waterCompIdx]) << std::endl;
- }
- else if (injectorType == InjectorType::GAS) {
- const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
- //cq_s[gasCompIdx] = - WI * (total_mob_dense * drawdown) * b_perfcells_dense[gasCompIdx];
- std::cout << "together" << (total_mob_dense * drawdown) / b_perfcells_dense[gasCompIdx] << std::endl;
- std::cout << "total_mob_dense" << total_mob_dense << std::endl;
- std::cout << "drawdown" << drawdown << std::endl;
- std::cout << "b_perfcells_dense[gasCompIdx]" <<
- b_perfcells_dense[gasCompIdx] << std::endl;
-
- cq_s[gasCompIdx] = - WI * drawdown;
- }
- } else {
- std::cout << "ML model not found";
-
- const Value cqt_i = - Tw * (total_mob_dense * drawdown);
-
- // compute volume ratio between connection at standard conditions
- Value volumeRatio = bhp * 0.0; // initialize it with the correct type
-;
- if (FluidSystem::enableVaporizedWater() && FluidSystem::enableDissolvedGasInWater()) {
- disOilVapWatVolumeRatio(volumeRatio, rvw, rsw, pressure,
- cmix_s, b_perfcells_dense, deferred_logger);
- } else {
- if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
- const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
- volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
- }
- }
-
- if constexpr (Indices::enableSolvent) {
- volumeRatio += cmix_s[Indices::contiSolventEqIdx] / b_perfcells_dense[Indices::contiSolventEqIdx];
- }
-
- if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
- gasOilVolumeRatio(volumeRatio, rv, rs, pressure,
- cmix_s, b_perfcells_dense, deferred_logger);
- }
- else {
- if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
- const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
- volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
- }
- if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
- const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
- volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
- }
- }
-
- // injecting connections total volumerates at standard conditions
- Value cqt_is = cqt_i / volumeRatio;
- for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
- cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
- }
- }
-
- // calculating the perforation solution gas rate and solution oil rates
- if (this->isProducer()) {
- if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
- gasOilPerfRateInj(cq_s, perf_rates,
- rv, rs, pressure, rvw, deferred_logger);
- }
- if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
- //no oil
- gasWaterPerfRateInj(cq_s, perf_rates, rvw, rsw,
- pressure, deferred_logger);
- }
- }
- }
- }
-
-
- template
- void
- StandardWell::
- assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
- const double dt,
- const Well::InjectionControls& inj_controls,
- const Well::ProductionControls& prod_controls,
- WellState& well_state,
- const GroupState& group_state,
- DeferredLogger& deferred_logger)
- {
- // TODO: only_wells should be put back to save some computation
- // for example, the matrices B C does not need to update if only_wells
- if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
-
- // clear all entries
- this->linSys_.clear();
-
- assembleWellEqWithoutIterationImpl(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
- }
-
-
-
-
- template
- void
- StandardWell::
- assembleWellEqWithoutIterationImpl(const Simulator& ebosSimulator,
- const double dt,
- const Well::InjectionControls& inj_controls,
- const Well::ProductionControls& prod_controls,
- WellState& well_state,
- const GroupState& group_state,
- DeferredLogger& deferred_logger)
- {
- // try to regularize equation if the well does not converge
- const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
- const double volume = 0.1 * unit::cubic(unit::feet) * regularization_factor;
-
- auto& ws = well_state.well(this->index_of_well_);
- ws.phase_mixing_rates.fill(0.0);
-
- const int np = this->number_of_phases_;
-
- std::vector connectionRates = this->connectionRates_; // Copy to get right size.
- auto& perf_data = ws.perf_data;
- auto& perf_rates = perf_data.phase_rates;
- for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
- // Calculate perforation quantities.
- std::vector cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
- EvalWell water_flux_s{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
- EvalWell cq_s_zfrac_effective{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
- calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
-
- // Equation assembly for this perforation.
- if constexpr (has_polymer && Base::has_polymermw) {
- if (this->isInjector()) {
- handleInjectivityEquations(ebosSimulator, well_state, perf, water_flux_s, deferred_logger);
- }
- }
- const int cell_idx = this->well_cells_[perf];
- for (int componentIdx = 0; componentIdx < this->num_components_; ++componentIdx) {
- // the cq_s entering mass balance equations need to consider the efficiency factors.
- const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
-
- connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
-
- StandardWellAssemble(*this).
- assemblePerforationEq(cq_s_effective,
- componentIdx,
- cell_idx,
- this->primary_variables_.numWellEq(),
- this->linSys_);
-
- // Store the perforation phase flux for later usage.
- if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
- auto& perf_rate_solvent = perf_data.solvent_rates;
- perf_rate_solvent[perf] = cq_s[componentIdx].value();
- } else {
- perf_rates[perf*np + this->ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
- }
- }
-
- if constexpr (has_zFraction) {
- StandardWellAssemble(*this).
- assembleZFracEq(cq_s_zfrac_effective,
- cell_idx,
- this->primary_variables_.numWellEq(),
- this->linSys_);
- }
- }
- // Update the connection
- this->connectionRates_ = connectionRates;
-
- // Accumulate dissolved gas and vaporized oil flow rates across all
- // ranks sharing this well (this->index_of_well_).
- {
- const auto& comm = this->parallel_well_info_.communication();
- comm.sum(ws.phase_mixing_rates.data(), ws.phase_mixing_rates.size());
- }
-
- // accumulate resWell_ and duneD_ in parallel to get effects of all perforations (might be distributed)
- this->linSys_.sumDistributed(this->parallel_well_info_.communication());
-
- // add vol * dF/dt + Q to the well equations;
- for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
- // TODO: following the development in MSW, we need to convert the volume of the wellbore to be surface volume
- // since all the rates are under surface condition
- EvalWell resWell_loc(this->primary_variables_.numWellEq() + Indices::numEq, 0.0);
- if (FluidSystem::numActivePhases() > 1) {
- assert(dt > 0);
- resWell_loc += (this->primary_variables_.surfaceVolumeFraction(componentIdx) -
- this->F0_[componentIdx]) * volume / dt;
- }
- resWell_loc -= this->primary_variables_.getQs(componentIdx) * this->well_efficiency_factor_;
- StandardWellAssemble(*this).
- assembleSourceEq(resWell_loc,
- componentIdx,
- this->primary_variables_.numWellEq(),
- this->linSys_);
- }
-
- const auto& summaryState = ebosSimulator.vanguard().summaryState();
- const Schedule& schedule = ebosSimulator.vanguard().schedule();
- StandardWellAssemble(*this).
- assembleControlEq(well_state, group_state,
- schedule, summaryState,
- inj_controls, prod_controls,
- this->primary_variables_,
- this->connections_.rho(),
- this->linSys_,
- deferred_logger);
-
-
- // do the local inversion of D.
- try {
- this->linSys_.invert();
- } catch( ... ) {
- OPM_DEFLOG_THROW(NumericalProblem, "Error when inverting local well equations for well " + name(), deferred_logger);
- }
- }
-
-
-
-
- template
- void
- StandardWell::
- calculateSinglePerf(const Simulator& ebosSimulator,
- const int perf,
- WellState& well_state,
- std::vector& connectionRates,
- std::vector& cq_s,
- EvalWell& water_flux_s,
- EvalWell& cq_s_zfrac_effective,
- DeferredLogger& deferred_logger) const
- {
- const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
- const EvalWell& bhp = this->primary_variables_.eval(Bhp);
- const int cell_idx = this->well_cells_[perf];
- const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
- std::vector mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
- getMobility(ebosSimulator, perf, mob, deferred_logger);
-
- PerforationRates perf_rates;
- double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(intQuants, cell_idx);
- const double Tw = this->wellIndex(perf) * trans_mult;
- const double elapsed_time = ebosSimulator.time();
-
- computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf, elapsed_time,
- cq_s, perf_rates, deferred_logger);
-
- auto& ws = well_state.well(this->index_of_well_);
- auto& perf_data = ws.perf_data;
- if constexpr (has_polymer && Base::has_polymermw) {
- if (this->isInjector()) {
- // Store the original water flux computed from the reservoir quantities.
- // It will be required to assemble the injectivity equations.
- const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
- water_flux_s = cq_s[water_comp_idx];
- // Modify the water flux for the rest of this function to depend directly on the
- // local water velocity primary variable.
- handleInjectivityRate(ebosSimulator, perf, cq_s);
- }
- }
-
- // updating the solution gas rate and solution oil rate
- if (this->isProducer()) {
- ws.phase_mixing_rates[ws.dissolved_gas] += perf_rates.dis_gas;
- ws.phase_mixing_rates[ws.dissolved_gas_in_water] += perf_rates.dis_gas_in_water;
- ws.phase_mixing_rates[ws.vaporized_oil] += perf_rates.vap_oil;
- ws.phase_mixing_rates[ws.vaporized_water] += perf_rates.vap_wat;
- }
-
- if constexpr (has_energy) {
- connectionRates[perf][Indices::contiEnergyEqIdx] =
- connectionRateEnergy(ebosSimulator.problem().maxOilSaturation(cell_idx),
- cq_s, intQuants, deferred_logger);
- }
-
- if constexpr (has_polymer) {
- std::variant polymerConcentration;
- if (this->isInjector()) {
- polymerConcentration = this->wpolymer();
- } else {
- polymerConcentration = this->extendEval(intQuants.polymerConcentration() *
- intQuants.polymerViscosityCorrection());
- }
-
- [[maybe_unused]] EvalWell cq_s_poly;
- std::tie(connectionRates[perf][Indices::contiPolymerEqIdx],
- cq_s_poly) =
- this->connections_.connectionRatePolymer(perf_data.polymer_rates[perf],
- cq_s, polymerConcentration);
-
- if constexpr (Base::has_polymermw) {
- updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state,
- perf, connectionRates, deferred_logger);
- }
- }
-
- if constexpr (has_foam) {
- std::variant foamConcentration;
- if (this->isInjector()) {
- foamConcentration = this->wfoam();
- } else {
- foamConcentration = this->extendEval(intQuants.foamConcentration());
- }
- connectionRates[perf][Indices::contiFoamEqIdx] =
- this->connections_.connectionRateFoam(cq_s, foamConcentration,
- FoamModule::transportPhase(),
- deferred_logger);
- }
-
- if constexpr (has_zFraction) {
- std::variant> solventConcentration;
- if (this->isInjector()) {
- solventConcentration = this->wsolvent();
- } else {
- solventConcentration = std::array{this->extendEval(intQuants.xVolume()),
- this->extendEval(intQuants.yVolume())};
- }
- std::tie(connectionRates[perf][Indices::contiZfracEqIdx],
- cq_s_zfrac_effective) =
- this->connections_.connectionRatezFraction(perf_data.solvent_rates[perf],
- perf_rates.dis_gas, cq_s,
- solventConcentration);
- }
-
- if constexpr (has_brine) {
- std::variant saltConcentration;
- if (this->isInjector()) {
- saltConcentration = this->wsalt();
- } else {
- saltConcentration = this->extendEval(intQuants.fluidState().saltConcentration());
- }
-
- connectionRates[perf][Indices::contiBrineEqIdx] =
- this->connections_.connectionRateBrine(perf_data.brine_rates[perf],
- perf_rates.vap_wat, cq_s,
- saltConcentration);
- }
-
- if constexpr (has_micp) {
- std::variant microbialConcentration;
- std::variant oxygenConcentration;
- std::variant ureaConcentration;
- if (this->isInjector()) {
- microbialConcentration = this->wmicrobes();
- oxygenConcentration = this->woxygen();
- ureaConcentration = this->wurea();
- } else {
- microbialConcentration = this->extendEval(intQuants.microbialConcentration());
- oxygenConcentration = this->extendEval(intQuants.oxygenConcentration());
- ureaConcentration = this->extendEval(intQuants.ureaConcentration());
- }
- std::tie(connectionRates[perf][Indices::contiMicrobialEqIdx],
- connectionRates[perf][Indices::contiOxygenEqIdx],
- connectionRates[perf][Indices::contiUreaEqIdx]) =
- this->connections_.connectionRatesMICP(cq_s,
- microbialConcentration,
- oxygenConcentration,
- ureaConcentration);
- }
-
- // Store the perforation pressure for later usage.
- perf_data.pressure[perf] = ws.bhp + this->connections_.pressure_diff(perf);
- }
-
-
-
- template
- template
- void
- StandardWell::
- getMobility(const Simulator& ebosSimulator,
- const int perf,
- std::vector& mob,
- DeferredLogger& deferred_logger) const
- {
- auto obtain = [this](const Eval& value)
- {
- if constexpr (std::is_same_v) {
- static_cast(this); // suppress clang warning
- return getValue(value);
- } else {
- return this->extendEval(value);
- }
- };
- WellInterface::getMobility(ebosSimulator, perf, mob,
- obtain, deferred_logger);
-
- // modify the water mobility if polymer is present
- if constexpr (has_polymer) {
- if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
- OPM_DEFLOG_THROW(std::runtime_error, "Water is required when polymer is active", deferred_logger);
- }
-
- // for the cases related to polymer molecular weight, we assume fully mixing
- // as a result, the polymer and water share the same viscosity
- if constexpr (!Base::has_polymermw) {
- if constexpr (std::is_same_v) {
- std::vector mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
- for (size_t i = 0; i < mob.size(); ++i)
- mob_eval[i].setValue(mob[i]);
- updateWaterMobilityWithPolymer(ebosSimulator, perf, mob_eval, deferred_logger);
- for (size_t i = 0; i < mob.size(); ++i) {
- mob[i] = getValue(mob_eval[i]);
- }
- } else {
- updateWaterMobilityWithPolymer(ebosSimulator, perf, mob, deferred_logger);
- }
- }
- }
-
- // if the injecting well has WINJMULT setup, we update the mobility accordingly
- if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
- const double bhp = this->primary_variables_.value(Bhp);
- const double perf_press = bhp + this->connections_.pressure_diff(perf);
- const double multiplier = this->getInjMult(perf, bhp, perf_press);
- for (size_t i = 0; i < mob.size(); ++i) {
- mob[i] *= multiplier;
- }
- }
- }
-
-
- template
- void
- StandardWell::
- updateWellState(const SummaryState& summary_state,
- const BVectorWell& dwells,
- WellState& well_state,
- DeferredLogger& deferred_logger)
- {
- if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
-
- const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
- updatePrimaryVariablesNewton(dwells, stop_or_zero_rate_target, deferred_logger);
-
- updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, summary_state, deferred_logger);
- Base::calculateReservoirRates(well_state.well(this->index_of_well_));
- }
-
-
-
-
-
- template
- void
- StandardWell::
- updatePrimaryVariablesNewton(const BVectorWell& dwells,
- const bool stop_or_zero_rate_target,
- DeferredLogger& deferred_logger)
- {
- const double dFLimit = this->param_.dwell_fraction_max_;
- const double dBHPLimit = this->param_.dbhp_max_rel_;
- this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit);
-
- // for the water velocity and skin pressure
- if constexpr (Base::has_polymermw) {
- this->primary_variables_.updateNewtonPolyMW(dwells);
- }
-
- this->primary_variables_.checkFinite(deferred_logger);
- }
-
-
-
-
-
- template
- void
- StandardWell::
- updateWellStateFromPrimaryVariables(const bool stop_or_zero_rate_target,
- WellState& well_state,
- const SummaryState& summary_state,
- DeferredLogger& deferred_logger) const
- {
- this->StdWellEval::updateWellStateFromPrimaryVariables(stop_or_zero_rate_target, well_state, summary_state, deferred_logger);
-
- // other primary variables related to polymer injectivity study
- if constexpr (Base::has_polymermw) {
- this->primary_variables_.copyToWellStatePolyMW(well_state);
- }
- }
-
-
-
-
-
- template
- void
- StandardWell::
- updateIPR(const Simulator& ebos_simulator, DeferredLogger& deferred_logger) const
- {
- // TODO: not handling solvent related here for now
-
- // initialize all the values to be zero to begin with
- std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
- std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
-
- for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
- std::vector mob(this->num_components_, 0.0);
- getMobility(ebos_simulator, perf, mob, deferred_logger);
-
- const int cell_idx = this->well_cells_[perf];
- const auto& int_quantities = ebos_simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
- const auto& fs = int_quantities.fluidState();
- // the pressure of the reservoir grid block the well connection is in
- double p_r = this->getPerfCellPressure(fs).value();
-
- // calculating the b for the connection
- std::vector b_perf(this->num_components_);
- for (size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
- if (!FluidSystem::phaseIsActive(phase)) {
- continue;
- }
- const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
- b_perf[comp_idx] = fs.invB(phase).value();
- }
- if constexpr (has_solvent) {
- b_perf[Indices::contiSolventEqIdx] = int_quantities.solventInverseFormationVolumeFactor().value();
- }
-
- // the pressure difference between the connection and BHP
- const double h_perf = this->connections_.pressure_diff(perf);
- const double pressure_diff = p_r - h_perf;
-
- // Let us add a check, since the pressure is calculated based on zero value BHP
- // it should not be negative anyway. If it is negative, we might need to re-formulate
- // to taking into consideration the crossflow here.
- if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
- deferred_logger.debug("CROSSFLOW_IPR",
- "cross flow found when updateIPR for well " + name()
- + " . The connection is ignored in IPR calculations");
- // we ignore these connections for now
- continue;
- }
-
- // the well index associated with the connection
- const double tw_perf = this->wellIndex(perf)*ebos_simulator.problem().template rockCompTransMultiplier(int_quantities, cell_idx);
-
- std::vector ipr_a_perf(this->ipr_a_.size());
- std::vector ipr_b_perf(this->ipr_b_.size());
- for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
- const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx];
- ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
- ipr_b_perf[comp_idx] += tw_mob;
- }
-
- // we need to handle the rs and rv when both oil and gas are present
- if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
- const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
- const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
- const double rs = (fs.Rs()).value();
- const double rv = (fs.Rv()).value();
-
- const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
- const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
-
- ipr_a_perf[gas_comp_idx] += dis_gas_a;
- ipr_a_perf[oil_comp_idx] += vap_oil_a;
-
- const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
- const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
-
- ipr_b_perf[gas_comp_idx] += dis_gas_b;
- ipr_b_perf[oil_comp_idx] += vap_oil_b;
- }
-
- for (size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
- this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
- this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
- }
- }
- this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
- this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
- }
-
-
- template
- void
- StandardWell::
- checkOperabilityUnderBHPLimit(const WellState& well_state, const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
- {
- const auto& summaryState = ebos_simulator.vanguard().summaryState();
- const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
- // Crude but works: default is one atmosphere.
- // TODO: a better way to detect whether the BHP is defaulted or not
- const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
- if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
- // if the BHP limit is not defaulted or the well does not have a THP limit
- // we need to check the BHP limit
- double total_ipr_mass_rate = 0.0;
- for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
- {
- if (!FluidSystem::phaseIsActive(phaseIdx)) {
- continue;
- }
-
- const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
- const double ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
-
- const double rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
- total_ipr_mass_rate += ipr_rate * rho;
- }
- if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
- this->operability_status_.operable_under_only_bhp_limit = false;
- }
-
- // checking whether running under BHP limit will violate THP limit
- if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
- // option 1: calculate well rates based on the BHP limit.
- // option 2: stick with the above IPR curve
- // we use IPR here
- std::vector well_rates_bhp_limit;
- computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
-
- this->adaptRatesForVFP(well_rates_bhp_limit);
- const double thp_limit = this->getTHPConstraint(summaryState);
- const double thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
- bhp_limit,
- this->connections_.rho(),
- this->getALQ(well_state),
- thp_limit,
- deferred_logger);
- if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
- this->operability_status_.obey_thp_limit_under_bhp_limit = false;
- }
- }
- } else {
- // defaulted BHP and there is a THP constraint
- // default BHP limit is about 1 atm.
- // when applied the hydrostatic pressure correction dp,
- // most likely we get a negative value (bhp + dp)to search in the VFP table,
- // which is not desirable.
- // we assume we can operate under defaulted BHP limit and will violate the THP limit
- // when operating under defaulted BHP limit.
- this->operability_status_.operable_under_only_bhp_limit = true;
- this->operability_status_.obey_thp_limit_under_bhp_limit = false;
- }
- }
-
-
-
-
-
- template
- void
- StandardWell