TNO Intern

Commit 86a6e171 authored by Hen Brett's avatar Hen Brett 🐔
Browse files

Tried doing some speedups

parent 89b16f61

src/pythermogis/thermogis_classes/doublet.py

+6 −8
Original line number Diff line number Diff line
@@ -6,9 +6,6 @@ def simulate_doublet(output_data: xr.Dataset, reservoir_properties: xr.Dataset, 
    # Calculate transmissivity scaled by ntg and converted to Dm

    output_data[f"transmissivity_with_ntg"] = (output_data[f"transmissivity"] * reservoir_properties.ntg) / 1e3



    # Instantiate ThermoGIS doublet

    doublet = instantiate_thermogis_doublet(utc_properties, rng_seed)



    # Calculate the doublet performance across all dimensions

    output_data_arrays = xr.apply_ufunc(calculate_performance_of_single_location,

                                        reservoir_properties.mask,
@@ -19,11 +16,11 @@ def simulate_doublet(output_data: xr.Dataset, reservoir_properties: xr.Dataset, 
                                        output_data.temperature,

                                        output_data.transmissivity,

                                        output_data.transmissivity_with_ntg,

                                        kwargs={"doublet": doublet, "utc_properties": utc_properties},

                                        input_core_dims=[[], [], [], [], [], [], [], []],

                                        rng_seed,

                                        kwargs={"utc_properties": utc_properties},

                                        input_core_dims=[[], [], [], [], [], [], [], [],[]],

                                        output_core_dims=[[], [], [], [], [], [], [], [], [], [], [], [], [], []],

                                        vectorize=True,

                                        dask="parallelized",

                                        )



    # Assign output DataArrays to the output_data object
@@ -44,8 +41,7 @@ def simulate_doublet(output_data: xr.Dataset, reservoir_properties: xr.Dataset, 
    return output_data





def calculate_performance_of_single_location(mask: float, depth: float, thickness: float, porosity: float, ntg: float, temperature: float, transmissivity: float, transmissivity_with_ntg: float, doublet: JClass = None ,

                                             utc_properties: JClass = None) -> float:

def calculate_performance_of_single_location(mask: float, depth: float, thickness: float, porosity: float, ntg: float, temperature: float, transmissivity: float, transmissivity_with_ntg: float, rng_seed: int, utc_properties: JClass = None) -> float:

    """

    Calculate the performance of a doublet at a single location.
@@ -96,6 +92,8 @@ def calculate_performance_of_single_location(mask: float, depth: float, thicknes 
    if not np.isnan(mask):

        return 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0



    # Instantiate ThermoGIS doublet

    doublet = instantiate_thermogis_doublet(utc_properties, rng_seed)

    set_doublet_parameters(doublet, transmissivity_with_ntg, depth, porosity, ntg, temperature, utc_properties)



    # The Java routine which calculates DoubletPerformance, for more detail on the simulation inspect the Java source code

tests/test_doc_examples.py

+30 −30
Original line number Diff line number Diff line
@@ -282,7 +282,7 @@ def test_example8(): 
    output_data_path.mkdir(parents=True, exist_ok=True)



    # generate simulation samples across desired reservoir properties

    Nsamples = 1000

    Nsamples = 100

    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)
@@ -309,34 +309,34 @@ def test_example8(): 
    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))

    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)

    plt.savefig(output_data_path / "reservoir_distributions.png")



    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])

    plt.savefig(output_data_path / "simulation_results.png")



    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

    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")

    # # plotting

    # 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)

    # plt.savefig(output_data_path / "reservoir_distributions.png")

    #

    # 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])

    # plt.savefig(output_data_path / "simulation_results.png")

    #

    # 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

    # 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")