NewFeature unequal segment lengths

[j-c-cook]

This builds on #133 and contains the commits from PR #134. This is related to #135 and should not be merged until after #134 is.

The segments function in boreholes.py had previously only accounted for the segment lengths to be equal on a specifi borehole. This modification provides the ability for the segments to have different lengths.

The gfunction.py file needed to address unequal segment lengths in the _map_axial_segment_pairs function so that the oreholes.segments function could properly create bore segments. The _baseSolver object had updates to the segment_lengths(), borehole_segments() and _check_inputs() functions.

• The summation of the segment lengths for each borehole are asserted to be within a tenth of a meter of the borehole height
• The borehole_segments() function accesses the new variable segmentLengths (if its not None) to split boreholes into segments
• The segment_lengths() function flattens the 2D segmentLengths list if segmentLengths is not None

[j-c-cook]

@MassimoCimmino I think this is ready to go. Can you take a look?

[MassimoCimmino]

I see potentially similar issues as in #133 / #134:

• The distribution of segment lengths should be compared when similarities are identified (like the number of segments are now checked);
• This might affect the visualization methods of the gFunction class.

I would suggest a list of arrays of segmentLengthRatios instead of segmentLengths to avoid the user have to regenerate the list any time the borehole length changes. I would expect any recommendation on the "optimal" discretization scheme to be in the form of ratios rather than lengths.

I will move on to this PR after #132, #138 and #146.

Also to be considered here is that the functionality for unequal segment lengths needs to be added to the pipe and network classes to preserve the full functionality of pygfunction. Depending on how complex the implementation is, we can consider the PR for v2.1 or v2.2.

[MassimoCimmino]

I will be working on this PR on my branch.

I have implemented segment length ratios, updated visualization methods, and fixed a few issues related to 'similarities' and the construction of the thermal response factor matrices. I will now be moving on to the pipe and network classes.

[j-c-cook]

Okay sounds good, I'm glad you're open for this being in v2.1, but I'll understand if it gets pushed back to v2.2. How is that looking? I've been busy with GLHEDT and it looks like I will be for the next couple months as I try to push out a report on the subject.

It looks like you are working hard and doing some rapid development. The progress is exciting. Thank you.

[MassimoCimmino]

I have also been busy with the start of the Fall semester. I have been looking into the pipes and networks module and it doesn't look like it should be too difficult to implement the variable segment lengths.

I am excited with the results so far with the 'UBWT' condition. It seems unequal segment lengths can get more accurate results faster and with fewer segments. I am looking forward to have this in v2.1. It should also be compatible with the recent 'equivalent' solver from #146.

[MassimoCimmino]

Thanks @j-c-cook. I just pushed my implementation. I do prefer your suggestion of imposing the end_length_ratio. I will update the examples and tests, and then do a final review of the PR.

For the optimal discretization parameters, I think this is a topic for another issue.

[MassimoCimmino]

Here is a first set of comparisons to verify the results of the current PR with the master branch. We also compare the accuracy and calculation time of the discretization method against a uniform discretization. Calculations are done with the 'similarities' solver.

Rectangular borefields

Comparison of g-functions using 12 segments between the PR and the master branch

It is shown here that there is no significant differences in g-function calculations when compared to the master.

Comparison of g-functions with a high-resolution solution (nSegments=64) from the master

It is shown here that the new discretization achieves better accuracy with fewer nSegments when compared to a uniform discretization.

Calculation time

This better accuracy also comes with a reduced calculation time. There is a small (but not so significant) increase in calculation time using uniform segments in the PR.

Field of 100 boreholes with two different lengths

Results also show better accuracy with lower calculation time in this case.

[MassimoCimmino]

Here is a second set of comparisons using the 'equivalent' solver. All cases are the same. Calculations are done with the 'equivalent' solver, except for the reference results using nSegments=64 which are generated using the 'similarities solver'.

Rectangular borefields

Comparison of g-functions using 12 segments between the PR and the master branch

It is shown here that there is no significant differences in g-function calculations when compared to the master.

Comparison of g-functions with a high-resolution solution (nSegments=64) from the master

It is shown here that the new discretization achieves better accuracy with fewer nSegments when compared to a uniform discretization.

Calculation time

This better accuracy also comes with a reduced calculation time. There is an increase in calculation time using uniform segments in the PR.

Field of 100 boreholes with two different lengths

Results also show better accuracy with lower calculation time in this case.

[MassimoCimmino]

Results are updated to 32394b6. This seems to solve the issues related to calculation time.

[j-c-cook]

@MassimoCimmino When you say

solve the issues related to calculation time

do you mean that for the same number of segments in the g-function calculation that unequal segment lengths and equal segment lengths now take the same amount of calculation time?

