diff --git a/.github/workflows/formatting.yml b/.github/workflows/formatting.yml index 480c3bf..2f02d52 100644 --- a/.github/workflows/formatting.yml +++ b/.github/workflows/formatting.yml @@ -1,5 +1,5 @@ name: Check code formatting -on: [push, pull_request] +on: [push, pull_request, workflow_dispatch] jobs: check-formatting: runs-on: ubuntu-latest diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 6b17bca..6ec1043 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -1,5 +1,5 @@ name: tests -on: [push, pull_request] +on: [push, pull_request, workflow_dispatch] jobs: test-linux: runs-on: ubuntu-latest diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..12ecc8f --- /dev/null +++ b/.gitmodules @@ -0,0 +1,6 @@ +[submodule "moose"] + path = moose + url = https://github.com/idaholab/moose.git +[submodule "nemlapp"] + path = nemlapp + url = https://github.com/sagarbhatt0904/nemlapp.git diff --git a/install.sh b/install.sh new file mode 100755 index 0000000..8d07d92 --- /dev/null +++ b/install.sh @@ -0,0 +1,81 @@ +#!/bin/bash +set -e + +# Parse arguments +BUILD_MOOSE=false +if [[ "$1" == "--with-moose" ]]; then + BUILD_MOOSE=true + echo "Installing srlife with MOOSE..." + # Initialize MOOSE submodule + git submodule update --init moose + cd moose + git checkout master + cd ../ +else + echo "Installing srlife only..." +fi + +# Get repo root +SRLIFE_DIR="$(pwd)" +ENV_NAME="srlifeMoose" +MOOSE_JOBS=32 + +# Create conda environment +echo "Creating conda environment..." +# ensure conda shell functions are available in this script +eval "$(conda shell.bash hook)" + +# Create the environment only if it doesn't already exist +if conda env list | awk '!/^#/ && NF>0 {print $1}' | grep -Fxq "$ENV_NAME"; then + echo "Conda environment '$ENV_NAME' already exists; skipping creation." +else + conda create -n "$ENV_NAME" moose-dev=2025.09.18=mpich + #TODO:install seacas --- maybe add this to requirements.txt later + conda install moose-seacas +fi + + +# Activate environment +conda activate $ENV_NAME + +# Install Python dependencies +echo "Installing Python dependencies..." +pip3 install --user wheel +pip3 install -r requirements.txt + + + + +# Build MOOSE (if specified) +if [[ "$BUILD_MOOSE" == "true" ]]; then + echo "Building MOOSE..." + cd $SRLIFE_DIR/moose/test + make -j$MOOSE_JOBS + + # Set environment variables + echo "Installing Thermohydraylics module..." + cd $SRLIFE_DIR/moose/modules/thermal_hydraulics/ + make -j$MOOSE_JOBS + + echo "Installing nemlapp for solid mechanics..." + cd $SRLIFE_DIR/nemlapp + make -j$MOOSE_JOBS +fi + +echo "Installation complete!" +echo "" +if [[ "$BUILD_MOOSE" == "true" ]]; then + echo "To use add the following to your rc file:" + echo " conda activate $ENV_NAME" + echo " export PYTHONPATH=$SRLIFE_DIR/" + echo " export NEML_DIR=/neml" + echo " export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/neml/lib" + echo " export MOOSE_DIR=$SRLIFE_DIR/moose/" + echo " export MOOSE_THM=$SRLIFE_DIR/moose/modules/thermal_hydraulics/thermal_hydraulics-opt" + echo " export NEMLAPP=$SRLIFE_DIR/nemlapp/nemlapp-opt" + echo " export MOOSE_MPI=$(which mpirun)" + echo " export MOOSE_NPROCS=" +else + echo "To use:" + echo " conda activate $ENV_NAME" +fi \ No newline at end of file diff --git a/moose b/moose new file mode 160000 index 0000000..296c9e8 --- /dev/null +++ b/moose @@ -0,0 +1 @@ +Subproject commit 296c9e817d13b35624be2423775b309a34d9336c diff --git a/nemlapp b/nemlapp new file mode 160000 index 0000000..edcf5c8 --- /dev/null +++ b/nemlapp @@ -0,0 +1 @@ +Subproject commit edcf5c8959ac46a3e69167aeefb54f8d83d5cb6f diff --git a/requirements.txt b/requirements.txt index d599151..731d36b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,5 @@ numpy -scipy==1.9.2 +scipy>=1.9.2 h5py vtk pytest diff --git a/srlife/damage.py b/srlife/damage.py index 4b9053e..56043e5 100644 --- a/srlife/damage.py +++ b/srlife/damage.py @@ -1033,7 +1033,7 @@ def calculate_volume_flaw_flattened_eq_stress( else: sigma_e[sigma_e < 0] = 0 g = ( - np.trapz( + np.trapezoid( (sigma_e / (sigma_e_max + self.tolerance)) ** Navg[..., None, None], time, axis=0, @@ -1196,7 +1196,7 @@ def calculate_surface_flaw_flattened_eq_stress( else: sigma_e[sigma_e < 0] = 0 g = ( - np.trapz( + np.trapezoid( (sigma_e / (sigma_e_max + self.tolerance)) ** Navg[..., None, None, None], time, @@ -1474,7 +1474,7 @@ def calculate_volume_flaw_element_log_reliability( # Calculating ratio of cyclic stress to max cyclic stress (in one cycle) # for time-independent and time-dependent cases g = ( - np.trapz( + np.trapezoid( (pstress / (pstress_max + 1.0e-14)) ** Navg[..., None], time, axis=0, @@ -1556,7 +1556,7 @@ def calculate_surface_time_dep_pstress( # Calculating ratio of cyclic stress to max cyclic stress (in one cycle) # for time-independent and time-dependent cases g = ( - np.trapz( + np.trapezoid( (surf_pstress / (surf_pstress_max + 1.0e-14)) ** Navg[..., None, None], time, axis=0, @@ -1743,7 +1743,7 @@ def calculate_volume_flaw_avg_normal_stress( sigma_n_0 = sigma_n else: g = ( - np.trapz( + np.trapezoid( (sigma_n / (sigma_n_max + self.tolerance)) ** Navg[..., None, None], time, axis=0, @@ -1890,7 +1890,7 @@ def calculate_surface_flaw_time_dep_normal_stress( sigma_n_0 = sigma_n else: g = ( - np.trapz( + np.trapezoid( (sigma_n / (sigma_n_max + self.tolerance)) ** Navg[..., None, None, None], time, diff --git a/srlife/data/damage/SiC.xml b/srlife/data/damage/SiC.xml index 64e6f87..0e38c92 100644 --- a/srlife/data/damage/SiC.xml +++ b/srlife/data/damage/SiC.xml @@ -46,12 +46,12 @@ - 298.15 1073.15 1273.15 1473.15 1673.15 1773.15 - 0.0 0.0 0.0 0.0 0.0 0.0 + 298.15 811.15 1089.15 1366.15 1566.15 + 0.0 0.0 0.0 0.0 0.0 - 298.15 1073.15 1273.15 1473.15 1673.15 1773.15 - 0.0 0.0 0.0 0.0 0.0 0.0 + 298.15 811.15 1089.15 1366.15 1566.15 + 0.0 0.0 0.0 0.0 0.0 298.15 811.15 1089.15 1366.15 1566.15 @@ -88,4 +88,100 @@ 520 520 750 210 210 + + + 298.15 873.15 1073.15 1273.15 + 0.0 0.0 0.0 0.0 + + + 298.15 873.15 1073.15 1273.15 + 0.0 0.0 0.0 0.0 + + + 298.15 873.15 1073.15 1273.15 + 515.3 515.3 515.3 515.3 + + + 298.15 873.15 1073.15 1273.15 + 711.8 711.8 711.8 711.8 + + + 298.15 873.15 1073.15 1273.15 + 8.29 8.29 8.29 8.29 + + + 298.15 873.15 1073.15 1273.15 + 6.75 6.75 6.75 6.75 + + 0.82 + 0.16 + + 298.15 873.15 1073.15 1273.15 + 63.65 63.65 63.65 63.65 + + + + 298.15 873.15 1073.15 1273.15 + 63.65 63.65 63.65 63.65 + + + + 298.15 873.15 1073.15 1273.15 + 14.29 2.84 2.84 2.84 + + + + 298.15 873.15 1073.15 1273.15 + 14.29 2.84 2.84 2.84 + + + + + + 298.15 873.15 1073.15 1273.15 + 252.2 252.2 252.2 252.2 + + + 298.15 873.15 1073.15 1273.15 + 191.9 191.9 191.9 191.9 + + + 298.15 873.15 1073.15 1273.15 + 1690.2 1690.2 1690.2 1690.2 + + + 298.15 873.15 1073.15 1273.15 + 1291.6 1291.6 1291.6 1291.6 + + + 298.15 873.15 1073.15 1273.15 + 1.11 1.11 1.11 1.11 + + + 298.15 873.15 1073.15 1273.15 + 2.02 2.02 2.02 2.02 + + 0.82 + 0.16 + + 298.15 873.15 1073.15 1273.15 + 63.65 63.65 63.65 63.65 + + + + 298.15 873.15 1073.15 1273.15 + 63.65 63.65 63.65 63.65 + + + + 298.15 873.15 1073.15 1273.15 + 14.29 2.84 2.84 2.84 + + + + 298.15 873.15 1073.15 1273.15 + 14.29 2.84 2.84 2.84 + + + diff --git a/srlife/interface.py b/srlife/interface.py new file mode 100644 index 0000000..3ac0a95 --- /dev/null +++ b/srlife/interface.py @@ -0,0 +1,1061 @@ +""" + This module has all of the helper functions needed to run an srlife analysis. + It is meant to allow a simple script interface to define a problem and run the + analysis. + + One such example is defined in test/test_interface.py +""" + +import numpy as np +from scipy.interpolate import RegularGridInterpolator + +# srlife specific functions +from srlife import ( + receiver, + solverparams, +) + + +# UNIT CONVERSION FUNCTIONS +def convert_mm_to_m(mm_in): + """ + Unit conversion from mm to m + + Args: mm_in (double): qty to convert + """ + return mm_in / 1000 + + +def convert_m_to_mm(m_in): + """ + Unit conversion from m to mm + + Args: m_in (double): qty to convert + """ + return m_in * 1000 + + +def convert_kWm2_to_Wmm2(kWm2_in): + """ + Unit conversion from kW/m^2 to W/mm^2 + + Args: kWm2_in (double): qty to convert + """ + return kWm2_in * 1000e-6 + + +def convert_Pa_to_MPa(Pa_in): + """ + Unit conversion from Pa to MPa + + Args: Pa_in (double): qty to convert + """ + return Pa_in * 1e-6 + + +def convert_MPa_to_Pa(MPa_in): + """ + Unit conversion from MPa to Pa + + Args: MPa_in (double): qty to convert + """ + return MPa_in * 1e6 + + +def convert_C_to_K(C_in): + """ + Unit conversion from C to K + + Args: C_in (double): qty to convert + """ + return C_in + 273.15 + + +def read_month_day_hour_flux_file(flux_file_name, data_shape, z_offset_in): + """ + Read and process the flux data from a solarPILOT file + + Args: + flux_file_name (string): path and name of flux file from solarpilot + data_shape (list[int]): list with n_rows and n_cols of flux data + z_offset_in (double): z offset between srlife model and solarpilot model + this will shift z data that is read up or down + + Returns: + z_data (np.array(double)): vector with z positions from flux sampling + theta_data (np.array(double)): vector with theta positions from flux sampling + flux_data (np.array(double)): matrix with flux data at sampled positions + """ + # read flux data + # massage into format + # output + row, col = data_shape + z_data = np.zeros([row, 1]) + theta_data = np.zeros([1, col]) + flux_data = np.zeros([row, col]) + try: + theta_data = np.loadtxt( + flux_file_name, delimiter=",", skiprows=6, usecols=range(1, col + 1) + )[0] + z_data = np.loadtxt( + flux_file_name, delimiter=",", skiprows=7, usecols=range(0, 1) + ) + z_data += z_offset_in + flux_data = np.loadtxt( + flux_file_name, delimiter=",", skiprows=7, usecols=range(1, col + 1) + ) + # Data is on cylinder, so 180 will repeat + theta_data = np.append(theta_data, -180) + flux_data = np.hstack((flux_data, flux_data[:, 0][:, None])) + except FileNotFoundError: + # if we cannot find the file + # print a message and just return all zeros + print( + f"File could not be found!!! Returning all zeros\n filename: {flux_file_name}" + ) + # convert height to mm + z_data = convert_m_to_mm(z_data) + return z_data, theta_data, flux_data + + +def get_flux_interpolators_from_data_files( + times, start_time, end_time, month, day, flux_data_dir, flux_data_shape, z_offset +): + """ + Read flux files and create interpolator functions for each month, day, hour + + Args: + times (list(int)): list of index of hours used for analysis + Could probably remove this if (end_time - start_time) == len(times) + start_time (int): first hour to read flux data + end_time (int): last hour to read flux data + month (int): month to read flux data + day (int): day to read flux data + flux_data_dir (string): path to flux data files + flux_data_shape (list(int)): n_rows and n_cols of flux data in flux files + z_offset_in (double): z offset between srlife model and solarpilot model + this will shift z data that is read up or down + + Returns: + flux_interpolators_by_hour (list(RegularGridInterpolators)): list with the interpolating + functions for each hour in times. + """ + print("Making flux interpolators from data") + # Time zero is NOT analyzed, so we do not make a fluxFn for that + flux_interpolators_by_hour = [0] * (len(times) - 1) + for i_count, i_time in enumerate(range(start_time, end_time + 1)): + flux_filename = f"{flux_data_dir}/data_{month}_d{day}_hr{int(i_time)}.csv" + z_pts, theta_pts, flux_pts = read_month_day_hour_flux_file( + flux_filename, flux_data_shape, z_offset + ) + interpFun = RegularGridInterpolator( + (z_pts, theta_pts), + flux_pts, + method="linear", + bounds_error=False, + fill_value=0.0, + ) + flux_interpolators_by_hour[i_count] = interpFun + return flux_interpolators_by_hour + + +def apply_tube_flux_bcs( + tube, + tube_absorbance, + tube_rec_theta, + flux_interpolators_by_hour, +): + """ + Actually apply the heat flux to the tube across width and length + + Args: + tube (srlife.Tube): srlife object for this tube + tube_z (list[double]): list of z coordinates of tube nodes + tube_theta (list[double]): list of theta coordinates of tube nodes + tube_absorbance (double): absorbance of tubes + tube_rec_theta (double): receiver theta value for this tube + flux_interpolators_by_hour (list[RegularGridInterpolator]): + list of interpolating functions for heat flux on tubes + """ + num_steps = len(flux_interpolators_by_hour) + tube_theta = np.linspace(0, 360, tube.nt + 1)[:-1] + tube_z = np.linspace(0, tube.h, tube.nz) + # add extra time step for time zero + tube_flux = np.zeros([num_steps + 1, len(tube_theta), len(tube_z)]) + for i_hour in range(num_steps): + # for each hour of flux data + # sample tube at tube_rec_theta and z + for i_z, z in enumerate(tube_z): + # for each zPt + pt_z_theta = [z, tube_rec_theta] + rec_flux = convert_kWm2_to_Wmm2( + flux_interpolators_by_hour[i_hour](pt_z_theta) + ) + # apply across crown so only sunSide gets flux + for i_theta, theta in enumerate(tube_theta): + # for each thetaPt around tube + if 0 < theta < 180: + tube_flux[i_hour + 1, i_theta, i_z] = rec_flux[0] * np.cos( + (theta - 90) * np.pi / 180 + ) + else: + tube_flux[i_hour + 1, i_theta, i_z] = 0.0 + # make bc object for tube with heat flux data + # NOTE: this is done here and not in cerate_tube because if we initialize + # bc data to np.zeros, flux will always be zero unless we reload from a file + tube.set_bc(apply_heatflux_bc(tube, tube_flux * tube_absorbance, "outer"), "outer") + # while doing this, also apply convective BCs to inner face + tube.set_bc(apply_convective_bc(tube, "inner"), "inner") + + +def apply_convective_bc(tube, loc): + """ + Defines and applies a convective BC to a surface of a tube + specificed by loc + + Args: + tube (srlife.Tube): srlife object for this tube + loc (string): "inner" or "outer" to specify which surface BC is applied to + + """ + nz = tube.nz + if loc == "outer": + r = tube.r + else: + r = tube.r - tube.t + return receiver.ConvectiveBC( + r, tube.h, nz, tube.times, np.zeros((len(tube.times), nz)) + ) + + +def apply_heatflux_bc(tube, data, loc): + """ + Defines and applied a heat flux BC to a surface of a tube object + specified by loc + + Args: + tube (srlife.Tube): srlife object for this tube + data (list[double]): 3d list of heat flux bcs [hour, theta, z] + loc (string): "inner" or "outer" to specify which surface BC is applied to + + """ + if loc == "outer": + r = tube.r + else: + r = tube.r - tube.t + nt = tube.nt + nz = tube.nz + return receiver.HeatFluxBC(r, tube.h, nt, nz, tube.times, data) + + +def calc_and_write_tube_flux_bcs( + rec, + ass_tube_per_panel, + num_panels, + flux_interpolators_by_hour, + tube_absorbance, +): + """ + Applies flux boundary conditions to analyzed tubes. Flux value is interpolated based on + centerline tube location. It is applied along the tube face as a cosine distribution + with max at sun-facing centerline to zero at 90 deg in either direction. + + Args: + rec (Receiver): the receiver object that the tubes belong to + ass_tube_per_panel (int): the assumed number of tubes per panel, or analyzed number of tubes + num_panels (int): number of panels in the receiver #BPMToDo: remove this and get from rec + flux_interpolators_by_hour (list(RegularGridInterpolators)): list of flux interpolators + on receiver by hour of analysis + tube_z (list(double)): z nodes points along tube + tube_theta (list(double)): theta nodes points along tube + tube_absorbance (double): absorbance of tubes + write_tube_flux_to_csv (bool): flag to dump flux data per tube to a csv file to debug + + Returns: None + """ + print("Assigning flux BCs to tubes") + for panel_key, panel in rec.panels.items(): + # in each panel of rec + for tube_key, tube in panel.tubes.items(): + print(f"Panel: {panel_key} Tube: {tube_key}") + # in each tube of panel + tube_rec_theta = 180 - (360 * int(panel_key) / num_panels) + if ass_tube_per_panel == 1: + # if only one tube, place at center of panel + tube_rec_theta += 360 / (2 * num_panels) + else: + # if more than one tube, evenly space them including edges + tube_rec_theta -= ( + 360 / num_panels * int(tube_key) / (ass_tube_per_panel - 1) + ) + if tube_rec_theta == -180: + # wrap data around tube + tube_rec_theta = 180 + apply_tube_flux_bcs( + tube, + tube_absorbance, + tube_rec_theta, + flux_interpolators_by_hour, + ) + + +def set_rec_flow_paths(rec, panel_flow_path, mass_flow_per_path, T_in_per_path): + """ + Set flowpaths within a receiver object + + Args: + rec (Receiver): the receiver object to set flowpaths for + panel_flow_path (list[string]): panel names indicating the flowpath. Can have a list + of lists if more that one flowpath + mass_flow_per_path (list[double]): list with entries of initial mass flow rate for + each path + T_in_per_path (list[double]): list with entries of inlet temp for each path + """ + print("Setting receiver flowpaths") + T_in_per_path = convert_C_to_K(T_in_per_path) + # add flowpaths to receiver + for i_path, flow_path_panels in enumerate(panel_flow_path): + # for each path in given list + print(f"Path {i_path}") + flow_path = [0] * len(flow_path_panels) + for i_panel, panel in enumerate(flow_path_panels): + # for each panel in path + print(f"Panel: {panel}") + flow_path[i_panel] = panel + panel = rec.panels[panel] + for tube in panel.tubes.values(): + # for each tube in panel + times = tube.times + # initialize mass flow for path + mass_flow = mass_flow_per_path[i_path] + # convert mass flow from kg/s to kg/hr + mass_flow *= 3600 + # set inlet temp + T_in_flow_path = T_in_per_path[i_path] * np.ones_like(times) + rec.add_flowpath(flow_path, times, mass_flow, T_in_flow_path) + + +def check_outlet_temps_and_step_mass_flow( + rec, path, path_key, T_out_target, pct_err_outlet_temp, breaker +): + """ + Check outlet temp of all tubes from last panel in flowpath. + Average their value and use it to step mass flow rate. + + Args: + rec (srlife.Receiver): Receiver object who owns the flowpath we are + optimizing + path (OrderedDict): dictionary of flowpath information + path_key (string): name of flowpath + T_out_target (double): target outlet temperature + pct_err_outler_temp (double): acceptable error % for outlet temp (for convergence) + breaker (list[bool]): check whether all flowpaths have converged + """ + # All tubes in last panel go to manifold, so estimate + # temp as avg of all tubes output + panels = path.get("panels") + T_in_path = path["inlet_temp"][1:] + last_panel = rec.panels[panels[-1]] + num_tubes = last_panel.ntubes + num_time_steps = next(iter(last_panel.tubes.values())).ntime + T_out_path = np.zeros([num_tubes, num_time_steps - 1]) + for i_tube, tube in enumerate(last_panel.tubes.values()): + # for tubes in the last panel of the flow path + # results for tube are (time, nz) shape + # want all time except 0 and the last zPos + T_out_path[i_tube] = tube.axial_results.get("fluid_temperature")[1:, -1] + T_out_path = T_out_path.mean(axis=0) + print(f"T_out = {T_out_path - 273.15}") + if ( + any(T_out_path > T_out_target * (1 + pct_err_outlet_temp / 100)) is False + and any(T_out_path < T_out_target * (1 - pct_err_outlet_temp / 100)) is False + ): + # we have converged this flowpath + breaker[int(path_key)] = True + else: + mass_flow = path.get("mass_flow") + mass_flow[1:] = ( + mass_flow[1:] * (T_out_path - T_in_path) / (T_out_target - T_in_path) + ) + mass_flow[0] = mass_flow[1] + rec.flowpaths[path_key]["mass_flow"] = mass_flow + print(f"Mass flow: {mass_flow/3600}") + + +def optimize_mass_flow_rate_per_path( + rec, rec_filename, T_out_target, pct_err_outlet_temp, solver, save_heat_to_vtu +): + """ + Iteratively calculate the ideal mass flow rate for each path such that the + outlet temperature matches a target temperature. + + Args: + rec (Receiver): srlife Receiver object we are optimizing flow for + rec_filename (string): file name for the receiver object to optimize + T_out_target (double): target output temperature + pct_err_outlet_temp (double): acceptable percentage of error for the outlet temp + solver (managers.SolutionManager): srlife solver object + save_heat_to_vtu (bool): whether or not to save vtu files of results + """ + print("Optimizing mass flow rate to match outlet temp") + T_out_target = convert_C_to_K(T_out_target) + # make sure receiver is up to date + solver.receiver = rec + # Mass flow optimization + # limit of error on outlet temp + max_mass_flow_iter = 5 + for i_opt in range(max_mass_flow_iter): + print(f"Iteration = {i_opt}") + breaker = [False] * len(rec.flowpaths) + if i_opt > 0: + for path_key, path in rec.flowpaths.items(): + print(f"FlowPath: {path_key}") + # check last panel tube outlet temps + check_outlet_temps_and_step_mass_flow( + rec, path, path_key, T_out_target, pct_err_outlet_temp, breaker + ) + for panel_key, panel in rec.panels.items(): + for tube_key, tube in panel.tubes.items(): + tube.write_vtk(f"ht-tube-{panel_key}-{tube_key}") + + if all(breaker): + print(f"Converged!!! solved in {i_opt} iterations") + if save_heat_to_vtu: + for panel_key, panel in rec.panels.items(): + for tube_key, tube in panel.tubes.items(): + tube.write_vtk(f"ht-tube-{panel_key}-{tube_key}") + break + # solve heat transfer problem + solver.solve_heat_transfer() + filename = rec_filename + "_ht_iter_" + str(i_opt) + ".hdf5" + rec.save(filename) + + +def calc_fluid_velocity(tube_mass_flow, fluid_rho, tube_Dh): + """ + Calculate the velocity of the fluid in a tube based on mass flow rate + + Args: + tube_mass_flow (double): mass flow rate in tube + fluid_rho (double): fluid density + tube_Dh (double): hydraulic diameter of tube + """ + return tube_mass_flow / fluid_rho * (4.0 / (np.pi * tube_Dh**2)) + + +def calc_reynolds_number(fluid_rho, fluid_vel, fluid_mu, tube_Dh): + """ + Calculate the Reynolds number for tube flow + + Args: + fluid_rho (double): fluid density + fluid_vel (double): fluid velocity + fluid_mu (double): fluid dynamic viscosity + tube_Dh (double): tube hydraulic diameter + + """ + return fluid_rho * fluid_vel * tube_Dh / fluid_mu + + +def calc_fd(Re, tube_eta, tube_Dh): + """ + Calculate Darcy friction factor from: + Zigrange and Sylvester 1985, A review of explicit friction factor equations. + J of Energy Resources Technology) + Equation 13 + + Args: + Re (double): Reynolds number + tube_eta (double): tube roughness + tube_Dh (double): tube hydraulic diameter + Returns: + fd (double): darcy friciton factor + """ + if Re < 4000: + # laminar flow + fd = 64 / Re + else: + fd = ( + 1 + / ( + -2 + * np.log10( + (tube_eta / tube_Dh) / 3.7 + - 5.02 / Re * np.log10((tube_eta / tube_Dh) / 3.7 + 13 / Re) + ) + ) + ** 2 + ) + return fd + + +def calc_friction_p_loss(fluid_rho, fluid_mu, tube_mass_flow, tube_dict, delta_l): + """ + Calculate pressure loss from friction in tubes + + Args: + fluid_rho (double): fluid density in tubes + fluid_mu (double): fluid dynamic visco in tubes + tube_mass_flow (double): mass flow rate in tube + tube_dict (dict): dictionary of tube specs + delta_l (double): length of tube + + Returns: + friction_p_loss (double): pressure loss from friction + """ + # tube roughness + tube_eta = tube_dict["eta"] + tube_Dh = convert_mm_to_m(tube_dict["od"] - 2 * tube_dict["t"]) + fluid_vel = calc_fluid_velocity(tube_mass_flow, fluid_rho, tube_Dh) + Re = calc_reynolds_number(fluid_rho, fluid_vel, fluid_mu, tube_Dh) + fd = calc_fd(Re, tube_eta, tube_Dh) + return fd * (0.5 * fluid_rho) * (fluid_vel**2) / tube_Dh * delta_l + + +def calc_manifold_bend_loss( + max_mass_flow, manifold_dict, manifold_rho, manifold_mu, num_bends +): + """ + Calculate pressure loss from flow through bends in manifold tubes + + Args: + max_mass_flow (double): mass flow in manifold + manifold_dict (dict): dictionary with manifold tube specs + manifold_rho (double): fluid density in manifold + manifold_mu (double): dynamic viscosity in manifold + num_bends (int): total number of bends in manifold tubes + """ + manifold_eta = manifold_dict["eta"] + manifold_pipe_t = convert_mm_to_m(manifold_dict["t"]) + manifold_pipe_id = convert_mm_to_m(manifold_dict["od"] - 2 * manifold_pipe_t) + manifold_pipe_bend_radius = convert_mm_to_m(manifold_dict["bend_radius"]) + manifold_eta = manifold_dict["eta"] + manifold_Dh = manifold_pipe_id + if manifold_eta / manifold_Dh > 0.0001: + K0 = 0.42 + elif manifold_eta / manifold_Dh == 0: + K0 = 0.21 + else: + K0 = 0.21 * (1 + 1000 * manifold_eta / manifold_Dh) + + manifold_vel = calc_fluid_velocity(max_mass_flow, manifold_rho, manifold_Dh) + Re = calc_reynolds_number(manifold_rho, manifold_vel, manifold_mu, manifold_Dh) + manifold_fd = calc_fd(Re, manifold_eta, manifold_Dh) + K = ( + K0 ** ((manifold_pipe_id / manifold_pipe_bend_radius) ** 0.5) + + 0.5 * np.pi * (manifold_pipe_bend_radius / manifold_pipe_id) * manifold_fd + ) + return num_bends * K * manifold_rho * manifold_vel**2 + + +def calc_path_total_p_loss( + rec, + path, + path_key, + delta_T, + tube_dict, + manifold_dict, + fluid, + rec_circumference, + num_panels, +): + """ + Calculate the total pressure loss in a flowpath from flow friction, + head loss, and flow through pipe bends. + + Args: + rec (Receiver): receiver object to analyze + path (dict): flowpath to get pressure loss through + path_key (string): name of flowpath + delta_T (double): temp step through flowpath + tube_dict (dict): dictionary of usefule tube info + manifold_dict (dict): dict of useful info for manifold tubes + fluid (thermalfluid.ThermalFluidMaterial): working fluid in receiver + rec_circumference (double): circumference of receiver + num_panels (int): number of panels in receiver + + Returns: + total_p_loss (double): total pressure loss in flowpath + """ + g = 9.81 # m/s^2 + # NOTE: this function works mostly with meters + act_tube_per_panel = tube_dict["ass_tube_per_panel"] * tube_dict["tube_mult"] + num_bends = num_panels * manifold_dict["bends_per_panel"] + + inlet_temp = path["inlet_temp"][0] + # get outlet avg outlet temp from each tube in last panel + # then average down again to single value + outlet_temp = np.mean( + [ + tube.axial_results.get("fluid_temperature")[:, -1][1:-1].mean() + for tube_key, tube in rec.panels[path["panels"][-1]].tubes.items() + ] + ) + max_mass_flow = path["mass_flow"].max() / 3600 # kg/s + tube_mass_flow = max_mass_flow / act_tube_per_panel + flow_path_length = convert_mm_to_m(tube_dict["h"]) * len(path["panels"]) + temp_steps = np.arange(inlet_temp, outlet_temp, delta_T) + rho_with_T = fluid.rho(temp_steps) * 1e9 # kg/m^3 + mu_with_T = fluid.mu(temp_steps) * 1000 / 3600 # Pa-s + + # Friction Loss + # pos spacing of temp change on tubes + delta_l = flow_path_length / np.size(temp_steps) + friction_p_loss = sum( + calc_friction_p_loss(rho, mu, tube_mass_flow, tube_dict, delta_l) + for rho, mu in zip(rho_with_T, mu_with_T) + ) + # Head Loss + avg_temp = 0.5 * (inlet_temp + outlet_temp) + head_p_loss = fluid.rho(avg_temp) * g * convert_mm_to_m(tube_dict["h"]) + # Manifold Loss + manifold_friction_p_loss = calc_friction_p_loss( + fluid.rho(avg_temp) * 1e9, + fluid.mu(avg_temp) * 1000 / 3600, + max_mass_flow, + manifold_dict, + (len(path["panels"]) - 1) * (rec_circumference / num_panels), + ) + manifold_bend_loss = calc_manifold_bend_loss( + max_mass_flow, + manifold_dict, + fluid.rho(avg_temp) * 1e9, + fluid.mu(avg_temp) * 1000 / 3600, + num_bends, + ) + return friction_p_loss + head_p_loss + manifold_bend_loss + manifold_friction_p_loss + + +def calc_p_loss_from_flows_temps(rec, tube_dict, manifold_dict, fluid, outlet_p): + """ + Calculate the pressure loss across the receiver based on the results + of a flow simulation. This captures friction loss, head loss, and bend losses + in the tubes and manifold pipes of the receiver + + Args: + rec (Receiver): the receiver object to calculate pressure loss + tube_dict (dict): a dictionary containing info about rec tubes + manifold_dict (dict): a dictionary containing info about manifold tubes + fluid (thermalfluid.ThermalFluidMaterial): thermal material of working fluid + outlet_p (double): pressure at flowpath outlets + """ + print("Calculating pressure loss from flow at temps") + # receiver details + num_panels = len(rec.panels) + # back out rec_diam + rec_diam = ( + tube_dict["ass_tube_per_panel"] + * tube_dict["tube_mult"] + * num_panels + * (tube_dict["od"] + tube_dict["spacing"]) + / np.pi + ) + rec_circumference = np.pi * convert_mm_to_m(rec_diam) + print(f"Rec diam = {convert_mm_to_m(rec_diam)}") + # This is number of tempsteps along tube + delta_T = 5 + flow_path_p_loss = np.zeros(len(rec.flowpaths)) + for path_key, path in rec.flowpaths.items(): + # for each flow path in rec + print(f"Flow path: {path_key}") + total_p_loss = calc_path_total_p_loss( + rec, + path, + path_key, + delta_T, + tube_dict, + manifold_dict, + fluid, + rec_circumference, + num_panels, + ) + flow_path_p_loss[int(path_key)] = total_p_loss + print(f"Total pressure loss= {(total_p_loss/1e6)} MPa") + print(f"Inlet pressure= {(outlet_p + total_p_loss/1e6)} MPa") + return flow_path_p_loss + + +def update_tube_pressure_bcs(rec, inlet_p_per_path, outlet_p): + """ + Set tube pressure BCs based on pressure loss calc after thm solve + + Args: + rec (Receiver): receiver object to set the pressure for + inlet_p_per_path (np.array): array of inlet pressures for each flowpath + outlet_p (double): outlet pressure, same for all paths + """ + # set tube pressures and temperatures + for path_key, path in rec.flowpaths.items(): + # for each flowpath in the model + tube_pressures = np.linspace( + inlet_p_per_path[int(path_key)], outlet_p, len(path["panels"]) + 1 + )[:-1] + for i_panel, panel_key in enumerate(path["panels"]): + # for each panel in flowpath + for tube in rec.panels[panel_key].tubes.values(): + # for each tube in panel + # get tube times + times = tube.times + pressure = np.ones_like(times) + pressure[0] = 0.0 + tube_pressure = tube_pressures[i_panel] * pressure + tube_pressure_bc = receiver.PressureBC(times, tube_pressure) + tube.set_pressure_bc(tube_pressure_bc) + + +def cycle_tube_pressure_bcs(tubes_dict, num_cycles, cyclic_times): + """ + Set pressure bcs on internal face of tube from flow + + Args: + tubes_dict (dict): rec.panels[panel].tubes dictionary tubes from a panel + tube_pressure (list[double]): pressure values for tubes in a given panel with height + pressure (list[double]): properly sized array of pressure data to define bc + times (list[int]): list of analysis times + + """ + for tube in tubes_dict.values(): + # NOTE: as soon as I cycle results + # ghost temp, fluid temp, fluid velocity results are no longer good + # so lets delete themo + try: + tube.quadrature_results.pop("ghost_temperature") + tube.axial_results.pop("fluid_temperature") + tube.axial_results.pop("fluid_velocity") + except KeyError: + # We do not need to do anything here, results are already gone + pass + # for each tube in panel + press_bc = tube.pressure_bc + pressure = press_bc.data + tube_pressure = np.tile(pressure[1:], num_cycles) + tube_pressure = np.append(0, tube_pressure) + tube_pressure_bc = receiver.PressureBC(cyclic_times, tube_pressure) + tube.set_pressure_bc(tube_pressure_bc) + + +def set_and_downsample_tube_temp_bcs( + tubes_dict, + set_init_T_to_inlet_T, + inlet_T, + cyclic_times, + num_cycles, + analysis_type, + loc, +): + """ + Handle specific details for temperature BCs in a tube. + Set the initial temperature + Find and set the appropriate temp for 2d analysis + + Args: + tubes_dict (dict): rec.panels[panel].tubes dictionary tubes from a panel + set_init_T_to_inlet_T (bool): whether to initialize tube temps as the + inlet temperature + inlet_T (double): inlet temperature for tubes + cyclic_times (list[int]): list of times across all cycles of analysis + num_cycles (int): number of cycles in analysis + analysis_type (string): "2d" or "3d" to handle turn on or off downsampling + of temperature results for simpler analysis + loc (string): "max_T" or "max_avg_T" to set where temp is sampled in 2d analysis + """ + for tube in tubes_dict.values(): + # for each tube in panel + T = tube.results["temperature"] + if set_init_T_to_inlet_T: + T[0] = inlet_T + tube.T0 = inlet_T + _, _, _, r4 = np.shape(T) + # this will tile stack all temp results num_cycles times + # It does not account for cycle heuristic + T_0 = np.array([T[0]]) + T = np.tile(T[1:], (num_cycles, 1, 1, 1)) + T = np.append(T_0, T, axis=0) + tube.results["temperature"] = T + tube.set_times(cyclic_times) + T_3d = tube.results["temperature"] + if analysis_type == "2d": + if loc == "max_T": + T_max = 0 + T_max_h_index = -1 + for i_temp in range(1, r4 - 1): + if np.max(T_3d[:, :, :, i_temp]) > T_max: + T_max = np.max(T_3d[:, :, :, i_temp]) + T_max_h_index = i_temp + tube.make_2D(tube.h / (r4 - 1) * T_max_h_index) + elif loc == "max_avg_T": + T_3d_avg = np.average(T_3d, axis=0) + T_max = 0 + T_max_h_index = -1 + for i_temp in range(1, r4 - 1): + if np.max(T_3d_avg[:, :, i_temp]) > T_max: + T_max = np.max(T_3d_avg[:, :, i_temp]) + T_max_h_index = i_temp + tube.make_2D(tube.h / (r4 - 1) * T_max_h_index) + tube.results["temperature"] = T_3d[:, :, :, T_max_h_index] + + +def process_single_panel_analysis(rec, struct_output_dict, single_panel_analysis_id): + """ + Remove all panels that will not be analyzed + + Args: + rec (Receiver): object we are analyzing + struct_output_dict (dict): dictionary with output information + save_struct_to_vtu (bool): whether or not to write files to vtu + st_fname (string): filename for structure + tube_fname (string): filename for tube + single_panel_analysis_is (string): name of panel to analyze + """ + # file name mods + st_fname = ( + f"sing_panel_{single_panel_analysis_id}_{struct_output_dict['st_filename']}" + ) + struct_output_dict["st_filename"] = st_fname + tube_fname = ( + f"sing_panel_{single_panel_analysis_id}_{struct_output_dict['tube_filename']}" + ) + struct_output_dict["tube_filename"] = tube_fname + # update model to use just one panel + single_panel_model = rec.panels[single_panel_analysis_id] + rec.panels.clear() + rec.flowpaths.clear() + rec.add_panel(single_panel_model) + + +def save_structural_results(rec_struct, st_fname, tube_fname): + """ + Save results from analysis to files + + Args: + rec_struct (Receiver): receiver object we analyzed + st_fname (string): hdf5 filename for structure + tube_fname (string): base filename for tube vtus + """ + # post process + rec_struct.save(st_fname + ".hdf5") + for panel_key, panel in rec_struct.panels.items(): + for tube_key, tube in panel.tubes.items(): + tube.write_vtk(f"{tube_fname}_{panel_key}_{tube_key}") + + +def run_struct_analysis( + rec_filename, + set_init_T_to_inlet_T, + loc, + analysis_type, + is_single_panel_analysis, + single_panel_analysis_id, + num_cycles, + solver, + struct_output_dict, +): + """ + Runs the structural analysis on the receiver with a given set of thermal + hydraulics results. + + Args: + rec_filename (string): filename of the receiver to analyze + set_init_T_to_inlet_T (bool): whether to initialize tube temps as the + inlet temperature + loc (string): "max_T" or "max_avg_T" used to select location to downsample + analysis domain if lower dimensional analysis is done + analysis_type (string): "2d" or "3d" used to decide whether full analysis + or downsampled analysis is used + is_single_panel_analysis (bool): whether a single panel is analyzed or the + full receiver + single_panel_analysis_id (string): name of single panel to analyze if + is_single_panel_analysis == True + num_cycles (double): number of times to repeat load cycles for the analysis + solver (managers.SolutionManager): system solver for the receiver + struct_output_dict (dict): dictionary with info about output for st files + save_to_vtu (bool): save files to vtu output + st_filename (string): hdf5 filename for structural output + tube_filename (string): filename for vtk output of tube analysis + + ToDos: + rename loc to better name, make it and analysis_type bools? + """ + # passing filename here because we may change receiver significantly + # i.e. removing many panels + rec_struct = receiver.Receiver.load(rec_filename + ".hdf5") + rec_struct.days *= num_cycles + inlet_T = rec_struct.flowpaths["0"]["inlet_temp"][0] + + # set tube pressures and temperatures + for path in rec_struct.flowpaths.values(): + # for each flowpath in the model + for panel_key in path["panels"]: + # for each panel in flowpath + # times we actually analyze receiver for thm + first_tube = next(iter(rec_struct.panels[panel_key].tubes.values())) + analysis_times = first_tube.times[1:] + # create cyclic time for life + cyclic_times = np.tile(analysis_times, num_cycles) + for i_cyc in range(num_cycles): + cyclic_times[ + i_cyc * len(analysis_times) : (i_cyc + 1) * len(analysis_times) + ] += (i_cyc * rec_struct.period) + # initial steps for pressure and time + cyclic_times = np.append([0], cyclic_times) + pressure = np.ones_like(cyclic_times) + pressure[0] = 0.0 + + cycle_tube_pressure_bcs( + rec_struct.panels[panel_key].tubes, + num_cycles, + cyclic_times, + ) + set_and_downsample_tube_temp_bcs( + rec_struct.panels[panel_key].tubes, + set_init_T_to_inlet_T, + inlet_T, + cyclic_times, + num_cycles, + analysis_type, + loc, + ) + # remove unnecessary panels if needed + if is_single_panel_analysis: + process_single_panel_analysis( + rec_struct, struct_output_dict, single_panel_analysis_id + ) + solver.receiver = rec_struct + # run structural problem + solver.solve_structural() + if struct_output_dict["save_to_vtu"]: + save_structural_results( + rec_struct, + struct_output_dict["st_filename"], + struct_output_dict["tube_filename"], + ) + rec_struct.save(rec_filename + "_struct.hdf5") + + +def create_receiver(tube_dict, num_days, times, period, panel_k, num_panels, results): + """ + Create an srlife.Receiver object based on user defined inputs + + Args: + tube_dict (dict): Python dictionary that holds useful tube definitions + "od", "t", "h", "nr", "nt", "nz", "spacing", "T0", "tube_k", + "eta", "tube_mult", "ass_tube_per_panel" + num_days (int): number of days to analyze + times (list[int]): list of index of hours to analyze + period (int): number of hours to analyze in a day + panel_k (double): panel stiffness used in SystemSolver object + num_panels (int): number of panels + results (list): set of results for tube objects + """ + # array of 1's with leading and trailing zeros multiplied by input pressure + rec = receiver.Receiver(period, num_days, panel_k) + for _ in range(num_panels): + rec.add_panel(create_panel(tube_dict, times, results)) + return rec + + +def create_panel(tube_dict, times, results): + """ + Create an srlife.Panel object to add to a receiver + + Args: + tube_dict (dictionary): Python dictionary that holds useful tube definitions + times (list[int]): list of index of hours to analyze + results (list): list of results for tube objects + + Returns: + panel (srlife.Panel): panel object + """ + panel = receiver.Panel(tube_dict["tube_k"]) + for _ in range(tube_dict["ass_tube_per_panel"]): + panel.add_tube(create_tube(tube_dict, times, results)) + return panel + + +def create_tube(tube_dict, times, results): + """ + Create an srlife.Tube object to add to panel + + Args: + tube_dict (dictionary): Python dictionary that holds useful tube definitions + times (list[int]): list of index of hours to analyze + results (list): list of results for tube object + + Returns: + tube (srlife.Tube): tube object + """ + tube = receiver.Tube( + 0.5 * tube_dict["od"], + tube_dict["t"], + tube_dict["h"], + tube_dict["nr"], + tube_dict["nt"], + tube_dict["nz"], + tube_dict["T0"], + ) + tube.set_times(times) + tube.multiplier_val = tube_dict["tube_mult"] + tube.T0 = tube_dict["T0"] + # initialize pressure data + pressure_data = np.ones(len(times)) * 1.0 + pressure_data[0] = 0.0 + pressure_data[-1] = 0.0 + pressure = receiver.PressureBC(times, pressure_data) + tube.set_pressure_bc(pressure) + for res in results: + tube.add_results( + res, + np.zeros((len(times), tube_dict["nr"], tube_dict["nt"], tube_dict["nz"])), + ) + return tube + + +def sample_parameters(num_threads, verbose, rtol, atol): + """ + Create a set of srlife.SolverParams to be used by srlife solvers + + Args: + num_threads (int): number of threads for analysis + verbose (bool): whether solvers should print verbose info + rtol (double): relative tol for solver convergence + atol (double): abs tol for solver convergence + Returns: None + """ + params = solverparams.ParameterSet() + + params["nthreads"] = num_threads + params["progress_bars"] = True + # If true store results on disk (slower, but less memory) + params["page_results"] = False + + params["thermal"]["miter"] = 200 + params["thermal"]["verbose"] = verbose + params["thermal"]["steady"] = True + params["thermal"]["substep"] = 2 # 10 + + params["thermal"]["solid"]["rtol"] = rtol + params["thermal"]["solid"]["atol"] = atol + params["thermal"]["solid"]["miter"] = 200 + params["thermal"]["solid"]["verbose"] = verbose + params["thermal"]["solid"]["substep"] = 2 # 10 + + params["thermal"]["fluid"]["rtol"] = rtol + params["thermal"]["fluid"]["atol"] = atol + params["thermal"]["fluid"]["miter"] = 200 + params["thermal"]["fluid"]["verbose"] = verbose + params["thermal"]["fluid"]["substep"] = 2 # 100 + + params["structural"]["rtol"] = rtol + params["structural"]["atol"] = atol + params["structural"]["miter"] = 50 + params["structural"]["verbose"] = verbose + + params["system"]["rtol"] = rtol + params["system"]["atol"] = atol + params["system"]["miter"] = 10 + params["system"]["verbose"] = verbose + + # If true store results on disk (slower, but less memory) + params["page_results"] = False + + return params diff --git a/srlife/moose_interface.py b/srlife/moose_interface.py new file mode 100644 index 0000000..8ecc4a9 --- /dev/null +++ b/srlife/moose_interface.py @@ -0,0 +1,580 @@ +""" +This module provides functions to interface between srlife and MOOSE. +It has functions to read MOOSE THM results, create MOOSE structural input files, +run the MOOSE structural model, and compute reliability from MOOSE structural outputs +using srlife's damage models. + +""" + +import os +import sys +import re +import subprocess +from pathlib import Path +import numpy as np +import pyhit # pylint: disable=import-error,wrong-import-position +from srlife.receiver import make_moose_hit_vector +from srlife.interface import convert_m_to_mm + +conda_env_dir = os.environ.get("CONDA_PREFIX") +ACCESS = os.getenv("ACCESS", f"{conda_env_dir}/seacas") +sys.path.append(os.path.join(ACCESS, "lib")) +import exodus as exo # pylint: disable=import-error,wrong-import-position + +SQRT2 = np.sqrt(2.0) + + +def find_pin_coords(mesh_file, tol=1e-6): + """ + Two diametrically-opposite outer-radius nodes on the bottom face. + This helps pin the bottom face without restricting the axial expansion + + Args: + mesh_file: path to the tube mesh exodus file + tol: tolerance for coordinate comparisons + + Returns: + (xa, ya, za), (xb, yb, zb): coordinates of the two pin nodes. + """ + model = exo.exodus(str(mesh_file), array_type="numpy") + x, y, z = model.get_coords() + model.close() + bot_mask = np.abs(z - z.min()) < tol + x_bot, y_bot = x[bot_mask], y[bot_mask] + cx, cy = x_bot.mean(), y_bot.mean() + dist = np.hypot(x_bot - cx, y_bot - cy) + outer_mask = np.abs(dist - dist.max()) < tol * 10 + x_outer, y_outer = x_bot[outer_mask], y_bot[outer_mask] + xa, ya = x_outer[0], y_outer[0] + idx_b = np.argmax(np.hypot(x_outer - xa, y_outer - ya)) + xb, yb = x_outer[idx_b], y_outer[idx_b] + z_bottom = float(z.min()) + return (float(xa), float(ya), z_bottom), (float(xb), float(yb), z_bottom) + + +PANEL_BLOCK_RE = re.compile(r"panel_(\d+)/") +FCH_TUBE_BLOCK_RE = re.compile(r"panel_(\d+)/fch_tube_(\d+)$") + + +def panels_in_thm_exodus(exo_path): + """ + Reads the THM exodus file and returns a sorted list of panel numbers present in the file. + + Args: + exo_path: path to the THM exodus file + Returns: + A sorted list of panel numbers (integers) found in the exodus file.""" + model = exo.exodus(str(exo_path), array_type="numpy") + names = [model.get_elem_blk_name(b) for b in model.get_elem_blk_ids()] + model.close() + panels = set() + for name in names: + m = PANEL_BLOCK_RE.match(name) + if m: + panels.add(int(m.group(1))) + return sorted(panels) + + +def read_time_axis(thm_exodus: Path): + """ + Reads the time axis from the THM exodus file. + + Args: + thm_exodus: Path to the THM exodus file + + Returns: + A list of time values. + """ + model = exo.exodus(str(thm_exodus), array_type="numpy") + times = list(model.get_times()) + model.close() + return times + +def discover_tubes(panel: int, out_dir: Path): + """ + Find tube IDs by globbing panel_{panel}_tube_*_in.e in out_dir + Args: + panel: int, panel number + out_dir: Path, directory to search for tube meshes + Returns: + A sorted list of tube IDs (integers) found in the directory. + """ + pattern = f"panel_{panel}_tube_*_in.e" + files = sorted(out_dir.glob(pattern)) + if not files: + raise FileNotFoundError(f"No tube meshes matching {pattern} in {out_dir}") + ids = [] + for f in files: + m = re.search(r"_tube_(\d+)_in\.e$", f.name) + if m: + ids.append(int(m.group(1))) + return sorted(ids) + + +def extract_thm_pressures_to_dat(thm_exo_path, out_dir): + """ For each panel_X/fch_tube_Y block in the THM exodus, write + panel_X_fch_tube_Y_p.dat (AXIS Z / AXIS T / DATA, pressure in MPa). + + Args: + thm_exo_path: path to the THM exodus file + out_dir: directory where the .dat files will be written + """ + out_dir = Path(out_dir) + model = exo.exodus(str(thm_exo_path), array_type="numpy") + times = model.get_times() + for blk_id in model.get_elem_blk_ids(): + blk_name = model.get_elem_blk_name(blk_id) + m = FCH_TUBE_BLOCK_RE.match(blk_name) + if not m: + continue + panel, tube = int(m.group(1)), int(m.group(2)) + conn, num_elem, num_nodes = model.get_elem_connectivity(blk_id) + conn = np.array(conn, dtype=int).reshape((num_elem, num_nodes)) + z_vals = np.array([ + np.mean([model.get_coord(nid)[2] for nid in elem_nodes]) + for elem_nodes in conn + ]) + order = np.argsort(z_vals) + out_path = out_dir / f"panel_{panel}_fch_tube_{tube}_p.dat" + with open(out_path, "w", encoding="utf-8") as f: + f.write(f"# Pressure from block {blk_name} variable p\n") + f.write("AXIS Z\n") + f.write(" ".join(f"{v:.16g}" for v in z_vals[order]) + "\n") + f.write("AXIS T\n") + f.write(" ".join(f"{t:.16g}" for t in times) + "\n") + f.write("DATA\n") + for step in range(len(times)): + p_vals = model.get_variable_values("EX_ELEM_BLOCK", blk_id, "p", step + 1) + f.write(" ".join(f"{v:.16g}" for v in np.array(p_vals)[order] / 1e6) + "\n") + model.close() + + +def create_moose_sm_inputs(moose_thm_filename, out_dir=None): + """ + Creates MOOSE Solid Mechanics input files for each receiver panel based on the THM results. + Expects the THM Exodus files to be named {moose_thm_filename}_flowpath_{fp}_exo.e + + Args: + moose_thm_filename (String): base filename of the MOOSE THM Exodus outputs + out_dir (String, optional): directory where THM exodus files are located. + This is also where the structural input files will be written. Defaults to + current working directory. + + Returns: + (input_paths, output_exodus_paths): input_paths are the .i files for moose sm + output_exodus_paths are the .e files for compute_moose_reliability. + """ + out_dir = Path(out_dir) if out_dir is not None else Path.cwd() + input_paths = [] + output_exodus_paths = [] + for fp in (0, 1): + exo_path = out_dir / f"{moose_thm_filename}_flowpath_{fp}_exo.e" + times = read_time_axis(exo_path) + extract_thm_pressures_to_dat(exo_path, out_dir) + for panel in panels_in_thm_exodus(exo_path): + tubes = discover_tubes(panel, out_dir) + root, i_name = build_structural_input( + panel, tubes, fp, times, out_dir, moose_thm_filename + ) + input_path = out_dir / i_name + pyhit.write(str(input_path), root) + input_paths.append(input_path) + # build_structural_input writes the [Outputs/exodus_out] file_base + # as panel_{P}_struct_from_fp{F}; MOOSE appends .e + output_exodus_paths.append( + out_dir / f"panel_{panel}_struct_from_fp{fp}.e" + ) + return input_paths, output_exodus_paths + +def run_moose_sm_model(moose_input_filename): + """ + Runs the MOOSE SolidMechanics module + Needs environment variables MOOSE_MPI, MOOSE_NPROCS, and DEER to be set. + + Args: + moose_input_filename (String): filename to call moose with + + """ + try: + print("Running MOOSE!") + mpirun = os.environ.get("MOOSE_MPI") + nprocs = os.environ.get("MOOSE_NPROCS") + moose_sm = os.environ.get("NEMLAPP") + argv = [mpirun, "-n", nprocs, moose_sm, "-i", moose_input_filename] + result = subprocess.run(argv, + check=True, capture_output=False, text=True) + except subprocess.CalledProcessError as e: + print(f"MOOSE returned error {e.returncode}") + print(f"stderr: {e.stderr}") + +# Disabling pylint warnings for local variables and statements. This is the main moose file writer +# function, its bound to be long and have many variables. +# pylint: disable=too-many-locals, too-many-statements +def build_structural_input(panel: int, tubes: list, flowpath: int, + times, out_dir: Path, moose_thm_filename: str): + """ + Build the structural solution file for receiver panel by panel. + + Args: + panel: int, panel number + tubes: list, list of tube IDs + flowpath: int, flow path number + times: list, time values + out_dir: Path, directory where the structural input file will be written + moose_thm_filename: str, base filename of the MOOSE THM Exodus outputs + + Returns: + root: pyhit.Node, used to write the MOOSE input file + i_name: str, the name of the structural input file + """ + thm_file = f"{moose_thm_filename}_flowpath_{flowpath}_exo.e" + output_base = f"panel_{panel}_struct_from_fp{flowpath}" + i_name = f"moose_structural_panel_{panel}_from_fp{flowpath}.i" + dt = float(times[1] - times[0]) if len(times) > 1 else 1.0 + end_time = float(times[-1]) + + root = pyhit.Node(parent=None, hitnode=None, offset=None) + + mesh = root.append("Mesh", construct_side_list_from_node_list="true") + for t in tubes: + rmin_id = 100 + t + blk_id = 10 + t + mesh.append(f"tube_{t}_mesh", + type="FileMeshGenerator", + file=f"panel_{panel}_tube_{t}_in.e") + mesh.append(f"tube_{t}_renum", + type="RenameBoundaryGenerator", + input=f"tube_{t}_mesh", + old_boundary=make_moose_hit_vector(["rmin"]), + new_boundary=make_moose_hit_vector([rmin_id])) + mesh.append(f"tube_{t}_rename", + type="RenameBoundaryGenerator", + input=f"tube_{t}_renum", + old_boundary=make_moose_hit_vector([rmin_id]), + new_boundary=make_moose_hit_vector([f"rmin_tube_{t}"])) + mesh.append(f"tube_{t}_blk_renum", + type="RenameBlockGenerator", + input=f"tube_{t}_rename", + old_block=make_moose_hit_vector(["0"]), + new_block=make_moose_hit_vector([blk_id])) + mesh.append(f"tube_{t}_blk_rename", + type="RenameBlockGenerator", + input=f"tube_{t}_blk_renum", + old_block=make_moose_hit_vector([blk_id]), + new_block=make_moose_hit_vector([f"tube_{t}"])) + mesh.append("combined", + type="CombinerGenerator", + inputs=make_moose_hit_vector([f"tube_{t}_blk_rename" for t in tubes])) + + prev_input = "combined" + for t in tubes: + coord_a, coord_b = find_pin_coords(out_dir / f"panel_{panel}_tube_{t}_in.e") + mesh.append(f"pin_xy_tube_{t}", + type="ExtraNodesetGenerator", + input=prev_input, + new_boundary=make_moose_hit_vector([f"pin_xy_tube_{t}"]), + coord=make_moose_hit_vector( + [f"{coord_a[0]:.16g}", f"{coord_a[1]:.16g}", f"{coord_a[2]:.16g}"])) + prev_input = f"pin_xy_tube_{t}" + mesh.append(f"pin_y_tube_{t}", + type="ExtraNodesetGenerator", + input=prev_input, + new_boundary=make_moose_hit_vector([f"pin_y_tube_{t}"]), + coord=make_moose_hit_vector( + [f"{coord_b[0]:.16g}", f"{coord_b[1]:.16g}", f"{coord_b[2]:.16g}"])) + prev_input = f"pin_y_tube_{t}" + + root.append("GlobalParams", + displacements=make_moose_hit_vector(["disp_x", "disp_y", "disp_z"])) + + aux_vars = root.append("AuxVariables") + aux_vars.append("temp", order="FIRST", family="LAGRANGE") + + user_objs = root.append("UserObjects") + user_objs.append("tube_temp_soln", + type="SolutionUserObject", + mesh=thm_file, + system_variables="T_solid", + execute_on=make_moose_hit_vector(["initial", "timestep_begin"])) + + functions = root.append("Functions") + for t in tubes: + functions.append(f"p_from_thm_tube_{t}", + type="PiecewiseMultilinear", + data_file=f"panel_{panel}_fch_tube_{t}_p.dat") + + physics = root.append("Physics") + sm = physics.append("SolidMechanics") + qs = sm.append("QuasiStatic") + qs.append("all", + strain="SMALL", + new_system="true", + eigenstrain_names="eigenstrain", + add_variables="true", + generate_output=make_moose_hit_vector([ + "cauchy_stress_xx", "cauchy_stress_yy", "cauchy_stress_zz", + "cauchy_stress_yz", "cauchy_stress_xz", "cauchy_stress_xy", + "mechanical_strain_xx", "mechanical_strain_yy", "mechanical_strain_zz", + "mechanical_strain_yz", "mechanical_strain_xz", "mechanical_strain_xy", + ])) + + aux_kernels = root.append("AuxKernels") + aux_kernels.append("temp_from_thm", + type="SolutionAux", + variable="temp", + solution="tube_temp_soln", + from_variable="T_solid", + execute_on=make_moose_hit_vector(["initial", "timestep_begin"])) + + bcs = root.append("BCs") + bcs.append("z_disp", + type="DirichletBC", + variable="disp_z", + boundary=make_moose_hit_vector(["bot"]), + value=0.0) + for t in tubes: + bcs.append(f"pin_x_tube_{t}", + type="DirichletBC", + variable="disp_x", + boundary=make_moose_hit_vector([f"pin_xy_tube_{t}"]), + value=0.0) + bcs.append(f"pin_y_tube_{t}", + type="DirichletBC", + variable="disp_y", + boundary=make_moose_hit_vector([f"pin_xy_tube_{t}"]), + value=0.0) + bcs.append(f"pin_y2_tube_{t}", + type="DirichletBC", + variable="disp_y", + boundary=make_moose_hit_vector([f"pin_y_tube_{t}"]), + value=0.0) + for t in tubes: + bcs.append(f"inner_pressure_x_tube_{t}", + type="Pressure", + variable="disp_x", + function=f"p_from_thm_tube_{t}", + boundary=f"rmin_tube_{t}") + bcs.append(f"inner_pressure_y_tube_{t}", + type="Pressure", + variable="disp_y", + function=f"p_from_thm_tube_{t}", + boundary=f"rmin_tube_{t}") + + constraints = root.append("Constraints") + constraints.append("ev_z", + type="EqualValueBoundaryConstraint", + variable="disp_z", + secondary=make_moose_hit_vector(["top"]), + penalty=1e7) + + materials = root.append("Materials") + materials.append("stress", + type="CauchyStressFromNEML", + database="../srlife/srlife/data/deformation/SiC.xml", + model="cares", + temperature="temp") + materials.append("thermal_strain", + type="ComputeThermalExpansionEigenstrainNEML", + database="../srlife/srlife/data/deformation/SiC.xml", + model="cares", + temperature="temp", + stress_free_temperature=300.0, + eigenstrain_name="eigenstrain") + + precond = root.append("Preconditioning") + precond.append("pc", type="SMP", full="true") + + root.append("Executioner", + type="Transient", + solve_type="NEWTON", + l_max_its=100, + l_tol=1e-10, + nl_max_its=15, + nl_rel_tol=1e-6, + nl_abs_tol=1e-8, + automatic_scaling="true", + compute_scaling_once="false", + resid_vs_jac_scaling_param=0.5, + petsc_options=make_moose_hit_vector([ + "-snes_converged_reason", "-ksp_converged_reason", + "-snes_linesearch_monitor"]), + petsc_options_iname=make_moose_hit_vector([ + "-pc_type", "-pc_factor_mat_solver_package"]), + petsc_options_value=make_moose_hit_vector(["lu", "superlu_dist"]), + line_search="none", + dt=dt, + start_time=0.0, + end_time=end_time) + + outputs = root.append("Outputs", print_linear_residuals="false") + outputs.append("exodus_out", + type="Exodus", + file_base=output_base, + sync_times=make_moose_hit_vector([f"{float(t):.16g}" for t in times]), + sync_only="true") + outputs.append("csv_out", type="CSV", file_base=output_base) + + return root, i_name + + +def get_element_temperatures(model, conn, i_step): + """ + Element-averaged temperature (K) at exodus 1-based step i_step + + Args: + model: exodus model object + conn: element connectivity matrix for the block (nelem, 8) + i_step: 1-based step index in the exodus file + Returns: + temperatures: element-averaged temperatures (nelem,) + """ + + temp_all = model.get_variable_values("EX_NODAL", 0, "temp", i_step) + return np.mean(temp_all[conn - 1], axis=1) + + +def read_tube_stress_and_temp(model, blk_id, conn, times): + """ Mandel-form Cauchy stress (ntime, nelem, 6) and element-averaged + temperatures (ntime, nelem) for one HEX8 tube block in a MOOSE + structural exodus + + Args: + model: exodus model object + blk_id: block ID for the tube in the exodus file + conn: element connectivity matrix for the block (nelem, 8) + times: list of time values in the exodus file + Returns: + mandel_stress: Cauchy stress in Mandel-form + temperatures: element-averaged temperatures + """ + n_times = len(times) + n_elem = conn.shape[0] + mandel_stress = np.zeros((n_times, n_elem, 6)) + temperatures = np.zeros((n_times, n_elem)) + stress_vars = [ + "cauchy_stress_xx", "cauchy_stress_yy", "cauchy_stress_zz", + "cauchy_stress_yz", "cauchy_stress_xz", "cauchy_stress_xy", + ] + mandel_mult = np.array([1.0, 1.0, 1.0, SQRT2, SQRT2, SQRT2]) + for t_idx in range(n_times): + i_step = t_idx + 1 + for s_idx, var_name in enumerate(stress_vars): + s_vals = model.get_variable_values( + "EX_ELEM_BLOCK", blk_id, var_name, i_step + ) + mandel_stress[t_idx, :, s_idx] = s_vals * mandel_mult[s_idx] + temperatures[t_idx] = get_element_temperatures(model, conn, i_step) + return mandel_stress, temperatures + + +def compute_moose_reliability(rec, mat_damage, damage_model, lifetime, + moose_sm_output_files, tube_multiplier): + """Read MOOSE structural exodus output and returns reliability using srlife's models. + + Args: + rec: srlife Receiver. + mat_damage: material damage model variant. + damage_model: srlife damage model (e.g. PIAModel). + lifetime: float, hours. + moose_sm_output_files: list of paths to MOOSE SM exodus files + tube_multiplier: float, scaling from analysis tubes to actual + tubes for per-panel + + Returns: + a dict that the driver script can use: + - "lifetime": input lifetime in hours + - "tube_volume": list of per-tube volume flaw reliabilities(VFR) + - "tube_surface": list of per-tube surface flaw reliabilities (SFR) + - "tube_combined": list of per-tube combined reliabilities (CR) + - "panel_volume": list of per-panel VFR + - "panel_surface": list of per-panel SFR + - "panel_combined": list of per-panel CR + - "overall_volume": overall VFR + - "overall_surface": overall SFR + - "overall_combined": overall CR + + """ + m3_to_mm3 = convert_m_to_mm(1.0) ** 3 + m2_to_mm2 = convert_m_to_mm(1.0) ** 2 + + # One representative tube + sample_tube = next(iter(next(iter(rec.panels.values())).tubes.values())) + nt, nz = sample_tube.nt, sample_tube.nz + + volumes = sample_tube.element_volumes() * m3_to_mm3 + surface, normals = sample_tube.surface_elements() + + # element_surface_areas() ships (z, t) per side; surface mask is (r, t, z) + # with z fastest -- transpose each side so areas align with the mask. + sa_raw = sample_tube.element_surface_areas() + half = (nz - 1) * nt + inner_sa = sa_raw[:half].reshape(nz - 1, nt).T.flatten() + outer_sa = sa_raw[half:].reshape(nz - 1, nt).T.flatten() + surface_areas = np.concatenate([inner_sa, outer_sa]) * m2_to_mm2 + + # Reorder so surface elements come first -- workaround for srlife's + # damage.py [:count_surface_elements] slicing in the surface-flaw call. + sort_order = np.argsort(~surface) + volumes_r = volumes[sort_order] + surface_r = surface[sort_order] + normals_r = normals[sort_order] + + per_panel_tube_results = [] + for sm_exo_path in moose_sm_output_files: + model = exo.exodus(str(sm_exo_path), array_type="numpy") + times = np.asarray(model.get_times()) + times_hr = times / 3600.0 # in moose solution the time is in seconds, here we need hours. + tube_results = [] + for blk_id in model.get_elem_blk_ids(): + conn_flat, num_elem, npe = model.get_elem_connectivity(blk_id) + if npe != 8: + continue + conn = np.array(conn_flat, dtype=int).reshape(num_elem, npe) + mandel_stress, temperatures = read_tube_stress_and_temp( + model, blk_id, conn, times + ) + mandel_stress = mandel_stress[:, sort_order] + temperatures = temperatures[:, sort_order] + vol_log_rel = damage_model.calculate_volume_flaw_element_log_reliability( + times_hr, mandel_stress, temperatures, volumes_r, + mat_damage, lifetime, + ) + surf_log_rel = damage_model.calculate_surface_flaw_element_log_reliability( + times_hr, mandel_stress, surface_r, normals_r, + temperatures, surface_areas, mat_damage, lifetime, + ) + combined_log_rel = vol_log_rel.copy() + combined_log_rel[np.where(surface_r)[0]] += surf_log_rel + tube_results.append({ + "volume": float(np.sum(vol_log_rel)), + "surface": float(np.sum(surf_log_rel)), + "combined": float(np.sum(combined_log_rel)), + }) + model.close() + per_panel_tube_results.append(tube_results) + + all_vol = np.array([r["volume"] for tubes in per_panel_tube_results for r in tubes]) + all_surf = np.array([r["surface"] for tubes in per_panel_tube_results for r in tubes]) + all_comb = np.array([r["combined"] for tubes in per_panel_tube_results for r in tubes]) + + panel_volume, panel_surface, panel_combined = [], [], [] + idx = 0 + for tubes in per_panel_tube_results: + n = len(tubes) + panel_volume.append(float(np.exp(np.sum(all_vol[idx:idx + n] * tube_multiplier)))) + panel_surface.append(float(np.exp(np.sum(all_surf[idx:idx + n] * tube_multiplier)))) + panel_combined.append(float(np.exp(np.sum(all_comb[idx:idx + n] * tube_multiplier)))) + idx += n + + return { + "lifetime": lifetime, + "tube_volume": [float(v) for v in np.exp(all_vol)], + "tube_surface": [float(v) for v in np.exp(all_surf)], + "tube_combined": [float(v) for v in np.exp(all_comb)], + "panel_volume": panel_volume, + "panel_surface": panel_surface, + "panel_combined": panel_combined, + "overall_volume": float(np.exp(np.sum(all_vol * tube_multiplier))), + "overall_surface": float(np.exp(np.sum(all_surf * tube_multiplier))), + "overall_combined": float(np.exp(np.sum(all_comb * tube_multiplier))), + } diff --git a/srlife/receiver.py b/srlife/receiver.py index fc43b01..1859eac 100644 --- a/srlife/receiver.py +++ b/srlife/receiver.py @@ -8,14 +8,44 @@ import itertools from collections import OrderedDict - +import subprocess +import os +import sys import numpy as np import scipy.interpolate as inter import h5py - from srlife import writers + +# Get absolute paths to moose python modules +current_dir = os.path.dirname(os.path.abspath(__file__)) +repo_root = os.path.dirname(current_dir) # go up one level +MOOSE_ROOT_PATH = os.path.join(repo_root, 'moose') # go to moose directory +# Add moose python paths for pyhit +moose_python_paths = [ + os.path.join(MOOSE_ROOT_PATH, 'python'), + os.path.join(MOOSE_ROOT_PATH, 'python', 'pyhit'), + os.path.join(MOOSE_ROOT_PATH, 'framework', 'contrib', 'hit'), +] + +for _p in moose_python_paths: + if os.path.exists(_p) and _p not in sys.path: + sys.path.insert(0, _p) + + +import pyhit # pylint: disable=import-error,wrong-import-position +from pyhit import moosetree # pylint: ddisable=import-error,wrong-import-position +conda_env_dir = os.environ.get("CONDA_PREFIX") +# Needed to use exodus.py +ACCESS = os.getenv("ACCESS", f"{conda_env_dir}/seacas") +sys.path.append(os.path.join(ACCESS, "lib")) +sys.path.append(os.path.join(ACCESS, "lib64")) +import exodus as exo # pylint: disable=import-error,wrong-import-position + + + + class Receiver: """Basic definition of the tubular receiver geometry. @@ -243,6 +273,682 @@ def load(cls, fobj): return res + def create_moose_thm_inlet( + self, moose_node, name, connectivity, m_dot_in, fluid_inlet_T + ): + """ + Create an inlet object for a moose THM analysis. + + Args: + node (pyhit.Node): the node object to append the inlet to + typically this would be a node for a flowpath in the Components + section of MOOSE input file + name (string): name for inlet, want to include panel info to make it searchable + connectivity (string): The moose path to the node the inlet connects to + m_dot_in (double): input mass flow rate to inlet + fluid_inlet_T (double): input fluid temp for inlet + + Return: None + """ + print("Creating Moose Inlet!!!") + moose_node.append( + name, + type="InletMassFlowRateTemperature1Phase", + input=connectivity, + m_dot=convert_kghr_to_kgs(m_dot_in), + T=fluid_inlet_T, + ) + + def create_moose_thm_outlet(self, moose_node, name, connectivity, outlet_p): + """ + Create an outlet object for a moose THM analysis. + + Args: + node (pyhit.Node): the node object to append the inlet to + typically this would be a node for a flowpath in the Components + section of MOOSE input file + name (string): name for outlet, want to include panel info to make it searchable + connectivity (string): The moose path to the node the outlet connects to + outlet_p (double): outlet pressure + + Return: None + """ + print("Creating Moose outlet!!!") + moose_node.append(name, type="Outlet1Phase", input=connectivity, p=outlet_p) + + def create_moose_thm_panel_to_panel_connection( + self, comp_node, panel_node, prev_panel_node, manifold_tube, tube_roughness, pos + ): + """ + Create a tube and connection between panels. This connects the + current panel being built to the previously built panel. + + Args: + comp_node (pyhit.Node): pyhit node for components section + panel_node (pyhit.Node): pyhit node for current panel + prev_panel_node (pyhit.Node): pyhit node for prev panel + manifold_tube (Tube): helper object for manifold tubes + tube_roughness (double): tube roughness used in flow simulation + pos (string): sets top or bottom connection + """ + # find in/out tubes from panels + prev_panel_tube_name = f"fch_{prev_panel_node.name}_" + panel_tube_name = f"fch_{panel_node.name}_" + if pos == "bot": + in_out_string = "in" + else: + in_out_string = "out" + prev_panel_tube_name += in_out_string + panel_tube_name += in_out_string + # Get in/out from prev panel + prev_panel_tube = moosetree.find( + prev_panel_node, func=lambda n: n.name == prev_panel_tube_name + ) + # Get in/out from this panel + panel_tube = moosetree.find( + panel_node, func=lambda n: n.name == panel_tube_name + ) + if prev_panel_tube is None or panel_tube is None: + print("CANNOT FIND PANEL in/out tube!!!") + print(prev_panel_tube_name) + print(panel_tube_name) + sys.exit() + # get height of this panels tube + in_out_tube_height = float(panel_tube["length"]) + # create tube between panels + # take hit vector to np.array of floats + tube_pos = str(prev_panel_tube["position"]) + tube_start = np.fromstring(tube_pos, dtype=float, sep=" ") + tube_pos = str(panel_tube["position"]) + tube_end = np.fromstring(tube_pos, dtype=float, sep=" ") + # adjust connector tubes to match correct height for top tubes + if pos == "top": + tube_start[2] += in_out_tube_height + tube_end[2] += in_out_tube_height + orientation = tube_end - tube_start + dist = np.sqrt(np.dot(orientation, orientation)) + tube_name = f"fch_{prev_panel_node.name}_to_{panel_node.name}" + manifold_tube_ir = convert_mm_to_m(manifold_tube.r - manifold_tube.t) + A_pipe = np.pi * manifold_tube_ir**2 + comp_node.append( + tube_name, + type="FlowChannel1Phase", + position=make_moose_hit_vector(tube_start), + orientation=make_moose_hit_vector(orientation), + n_elems=manifold_tube.nz, + length=dist, + A=A_pipe, + D_h=2.0 * manifold_tube_ir, + roughness=tube_roughness, + fp="fp", + ) + # create jcts between tubes + connectivity = [ + f"{prev_panel_node.name}/{prev_panel_tube_name}:{in_out_string}", + f"{tube_name}:in", + ] + comp_node.append( + f"jct_{prev_panel_node.name}_to_connector", + type="VolumeJunction1Phase", + position=make_moose_hit_vector(tube_start), + volume=A_pipe * 2.0, + connections=make_moose_hit_vector(connectivity), + ) + connectivity = [ + f"{tube_name}:out", + f"{panel_node.name}/{panel_tube_name}:{in_out_string}", + ] + comp_node.append( + f"jct_connector_to_{panel_node.name}", + type="VolumeJunction1Phase", + position=make_moose_hit_vector(tube_end), + volume=A_pipe * 2.0, + connections=make_moose_hit_vector(connectivity), + ) + + def create_moose_thm_front_matter(self, moose_root, press, fluid_inlet_T): + """ + Create front matter before components in MOOSE input file + This includes GlobalParams, materials, functions, closures, etc + + Args: + moose_root (pyhit.Node): root node for pyhit moose interface + press (double): outlet pressure for flowpath + fluid_inlet_T (double): inlet fluid temp, used to set initial_T + of all components and fluid + + ToDo: Use solver, fluid, material, etc objects + """ + # global params + moose_root.append( + "GlobalParams", + initial_p=press, + initial_vel=0.0, + initial_T=fluid_inlet_T, + initial_vel_x=0.0, + initial_vel_y=0.0, + initial_vel_z=0.0, + gravity_vector=make_moose_hit_vector([0, 0, -9.81]), + rdg_slope_reconstruction="full", + scaling_factor_1phase=make_moose_hit_vector([1, 1e-2, 1e-4]), + closures="thm_closure", + ) + # fp = "fp") + # material properties + mat_node = moose_root.append("Materials") + + mat_node.append( + "tube_mat", + type="ADGenericConstantMaterial", + prop_names=make_moose_hit_vector( + ["density", "specific_heat", "thermal_conductivity"] + ), + prop_values=make_moose_hit_vector([3200, 947, 60]), + ) + # fluid properties + fluid_node = moose_root.append("FluidProperties") + fluid_node.append( + "sco2", type="CO2FluidProperties", allow_imperfect_jacobians="true" + ) + fluid_node.append( + "fp", + type="TabulatedBicubicFluidProperties", + fp="sco2", + interpolated_properties=make_moose_hit_vector( + [ + "density", + "enthalpy", + "viscosity", + "internal_energy", + "k", + "cv", + "cp", + "c", + "entropy", + ] + ), + # out_of_bounds_behavior="declare_invalid", + temperature_min=240, + temperature_max=999, + pressure_min=1e7, + pressure_max=8e7, + num_T=100, + num_p=100, + tolerance=1e-5, + T_initial_guess=773, + p_initial_guess=2e7, + construct_pT_from_ve="true", + construct_pT_from_vh="true", + allow_imperfect_jacobians="false", + ) + # fluid_node.append("fp", + # type="SimpleFluidProperties", + # density0=1500, + # bulk_modulus=2.0e9, + # porepressure_coefficient=0, + # molar_mass=0.0812, + # specific_entropy=0.0, + # thermal_expansion=1.0e-12, + # cv=1008, + # cp=1010, + # viscosity=4.0e-3, + # thermal_conductivity=0.44) + + # Closures + closure_node = moose_root.append("Closures") + closure_node.append("thm_closure", type="Closures1PhaseTHM") + + # Functions + moose_root.append("Functions") + # user objects needed for flux SolutionFunctions + moose_root.append("UserObjects") + + def create_moose_thm_back_matter( + self, + moose_root, + flow_path_name, + target_outlet_T, + mass_flow_with_t, + inlet_name, + inlet_pipe_name, + outlet_pipe_name, + start_time, + end_time, + dt, + dtmin, + dtmax, + nl_rel_tol, + nl_abs_tol, + nl_max_its, + ): + """ + Add back matter to moose_root node. + BCs, controls, postProc, precon, execs, outputs + + Args: + moose_root (pyhit.Node): pyhit node for the root of moose sim + flow_path_name (string): flowpath name used to set filename for csv outs + target_outlet_T (double): target temperature for outlet, used in controls + mass_flow_with_t (np.array): mass flow rate with time + inlet_name (string): name of pyhit node for inlet to flowpath + inlet_pipe_name (string): name of pipe after inlet + outlet_pipe_name (string): name of pipe before outlet + start_time (double): sim start time (sec) + end_time (double): sim end time (sec) + dt (double): init time step (sec) + dtmin (double): min time step (sec) + dtmax (double): max time step (sec) + nl_rel_tol (double): relative tolerance for newton iter + nl_abs_tol (double): absolute tolerance for newton iter + nl_max_its (int): maximum number of iterations for newton solver + + """ + # Controls + useControls = True + if useControls: + m_dot_fun = "m_dot_time_fun" + func_node = moosetree.find(moose_root, func=lambda n: n.name == "Functions") + if func_node is None: + func_node = moose_root.append("Functions") + times = np.arange(start_time, end_time, dtmax) + times = np.append(times, end_time) + func_node.append( + m_dot_fun, + type="PiecewiseLinear", + x=make_moose_hit_vector(times), + y=make_moose_hit_vector(mass_flow_with_t), + ) + control_node = moose_root.append("ControlLogic") + control_node.append( + "set_inlet_mass_flow", + type="TimeFunctionComponentControl", + component=inlet_name, + parameter="m_dot", + function=m_dot_fun, + ) + + # post processors + post_proc_node = moose_root.append("Postprocessors") + post_proc_node.append( + "m_dot_inlet", + type="RealComponentParameterValuePostprocessor", + component=inlet_name, + parameter="m_dot", + ) + post_proc_node.append( + "path_T_in", type="SideAverageValue", boundary=inlet_pipe_name, variable="T" + ) + post_proc_node.append( + "path_T_out", + type="SideAverageValue", + boundary=outlet_pipe_name, + variable="T", + ) + post_proc_node.append( + "core_p_in", type="SideAverageValue", boundary=inlet_pipe_name, variable="p" + ) + post_proc_node.append( + "core_p_out", + type="SideAverageValue", + boundary=outlet_pipe_name, + variable="p", + ) + post_proc_node.append( + "core_delta_p", + type="ParsedPostprocessor", + pp_names=make_moose_hit_vector(["core_p_in", "core_p_out"]), + function=make_moose_hit_vector(["core_p_in - core_p_out"]), + ) + + # PRECONDITIONER + precon_node = moose_root.append("Preconditioning") + precon_node.append( + "pc", + type="SMP", + full="true", + petsc_options_iname=make_moose_hit_vector(["-snes_test_err"]), + petsc_options_value=make_moose_hit_vector([1.0e-9]), + ) + + # EXECUTIONER + exec_node = moose_root.append( + "Executioner", + type="Transient", + start_time=start_time, + dtmin=dtmin, + dtmax=dtmax, + end_time=end_time, + line_search="basic", + solve_type="NEWTON", + petsc_options=make_moose_hit_vector( + [ + "-snes_converged_reason", + "-ksp_converged_reason", + "-snes_linesearch_monitor", + ] + ), + petsc_options_iname=make_moose_hit_vector( + ["-pc_type", "-pc_factor_mat_solver_package"] + ), + petsc_options_value=make_moose_hit_vector(["lu", "superlu_dist"]), + l_max_its=200, + l_tol=1.0e-10, + nl_rel_tol=nl_rel_tol, + nl_abs_tol=nl_abs_tol, + nl_max_its=nl_max_its, + ) + exec_node.append( + "TimeStepper", type="IterationAdaptiveDT", dt=dt, growth_factor=4 + ) + exec_node.append("TimeIntegrator", type="BDF2") + # OUTPUTS + sync_times = np.arange(start_time, end_time + dtmax, dtmax) + output_node = moose_root.append("Outputs", print_linear_residuals="false") + output_node.append( + "exo", + type="Exodus", + sync_times=make_moose_hit_vector(sync_times), + sync_only="true", + ) + output_node.append( + "console", type="Console", max_rows=1, outlier_variable_norms="false" + ) + output_node.append( + "csv", + type="CSV", + file_base=f"{flow_path_name}_outs", + delimiter=",", + show=make_moose_hit_vector(["core_delta_p", "m_dot_inlet", "path_T_out"]), + ) + + def create_moose_thm_model( + self, + moose_filename, + rec_diam, + tube_od, + tube_spacing, + tube_roughness, + manifold_tube, + outlet_p, + target_outlet_T, + start_time, + end_time, + dt, + nl_rel_tol, + nl_abs_tol, + ): + """ + Create a moose THM input file for the receiver. This file + can then be used with MOOSE to run the Thermal Hydraulics analysis. + + This assumes the receiver is fully defined and ready for analysis + + NOTE: To work here we must assume panel ordering. This allows + us to put the panels/tubes in correct positions + + Args: + moose_filename (string): filename used for the moose input + rec_diam (double): diameter of receiver (THIS SHOULD BE MEMBER OF REC) + tube_od (double): diameter of receiver tubes + tube_spacing (double): space between tubes, for edge spacing to tubes + tube_roughness (double): tube roughness used in flow simulation + manifold_tube (Tube): srlife tube objects for manifold pipes + outlet_p (double): outlet pressure for flow + target_outlet_T (double): target fluid temp at outlet + start_time (double): sim start time + end_time (double): sim end time + dt (double): init time step + nl_rel_tol (double): relative tolerance for newton iter + nl_abs_tol (double): absolute tolerance for newton iter + + Returns: None + """ + print("Creating Moose Model!!!") + # simulation specifics + dtmax = 3600 + dtmin = 1 + nl_max_its = 40 + # length of in/out tubes to panels + panel_in_out_length = 0.1 + + # This loop will create a moose input file for each flowpath in receiver + panel_delta_theta = 2.0 * np.pi / self.npanels + rec_radius = 0.5 * convert_mm_to_m(rec_diam) + edge_spacing_angle = convert_mm_to_m(tube_spacing + 0.5 * tube_od) / ( + rec_radius + ) + filenames = [] + for path_key, flowpath in self.flowpaths.items(): + # for each flowpath in receiver + # make pyhit root + # loop over panels and call their functions + moose_root = pyhit.Node(parent=None, hitnode=None, offset=None) + # add moose front matter: Global params, fluid props, closures, functions, etc + self.create_moose_thm_front_matter( + moose_root, outlet_p, flowpath["inlet_temp"][0] + ) + # Component loop: This will create all components needed in MOOSE + comp_node = moose_root.append("Components") + # NOTE: MOOSE THM cannot do groups of groups + # so I cannot create a flowpath group + print(f"path: {path_key}") + connectivity = ( + f"panel_{flowpath['panels'][0]}/fch_panel_{flowpath['panels'][0]}_in:in" + ) + inlet_name = f"inlet_panel_{flowpath['panels'][0]}" + self.create_moose_thm_inlet( + comp_node, + inlet_name, + connectivity, + flowpath["mass_flow"][0], + flowpath["inlet_temp"][0], + ) + for iPanel, panel_name in enumerate(flowpath["panels"]): + print(f"panel: {panel_name}") + # for each panel in flowpath + # get theta start and delta theta + panel_node = comp_node.append(f"panel_{panel_name}") + panel = self.panels[panel_name] + panel_theta_start = ( + int(panel_name) * panel_delta_theta + edge_spacing_angle + ) + panel_theta_end = ( + int(panel_name) + 1 + ) * panel_delta_theta - edge_spacing_angle + panel.create_panel_components( + panel_node, + panel_theta_start, + panel_theta_end, + rec_radius, + tube_roughness, + manifold_tube, + panel_in_out_length, + ) + if iPanel != 0: + # if we are not in first panel, connect panels + prev_panel_name = f"panel_{flowpath['panels'][iPanel-1]}" + prev_panel_node = moosetree.find( + comp_node, func=lambda n, name=prev_panel_name: n.name == name + ) + if iPanel % 2 == 0: + # even panel, connect bot + pos = "bot" + self.create_moose_thm_panel_to_panel_connection( + comp_node, + panel_node, + prev_panel_node, + manifold_tube, + tube_roughness, + pos, + ) + else: + # odd panel, connect top + pos = "top" + self.create_moose_thm_panel_to_panel_connection( + comp_node, + panel_node, + prev_panel_node, + manifold_tube, + tube_roughness, + pos, + ) + connectivity = ( + f"panel_{flowpath['panels'][-1]}/fch_panel_{flowpath['panels'][-1]}_" + ) + if len(flowpath["panels"]) % 2: + # if there are even number of panels in path + outlet_pos = "out:out" + else: + # odd number of panels in path + outlet_pos = "in:in" + connectivity += outlet_pos + outlet_name = f"outlet_panel_{flowpath['panels'][-1]}" + self.create_moose_thm_outlet(comp_node, outlet_name, connectivity, outlet_p) + # create moose back matter: BCs, controls, postProcessors, + # precon, execs, outputs + inlet_pipe_name = ( + f"panel_{flowpath['panels'][0]}/fch_panel_{flowpath['panels'][0]}_in:in" + ) + outlet_pipe_name = connectivity + flow_path_name = f"flowpath_{path_key}" + self.create_moose_thm_back_matter( + moose_root, + flow_path_name, + target_outlet_T, + convert_kghr_to_kgs(flowpath["mass_flow"]), + inlet_name, + inlet_pipe_name, + outlet_pipe_name, + start_time, + end_time, + dt, + dtmin, + dtmax, + nl_rel_tol, + nl_abs_tol, + nl_max_its, + ) + moose_flow_path_filename = f"{moose_filename}_{flow_path_name}.i" + pyhit.write(moose_flow_path_filename, moose_root) + filenames.append(moose_flow_path_filename) + return filenames + + def run_moose_thm_model(self, moose_input_filename): + """ + Runs the MOOSE THM model using inputs + + Args: + moose_exec (String): moose executable path and name + moose_input_filename (String): filename to call moose with + + """ + try: + print("Running MOOSE!") + mpirun = os.environ.get("MOOSE_MPI") + nprocs = os.environ.get("MOOSE_NPROCS") + moose_thm = os.environ.get("MOOSE_THM") + argv = [mpirun, "-n", nprocs, moose_thm, "-i", moose_input_filename] + subprocess.run(argv,check=True, capture_output=False, text=True) + except subprocess.CalledProcessError as e: + print(f"MOOSE returned error {e.returncode}") + print(f"stderr: {e.stderr}") + + def get_moose_thm_results(self, moose_input_filenames): + """ + This will hop into MOOSE exodus file and get + flow tube nodal values of pressure and + heat tube nodal values of temp + Then load them into srlife Tube objects + + Agrs: + moose_input_filename (list(string)): name of moose input file used to get output + There will be one filename per flowpath + """ + # Steps: + # 1) read MOOSE data and get dicts + # using exodus.py from SEACAS + # 2) load MOOSE data into Tubes + # using add_results + # probably better to do this through tube objects + # loop over tubes, get results from name/dict + # have tube add results + # can loop over panels in rec (get panel_node) + # then tube in panels, find heat_tube? or just use name from there? + # name = {panel_node.name}/heat_tube_i for i in panel.tubes + for iPath, flowpath in enumerate(self.flowpaths.values()): + # for each flowpath in model + input_filename = moose_input_filenames[iPath] + output_filename = input_filename[0:-2] + "_exo.e" + # read temp and pressure results from output exodus + times, press_results, temp_results = read_moose_thm_exodus_file( + output_filename + ) + for panel_name in flowpath["panels"]: + # for each panel in receiver + panel_str = f"panel_{panel_name}/" + panel = self.panels[panel_name] + for iTube, tube in enumerate(panel.tubes.values()): + # for each tube in panel + flow_tube_name = f"{panel_str}fch_tube_{iTube}" + press_res_time_space = press_results[flow_tube_name] + heat_tube_name = f"{panel_str}heat_tube_{iTube}:0" + heat_res_time_space = temp_results[heat_tube_name] + # load MOOSE results into Tube object + tube.load_moose_thm_results( + times, press_res_time_space, heat_res_time_space + ) + + +def read_moose_thm_exodus_file(moose_output_filename): + """ + Read an exodus file and return results + In particular we get dictionary of elem block name and results. + Flow tube nodal pressure results and + Heat tube nodal temperature results + + Args: + moose_output_filename (string): name of moose exodus output file + + Returns: + times (list): list of analysis times + press_results (dict): {elem_block_name, nodal_pressures} pressure results + temp_results (dict): {elem_block_name, nodal_temp} temp results + """ + model = exo.exodus(moose_output_filename, array_type="numpy") + times = model.get_times() + pressure_dict = {} + temp_dict = {} + + for iElemBlk in model.get_elem_blk_ids(): + name = model.get_elem_blk_name(iElemBlk) + if name[0:5] == "panel": + if "fch_tube" in name: + pressure_results = [] + for iStep in range(len(times)): + # NOTE: exodus.py indexes times steps from 1 + press_data = model.get_variable_values( + "EX_ELEM_BLOCK", iElemBlk, "p", iStep + 1 + ) + pressure_results.append(press_data) + pressure_dict[name] = np.array(pressure_results) + if "heat_tube" in name: + connect, _, _ = model.get_elem_connectivity(iElemBlk) + nodes = np.unique(connect) + temp_results = [] + nodeCoord = [] + for node in nodes: + x, y, z = model.get_coord(node) + nodeCoord.append([x, y, z]) + + for iStep in range(len(times)): + # NOTE: exodus.py indexes times steps from 1 + temp_data = model.get_variable_values( + "EX_NODAL", iElemBlk, "T_solid", iStep + 1 + ) + temp_results.append(temp_data[nodes - 1]) + temp_dict[name] = [np.array(nodeCoord), np.array(temp_results)] + return times, pressure_dict, temp_dict + class Panel: """Basic definition of a panel in a tubular receiver. @@ -355,6 +1061,422 @@ def load(cls, fobj): return res + def create_moose_thm_split_connector_tube( + self, + panel_node, + iTube, + manifold_tube, + tube_roughness, + x_tube_prev, + y_tube_prev, + x_tube, + y_tube, + tube_height, + panel_in_out_length, + ): + """ + Create thm tubes to connect tubes within a panel when the connector + tube must be split for panel in/out tubes. This also creates the jcts + and panel in/out pipes + + Args: + panel_node (pyhit.Node): panel node that tubes live on + manifold_tube (Tube): tube object describing the manifold pipes + tube_roughness (double): roughness of tube for flow simulation + x_tube_prev (double): x position of tube before current + y_tube_prev (double): x position of tube before current + x_tube (double): x position of current tube + y_tube (double): x position of current tube + tube_height (double): height of current tube to make connection at + top and bot + is_center_tube (bool): whether or not this tube is in center of panel. + if it is, we will create and join panel in/out tubes + panel_in_out_length (double): length of panel in/out tubes + """ + # Assuming straight path between tubes + tube_to_tube_vect = np.array([x_tube - x_tube_prev, y_tube - y_tube_prev]) + midpoint = 0.5 * np.array([x_tube + x_tube_prev, y_tube + y_tube_prev]) + dist = np.sqrt(np.dot(tube_to_tube_vect, tube_to_tube_vect)) + manifold_tube_ir = convert_mm_to_m(manifold_tube.r - manifold_tube.t) + # bottom tube + bot_connector_tube_1_name = f"fch_tube_{iTube-1}_to_in_bot" + panel_node.append( + bot_connector_tube_1_name, + type="FlowChannel1Phase", + position=make_moose_hit_vector([x_tube_prev, y_tube_prev, 0.0]), + orientation=make_moose_hit_vector( + [x_tube - x_tube_prev, y_tube - y_tube_prev, 0.0] + ), + n_elems=manifold_tube.nz, + length=0.5 * dist, + A=np.pi * manifold_tube_ir**2, + D_h=2.0 * manifold_tube_ir, + roughness=tube_roughness, + fp="fp", + ) + # find and add to jct for prev tube + jct_name = f"jct_tube_{iTube-1}_bot" + add_tube_to_moose_thm_jct( + panel_node, + jct_name, + panel_node.name + "/" + bot_connector_tube_1_name, + "in", + ) + # second half of bottom tube + bot_connector_tube_2_name = f"fch_tube_in_to_{iTube}_bot" + panel_node.append( + bot_connector_tube_2_name, + type="FlowChannel1Phase", + position=make_moose_hit_vector([midpoint[0], midpoint[1], 0.0]), + orientation=make_moose_hit_vector( + [x_tube - x_tube_prev, y_tube - y_tube_prev, 0.0] + ), + n_elems=manifold_tube.nz, + length=0.5 * dist, + A=np.pi * manifold_tube_ir**2, + D_h=2.0 * manifold_tube_ir, + roughness=tube_roughness, + fp="fp", + ) + # find and add to jct for this tube + jct_name = f"jct_tube_{iTube}_bot" + add_tube_to_moose_thm_jct( + panel_node, + jct_name, + panel_node.name + "/" + bot_connector_tube_2_name, + "out", + ) + # create panel in tube and junction + panel_in_tube_name = f"fch_{panel_node.name}_in" + A_pipe = np.pi * manifold_tube_ir**2 + panel_node.append( + panel_in_tube_name, + type="FlowChannel1Phase", + position=make_moose_hit_vector( + [midpoint[0], midpoint[1], -panel_in_out_length] + ), + orientation=make_moose_hit_vector([0.0, 0.0, 1.0]), + n_elems=manifold_tube.nz, + length=panel_in_out_length, + A=A_pipe, + D_h=2.0 * manifold_tube_ir, + roughness=tube_roughness, + fp="fp", + ) + connectivity = [ + f"{panel_node.name}/{panel_in_tube_name}:out", + f"{panel_node.name}/{bot_connector_tube_1_name}:out", + f"{panel_node.name}/{bot_connector_tube_2_name}:in", + ] + panel_node.append( + f"jct_{panel_node.name}_t_in", + type="VolumeJunction1Phase", + position=make_moose_hit_vector([midpoint[0], midpoint[1], 0]), + volume=A_pipe * 2.0, + connections=make_moose_hit_vector(connectivity), + ) + + # first half of top tube + top_connector_tube_1_name = f"fch_tube_{iTube-1}_to_out_top" + panel_node.append( + top_connector_tube_1_name, + type="FlowChannel1Phase", + position=make_moose_hit_vector([x_tube_prev, y_tube_prev, tube_height]), + orientation=make_moose_hit_vector( + [x_tube - x_tube_prev, y_tube - y_tube_prev, 0.0] + ), + n_elems=manifold_tube.nz, + length=0.5 * dist, + A=A_pipe, + D_h=2.0 * manifold_tube_ir, + roughness=tube_roughness, + fp="fp", + ) + # find and add to jct for prev tube + jct_name = f"jct_tube_{iTube-1}_top" + add_tube_to_moose_thm_jct( + panel_node, + jct_name, + panel_node.name + "/" + top_connector_tube_1_name, + "in", + ) + # second half of top tube + top_connector_tube_2_name = f"fch_tube_out_to_{iTube}_top" + panel_node.append( + top_connector_tube_2_name, + type="FlowChannel1Phase", + position=make_moose_hit_vector([midpoint[0], midpoint[1], tube_height]), + orientation=make_moose_hit_vector( + [x_tube - x_tube_prev, y_tube - y_tube_prev, 0.0] + ), + n_elems=manifold_tube.nz, + length=0.5 * dist, + A=A_pipe, + D_h=2.0 * manifold_tube_ir, + roughness=tube_roughness, + fp="fp", + ) + # find and add to jct for this tube + jct_name = f"jct_tube_{iTube}_top" + add_tube_to_moose_thm_jct( + panel_node, + jct_name, + panel_node.name + "/" + top_connector_tube_2_name, + "out", + ) + # create panel out tube and junction + panel_out_tube_name = f"fch_{panel_node.name}_out" + panel_node.append( + panel_out_tube_name, + type="FlowChannel1Phase", + position=make_moose_hit_vector([midpoint[0], midpoint[1], tube_height]), + orientation=make_moose_hit_vector([0.0, 0.0, 1.0]), + n_elems=manifold_tube.nz, + length=panel_in_out_length, + A=A_pipe, + D_h=2.0 * manifold_tube_ir, + roughness=tube_roughness, + fp="fp", + ) + connectivity = [ + f"{panel_node.name}/{top_connector_tube_1_name}:out", + f"{panel_node.name}/{top_connector_tube_2_name}:in", + f"{panel_node.name}/{panel_out_tube_name}:in", + ] + panel_node.append( + f"jct_{panel_node.name}_t_out", + type="VolumeJunction1Phase", + position=make_moose_hit_vector([midpoint[0], midpoint[1], tube_height]), + volume=A_pipe * 2.0, + connections=make_moose_hit_vector(connectivity), + ) + + def create_moose_thm_connector_tube( + self, + panel_node, + iTube, + manifold_tube, + tube_roughness, + x_tube_prev, + y_tube_prev, + x_tube, + y_tube, + tube_height, + panel_in_out_length, + is_center_tube, + ): + """ + Create thm tubes to connect tubes within a panel. It also creates in/out + tubes when the tube is at panel centerline. + + Args: + panel_node (pyhit.Node): panel node that tubes live on + manifold_tube (Tube): tube object describing the manifold pipes + tube_roughness (double): roughness of tube for flow simulation + x_tube_prev (double): x position of tube before current + y_tube_prev (double): x position of tube before current + x_tube (double): x position of current tube + y_tube (double): x position of current tube + tube_height (double): height of current tube to make connection at + top and bot + is_center_tube (bool): whether or not this tube is in center of panel. + if it is, we will create and join panel in/out tubes + panel_in_out_length (double): length of panel in/out tubes + """ + # Assuming straight path between tubes + tube_to_tube_vect = np.array([x_tube - x_tube_prev, y_tube - y_tube_prev]) + dist = np.sqrt(np.dot(tube_to_tube_vect, tube_to_tube_vect)) + manifold_tube_ir = convert_mm_to_m(manifold_tube.r - manifold_tube.t) + # bottom tube + bot_connector_tube_name = f"fch_tube_{iTube-1}_to_{iTube}_bot" + panel_node.append( + bot_connector_tube_name, + type="FlowChannel1Phase", + position=make_moose_hit_vector([x_tube_prev, y_tube_prev, 0.0]), + orientation=make_moose_hit_vector( + [x_tube - x_tube_prev, y_tube - y_tube_prev, 0.0] + ), + n_elems=manifold_tube.nz, + length=dist, + A=np.pi * manifold_tube_ir**2, + D_h=2.0 * manifold_tube_ir, + roughness=tube_roughness, + fp="fp", + ) + # find and add to jct for prev tube + jct_name = f"jct_tube_{iTube-1}_bot" + add_tube_to_moose_thm_jct( + panel_node, jct_name, panel_node.name + "/" + bot_connector_tube_name, "in" + ) + # find and add to jct for this tube + jct_name = f"jct_tube_{iTube}_bot" + add_tube_to_moose_thm_jct( + panel_node, jct_name, panel_node.name + "/" + bot_connector_tube_name, "out" + ) + + # top tube + top_connector_tube_name = f"fch_tube_{iTube-1}_to_{iTube}_top" + panel_node.append( + top_connector_tube_name, + type="FlowChannel1Phase", + position=make_moose_hit_vector([x_tube_prev, y_tube_prev, tube_height]), + orientation=make_moose_hit_vector( + [x_tube - x_tube_prev, y_tube - y_tube_prev, 0.0] + ), + n_elems=manifold_tube.nz, + length=dist, + A=np.pi * manifold_tube_ir**2, + D_h=2.0 * manifold_tube_ir, + roughness=tube_roughness, + fp="fp", + ) + # find and add to jct for prev tube + jct_name = f"jct_tube_{iTube-1}_top" + add_tube_to_moose_thm_jct( + panel_node, jct_name, panel_node.name + "/" + top_connector_tube_name, "in" + ) + # find and add to jct for this tube + jct_name = f"jct_tube_{iTube}_top" + add_tube_to_moose_thm_jct( + panel_node, jct_name, panel_node.name + "/" + top_connector_tube_name, "out" + ) + + if is_center_tube: + panel_in_tube_name = f"fch_{panel_node.name}_in" + panel_node.append( + panel_in_tube_name, + type="FlowChannel1Phase", + position=make_moose_hit_vector([x_tube, y_tube, -panel_in_out_length]), + orientation=make_moose_hit_vector([0.0, 0.0, 1.0]), + n_elems=manifold_tube.nz, + length=panel_in_out_length, + A=np.pi * manifold_tube_ir**2, + D_h=2.0 * manifold_tube_ir, + roughness=tube_roughness, + fp="fp", + ) + # find and add to jct for this tube + jct_name = f"jct_tube_{iTube}_bot" + add_tube_to_moose_thm_jct( + panel_node, jct_name, panel_node.name + "/" + panel_in_tube_name, "out" + ) + panel_out_tube_name = f"fch_{panel_node.name}_out" + panel_node.append( + panel_out_tube_name, + type="FlowChannel1Phase", + position=make_moose_hit_vector([x_tube, y_tube, tube_height]), + orientation=make_moose_hit_vector([0.0, 0.0, 1.0]), + n_elems=manifold_tube.nz, + length=panel_in_out_length, + A=np.pi * manifold_tube_ir**2, + D_h=2.0 * manifold_tube_ir, + roughness=tube_roughness, + fp="fp", + ) + # find and add to jct for this tube + jct_name = f"jct_tube_{iTube}_top" + add_tube_to_moose_thm_jct( + panel_node, jct_name, panel_node.name + "/" + panel_out_tube_name, "in" + ) + + def create_panel_components( + self, + panel_node, + panel_theta_start, + panel_theta_end, + rec_radius, + tube_roughness, + manifold_tube, + panel_in_out_length, + ): + """ + Create the create components of the panel in a MOOSE THM + input file. Uses pyhit to create MOOSE .i files + + NOTE: We have assumed the location of the panel based on its name + + Args: + panel_node (pyhit.Node): the node for this panel in moose components + panel_theta_start (double): rec theta at start of panel (radians) + panel_theta_end (double): rec_theta at end of panel (radians) + rec_radius (double): radius of receiver + tube_roughness (double): roughness of tube for flow simulation + manifold_tube (Tube): tube object created for manifold tubes + panel_in_out_length (double): length of panel in/out tubes + Returns: None + """ + print("Creating Moose Panel!!!") + # set tube thetas based on number of tubes analyzed + thetas = np.linspace(panel_theta_start, panel_theta_end, self.ntubes) + tube_xs = rec_radius * np.cos(thetas) + tube_ys = rec_radius * np.sin(thetas) + panel_center_theta = 0.5 * (panel_theta_start + panel_theta_end) + # This function will create a 3D mesh for the tubes + # to be used for heat transfer + # This is done once per panel, but could be done once per receiver + # if we assume that the tubes are the same on all panels + # and the mesh is also the same + for iTube, tube in enumerate(self.tubes.values()): + x_tube = tube_xs[iTube] + y_tube = tube_ys[iTube] + tube.create_moose_thm_tube_and_jcts( + panel_node, iTube, x_tube, y_tube, tube_roughness + ) + if iTube != 0: + # if this is not the first tube + # connect it to previous tube + tube_height = convert_mm_to_m(tube.h) + x_tube_prev = tube_xs[iTube - 1] + y_tube_prev = tube_ys[iTube - 1] + if ( + thetas[iTube - 1] < panel_center_theta < thetas[iTube] + ): + # this connector crosses centerline + self.create_moose_thm_split_connector_tube( + panel_node, + iTube, + manifold_tube, + tube_roughness, + x_tube_prev, + y_tube_prev, + x_tube, + y_tube, + tube_height, + panel_in_out_length, + ) + + elif thetas[iTube] == panel_center_theta: + # this tube is at panel centerline + # this is normal tubeToTube connection + self.create_moose_thm_connector_tube( + panel_node, + iTube, + manifold_tube, + tube_roughness, + x_tube_prev, + y_tube_prev, + x_tube, + y_tube, + tube_height, + panel_in_out_length, + True, + ) + else: + # this is normal tubeToTube connection + self.create_moose_thm_connector_tube( + panel_node, + iTube, + manifold_tube, + tube_roughness, + x_tube_prev, + y_tube_prev, + x_tube, + y_tube, + tube_height, + panel_in_out_length, + False, + ) + def next_name(names): """Determine the next numeric string name based on a list @@ -374,6 +1496,84 @@ def next_name(names): return str(max(curr_ints) + 1) +def make_moose_hit_vector(list_in): + """ + Convert a python list into a moose hit vector + + Args: list_in (list): the python list to convert + Returns: a vector in moose hit form + """ + return "'" + " ".join(str(v) for v in list_in) + "'" + + +def convert_mm_to_m(mm_in): + """ + Unit conversion from mm to m + + Args: mm_in (double): qty to convert + """ + return mm_in / 1000 + + +def convert_m_to_mm(m_in): + """ + Unit conversion from m to mm + + Args: m_in (double): qty to convert + """ + return m_in * 1000 + + +def convert_kghr_to_kgs(kghr_in): + """ + Unit conversion from kg/hr to kg/sec + + Args: kghr_in (double): qty to convert + """ + return kghr_in / 3600 + + +def convert_Pa_to_MPa(Pa_in): + """ + Unit conversion from Pa to MPa + + Args: Pa_in (double): qty to convert + """ + return Pa_in / 10**6 + + +def convert_Wmm2_to_Wm2(Wmm2_in): + """ + Unit conversion from W/mm^2 to W/m^2 + + Args: Wmm2_in (double): qty to convert + """ + return Wmm2_in * 1000**2 + + +def add_tube_to_moose_thm_jct(panel_node, jct_name, tube_name, tube_in_out): + """ + Add a tube to a junctions connection attribute in moose THM model. + Find the jct on given node, append tube_name to its connections attr + + NOTE: may want to make this just be generic + + Args: + panel_node (pyhit.Node): node to search for junction within + jct_name (string): name of jct to add to + tube_name (string): name of tube to add to jct connections + tube_in_out (string): specifies tube_name:in or :out connection + """ + jct_node = moosetree.find(panel_node, func=lambda n: n.name == jct_name) + if jct_node is None: + print(f"COULD NOT FIND JCT!!! {jct_name}") + # BPMToDo: better error handling here + sys.exit() + connect = jct_node["connections"] + connect = make_moose_hit_vector([connect, tube_name + ":" + tube_in_out]) + jct_node["connections"] = connect + + class Tube: """Geometry, boundary conditions, and results for a single tube. @@ -448,6 +1648,10 @@ def __init__( self.multiplier_val = multiplier + # BPM: for use with MOOSE (default) + self.center_point = [0, 0, 0] + self.tube_yaw = 0.0 + @property def multiplier(self): """ @@ -862,7 +2066,7 @@ def add_quadrature_results(self, name, data): name (str): parameter set name data (np.array): actual results data """ - if data.shape[0] != self.ntime: + if data.shape[0] != self.ntime and name != "ghost_temperature": raise ValueError("Quadrature data must have time axis first!") self.quadrature_results[name] = self._setup_memmap(name + "_quad", data.shape) self.quadrature_results[name][:] = data[:] @@ -885,7 +2089,7 @@ def add_axial_results(self, name, data): name (str): results name data (np.array): data to store """ - if data.shape != (self.ntime, self.nz): + if data.shape != (self.ntime, self.nz) and name[0:5] != "fluid": raise ValueError("Axial result field must have shape (ntime, nz)!") self.axial_results[name] = self._setup_memmap(name + " _axial", data.shape) self.axial_results[name][:] = data[:] @@ -1081,6 +2285,446 @@ def load(cls, fobj): return res + def create_moose_thm_3D_tube_mesh(self, panel_node, tube_num, x_tube, y_tube): + """ + Creates a moose input file for tubes on this panel. + This means making an input file and executing moose + with --mesh-only flag + + Args: + panel_node (pyhit.Node): pyhit node object for this panel + tube_num (int): number designating tube within panel + x_tube (double): x coordinate of tube base + y_tube (double): y coordinate of tube base + + Returns: + tube_mesh_filename (string): filename of this panel's tube mesh + """ + # angle to rotate tube and apply flux BCs to outer face + # tube centerline should be perp to its x,y vector + # this generator splits tube outer face at x = 0 + # so we want to rotate another 90 deg to get sunFaceing + # ASSUMES that center of rec is at 0,0,0 + tube_angle = np.arctan2(y_tube, x_tube) * 180 / np.pi + # set tube center point and yaw + self.center_point = [x_tube, y_tube, 0.0] + # but the yaw of the tube should be just the perp angle + self.tube_yaw = tube_angle * np.pi / 180 - np.pi / 2 + + # create new root for this file + tube_mesh_root = pyhit.Node(parent=None, hitnode=None, offset=None) + mesh_node = tube_mesh_root.append("Mesh") + # add ring mesh + ring_name = "ring_2d" + r_outer = convert_mm_to_m(self.r) + h = convert_mm_to_m(self.h) + mesh_node.append( + ring_name + "_mesh", + type="AnnularMeshGenerator", + nr=self.nr - 1, + nt=self.nt, + rmin=convert_mm_to_m(self.r - self.t), + rmax=r_outer, + ) + # split outer boundary of tube + mesh_node.append( + ring_name, + type="PatchSidesetGenerator", + boundary="rmax", + n_patches=2, + input=ring_name + "_mesh", + ) + # turn ring mesh into tube + mesh_node.append( + "tube", + type="AdvancedExtruderGenerator", + input=ring_name, + heights=make_moose_hit_vector([1, 1, h]), + num_layers=make_moose_hit_vector([0, 0, (self.nz - 1)]), + direction=make_moose_hit_vector([0, 0, 1]), + bottom_boundary="bot", + top_boundary="top", + ) + + # rotate tube + mesh_node.append( + "tube_rot", + type="TransformGenerator", + input="tube", + transform="ROTATE", + vector_value=make_moose_hit_vector([0, 0, tube_angle]), + ) + # transform tube + mesh_node.append( + "tube_pos", + type="TransformGenerator", + input="tube_rot", + transform="TRANSLATE", + vector_value=make_moose_hit_vector([x_tube, y_tube, 0.0]), + ) + tube_mesh_root.append("Outputs", exodus="true") + # write tube input to a file + tube_mesh_moose_input = f"{panel_node.name}_tube_{tube_num}" + pyhit.write(tube_mesh_moose_input + ".i", tube_mesh_root) + # run moose to generate tube mesh + try: + moose_exec = os.environ.get("MOOSE_THM", "MOOSE_THM") + result = subprocess.run([moose_exec, "-i", tube_mesh_moose_input + ".i", "--mesh-only"], + check=True, capture_output=True, text=True) + except subprocess.CalledProcessError as e: + print(f"MOOSE returned error {e.returncode}") + print(f"stderr: {e.stderr}") + + # if this runs, then output file will be below + tube_mesh_file = f"{tube_mesh_moose_input}_in.e" + # now I want to write flux bc data to this exodus file + self.write_flux_bc_to_3d_tube_mesh(tube_mesh_file) + # then create a solution user object and a function to use it + moose_root = (panel_node.parent).parent + user_obj_node = moosetree.find( + moose_root, func=lambda n: n.name == "UserObjects" + ) + if user_obj_node is None: + # is we could not find it, make it + user_obj_node = moose_root.append("UserObjects") + user_obj_node.append( + f"{tube_mesh_moose_input}_SolObj", + type="SolutionUserObject", + mesh=tube_mesh_file, + system_variables="flux", + ) + func_node = moosetree.find(moose_root, func=lambda n: n.name == "Functions") + if func_node is None: + # is we could not find it, make it + func_node = moose_root.append("Functions") + func_node.append( + f"flux_{tube_mesh_moose_input}_SolFn", + type="SolutionFunction", + solution=f"{tube_mesh_moose_input}_SolObj", + ) + return tube_mesh_file + + def write_flux_bc_to_3d_tube_mesh(self, tube_mesh_file): + """ + Writes flux bcs defined in srlife to 3D tube mesh in exodus + This allows moose to generate a solution function to apply time + and space dependent flux bcs + + Args: + tube_mesh_file: exodus filename of tube to write to + """ + model = exo.exodus(tube_mesh_file, array_type="numpy", mode="a") + exo.add_variables(model, nodal_vars=["flux"]) + # then loop over side set nodes with on rmax_0 + rmax_0_id = 2 + _, flux_bc_nodes = model.get_side_set_node_list(rmax_0_id) + if model.get_side_set_name(rmax_0_id) != "rmax_0": + sys.exit("PROBLEM!!!") + # we only want to operate on each node once + flux_bc_nodes = np.unique(flux_bc_nodes) + flux_data = np.zeros((len(self.times), model.num_nodes())) + for node in flux_bc_nodes: + # get coords in MOOSE x,y,z coordinates (expects node index) + x, y, z = model.get_coord(node) + coord = np.array([x, y, z]) + # get srlife tube coords (note: could just map coords) + _, theta, z = self.map_moose_coords_to_srlife_tube_coords(coord) + z_tube = convert_m_to_mm(z) + for iTime, time in enumerate(self.times): + # Get the correct value of flux + flux_data[iTime, node - 1] = convert_Wmm2_to_Wm2( + self.outer_bc.flux(time, theta, z_tube) + ) + # make all times and variables in exodus + for iTime, time in enumerate(self.times): + model.put_time(iTime + 1, time * 3600) + model.put_node_variable_values("flux", iTime + 1, flux_data[iTime]) + # print(model.get_node_variable_values("flux", 2)) + model.close() + + def map_moose_coords_to_srlife_tube_coords(self, coords): + """ + Check if a set of moose node coords is on the outer boundary + of the tube + + Args: + coodrs (list): [x, y, z] in moose analysis coordinates + Returns: + r (double): tube radial coord + theta (double): tube theta coord + z (double): tube z coord + """ + # transform to tube coords + # translate via self.center_point + # rotate by self.tube_yaw + c_yaw = np.cos(self.tube_yaw) + s_yaw = np.sin(self.tube_yaw) + R = np.array([[c_yaw, s_yaw, 0], [-s_yaw, c_yaw, 0], [0, 0, 1]]) + coord_tube = R @ np.array(coords - self.center_point) + x_tube = coord_tube[0] + y_tube = coord_tube[1] + z_tube = coord_tube[2] + # radial coords + r = np.sqrt(x_tube**2 + y_tube**2) + theta = np.arctan2(y_tube, x_tube) + if theta < 0.0: + # we want 0,2pi, not -pi,pi + theta += 2.0 * np.pi + return r, theta, z_tube + + def create_moose_thm_tube_and_jcts( + self, panel_node, tube_num, x_tube, y_tube, tube_roughness + ): + """ + Create a tube node in a moose input file using pyhit. + + Args: + panel_node (pyhit.Node): panel node in moose input file + tube_num (int): number designating tube within panel + x_tube (double): x position of tube + y_tube (double): y position of tube + tube_roughness (double): roughness of tube for flow simulation + """ + tube_ir = convert_mm_to_m(self.r - self.t) + tube_fch_name = f"fch_tube_{tube_num}" + A_pipe = np.pi * (tube_ir) ** 2 + panel_node.append( + tube_fch_name, + type="FlowChannel1Phase", + position=make_moose_hit_vector([x_tube, y_tube, 0]), + orientation=make_moose_hit_vector([0, 0, 1]), + n_elems=(self.nz - 1), + length=convert_mm_to_m(self.h), + A=A_pipe, + D_h=2.0 * (tube_ir), + roughness=tube_roughness, + fp="fp", + ) + # we are just initializing these now + # will tie in with connectors + panel_node.append( + f"jct_tube_{tube_num}_bot", + type="VolumeJunction1Phase", + position=make_moose_hit_vector([x_tube, y_tube, 0]), + volume=A_pipe * 2.0, + connections=make_moose_hit_vector( + [f"{panel_node.name}/fch_tube_{tube_num}:in"] + ), + ) + panel_node.append( + f"jct_tube_{tube_num}_top", + type="VolumeJunction1Phase", + position=make_moose_hit_vector([x_tube, y_tube, convert_mm_to_m(self.h)]), + volume=A_pipe * 2.0, + connections=make_moose_hit_vector( + [f"{panel_node.name}/fch_tube_{tube_num}:out"] + ), + ) + # make 3D tube for heat transfer + panel_tube_mesh_filename = self.create_moose_thm_3D_tube_mesh( + panel_node, tube_num, x_tube, y_tube + ) + # add tube heat structure and heat transfer + # NOTE: since tube mesh is at correct pos, do not need to set position here + heat_tube_name = f"heat_tube_{tube_num}" + panel_node.append( + heat_tube_name, + type="HeatStructureFromFile3D", + file=panel_tube_mesh_filename, + position=make_moose_hit_vector([0.0, 0.0, 0.0]), + ) + panel_node.append( + f"heat_transfer_tube_{tube_num}", + type="HeatTransferFromHeatStructure3D1Phase", + flow_channels=f"{panel_node.name}/{tube_fch_name}", + hs=f"{panel_node.name}/{heat_tube_name}", + boundary=f"{panel_node.name}/{heat_tube_name}:rmin", + P_hf=2 * np.pi * tube_ir, + ) + # add flux BCs for 3D tube + moose_root = (panel_node.parent).parent + bc_node_name = "BCs" + bc_node = moosetree.find(moose_root, func=lambda n: n.name == bc_node_name) + if bc_node is None: + # is we could not find it, make it + bc_node = moose_root.append("BCs") + bc_func_name = f"flux_{panel_node.name}_tube_{tube_num}_SolFn" + bc_node.append( + f"flux_{panel_node.name}_tube_{tube_num}", + type="FunctionNeumannBC", + variable="T_solid", + boundary=f"{panel_node.name}/{heat_tube_name}:rmax_0", + function=bc_func_name, + ) + + def heat_flux_data_to_moose_thm_data(self, func_name): + """ + Take tube HeatFluxBC.data and parse into PiecewiseMultilinear + compliant file for use with MOOSE + + Args: + func_name (String): name of function for this tube on this panel + + Returns: + axis_t (list): list of times in seconds + axis_x (list): list of x coordinates of data in MOOSE coords + axis_y (list): list of y coordinates of data in MOOSE coords + axis_z (list): list of z coordinates of data in MOOSE coords + data (list(list)): array of heat flux values in W/m^2 + [time, nx, ny, nz] + """ + bc = self.outer_bc + if bc is None: + print("NO FLUX BC ON TUBE???") + t = self.times * 3600 + # map srlife node index to tube coords + r = convert_mm_to_m(self.r) + thetas = np.linspace(np.pi, 0, self.nt) + h = convert_mm_to_m(self.h) + x_hats = np.array([r * np.cos(theta) for theta in thetas]) + y_hats = np.array([r * np.sin(theta) for theta in thetas]) + z_hats = np.linspace(0, h, self.nz) + # map tube coords to MOOSE THM coords + x = np.zeros_like(x_hats) + data = np.zeros([len(t), len(x), len(z_hats)]) + for iTime, time in enumerate(t): + hour_data = np.zeros([len(x), len(z_hats)]) + for iX, x_hat in enumerate(x_hats): + theta = thetas[iX] + y_hat = y_hats[iX] + for iZ, z_hat in enumerate(z_hats): + x[iX], _, z = self.map_tube_coords_to_moose_thm_coords( + x_hat, y_hat, z_hat + ) + hour_data[iX, iZ] = convert_Wmm2_to_Wm2( + bc.flux(time / 3600, theta, convert_m_to_mm(z)) + ) + # now I need to sort everything to be monotonically increasing + x_order = np.argsort(x) + x = x[x_order] + z_order = np.argsort(z_hats) + z_hats = z_hats[z_order] + hour_data = hour_data[x_order, :] + hour_data = hour_data[:, z_order] + data[iTime] = hour_data + return t, x, z_hats, data + + def map_tube_coords_to_moose_thm_coords(self, x_hat, y_hat, z_hat): + """ + Takes tube coordinates and maps them to MOOSE THM simulation coords + x = xc + R@x_hat + + Args: + x_hat (double): tube x coordinates + y_hat (double): tube y coordinates + z_hat (double): tube z coordinates + Returns: + x (double): moose thm x coords + y (double): moose thm y coords + z (double): moose thm z coords + """ + coord_tube = np.array([x_hat, y_hat, z_hat]) + # rotation matrix + c_yaw = np.cos(self.tube_yaw) + s_yaw = np.sin(self.tube_yaw) + R = np.array([[c_yaw, -s_yaw, 0], [s_yaw, c_yaw, 0], [0, 0, 1]]) + # transformation + coord = R @ coord_tube + np.array(self.center_point) + x = coord[0] + y = coord[1] + z = coord[2] + return x, y, z + + def map_nodes_moose_to_srlife(self, coords): + """ + Map nodal coordinates to nodal indices for result data coming from + MOOSE to be loaded into srlife Tube objects. + Tube stores results as [time, nr, nt, nz] structs + Data from MOOSE is [time, all_tube_nodes] + + Args: + coords (list[double]): nodal coordinates + + Returns: + r_ind, t_ind, z_ind (int): indices for current node from MOOSE + into srlife data structure + """ + # transform to tube coords + # translate via self.center_point + # rotate by self.tube_yaw + c_yaw = np.cos(self.tube_yaw) + s_yaw = np.sin(self.tube_yaw) + R = np.array([[c_yaw, s_yaw, 0], [-s_yaw, c_yaw, 0], [0, 0, 1]]) + coord_tube = R @ np.array(coords - self.center_point) + x = coord_tube[0] + y = coord_tube[1] + z = coord_tube[2] + # radial coords + r = np.sqrt(x**2 + y**2) + # correct with tube rotation + theta = np.pi + np.arctan2(y, x) + # helper q's + t = convert_mm_to_m(self.t) + ir = convert_mm_to_m(self.r - self.t) + h = convert_mm_to_m(self.h) + r_hat = r - ir + + # indices + r_ind = int(round((self.nr - 1) / t * r_hat)) + t_ind = int(round((self.nt) / (2 * np.pi) * theta)) + z_ind = int(round((self.nz - 1) / h * z)) + if t_ind == self.nt: + # have to correct for wrapped angle + t_ind = 0 + return r_ind, t_ind, z_ind + + def load_moose_thm_results(self, times, press_data, temp_data): + """ + Load solution data from MOOSE THM simulation into tube result data structures + + Args: + times (list): list of analysis times + press_data (np.array): Array nodal values for pressure along tube length + [time, node] + temp_data (list(np.array)): list of two arrays + 1 - nodal coordinates by index [node] + 2 - nodal values for temperature in tube with time [time, node] + """ + # take times from seconds to hours + times = np.array([it / 3600 for it in times]) + # set tube pressure BC + # downsample pressure data to a single value per time + # also convert from Pa to MPa + press = np.array( + [convert_Pa_to_MPa(np.average(press_step)) for press_step in press_data] + ) + + bc = PressureBC(times, press) + self.set_pressure_bc(bc) + # add temp results to tube + coords = temp_data[0] + temps = temp_data[1] + # loop through data + self.set_times(times) + tube_nod_temps = np.zeros([len(times), self.nr, self.nt, self.nz]) + repeats = 0 + coord_catch = set() + for iCoord, coord in enumerate(coords): + # for each node, get index mapping + r_ind, t_ind, z_ind = self.map_nodes_moose_to_srlife(coord) + coord_tup = (r_ind, t_ind, z_ind) + if coord_tup in coord_catch: + print("ERROR!!!!") + print(coord_tup) + repeats += 1 + print(f"{repeats} repeated coords!!!") + else: + coord_catch.add(coord_tup) + for iStep, temp in enumerate(temps): + # for each timestep, get temp of this node + tube_nod_temps[iStep, r_ind, t_ind, z_ind] = temp[iCoord] + self.add_results("temperature", tube_nod_temps) + def _vector_interpolate(base, data): """Interpolate as a vector @@ -1238,7 +2882,6 @@ def _generate_surface_mesh(self): """ ts = np.linspace(0, 2 * np.pi, self.nt + 1)[:-1] zs = np.linspace(0, self.h, self.nz) - return self.times, ts, zs def _generate_ifn(self, data): @@ -1257,7 +2900,6 @@ def _generate_ifn(self, data): bounds_error=False, fill_value=None, ) - return _make_ifn(base) @@ -1320,7 +2962,9 @@ def flux(self, t, theta, z): Returns: float: the flux value at this time and location """ - return self.ifn([t, theta, z]) + + res = self.ifn([t, theta, z]) + return res def save(self, fobj): """Save to an HDF5 file @@ -1509,7 +3153,7 @@ def __init__(self, radius, height, nz, fluid_T, film): self.fluid_T = fluid_T self.film = film - if fluid_T.shape != (nz,) or film.shape != (nz): + if fluid_T.shape != (nz,) or film.shape != (nz,): raise ValueError( "Film coefficient and fluid temperature data must have size (nz,)" ) diff --git a/srlife/thermal.py b/srlife/thermal.py index 14f21c6..3c4243a 100644 --- a/srlife/thermal.py +++ b/srlife/thermal.py @@ -86,6 +86,7 @@ def __init__(self, pset=solverparams.ParameterSet()): self.solid_params = pset.get_default("solid", solverparams.ParameterSet()) self.thermo_params = pset.get_default("fluid", solverparams.ParameterSet()) + self.conv_fluid_nodal_T = np.array([]) def solve_receiver( self, @@ -311,7 +312,9 @@ def solve_fluid(self, i, time, dt): ) # Solve for the fluid temperatures - nodal_temps = model.solve(time) + nodal_temps = model.solve(time, self.conv_fluid_nodal_T) + # save these values to use as initial guess in next solve + self.conv_fluid_nodal_T = nodal_temps # Update the stored fluid temperature and flow velocities flow_rates, tube_temperatures = model.recover_tube_results( diff --git a/srlife/thermohydraulics/flowpath.py b/srlife/thermohydraulics/flowpath.py index 40a17b0..1b11bf0 100644 --- a/srlife/thermohydraulics/flowpath.py +++ b/srlife/thermohydraulics/flowpath.py @@ -275,6 +275,11 @@ def __init__( self.miter = 100 self.verbose = verbose + # Save previous solution for use as initial guess + # to start, this is empty + self.T_prev_sol = np.array([]) + self.has_prev_conv_T = False + def add_panel(self, weights, ri, h, metal_temp, material): """ Construct and add the standard panel -> manifold link @@ -342,7 +347,7 @@ def _setup(self): self.dof_map.append(list(range(self.nvals, self.nvals + obj.size))) self.nvals += obj.size - def solve(self, t): + def solve(self, t, T_init_guess=np.array([])): """ Solve for the current fluid temperatures @@ -353,7 +358,10 @@ def solve(self, t): self.t = t # Significant decision... - T = np.zeros((self.nvals)) + if T_init_guess.size != 0: + T = T_init_guess + else: + T = np.zeros((self.nvals)) # Initial residual R, J = self.RJ(T) diff --git a/test/Rec_Gold.hdf5 b/test/Rec_Gold.hdf5 new file mode 100644 index 0000000..23ef1c7 Binary files /dev/null and b/test/Rec_Gold.hdf5 differ diff --git a/test/data_6_d20_hr6.csv b/test/data_6_d20_hr6.csv new file mode 100644 index 0000000..314e01e --- /dev/null +++ b/test/data_6_d20_hr6.csv @@ -0,0 +1,32 @@ +Flux intensity plot (units kW/m2) +Max. flux,845.368290,kW/m2 +Min. flux,5.230420,kW/m2 +Ave. flux,317.202681,kW/m2 + +,Receiver circumferential position (0=North; +East) [deg] --> +Receiver vertical position [m],180.0,165.6,151.2,136.8,122.4,108.0,93.6,79.2,64.8,50.4,36.0,21.6,7.2,-7.2,-21.6,-36.0,-50.4,-64.8,-79.2,-93.6,-108.0,-122.4,-136.8,-151.2,-165.6, +11.5,5.865428,11.283991,13.880661,16.842514,22.815152,32.645430,43.186667,58.565392,66.347012,72.230932,85.402586,97.840197,109.429092,121.305551,141.017264,153.904936,151.934404,162.649725,138.853229,117.681587,74.587489,43.714894,27.364659,20.800167,6.491756, +11.0,17.658028,29.897828,34.924294,38.709276,45.014505,56.841670,73.787590,98.083153,110.704619,123.382102,152.472569,186.575658,219.479074,254.866113,308.945140,333.658735,333.745597,368.698913,315.307864,273.352880,179.640739,110.583814,70.516157,54.587778,27.069226, +10.5,36.674288,53.725361,61.776878,66.756774,72.478764,85.964130,111.261308,144.960179,163.938502,186.376298,235.371592,297.652387,357.103579,420.466061,513.749207,548.871754,543.427509,563.915575,478.943700,407.648136,272.732278,164.844613,104.091638,79.118222,59.918282, +10.1,58.385827,75.340911,86.230176,94.034261,99.839781,115.841957,150.200820,191.523915,218.579634,253.160567,321.114763,408.338732,488.118858,568.696977,680.669811,724.857991,707.606809,678.833907,585.772046,491.729620,344.654885,204.458647,127.979380,92.955093,89.964197, +9.6,72.547265,89.155024,104.661356,116.229239,123.349889,143.254884,184.851984,230.298822,267.376752,315.788281,398.582772,501.378057,590.364240,668.949451,769.061696,824.309470,794.997917,698.053072,611.227515,502.933148,360.085565,211.882003,130.134056,94.298304,104.189166, +9.1,73.377691,92.286521,116.005702,132.032528,142.310761,166.880555,211.843447,257.966925,306.702048,369.585016,462.757405,570.680349,657.203287,715.064570,778.664687,840.188562,803.299265,678.705748,599.267354,483.599009,346.786118,204.969999,125.356205,92.589900,98.652784, +8.6,68.126975,88.699430,118.705590,140.683250,156.752226,185.942657,231.035091,276.232795,336.840817,412.816882,511.412373,614.727913,685.250575,718.887212,751.434871,810.145711,772.694578,659.723400,585.398800,463.408221,333.813909,199.803179,124.830549,91.094809,91.101483, +8.1,63.940598,85.630442,116.109772,143.114612,166.228001,199.362980,244.063480,289.067585,359.612452,444.521903,542.487606,635.177725,682.745331,704.719091,734.388058,785.658757,755.393345,657.349880,579.319679,449.218987,329.253533,188.562395,114.133972,84.836870,89.174000, +7.7,65.938220,87.002191,114.267323,142.468584,170.958670,206.581963,252.037516,299.190928,375.814454,463.690297,555.999563,637.906126,671.859390,694.371924,738.305565,794.258378,769.387549,671.946760,588.653729,467.603572,344.652994,194.437526,116.424957,90.995471,93.698518, +7.2,72.749707,90.784625,114.938310,142.247289,172.041149,208.226388,255.404607,307.232328,385.451697,470.630232,555.793486,631.760997,669.048445,694.993634,748.389640,819.531479,788.565498,685.087510,601.891633,496.902157,360.697599,208.416296,128.624520,104.597274,98.665569, +6.7,75.284866,92.306701,116.572879,142.513004,170.562028,206.046851,255.191811,312.701706,389.458540,468.575243,547.455033,622.168154,672.309036,702.227319,750.151506,827.207277,787.296141,686.564217,600.968183,498.867845,354.642957,214.084027,132.532529,105.996993,99.333239, +6.2,70.965130,90.513041,117.968740,141.772216,167.594225,202.900224,254.045631,316.943235,391.111456,463.535258,537.986543,612.409394,673.077829,709.334630,740.888804,807.798666,764.351904,675.427188,582.730503,469.032813,332.308408,202.373653,120.507249,91.521556,94.799816, +5.7,65.369504,87.770128,118.570521,140.517851,165.446535,201.413718,254.932867,321.221018,393.579647,461.522886,533.271927,608.270106,672.165140,712.171492,736.292361,788.051437,755.380651,664.216106,562.734177,445.375915,320.775324,191.000050,112.753438,80.974581,90.969674, +5.3,65.242545,87.630401,117.939154,140.059561,165.567066,202.408342,258.663648,324.952294,397.262377,465.003605,535.895485,612.184177,673.197282,709.606967,745.951442,798.345822,773.845989,670.894856,571.483626,466.254142,338.391464,198.446811,125.640660,89.312089,95.177955, +4.8,70.196015,90.199699,116.108580,140.232651,167.285845,205.241671,262.334833,324.229833,398.012697,471.886049,544.303563,620.755932,671.838768,701.648864,754.459887,818.378749,785.713573,689.098141,592.352752,496.439935,352.445625,210.427393,141.827345,99.332781,104.781586, +4.3,73.220712,91.635071,114.157709,140.770533,169.512385,207.393869,260.669920,313.829163,389.590378,474.782080,554.314978,629.917527,667.109611,692.568176,748.590553,814.474930,771.477045,691.755474,595.052072,497.307045,346.584656,207.445635,140.437135,96.522706,107.539481, +3.8,70.409957,90.916623,114.300921,142.161480,170.575333,205.693454,249.149957,291.363464,367.493169,464.898938,558.319548,638.446104,669.411450,691.472540,736.063023,792.200397,744.409295,675.053970,581.309402,471.844117,329.698648,201.432119,121.422787,84.903556,98.904763, +3.3,66.780681,89.950648,115.486836,142.499273,167.874759,197.647559,227.030377,258.131103,331.951006,436.722451,545.794823,639.946223,682.582719,702.179375,734.972898,783.185103,743.191115,655.591600,564.161077,455.173329,322.085878,201.880120,117.572959,83.055553,93.472703, +2.9,69.395461,90.794585,116.332867,138.255296,159.014698,182.116590,197.061057,218.908629,286.921045,390.521892,509.117976,620.450519,688.392962,715.721760,754.817587,808.516245,776.752821,660.333729,571.716366,467.782506,332.096122,208.290729,130.425582,90.604223,102.996816, +2.4,76.096444,93.912987,115.017278,127.997929,143.141059,159.085041,162.804739,178.222814,237.049681,330.724159,447.850780,569.836090,660.615014,710.692050,779.728021,845.368290,807.663824,685.291615,590.610631,485.669678,335.146460,201.751498,131.093149,90.737506,114.298249, +1.9,76.145988,92.622307,105.890143,111.298611,120.551093,129.631886,127.064487,138.460137,185.790386,263.167141,367.977154,488.420249,590.661192,662.443853,765.339580,833.682320,789.255551,705.651209,606.208101,503.061993,344.327190,207.680376,129.159426,88.261077,107.457543, +1.4,63.130895,80.433153,85.132597,87.938759,92.506930,96.425126,92.207886,101.001451,136.037633,194.051011,278.519378,384.219591,485.234091,559.577868,674.698926,729.640393,693.747841,679.421832,588.045764,492.662784,339.485696,217.604536,119.520737,85.949469,77.003034, +0.9,41.958724,59.259016,56.733468,59.984219,62.544511,63.997906,60.940287,67.682370,91.338342,130.400483,190.666579,270.068952,353.606763,411.888101,509.713656,544.559336,528.545285,554.225394,484.146205,404.930489,279.459182,187.582891,94.763647,74.438300,43.539273, +0.5,20.874981,34.346603,29.849904,33.548026,36.352097,37.111611,35.968006,40.788031,55.310048,78.728883,116.325708,165.596506,218.557759,252.014828,310.135687,325.431913,323.981447,358.292285,317.543371,267.717581,185.578683,129.770336,62.883900,51.794553,19.053695, +-0.0,6.838274,13.739533,11.878137,14.792917,17.984004,18.461297,18.648988,21.709624,29.880444,42.480334,62.727140,88.124764,112.684268,125.552273,146.843451,150.444623,149.552792,156.682002,138.905323,112.543353,75.687404,50.695634,23.663315,19.962505,5.230420, diff --git a/test/data_6_d20_hr7.csv b/test/data_6_d20_hr7.csv new file mode 100644 index 0000000..1beb8a5 --- /dev/null +++ b/test/data_6_d20_hr7.csv @@ -0,0 +1,32 @@ +Flux intensity plot (units kW/m2) +Max. flux,983.427384,kW/m2 +Min. flux,8.922926,kW/m2 +Ave. flux,396.461019,kW/m2 + +,Receiver circumferential position (0=North; +East) [deg] --> +Receiver vertical position [m],180.0,165.6,151.2,136.8,122.4,108.0,93.6,79.2,64.8,50.4,36.0,21.6,7.2,-7.2,-21.6,-36.0,-50.4,-64.8,-79.2,-93.6,-108.0,-122.4,-136.8,-151.2,-165.6, +11.5,8.922926,15.235459,15.685769,21.794866,31.948658,45.217232,54.469739,77.368883,83.637944,99.849242,119.121781,128.287710,151.895427,170.283334,164.679445,207.063765,211.601128,202.007188,170.970138,148.781458,95.484908,56.229191,34.401579,15.069059,11.714641, +11.0,29.528849,41.431577,39.220303,46.098350,61.591009,82.280134,98.675091,135.926759,149.747750,178.183071,220.106817,249.931875,311.494423,367.480888,359.621292,466.268497,489.703524,466.166995,392.716792,350.138107,232.988908,147.912204,96.486663,47.187907,39.618773, +10.5,60.938384,76.785420,71.634903,77.705316,99.077692,126.279086,152.031612,204.699260,230.785744,274.426236,343.744364,404.745647,505.332961,604.958286,603.398440,732.057587,746.701333,698.254957,587.610869,510.333073,340.939173,219.483638,155.357787,97.907739,83.279784, +10.1,87.588346,109.179247,103.311393,110.488520,137.335142,168.011649,204.633129,270.291947,311.761649,372.172976,464.858497,558.507991,677.500809,796.079043,820.935327,902.933964,898.445635,835.036405,712.882125,607.175918,410.498907,262.293785,183.012993,154.015067,122.510062, +9.6,96.287036,127.714499,127.011989,139.574885,169.231515,200.762637,248.367286,321.850411,379.415697,456.763456,563.392400,679.621901,791.182740,895.914871,958.689898,956.132834,916.681507,847.301445,731.708069,603.305992,409.601966,254.966704,172.228053,175.158471,131.236359, +9.1,90.143800,129.039677,140.242220,160.077581,190.294213,224.047020,281.724333,357.231639,429.676482,521.138343,629.886225,750.983334,835.693442,904.219990,983.427384,939.968280,883.380317,812.272659,711.540924,575.027630,396.005912,242.918301,157.120404,163.854282,117.177716, +8.6,83.518119,122.080194,142.697478,167.944854,200.294942,241.236328,306.475531,379.959550,464.564908,564.062626,662.694268,771.383697,832.136784,870.824521,936.427630,903.316042,848.285491,778.114304,694.105988,563.496718,397.992773,255.916526,177.988589,149.415893,111.110618, +8.1,84.844018,118.463194,138.658459,165.599334,202.250879,253.897090,322.203681,392.994969,485.475046,586.946228,669.439170,758.922648,816.250320,851.797831,891.460919,882.148651,841.094116,768.139258,676.133189,549.975923,387.048359,255.158521,190.101168,160.930912,126.711105, +7.7,88.606675,121.368241,134.598874,160.105816,201.085459,261.182535,326.877128,397.053303,491.820843,593.060860,664.836097,741.967334,808.589537,860.190711,899.164891,896.083232,868.887095,802.126067,693.041056,562.966183,398.573746,251.822365,183.396783,176.402568,137.802650, +7.2,88.710764,125.769425,133.170083,157.922124,200.522821,261.621018,321.133330,393.126654,485.161978,587.240577,661.752582,737.851675,811.417253,876.622836,935.762797,923.708858,892.948306,836.065051,723.803863,591.259520,419.344867,242.410757,166.677255,168.130270,128.310612, +6.7,86.065741,124.753428,134.301412,160.309443,200.867271,255.742796,310.372835,384.517840,471.875633,574.340288,662.382007,745.739866,820.111667,881.109089,947.579867,929.113098,888.759499,832.776927,739.211230,587.313809,433.472993,244.231068,177.563916,154.157871,116.069420, +6.2,86.851176,119.416419,137.085627,164.416711,200.750613,247.559469,301.901607,376.635029,460.551694,559.939780,662.615454,756.177377,828.487583,870.598640,922.294510,905.099617,854.867439,796.788237,718.561180,558.587341,412.199248,263.844891,199.314350,161.916323,122.734418, +5.7,89.944198,117.985859,139.770124,166.781516,199.735027,242.690870,301.664359,374.567335,457.908739,550.657663,660.314397,762.249796,832.090315,862.507878,898.352879,878.265216,838.176828,783.117760,690.109151,526.153235,387.453526,260.598841,194.389202,163.884132,130.553481, +5.3,89.997151,122.085845,140.100728,166.342130,198.567216,245.294549,311.857081,381.799322,466.912137,551.666036,656.961141,760.431916,830.252376,873.655748,915.045446,887.889543,865.355925,808.681142,694.228783,549.774008,393.866204,256.784323,169.279408,147.322551,121.243026, +4.8,89.016938,125.803005,137.797867,163.431607,198.236588,253.544994,328.770256,396.809319,484.309183,563.091723,654.056414,749.780361,821.193348,888.449898,944.178123,920.186307,893.232929,835.925933,720.447840,586.505016,411.281675,240.192916,164.269971,136.802472,106.931973, +4.3,89.585605,125.020490,135.508295,161.239002,199.411261,261.879500,343.723023,410.684704,501.928364,580.298203,655.398990,738.531490,807.999827,885.633338,941.081684,926.871260,891.109574,832.541985,723.872008,593.237626,413.601192,238.403508,190.497411,159.829692,113.365312, +3.8,89.184180,119.888233,136.709522,163.201599,201.690289,264.297254,345.662312,411.045793,507.533261,594.457950,664.975344,741.203071,801.621025,870.108223,910.592801,904.258621,858.150510,807.080057,699.575834,579.324819,400.026988,257.210962,196.531796,172.245066,126.789928, +3.3,88.066721,115.749899,141.537761,168.288468,202.118211,257.171147,328.525994,390.148364,490.999016,592.495575,676.400264,761.210249,810.063869,859.240873,897.831413,875.389555,835.998638,792.386018,673.780323,561.338860,393.633684,259.267434,179.446006,153.165240,123.084747, +2.9,89.637651,118.346971,146.327505,170.610268,196.608815,239.559905,294.101691,350.298785,449.723285,562.537575,669.778693,779.674951,830.921057,874.765262,930.612874,880.033816,849.270790,812.020584,679.540443,570.250074,411.013310,258.931136,157.142770,131.240173,109.023860, +2.4,95.625073,126.155695,143.126580,162.824216,181.968365,212.401973,249.688600,300.479925,389.101990,501.243356,625.314710,762.514919,840.506703,908.864416,970.330468,914.828105,882.660814,839.002695,703.170130,575.029215,411.689479,245.738549,174.539984,150.499486,109.352259, +1.9,98.987652,125.427274,123.163849,139.814752,157.063236,177.196743,201.199731,246.809825,316.670459,414.383581,537.723478,685.399883,796.632680,900.007702,949.631834,937.812431,917.069825,860.379968,731.356639,597.231557,420.878011,260.531246,197.470836,170.412804,119.088469, +1.4,88.497300,105.834032,88.435603,104.431550,122.911758,135.594408,151.351234,190.013396,238.658210,313.131211,418.000464,550.699542,676.507731,793.646209,826.210678,894.664366,902.993490,832.705332,720.638009,599.666548,414.077457,276.754368,179.212555,153.595975,114.637919, +0.9,63.003641,73.883265,51.880555,66.051805,84.542109,92.352880,103.037852,131.823296,162.457595,211.949464,287.643763,384.811562,496.902131,595.715128,619.813791,727.392819,755.892001,682.046833,597.090584,501.748379,337.279137,237.624803,118.780378,91.264270,84.510060, +0.5,32.955955,40.059096,24.578039,34.593556,49.616879,54.264599,61.490290,79.176869,97.472717,126.205968,172.268150,228.552909,303.243562,360.032313,377.496057,466.057780,497.897431,449.949907,397.557133,339.144048,226.251746,160.985653,61.623834,35.757786,44.002549, +-0.0,11.499806,14.719394,9.367245,14.635576,24.200211,26.711428,31.423467,39.928058,50.755003,65.652015,89.589039,115.675556,150.847196,170.674352,179.282804,209.431456,220.527056,196.411236,171.605219,144.060590,90.502142,60.620801,20.447250,9.534841,13.907637, diff --git a/test/test_interface.py b/test/test_interface.py new file mode 100644 index 0000000..4b186a6 --- /dev/null +++ b/test/test_interface.py @@ -0,0 +1,555 @@ +import unittest +import tempfile + +import numpy as np +from scipy.interpolate import RegularGridInterpolator +import multiprocess +import subprocess +import h5py + +from srlife import ( + receiver, + thermal, + structural, + system, + damage, + library, + managers, + interface, +) + +header = b"\n\n\n\n\n\n" + +fake_data = b"""-,180,135,90,45,0,-45,-90,-135, +0,1,1,1,1,1,1,1,1 +1,1,1,1,1,1,1,1,1 +2,1,1,1,1,1,1,1,1 +3,1,1,1,1,1,1,1,1 +""" + +nr_unit = 2 +nt_unit = 4 +nz_unit = 4 +tube_dict_unit = { + "od": 21.3, + "t": 1.25, + "h": 3, + "nr": nr_unit, + "nt": nt_unit, + "nz": nz_unit, + "spacing": 1, + "T0": 300, + "tube_k": "rigid", + "eta": 0.1, + "tube_mult": 58.5, + "ass_tube_per_panel": 2, +} + + +def calc_fd(Re, tube_eta, tube_Dh): + """ + Calculate Darcy friction factor from: + Zigrange and Sylvester 1985, A review of explicit friction factor equations. + J of Energy Resources Technology) + Equation 13 + + Args: + Re (double): Reynolds number + tube_eta (double): tube roughness + tube_Dh (double): tube hydraulic diameter + Returns: + fd (double): darcy friciton factor + """ + if Re < 4000: + # laminar flow + fd = 64 / Re + else: + fd = ( + 1 + / ( + -2 + * np.log10( + (tube_eta / tube_Dh) / 3.7 + - 5.02 / Re * np.log10((tube_eta / tube_Dh) / 3.7 + 13 / Re) + ) + ) + ** 2 + ) + return fd + + +class InterfaceUnitTests(unittest.TestCase): + """ + Run unit tests on some funcitons of interface.py + """ + + def test_file_not_found(self): + flux_filename = "no_file_existing.csv" + r = 10 + c = 17 + data_shape = [r, c] + z_offset_in = 0.0 + z_data, theta_data, flux_data = interface.read_month_day_hour_flux_file( + flux_filename, data_shape, z_offset_in + ) + print(theta_data) + self.assertTrue(np.array_equal(z_data, np.zeros([r, 1]))) + self.assertTrue(np.array_equal(theta_data, np.zeros([1, c]))) + self.assertTrue(np.array_equal(flux_data, np.zeros([r, c]))) + + def test_fake_file(self): + # Create a mock CSV file content to test the parsing + with tempfile.NamedTemporaryFile(delete=True) as f: + f.write(header) + f.write(fake_data) + f.seek(0) + z_data, theta_data, flux_data = interface.read_month_day_hour_flux_file( + f.name, [4, 8], 0 + ) + # spot check data + print(z_data) + self.assertEqual(z_data[2], 2000) + self.assertEqual(theta_data[-1], -180) + self.assertEqual(flux_data[2, 5], 1) + + def test_apply_tube_flux_bcs(self): + # test applying BCs via interpolator + times = np.array([0, 1]) + tube = interface.create_tube(tube_dict_unit, times, []) + tube_absorbance = 1.0 + tube_rec_theta = 12 + z_pts = np.linspace(10, 20, nz_unit) + theta_pts = np.linspace(15, -15, nt_unit) + data = np.array([[(2 * ix + 4 * iy) for ix in z_pts] for iy in theta_pts]) + interpFun = RegularGridInterpolator( + (z_pts, theta_pts), + data, + method="linear", + bounds_error=False, + fill_value=0.0, + ) + tube_rec_theta = 4 + interface.apply_tube_flux_bcs( + tube, tube_absorbance, tube_rec_theta, [interpFun] + ) + # spot check tube BCs + self.assertAlmostEqual(tube.outer_bc.data[0, 0, 2], interpFun([0.0, 3.0])) + self.assertEqual(tube.outer_bc.data[0, 3, 3], 0.0) + + def test_calc_fluid_velocity(self): + mfr = 10 + rho = 5 + Dh = 3 + vel = mfr / rho * (4.0 / (np.pi * Dh**2)) + self.assertAlmostEqual(interface.calc_fluid_velocity(mfr, rho, Dh), vel) + + def test_calc_reynolds_number(self): + mfr = 28 + rho = 5 + mu = 5 + Dh = 3 + vel = mfr / rho * (4.0 / (np.pi * Dh**2)) + Re = rho * vel * Dh / mu + self.assertAlmostEqual(interface.calc_reynolds_number(rho, vel, mu, Dh), Re) + + def test_calc_fd(self): + Re = [2000, 10000] + eta = 0.2 + Dh = 0.012 + for iRe in Re: + fd = interface.calc_fd(iRe, eta, Dh) + fd_test = calc_fd(iRe, eta, Dh) + self.assertAlmostEqual(fd, fd_test, places=12) + + def test_calc_friction_p_loss(self): + mfr = 0.001 + rho = 50 + mu = 0.004 + delta_l = 1 + p_loss = interface.calc_friction_p_loss(rho, mu, mfr, tube_dict_unit, delta_l) + test_val = 26.09267992590423 + self.assertAlmostEqual(p_loss, test_val) + + def test_calc_manifold_bend_loss(self): + mfr = 0.001 + rho = 1550 + mu = 0.004 + manifold_tube_dict = { + "od": 21.3, + "t": 1.25, + "bend_radius": 18, + "bends_per_panel": 4, + "eta": 0.01, + } + num_bends = 4 + + p_loss = interface.calc_manifold_bend_loss( + mfr, manifold_tube_dict, rho, mu, num_bends + ) + test_val = 0.12351409567895541 + self.assertAlmostEqual(p_loss, test_val) + + def test_create_receiver(self): + num_days = 1 + times = np.array([0, 1, 2]) + period = 2 + panel_k = "disconnect" + num_panels = 12 + results = [] + rec = interface.create_receiver( + tube_dict_unit, num_days, times, period, panel_k, num_panels, results + ) + self.assertEqual(num_panels, rec.npanels) + for panel in rec.panels.values(): + self.assertEqual(tube_dict_unit["ass_tube_per_panel"], panel.ntubes) + + def test_cycle_tube_pressure_bcs(self): + times = np.array([0, 1]) + tube_pressure = 12 + pressure = np.ones_like(times) * tube_pressure + pressure[0] = 0.0 + tube = interface.create_tube(tube_dict_unit, times, []) + tube_pressure_bc = receiver.PressureBC(times, pressure) + tube.set_pressure_bc(tube_pressure_bc) + + tubes_dict = {"0": tube} + num_cycles = 3 + cyclic_times = np.array([0, 1, 2, 3]) + interface.cycle_tube_pressure_bcs(tubes_dict, num_cycles, cyclic_times) + gold = np.array([0, 12, 12, 12]) + self.assertTrue(np.array_equal(gold, tube.pressure_bc.data)) + + def test_set_and_downsample_tube_temp_bcs(self): + times = np.array([0, 1, 2]) + cyclic_times = np.array([0, 1, 2]) + num_cycles = 1 + tube = interface.create_tube(tube_dict_unit, times, []) + tubes_dict = {"0": tube} + data = 245 * np.ones([tube.ntime, tube.nr, tube.nt, tube.nz]) + tube.add_results("temperature", data) + # Test 1: inlet T + interface.set_and_downsample_tube_temp_bcs( + tubes_dict, + True, + 245, + cyclic_times, + num_cycles, + "3d", + "", + ) + self.assertEqual(245, tube.T0) + analysis_type = "2d" + locs = ["max_T", "max_avg_T"] + for i_loc, loc in enumerate(locs): + tube = interface.create_tube(tube_dict_unit, times, []) + tubes_dict = {"0": tube} + data = (i_loc + 245) * np.ones([tube.ntime, tube.nr, tube.nt, tube.nz]) + tube.add_results("temperature", data) + interface.set_and_downsample_tube_temp_bcs( + tubes_dict, + False, + 245, + cyclic_times, + num_cycles, + analysis_type, + loc, + ) + self.assertTrue( + np.array_equal(tube.results["temperature"], data[:, :, :, 0]) + ) + + +class InterfaceRegressionTest(unittest.TestCase): + def test_simple_model(self): + multiprocess.set_start_method("spawn", force=True) + # Receiver filename to write to + rec_filename = "./test/Rec_Test" + # file extension used by srlife + hdf5_ext = ".hdf5" + + # GEOMETRY + # rec geom + num_panels = 2 # number of panels in receiver + rec_diam = 10000 # mm + rec_height = 12000 # mm + + # rec tube geom + tube_od = 21.3 # mm (Custom, Schedule 5) + tube_t = 1.24 # mm (Custom, Schedule 5) + tube_h = 12000 # mm (also rec height) + tube_spacing = 1 # mm + tube_eta = 0.01e-3 # tube roughness + tube_absorbance = 0.98 # absorbance of tube + z_offset = 0.25 # z offset for flux data + + # manifold tube geom + manifold_tube_od = 502.15 # mm + manifold_tube_t = 45.24 # mm + manifold_bend_radius = 2 * manifold_tube_od # used in p loss calc + manifold_tube_eta = tube_eta # roughness + manifold_bends_per_panel = 4 + # END GEOMETRY + + # ASSUMPTIONS + ass_tube_per_panel = 2 # number of tubes analyzed at each panel + outlet_p = 0.5 # MPa + + # tube assumption work + act_tube_per_panel = int( + np.pi * rec_diam / (num_panels * (tube_od + tube_spacing)) + ) + tube_multiplier = act_tube_per_panel / ass_tube_per_panel + print(f"Actual tubes per panel = {act_tube_per_panel}") + print(f"Tube Multiplier = {tube_multiplier}") + + # Structure assumptions + panel_k = "disconnect" # spring stiffness panel connect + tube_k = "rigid" # stiffness of tube connect + set_init_T_to_inlet_T = True + loc = "max_avg_T" # "max_T" + analysis_type = "2d" # "3d" + is_single_panel_analysis = True + single_panel_analysis_id = "0" + # END ASSUMPTIONS + + # BC/IC INFO + # flux data input and shape + flux_data_dir = "./test" + flux_data_shape = [25, 25] + tube_initial_temp = 300.0 # K + # END BC/IC INFO + + # ANALYSIS TIMES and CYCLE INFO + num_days = 1 + start_time = 6 + end_time = 7 + time_step = 1 # hour, can do nonInt + # stimes is index of hours to analyze, starting with first hour of flux BCs + stimes = np.arange(1, (end_time + 1) - start_time + time_step, time_step) + period = len(stimes) + # array of sample times of across number of days + times = np.tile(stimes, num_days) + # add time zero for initial conditions + times = np.append(np.array([0]), times) + # month and date from flux data + month = 6 + day = 20 + print("Times analyzed: \n", times) + num_cycles = 2 + # END ANALYSIS TIMES AND CYCLE INFO + + # TUBE DISC INFO + nz = int(interface.convert_mm_to_m(rec_height) * 0.5 + 1) + nt = 8 + nr = 2 + # END TUBE DISC INFO + + # SOLVERS # + # setup solver parameters + num_threads = 1 + # NOTE: using large tolerances to assure quick convergence + rtol = 1.0e-2 + atol = 1.0e-3 + params = interface.sample_parameters(num_threads, False, rtol, atol) + # thermal + thermal_solver = thermal.ThermohydraulicsThermalSolver(params["thermal"]) + # Structural solver + structural_solver = structural.PythonTubeSolver(params["structural"]) + # receiver system solver + system_solver = system.SpringSystemSolver(params["system"]) + # material damage model + damage_model = damage.TimeFractionInteractionDamage(params["damage"]) + # END SOLVERS # + + # MATERIALS + # fluid model + mat_fluid = library.load_thermal_fluid("32MgCl2-68KCl", "base") + # material models + st_mat_model = "elastic_creep" # "elastic_model" "base" + mat_thermal, mat_deformation, mat_damage = library.load_material( + "740H", "base", st_mat_model, "base" + ) + # END MATERIALS + + # structure output + save_struct_files_to_vtu = False + st_filename = "structure_" + st_mat_model + tube_filename = "tube_" + st_mat_model + + # FLOW SYSTEM SETUP + T_out_target = 720 # target outlet temperature, unit: C + # tolerance for mass flow optimization + # NOTE: making very large for test for single iter convergence + pct_err_outlet_temp = 100.25 + panel_flow_path = [["1"], ["0"]] + mass_flow_per_path = np.ones((len(panel_flow_path), len(times))) * 650 # kg/s + T_in_per_path = np.array([500, 500]) # Celcius + use_cycle_reset_heuristic = False + save_heat_to_vtu = False + # END FLOW SYSTEM SETUP + + # START ANALYSIS WORK + # rec tube data structure + tube_dict = { + "od": tube_od, + "t": tube_t, + "h": tube_h, + "nr": nr, + "nt": nt, + "nz": nz, + "spacing": tube_spacing, + "T0": tube_initial_temp, + "tube_k": tube_k, + "eta": tube_eta, + "tube_mult": tube_multiplier, + "ass_tube_per_panel": ass_tube_per_panel, + } + # manifold tube data structure + manifold_tube_dict = { + "od": manifold_tube_od, + "t": manifold_tube_t, + "bend_radius": manifold_bend_radius, + "bends_per_panel": manifold_bends_per_panel, + "eta": manifold_tube_eta, + } + + # STEP 1: create a receiver + rec = interface.create_receiver( + tube_dict, num_days, times, period, panel_k, num_panels, results=[] + ) + + # STEP 2: read and assign flux BCs to receiver tubes + # read flux files, make interpolator functions by hour and save + # NOTE: not doing multiple days right now + # NOTE: can we assume zData and thetaData are constant in time??? + flux_interpolators_by_hour = interface.get_flux_interpolators_from_data_files( + times, + start_time, + end_time, + month, + day, + flux_data_dir, + flux_data_shape, + z_offset, + ) + # Take flux interpolators and write flux BCs for each tube + # flux is sampled at tube centerline and applied smoothly + # across sunSideSurface via cosine + interface.calc_and_write_tube_flux_bcs( + rec, + ass_tube_per_panel, + num_panels, + flux_interpolators_by_hour, + tube_absorbance, + ) + + rec.save(rec_filename + hdf5_ext) + + # STEP 3: add flow paths to receiver + # add flow paths to receiver + interface.set_rec_flow_paths( + rec, panel_flow_path, mass_flow_per_path, T_in_per_path + ) + rec.save(rec_filename + hdf5_ext) + + # make solver + solver = managers.SolutionManager( + rec, + thermal_solver, + mat_thermal, + mat_fluid, + structural_solver, + mat_deformation, + mat_damage, + system_solver, + damage_model, + pset=params, + ) + + # STEP 4: optimize mass flow rate to match target outlet temp + # now optimize mass flow rate for each path + # this takes a long time, can skip if it has already been done + run_mass_flow_opt = True + # Heuristics + if use_cycle_reset_heuristic: + solver.add_heuristic(managers.CycleResetHeuristic()) + if run_mass_flow_opt: + interface.optimize_mass_flow_rate_per_path( + rec, + rec_filename, + T_out_target, + pct_err_outlet_temp, + solver, + save_heat_to_vtu, + ) + # Saving optimized results to base filename + rec.save(rec_filename + hdf5_ext) + else: + rec_filename = "./Rec_Test_ht_iter_2" + rec = receiver.Receiver.load(rec_filename + hdf5_ext) + + mass_flow = np.array([path["mass_flow"] for path in rec.flowpaths.values()]) + print(f"Mass flow: {mass_flow/3600}") + + # STEP 5: calc pressure loss from flow + # calculate pressure loss from + # pipe friction + # head loss + # manifold loss + flow_path_p_loss = interface.calc_p_loss_from_flows_temps( + rec, tube_dict, manifold_tube_dict, mat_fluid, outlet_p + ) + inlet_p = outlet_p + interface.convert_Pa_to_MPa(flow_path_p_loss) + interface.update_tube_pressure_bcs(rec, inlet_p, outlet_p) + rec.save(rec_filename + hdf5_ext) + + # STEP 6: solve structure and life receiver + struct_output_dict = { + "save_to_vtu": save_struct_files_to_vtu, + "st_filename": st_filename, + "tube_filename": tube_filename, + } + interface.run_struct_analysis( + rec_filename, + set_init_T_to_inlet_T, + loc, + analysis_type, + is_single_panel_analysis, + single_panel_analysis_id, + num_cycles, + solver, + struct_output_dict, + ) + # END ANALYSIS WORK + # do some checks on hdf5 files + file_gold = h5py.File("./test/Rec_Gold.hdf5", "r+") + file_test = h5py.File(rec_filename + hdf5_ext, "r+") + # flux bcs for a panel/tube + data_name = "/panels/1/tubes/0/outer_bc/data" + print("Flux BC") + self.assertTrue( + np.allclose( + file_test[data_name], file_gold[data_name], rtol=rtol, atol=atol + ) + ) + # fluid temps for two tubes + print("Fluid Temp") + data_name = "/panels/1/tubes/1/axial_results/fluid_temperature" + self.assertTrue( + np.allclose( + file_test[data_name], file_gold[data_name], rtol=rtol, atol=atol + ) + ) + # temps for a tube + print("Temp") + data_name = "/panels/1/tubes/0/results/temperature" + self.assertTrue( + np.allclose( + file_test[data_name], file_gold[data_name], rtol=rtol, atol=atol + ) + ) + # NOTE: in interface.run_struct_analysis we make this a + # single panel model, and don;t save it in this run + # so structure results are not available + # all else matches and run_struct_analysis finished, + # so call it cleared