diff --git a/Example/3WELL.mako b/Example/3WELL/3WELL.mako similarity index 100% rename from Example/3WELL.mako rename to Example/3WELL/3WELL.mako diff --git a/Example/README.ipynb b/Example/3WELL/README.ipynb similarity index 100% rename from Example/README.ipynb rename to Example/3WELL/README.ipynb diff --git a/Example/README.md b/Example/3WELL/README.md similarity index 100% rename from Example/README.md rename to Example/3WELL/README.md diff --git a/Example/README_11_0.png b/Example/3WELL/README_11_0.png similarity index 100% rename from Example/README_11_0.png rename to Example/3WELL/README_11_0.png diff --git a/Example/README_12_0.png b/Example/3WELL/README_12_0.png similarity index 100% rename from Example/README_12_0.png rename to Example/3WELL/README_12_0.png diff --git a/Example/README_13_0.png b/Example/3WELL/README_13_0.png similarity index 100% rename from Example/README_13_0.png rename to Example/3WELL/README_13_0.png diff --git a/Example/README_14_0.png b/Example/3WELL/README_14_0.png similarity index 100% rename from Example/README_14_0.png rename to Example/3WELL/README_14_0.png diff --git a/Example/README_15_0.png b/Example/3WELL/README_15_0.png similarity index 100% rename from Example/README_15_0.png rename to Example/3WELL/README_15_0.png diff --git a/Example/TRUEPERMX.INC b/Example/3WELL/TRUEPERMX.INC similarity index 100% rename from Example/TRUEPERMX.INC rename to Example/3WELL/TRUEPERMX.INC diff --git a/Example/TINY/PERMX.npy b/Example/TINY/PERMX.npy new file mode 100644 index 0000000..cef50c2 Binary files /dev/null and b/Example/TINY/PERMX.npy differ diff --git a/Example/TINY/RUNFILE.mako b/Example/TINY/RUNFILE.mako new file mode 100644 index 0000000..49fbd67 --- /dev/null +++ b/Example/TINY/RUNFILE.mako @@ -0,0 +1,175 @@ +<%! +import numpy as np +%> + +-- *------------------------------------------* +-- * * +-- * base grid model with input parameters * +-- * * +-- *------------------------------------------* +RUNSPEC + +TITLE + TINY BOX MODEL + +--DIMENS +-- NDIVIX NDIVIY NDIVIZ +-- 40 20 5 / + +--BLACKOIL +OIL +WATER +GAS +DISGAS + +METRIC + +TABDIMS +-- NTSFUN NTPVT NSSFUN NPPVT NTFIP NRPVT NTENDP + 1 1 35 30 5 30 1 / + +EQLDIMS +-- NTEQUL NDRXVD NDPRVD + 1 5 100 / + +WELLDIMS +-- NWMAXZ NCWMAX NGMAXZ MWGMAX + 15 15 2 20 / + +VFPPDIMS +-- MXMFLO MXMTHP MXMWFR MXMGFR MXMALQ NMMVFT + 10 10 10 10 1 1 / + +VFPIDIMS +-- MXSFLO MXSTHP NMSVFT + 10 10 1 / + +AQUDIMS +-- MXNAQN MXNAQC NIFTBL NRIFTB NANAQU NCAMAX + 0 0 1 36 2 200/ + +FAULTDIM + 500 / + +START + 01 JAN 2022 / + +NSTACK + 25 / + + +NOECHO + +GRID +INIT + +INCLUDE +'../include/Grid.grdecl' / +/ + +PERMX +% for i in range(0, len(log_permx)): +% if log_permx[i] < 6: +${"%.3f" %(np.exp(log_permx[i]))} +% else: +${"%.3f" %(np.exp(6))} +% endif +% endfor +/ + +COPY + 'PERMX' 'PERMY' / + 'PERMX' 'PERMZ' / +/ + + +PROPS =============================================================== + +INCLUDE + '../include/pvt.txt' / +/ + +REGIONS =============================================================== + +ENDBOX + +SOLUTION =============================================================== + + +-- DATUM DATUM OWC OWC GOC GOC RSVD RVVD SOLN +-- DEPTH PRESS DEPTH PCOW DEPTH PCOG TABLE TABLE METH +EQUIL + 2355.00 200.46 3000 0.00 2355.0 0.000 / + + +RPTSOL +'PRES' 'SWAT' / + +RPTRST + BASIC=2 / + + + +SUMMARY ================================================================ + +RUNSUM + + +RPTONLY + +WWIR + 'INJ1' + 'INJ2' + 'INJ3' +/ + +WOPR + 'PRO1' + 'PRO2' + 'PRO3' +/ + +WWPR + 'PRO1' + 'PRO2' + 'PRO3' +/ + +ELAPSED + +SCHEDULE ============================================================= + + +RPTSCHED + 'NEWTON=2' / + +RPTRST + BASIC=2 FIP RPORV / + +------------------- WELL SPECIFICATION DATA -------------------------- + +INCLUDE +'../include/Schdl.sch' / +/ + + +WCONINJE +'INJ1' WATER 'OPEN' BHP 2* 300/ +'INJ2' WATER 'OPEN' BHP 2* 300/ +'INJ3' WATER 'OPEN' BHP 2* 300/ +/ + +WCONPROD + 'PRO1' 'OPEN' BHP 5* 90 / + 'PRO2' 'OPEN' BHP 5* 90 / + 'PRO3' 'OPEN' BHP 5* 90 / +/ +--------------------- PRODUCTION SCHEDULE ---------------------------- + + + +TSTEP +10*400 / +/ + +END diff --git a/Example/TINY/gradient_log_permx.png b/Example/TINY/gradient_log_permx.png new file mode 100644 index 0000000..c0a203b Binary files /dev/null and b/Example/TINY/gradient_log_permx.png differ diff --git a/Example/TINY/gradient_permx.png b/Example/TINY/gradient_permx.png new file mode 100644 index 0000000..d3c4157 Binary files /dev/null and b/Example/TINY/gradient_permx.png differ diff --git a/Example/TINY/include/Grid.grdecl b/Example/TINY/include/Grid.grdecl new file mode 100644 index 0000000..17c0611 --- /dev/null +++ b/Example/TINY/include/Grid.grdecl @@ -0,0 +1,1936 @@ +-- Generated [ +-- Exported by: PyReSiTo (single) v1.0 +-- Generated ] +SPECGRID +10 10 2 1 F +/ + +COORD + 0.00 0.00 2355.00 0.00 0.00 2505.00 + 180.00 0.00 2355.00 180.00 0.00 2505.00 + 360.00 0.00 2355.00 360.00 0.00 2505.00 + 540.00 0.00 2355.00 540.00 0.00 2505.00 + 720.00 0.00 2355.00 720.00 0.00 2505.00 + 900.00 0.00 2355.00 900.00 0.00 2505.00 + 1080.00 0.00 2355.00 1080.00 0.00 2505.00 + 1260.00 0.00 2355.00 1260.00 0.00 2505.00 + 1440.00 0.00 2355.00 1440.00 0.00 2505.00 + 1620.00 0.00 2355.00 1620.00 0.00 2505.00 + 1800.00 0.00 2355.00 1800.00 0.00 2505.00 + 0.00 180.00 2355.00 0.00 180.00 2505.00 + 180.00 180.00 2355.00 180.00 180.00 2505.00 + 360.00 180.00 2355.00 360.00 180.00 2505.00 + 540.00 180.00 2355.00 540.00 180.00 2505.00 + 720.00 180.00 2355.00 720.00 180.00 2505.00 + 900.00 180.00 2355.00 900.00 180.00 2505.00 + 1080.00 180.00 2355.00 1080.00 180.00 2505.00 + 1260.00 180.00 2355.00 1260.00 180.00 2505.00 + 1440.00 180.00 2355.00 1440.00 180.00 2505.00 + 1620.00 180.00 2355.00 1620.00 180.00 2505.00 + 1800.00 180.00 2355.00 1800.00 180.00 2505.00 + 0.00 360.00 2355.00 0.00 360.00 2505.00 + 180.00 360.00 2355.00 180.00 360.00 2505.00 + 360.00 360.00 2355.00 360.00 360.00 2505.00 + 540.00 360.00 2355.00 540.00 360.00 2505.00 + 720.00 360.00 2355.00 720.00 360.00 2505.00 + 900.00 360.00 2355.00 900.00 360.00 2505.00 + 1080.00 360.00 2355.00 1080.00 360.00 2505.00 + 1260.00 360.00 2355.00 1260.00 360.00 2505.00 + 1440.00 360.00 2355.00 1440.00 360.00 2505.00 + 1620.00 360.00 2355.00 1620.00 360.00 2505.00 + 1800.00 360.00 2355.00 1800.00 360.00 2505.00 + 0.00 540.00 2355.00 0.00 540.00 2505.00 + 180.00 540.00 2355.00 180.00 540.00 2505.00 + 360.00 540.00 2355.00 360.00 540.00 2505.00 + 540.00 540.00 2355.00 540.00 540.00 2505.00 + 720.00 540.00 2355.00 720.00 540.00 2505.00 + 900.00 540.00 2355.00 900.00 540.00 2505.00 + 1080.00 540.00 2355.00 1080.00 540.00 2505.00 + 1260.00 540.00 2355.00 1260.00 540.00 2505.00 + 1440.00 540.00 2355.00 1440.00 540.00 2505.00 + 1620.00 540.00 2355.00 1620.00 540.00 2505.00 + 1800.00 540.00 2355.00 1800.00 540.00 2505.00 + 0.00 720.00 2355.00 0.00 720.00 2505.00 + 180.00 720.00 2355.00 180.00 720.00 2505.00 + 360.00 720.00 2355.00 360.00 720.00 2505.00 + 540.00 720.00 2355.00 540.00 720.00 2505.00 + 720.00 720.00 2355.00 720.00 720.00 2505.00 + 900.00 720.00 2355.00 900.00 720.00 2505.00 + 1080.00 720.00 2355.00 1080.00 720.00 2505.00 + 1260.00 720.00 2355.00 1260.00 720.00 2505.00 + 1440.00 720.00 2355.00 1440.00 720.00 2505.00 + 1620.00 720.00 2355.00 1620.00 720.00 2505.00 + 1800.00 720.00 2355.00 1800.00 720.00 2505.00 + 0.00 900.00 2355.00 0.00 900.00 2505.00 + 180.00 900.00 2355.00 180.00 900.00 2505.00 + 360.00 900.00 2355.00 360.00 900.00 2505.00 + 540.00 900.00 2355.00 540.00 900.00 2505.00 + 720.00 900.00 2355.00 720.00 900.00 2505.00 + 900.00 900.00 2355.00 900.00 900.00 2505.00 + 1080.00 900.00 2355.00 1080.00 900.00 2505.00 + 1260.00 900.00 2355.00 1260.00 900.00 2505.00 + 1440.00 900.00 2355.00 1440.00 900.00 2505.00 + 1620.00 900.00 2355.00 1620.00 900.00 2505.00 + 1800.00 900.00 2355.00 1800.00 900.00 2505.00 + 0.00 1080.00 2355.00 0.00 1080.00 2505.00 + 180.00 1080.00 2355.00 180.00 1080.00 2505.00 + 360.00 1080.00 2355.00 360.00 1080.00 2505.00 + 540.00 1080.00 2355.00 540.00 1080.00 2505.00 + 720.00 1080.00 2355.00 720.00 1080.00 2505.00 + 900.00 1080.00 2355.00 900.00 1080.00 2505.00 + 1080.00 1080.00 2355.00 1080.00 1080.00 2505.00 + 1260.00 1080.00 2355.00 1260.00 1080.00 2505.00 + 1440.00 1080.00 2355.00 1440.00 1080.00 2505.00 + 1620.00 1080.00 2355.00 1620.00 1080.00 2505.00 + 1800.00 1080.00 2355.00 1800.00 1080.00 2505.00 + 0.00 1260.00 2355.00 0.00 1260.00 2505.00 + 180.00 1260.00 2355.00 180.00 1260.00 2505.00 + 360.00 1260.00 2355.00 360.00 1260.00 2505.00 + 540.00 1260.00 2355.00 540.00 1260.00 2505.00 + 720.00 1260.00 2355.00 720.00 1260.00 2505.00 + 900.00 1260.00 2355.00 900.00 1260.00 2505.00 + 1080.00 1260.00 2355.00 1080.00 1260.00 2505.00 + 1260.00 1260.00 2355.00 1260.00 1260.00 2505.00 + 1440.00 1260.00 2355.00 1440.00 1260.00 2505.00 + 1620.00 1260.00 2355.00 1620.00 1260.00 2505.00 + 1800.00 1260.00 2355.00 1800.00 1260.00 2505.00 + 0.00 1440.00 2355.00 0.00 1440.00 2505.00 + 180.00 1440.00 2355.00 180.00 1440.00 2505.00 + 360.00 1440.00 2355.00 360.00 1440.00 2505.00 + 540.00 1440.00 2355.00 540.00 1440.00 2505.00 + 720.00 1440.00 2355.00 720.00 1440.00 2505.00 + 900.00 1440.00 2355.00 900.00 1440.00 2505.00 + 1080.00 1440.00 2355.00 1080.00 1440.00 2505.00 + 1260.00 1440.00 2355.00 1260.00 1440.00 2505.00 + 1440.00 1440.00 2355.00 1440.00 1440.00 2505.00 + 1620.00 1440.00 2355.00 1620.00 1440.00 2505.00 + 1800.00 1440.00 2355.00 1800.00 1440.00 2505.00 + 0.00 1620.00 2355.00 0.00 1620.00 2505.00 + 180.00 1620.00 2355.00 180.00 1620.00 2505.00 + 360.00 1620.00 2355.00 360.00 1620.00 2505.00 + 540.00 1620.00 2355.00 540.00 1620.00 2505.00 + 720.00 1620.00 2355.00 720.00 1620.00 2505.00 + 900.00 1620.00 2355.00 900.00 1620.00 2505.00 + 1080.00 1620.00 2355.00 1080.00 1620.00 2505.00 + 1260.00 1620.00 2355.00 1260.00 1620.00 2505.00 + 1440.00 1620.00 2355.00 1440.00 1620.00 2505.00 + 1620.00 1620.00 2355.00 1620.00 1620.00 2505.00 + 1800.00 1620.00 2355.00 1800.00 1620.00 2505.00 + 0.00 1800.00 2355.00 0.00 1800.00 2505.00 + 180.00 1800.00 2355.00 180.00 1800.00 2505.00 + 360.00 1800.00 2355.00 360.00 1800.00 2505.00 + 540.00 1800.00 2355.00 540.00 1800.00 2505.00 + 720.00 1800.00 2355.00 720.00 1800.00 2505.00 + 900.00 1800.00 2355.00 900.00 1800.00 2505.00 + 1080.00 1800.00 2355.00 1080.00 1800.00 2505.00 + 1260.00 1800.00 2355.00 1260.00 1800.00 2505.00 + 1440.00 1800.00 2355.00 1440.00 1800.00 2505.00 + 1620.00 1800.00 2355.00 1620.00 1800.00 2505.00 + 1800.00 1800.00 2355.00 1800.00 1800.00 2505.00 +/ +ZCORN +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2355.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2430.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +2505.00 +/ +ACTNUM +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +/ +PORO +200*0.2 / diff --git a/Example/TINY/include/Schdl.sch b/Example/TINY/include/Schdl.sch new file mode 100644 index 0000000..34c51f0 --- /dev/null +++ b/Example/TINY/include/Schdl.sch @@ -0,0 +1,26 @@ + +------------------- WELL SPECIFICATION DATA -------------------------- +WELSPECS +INJ1 G 1 1 2357 WATER 1* STD 1* 1* 1* / +INJ2 G 5 1 2357 WATER 1* STD 1* 1* 1* / +INJ3 G 10 1 2357 WATER 1* STD 1* 1* 1* / +PRO1 G 1 10 2357 WATER 1* STD 1* 1* 1* / +PRO2 G 5 10 2357 WATER 1* STD 1* 1* 1* / +PRO3 G 10 10 2357 WATER 1* STD 1* 1* 1* / +/ + +COMPDAT +INJ1 1 1 1 1 OPEN 1* 1* 0.15 1* 5.0 / +INJ1 1 1 2 2 OPEN 1* 1* 0.15 1* 5.0 / +INJ2 5 1 1 1 OPEN 1* 1* 0.15 1* 5.0 / +INJ2 5 1 2 2 OPEN 1* 1* 0.15 1* 5.0 / +INJ3 10 1 1 1 OPEN 1* 1* 0.15 1* 5.0 / +INJ3 10 1 2 2 OPEN 1* 1* 0.15 1* 5.0 / +PRO1 1 10 1 1 OPEN 1* 1* 0.15 1* 5.0 / +PRO1 1 10 2 2 OPEN 1* 1* 0.15 1* 5.0 / +PRO2 5 10 1 1 OPEN 1* 1* 0.15 1* 5.0 / +PRO2 5 10 2 2 OPEN 1* 1* 0.15 1* 5.0 / +PRO3 10 10 1 1 OPEN 1* 1* 0.15 1* 5.0 / +PRO3 10 10 2 2 OPEN 1* 1* 0.15 1* 5.0 / +/ + diff --git a/Example/TINY/include/Schdl_fixed_wi.sch b/Example/TINY/include/Schdl_fixed_wi.sch new file mode 100644 index 0000000..5ffa971 --- /dev/null +++ b/Example/TINY/include/Schdl_fixed_wi.sch @@ -0,0 +1,26 @@ + +------------------- WELL SPECIFICATION DATA -------------------------- +WELSPECS +INJ1 G 1 1 2357 WATER 1* STD 1* 1* 1* / +INJ2 G 5 1 2357 WATER 1* STD 1* 1* 1* / +INJ3 G 10 1 2357 WATER 1* STD 1* 1* 1* / +PRO1 G 1 10 2357 WATER 1* STD 1* 1* 1* / +PRO2 G 5 10 2357 WATER 1* STD 1* 1* 1* / +PRO3 G 10 10 2357 WATER 1* STD 1* 1* 1* / +/ + +COMPDAT +INJ1 1 1 1 1 OPEN 1* 4.51401363e+01 0.15 1* 5.0 / +INJ1 1 1 2 2 OPEN 1* 4.52884317e+01 0.15 1* 5.0 / +INJ2 5 1 1 1 OPEN 1* 4.17800918e+01 0.15 1* 5.0 / +INJ2 5 1 2 2 OPEN 1* 3.59314189e+01 0.15 1* 5.0 / +INJ3 10 1 1 1 OPEN 1* 1.55353913e+01 0.15 1* 5.0 / +INJ3 10 1 2 2 OPEN 1* 1.25381669e+01 0.15 1* 5.0 / +PRO1 1 10 1 1 OPEN 1* 9.87532715e+00 0.15 1* 5.0 / +PRO1 1 10 2 2 OPEN 1* 9.67052102e+00 0.15 1* 5.0 / +PRO2 5 10 1 1 OPEN 1* 8.36501688e+00 0.15 1* 5.0 / +PRO2 5 10 2 2 OPEN 1* 7.23084441e+00 0.15 1* 5.0 / +PRO3 10 10 1 1 OPEN 1* 1.19734202e+01 0.15 1* 5.0 / +PRO3 10 10 2 2 OPEN 1* 1.14788620e+01 0.15 1* 5.0 / +/ + diff --git a/Example/TINY/include/pvt.txt b/Example/TINY/include/pvt.txt new file mode 100644 index 0000000..145fdce --- /dev/null +++ b/Example/TINY/include/pvt.txt @@ -0,0 +1,319 @@ +-- +-- Water +SWFN + 0.000000 0.000000 0.000000 + 0.006819 0.003554 0.000000 + 0.027091 0.013908 0.000000 + 0.060263 0.025952 0.000000 + 0.105430 0.044108 0.000000 + 0.161359 0.068765 0.000000 + 0.226526 0.101539 0.000000 + 0.299152 0.145360 0.000000 + 0.377257 0.193942 0.000000 + 0.458710 0.244954 0.000000 + 0.541290 0.308844 0.000000 + 0.622743 0.394668 0.000000 + 0.700848 0.491067 0.000000 + 0.773474 0.593135 0.000000 + 0.838641 0.699496 0.000000 + 0.894570 0.800074 0.000000 + 0.939737 0.888146 0.000000 + 0.972909 0.949655 0.000000 + 0.993181 0.987216 0.000000 + 1.000000 1.000000 0.000000 +/ +-- Oil +SOF3 + 0.000000 0.000000 0.000000 + 0.006819 0.000686 0.002515 + 0.027091 0.002746 0.013164 + 0.060263 0.006794 0.034361 + 0.105430 0.014497 0.067229 + 0.161359 0.026915 0.112035 + 0.226526 0.045924 0.168322 + 0.299152 0.073425 0.235001 + 0.377257 0.111978 0.310431 + 0.458710 0.166662 0.392506 + 0.541290 0.228933 0.478758 + 0.622743 0.297987 0.566462 + 0.700848 0.383686 0.652753 + 0.773474 0.495582 0.734742 + 0.838641 0.617999 0.809639 + 0.894570 0.732288 0.874857 + 0.939737 0.839305 0.928127 + 0.972909 0.926725 0.967580 + 0.993181 0.980640 0.991823 + 1.000000 1.000000 1.000 +/ +-- Gas +SGFN + 0.0 0.0 0.000 + 0.006819 0.000563 0.000 + 0.027091 0.004459 0.000 + 0.060263 0.014794 0.000 + 0.105430 0.034233 0.000 + 0.161359 0.064817 0.000 + 0.226526 0.107814 0.000 + 0.299152 0.163621 0.000 + 0.377257 0.231716 0.000 + 0.458710 0.310676 0.000 + 0.541290 0.398240 0.000 + 0.622743 0.491432 0.000 + 0.700848 0.586727 0.000 + 0.773474 0.680250 0.000 + 0.838641 0.768005 0.000 + 0.894570 0.846100 0.000 + 0.939737 0.910981 0.000 + 0.972909 0.959640 0.000 + 0.993181 0.989789 0.000 + 1.0 1.000000 0.000 +/ + +-- *** G A S D A T A *** +-- BAR (CP) +-- +PVTG +-- + 20.00 0.00002448 0.061895 0.01299 + 0.00001224 0.061810 0.01300 + 0.00000000 0.061725 0.01300 / + 40.00 0.00000628 0.030252 0.01383 + 0.00000314 0.030249 0.01383 + 0.00000000 0.030245 0.01383 / + 60.00 0.00000585 0.019844 0.01450 + 0.00000292 0.019845 0.01450 + 0.00000000 0.019846 0.01449 / + 80.00 0.00000728 0.014686 0.01520 + 0.00000364 0.014689 0.01519 + 0.00000000 0.014692 0.01518 / + 100.00 0.00001017 0.011627 0.01596 + 0.00000509 0.011633 0.01595 + 0.00000000 0.011638 0.01593 / + 120.00 0.00001485 0.009619 0.01682 + 0.00000743 0.009627 0.01679 + 0.00000000 0.009635 0.01676 / + 140.00 0.00002182 0.008213 0.01780 + 0.00001091 0.008224 0.01774 + 0.00000000 0.008235 0.01767 / + 160.00 0.00003155 0.007184 0.01890 + 0.00001577 0.007198 0.01878 + 0.00000000 0.007212 0.01866 / + 197.66 0.00006327 0.005820 0.02160 + 0.00003164 0.005840 0.02122 + 0.00000000 0.005860 0.02086 / + 231.13 0.00010861 0.005042 0.02477 + 0.00005431 0.005061 0.02389 + 0.00000000 0.005082 0.02306 / + 261.31 0.00016781 0.004561 0.02844 + 0.00008391 0.004571 0.02672 + 0.00000000 0.004584 0.02515 / + 288.87 0.00024205 0.004255 0.03272 + 0.00012103 0.004243 0.02976 + 0.00000000 0.004241 0.02711 / + 314.34 0.00033405 0.004062 0.03783 + 0.00016703 0.004017 0.03311 + 0.00000000 0.003990 0.02893 / + 338.20 0.00044866 0.003953 0.04410 + 0.00022433 0.003860 0.03693 + 0.00000000 0.003797 0.03065 / + 360.83 0.00059341 0.003915 0.05210 + 0.00029670 0.003756 0.04150 + 0.00000000 0.003644 0.03227 / + 382.58 0.00077814 0.003947 0.06273 + 0.00038907 0.003698 0.04725 + 0.00000000 0.003518 0.03382 / + 403.60 0.00100943 0.004048 0.07723 + 0.00050471 0.003683 0.05472 + 0.00000000 0.003413 0.03529 / + 423.77 0.00127517 0.004207 0.09631 + 0.00063758 0.003705 0.06418 + 0.00000000 0.003325 0.03664 / +/ + +-- *** O I L D A T A *** + +-- RSO PRESSURE B-OIL VISCOSITY +-- (BAR) (CP) +PVTO +-- + 19.6 20.00 1.12324 0.96519 + 55.00 1.11698 1.03237 + 90.00 1.11127 1.10051 + 125.00 1.10602 1.16942 + 160.00 1.10119 1.23890 + 195.00 1.09672 1.30876 + 230.00 1.09256 1.37884 + 265.00 1.08868 1.44899 + 300.00 1.08504 1.51908 + 335.00 1.08164 1.58903 + 370.00 1.07843 1.65876 / + + 31.5 40.00 1.15981 0.85738 + 75.00 1.15288 0.91402 + 110.00 1.14657 0.97137 + 145.00 1.14079 1.02927 + 180.00 1.13546 1.08759 + 215.00 1.13053 1.14617 + 250.00 1.12595 1.20488 + 285.00 1.12168 1.26360 + 320.00 1.11768 1.32224 + 355.00 1.11394 1.38073 + 390.00 1.11042 1.43898 / + + 42.4 60.00 1.191 0.7868 + 95.00 1.184 0.8364 + 130.00 1.177 0.8866 + 165.00 1.171 0.9371 + 200.00 1.165 0.9880 + 235.00 1.160 1.0390 + 270.00 1.155 1.0902 + 305.00 1.150 1.1413 + 340.00 1.146 1.1922 + 375.00 1.142 1.2431 + 410.00 1.138 1.2936 / +-- + 53.4 80.00 1.222 0.7175 + 115.00 1.214 0.7608 + 150.00 1.206 0.8045 + 185.00 1.200 0.8485 + 220.00 1.194 0.8928 + 255.00 1.188 0.9371 + 290.00 1.183 0.9815 + 325.00 1.178 1.0258 + 360.00 1.173 1.0700 + 395.00 1.169 1.1141 + 430.00 1.165 1.1579 / +-- + 64.6 100.00 1.252 0.6544 + 135.00 1.244 0.6923 + 170.00 1.236 0.7305 + 205.00 1.229 0.7689 + 240.00 1.222 0.8075 + 275.00 1.216 0.8461 + 310.00 1.211 0.8847 + 345.00 1.205 0.9233 + 380.00 1.200 0.9618 + 415.00 1.196 1.0000 + 450.00 1.191 1.0381 / +-- + 76.3 120.00 1.284 0.5978 + 155.00 1.275 0.6312 + 190.00 1.266 0.6648 + 225.00 1.259 0.6985 + 260.00 1.252 0.7323 + 295.00 1.245 0.7661 + 330.00 1.239 0.7999 + 365.00 1.234 0.8337 + 400.00 1.229 0.8673 + 435.00 1.224 0.9007 + 470.00 1.219 0.9340 / +-- + 88.5 140.00 1.316 0.5477 + 175.00 1.307 0.5749 + 210.00 1.298 0.6020 + 245.00 1.290 0.6290 + 280.00 1.282 0.6559 + 315.00 1.276 0.6827 + 350.00 1.269 0.7095 + 385.00 1.263 0.7362 + 420.00 1.258 0.7629 + 455.00 1.253 0.7895 + 490.00 1.248 0.8161 / +-- + 101.3 160.00 1.350 0.5020 + 195.00 1.340 0.5227 + 230.00 1.331 0.5432 + 265.00 1.322 0.5635 + 300.00 1.314 0.5835 + 335.00 1.307 0.6034 + 370.00 1.300 0.6231 + 405.00 1.294 0.6426 + 440.00 1.288 0.6620 + 475.00 1.283 0.6813 / +-- + 114.7 180.00 1.385 0.4636 + 215.00 1.375 0.4820 + 250.00 1.365 0.5003 + 285.00 1.356 0.5183 + 320.00 1.347 0.5362 + 355.00 1.340 0.5538 + 390.00 1.333 0.5712 + 425.00 1.326 0.5885 + 460.00 1.320 0.6055 + 495.00 1.314 0.6222 / +-- + 128.9 200.00 1.422 0.4290 + 235.00 1.411 0.4455 + 270.00 1.401 0.4618 + 305.00 1.391 0.4779 + 340.00 1.382 0.4938 + 375.00 1.374 0.5096 + 410.00 1.367 0.5252 + 445.00 1.360 0.5406 + 480.00 1.353 0.5558 / +-- + 143.8 220.00 1.461 0.3977 + 255.00 1.449 0.4125 + 290.00 1.438 0.4271 + 325.00 1.428 0.4415 + 360.00 1.419 0.4558 + 395.00 1.410 0.4699 + 430.00 1.402 0.4839 + 465.00 1.395 0.4977 / +-- + 159.5 240.00 1.502 0.3692 + 275.00 1.489 0.3825 + 310.00 1.478 0.3956 + 345.00 1.467 0.4086 + 380.00 1.458 0.4214 + 415.00 1.449 0.4341 + 450.00 1.440 0.4466 + 485.00 1.432 0.4590 / +-- + 184.0 268.79 1.565 0.3324 + 303.79 1.551 0.3438 + 338.79 1.539 0.3551 + 373.79 1.528 0.3663 + 408.79 1.517 0.3774 + 443.79 1.508 0.3883 + 478.79 1.499 0.3991 / +-- + 226.3 306.18 1.679 0.2855 + 341.18 1.664 0.2949 + 376.18 1.650 0.3041 + 411.18 1.637 0.3132 + 446.18 1.625 0.3222 + 481.18 1.614 0.3311 / +-- + 268.6 339.93 1.792 0.2517 + 374.93 1.775 0.2597 + 409.93 1.760 0.2675 + 444.93 1.746 0.2751 + 479.93 1.732 0.2827 / +-- + 310.9 371.44 1.903 0.2265 + 406.44 1.885 0.2333 + 441.44 1.868 0.2401 + 476.44 1.853 0.2468 / +-- + 353.3 401.66 2.013 0.2071 + 436.66 1.993 0.2132 + 471.66 1.975 0.2192 / +/ +-- +--****************************************************** +--* PVT WATER +--****************************************************** +PVTW + 344.83 1.0292 4.002E-05 0.36000 0.00E+00 / + +-------------------------------------------------------- +-- +--P(DATUM) CR +ROCK + 383.0 4.12E-05 / +-- + + +DENSITY + 842.3 1001.1 0.900 / diff --git a/Example/TINY/permx.png b/Example/TINY/permx.png new file mode 100644 index 0000000..9655b27 Binary files /dev/null and b/Example/TINY/permx.png differ diff --git a/Example/TINY/rates.png b/Example/TINY/rates.png new file mode 100644 index 0000000..40df74c Binary files /dev/null and b/Example/TINY/rates.png differ diff --git a/Example/TINY/run_sim.py b/Example/TINY/run_sim.py new file mode 100644 index 0000000..7fc0a70 --- /dev/null +++ b/Example/TINY/run_sim.py @@ -0,0 +1,287 @@ +import time +import numpy as np +import matplotlib.pyplot as plt +from datetime import datetime +from subsurface.multphaseflow.jutul_darcy import JutulDarcy + +# Report steps +datetimes = [ + datetime(2023, 2, 5), + datetime(2024, 3, 11), + datetime(2025, 4, 15), + datetime(2026, 5, 20), + datetime(2027, 6, 24), + datetime(2028, 7, 28), + datetime(2029, 9, 1), + datetime(2030, 10, 6), + datetime(2031, 11, 10), + datetime(2032, 12, 14), +] + +# Data types to report +datatype = [ + 'WOPR:PRO1', + 'WOPR:PRO2', + 'WOPR:PRO3', + 'WWPR:PRO1', + 'WWPR:PRO2', + 'WWPR:PRO3', +] + +# Simulator settings and adjoint configuration +kwargs = { + 'reporttype': 'dates', + 'reportpoint': datetimes, + 'runfile': 'RUNFILE.mako', + 'datatype': datatype, + 'adjoint_pbar': True, + 'perm_copied': True, # Include total derivative when PERMX is copied to PERMY and PERMZ + 'adjoints': {'WOPR': + {'steps': [datetime(2032, 12, 14)], 'wellID': 'PRO2', 'parameters': ['log_permx', 'permx']}, + }, +} + +def simulate(): + log_permx = np.log(np.load('PERMX.npy')) + simulator = JutulDarcy(kwargs) + res, grad = simulator({'log_permx': log_permx}) + return res, grad + +def plot(res): + """Plot production rates with enhanced styling.""" + plt.style.use('seaborn-v0_8-darkgrid') + + fig, ax = plt.subplots(ncols=3, figsize=(12, 3), sharey=True) + fig.patch.set_facecolor('white') + + colors = {'WOPR': '#2E86AB', 'WWPR': '#A23B72'} + linestyles = {'WOPR': '-', 'WWPR': '--'} + linewidths = {'WOPR': 2.5, 'WWPR': 2.2} + + for col in range(3): + producer_num = col + 1 + + # Plot well production rates + ax[col].plot( + res.index, res[f'WOPR:PRO{producer_num}'], + label='Oil Production', color=colors['WOPR'], + linestyle=linestyles['WOPR'], linewidth=linewidths['WOPR'], + marker='o', markersize=5, alpha=0.85 + ) + ax[col].plot( + res.index, res[f'WWPR:PRO{producer_num}'], + label='Water Production', color=colors['WWPR'], + linestyle=linestyles['WWPR'], linewidth=linewidths['WWPR'], + marker='s', markersize=4, alpha=0.85 + ) + + # Enhanced styling + ax[col].set_title(f'Producer {producer_num}', fontsize=14, fontweight='bold', pad=15) + ax[col].set_ylabel('Rate (m³/day)', fontsize=11, fontweight='500') + ax[col].legend(loc='best', framealpha=0.95, fontsize=10, edgecolor='gray') + ax[col].grid(True, alpha=0.7, linewidth=0.8) + ax[col].spines['top'].set_visible(False) + ax[col].spines['right'].set_visible(False) + + # Rotate x-axis labels for better readability + ax[col].tick_params(axis='x', rotation=45) + + plt.tight_layout() + fig.savefig('rates.png', dpi=300, bbox_inches='tight') + + +def plot_permx(): + """Plot the log permeability field as a 3D voxel visualization, including wells as small sticks above the grid.""" + nx, ny, nz = 10, 10, 2 + permx = np.load('PERMX.npy').reshape((nx, ny, nz), order='F') + + wells = { + "INJ1": {"ij": (0, 0), "color": "deepskyblue"}, + "INJ2": {"ij": (4, 0), "color": "deepskyblue"}, + "INJ3": {"ij": (9, 0), "color": "deepskyblue"}, + "PRO1": {"ij": (0, 9), "color": "crimson"}, + "PRO2": {"ij": (4, 9), "color": "crimson"}, + "PRO3": {"ij": (9, 9), "color": "crimson"}, + } + + permx_norm = (permx - permx.min()) / (permx.max() - permx.min()) + filled = np.ones((nx, ny, nz), dtype=bool) + + cmap = plt.get_cmap('YlGnBu_r') + facecolors = cmap(permx_norm) + edgecolor_intensity = permx_norm * 0.5 + 0.5 + edgecolors = plt.cm.gray(edgecolor_intensity) + + fig = plt.figure(figsize=(10, 6)) + ax = fig.add_subplot(111, projection='3d') + ax.computed_zorder = False # allow manual zorder in 3D + fig.subplots_adjust(left=0.05, right=0.95, top=0.95, bottom=0.12) + + x, y, z = np.indices(np.array(filled.shape) + 1).astype(float) + x = x / nx + y = y / ny + z = z / nz + + ax.voxels( + x, y, z, filled, + facecolors=facecolors, + edgecolors=edgecolors, + linewidth=0.5, + alpha=1.0, + zsort='max' + ) + + # Very thin, taller well sticks: centered in each cell + stick_extra_above = 0.75 # taller above z=1.0 (was 0.30) + stick_size_x = 0.15 / nx # thinner + stick_size_y = 0.15 / ny # thinner + + for name, w in wells.items(): + i, j = w["ij"] + color = w.get("color", "black") + + # exact cell center in normalized coordinates + cx = (i + 0.5) / nx + cy = (j + 0.5) / ny + + # bar3d expects lower-left corner, so shift by half size to keep centered + ax.bar3d( + cx - 0.5 * stick_size_x, cy - 0.5 * stick_size_y, 1.0, + stick_size_x, stick_size_y, 1.0 + stick_extra_above, + color=color, edgecolor=None, linewidth=0.8, shade=True, alpha=0.7, zsort='max', + ) + ax.text(cx, cy, 1.0 + stick_extra_above + 1.2, name, color=color, fontsize=9, ha='center', zorder=1) + + ax.set_zlim(0.0, 1.0 + stick_extra_above + 0.05) + + sm = plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=permx.min(), vmax=permx.max())) + sm.set_array([]) + cbar_ax = fig.add_axes([0.15, 0.2, 0.7, 0.03]) + cbar = plt.colorbar(sm, cax=cbar_ax, orientation='horizontal') + cbar.set_label('Permeability (mD)', fontsize=11, fontweight='bold') + + ax.set_xticklabels([]) + ax.set_yticklabels([]) + ax.set_zticklabels([]) + + ax.view_init(elev=20, azim=45) + ax.set_box_aspect([nx/10, ny/10, nz/10]) + fig.savefig('permx.png', dpi=300, bbox_inches='tight') + + +def plot_gradient(grad, savename='gradient_log_permx.png'): + nx, ny, nz = 10, 10, 2 + grad = grad.reshape((nx, ny, nz), order='F') + + wells = { + "INJ1": {"ij": (0, 0), "color": "deepskyblue"}, + "INJ2": {"ij": (4, 0), "color": "deepskyblue"}, + "INJ3": {"ij": (9, 0), "color": "deepskyblue"}, + "PRO1": {"ij": (0, 9), "color": "crimson"}, + "PRO2": {"ij": (4, 9), "color": "crimson"}, + "PRO3": {"ij": (9, 9), "color": "crimson"}, + } + + # Normalize gradient for color mapping (0-1) + grad_norm = (grad - grad.min()) / (grad.max() - grad.min()) + + # Create voxel array (all filled) + filled = np.ones((nx, ny, nz), dtype=bool) + + # Create color map: low gradient (dark blue) to high gradient (bright yellow) + cmap = plt.get_cmap('YlGnBu_r') + rgba_colors = cmap(grad_norm) + facecolors = rgba_colors + + # Edge colors based on intensity + edgecolor_intensity = grad_norm * 0.5 + 0.5 # Scale to [0.5, 1] + edgecolors = plt.cm.gray(edgecolor_intensity) + + # Create figure with 3D axes and space for colorbar + fig = plt.figure(figsize=(10, 6)) + ax = fig.add_subplot(111, projection='3d') + ax.computed_zorder = False # allow manual zorder in 3D + + # Minimize whitespace by adjusting margins + fig.subplots_adjust(left=0.05, right=0.95, top=0.95, bottom=0.12) + + # Create coordinate arrays + x, y, z = np.indices(np.array(filled.shape) + 1).astype(float) + x = x / nx + y = y / ny + z = z / nz + + # Plot voxels + ax.voxels( + x, y, z, filled, + facecolors=facecolors, + edgecolors=edgecolors, + linewidth=0.5, + alpha=1.0, + zsort='max' + ) + + # Very thin, taller well sticks: centered in each cell + stick_extra_above = 0.75 # taller above z=1.0 (was 0.30) + stick_size_x = 0.15 / nx # thinner + stick_size_y = 0.15 / ny # thinner + + for name, w in wells.items(): + i, j = w["ij"] + color = w.get("color", "black") + + # exact cell center in normalized coordinates + cx = (i + 0.5) / nx + cy = (j + 0.5) / ny + + # bar3d expects lower-left corner, so shift by half size to keep centered + ax.bar3d( + cx - 0.5 * stick_size_x, cy - 0.5 * stick_size_y, 1.0, + stick_size_x, stick_size_y, 1.0 + stick_extra_above, + color=color, edgecolor=None, linewidth=0.8, shade=True, alpha=0.7, zsort='max', + ) + ax.text(cx, cy, 1.0 + stick_extra_above + 1.2, name, color=color, fontsize=9, ha='center', zorder=1) + + ax.set_zlim(0.0, 1.0 + stick_extra_above + 0.05) + + + # Add horizontal colorbar below the plot + sm = plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=grad.min(), vmax=grad.max())) + sm.set_array([]) + cbar_ax = fig.add_axes([0.15, 0.2, 0.7, 0.03]) + cbar = plt.colorbar(sm, cax=cbar_ax, orientation='horizontal') + cbar.set_label('Gradient of log-PERMX (Sm3/day)', fontsize=11, fontweight='bold') + + ax.set_xticklabels([]) + ax.set_yticklabels([]) + ax.set_zticklabels([]) + + # View angle + ax.view_init(elev=20, azim=45) + ax.set_box_aspect([nx/10, ny/10, nz/10]) + plt.show() + fig.savefig(savename, dpi=300, bbox_inches='tight') + + +if __name__ == "__main__": + + # Visualize permeability field + plot_permx() + + # Run simulation + results, gradient = simulate() + print(results) + print(gradient) + + # Plot gradient field for report time step + grad = gradient.loc[datetime(2032, 12, 14), ('WOPR:PRO2', 'log_permx')] + plot_gradient(grad, savename='gradient_log_permx.png') + + # Plot production rates + plot(results) + + # Test that gradient with respect to log_permx is consistent with gradient with respect to permx + permx = np.load('PERMX.npy') + grad_permx = gradient.loc[datetime(2032, 12, 14), ('WOPR:PRO2', 'permx')] + grad_log_permx = gradient.loc[datetime(2032, 12, 14), ('WOPR:PRO2', 'log_permx')] + np.testing.assert_almost_equal(grad_log_permx, grad_permx * permx) diff --git a/Example/TINY/test_grad.py b/Example/TINY/test_grad.py new file mode 100644 index 0000000..f6ecb96 --- /dev/null +++ b/Example/TINY/test_grad.py @@ -0,0 +1,286 @@ +import os +import numpy as np +import matplotlib.pyplot as plt + +from datetime import datetime +from subsurface.multphaseflow.jutul_darcy import JutulDarcy +from misc.structures import PETDataFrame, PETStateArray + + +def test_finite_difference_gradient(run_fda=True, run_adjoint=True, folder='TEST/FINITE_DIFF'): + + datapoint = 'WOPR:PRO2' + date = datetime(2032, 12, 14) + options = { + 'parallel': 5, + 'datatype': [datapoint], + 'reporttype': 'dates', + 'reportpoint': [date], + 'runfile': 'RUNFILE.mako', + 'startdate': datetime(2022, 1, 1), + 'adjoint_pbar': False, + 'adjoints': {'WOPR': + {'steps': [date], 'wellID': 'PRO2', 'parameters': 'poro'} + }, + 'perm_copied': True, + } + + os.makedirs(folder, exist_ok=True) + + # Load PERMX + permx = np.load('PERMX.npy') + + reps = 0.01 # Relative perturbation size for finite difference approximation (0.5%) + if run_fda: + # Calculate finite difference gradient + kw = dict(options) + kw.pop('adjoints') + + # Pertubed input + permx_p = np.tile(permx[:, None], (1, permx.size)) + permx_m = np.tile(permx[:, None], (1, permx.size)) + for i in range(permx.size): + permx_p[i, i] += reps*permx[i] + permx_m[i, i] -= reps*permx[i] + + # Run simulator on perturbed inputs + simulator = JutulDarcy(kw) + inputs_p = [{'log_permx': np.log(permx_p[:, i])} for i in range(permx.size)] + inputs_m = [{'log_permx': np.log(permx_m[:, i])} for i in range(permx.size)] + results_p = simulator(inputs_p) + results_m = simulator(inputs_m) + + # Convert results to PETDataFrame and save + results_p = PETDataFrame.merge_dataframes(results_p) + results_m = PETDataFrame.merge_dataframes(results_m) + results_p.to_pickle(os.path.join(folder, 'results_p.pkl')) + results_m.to_pickle(os.path.join(folder, 'results_m.pkl')) + + if run_adjoint: + # Run simulator with adjoint to get gradient + simulator = JutulDarcy(options) + inputs = [{'log_permx': np.log(permx)}] + results, adjoint = simulator(inputs) + adjoint.to_pickle(os.path.join(folder, 'adjoint.pkl')) + + # Analyze results + results_p = PETDataFrame.from_pickle(os.path.join(folder, 'results_p.pkl')).loc[date, datapoint] + results_m = PETDataFrame.from_pickle(os.path.join(folder, 'results_m.pkl')).loc[date, datapoint] + grad_fda = (results_p - results_m)/(2*reps*permx) # Finite difference approximation of gradient + grad_adj = PETDataFrame.from_pickle(os.path.join(folder, 'adjoint.pkl')).loc[date, (datapoint, 'permx')] + + # Plot FDA vs adjoint gradient on grid + nx = 10 + ny = 10 + nz = 2 + grad_fda_grid = grad_fda.reshape((nx, ny, nz), order='F') + grad_adj_grid = grad_adj.reshape((nx, ny, nz), order='F') + + fig, axes = plt.subplots(2, 2, figsize=(10, 8), dpi=140) + cmap = 'seismic' + + # im0: Finite difference gradient for layer 1 + maxv = np.max(np.abs(grad_adj_grid[:, :, 0])) + im0 = axes[0, 0].imshow(grad_fda_grid[:, :, 0], cmap=cmap, vmin=-maxv, vmax=maxv) + axes[0, 0].set_title('FDA Gradient (Layer 1)', fontsize=11) + fig.colorbar(im0, ax=axes[0, 0], fraction=0.046, pad=0.04) + + # im1: Adjoint gradient for layer 1 + im1 = axes[0, 1].imshow(grad_adj_grid[:, :, 0], cmap=cmap, vmin=-maxv, vmax=maxv) + axes[0, 1].set_title('Adjoint Gradient (Layer 1)', fontsize=11) + fig.colorbar(im1, ax=axes[0, 1], fraction=0.046, pad=0.04) + + + # im2: Finite difference gradient for layer 2 + maxv = np.max(np.abs(grad_adj_grid[:, :, 1])) + im2 = axes[1, 0].imshow(grad_fda_grid[:, :, 1], cmap=cmap, vmin=-maxv, vmax=maxv) + axes[1, 0].set_title('FDA Gradient (Layer 2)', fontsize=11) + fig.colorbar(im2, ax=axes[1, 0], fraction=0.046, pad=0.04) + + # im3: Adjoint gradient for layer 2 + im3 = axes[1, 1].imshow(grad_adj_grid[:, :, 1], cmap=cmap, vmin=-maxv, vmax=maxv) + axes[1, 1].set_title('Adjoint Gradient (Layer 2)', fontsize=11) + fig.colorbar(im3, ax=axes[1, 1], fraction=0.046, pad=0.04) + + for ax in axes.flatten(): + ax.set_xticks([]) + ax.set_yticks([]) + + fig.tight_layout() + fig.savefig(os.path.join(folder, 'fda_vs_adjoint_gradient.png'), dpi=300) + #plt.show() + + grad_fda_vec = grad_fda_grid.flatten(order='F') + grad_adj_vec = grad_adj_grid.flatten(order='F') + + relative_error = np.linalg.norm(grad_fda_vec - grad_adj_vec) / np.linalg.norm(grad_adj_vec) + print(f"Relative error between FDA and adjoint gradients: {relative_error:.4e}") + print(f"Norm of FDA gradient: {np.linalg.norm(grad_fda_vec):.4e}") + print(f"Norm of adjoint gradient: {np.linalg.norm(grad_adj_vec):.4e}") + + + rtol = grad_adj_vec.size * np.finfo(float).eps * 100 # Relative tolerance scaled by size of gradient vector + atol = 1e-2 * np.linalg.norm(grad_adj_vec) # + print(f"Using rtol={rtol:.4e} and atol={atol:.4e} for np.testing.assert_allclose") + np.testing.assert_allclose(grad_fda_vec, grad_adj_vec, rtol=rtol, atol=atol) + + +def test_sens_matrix_of_log_permx(run=True, folder='TEST'): + + datapoint = 'WOPR:PRO2' + date = datetime(2032, 12, 14) + options = { + 'parallel': 5, + 'datatype': [datapoint], + 'reporttype': 'dates', + 'reportpoint': [date], + 'runfile': 'RUNFILE.mako', + 'startdate': datetime(2022, 1, 1), + 'adjoint_pbar': False, + 'adjoints': {'WOPR': + {'steps': [date], 'wellID': 'PRO2', 'parameters': 'log_permx'} + }, + } + + os.makedirs(folder, exist_ok=True) + + ne = 10_000 + if run: + np.random.seed(29_01_1983) + + pinfo = { + 'nx': 10, + 'ny': 10, + 'nz': 2, + 'vario': ['sph', 'sph'], + 'mean': 200*[4.0], + 'variance': [1.0, 1.0], + 'corr_length': [10.0, 10.0], + 'aniso': [1.0, 1.0], + 'angle': [0.0, 0.0], + } + prior_log_permx_ensemble = PETStateArray.generate_from_prior_info( + prior_info = {'log_permx': pinfo}, + ne=ne, + save=False + ) + np.save(os.path.join(folder, 'prior_log_permx_ensemble.npy'), prior_log_permx_ensemble) + + # Run simulator on ensemble + simulator = JutulDarcy(options) + inputs = [{'log_permx': prior_log_permx_ensemble[:, i]} for i in range(ne)] + results, adjoint = simulator(inputs) + results = PETDataFrame.merge_dataframes(results) + adjoint = PETDataFrame.merge_dataframes(adjoint) + results.to_pickle(os.path.join(folder, 'results.pkl')) + adjoint.to_pickle(os.path.join(folder, 'adjoint.pkl')) + + + # Analyze results + col = datapoint + idx = date + results = PETDataFrame.from_pickle(os.path.join(folder, 'results.pkl')).loc[idx, col] + adjoint = PETDataFrame.from_pickle(os.path.join(folder, 'adjoint.pkl')).loc[idx, (col, 'log_permx')] + + nx = 200 + ny = 1 + enX = np.load(os.path.join(folder, 'prior_log_permx_ensemble.npy')) + enY = results[np.newaxis, :] + enG = adjoint[np.newaxis, :, :] + + assert enX.shape == (nx, ne) + assert enY.shape == (ny, ne) + assert enG.shape == (ny, nx, ne) + + # Compute sensitivity matrix using ensemble gradients + P = (np.eye(ne) - np.ones((ne,ne))/ne)/np.sqrt(ne-1) + A = enX @ P + Y = enY @ P + Gbar = np.mean(enG, axis=-1) + + Cyx = Y @ A.T + GbarCxx = Gbar @ A @ A.T + + # ------------------------------------------------------------------- + # Stein's lemma tells us that: + # E[G]Cxx = Cyx ---> Gbar @ A @ A.T ≈ Y @ A.T (for large ne) + # ------------------------------------------------------------------- + + # Loop over different ensemble sizes and compute norms of Cyx, GbarCxx, and their difference + Cxy_norm = [] + GbarCxx_norm = [] + err = [] + ens = np.logspace(1, np.log10(ne), num=10, dtype=int) + for n in ens: + P_n = (np.eye(n) - np.ones((n,n))/n)/np.sqrt(n-1) + A_n = enX[:, :n] @ P_n + Y_n = enY[:, :n] @ P_n + Gbar_n = np.mean(enG[:, :, :n], axis=-1) + + GbarCxx_n = Gbar_n @ A_n @ A_n.T + Cyx_n = Y_n @ A_n.T + + err.append(np.linalg.norm(GbarCxx_n-Cyx_n)) + Cxy_norm.append(np.linalg.norm(Cyx_n)) + GbarCxx_norm.append(np.linalg.norm(GbarCxx_n)) + + + # Plot norm of Cyx vs ensemble size + plt.style.use('seaborn-v0_8-whitegrid') + fig, ax = plt.subplots(figsize=(7, 4), dpi=140) + + ax.plot( + ens, + Cxy_norm, + color='#1f77b4', + linewidth=2.5, + marker='o', + markersize=6, + markerfacecolor='white', + markeredgewidth=1.3, + label=r'$\|C_{yx}\|$' + ) + ax.plot( + ens, + GbarCxx_norm, + color='#2ca02c', + linewidth=2.5, + marker='s', + markersize=6, + markerfacecolor='white', + markeredgewidth=1.3, + label=r'$\|\bar{G}C_{xx}\|$' + ) + ax.plot( + ens, + err, + color='#d62728', + linewidth=2.5, + linestyle='--', + marker='^', + markersize=6, + markerfacecolor='white', + markeredgewidth=1.3, + label=r'$\|\bar{G}C_{xx} - C_{yx}\|$' + ) + + ax.set_xscale('log') + ax.set_ylim(0, None) + ax.set_xlabel('Ensemble size', fontsize=11) + ax.set_ylabel(r'$L_2$-norm', fontsize=11) + + ax.grid(True, which='major', linestyle='-', linewidth=0.7, alpha=0.45) + ax.grid(True, which='minor', linestyle=':', linewidth=0.6, alpha=0.25) + ax.spines['top'].set_visible(False) + ax.spines['right'].set_visible(False) + + ax.legend(loc='best', frameon=True, fancybox=True, framealpha=0.92) + fig.tight_layout() + fig.savefig(os.path.join(folder, 'sensitivity_matrix_convergence.png'), dpi=300) + #plt.show() + + + +if __name__ == "__main__": + test_finite_difference_gradient(run_fda=True, run_adjoint=True, folder='TEST/TEMP') + #test_sens_matrix_of_log_permx(run=True, folder='TEST/SENS_WOPR_LOG_PERMX_SUBSTATES') \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index db1ea66..d0690e3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,3 +28,12 @@ dependencies = [ [project.urls] Homepage = "https://github.com/Python-Ensemble-Toolbox/SimulatorWrap" Issues = "https://github.com/Python-Ensemble-Toolbox/SimulatorWrap/ issues" + +[tool.pytest.ini_options] +# Suppress the DeprecationWarning from `multiprocess` calling os.fork() in a +# multi-threaded process (triggered by p_tqdm's parallel worker launch). +# This is a known limitation of the fork-based multiprocessing backend and +# does not affect correctness; the warning is harmless in this context. +filterwarnings = [ + "ignore::DeprecationWarning:multiprocess", +] diff --git a/src/subsurface/multphaseflow/jutul_darcy.py b/src/subsurface/multphaseflow/jutul_darcy.py index 9b0b1c2..f5e1c28 100644 --- a/src/subsurface/multphaseflow/jutul_darcy.py +++ b/src/subsurface/multphaseflow/jutul_darcy.py @@ -1,690 +1,1232 @@ -''' +""" Simulator wrapper for the JutulDarcy simulator. -This module provides a wrapper interface for running JutulDarcy simulations -with support for ensemble-based workflows and flexible output formatting. -''' +This module provides a Python interface (:class:`JutulDarcy`) for running +JutulDarcy reservoir simulations from a single configuration dictionary. +It supports: + +* Single or ensemble forward simulations (parallel via :mod:`p_tqdm`). +* Input through either a static ``.DATA`` file or a Mako-templated ``.mako`` + file rendered per ensemble member. +* Flexible result extraction (field and per-well summary keywords) with + ``list`` / ``dict`` / ``DataFrame`` output formats. +* Adjoint sensitivity computation for well-based objectives with respect to + reservoir parameters (porosity, permeability, optionally log-scaled or + copied across PERMX/Y/Z). + +The wrapper communicates with Julia via :mod:`juliacall`. Heavy Julia state +is created lazily inside worker processes so the class can be pickled for +multiprocessing. +""" import os import shutil import warnings - +import datetime as dt import numpy as np import pandas as pd -import datetime as dt +from collections import deque +from contextlib import contextmanager +from dataclasses import dataclass +from pathlib import Path +from typing import Any, Iterable from mako.template import Template from p_tqdm import p_map from tqdm import tqdm - -__author__ = "Mathias Methlie Nilsen" +__author__ = ["Mathias Methlie Nilsen"] # With help from "Claude Opus 4.7" __all__ = ["JutulDarcy"] -os.environ["PYTHON_JULIACALL_HANDLE_SIGNALS"] = "yes" -os.environ["PYTHON_JULIACALL_THREADS"] = "1" -os.environ["PYTHON_JULIACALL_OPTLEVEL"] = "3" +# ============================================================================ # +# Environment configuration +# ============================================================================ # +# These environment variables must be set BEFORE juliacall is imported. We use +# `setdefault` so users can still override them externally. +os.environ.setdefault("PYTHON_JULIACALL_HANDLE_SIGNALS", "yes") +os.environ.setdefault("PYTHON_JULIACALL_THREADS", "1") +os.environ.setdefault("PYTHON_JULIACALL_OPTLEVEL", "3") + +# Silence the noisy "juliacall module already imported" warning that appears +# in worker processes when Julia is re-imported. warnings.filterwarnings("ignore", message=".*juliacall module already imported.*") -PBAR_OPTS = { +# ============================================================================ # +# Module-level constants +# ============================================================================ # +SECONDS_PER_DAY: int = 86_400 + +#: Shared progress-bar styling used by both forward and adjoint loops. +PBAR_OPTS: dict[str, Any] = { "ncols": 110, "colour": "#285475", - "bar_format": "{desc}: {percentage:3.0f}% [{bar}] {n_fmt}/{total_fmt} │ ⏱ {elapsed}<{remaining} │ {rate_fmt}", + "bar_format": ( + "{desc}: {percentage:3.0f}% [{bar}] {n_fmt}/{total_fmt} " + "│ ⏱ {elapsed}<{remaining} │ {rate_fmt}" + ), "ascii": "-◼", } +#: Summary keywords whose values are cumulative (totals), as opposed to rates. +CUMULATIVE_KEYS: frozenset[str] = frozenset({ + "FOPT", "FGPT", "FWPT", "FWLT", "FWIT", + "WOPT", "WGPT", "WWPT", "WLPT", +}) + +VALID_OUTPUT_FORMATS: frozenset[str] = frozenset({"list", "dict", "dataframe"}) +VALID_ADJOINT_MODES: frozenset[str] = frozenset({"sensitivities", "optimization"}) +VALID_REPORT_TYPES: frozenset[str] = frozenset({"days", "dates"}) + +#: Map summary keyword -> (phase string, is_rate). Cumulative totals +#: (``WxPT``) are marked ``is_rate=False``; instantaneous rates (``WxPR``, +#: ``WxIR``) are marked ``is_rate=True``. +PHASE_MAP: dict[str, tuple[str, bool]] = { + "WOPT": ("oil", False), "WGPT": ("gas", False), + "WWPT": ("water", False), "WLPT": ("liquid", False), + "WOPR": ("oil", True), "WGPR": ("gas", True), + "WWPR": ("water", True), "WWIR": ("water", True), + "WLPR": ("liquid", True), +} -class JutulDarcy: +#: Map a phase name to the corresponding JutulDarcy rate-target type name. +RATE_ID_MAP: dict[str, str] = { + "mass": "TotalSurfaceMassRate", + "liquid": "SurfaceLiquidRateTarget", + "water": "SurfaceWaterRateTarget", + "oil": "SurfaceOilRateTarget", + "gas": "SurfaceGasRateTarget", + "rate": "TotalRateTarget", +} - def __init__(self, options: dict): - """ - JutulDarcy simulation wrapper. +#: Display units for known summary keywords (metric system). +UNIT_MAP: dict[str, str] = { + "PORO": "", "PERMX": "mD", "PERMY": "mD", "PERMZ": "mD", + "FOPT": "Sm3", "FGPT": "Sm3", "FWPT": "Sm3", "FWLT": "Sm3", "FWIT": "Sm3", + "FOPR": "Sm3/day", "FGPR": "Sm3/day", "FWPR": "Sm3/day", + "FLPR": "Sm3/day", "FWIR": "Sm3/day", + "WOPR": "Sm3/day", "WGPR": "Sm3/day", "WWPR": "Sm3/day", + "WLPR": "Sm3/day", "WWIR": "Sm3/day", +} - Parameters - ---------- - options : dict - Configuration dictionary controlling input files, reporting, output - formatting, parallel execution, and optional adjoint computation. - Supported keys: +#: Permeability keywords in canonical (x, y, z) order. +PERM_KEYS: tuple[str, ...] = ("PERMX", "PERMY", "PERMZ") + +#: Map lower-case permeability parameter substring -> axis index in PERM_KEYS. +PERM_INDEX: dict[str, int] = {"permx": 0, "permy": 1, "permz": 2} + + +# ============================================================================ # +# Configuration dataclasses +# ============================================================================ # +@dataclass +class AdjointObjective: + """ + One adjoint objective specification (one well + phase). + + Attributes + ---------- + wellID : str + Name of the well associated with this objective. + phase : str + Phase string used to pick the rate target (see :data:`RATE_ID_MAP`). + is_rate : bool + True for instantaneous rate objectives, False for cumulative totals. + parameters : list[str] + Parameters w.r.t. which the gradient should be evaluated, e.g. + ``["PORO", "PERMX"]`` (case-insensitive; ``"log"`` enables log-scaling). + steps : Any + Either the string ``"all"`` (use the wrapper's global report points) + or an explicit list of ``int`` days or :class:`datetime.datetime` + objects at which to evaluate the objective. + """ + wellID: str + phase: str + is_rate: bool + parameters: list[str] + steps: Any # 'all' | list[int] | list[datetime] + + +# ============================================================================ # +# Helper functions (pure / stateless) +# ============================================================================ # +def get_metric_unit(key: str) -> str: + """ + Return the metric-system unit string for a summary keyword. + + Parameters + ---------- + key : str + Summary keyword (case-insensitive), e.g. ``"FOPT"``. + + Returns + ------- + str + Unit string from :data:`UNIT_MAP`, or ``"Unknown"`` if not found. + """ + return UNIT_MAP.get(key.upper(), "Unknown") + + +def _process_datatype_info(datatypes: Iterable[str]) -> list[str]: + """ + Expand grouped datatype specs into one entry per well. + + A spec like ``"WOPR:W1:W2"`` is expanded to ``["WOPR:W1", "WOPR:W2"]``. + Specs without a colon (field-level) are passed through unchanged. + + Parameters + ---------- + datatypes : iterable of str + User-supplied datatype identifiers. + + Returns + ------- + list of str + Normalised, one-well-per-entry list. + """ + out: list[str] = [] + for d in datatypes: + if ":" in d: + base, *wells = d.split(":") + out.extend(f"{base}:{w}" for w in wells) + else: + out.append(d) + return out + + +def _process_adjoint_info(adjoint_info: dict) -> dict[str, AdjointObjective]: + """ + Normalise the ``adjoints`` config dict into :class:`AdjointObjective` items. + + Each input entry may specify a single well or a list of wells; the output + contains one :class:`AdjointObjective` per (datatype, well) pair, keyed by + ``":"``. + + Parameters + ---------- + adjoint_info : dict + Raw mapping ``{datatype: {"wellID": ..., "parameters": ..., "steps": ...}}``. + + Returns + ------- + dict[str, AdjointObjective] + Flattened objective specifications. + """ + info: dict[str, AdjointObjective] = {} + for dataID, spec in adjoint_info.items(): + # Normalise scalar values to lists for uniform iteration. + wells = spec["wellID"] + if isinstance(wells, str): + wells = [wells] + params = spec["parameters"] + if isinstance(params, str): + params = [params] + + phase, is_rate = PHASE_MAP[dataID] + for w in wells: + info[f"{dataID}:{w}"] = AdjointObjective( + wellID=w, phase=phase, is_rate=is_rate, + parameters=list(params), steps=spec["steps"], + ) + return info - - ``runfile`` : str, optional - Path to either a ``.mako`` template or a ``.DATA`` run file. - - ``reporttype`` : {"days", "dates"}, optional - Type of report index. Default is ``"days"``. +def _active_to_full_grid(vec: np.ndarray, actnum_vec: np.ndarray, + fill_value: float = 0.0) -> np.ndarray: + """ + Embed an active-cell vector into the full grid layout. + + Parameters + ---------- + vec : np.ndarray + Either a vector defined only on active cells (length ``actnum_vec.sum()``) + or already on the full grid (length ``len(actnum_vec)``). + actnum_vec : np.ndarray + Flat ACTNUM array (1 = active, 0 = inactive), Fortran-ordered. + fill_value : float, optional + Value placed at inactive cells. Defaults to ``0.0``. + + Returns + ------- + np.ndarray + Full-grid vector of length ``len(actnum_vec)`` and dtype ``float64``. + + Raises + ------ + ValueError + If ``vec`` matches neither the active nor the full grid size. + """ + n_active = int(actnum_vec.sum()) + if len(vec) == n_active: + full = np.full(actnum_vec.shape, fill_value, dtype=np.float64, order="F") + full[actnum_vec == 1] = vec + return full + if len(vec) == len(actnum_vec): + return np.asarray(vec, dtype=np.float64) + raise ValueError("Parameter length does not match number of active cells") - - ``reportpoint`` : list, optional - Report points used for output indexing. For ``"days"`` this is - a list of numeric day values; for ``"dates"`` a list of - ``datetime`` objects. - - ``datatype`` : list[str], optional - Result keywords to extract. Supports field keys (for example - ``"FOPT"``) and well keys (for example ``"WOPR:PROD1"``). +@contextmanager +def _chdir(path: str | os.PathLike): + """ + Context manager that temporarily changes the working directory. + + The previous working directory is always restored, including when the + wrapped block raises an exception. + + Parameters + ---------- + path : str or path-like + Directory to switch into. + + Yields + ------ + None + """ + prev = os.getcwd() + os.chdir(path) + try: + yield + finally: + os.chdir(prev) + + +def _suppress_julia(julia, code: str): + """ + Evaluate a Julia expression with ``stdout`` and ``stderr`` redirected to devnull. + + Used to silence verbose JutulDarcy setup/solver messages. + + Parameters + ---------- + julia : juliacall.Main + The Julia main module reference. + code : str + A single Julia expression (will be wrapped, not multi-line statements). + + Returns + ------- + Any + Result of the inner Julia expression. + """ + return julia.seval(f""" + redirect_stdout(devnull) do + redirect_stderr(devnull) do + {code} + end + end + """) + + +def _get_mapping_value(obj, key_name: str, julia): + """ + Look up ``key_name`` in a Julia mapping using both Symbol and string keys. + + JutulDarcy mappings sometimes use ``Symbol`` keys and sometimes plain + strings; this helper tries both so callers don't have to. + + Parameters + ---------- + obj : Any + A Julia mapping-like object supporting ``haskey`` / ``getindex``. + key_name : str + Candidate key as a Python string. + julia : juliacall.Main + Julia main module (used to construct ``Symbol``). + + Returns + ------- + Any or None + The value if found; ``None`` otherwise. + """ + for candidate in (julia.Symbol(key_name), key_name): + try: + if julia.haskey(obj, candidate): + return obj[candidate] + except Exception: + # `haskey` may not be defined for every object type encountered. + continue + return None + + +def _extract_key_value(root_object, keys, julia): + """ + Breadth-first search for the first matching key in nested Julia mappings. + + Useful because JutulDarcy gradient structures can nest gradients several + layers deep (e.g. ``grad[:model][:reservoir][:porosity]``). + + Parameters + ---------- + root_object : Any or list + Starting Julia object(s) to search. + keys : str or list of str + Candidate key name(s). The first key found anywhere in the tree wins. + julia : juliacall.Main + Julia main module. + + Returns + ------- + Any or None + Value associated with the first matching key, or ``None``. + """ + if not isinstance(keys, list): + keys = [keys] + if not isinstance(root_object, list): + root_object = [root_object] + + queue = deque(root_object) + visited: set[int] = set() # avoid revisiting the same Julia object + + while queue: + cur = queue.popleft() + oid = id(cur) + if oid in visited: + continue + visited.add(oid) + + # Try to match at the current level first. + for k in keys: + val = _get_mapping_value(cur, k, julia) + if val is not None: + return val + + # Otherwise, enqueue children for further exploration. + try: + for ck in list(julia.keys(cur)): + try: + queue.append(cur[ck]) + except Exception: + continue + except Exception: + # `cur` is a leaf (no `keys` method); skip it. + continue + + return None + + +def _extract_adjoint(jlgrad, jlcase, parameter: str, actnum: np.ndarray, + perm_copied: bool, julia) -> np.ndarray: + """ + Extract and post-process an adjoint gradient for a single parameter. + + Handles unit scaling (mD → SI), optional log-scaling, embedding back onto + the full grid, and the "copied permeability" convention where PERMY/PERMZ + are duplicates of PERMX in the input deck. + + Parameters + ---------- + jlgrad : Any + Julia gradient object returned by the adjoint solver. + jlcase : Any + Julia case object (needed to read PERMX/Y/Z values when log-scaling). + parameter : str + Parameter name. May contain ``"poro"``, ``"perm"``, ``"permx"``, + ``"permy"``, ``"permz"`` (case-insensitive). Including ``"log"`` + enables log-scaling for permeability. + actnum : np.ndarray + Flat ACTNUM vector (1/0) for embedding into the full grid. + perm_copied : bool + If ``True``, treat PERMX as the master and sum contributions from all + three permeability directions into the returned gradient. + julia : juliacall.Main + Julia main module. + + Returns + ------- + np.ndarray + Full-grid gradient vector. + + Raises + ------ + ValueError + If the gradient could not be located or the parameter is unsupported. + """ + p = parameter.lower() + + # ---- Porosity -------------------------------------------------------- # + if "poro" in p: + grad = _extract_key_value(jlgrad, "porosity", julia) + if grad is None: + raise ValueError(f"Could not find porosity gradient for '{parameter}'") + return _active_to_full_grid(np.asarray(grad), actnum) + + # ---- Permeability ---------------------------------------------------- # + if "perm" in p: + log_scale = "log" in p + grad = _extract_key_value(jlgrad, "permeability", julia) + if grad is None: + raise ValueError(f"Could not find permeability gradient for '{parameter}'") + + full = np.asarray(grad) # shape: (3, n_active) + mdarcy = julia.seval("si_unit(:milli)*si_unit(:darcy)") # mD → SI factor + + def _scale(adj: np.ndarray, axis: int) -> np.ndarray: + """Apply log-scale (∂J/∂log(k) = k·∂J/∂k) or unit conversion.""" + if log_scale: + perm = np.array(jlcase.input_data["GRID"][PERM_KEYS[axis]]) + return adj * perm.flatten(order="F") + return adj * mdarcy + + # Copied-permeability case: aggregate gradients from all 3 axes. + if perm_copied: + out = np.zeros(actnum.shape, dtype=np.float64) + for i in range(3): + adj = _active_to_full_grid(full[i], actnum) + out += _scale(adj, i) + return out + + # Otherwise pick the axis encoded in the parameter name. + idx = next((v for k, v in PERM_INDEX.items() if k in p), None) + if idx is None: + raise ValueError(f"Adjoint not implemented for '{parameter}'") + return _scale(_active_to_full_grid(full[idx], actnum), idx) + + raise ValueError(f"Adjoint not implemented for parameter '{parameter}'") + + +def well_QOI_objective(wellID: str, phaseID: str, time: Iterable[float], + step_index=None, is_rate: bool = True, julia=None): + """ + Build per-timestep Julia QOI closures for a single well/phase. + + Each closure returns either a daily-volume contribution (``is_rate=True``, + one time only) or the cumulative integrand over time up to a horizon + (``is_rate=False``). The sign is flipped automatically for producers so + objectives are positive for "good" outcomes (e.g. produced oil). + + Parameters + ---------- + wellID : str + Name of the well as it appears in the case. + phaseID : str + Phase identifier (see :data:`RATE_ID_MAP`). + time : iterable of float + Sequence of horizon times in seconds. + step_index : Any, optional + Currently unused; retained for backwards compatibility. + is_rate : bool, optional + Selects rate (instantaneous) vs cumulative objective formulation. + julia : juliacall.Main, optional + Existing Julia module; one is imported on demand if omitted. + + Returns + ------- + list[Any] + List of Julia closures, one per entry in ``time``. + + Raises + ------ + ValueError + If ``phaseID`` is not in :data:`RATE_ID_MAP`. + """ + if julia is None: + from juliacall import Main as julia + julia.seval("using JutulDarcy") - - ``adjoints`` : dict, optional - Objective and parameter configuration for adjoint sensitivities. + if phaseID not in RATE_ID_MAP: + raise ValueError(f"Unknown rate type: {phaseID}") + rate_sym = RATE_ID_MAP[phaseID] - - ``output_format`` : {"list", "dict", "dataframe"}, optional - Output representation for forward results. Default is ``"list"``. - (will probably be dataframe in the future) + qois = [] + for i, sec in enumerate(time): + if is_rate: + # Instantaneous rate: contribute only at the exact step matching `sec`. + obj = julia.seval( + f""" + function well_QOI_{i}(model, state, dt, step_info, forces) + if step_info[:time] != {sec} + return 0.0 + end + ctrl = forces[:Facility].control[Symbol("{wellID}")] + sign = ctrl isa JutulDarcy.ProducerControl ? -1.0 : 1.0 + rate = JutulDarcy.compute_well_qoi(model, state, forces, Symbol("{wellID}"), {rate_sym}) + return sign * rate * si_unit(:day) # rate is in Sm3/s: convert to Sm3/day + end + """ + ) + else: + # Cumulative: integrate (dt * rate) up to the horizon `sec`. + obj = julia.seval( + f""" + function well_QOI_{i}(model, state, dt, step_info, forces) + if step_info[:time] > {sec} + return 0.0 + end + ctrl = forces[:Facility].control[Symbol("{wellID}")] + sign = ctrl isa JutulDarcy.ProducerControl ? -1.0 : 1.0 + rate = JutulDarcy.compute_well_qoi(model, state, forces, Symbol("{wellID}"), {rate_sym}) + return sign * dt * rate + end + """ + ) + qois.append(obj) + return qois - - ``adjoint_pbar`` : bool, optional - If ``True``, display progress bars during adjoint solves. - Not recommended for large ensembles due to many progress bars. - Default is ``True``. - - ``parallel`` : int, optional - Number of processes for ensemble simulations. Default is ``1``. - """ - # Check runfile - runfile = options.get('runfile') - if runfile: - self.makofile = runfile if runfile.endswith('.mako') else None - self.datafile = runfile if runfile.endswith('.DATA') else None - - # Report Options - self.report_type = options.get('reporttype', 'days') # days or dates - self.report = options.get('reportpoint', None) # list of days or dates (dates as datetime objects) - self.index = [self.report_type, self.report] - - # Process datatypes - datatype = options.get('datatype', ['FOPT', 'FGPT', 'FWPT', 'FWIT']) - self.datatype = _process_datatype_info(datatype) - - # Adjoint information - if 'adjoints' in options: +# ============================================================================ # +# Main wrapper +# ============================================================================ # +class JutulDarcy: + """ + Python wrapper around a JutulDarcy reservoir simulation. + + The class is instantiated once with an options dictionary describing the + simulation setup, then called like a function on either a single input + (``dict`` or ``.DATA`` path) or a list of inputs (one per ensemble member). + + See the module docstring for the list of supported options. + + Parameters + ---------- + options : dict + Configuration dictionary. Recognised keys: + + - ``runfile`` : str + Path to a ``.mako`` template or ``.DATA`` file. + - ``reporttype`` : {"days", "dates"} + How report points are interpreted. Default ``"days"``. + - ``reportpoint`` : list + Report points (numeric days or :class:`datetime.datetime`). + - ``datatype`` : list[str] + Summary keywords to extract. Default + ``["FOPT", "FGPT", "FWPT", "FWIT"]``. + - ``adjoints`` : dict + Objective/parameter spec; enables adjoint computation. + - ``output_format`` : {"list", "dict", "dataframe"} + Output container type. Default ``"dataframe"``. + - ``adjoint_pbar`` : bool + Show a per-objective progress bar. Default ``False``. + - ``parallel`` : int + Number of worker processes. Default ``1``. + - ``perm_copied`` : bool + See :func:`_extract_adjoint`. Default ``False``. + - ``adjoint_mode`` : {"sensitivities", "optimization"} + Selects the JutulDarcy gradient pathway. Default + ``"sensitivities"``. + - ``optimization_targets`` : Any + Reserved for future use. + - ``eval_adjoint_funcs`` : bool + If True, store objective function values in ``self.adjoint_funcs``. + + Raises + ------ + ValueError + If ``reporttype``, ``output_format`` or ``adjoint_mode`` is invalid. + """ + + def __init__(self, options: dict): + # ---- Runfile (template vs static deck) --------------------------- # + runfile = options.get("runfile") + self.makofile = runfile if runfile and runfile.endswith(".mako") else None + self.datafile = runfile if runfile and runfile.endswith(".DATA") else None + + # ---- Report-point configuration ---------------------------------- # + self.report_type = options.get("reporttype", "days") + if self.report_type not in VALID_REPORT_TYPES: + raise ValueError( + f"Invalid reporttype '{self.report_type}'. " + f"Must be one of {sorted(VALID_REPORT_TYPES)}" + ) + self.report = options.get("reportpoint") + self.index: list = [self.report_type, self.report] + + # Computed lazily during the first forward run. + self.report_seconds: np.ndarray | None = None + self.start_date: dt.datetime | None = None + + # ---- Datatypes to extract ---------------------------------------- # + self.datatype = _process_datatype_info( + options.get("datatype", ["FOPT", "FGPT", "FWPT", "FWIT"]) + ) + + # ---- Adjoint configuration --------------------------------------- # + if "adjoints" in options: self.compute_adjoints = True - self.adjoint_info = _process_adjoint_info(options['adjoints']) - self.eval_adjoint_funcs = options.get('eval_adjoint_funcs', True) - self.adjoint_funcs = None + self.adjoint_info = _process_adjoint_info(options["adjoints"]) + self.eval_adjoint_funcs = options.get("eval_adjoint_funcs", False) else: self.compute_adjoints = False self.eval_adjoint_funcs = False - self.adjoint_funcs = None - - # Other options - self.output_format = options.get('output_format', 'dataframe') # list, dict or dataframe - self.adjoint_pbar = options.get('adjoint_pbar', True) - self.parallel = options.get('parallel', 1) + self.adjoint_funcs: pd.DataFrame | None = None + + # ---- Output / execution options ---------------------------------- # + self.output_format = options.get("output_format", "dataframe") + if self.output_format not in VALID_OUTPUT_FORMATS: + raise ValueError( + f"Invalid output_format '{self.output_format}'. " + f"Must be one of {sorted(VALID_OUTPUT_FORMATS)}" + ) + self.adjoint_pbar = options.get("adjoint_pbar", False) + self.parallel = int(options.get("parallel", 1)) + self.perm_copied = options.get("perm_copied", False) + self.adjoint_mode = options.get("adjoint_mode", "sensitivities") + if self.adjoint_mode not in VALID_ADJOINT_MODES: + raise ValueError( + f"Invalid adjoint_mode '{self.adjoint_mode}'. " + f"Must be one of {sorted(VALID_ADJOINT_MODES)}" + ) + self.optimization_targets = options.get("optimization_targets") - # This is for PET to work properly + # ---- PET compatibility attributes -------------------------------- # + # `input_dict` and `true_order` are consumed by the surrounding PET + # framework; we keep them as direct passthroughs. self.input_dict = options self.true_order = self.index - - def __call__(self, inputs: list[dict]|dict|str): + # ---------------------------------------------------------------- # + # Public entry point + # ---------------------------------------------------------------- # + def __call__(self, inputs: list[dict] | dict | str): """ - Execute parallel forward simulations for all ensemble members. - - This method runs the forward simulation for all input parameter sets in parallel, - optionally computing adjoint sensitivities if configured. + Run the configured simulation for one or many ensemble members. Parameters ---------- - inputs : list of dict, dict, or str - List containing input parameter sets, indexed by ensemble member, or a single input dictionary, - or a string representing a file path to input datafile. - - number (0, 1, 2, ...). Each value is a dict of parameters for that member. + inputs : list[dict] or dict or str + Either a list of inputs (one per ensemble member), a single dict + of Mako template parameters, or a path to a ``.DATA`` deck. Returns ------- - output : list of (dict, list, or pd.DataFrame) - Forward simulation results for all ensemble members, formatted according - to the 'out_format' option specified during initialization. - adjoints : pd.DataFrame, optional - Adjoint sensitivities for all objectives and parameters. Only returned if - 'adjoints' configuration was provided during initialization. Contains - multi-indexed columns (objective, parameter) and index based on reporttype. - - Notes - ----- - - Existing simulation folders (En_*) are deleted before execution to prevent - conflicts from previous runs. - - Parallel execution uses the number of CPUs specified in the 'parallel' option. - - Progress is displayed via a progress bar during execution. - - Examples - -------- - >>> simulator = JutulDarcyWrapper(options) - >>> inputs = [{'param1': [1, 2, 3]}, {'param1': [1.1, 2.1, 3.1]}] - >>> results = simulator(inputs) + Forward-only mode + Single input → one result; list input → list of results. The + result type follows ``self.output_format``. + Adjoint mode + Same shape as above, but each result is a ``(forward, adjoint)`` + tuple where ``adjoint`` is a :class:`pandas.DataFrame` with a + ``(objective, parameter)`` MultiIndex on the columns. """ + # Normalise to a list so the parallel path is uniform. if isinstance(inputs, (dict, str)): inputs = [inputs] - # Delet all existing En_* folders - for item in os.listdir('.'): - if os.path.isdir(item) and item.startswith('En_'): - shutil.rmtree(item) - - # simulate all inputs in parallel + # Clean up any stale simulation folders from a previous (possibly failed) run. + self._cleanup_simulation_folders() + + n = len(inputs) outputs = p_map( - self.run_fwd_sim, - [inputs[n] for n in range(len(inputs))], - list(range(len(inputs))), + self.run_fwd_sim, + list(inputs), + list(range(n)), num_cpus=self.parallel, - unit='sim', - desc='Simulations', + unit="sim", + desc="Simulations", leave=False, - **PBAR_OPTS + **PBAR_OPTS, ) + # In adjoint mode each worker returns a (result, adjoint) tuple; split + # them out so callers see two parallel collections rather than a list + # of tuples. if self.compute_adjoints: results, adjoints = zip(*outputs) - if len(inputs) == 1: - results = results[0] - adjoints = adjoints[0] + if n == 1: + return results[0], adjoints[0] + return list(results), list(adjoints) - return results, adjoints - else: - return outputs - + return outputs[0] if n == 1 else outputs + + # ---------------------------------------------------------------- # + # Per-member forward simulation + # ---------------------------------------------------------------- # + def run_fwd_sim(self, input: dict | str, idn: int = 0, + delete_folder: bool = True): + """ + Run a single forward simulation (and optionally adjoints) in isolation. - def run_fwd_sim(self, input: dict|str, idn: int=0, delete_folder: bool=True): + Each call creates and operates inside its own ``En_`` folder so + multiple workers cannot collide on intermediate files written by + JutulDarcy. + + Parameters + ---------- + input : dict or str + Mako template parameters or a path to a ``.DATA`` deck. + idn : int, optional + Ensemble member index (used for folder naming and pbar position). + delete_folder : bool, optional + If True (default), remove the ``En_`` folder when done. - # Import Julia and JutulDarcy (this needs to be done here for multiprocessing to work properly) + Returns + ------- + Same as :meth:`__call__` for a single member. + """ + # Julia is imported lazily inside the worker so it picks up the + # process-local thread/state. Importing at module scope would break + # multiprocessing's fork/spawn semantics. from juliacall import Main as julia - julia.seval('using JutulDarcy, Jutul') - - # Make simulation folder - folder = f'En_{idn}' - os.makedirs(folder) + julia.seval("using JutulDarcy, Jutul") + + folder = Path(f"En_{idn}") + folder.mkdir(exist_ok=False) + + try: + # Stage the deck (render template or copy file). + datafile = self._stage_input(input, folder, idn) + + # All Julia file I/O must happen relative to the deck location. + with _chdir(folder): + case = self._setup_case(datafile, julia) + julia.case = case # expose to Julia-side eval'd expressions + + units = self._detect_units(case, julia) + actnum_vec = self._extract_actnum(case) + + # Forward solve. `output_substates=True` keeps intermediate + # states needed for adjoint reconstruction. + jlres = julia.simulate_reservoir( + case, info_level=-1, output_substates=True + ) + julia.res = jlres + + pyres = self.extract_datatypes(jlres, case, units, julia) + output = self._format_output(pyres) + + if self.compute_adjoints: + adjoints = self._compute_adjoints( + case, jlres, pyres, units, actnum_vec, idn, julia + ) + return output, adjoints + + return output + finally: + # Always clean up, even on failure, to keep the workspace tidy. + if delete_folder and folder.exists(): + shutil.rmtree(folder, ignore_errors=True) + + # ---------------------------------------------------------------- # + # Setup helpers + # ---------------------------------------------------------------- # + @staticmethod + def _cleanup_simulation_folders() -> None: + """Remove any ``En_*`` directories left in the current directory.""" + for item in os.listdir("."): + if item.startswith("En_") and os.path.isdir(item): + shutil.rmtree(item, ignore_errors=True) + + def _stage_input(self, input: dict | str, folder: Path, idn: int) -> str: + """ + Render a Mako template or copy a static ``.DATA`` file into ``folder``. + + Parameters + ---------- + input : dict or str + Either Mako parameters or a path to a ``.DATA`` deck. + folder : Path + Destination directory. + idn : int + Ensemble index (injected into the Mako context as ``member``). - # Check input type (datafile or input for makofile) + Returns + ------- + str + Basename of the deck file inside ``folder``. + + Raises + ------ + FileNotFoundError + If a string input does not point to an existing file. + ValueError + If a string input is not a ``.DATA`` file. + TypeError + For unsupported input types. + """ if isinstance(input, dict): - input['member'] = idn - datafile = self.render_makofile(self.makofile, folder, input) - elif isinstance(input, str): - assert os.path.isfile(input), f'Input file {input} not found' - assert input.endswith('.DATA'), 'Input string must be a path to a .DATA file' - - # Copy datafile to simulation folder + input["member"] = idn + return self.render_makofile(self.makofile, str(folder), input) + + if isinstance(input, str): + if not os.path.isfile(input): + raise FileNotFoundError(f"Input file {input} not found") + if not input.endswith(".DATA"): + raise ValueError("Input string must be a path to a .DATA file") datafile = os.path.basename(input) shutil.copy(input, folder) + return datafile - # Enter simulation folder and run simulation - os.chdir(folder) + raise TypeError(f"Unsupported input type: {type(input).__name__}") - # Setup case from datafile - suppress = True - if suppress: - case = julia.seval(f""" - redirect_stdout(devnull) do - redirect_stderr(devnull) do - setup_case_from_data_file("{datafile}") - end - end - """) - else: - case = julia.setup_case_from_data_file(datafile) - - # Get Units - if julia.haskey(case.input_data["RUNSPEC"], "METRIC"): - units = 'metric' - elif julia.haskey(case.input_data["RUNSPEC"], "SI"): - units = 'si' - elif julia.haskey(case.input_data["RUNSPEC"], "FIELD"): - units = 'field' - else: - units = julia.missing + @staticmethod + def _setup_case(datafile: str, julia): + """Invoke JutulDarcy's case-setup routine with output suppressed.""" + return _suppress_julia(julia, f'setup_case_from_data_file("{datafile}")') + + @staticmethod + def _detect_units(case, julia) -> str | Any: + """ + Detect the unit system from the deck's RUNSPEC section. - # Get some grid info + Returns + ------- + str + One of ``"metric"``, ``"si"``, ``"field"`` if specified. + Any + ``julia.missing`` if no unit keyword is present. + """ + runspec = case.input_data["RUNSPEC"] + for key, name in (("METRIC", "metric"), ("SI", "si"), ("FIELD", "field")): + if julia.haskey(runspec, key): + return name + return julia.missing + + @staticmethod + def _extract_actnum(case) -> np.ndarray: + """ + Return the flat ACTNUM array (1/0), defaulting to all-active. + + Parameters + ---------- + case : Any + Julia case object. + + Returns + ------- + np.ndarray + Flat (Fortran-ordered) ACTNUM vector. + """ nx, ny, nz = case.input_data["GRID"]["cartDims"] try: - actnum = np.array(case.input_data["GRID"]["ACTNUM"]) # Shape (nx, ny, nz) - actnum_vec = actnum.flatten(order='F') # Fortran order flattening - except: - actnum_vec = np.ones(nx*ny*nz) - - # Simulate and get results - jlres = julia.simulate_reservoir(case, info_level=-1, output_substates=True) - pyres, daysIDX = self.extract_datatypes(jlres, case, units, julia) - - # Convert output to requested format - if not self.output_format == 'dataframe': - if self.output_format == 'dict': - output = pyres.to_dict(orient='list') - elif self.output_format == 'list': - output = [] - for _, row in pyres.iterrows(): - row_dict = row.to_dict() - row_dict = {k: np.atleast_1d(v) for k, v in row_dict.items()} - output.append(row_dict) - else: - output = pyres + actnum = np.array(case.input_data["GRID"]["ACTNUM"]) + return actnum.flatten(order="F") + except (KeyError, AttributeError): + # No ACTNUM keyword in the deck => every cell is active. + return np.ones(nx * ny * nz) + + def _format_output(self, pyres: pd.DataFrame): + """Convert the per-member DataFrame to the user-requested container.""" + if self.output_format == "dataframe": + return pyres + if self.output_format == "dict": + return pyres.to_dict(orient="list") + return pyres.to_dict(orient="records") # 'list' format → list of records + + @staticmethod + def render_makofile(makofile: str, folder: str, input: dict) -> str: + """ + Render a Mako template into a ``.DATA`` file in ``folder``. + Parameters + ---------- + makofile : str + Path to the ``.mako`` template. + folder : str + Destination directory. + input : dict + Template context variables. - # ---------------------------------------------------------------------------------------------- - # Compute adjoints - # ---------------------------------------------------------------------------------------------- - if self.compute_adjoints: + Returns + ------- + str + Basename of the rendered ``.DATA`` file. + """ + datafile = os.path.basename(makofile).replace(".mako", ".DATA") + template = Template(filename=makofile) + with open(os.path.join(folder, datafile), "w") as f: + f.write(template.render(**input)) + return datafile - # Initialize adjoint results storage - info = self.adjoint_info - columns = info.keys() - grad_dict = {(col, param): [] for col in columns for param in info[col]['parameters']} + # ---------------------------------------------------------------- # + # Datatype extraction + # ---------------------------------------------------------------- # + def extract_datatypes(self, jlres, jlcase, units, julia) -> pd.DataFrame: + """ + Extract requested summary keywords from a finished simulation. + + Parameters + ---------- + jlres : Any + Julia result object from ``simulate_reservoir``. + jlcase : Any + Julia case object (used for the start date). + units : str or Any + Unit-system identifier (see :meth:`_detect_units`). + julia : juliacall.Main + Julia main module. - if self.eval_adjoint_funcs: - func_dict = {col: [] for col in columns} # For storing functions for each objective + Returns + ------- + pd.DataFrame + One row per report point; one column per datatype. ``df.attrs`` + holds per-column unit strings. + + Raises + ------ + ValueError + If a requested datatype or well is missing from the results. + """ + jl_units = julia.Symbol(units) if isinstance(units, str) else units + smry = julia.JutulDarcy.summary_result(jlcase, jlres, jl_units) - # Setup progress bar (iterate over adjoint objectives, not DataFrame columns) - if self.adjoint_pbar: - PBAR_OPTS.pop('colour', None) - pbar = tqdm( - self.adjoint_info.items(), - desc=f'Solving adjoints', - position=idn+1, - leave=False, - unit='obj', - dynamic_ncols=True, - colour="#713996", - **PBAR_OPTS - ) + # JutulDarcy reports time in integer seconds; match against our + # requested report points by intersecting the two integer arrays. + sim_seconds = np.array(list(smry["TIME"].seconds), dtype=np.int64) + self.start_date = jlcase.input_data["RUNSPEC"]["START"] + self.report_seconds = self._compute_report_seconds() + + idx = np.flatnonzero(np.isin(sim_seconds, self.report_seconds)) + + res: dict[str, np.ndarray] = {} + attrs: dict[str, str] = {} + + wells = smry["VALUES"]["WELLS"] + field = smry["VALUES"]["FIELD"] + + for datatype in self.datatype: + if ":" in datatype: + # Well-level datatype: ":" + baseID, wellID = datatype.split(":") + if wellID not in wells: + raise ValueError( + f"Well ID '{wellID}' not found for datatype '{baseID}'" + ) + well_data = wells[wellID] + if baseID not in well_data: + raise ValueError( + f"Datatype '{baseID}' not found for well '{wellID}'" + ) + res[datatype] = np.array(well_data[baseID])[idx] + attrs[datatype] = get_metric_unit(baseID) else: - pbar = self.adjoint_info.items() - - - # Loop over adjoint objectives - for col, info in pbar: - - if info['steps'] == 'all': # If 'all', use same steps as forward simulation results - stepIDX = daysIDX - elif info['steps'] == 'acc': - stepIDX = None - else: - smry = julia.JutulDarcy.summary_result(case, jlres, units) - sim_days = np.array(list(smry["TIME"].seconds)) / (24*60*60) - sim_days = sim_days.astype(int) - - if isinstance(info['steps'][0], int): - if not np.all(np.isin(info['steps'], sim_days)): - raise ValueError(f'Steps {info["steps"]} not found in simulation results for objective {col}, Available steps: {sim_days}') - stepIDX = np.argwhere(np.isin(sim_days, info['steps'])).flatten() - - elif isinstance(info['steps'][0], dt.datetime): - start_date = case.input_data["RUNSPEC"]['START'] - sim_dates = np.array([start_date + dt.timedelta(days=int(d)) for d in sim_days]) - if not np.all(np.isin(info['steps'], sim_dates)): - raise ValueError(f'Dates {info["steps"]} not found in simulation results for objective {col}, Available dates: {sim_dates}') - stepIDX = np.argwhere(np.isin(sim_dates, info['steps'])).flatten() - - - # Get QOI function for this objective - funcs = well_QOI_objective( - wellID=info['wellID'], - phaseID=info['phase'], - stepID=stepIDX, - time=np.array(jlres.time), - rate=info['is_rate'], - julia=julia - ) + # Field-level datatype. + if datatype not in field: + raise ValueError( + f"Datatype '{datatype}' not found in field results" + ) + res[datatype] = np.array(field[datatype])[idx] + attrs[datatype] = get_metric_unit(datatype) + + df = pd.DataFrame(res, index=self.index[1]) + df.index.name = self.index[0] + df.attrs = attrs + return df + + def _compute_report_seconds(self) -> np.ndarray: + """ + Convert the configured report points to integer seconds. - # Update progress bar - if self.adjoint_pbar: - update_desc = f'Adjoints for {col}' - pbar.set_description_str(update_desc) + Returns + ------- + np.ndarray + Seconds elapsed from the simulation start for each report point. + + Raises + ------ + ValueError + If ``self.report_type`` is not recognised (defensive; should be + blocked by :meth:`__init__` validation). + """ + rtype, rpoints = self.index + if rtype == "days": + return np.array(rpoints, dtype=np.int64) * SECONDS_PER_DAY + if rtype == "dates": + return np.array( + [(d - self.start_date).total_seconds() for d in rpoints], + dtype=np.int64, + ) + raise ValueError(f"Invalid report type: {rtype}") - # Comute adjoint sensitivities - julia.case = case - julia.res = jlres - for i, func in enumerate(funcs): - julia.func = func - grad = julia.seval(""" - redirect_stdout(devnull) do - redirect_stderr(devnull) do - JutulDarcy.reservoir_sensitivities( - case, res, func, - include_parameters=true, - ) - end - end - """) - - # Evaluate functions - if self.eval_adjoint_funcs: - func_val = julia.Jutul.evaluate_objective(func, case, jlres.result) - func_val = func_val*julia.seval('si_unit(:day)') if info['is_rate'] else func_val - func_dict[col].append(func_val) - - - # Extract parameter sensitivities for this objective and store in dict - for param in info['parameters']: - grad_param = _extract_adjoint( - grad, - case, - param, - actnum_vec, - info['is_rate'], - julia - ) - grad_dict[(col, param)].append(grad_param) + # ---------------------------------------------------------------- # + # Adjoint computation + # ---------------------------------------------------------------- # + def _compute_adjoints(self, case, jlres, pyres, units, actnum_vec, + idn: int, julia) -> pd.DataFrame: + """ + Compute gradients for all configured adjoint objectives. + + Parameters + ---------- + case : Any + Julia case object. + jlres : Any + Julia simulation result. + pyres : pd.DataFrame + Forward-extracted results, used as a sanity check against + Jutul's own objective evaluation. + units : Any + Unit system identifier (stored in returned ``df.attrs``). + actnum_vec : np.ndarray + Flat ACTNUM vector for full-grid embedding. + idn : int + Ensemble index (used only for pbar positioning). + julia : juliacall.Main + Julia main module. + + Returns + ------- + pd.DataFrame + Index: requested adjoint evaluation points; + columns: ``(objective, parameter)`` MultiIndex. + """ + # Optimization mode requires setting up a parameter dictionary and + # freeing the parameters Jutul will differentiate w.r.t. + if self.adjoint_mode == "optimization": + julia.grad_case = julia.seval( + "JutulDarcy.setup_reservoir_dict_optimization(case, verbose=false)" + ) + julia.seval("JutulDarcy.free_optimization_parameters!(grad_case)") + + # Pre-allocate output containers. + grad_dict: dict[tuple[str, str], list] = { + (col, p): [] + for col, obj in self.adjoint_info.items() + for p in obj.parameters + } + func_dict: dict[str, list] = ( + {col: [] for col in self.adjoint_info} if self.eval_adjoint_funcs else {} + ) + + pbar = self._make_adjoint_pbar(idn) + sim_times = np.array(jlres.time) # Unit: seconds + adjoint_index_final = None # captured from the last objective + + for col, info in pbar: + # Resolve the evaluation points for this objective. + adj_seconds, adj_index = self._resolve_adjoint_steps(info) + adjoint_index_final = adj_index + + # Closest simulation step to each requested time (Julia is 1-indexed). + adj_step_idx = [ + int(np.argmin(np.abs(sim_times - s))) + 1 for s in adj_seconds + ] + + funcs = well_QOI_objective( + info.wellID, info.phase, adj_seconds, + adj_step_idx, info.is_rate, julia=julia, + ) if self.adjoint_pbar: - pbar.close() + pbar.set_description_str(f"Adjoints for {col}") + + for i, func in enumerate(funcs): + julia.func = func + grad = self._solve_adjoint(julia) + + # Cross-check Julia's objective evaluation against the value + # we already obtained from the forward extraction. + func_val = julia.Jutul.evaluate_objective(func, case, jlres.result) + expected = pyres.loc[adj_index[i]][col] + assert np.isclose(func_val, expected), ( + f"func_val: {func_val:.3e}, pyres: {expected:.3e}" + ) + if self.eval_adjoint_funcs: + func_dict[col].append(func_val) + + # Post-process the raw Julia gradient for each requested parameter. + for param in info.parameters: + grad_dict[(col, param)].append( + _extract_adjoint(grad, case, param, actnum_vec, + self.perm_copied, julia) + ) - # Create adjoint dataframe with MultiIndex columns - cols = pd.MultiIndex.from_tuples(grad_dict.keys()) - adjoints = pd.DataFrame(grad_dict, columns=cols, index=self.index[1]) - adjoints.index.name = self.index[0] - adjoints.attrs = {'units': units} - # ---------------------------------------------------------------------------------------------- + if self.adjoint_pbar: + pbar.close() - os.chdir('..') + # Wrap the gradients into a multi-indexed DataFrame. + cols = pd.MultiIndex.from_tuples(grad_dict.keys()) + adjoints = pd.DataFrame(grad_dict, columns=cols, index=adjoint_index_final) + adjoints.index.name = self.index[0] + adjoints.attrs = {"units": units} - # Function values dataframe if self.eval_adjoint_funcs: fun = pd.DataFrame(func_dict, index=adjoints.index) fun.index.name = adjoints.index.name self.adjoint_funcs = fun - - # Delete simulation folder - if delete_folder: - shutil.rmtree(folder) - - if self.compute_adjoints: - return output, adjoints - else: - return output + return adjoints - - def extract_datatypes(self, jlres, jlcase, units, julia) -> tuple[pd.DataFrame, np.ndarray]: - res = {} - attrs = {} # For datatype units + def _solve_adjoint(self, julia): + """ + Invoke the appropriate JutulDarcy adjoint solver for the current mode. - if isinstance(units, str): - jl_units = julia.Symbol(units) - else: - jl_units = units # julia.missing or already a Julia Symbol + Parameters + ---------- + julia : juliacall.Main + Julia main module. The names ``case``, ``res``, ``grad_case`` and + ``func`` are expected to be already bound in the Julia namespace. - # Summary - smry = julia.JutulDarcy.summary_result(jlcase, jlres, jl_units) + Returns + ------- + Any + Raw Julia gradient object. + """ + if self.adjoint_mode == "sensitivities": + return _suppress_julia( + julia, + "JutulDarcy.reservoir_sensitivities(" + "case, res, func, include_parameters=true)" + ) + return _suppress_julia( + julia, + "JutulDarcy.parameters_gradient_reservoir(grad_case, func, deps=:case)" + ) - # Start date - start_date = jlcase.input_data["RUNSPEC"]['START'] - - # Get report steps - report_days = [] - sim_days = np.array(list(smry["TIME"].seconds), dtype=np.int64) / (24*60*60) - for d, sday in enumerate(sim_days): - if self.index[0] == 'days': - if sday in self.index[1]: - report_days.append(int(sday)) - elif self.index[0] == 'dates': - sim_date = start_date + dt.timedelta(days=sday) - if sim_date in self.index[1]: - report_days.append(int(sday)) - else: - raise ValueError(f"Invalid report type: {self.index[0]}. Must be 'days' or 'dates'.") - - report_daysIDX = np.argwhere(np.isin(sim_days, report_days)).flatten() - - # Extract datatypes for step - for datatype in self.datatype: - - # Well results - if ':' in datatype: - baseID, wellID = datatype.split(':') - jlres_wells = smry["VALUES"]["WELLS"] - - if wellID in jlres_wells: - data = jlres_wells[wellID] - else: - raise ValueError(f"Well ID '{wellID}' not found in simulation results for datatype '{baseID}'") - - if baseID in data: - res[datatype] = np.array(data[baseID])[report_daysIDX] - attrs[datatype] = get_metric_unit(baseID) - else: - raise ValueError(f"Datatype '{baseID}' not found for well '{wellID}' in simulation results") - - # Field results - else: - jlres_field = smry["VALUES"]["FIELD"] - if datatype in jlres_field: - data = np.array(jlres_field[datatype]) - res[datatype] = data[report_daysIDX] - attrs[datatype] = get_metric_unit(datatype) - else: - raise ValueError(f"Datatype '{datatype}' not found in field results of simulation") - - # TODO: Add support for other datatypes (saturation, pressure, PERMX, etc.) - - # Make DataFrame - res = pd.DataFrame(res, index=self.index[1]) - res.index.name = self.index[0] - res.attrs = attrs - - return res, report_daysIDX - - - def render_makofile(self, makofile: str, folder: str, input: dict): - datafile = makofile.replace('.mako', '.DATA') - template = Template(filename=makofile) - with open(os.path.join(folder, datafile), 'w') as f: - f.write(template.render(**input)) - - return datafile + def _resolve_adjoint_steps(self, info: AdjointObjective): + """ + Convert an objective's step specification into seconds + index labels. + Parameters + ---------- + info : AdjointObjective + Spec containing ``steps`` (``"all"``, list of ints, or list of + :class:`datetime.datetime`). - -def _extract_adjoint(jlgrad, jlcase, parameter, actnum, is_rate, julia): - - if 'poro' in parameter.lower(): - adjoint = np.asarray(jlgrad[julia.Symbol("porosity")]) # F-order array (only active cells) - adjoint = _active_to_full_grid(adjoint, actnum) - - elif 'perm' in parameter.lower(): - adjoint = np.asarray(jlgrad[julia.Symbol("permeability")]) # F-order array (only active cells) - - if 'permx' in parameter.lower(): index = 0 - if 'permy' in parameter.lower(): index = 1 - if 'permz' in parameter.lower(): index = 2 - - adjoint = _active_to_full_grid(adjoint[index], actnum) # SI: per m2, convert to per mD - - if 'log' in parameter.lower(): - perm = np.array(jlcase.input_data['GRID'][['PERMX', 'PERMY', 'PERMZ'][index]]) - adjoint = adjoint * perm.flatten(order='F') - # Note: For dJ/dlog(perm) = dJ/dperm * perm, no need to convert from m2 to mD. - else: - adjoint = adjoint * julia.seval('si_unit(:milli)*si_unit(:darcy)') - else: - raise ValueError(f"Adjoint not implemented for parameter '{parameter}'") - - if is_rate: - adjoint = adjoint * julia.seval('si_unit(:day)') # Convert from per sec to per day - - return adjoint - - -def _active_to_full_grid(vec, actnum_vec, fill_value=0.0): - if len(vec) == actnum_vec.sum(): - full_vec = np.full(actnum_vec.shape, fill_value, dtype=np.float64, order='F') - full_vec[actnum_vec == 1] = vec - return full_vec - if len(vec) == len(actnum_vec): - return vec.astype(np.float64, copy=False) + Returns + ------- + tuple + ``(seconds_array, index_labels)`` ready for adjoint evaluation. - raise ValueError("Parameter length does not match number of active cells") - - -def _process_datatype_info(datatypes): - processed = [] - for dataID in datatypes: - if ':' in dataID: - s = dataID.split(':') - baseID = s[0] - wellIDs = s[1:] - for wID in wellIDs: - processed.append(f'{baseID}:{wID}') - else: - processed.append(dataID) - return processed - -def _process_adjoint_info(adjoint_info): - phase_map = { - 'WOPT': ('oil', False), - 'WGPT': ('gas', False), - 'WWPT': ('water', False), - 'WLPT': ('liquid', False), - 'WOPR': ('oil', True), - 'WGPR': ('gas', True), - 'WWPR': ('water', True), - 'WWIR': ('water', True), - 'WLPR': ('liquid', True), - } - info = {} - for dataID in adjoint_info: - wellID = adjoint_info[dataID]['wellID'] - if isinstance(wellID, str): wellID = [wellID] - - parameters = adjoint_info[dataID]['parameters'] - if isinstance(parameters, str): - parameters = [parameters] - - steps = adjoint_info[dataID]['steps'] # 'all' or 'acc' or list of steps (days or dates) - phaseID, is_rate = phase_map[dataID] - for wID in wellID: - info[f'{dataID}:{wID}'] = { - 'wellID': wID, - 'phase': phaseID, - 'is_rate': is_rate, - 'parameters': parameters, - 'steps': steps, - } - return info - - -def get_metric_unit(key: str) -> str: - unit_map = { - "PORO": "", - "PERMX": "mD", - "PERMY": "mD", - "PERMZ": "mD", - "FOPT": "Sm3", - "FGPT": "Sm3", - "FWPT": "Sm3", - "FWLT": "Sm3", - "FWIT": "Sm3", - "FOPR": "Sm3/day", - "FGPR": "Sm3/day", - "FWPR": "Sm3/day", - "FLPR": "Sm3/day", - "FWIR": "Sm3/day", - "WOPR": "Sm3/day", - "WGPR": "Sm3/day", - "WWPR": "Sm3/day", - "WLPR": "Sm3/day", - "WWIR": "Sm3/day", - } - return unit_map.get(key.upper(), "Unknown") - - -def well_QOI_objective(wellID, phaseID, stepID=None, time=None, rate=True, julia=None): - if julia is None: - from juliacall import Main as julia - julia.seval("using JutulDarcy") + Raises + ------ + TypeError + If the step list contains an unsupported element type. + """ + if info.steps == "all": + return self.report_seconds, self.index[1] + + first = info.steps[0] + if isinstance(first, int): + return ( + np.array(info.steps, dtype=np.int64) * SECONDS_PER_DAY, + info.steps, + ) + if isinstance(first, dt.datetime): + return ( + np.array( + [(d - self.start_date).total_seconds() for d in info.steps], + dtype=np.int64, + ), + info.steps, + ) + raise TypeError(f"Unsupported adjoint step type: {type(first).__name__}") - dt_factor = "" if rate else "dt*" - - rateID_map = { - "mass" : "TotalSurfaceMassRate", - "liquid" : "SurfaceLiquidRateTarget", - "water" : "SurfaceWaterRateTarget", - "oil" : "SurfaceOilRateTarget", - "gas" : "SurfaceGasRateTarget", - "rate" : "TotalRateTarget", - } - if phaseID not in rateID_map: - raise ValueError(f"Unknown rate type: {phaseID}") - rateID_symbol = rateID_map[phaseID] - - if stepID is None: - julia.seval( - f""" - function well_QOI(model, state, dt, step_i, forces) - ctrl = forces[:Facility].control[Symbol("{wellID}")] - if ctrl isa JutulDarcy.DisabledControl - return 0.0 - elseif ctrl isa JutulDarcy.ProducerControl - sgn = -1.0 - else - sgn = 1.0 - end - rate = JutulDarcy.compute_well_qoi( - model, - state, - forces, - Symbol("{wellID}"), - {rateID_symbol} - ) - return sgn*{dt_factor}rate - end - """ - ) - return julia.well_QOI + def _make_adjoint_pbar(self, idn: int): + """ + Build an iterator over ``self.adjoint_info``, optionally with a pbar. - if isinstance(stepID, (list, np.ndarray)): - qois = [] - for sID in stepID: - julia.seval( - f""" - function well_QOI_{sID}(model, state, dt, step_i, forces) - if (step_i[:step] != {sID+1}) && (step_i[:time] != {time[sID+1]}) - return 0.0 - else - ctrl = forces[:Facility].control[Symbol("{wellID}")] - if ctrl isa JutulDarcy.DisabledControl - return 0.0 - elseif ctrl isa JutulDarcy.ProducerControl - sgn = -1.0 - else - sgn = 1.0 - end - rate = JutulDarcy.compute_well_qoi( - model, - state, - forces, - Symbol("{wellID}"), - {rateID_symbol} - ) - return sgn*{dt_factor}rate - end - end - """ - ) - qois.append(julia.seval(f"well_QOI_{sID}")) - return qois - - julia.seval( - f""" - function well_QOI(model, state, dt, step_i, forces) - if (step_i[:step] != {stepID+1}) && (step_i[:time] != {time[stepID+1]}) - return 0.0 - else - ctrl = forces[:Facility].control[Symbol("{wellID}")] - if ctrl isa JutulDarcy.DisabledControl - return 0.0 - elseif ctrl isa JutulDarcy.ProducerControl - sgn = -1.0 - else - sgn = 1.0 - end - rate = JutulDarcy.compute_well_qoi( - model, - state, - forces, - Symbol("{wellID}"), - {rateID_symbol} - ) - return sgn*{dt_factor}rate - end - end + Parameters + ---------- + idn : int + Ensemble index (used to stagger pbar positions in parallel runs). + + Returns + ------- + Iterable + Either ``self.adjoint_info.items()`` or a wrapped :class:`tqdm`. """ - ) - return julia.well_QOI \ No newline at end of file + if not self.adjoint_pbar: + return self.adjoint_info.items() + + # Drop the default colour so we can override it for the adjoint bars. + opts = {k: v for k, v in PBAR_OPTS.items() if k != "colour"} + return tqdm( + self.adjoint_info.items(), + desc="Solving adjoints", + position=idn + 1, + leave=False, + unit="obj", + dynamic_ncols=False, + colour="#713996", + **opts, + ) \ No newline at end of file diff --git a/tests/test_jutul_darcy.py b/tests/test_jutul_darcy.py new file mode 100644 index 0000000..1105c4e --- /dev/null +++ b/tests/test_jutul_darcy.py @@ -0,0 +1,158 @@ +"""Integration tests for the JutulDarcy wrapper. + +These tests execute a real simulation and are intentionally heavier than unit tests. +Run selectively when validating simulator integration. +""" + +from contextlib import contextmanager +from datetime import datetime +from pathlib import Path + +import os +import shutil +import numpy as np +import pytest + +from subsurface.multphaseflow.jutul_darcy import JutulDarcy + + +REPORT_DATES = [ + datetime(2023, 2, 5), + datetime(2024, 3, 11), + datetime(2025, 4, 15), + datetime(2026, 5, 20), + datetime(2027, 6, 24), + datetime(2028, 7, 28), + datetime(2029, 9, 1), + datetime(2030, 10, 6), + datetime(2031, 11, 10), + datetime(2032, 12, 14), +] + +DATA_TYPES = [ + "WOPR:PRO1", + "WOPR:PRO2", + "WOPR:PRO3", + "WWPR:PRO1", + "WWPR:PRO2", + "WWPR:PRO3", + "WWIR:INJ1", +] + +GRADIENT_STEP = datetime(2032, 12, 14) +GRADIENT_TARGET = "WOPR:PRO2" + + +def _tiny_folder() -> Path: + """Return absolute path to the `Example/TINY` input case directory.""" + tiny_path = Path(__file__).resolve().parents[1] / "Example" / "TINY" + if not tiny_path.exists(): + raise FileNotFoundError(f"TINY folder not found at: {tiny_path}") + return tiny_path + + +@contextmanager +def _working_directory(path: Path): + """Temporarily change current working directory.""" + previous = Path.cwd() + os.chdir(path) + try: + yield + finally: + os.chdir(previous) + + +def _copy_case_folder(target: Path, source: Path, name: str) -> Path: + """Copy the TINY case folder into a per-test destination.""" + case_path = target / name + shutil.copytree(source, case_path) + return case_path + + +def _run_case(case_path: Path, options: dict): + """Run JutulDarcy for the given case path and options.""" + log_permx = np.log(np.load(case_path / "PERMX.npy")) + simulator = JutulDarcy(options) + with _working_directory(case_path): + return simulator({"log_permx": log_permx}) + + +@pytest.fixture(scope="module") +def options() -> dict: + """Base simulation options shared by tests.""" + return { + "reporttype": "dates", + "reportpoint": REPORT_DATES, + "runfile": "RUNFILE.mako", + "datatype": DATA_TYPES, + } + + +@pytest.fixture(scope="module") +def tiny_folder() -> Path: + """Path to the immutable source TINY case in repository.""" + return _tiny_folder() + + +@pytest.fixture(scope="module") +def run_simulation_with_adjoint(tmp_path_factory, options, tiny_folder): + """Run one adjoint-enabled simulation and share output across gradient tests.""" + base_tmp = tmp_path_factory.mktemp("adjoint") + case_path = _copy_case_folder(base_tmp, tiny_folder, "TINY_ADJOINT") + + case_options = { + **options, + 'perm_copied': True, # Include total derivative when PERMX is copied to PERMY and PERMZ + "adjoints": { + "WOPR": { + "steps": [GRADIENT_STEP], + "wellID": "PRO2", + "parameters": ["log_permx", "permx"], + } + }, + } + return _run_case(case_path, case_options) + + +def test_simulation_runs_and_matches_requested_outputs(tmp_path, options, tiny_folder): + """Simulation returns a non-empty table with all requested columns and dates.""" + case_path = _copy_case_folder(tmp_path, tiny_folder, "TINY") + results = _run_case(case_path, options) + + assert not results.empty, "Simulation result table is empty" + + missing_columns = sorted(set(options["datatype"]) - set(results.columns)) + assert not missing_columns, f"Missing expected result columns: {missing_columns}" + + missing_dates = [date for date in options["reportpoint"] if date not in results.index] + assert not missing_dates, f"Missing expected report dates: {missing_dates}" + + +def test_gradient_contains_expected_structure(run_simulation_with_adjoint, options): + """Adjoint run returns expected gradient columns and valid report dates.""" + _, gradient = run_simulation_with_adjoint + + assert not gradient.empty, "Gradient table is empty" + + expected_columns = { + (GRADIENT_TARGET, "log_permx"), + (GRADIENT_TARGET, "permx"), + } + missing_columns = sorted(expected_columns - set(gradient.columns)) + assert not missing_columns, f"Missing expected gradient columns: {missing_columns}" + + unexpected_dates = [date for date in gradient.index if date not in options["reportpoint"]] + assert not unexpected_dates, f"Gradient has unexpected report dates: {unexpected_dates}" + + +def test_gradient_log_permx_matches_chain_rule(run_simulation_with_adjoint, tiny_folder): + """Validate dF/d(log(k)) = dF/d(k) * k for the configured objective and step.""" + _, gradient = run_simulation_with_adjoint + + grad_log_permx = gradient.loc[GRADIENT_STEP, (GRADIENT_TARGET, "log_permx")] + grad_permx = gradient.loc[GRADIENT_STEP, (GRADIENT_TARGET, "permx")] + permx = np.load(tiny_folder / "PERMX.npy") + + np.testing.assert_allclose(grad_log_permx, grad_permx * permx) + +