[MassimoCimmino]

@j-c-cook Not exactly. We saw an increase in the calculation time before 32394b6 when compared to the master for the same number of segments with a uniform discretization. This was much more visible in the equivalent solver. See the figures below for results before 32394b6.

For the calculation time when comparing uniform and non-uniform discretizations, the increase in calculation time is simply due to more thermal response factors being evaluated (due to loss of similarities). As an example, for a pair of boreholes with 12 uniform segments, we have 78 FLS evaluations. For 12 non-uniform segments, it increases to 138. There are 60 evaluations for 8 non-uniform segments.

Similarities (before fix)

Equivalent (before fix)

[j-c-cook]

That equivalence method breaks the performance barrier for computing g-functions. Really incredible work by Prieto and Cimmino (2021).

[j-c-cook]

@MassimoCimmino I've looked closer at #136 (comment) and it appears that this PR is still slower for equivalent methods with 12 uniform segments. Is that an updated plot, or an old one before your fix?

[MassimoCimmino]

#136 (comment) shows the old results. I updated the comment to clarify.
#136 (comment) and #136 (comment) show the latest results.

[j-c-cook]

I just ran 28 cases ranging from 2x2 to 26x26 with the equivalent solver combined with the unequal segments introduced in this PR. All of the cases ran fine, except for one.

I'm not sure how its possible that just one case can break. I also ran these exact cases yesterday and they all ran fine. Can you take a look at the line that threw the error?

Edit: This was run on a high performance computing cluster, and could very well be a fluke. Sometimes those computers throw strange errors that don't make sense, but I think its worth considering nonetheless.

Error Message

  File "pyg2_uift.py", line 156, in <module>
    main(sys.argv)
  File "pyg2_uift.py", line 128, in main
    options=options, method='equivalent')
  File "/scratch/jackcook/Masters-Thesis/pygfunction/v07/pygfunction/pygfunctionv2_0/gfunction.py", line 202, in __init__
    self.gFunc = self.evaluate_g_function(self.time)
  File "/scratch/jackcook/Masters-Thesis/pygfunction/v07/pygfunction/pygfunctionv2_0/gfunction.py", line 233, in evaluate_g_function
    self.gFunc = self.solver.solve(time, self.alpha)
  File "/scratch/jackcook/Masters-Thesis/pygfunction/v07/pygfunction/pygfunctionv2_0/gfunction.py", line 1428, in solve
    segment_ratios=self.segment_ratios)
  File "/scratch/jackcook/Masters-Thesis/pygfunction/v07/pygfunction/pygfunctionv2_0/networks.py", line 763, in coefficients_borehole_heat_extraction_rate
    m_flow_network, cp_f, nSegments, segment_ratios=segment_ratios)
  File "/scratch/jackcook/Masters-Thesis/pygfunction/v07/pygfunction/pygfunctionv2_0/networks.py", line 483, in coefficients_inlet_temperature
    for i in range(self.nBoreholes)]))
  File "/scratch/jackcook/Masters-Thesis/pygfunction/v07/pygfunction/pygfunctionv2_0/networks.py", line 483, in <listcomp>
    for i in range(self.nBoreholes)]))
  File "/scratch/jackcook/Masters-Thesis/pygfunction/v07/pygfunction/pygfunctionv2_0/pipes.py", line 391, in coefficients_outlet_temperature
    segment_ratios=segment_ratios)
  File "/scratch/jackcook/Masters-Thesis/pygfunction/v07/pygfunction/pygfunctionv2_0/pipes.py", line 1010, in _continuity_condition_base
    a_b = np.zeros((self.nOutlets, nSegments))
TypeError: 'list' object cannot be interpreted as an integer

Proof that all other cases run fine

Input file for cluster that accepts a path to a json file containing inputs

This is basically the same as discretize_boreholes.py.

# Jack C. Cook
# Tuesday, August 31, 2021

# Compute g-functions using pygfunction v2.0

from pygfunction import pygv2 as gt
import sys
import json
import time as tyme
import numpy as np
from numpy import pi


def js_r(filename: str):
    with open(filename) as f_in:
        return json.load(f_in)


