From 2a742b3aa91c28b9e920ebf8edb834cea0251e4e Mon Sep 17 00:00:00 2001 From: bretth Date: Fri, 21 Nov 2025 12:28:26 +0100 Subject: [PATCH 1/3] Adding an example to the test_doc_examples --- .../deterministic_doublet.py | 2 +- src/pythermogis/thermogis_classes/doublet.py | 3 + tests/test_doc_examples.py | 445 +++++++++++++----- 3 files changed, 337 insertions(+), 113 deletions(-) diff --git a/src/pythermogis/doublet_simulation/deterministic_doublet.py b/src/pythermogis/doublet_simulation/deterministic_doublet.py index 02c75b5..4e560b0 100644 --- a/src/pythermogis/doublet_simulation/deterministic_doublet.py +++ b/src/pythermogis/doublet_simulation/deterministic_doublet.py @@ -132,7 +132,7 @@ def validate_input(reservoir_properties: xr.Dataset): # check all necessary variables are present missing_variables = [] for variable in ["thickness", "porosity", "ntg", "depth"]: - if variable not in reservoir_properties: + if variable not in reservoir_properties and variable not in reservoir_properties.dims: missing_variables.append(variable) if len(missing_variables) > 0: diff --git a/src/pythermogis/thermogis_classes/doublet.py b/src/pythermogis/thermogis_classes/doublet.py index 3360a82..2bea3e9 100644 --- a/src/pythermogis/thermogis_classes/doublet.py +++ b/src/pythermogis/thermogis_classes/doublet.py @@ -122,6 +122,9 @@ def calculate_performance_of_single_location(mask: float, depth: float, thicknes results = doublet.calculateDoubletPerformance(input) # If calculation was not successful, return mask value + if results is None: + return (np.nan,) * 14 + if results.utc() == -9999.0: return (mask_value,) * 14 diff --git a/tests/test_doc_examples.py b/tests/test_doc_examples.py index faa28fd..b3761bf 100644 --- a/tests/test_doc_examples.py +++ b/tests/test_doc_examples.py @@ -2,102 +2,130 @@ def test_example_1(): from pythermogis import calculate_doublet_performance import xarray as xr - input_data = xr.Dataset({ - "thickness": 300, - "ntg": 0.5, - "porosity": 0.5, - "depth": 2000, - "permeability": 300, - }) + input_data = xr.Dataset( + { + "thickness": 300, + "ntg": 0.5, + "porosity": 0.5, + "depth": 2000, + "permeability": 300, + } + ) results = calculate_doublet_performance(input_data) print(results) + def test_example_2(): from pythermogis import calculate_doublet_performance import xarray as xr import numpy as np - input_data = xr.Dataset({ - "thickness": (("x", "y"), np.array([[300, 300], [200, 200]])), - "ntg": (("x", "y"), np.array([[0.5, 0.5], [0.25, 0.25]])), - "porosity": (("x", "y"), np.array([[0.5, 0.5], [0.75, 0.7]])), - "depth": (("x", "y"), np.array([[5000, 5000], [4500, 4500]])), - "permeability": (("x", "y"), np.array([[300, 300], [450, 450]])), - }, coords={"x": [0, 1], "y": [10, 20]}) + input_data = xr.Dataset( + { + "thickness": (("x", "y"), np.array([[300, 300], [200, 200]])), + "ntg": (("x", "y"), np.array([[0.5, 0.5], [0.25, 0.25]])), + "porosity": (("x", "y"), np.array([[0.5, 0.5], [0.75, 0.7]])), + "depth": (("x", "y"), np.array([[5000, 5000], [4500, 4500]])), + "permeability": (("x", "y"), np.array([[300, 300], [450, 450]])), + }, + coords={"x": [0, 1], "y": [10, 20]}, + ) results = calculate_doublet_performance(input_data) print(results) + def test_example_3(): - from pythermogis import calculate_doublet_performance, instantiate_utc_properties_builder + from pythermogis import ( + calculate_doublet_performance, + instantiate_utc_properties_builder, + ) import xarray as xr - input_data = xr.Dataset({ - "thickness": 300, - "ntg": 0.5, - "porosity": 0.5, - "depth": 2000, - "permeability": 300, - }) + input_data = xr.Dataset( + { + "thickness": 300, + "ntg": 0.5, + "porosity": 0.5, + "depth": 2000, + "permeability": 300, + } + ) utc_properties = instantiate_utc_properties_builder().setUseHeatPump(True).build() results = calculate_doublet_performance(input_data, utc_properties=utc_properties) print(results) + def test_example_4(): from pythermogis import calculate_doublet_performance_stochastic import xarray as xr - input_data = xr.Dataset({ - "thickness_mean": 300, - "thickness_sd": 50, - "ntg": 0.5, - "porosity": 0.5, - "depth": 2000, - "ln_permeability_mean": 4, - "ln_permeability_sd": 0.5, - }) + input_data = xr.Dataset( + { + "thickness_mean": 300, + "thickness_sd": 50, + "ntg": 0.5, + "porosity": 0.5, + "depth": 2000, + "ln_permeability_mean": 4, + "ln_permeability_sd": 0.5, + } + ) results = calculate_doublet_performance_stochastic(input_data) print(results) + def test_example_5(): from pythermogis import calculate_doublet_performance_stochastic import xarray as xr from matplotlib import pyplot as plt from pathlib import Path from os import path - output_data_path = Path(path.dirname(__file__), "resources") / "test_output" / "example_data" + + output_data_path = ( + Path(path.dirname(__file__), "resources") / "test_output" / "example_data" + ) output_data_path.mkdir(parents=True, exist_ok=True) - input_data = xr.Dataset({ - "thickness_mean": 300, - "thickness_sd": 50, - "ntg": 0.5, - "porosity": 0.5, - "depth": 2000, - "ln_permeability_mean": 5, - "ln_permeability_sd": 0.5, - }) - results = calculate_doublet_performance_stochastic(input_data, p_values=[1, 10, 20, 30, 40, 50, 60, 70, 80, 90, 99]) + input_data = xr.Dataset( + { + "thickness_mean": 300, + "thickness_sd": 50, + "ntg": 0.5, + "porosity": 0.5, + "depth": 2000, + "ln_permeability_mean": 5, + "ln_permeability_sd": 0.5, + } + ) + results = calculate_doublet_performance_stochastic( + input_data, p_values=[1, 10, 20, 30, 40, 50, 60, 70, 80, 90, 99] + ) fig, axe = plt.subplots(figsize=(10, 5)) kh = results.transmissivity_with_ntg * 1.0 kh.plot(y="p_value", ax=axe) axe.set_title(f"Aquifer Net transmissivity\n kH") - temp = input_data.temperature + temp = input_data.temperature inj_temp = results.inj_temp.sel(p_value=50, method="nearest") prd_temp = results.prd_temp.sel(p_value=50, method="nearest") - axe.axhline(50.0, label=f"TEMP res: {temp:.1f} inj: {inj_temp:.1f} prd: {prd_temp:.1f} ", - ls="--", c="tab:orange") + axe.axhline( + 50.0, + label=f"TEMP res: {temp:.1f} inj: {inj_temp:.1f} prd: {prd_temp:.1f} ", + ls="--", + c="tab:orange", + ) # axes.axhline(pos, label=f"probability of success: {pos:.1f}%",ls="--", c="tab:orange") kh50 = kh.sel(p_value=50, method="nearest") - axe.axvline(kh50, ls="--", c="tab:orange" ) + axe.axvline(kh50, ls="--", c="tab:orange") axe.legend() axe.set_xlabel("kH [Dm]") axe.set_ylabel("p-value [%]") plt.savefig(output_data_path / "kh.png") + def test_example_6(): from pythermogis import calculate_doublet_performance_stochastic, calculate_pos from pygridsio import read_grid, plot_grid, resample_grid @@ -108,10 +136,14 @@ def test_example_6(): from os import path # the location of the input data: the data can be found in the resources/example_data directory of the repo - input_data_path = Path(path.dirname(__file__), "resources") / "test_input" / "example_data" + input_data_path = ( + Path(path.dirname(__file__), "resources") / "test_input" / "example_data" + ) # create a directory to write the output files to - output_data_path = Path(path.dirname(__file__), "resources") / "test_output" / "example_data" + output_data_path = ( + Path(path.dirname(__file__), "resources") / "test_output" / "example_data" + ) output_data_path.mkdir(parents=True, exist_ok=True) # if set to True then simulation is always run, otherwise pre-calculated results are read (if available) @@ -119,51 +151,96 @@ def test_example_6(): # grids can be in .nc, .asc, .zmap or .tif format: see pygridsio package new_cellsize = 20000 # in m, this sets the resolution of the model; to speed up calcualtions you can increase the cellsize - reservoir_properties = resample_grid(read_grid(input_data_path / "ROSL_ROSLU__thick.zmap"), new_cellsize=new_cellsize).to_dataset(name="thickness_mean") - reservoir_properties["thickness_sd"] = resample_grid(read_grid(input_data_path / "ROSL_ROSLU__thick_sd.zmap"), new_cellsize=new_cellsize) - reservoir_properties["ntg"] = resample_grid(read_grid(input_data_path / "ROSL_ROSLU__ntg.zmap"), new_cellsize=new_cellsize) - reservoir_properties["porosity"] = resample_grid(read_grid(input_data_path / "ROSL_ROSLU__poro.zmap"), new_cellsize=new_cellsize) / 100 - reservoir_properties["depth"] = resample_grid(read_grid(input_data_path / "ROSL_ROSLU__top.zmap"), new_cellsize=new_cellsize) - reservoir_properties["ln_permeability_mean"] = np.log(resample_grid(read_grid(input_data_path / "ROSL_ROSLU__perm.zmap"), new_cellsize=new_cellsize)) - reservoir_properties["ln_permeability_sd"] = resample_grid(read_grid(input_data_path / "ROSL_ROSLU__ln_perm_sd.zmap"), new_cellsize=new_cellsize) + reservoir_properties = resample_grid( + read_grid(input_data_path / "ROSL_ROSLU__thick.zmap"), new_cellsize=new_cellsize + ).to_dataset(name="thickness_mean") + reservoir_properties["thickness_sd"] = resample_grid( + read_grid(input_data_path / "ROSL_ROSLU__thick_sd.zmap"), + new_cellsize=new_cellsize, + ) + reservoir_properties["ntg"] = resample_grid( + read_grid(input_data_path / "ROSL_ROSLU__ntg.zmap"), new_cellsize=new_cellsize + ) + reservoir_properties["porosity"] = ( + resample_grid( + read_grid(input_data_path / "ROSL_ROSLU__poro.zmap"), + new_cellsize=new_cellsize, + ) + / 100 + ) + reservoir_properties["depth"] = resample_grid( + read_grid(input_data_path / "ROSL_ROSLU__top.zmap"), new_cellsize=new_cellsize + ) + reservoir_properties["ln_permeability_mean"] = np.log( + resample_grid( + read_grid(input_data_path / "ROSL_ROSLU__perm.zmap"), + new_cellsize=new_cellsize, + ) + ) + reservoir_properties["ln_permeability_sd"] = resample_grid( + read_grid(input_data_path / "ROSL_ROSLU__ln_perm_sd.zmap"), + new_cellsize=new_cellsize, + ) # simulate # if doublet simulation has already been run then read in results, or run the simulation and write results out if (output_data_path / "output_results.nc").exists and not run_simulation: results = xr.load_dataset(output_data_path / "output_results.nc") else: - results = calculate_doublet_performance_stochastic(reservoir_properties, p_values=[10, 20, 30, 40, 50, 60, 70, 80, 90]) - results.to_netcdf(output_data_path / "output_results.nc") # write doublet simulation results to file + results = calculate_doublet_performance_stochastic( + reservoir_properties, p_values=[10, 20, 30, 40, 50, 60, 70, 80, 90] + ) + results.to_netcdf( + output_data_path / "output_results.nc" + ) # write doublet simulation results to file # plot results as maps - variables_to_plot = ["transmissivity","power", "pres", "npv"] + variables_to_plot = ["transmissivity", "power", "pres", "npv"] p_values_to_plot = [10, 50, 90] - fig, axes = plt.subplots(nrows=len(variables_to_plot), ncols=len(p_values_to_plot), figsize=(5*len(variables_to_plot), 5*len(p_values_to_plot)), sharex=True, sharey=True) + fig, axes = plt.subplots( + nrows=len(variables_to_plot), + ncols=len(p_values_to_plot), + figsize=(5 * len(variables_to_plot), 5 * len(p_values_to_plot)), + sharex=True, + sharey=True, + ) for j, p_value in enumerate(p_values_to_plot): results_p_value = results.sel(p_value=p_value) for i, variable in enumerate(variables_to_plot): - plot_grid(results_p_value[variable], axes=axes[i,j], add_netherlands_shapefile=True) # See documentation on plot_grid in pygridsio, you can also provide your own shapefile + plot_grid( + results_p_value[variable], + axes=axes[i, j], + add_netherlands_shapefile=True, + ) # See documentation on plot_grid in pygridsio, you can also provide your own shapefile plt.tight_layout() # ensure there is enough spacing plt.savefig(output_data_path / "maps.png") # plot net-present value at a single location as a function of p-value and find the probability of success - results["pos"] = calculate_pos(results) # calculate probability of success across whole aquifer + results["pos"] = calculate_pos( + results + ) # calculate probability of success across whole aquifer x, y = 125e3, 525e3 # define location - results_loc = results.sel(x=x,y=y,method="nearest") # select only the location of interest - pos = results_loc.pos.data # get probability of success at location of interest as a single value + results_loc = results.sel( + x=x, y=y, method="nearest" + ) # select only the location of interest + pos = ( + results_loc.pos.data + ) # get probability of success at location of interest as a single value # plot npv versus p-value and a map showing the location of interest - fig, axes = plt.subplots(ncols=2, figsize=(10,5)) + fig, axes = plt.subplots(ncols=2, figsize=(10, 5)) results_loc.npv.plot(y="p_value", ax=axes[0]) axes[0].set_title(f"Aquifer: ROSL_ROSLU\nlocation: [{x:.0f}m, {y:.0f}m]") - axes[0].axhline(pos, label=f"probability of success: {pos:.1f}%",ls="--", c="tab:orange") - axes[0].axvline(0.0,ls="--", c="tab:orange") + axes[0].axhline( + pos, label=f"probability of success: {pos:.1f}%", ls="--", c="tab:orange" + ) + axes[0].axvline(0.0, ls="--", c="tab:orange") axes[0].legend() axes[0].set_xlabel("net-present-value [Million €]") axes[0].set_ylabel("p-value [%]") - plot_grid(results.pos,axes=axes[1], add_netherlands_shapefile=True) - axes[1].scatter(x,y,marker="x",color="tab:red") + plot_grid(results.pos, axes=axes[1], add_netherlands_shapefile=True) + axes[1].scatter(x, y, marker="x", color="tab:red") plt.tight_layout() # ensure there is enough spacing plt.savefig(output_data_path / "pos.png") @@ -178,21 +255,41 @@ def test_example_7(): from pathlib import Path from os import path - input_data_path = Path(path.dirname(__file__), "resources") / "test_input" / "example_data" - output_data_path = Path(path.dirname(__file__), "resources") / "test_output" / "example_data" + input_data_path = ( + Path(path.dirname(__file__), "resources") / "test_input" / "example_data" + ) + output_data_path = ( + Path(path.dirname(__file__), "resources") / "test_output" / "example_data" + ) output_data_path.mkdir(parents=True, exist_ok=True) # read in reservoir property maps - reservoir_properties = read_grid(input_data_path / "ROSL_ROSLU__thick.zmap").to_dataset(name="thickness_mean") - reservoir_properties["thickness_sd"] = read_grid(input_data_path / "ROSL_ROSLU__thick_sd.zmap") + reservoir_properties = read_grid( + input_data_path / "ROSL_ROSLU__thick.zmap" + ).to_dataset(name="thickness_mean") + reservoir_properties["thickness_sd"] = read_grid( + input_data_path / "ROSL_ROSLU__thick_sd.zmap" + ) reservoir_properties["ntg"] = read_grid(input_data_path / "ROSL_ROSLU__ntg.zmap") - reservoir_properties["porosity"] = 0.01* read_grid(input_data_path / "ROSL_ROSLU__poro.zmap") + reservoir_properties["porosity"] = 0.01 * read_grid( + input_data_path / "ROSL_ROSLU__poro.zmap" + ) reservoir_properties["depth"] = read_grid(input_data_path / "ROSL_ROSLU__top.zmap") - reservoir_properties["ln_permeability_mean"] = np.log( read_grid(input_data_path / "ROSL_ROSLU__perm.zmap")) - reservoir_properties["ln_permeability_sd"] = read_grid(input_data_path / "ROSL_ROSLU__ln_perm_sd.zmap") + reservoir_properties["ln_permeability_mean"] = np.log( + read_grid(input_data_path / "ROSL_ROSLU__perm.zmap") + ) + reservoir_properties["ln_permeability_sd"] = read_grid( + input_data_path / "ROSL_ROSLU__ln_perm_sd.zmap" + ) # Select locations of interest from the reservoir properties and simulate doublet performance - portfolio_coords = [(125e3, 525e3), (100e3, 525e3), (110e3, 525e3), (125e3, 515e3), (125e3, 520e3)] + portfolio_coords = [ + (125e3, 525e3), + (100e3, 525e3), + (110e3, 525e3), + (125e3, 515e3), + (125e3, 520e3), + ] # Collect selected locations, create a new Dataset with the coordinate 'location' locations = [] @@ -202,18 +299,19 @@ def test_example_7(): locations.append(point) portfolio_reservoir_properties = xr.concat(locations, dim="location") - results_portfolio = calculate_doublet_performance_stochastic(portfolio_reservoir_properties, p_values=[10, 20, 30, 40, 50, 60, 70, 80, 90]) + results_portfolio = calculate_doublet_performance_stochastic( + portfolio_reservoir_properties, p_values=[10, 20, 30, 40, 50, 60, 70, 80, 90] + ) # post process: clip the npv and calculate the probability of success AEC = -3.5 - results_portfolio['npv'] = results_portfolio.npv.clip(min=AEC) + results_portfolio["npv"] = results_portfolio.npv.clip(min=AEC) results_portfolio["pos"] = calculate_pos(results_portfolio) # plot results fig, axes = plt.subplots(ncols=2, nrows=2, figsize=(10, 10)) colors = plt.cm.tab10.colors for i_loc, loc in enumerate(portfolio_coords): - # select single portfolio location to plot results_loc = results_portfolio.isel(location=i_loc) pos = results_loc.pos.data @@ -223,7 +321,12 @@ def test_example_7(): ax = axes[0, 0] results_loc.npv.plot(y="p_value", ax=ax, color=colors[i_loc]) ax.set_title(f"Aquifer: ROSL_ROSLU\n POS ") - ax.axhline(pos, label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", ls="--", c=colors[i_loc]) + ax.axhline( + pos, + label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", + ls="--", + c=colors[i_loc], + ) ax.axvline(0.0, ls="--", c=colors[i_loc]) ax.legend() ax.set_xlabel("net-present-value [Million €]") @@ -236,8 +339,12 @@ def test_example_7(): temp = results_loc.temperature.sel(p_value=50) inj_temp = results_loc.inj_temp.sel(p_value=50) prd_temp = results_loc.prd_temp.sel(p_value=50) - ax.axhline(50.0, label=f"POS: {pos:.1f}% TEMP res: {temp:.1f} inj: {inj_temp:.1f} prd: {prd_temp:.1f} ", - ls="--", c=colors[i_loc]) + ax.axhline( + 50.0, + label=f"POS: {pos:.1f}% TEMP res: {temp:.1f} inj: {inj_temp:.1f} prd: {prd_temp:.1f} ", + ls="--", + c=colors[i_loc], + ) ax.axvline(kh.sel(p_value=50), ls="--", c=colors[i_loc]) ax.legend() ax.set_xlabel("kH [Dm]") @@ -246,7 +353,12 @@ def test_example_7(): ax = axes[1, 0] results_loc.power.plot(y="p_value", ax=ax, color=colors[i_loc]) ax.set_title(f"Aquifer: ROSL_ROSLU\n Power") - ax.axhline(pos, label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", ls="--", c=colors[i_loc]) + ax.axhline( + pos, + label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", + ls="--", + c=colors[i_loc], + ) valp50 = results_loc.power.sel(p_value=50) ax.axvline(valp50, ls="--", c=colors[i_loc]) ax.legend() @@ -256,7 +368,12 @@ def test_example_7(): ax = axes[1, 1] results_loc.cop.plot(y="p_value", ax=ax, color=colors[i_loc]) ax.set_title(f"Aquifer: ROSL_ROSLU\n COP") - ax.axhline(pos, label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", ls="--", c=colors[i_loc]) + ax.axhline( + pos, + label=f"POS: {pos:.1f}% [ {x:.0f}m, {y:.0f}m]", + ls="--", + c=colors[i_loc], + ) valp50 = results_loc.cop.sel(p_value=50) ax.axvline(valp50, ls="--", c=colors[i_loc]) ax.legend() @@ -266,17 +383,23 @@ def test_example_7(): plt.tight_layout() # ensure there is enough spacing plt.savefig(output_data_path / "npv.png") + def test_example8(): - from pythermogis import calculate_doublet_performance, calculate_pos, plot_exceedance + from pythermogis import ( + calculate_doublet_performance, + calculate_pos, + plot_exceedance, + ) import xarray as xr from matplotlib import pyplot as plt from pathlib import Path from os import path import numpy as np - # output directory to write output to - output_data_path = Path(path.dirname(__file__), "resources") / "test_output" / "exceedance_example" + output_data_path = ( + Path(path.dirname(__file__), "resources") / "test_output" / "exceedance_example" + ) output_data_path.mkdir(parents=True, exist_ok=True) # generate simulation samples across desired reservoir properties @@ -284,7 +407,7 @@ def test_example8(): thickness_samples = np.random.uniform(low=150, high=300, size=Nsamples) porosity_samples = np.random.uniform(low=0.5, high=0.8, size=Nsamples) ntg_samples = np.random.uniform(low=0.25, high=0.5, size=Nsamples) - depth_samples = np.random.uniform(low=1800, high=2200, size=Nsamples) + depth_samples = np.random.uniform(low=1800, high=2200, size=Nsamples) permeability_samples = np.random.uniform(low=200, high=800, size=Nsamples) reservoir_properties = xr.Dataset( { @@ -294,48 +417,146 @@ def test_example8(): "depth": (["sample"], depth_samples), "permeability": (["sample"], permeability_samples), }, - coords={"sample": np.arange(Nsamples)} + coords={"sample": np.arange(Nsamples)}, ) # For every sample, run a doublet simulation store the output values - simulations = calculate_doublet_performance(reservoir_properties, print_execution_duration=True) + simulations = calculate_doublet_performance( + reservoir_properties, print_execution_duration=True + ) # post process the samples to calculate the percentiles of each variable - percentiles = np.arange(1,99) - results = simulations.quantile(q=percentiles/100, dim="sample") - results = results.assign_coords(quantile=("quantile", percentiles)) # Overwrite the 'quantile' coordinate with integer percentiles - results = results.rename({"quantile": "p_value"}) # Rename dimension to 'p_value' to match nomenclature in the rest of pythermogis + percentiles = np.arange(1, 99) + results = simulations.quantile(q=percentiles / 100, dim="sample") + results = results.assign_coords( + quantile=("quantile", percentiles) + ) # Overwrite the 'quantile' coordinate with integer percentiles + results = results.rename( + {"quantile": "p_value"} + ) # Rename dimension to 'p_value' to match nomenclature in the rest of pythermogis results["pos"] = calculate_pos(results) # plotting - fig, axes = plt.subplots(3,2, figsize=(10,15)) + fig, axes = plt.subplots(3, 2, figsize=(10, 15)) plt.suptitle("Reservoir Property Distributions") - reservoir_properties.thickness.plot.hist(ax=axes[0,0]) - reservoir_properties.porosity.plot.hist(ax=axes[0,1]) - reservoir_properties.ntg.plot.hist(ax=axes[1,0]) - reservoir_properties.depth.plot.hist(ax=axes[1,1]) - reservoir_properties.permeability.plot.hist(ax=axes[2,0]) - axes[2,1].set_visible(False) + reservoir_properties.thickness.plot.hist(ax=axes[0, 0]) + reservoir_properties.porosity.plot.hist(ax=axes[0, 1]) + reservoir_properties.ntg.plot.hist(ax=axes[1, 0]) + reservoir_properties.depth.plot.hist(ax=axes[1, 1]) + reservoir_properties.permeability.plot.hist(ax=axes[2, 0]) + axes[2, 1].set_visible(False) plt.savefig(output_data_path / "reservoir_distributions.png") - fig, axes = plt.subplots(2,2, figsize=(10,10)) + fig, axes = plt.subplots(2, 2, figsize=(10, 10)) plt.suptitle(f"Simulation results\nNo. of Samples: {Nsamples}") - simulations.plot.scatter(x="permeability", y="transmissivity", ax=axes[0,0]) - simulations.plot.scatter(x="permeability", y="utc", ax=axes[0,1]) - simulations.plot.scatter(x="permeability", y="power", ax=axes[1,0]) - simulations.plot.scatter(x="permeability", y="npv", ax=axes[1,1]) + simulations.plot.scatter(x="permeability", y="transmissivity", ax=axes[0, 0]) + simulations.plot.scatter(x="permeability", y="utc", ax=axes[0, 1]) + simulations.plot.scatter(x="permeability", y="power", ax=axes[1, 0]) + simulations.plot.scatter(x="permeability", y="npv", ax=axes[1, 1]) plt.savefig(output_data_path / "simulation_results.png") - fig, axes = plt.subplots(2,2, figsize=(10,10)) + fig, axes = plt.subplots(2, 2, figsize=(10, 10)) plt.suptitle(f"Exceedance probability curves\nNo. of Samples: {Nsamples}") - plot_exceedance(results["permeability"], ax=axes[0,0]) - plot_exceedance(results["utc"], ax=axes[0,1]) - plot_exceedance(results["power"], ax=axes[1,0]) - plot_exceedance(results["npv"], ax=axes[1,1]) - axes[1, 1].axhline(results["pos"], ls="--", c="tab:orange", label=f"POS: {results["pos"]:.1f}%") # add the probability of success + plot_exceedance(results["permeability"], ax=axes[0, 0]) + plot_exceedance(results["utc"], ax=axes[0, 1]) + plot_exceedance(results["power"], ax=axes[1, 0]) + plot_exceedance(results["npv"], ax=axes[1, 1]) + axes[1, 1].axhline( + results["pos"], ls="--", c="tab:orange", label=f"POS: {results['pos']:.1f}%" + ) # add the probability of success axes[1, 1].axvline(0, ls="--", c="tab:orange") # add the probability of success axes[1, 1].legend() plt.savefig(output_data_path / "exceedance_probabilities.png") +def test_example9(): + from pathlib import Path + + import numpy as np + import xarray as xr + from matplotlib import pyplot as plt + + from pythermogis import calculate_doublet_performance + + # output directory to write output to + output_data_path = Path(__file__).parent / "resources" / "test_output" / "example9" + output_data_path.mkdir(parents=True, exist_ok=True) + + # define the poro-depth and the poro-perm relationships: + # porosity = a * exp(-b * depth) + c + # -with depth in km + poro_a = 28.97 + poro_b = 0.21 + poro_c = 0.0 + poro_std = 1 + + # ln(permeability) = a + b * poro + perm_a = 0.334 + perm_b = -1.655 + perm_std = 1.556 + + # define range of depth interest: + n_samples = 500 + n_depths = 10 + + depths = np.linspace(1e3, 5e3, n_depths) + thickness = 100 + ntg = 1.0 + + porosity = poro_a * np.exp(-poro_b * (depths * 1e-3)) + poro_c + porosity = np.random.normal(porosity[:, None], poro_std, size=(n_depths, n_samples)) + + permeability = np.exp(perm_a * porosity + perm_b) + permeability = np.random.normal(permeability, perm_std) + permeability = np.clip(permeability, 1.0, 1e4) + + reservoir_properties = xr.Dataset( + { + "thickness": thickness, + "porosity": (["depth", "samples"], porosity / 100), + "ntg": ntg, + "permeability": (["depth", "samples"], permeability), + }, + coords={ + "samples": np.arange(n_samples), + "depth": depths, + }, + ) + + # For every sample, run a doublet simulation store the output values + simulations = calculate_doublet_performance( + reservoir_properties, + chunk_size=100, + ) + simulations_stoch = simulations.quantile([0.1,0.5,0.9], dim="samples") + print(simulations) + print(simulations_stoch) + + # make a plot showing there is a sweet spot between Doublet power, Cost of energy + # and Transmissivity + fig, ax = plt.subplots(figsize=(10, 6)) + ax2 = ax.twiny() + + p10 = simulations_stoch.sel(quantile=0.9) + p50 = simulations_stoch.sel(quantile=0.5) + p90 = simulations_stoch.sel(quantile=0.1) + + ax2.fill_betweenx(p50.depth, p10.utc, p90.utc, color="lightgrey", alpha=0.5) + ax2.plot(p50.utc, p50.depth, label="UTC [€/kWh]", color="lightgrey") + + ax2.fill_betweenx(p50.depth, p10.power, p90.power, color="black", alpha=0.5) + ax2.plot(p50.power, p50.power.depth, label="Power [MW]", color="black") + + ax.fill_betweenx(p50.depth, p10.transmissivity_with_ntg, p90.transmissivity_with_ntg, color="tab:green", alpha=0.5) + ax.plot(p50.transmissivity_with_ntg, p50.transmissivity_with_ntg.depth, label="Transmissivity [Dm]", color="tab:green") + + ax.set_ylabel("Depth [m]") + ax2.set_xlabel("Power [MW] & UTC [€/kWh]") + ax.set_xlabel("Transmissivity [Dm]", color="tab:green") + ax.set_ylim([depths[-1], 0.0]) + plt.show() + + + + -- GitLab From 46afc6e90afe8dcb9f1e8e6bea577d920fe1980f Mon Sep 17 00:00:00 2001 From: bretth Date: Fri, 21 Nov 2025 15:49:02 +0100 Subject: [PATCH 2/3] fixing the example --- tests/test_doc_examples.py | 52 +++++++++++++++++++++++++------------- 1 file changed, 35 insertions(+), 17 deletions(-) diff --git a/tests/test_doc_examples.py b/tests/test_doc_examples.py index b3761bf..4b108f3 100644 --- a/tests/test_doc_examples.py +++ b/tests/test_doc_examples.py @@ -483,19 +483,16 @@ def test_example9(): output_data_path.mkdir(parents=True, exist_ok=True) # define the poro-depth and the poro-perm relationships: - # porosity = a * exp(-b * depth) + c - # -with depth in km + # porosity = a * exp(-b * depth), with depth in km poro_a = 28.97 poro_b = 0.21 - poro_c = 0.0 - poro_std = 1 + poro_std = 5 # 1% uncertainty - # ln(permeability) = a + b * poro + # ln(permeability) = a * poro + b perm_a = 0.334 perm_b = -1.655 - perm_std = 1.556 - # define range of depth interest: + # number of depths and number of samples per depth n_samples = 500 n_depths = 10 @@ -503,12 +500,13 @@ def test_example9(): thickness = 100 ntg = 1.0 - porosity = poro_a * np.exp(-poro_b * (depths * 1e-3)) + poro_c + porosity = poro_a * np.exp(-poro_b * (depths * 1e-3)) + + # expand porosity to have n_samples per depth, with uncertainty in porosity porosity = np.random.normal(porosity[:, None], poro_std, size=(n_depths, n_samples)) + porosity = np.clip(porosity, 0.0, 100.0) permeability = np.exp(perm_a * porosity + perm_b) - permeability = np.random.normal(permeability, perm_std) - permeability = np.clip(permeability, 1.0, 1e4) reservoir_properties = xr.Dataset( { @@ -523,14 +521,12 @@ def test_example9(): }, ) - # For every sample, run a doublet simulation store the output values + # For every sample, run a doublet simulation simulations = calculate_doublet_performance( reservoir_properties, chunk_size=100, ) simulations_stoch = simulations.quantile([0.1,0.5,0.9], dim="samples") - print(simulations) - print(simulations_stoch) # make a plot showing there is a sweet spot between Doublet power, Cost of energy # and Transmissivity @@ -541,20 +537,42 @@ def test_example9(): p50 = simulations_stoch.sel(quantile=0.5) p90 = simulations_stoch.sel(quantile=0.1) + # UTC ax2.fill_betweenx(p50.depth, p10.utc, p90.utc, color="lightgrey", alpha=0.5) ax2.plot(p50.utc, p50.depth, label="UTC [€/kWh]", color="lightgrey") + # power ax2.fill_betweenx(p50.depth, p10.power, p90.power, color="black", alpha=0.5) ax2.plot(p50.power, p50.power.depth, label="Power [MW]", color="black") - ax.fill_betweenx(p50.depth, p10.transmissivity_with_ntg, p90.transmissivity_with_ntg, color="tab:green", alpha=0.5) - ax.plot(p50.transmissivity_with_ntg, p50.transmissivity_with_ntg.depth, label="Transmissivity [Dm]", color="tab:green") + # transmissivity Dm + ax.fill_betweenx( + p50.depth, + p10.transmissivity_with_ntg, + p90.transmissivity_with_ntg, + color="tab:green", + alpha=0.5, + ) + ax.plot( + p50.transmissivity_with_ntg, + p50.transmissivity_with_ntg.depth, + label="Transmissivity [Dm]", + color="tab:green", + ) ax.set_ylabel("Depth [m]") ax2.set_xlabel("Power [MW] & UTC [€/kWh]") ax.set_xlabel("Transmissivity [Dm]", color="tab:green") - ax.set_ylim([depths[-1], 0.0]) - plt.show() + ax.set_ylim([depths[-1], depths[0]]) + + # Collect handles and labels from both axes + handles1, labels1 = ax.get_legend_handles_labels() + handles2, labels2 = ax2.get_legend_handles_labels() + handles = handles1 + handles2 + labels = labels1 + labels2 + ax.legend(handles, labels) + + plt.savefig(output_data_path / "power_utc_transmissivity.png") -- GitLab From 8bf8e80dc62263af2453dfaa764ec8b7d7c98cb2 Mon Sep 17 00:00:00 2001 From: bretth Date: Wed, 26 Nov 2025 15:05:21 +0100 Subject: [PATCH 3/3] fixing the example --- tests/test_doc_examples.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_doc_examples.py b/tests/test_doc_examples.py index 4b108f3..eabbf46 100644 --- a/tests/test_doc_examples.py +++ b/tests/test_doc_examples.py @@ -486,7 +486,7 @@ def test_example9(): # porosity = a * exp(-b * depth), with depth in km poro_a = 28.97 poro_b = 0.21 - poro_std = 5 # 1% uncertainty + poro_std = 1 # 1% uncertainty, if this goes larger it blows up massively... # ln(permeability) = a * poro + b perm_a = 0.334 @@ -545,7 +545,7 @@ def test_example9(): ax2.fill_betweenx(p50.depth, p10.power, p90.power, color="black", alpha=0.5) ax2.plot(p50.power, p50.power.depth, label="Power [MW]", color="black") - # transmissivity Dm + # transmissivity ax.fill_betweenx( p50.depth, p10.transmissivity_with_ntg, -- GitLab