TNO Intern

Commit c90a084a authored by Florian Knappers's avatar Florian Knappers
Browse files

added loop over p values

parent a7e46c4f

src/pythermogis/workflow/utc/utc.py

+52 −72
Original line number Diff line number Diff line
@@ -51,7 +51,7 @@ def run_utc_workflow(config: UTCConfiguration) -> xr.DataTree: 
        start = time.perf_counter()



        aquifer_ds = root[aquifer].to_dataset()

        updated_ds, n_cells = compute_results_for_aquifer(

        updated_ds = compute_results_for_aquifer(

            aquifer_ds, aquifer, config, temperature_model

        )

        root[aquifer] = xr.DataTree(dataset=updated_ds, name=aquifer)
@@ -59,12 +59,7 @@ def run_utc_workflow(config: UTCConfiguration) -> xr.DataTree: 
        end = time.perf_counter()

        elapsed = end - start



        cells_per_sec = n_cells / elapsed



        print(

            f"Time required for aquifer {aquifer}: {elapsed:.2f} seconds. "

            f"{cells_per_sec} cells per second"

        )

        print(f"Time required for aquifer {aquifer}: {elapsed:.2f} seconds.")



    if config.results_dir != "":

        zarr_path = f"{config.results_dir}/ROSL_ROSLU_test.zarr"
@@ -149,15 +144,13 @@ def compute_results_for_aquifer( 
            tuple(result) if result is not None else tuple(np.nan for _ in output_names)

        )



    # TODO: loop over pvalues and save grids with pvalue noted

    # TODO: Maybe save the results as datasets iso dataarrays & have p value be a dim?



    for p_value in config.p_values:

        # calc transmissivity(_with_ntg)

    hydro = compute_hydro_properties_mc(aquifer_ds, 90)

    transmissivity = hydro["transmissivity"]

    transmissivity_with_ntg = hydro["transmissivity_with_ntg"]

    aquifer_ds["transmissivity"] = transmissivity

    aquifer_ds["transmissivity_with_ntg"] = transmissivity_with_ntg

        hydro = compute_hydro_properties_mc(aquifer_ds, p_value)

        aquifer_ds["transmissivity"] = hydro["transmissivity"]

        aquifer_ds["transmissivity_with_ntg"] = hydro["transmissivity_with_ntg"]

        aquifer_ds["thickness"] = hydro["thickness"]

        aquifer_ds["permeability"] = hydro["permeability"]



        # calc temperature grid

        mid_depth = -(aquifer_ds["depth"] + 0.5 * aquifer_ds["thickness"])
@@ -168,19 +161,6 @@ def compute_results_for_aquifer( 


        aquifer_ds["temperature"] = temperature



    valid_mask = (

        aquifer_ds["thickness"].notnull()

        & aquifer_ds["transmissivity"].notnull()

        & aquifer_ds["transmissivity_with_ntg"].notnull()

        & aquifer_ds["ntg"].notnull()

        & aquifer_ds["depth"].notnull()

        & aquifer_ds["porosity"].notnull()

        & aquifer_ds["temperature"].notnull()

        & (aquifer_ds["transmissivity_with_ntg"] >= config.kh_cutoff)

    )



    n_cells = int(valid_mask.sum().compute())



        # use hc_accum mask..

        hc_mask = aquifer_ds["hc_accum"] != 0

        vars_to_mask = [
@@ -213,13 +193,13 @@ def compute_results_for_aquifer( 
        )



        for name, result in zip(output_names, results, strict=True):

        aquifer_ds[name] = result

            aquifer_ds[f"{name}_p{p_value}"] = result



    aquifer_ds.compute()



    print(f"Computed results for {aquifer_name}")



    return aquifer_ds, n_cells

    return aquifer_ds





def compute_hydro_properties_mc(

tests/utc/test_rosl_roslu.py

+2 −1
Original line number Diff line number Diff line
@@ -14,10 +14,11 @@ def test_rosl_roslu(): 
        results_dir="C:/Users/knappersfy/work/thermogis/pythermogis/workdir/ROSL_ROSLU/output",

        aquifers=["ROSL_ROSLU"],

        segment_length=1,

        p_values=[95]

    )

    result = run_utc_workflow(config)



    result.ROSL_ROSLU.utc.plot()

    result.ROSL_ROSLU.utc_p95.plot()

    plt.show()



    print(result)
 
 No newline at end of file