Merge branch 'NREL:main' into main

Author: softwareengineerprogrammer

setup.py

@@ -77,7 +77,10 @@ def read(*names, **kwargs):
         'h5py==3.10.0',  # TODO resolve apparent h5py==3.11.0 build compatibility issues on some platforms
         'scipy',
         'iapws',
-        'coolprop',
+        'coolprop; python_version > "3.8"',
+        # coolprop > 6.7.0 fails on Python 3.8, see
+        # https://github.com/softwareengineerprogrammer/GEOPHIRES/actions/runs/14599234458/job/40952924048?pr=69#step:5:302
+        'coolprop==6.7.0; python_version == "3.8"',
         'rich',
         'pylocker',
     ],

src/geophires_x/SBTReservoir.py

@@ -5,16 +5,18 @@
 import pandas as pd
 from scipy.special import erf, erfc, jv, yv, exp1
 from scipy.interpolate import interp1d
+from scipy.integrate import trapezoid
 import scipy.io as sio
 import matplotlib.pyplot as plt
 
+import CoolProp.CoolProp as CP
+
 import geophires_x.Model as Model
 from .CylindricalReservoir import CylindricalReservoir
 from .OptionList import FlowrateModel, InjectionTemperatureModel, Configuration
 from .Parameter import intParameter, floatParameter, OutputParameter, ReadParameter, strParameter, boolParameter
 from .Reservoir import Reservoir
 from .Units import *
-import sys
 from functools import lru_cache
 
 