def main(argv):
    input_path = argv[1]

    d = js_r(input_path)
    
    # Pipe dimensions
    # ---------------
    r_out = 26.67 / 1000. / 2.  # Pipe outer radius (m)
    r_in = 21.6 / 1000. / 2.  # Pipe inner radius (m)
    s = 32.3 / 1000.  # Inner-tube to inner-tube Shank spacing (m)
    epsilon = 1.0e-6  # Pipe roughness (m)

    pos = [(-0.029484999999999997, 3.610871087285971e-18),
           (0.029484999999999997, -7.221742174571942e-18)]

    m_flow_borehole = 0.2

    # Thermal conductivities
    # ----------------------
    k_p = 0.4  # Pipe thermal conductivity (W/m.K)
    k_s = 2.0  # Ground thermal conductivity (W/m.K)
    k_g = 1.0  # Grout thermal conductivity (W/m.K)

    # Fluid properties
    mixer = 'MEG'  # Ethylene glycol mixed with water
    percent = 0.  # Percentage of ethylene glycol added in
    fluid = gt.media.Fluid(mixer=mixer, percent=percent)

    # Volumetric heat capacities
    # --------------------------
    rhoCp_p = 1542. * 1000.  # Pipe volumetric heat capacity (J/K.m3)
    rhoCp_s = 2343.493 * 1000.  # Soil volumetric heat capacity (J/K.m3)
    rhoCp_g = 3901. * 1000.  # Grout volumetric heat capacity (J/K.m3)
    
    # g-Function calculation options
    nSegments = 12
    options = {'nSegments':nSegments, 'disp':False}

    # nSegments = 8
    # segment_ratios = \
    #     gt.utilities.segment_ratios(nSegments, end_length_ratio=0.02)
    # options = \
    #     {'nSegments': nSegments, 'segment_ratios': segment_ratios,
    #      'disp':False}
    
    alpha = k_s / rhoCp_s

    # -------------------------------------------------------------------------
    # Simulation parameters
    # -------------------------------------------------------------------------

    for i in range(len(d)):
        parameters = d[i]

        # Borehole dimensions
        D = parameters['D'][0]             # Borehole buried depth (m)
        H = parameters['H'][0]           # Borehole length (m)
        r_b = parameters['r_b'][0]         # Borehole radius (m)
        B = parameters['B'][0]             # Borehole spacing (m)

        # g-Function calculation options

        # Geometrically expanding time vector.
        ts = H**2/(9.*alpha)            # Bore field characteristic time
        # time = gt.utilities.time_geometric(dt, tmax, Nt)
        log_time = [-8.5, -7.8, -7.2, -6.5, -5.9, -5.2, -4.5,
                    -3.963, -3.27, -2.864,-2.577, -2.171, -1.884,
                    -1.191, -0.497, -0.274, -0.051, 0.196, 0.419,
                    0.642, 0.873, 1.112, 1.335, 1.679, 2.028,
                    2.275, 3.003]
        time = np.exp(log_time) * ts

        coordinates = parameters['bore_locations']
        boreField = []
        for i in range(len(coordinates)):
            _x = coordinates[i][0]
            _y = coordinates[i][1]
            borehole = gt.boreholes.Borehole(H, D, r_b, _x, _y)
            boreField.append(borehole)

        nbh = len(coordinates)
        nBoreholes = len(boreField)
        
        # Pipe thermal resistance
        R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(r_in, r_out, k_p)
        # Fluid to inner pipe wall thermal resistance (Single U-tube)
        m_flow_pipe = m_flow_borehole
        h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
            m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp, epsilon)
        R_f = 1.0/(h_f*2*pi*r_in)

        # Single U-tube, same for all boreholes in the bore field
        UTubes = []
        for borehole in boreField:
            SingleUTube = gt.pipes.SingleUTube(pos, r_in, r_out,
                                               borehole, k_s, k_g, R_f + R_p)
            UTubes.append(SingleUTube)
        m_flow_network = m_flow_borehole*nBoreholes
        network = gt.networks.Network(
            boreField, UTubes, m_flow_network=m_flow_network, cp_f=fluid.cp,
            nSegments=nSegments)

        nSources = nbh * nSegments
        tic = tyme.time()
        # gfunc = gt.gfunction.gFunction(
        #     network, alpha, time=time, boundary_condition='MIFT',
        #     options=options, method='similarities')
        gfunc = gt.gfunction.gFunction(
            network, alpha, time=time, boundary_condition='MIFT',
            options=options, method='equivalent')
        toc = tyme.time()

        total = toc - tic

        print('nbh\t{}\tH\t{}\tnSegments\t{}\ttime\t{}\tTask\t{}'.
              format(nbh, H, nSegments, total, parameters['name']))

        # Export g-Function to the output folder
        d_out = {'g': {}}
        key = '{}_{}_{}_{}'.format(B, H, r_b, D)
        d_out['g'][key] = gfunc.gFunc.tolist()

        d_out['bore_locations'] = coordinates

        d_out['log_time'] = log_time

        js_out('output/' + parameters['name'], d_out)

    return


def js_out(name, d):
    with open(name + '.json', 'w+') as fp:
        json.dump(d, fp)


if __name__ == '__main__':
    main(sys.argv)

[j-c-cook]

The error is reproducible. The field is 2x2. The attached file throws the error.
recreate_error.zip