Hi Wouter, I'm a big fan of GHETool and the flexible Python backend! 😊

While simulating energy piles with monthly loads, I ran into the same non-physical “wiggles” in the gfunction that you described in MassimoCimmino/pygfunction#309; with resulting effects on the fluid temperatures, especially when the peak_load_duration is smaller than 730 hours.

I think I found a solution to avoid this: It seems tied to the cylindrical correction: in GHETool/VariableClasses/Cylindrical_correction.py Line 224, the CHS is evaluated exactly at the borehole wall, which seems to lead to numerical instabilities somewhere in the CHS:

h_chs = cylindrical_heat_source(time, alpha, r_b, r_b)

A simple (and very effective) workaround was for me to evaluate slightly off the wall, e.g. r = r_b + 0.01 m:

h_chs = cylindrical_heat_source(time, alpha, r_b + 0.01, r_b)

This removes the artifacts in the g-function without any other code changes. The difference in g between using r_b ± 0.01 is in the order of ~1e-5, so the impact seems negligible in practice. I know this is a minimal tweak rather than a principled fix, but it seems robust and the results physically reasonable.

Happy to implement this little change to the code if helpful. Thanks again for the great tool!