@@ -378,6 +380,8 @@ def Calculate_Coaxial(self, model):
         """
         model.logger.info(f'Init {str(__class__)}: {sys._getframe().f_code.co_name}')
 
+        raise NotImplementedError('SBT with coaxial configuration is not implemented at this time.')
+
         # Clear all equivalent: Initialize variables and import necessary libraries
 
         # SBT v2 for co-axial heat exchanger with high-temperature capability
@@ -458,28 +462,30 @@ def Calculate_Coaxial(self, model):
         eps_annulus = Dh_annulus * piperoughness  # Relative roughness annulus [-]
         eps_centerpipe = 2 * radiuscenterpipe * piperoughness  # Relative roughness inner pipe [-]
 
+        # These variables are set but never used, and sio.loadmat is causing a "No such file" error.
+        # (see https://github.com/NREL/GEOPHIRES-X/issues/373)
         # Variable fluid properties logic
-        if variablefluidproperties == 0:
-            Pvector = [1, 1e9]
-            Tvector = [1, 1e4]
-            density = np.full((2, 2), rho_f)
-            heatcapacity = np.full((2, 2), cp_f)
-            thermalconductivity = np.full((2, 2), k_f)
-            viscosity = np.full((2, 2), mu_f)
-            thermalexpansion = np.zeros((2, 2))
-        else:
-            print('Loading fluid properties ...')
-            if fluid == 1:
-                # Load properties for water from pre-generated CoolProp data
-                properties = sio.loadmat('properties_H2O_HT_v3.mat')
-                print('Fluid properties for water loaded successfully')
-            elif fluid == 2:
-                # Load properties for CO2 from pre-generated CoolProp data
-                properties = sio.loadmat('properties_CO2.mat')
-                print('Fluid properties for CO2 loaded successfully')
-            else:
-                print('No valid fluid selected')
-                exit()
+        #if variablefluidproperties == 0:
+        #    Pvector = [1, 1e9]
+        #    Tvector = [1, 1e4]
+        #    density = np.full((2, 2), rho_f)
+        #    heatcapacity = np.full((2, 2), cp_f)
+        #    thermalconductivity = np.full((2, 2), k_f)
+        #    viscosity = np.full((2, 2), mu_f)
+        #    thermalexpansion = np.zeros((2, 2))
+        #else:
+        #    print('Loading fluid properties ...')
+        #    if fluid == 1:
+        #        # Load properties for water from pre-generated CoolProp data
+        #        properties = sio.loadmat('properties_H2O_HT_v3.mat')
+        #        print('Fluid properties for water loaded successfully')
+        #    elif fluid == 2:
+        #        # Load properties for CO2 from pre-generated CoolProp data
+        #        properties = sio.loadmat('properties_CO2.mat')
+        #        print('Fluid properties for CO2 loaded successfully')
+        #    else:
+        #        print('No valid fluid selected')
+        #        exit()
 
         # Length of each segment [m]
         Deltaz = np.sqrt(np.diff(x) ** 2 + np.diff(y) ** 2 + np.diff(z) ** 2)
@@ -645,9 +651,10 @@ def Calculate_Coaxial(self, model):
         # MIR N = len(x) - 1
         N = len(x) - 1
         Nboiler = len(boilerelements)
-        # MIR Nreg = N - Nboiler
-        Nreg = 1 + N - Nboiler
+        Nreg = N - Nboiler
+        #Nreg = 1 + N - Nboiler
         self.krock.value_vector = np.full(N, self.krock.value)
+        k_m_vector = np.zeros(N)
         k_m_vector[Nreg:] = k_m_boiler
         alpha_m_vector = k_m_vector / self.rhorock.value / self.cprock.value
 
@@ -662,24 +669,26 @@ def Calculate_Coaxial(self, model):
         SoverLSorted = SMatrixSorted / (np.ones((N, 1)) * Deltaz)
         mindexNPCP = np.argmax(np.min(SoverLSorted, axis=1) < LimitSoverL)
 
-        # MIR midpointsx = 0.5 * x[1:] + 0.5 * x[:-1]
-        # MIR midpointsy = 0.5 * y[1:] + 0.5 * y[:-1]
-        # MIR midpointsz = 0.5 * z[1:] + 0.5 * z[:-1]
-        midpointsx = 0.5 * x + 0.5 * x
-        midpointsy = 0.5 * y + 0.5 * y
-        midpointsz = 0.5 * z + 0.5 * z
+        midpointsx = 0.5 * x[1:] + 0.5 * x[:-1]
+        midpointsy = 0.5 * y[1:] + 0.5 * y[:-1]
+        midpointsz = 0.5 * z[1:] + 0.5 * z[:-1]
+        #midpointsx = 0.5 * x + 0.5 * x
+        #midpointsy = 0.5 * y + 0.5 * y
+        #midpointsz = 0.5 * z + 0.5 * z
         verticalchange = np.diff(z)
         # MIR
-        verticalchange = np.append(verticalchange, verticalchange[-1])
+        #verticalchange = np.append(verticalchange, verticalchange[-1])
 
         if initialtemperatureprofile == 0:
             BBinitial = Tsurf - GeoGradient * midpointsz
             Tfluidupnodes = Tsurf - GeoGradient * z
             Tfluiddownnodes = Tsurf - GeoGradient * z
         elif initialtemperatureprofile == 1:
             BBinitial = np.interp(midpointsz, initialtemperaturedata[:, 0], initialtemperaturedata[:, 1])
-            Tfluidupnodes = np.interp(z, initialtemperaturedata[:, 0], initialtemperaturedata[:, 1])
-            Tfluiddownnodes = np.interp(z, initialtemperaturedata[:, 0], initialtemperaturedata[:, 1])
+            #Tfluidupnodes = np.interp(z, initialtemperaturedata[:, 0], initialtemperaturedata[:, 1])
+            #Tfluiddownnodes = np.interp(z, initialtemperaturedata[:, 0], initialtemperaturedata[:, 1])
+            Tfluidupnodes = np.interp(z[:-1], initialtemperaturedata[:, 0], initialtemperaturedata[:, 1])
+            Tfluiddownnodes = np.interp(z[:-1], initialtemperaturedata[:, 0], initialtemperaturedata[:, 1])
 
         # MIR Tfluiddownmidpoints = 0.5 * Tfluiddownnodes[1:] + 0.5 * Tfluiddownnodes[:-1]
         # MIR Tfluidupmidpoints = 0.5 * Tfluidupnodes[1:] + 0.5 * Tfluidupnodes[:-1]
@@ -806,26 +815,38 @@ def Calculate_Coaxial(self, model):
         self.Tresoutput.value[0] = Tsurf
         Poutput[0] = Pin * 1e5
 
-        Tfluidupnodesstore = np.zeros((N + 1, len(times)))
-        Tfluiddownnodesstore = np.zeros((N + 1, len(times)))
-        # MIR Tfluidupmidpointsstore = np.zeros((N, len(times)))
-        Tfluidupmidpointsstore = np.zeros((N + 1, len(times)))
-        # MIR Tfluiddownmidpointsstore = np.zeros((N, len(times)))
-        Tfluiddownmidpointsstore = np.zeros((N + 1, len(times)))
-        Pfluidupnodesstore = np.zeros((N + 1, len(times)))
-        Pfluiddownnodesstore = np.zeros((N + 1, len(times)))
-        # MIR Pfluidupmidpointsstore = np.zeros((N, len(times)))
-        # MIR Pfluiddownmidpointsstore = np.zeros((N, len(times)))
-        Pfluidupmidpointsstore = np.zeros((N + 1, len(times)))
-        Pfluiddownmidpointsstore = np.zeros((N + 1, len(times)))
-        Qfluidupnodesstore = np.zeros((N + 1, len(times)))
-        Qfluiddownnodesstore = np.zeros((N + 1, len(times)))
-        Phasefluidupnodesstore = np.zeros((N + 1, len(times)))
-        Phasefluiddownnodesstore = np.zeros((N + 1, len(times)))
-        Hfluidupnodesstore = np.zeros((N + 1, len(times)))
-        Hfluiddownnodesstore = np.zeros((N + 1, len(times)))
+        #Tfluidupnodesstore = np.zeros((N + 1, len(times)))
+        #Tfluiddownnodesstore = np.zeros((N + 1, len(times)))
+        Tfluidupnodesstore = np.zeros((N, len(times)))
+        Tfluiddownnodesstore = np.zeros((N, len(times)))
+        Tfluidupmidpointsstore = np.zeros((N, len(times)))
+        #Tfluidupmidpointsstore = np.zeros((N + 1, len(times)))
+        Tfluiddownmidpointsstore = np.zeros((N, len(times)))
+        #Tfluiddownmidpointsstore = np.zeros((N + 1, len(times)))
+        #Pfluidupnodesstore = np.zeros((N + 1, len(times)))
+        #Pfluiddownnodesstore = np.zeros((N + 1, len(times)))
+        Pfluidupnodesstore = np.zeros((N, len(times)))
+        Pfluiddownnodesstore = np.zeros((N, len(times)))
+        Pfluidupmidpointsstore = np.zeros((N, len(times)))
+        Pfluiddownmidpointsstore = np.zeros((N, len(times)))
+        #Pfluidupmidpointsstore = np.zeros((N + 1, len(times)))
+        #Pfluiddownmidpointsstore = np.zeros((N + 1, len(times)))
+        #Qfluidupnodesstore = np.zeros((N + 1, len(times)))
+        #Qfluiddownnodesstore = np.zeros((N + 1, len(times)))
+        #Phasefluidupnodesstore = np.zeros((N + 1, len(times)))
+        #Phasefluiddownnodesstore = np.zeros((N + 1, len(times)))
+        #Hfluidupnodesstore = np.zeros((N + 1, len(times)))
+        #Hfluiddownnodesstore = np.zeros((N + 1, len(times)))
+        Qfluidupnodesstore = np.zeros((N, len(times)))
+        Qfluiddownnodesstore = np.zeros((N, len(times)))
+        Phasefluidupnodesstore = np.zeros((N, len(times)))
+        Phasefluiddownnodesstore = np.zeros((N, len(times)))
+        Hfluidupnodesstore = np.zeros((N, len(times)))
+        Hfluiddownnodesstore = np.zeros((N, len(times)))
         Qinterexchangestore = np.zeros((N, len(times)))
         QinterexchangeUp = np.zeros(N)
+        velocityfluiddownmidpointsstore = np.zeros((N, len(times)))
+        heatcapacityfluidupmidpointsstore = np.zeros((N, len(times)))
 
         # Store initial values
         Tfluidupnodesstore[:, 0] = Tfluidupnodes
@@ -841,6 +862,8 @@ def Calculate_Coaxial(self, model):
         Phasefluidupnodesstore[:, 0] = Phasefluidupnodes
         Phasefluiddownnodesstore[:, 0] = Phasefluiddownnodes
         Qinterexchangestore[:, 0] = np.zeros(N)
+        velocityfluiddownmidpointsstore[:, 0] = velocityfluiddownmidpoints
+        heatcapacityfluidupmidpointsstore[:, 0] = heatcapacityfluidupmidpoints
 
         print('Pre-processing completed successfully. Starting simulation ...')
 
@@ -937,7 +960,8 @@ def Calculate_Coaxial(self, model):
             else:
                 maxindextoconsider = np.where(maxspacingtest > 1)[0][-1]
 
-            if mindexNPCP < maxindextoconsider + 1:
+            #if mindexNPCP < maxindextoconsider + 1:
+            if mindexNPCP < maxindextoconsider:
                 indicestocalculate = SortedIndices[:, mindexNPCP + 1:maxindextoconsider + 1]
                 NPCP[range(N), indicestocalculate] = Deltaz[indicestocalculate] / (
                         4 * np.pi * k_m_vector[indicestocalculate] * SMatrix[range(N), indicestocalculate]) * erfc(
@@ -1022,40 +1046,53 @@ def Calculate_Coaxial(self, model):
 
             Q[:, i] = BB
 
-            R[0:N] = Vvector * Tfluiddownmidpointsstore[:, i - 1] / Deltat
-            R[0:N] += Q[:, i]
-            R[N:2 * N] = Vvector * Tfluidupmidpointsstore[:, i - 1] / Deltat
-            R[N:2 * N] += Q[:, i]
-
-            R[2 * N:3 * N] = Pfluidupmidpointsstore[:, i - 1] - Pfluiddownmidpointsstore[:,
+            #R[0:N] = Vvector * Tfluiddownmidpointsstore[:, i - 1] / Deltat
+            #R[0:N] += Q[:, i]
+            #R[N:2 * N] = Vvector * Tfluidupmidpointsstore[:, i - 1] / Deltat
+            #R[N:2 * N] += Q[:, i]
+            R[0:N, 0] = Vvector * Tfluiddownmidpointsstore[:, i - 1] / Deltat
+            R[0:N, 0] += Q[:, i]
+            R[N:2 * N, 0] = Vvector * Tfluidupmidpointsstore[:, i - 1] / Deltat
+            R[N:2 * N, 0] += Q[:, i]
+
+            #R[2 * N:3 * N] = Pfluidupmidpointsstore[:, i - 1] - Pfluiddownmidpointsstore[:,
+            #                                                    i - 1] + velocityfluiddownmidpointsstore[:,
+            #                                                             i - 1] * Pfluiddownmidpointsstore[:,
+            #                                                                      i - 1] * Deltat
+            R[2 * N:3 * N, 0] = Pfluidupmidpointsstore[:, i - 1] - Pfluiddownmidpointsstore[:,
                                                                 i - 1] + velocityfluiddownmidpointsstore[:,
                                                                          i - 1] * Pfluiddownmidpointsstore[:,
                                                                                   i - 1] * Deltat
 
-            R[3 * N:4 * N] = heatcapacityfluidupmidpointsstore[:, i - 1] * Tfluidupmidpointsstore[:, i - 1]
+            #R[3 * N:4 * N] = heatcapacityfluidupmidpointsstore[:, i - 1] * Tfluidupmidpointsstore[:, i - 1]
+            R[3 * N:4 * N, 0] = heatcapacityfluidupmidpointsstore[:, i - 1] * Tfluidupmidpointsstore[:, i - 1]
 
             try:
                 solutions = np.linalg.solve(L, R)
             except np.linalg.LinAlgError:
                 print(f'Simulation terminated prematurely due to linear algebra error at time step {i}.')
                 break
 
-            BB = solutions[0:N]
+            #BB = solutions[0:N]
+            BB = solutions[0:N, 0]
             Tfluidupmidpointsstore[:, i] = BB
             Tfluidupmidpoints = BB
             Tfluidupmidpoints = Tfluidupmidpoints
 
-            BB = solutions[N:2 * N]
+            #BB = solutions[N:2 * N]
+            BB = solutions[N:2 * N, 0]
             Tfluiddownmidpointsstore[:, i] = BB
             Tfluiddownmidpoints = BB
             Tfluiddownmidpoints = Tfluiddownmidpoints
 
-            BB = solutions[2 * N:3 * N]
+            #BB = solutions[2 * N:3 * N]
+            BB = solutions[2 * N:3 * N, 0]
             Pfluidupmidpointsstore[:, i] = BB
             Pfluidupmidpoints = BB
             Pfluidupmidpoints = Pfluidupmidpoints
 
-            BB = solutions[3 * N:4 * N]
+            #BB = solutions[3 * N:4 * N]
+            BB = solutions[3 * N:4 * N, 0]
             Pfluiddownmidpointsstore[:, i] = BB
             Pfluiddownmidpoints = BB
             Pfluiddownmidpoints = Pfluiddownmidpoints

tests/test_geophires_x.py

@@ -913,3 +913,15 @@ def test_negative_electricity_production_raises_error(self):
             )
             client.get_geophires_result(params)
         self.assertIn('Electricity production calculated as negative', str(e.exception))
+
+    def test_sbt_coaxial_raises_error(self):
+        client = GeophiresXClient()
+        with self.assertRaises(RuntimeError) as e:
+            params = GeophiresInputParameters(
+                {
+                    'Reservoir Model': 8,
+                    'Well Geometry Configuration': 2,
+                }
+            )
+            client.get_geophires_result(params)
+        self.assertIn('SBT with coaxial configuration is not implemented', str(e.exception))