From fd5c2e91c17cfa6742f24f9ccc1044fa96dd859c Mon Sep 17 00:00:00 2001 From: korevaarzer Date: Wed, 26 Aug 2026 13:07:30 +0200 Subject: [PATCH 1/9] fix: give each MC sample independent copy of locked list params --- src/geoloop/utils/RunManager.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/geoloop/utils/RunManager.py b/src/geoloop/utils/RunManager.py index d218991..377b228 100644 --- a/src/geoloop/utils/RunManager.py +++ b/src/geoloop/utils/RunManager.py @@ -81,10 +81,15 @@ def sample_parameter_space(config: SingleRunConfig) -> tuple[dict, list[str]]: sample_space[param_name] = np.full(n_samples, float(fixed_value)) # broadcast locked (non-varying) parameters + # Use list() to ensure each sample gets an independent copy of list parameters, + # preventing shared-reference issues when optimization modifies values in-place + # (e.g. r_out[-1] for COAXIAL pipes) for param_name in locked_param_names: - sample_space[param_name] = [ - getattr(config, param_name) for _ in range(n_samples) - ] + val = getattr(config, param_name) + if isinstance(val, list): + sample_space[param_name] = [list(val) for _ in range(n_samples)] + else: + sample_space[param_name] = [val for _ in range(n_samples)] return sample_space, varying_parameters -- GitLab From 6fac50f8377914c26e3ad01cb7b8b594501847d7 Mon Sep 17 00:00:00 2001 From: korevaarzer Date: Wed, 26 Aug 2026 13:08:04 +0200 Subject: [PATCH 2/9] fix: correct r_in/r_out mapping and guard Reff division by zero --- src/geoloop/geoloopcore/CoaxialPipe.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/geoloop/geoloopcore/CoaxialPipe.py b/src/geoloop/geoloopcore/CoaxialPipe.py index efec838..044c35b 100644 --- a/src/geoloop/geoloopcore/CoaxialPipe.py +++ b/src/geoloop/geoloopcore/CoaxialPipe.py @@ -198,8 +198,11 @@ class CoaxialPipe(CustomPipe): R_f_out_out = 1.0 / (h_f_a_out * 2 * np.pi * r_out_in) R_fp = R_p_out + R_f_out_out - r_inner = np.roll(self.r_in, 1) - r_outer = np.roll(self.r_out, 1) + # r_inner = np.roll(self.r_in, 1) + # r_outer = np.roll(self.r_out, 1) + + r_inner = self.r_in + r_outer = self.r_out self.coaxial = gt.pipes.Coaxial( (0, 0), r_inner, r_outer, self.b, k_s, k_g[0], R_ff, R_fp, J=self.J @@ -498,6 +501,11 @@ class CoaxialPipe(CustomPipe): if abs(np.sum(qbz) - power) > 1: print("power is not sumq (sumq, power):", np.sum(qbz), ",", power) - Reff = -np.average(Tf - Tb) * bh.H / power + # Reff = -np.average(Tf - Tb) * bh.H / (power + 1e-6) + + if abs(power) > 1e-3: + Reff = -np.average(Tf - Tb) * bh.H / power + else: + Reff = -1 return T_f_out, power, Reff, ptemp, Tb, qbz -- GitLab From 6b57e63acb7f575990c582fd81d10c26e58a4255 Mon Sep 17 00:00:00 2001 From: korevaarzer Date: Wed, 26 Aug 2026 13:08:25 +0200 Subject: [PATCH 3/9] fix: use Nelder-Mead optimizer with fatol=5 instead of Powell/ftol --- src/geoloop/bin/SingleRunSim.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/geoloop/bin/SingleRunSim.py b/src/geoloop/bin/SingleRunSim.py index 0be7569..6f9aa3d 100644 --- a/src/geoloop/bin/SingleRunSim.py +++ b/src/geoloop/bin/SingleRunSim.py @@ -61,7 +61,7 @@ def optimize_forkeys( x0[i] = 0.5 * sum(boundsval[i]) method = "Nelder-Mead" - method = "Powell" + # method = "Powell" # ftol is tolerance for output power error of optimized point result = minimize( @@ -70,7 +70,7 @@ def optimize_forkeys( args=(optimize_keys, config, copcrit, isample), bounds=boundsval, method=method, - options={"maxiter": 50, "ftol": 100}, + options={"maxiter": 50, "fatol": 5}, ) return result, config -- GitLab From 5740ee0430823f0b2f15143c06ac562d718ebaaa Mon Sep 17 00:00:00 2001 From: korevaarzer Date: Wed, 26 Aug 2026 13:09:05 +0200 Subject: [PATCH 4/9] feat: store r_out as per-sample output when varied by optimization --- src/geoloop/utils/helpers.py | 35 ++++++++++++++++++++++++++++++----- 1 file changed, 30 insertions(+), 5 deletions(-) diff --git a/src/geoloop/utils/helpers.py b/src/geoloop/utils/helpers.py index 6fc93ec..3c22f4e 100644 --- a/src/geoloop/utils/helpers.py +++ b/src/geoloop/utils/helpers.py @@ -625,6 +625,16 @@ def save_MCrun_results( variable_param_dict = {key: param_dict[key] for key in variable_param_names} locked_param_dict = {key: param_dict[key] for key in locked_param_names} + # Check if r_out varies across samples (e.g. due to optimization) + # by computing the std of the last pipe value across all samples + r_out_varies = False + if "r_out" in locked_param_dict: + r_out_last_values = np.array( + [np.asarray(locked_param_dict["r_out"][s])[-1] for s in range(n_samples)] + ) + if round(float(np.std(r_out_last_values)), 7) != 0: + r_out_varies = True + # Variable parameters: Create dataset with [param, samples] dimensions variable_param_da = xr.DataArray( list(variable_param_dict.values()), @@ -662,16 +672,31 @@ def save_MCrun_results( } pipe_pos = np.asarray(locked_param_x["pos"]) - r_out = locked_param_x["r_out"] pipe_pos_ds = xr.Dataset( data_vars=dict(pipe_pos=(["nPipes", "xy"], pipe_pos)), coords=dict(nPipes=(["nPipes"], pipes_coord)), ) - r_out_ds = xr.Dataset( - data_vars=dict(r_out=(["nPipes"], r_out)), - coords=dict(nPipes=(["nPipes"], pipes_coord)), - ) + + # Store r_out: if it varies across samples (optimization), store with [samples, nPipes] dims + # otherwise store as locked with [nPipes] dim only + if r_out_varies: + r_out_per_sample = np.array( + [np.asarray(locked_param_dict["r_out"][s]) for s in range(n_samples)] + ) + r_out_ds = xr.Dataset( + data_vars=dict(r_out=(["samples", "nPipes"], r_out_per_sample)), + coords=dict(samples=(["samples"], list(sample_coord)), nPipes=(["nPipes"], pipes_coord)), + ) + else: + r_out = locked_param_x.get("r_out", locked_param_single_dict.get("r_out")) + if r_out is not None and isinstance(r_out, list): + r_out_ds = xr.Dataset( + data_vars=dict(r_out=(["nPipes"], r_out)), + coords=dict(nPipes=(["nPipes"], pipes_coord)), + ) + else: + r_out_ds = None # Then, create dataset from locked parameters with dim(len(parnames), nlayers_k_s) locked_param_z_k_s = { -- GitLab From 3b46ce5b76eb330c6cc830eeb31d97d54ebf9a1d Mon Sep 17 00:00:00 2001 From: korevaarzer Date: Wed, 26 Aug 2026 13:09:33 +0200 Subject: [PATCH 5/9] feat: add r_out entries to units_dict and format_dict --- src/geoloop/constants.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/geoloop/constants.py b/src/geoloop/constants.py index c33539f..4cf84c2 100644 --- a/src/geoloop/constants.py +++ b/src/geoloop/constants.py @@ -35,6 +35,7 @@ units_dict = { "insu_k": "W/mK", "insu_dr": "%", "k_g_res": "W/mK", + "r_out": "m", } # Dictionary of formatted names for various parameters in the model. @@ -74,6 +75,7 @@ format_dict = { "insu_k": "$insu_{k}$", "insu_dr": "$insu_{dr}$", "k_g_res": "$k_{g-res}$", + "r_out": "$r_{out}$", } # Dictionary with colors for lithology plotting -- GitLab From a97403f09f08762607e19457337411664b68691d Mon Sep 17 00:00:00 2001 From: korevaarzer Date: Wed, 26 Aug 2026 13:10:16 +0200 Subject: [PATCH 6/9] feat: add r_out scatter plots to stochastic output --- src/geoloop/plotting/create_plots.py | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/src/geoloop/plotting/create_plots.py b/src/geoloop/plotting/create_plots.py index 937132c..cb966d5 100644 --- a/src/geoloop/plotting/create_plots.py +++ b/src/geoloop/plotting/create_plots.py @@ -84,16 +84,23 @@ class PlotResults: ] variable_param_names, locked_param_names = get_param_names() - x_variables = result_variables + variable_param_names + x_variables = result_variables + variable_param_names + ["r_out"] # Loop through variables and create scatterplots for x_variable in x_variables: if x_variable != y_variable and x_variable in results_df.columns: # only plot the variables that are truly variable (stdev != 0) if results_df[x_variable].std().round(7) != 0: - title = f"{format_dict[x_variable]} ({units_dict[x_variable]}) after {results_df.loc[0, 'hours']} hours production vs. {format_dict[y_variable]} ({units_dict[y_variable]})" + # Handle r_out unit conversion + plot_df = results_df.copy() + x_unit = units_dict[x_variable] + if x_variable == "r_out": + plot_df[x_variable] = plot_df[x_variable] * 1000 + x_unit = "mm" + + title = f"{format_dict[x_variable]} ({x_unit}) after {plot_df.loc[0, 'hours']} hours production vs. {format_dict[y_variable]} ({units_dict[y_variable]})" g = sns.JointGrid( - results_df, x=x_variable, y=y_data[y_variable], hue="file_name" + plot_df, x=x_variable, y=y_data[y_variable], hue="file_name" ) g.plot_joint( sns.scatterplot, alpha=0.7, edgecolor=".2", linewidth=0.5 @@ -118,7 +125,7 @@ class PlotResults: ) g.ax_joint.legend() # Ensure the legend is updated g.set_axis_labels( - xlabel=f"{format_dict[x_variable]} ({units_dict[x_variable]})", + xlabel=f"{format_dict[x_variable]} ({x_unit})", ylabel=f"{format_dict[y_variable]} ({units_dict[y_variable]})", ) g.fig.suptitle(title) @@ -135,9 +142,16 @@ class PlotResults: in params_df.select_dtypes(include=["float64", "int64"]).columns ): if params_df[x_variable].std().round(7) != 0: - title = f"{format_dict[x_variable]} ({units_dict[x_variable]}) after {params_df.loc[0, 'nyear']} year production vs. {format_dict[y_variable]} ({units_dict[y_variable]})" + # Handle r_out unit conversion + plot_params_df = params_df.copy() + x_unit = units_dict[x_variable] + if x_variable == "r_out": + plot_params_df[x_variable] = plot_params_df[x_variable] * 1000 + x_unit = "mm" + + title = f"{format_dict[x_variable]} ({x_unit}) after {plot_params_df.loc[0, 'nyear']} year production vs. {format_dict[y_variable]} ({units_dict[y_variable]})" g = sns.JointGrid( - params_df, x=x_variable, y=y_data[y_variable], hue="file_name" + plot_params_df, x=x_variable, y=y_data[y_variable], hue="file_name" ) g.plot_joint( sns.scatterplot, alpha=0.7, edgecolor=".2", linewidth=0.5 @@ -156,7 +170,7 @@ class PlotResults: g.ax_joint.invert_yaxis() g.set_axis_labels( - xlabel=f"{format_dict[x_variable]} ({units_dict[x_variable]})", + xlabel=f"{format_dict[x_variable]} ({x_unit})", ylabel=f"{format_dict[y_variable]} ({units_dict[y_variable]})", ) g.fig.suptitle(title) -- GitLab From 062936cc5dc4c06fb16f0023ca0bc7e3f9cf623e Mon Sep 17 00:00:00 2001 From: korevaarzer Date: Wed, 26 Aug 2026 13:12:01 +0200 Subject: [PATCH 7/9] docs: mention r_out optimization and per-sample output --- docs/manual/simulation_configuration.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/docs/manual/simulation_configuration.md b/docs/manual/simulation_configuration.md index f3cc2cc..ff4d0a6 100644 --- a/docs/manual/simulation_configuration.md +++ b/docs/manual/simulation_configuration.md @@ -229,6 +229,7 @@ Supported parameters include: - Tgrad - Q - fluid_percent +- r_out Supported sampling distributions and their corresponding arguments for the configuration JSON are: @@ -256,6 +257,10 @@ optimization simulation, see the example for [Optimization of a coaxial BHE](../ > and upwards in the inner pipe. Optimization of the pipe radius for > this design is only supported for the inner pipe. +> When `r_out` is optimized, its resulting value can differ between samples. +> In that case, `r_out` is stored per-sample (dimensions `samples`; `nPipes`) +> in the results file instead of the usual single locked-parameter value. + --- ## Simulation results -- GitLab From 972a90fc13ac3c67259ea23f0e9c4b9da6c70f84 Mon Sep 17 00:00:00 2001 From: korevaarzer Date: Wed, 26 Aug 2026 13:12:15 +0200 Subject: [PATCH 8/9] chore: bump version to 1.1.0 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index b586f53..c653d0b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "geoloop" -version = "1.0.7" +version = "1.1.0" description = "This is a Python package for simulation of (deep) vertical Borehole Heat Exchanger (BHE) systems" authors = [ { name = "Zanne Korevaar", email = "zanne.korevaar@tno.nl" }, -- GitLab From 81e0a8c7e09d10aca74291e6870b99c0fb1adc22 Mon Sep 17 00:00:00 2001 From: korevaarzer Date: Wed, 26 Aug 2026 13:29:32 +0200 Subject: [PATCH 9/9] update uv.lock --- uv.lock | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/uv.lock b/uv.lock index 4815ea6..2265c78 100644 --- a/uv.lock +++ b/uv.lock @@ -340,7 +340,7 @@ wheels = [ [[package]] name = "geoloop" -version = "1.0.2" +version = "1.1.0" source = { editable = "." } dependencies = [ { name = "h5netcdf" }, -- GitLab