Commit 61ea3e8d authored by William Gurecky's avatar William Gurecky
Browse files

update pod rom model to simplify use

parent b11511ec

dthydro_cfd/fv_domain.py

+49 −2
Original line number Diff line number Diff line
@@ -4,6 +4,7 @@ Defines domain and time integration 
import numpy as np

from typing import List, Union, Dict, Tuple

from copy import deepcopy

from collections import deque

#from abc import ABC

import abc

from dataclasses import dataclass
@@ -95,6 +96,9 @@ class PipeNetwork1D(object): 
        self.ode_w = np.array([])

        self.control_vars = {}

        self.transient_coeffs = {}

        # storage for past time derivatives of system variables

        self._df_dt_hist = deque(maxlen=kwargs.get("df_dt_hist_len", 10))

        self._df_dt_hist_dt = deque(maxlen=kwargs.get("df_dt_hist_len", 10))



    def set_control_var(self, control_var_name, val):

        """
@@ -174,6 +178,23 @@ class PipeNetwork1D(object): 
                term_name, u_var_name, flux_fn, b_fn,

                flux_fn_kwargs, b_fn_kwargs, scale)



    def add_vol_dt_source_term(

            self, term_name: str, u_var_name: List[str], source_fn: callable,

            source_fn_kwargs={}, scale=1.0):

        """

        Adds volumetric src term of the form:



        .. math::



            (1/\Delta t)\in_{t=0}^{\Delta t} f(U) \frac{dU(\tau)}{dt} d\tau



        Certain friction and dampening terms depend on

        the time driviative of the flow rate or pressure.

        """

        self.terms_dict[term_name] = SourceTermDuDt(

                term_name, u_var_name, source_fn,

                source_fn_kwargs=source_fn_kwargs, scale=scale)



    def add_vol_source_term(

            self, term_name: str, u_var_name: List[str], source_fn: callable,

            source_fn_kwargs={}, scale=1.0):
@@ -329,8 +350,15 @@ class PipeNetwork1D(object): 
        # compute contribution from flux, diffusion, and source terms

        f_t = np.zeros(my_fields.cc_fields.shape)

        for term_name, field_term  in iteritems(self.terms_dict):

            #assert isinstance(bc, FieldTerm1D)

            f_t += field_term.apply(my_fields, control_vars=self.control_vars)

            if isinstance(field_term, SourceTermDuDt):

                df_dt_hist_dt, df_dt_hist = self.get_df_dt_hist()

                f_t += field_term.apply(

                        my_fields,

                        df_dt_hist_dt, df_dt_hist,

                        control_vars=self.control_vars)

            else:

                f_t += field_term.apply(

                        my_fields, control_vars=self.control_vars)

        # Add volume terms from b cdot dF(u)/ds divergence theorem

        pass

        # eval coupled ODEs
@@ -353,6 +381,19 @@ class PipeNetwork1D(object): 
        self.fields.cc_fields[:] = self.old_cc_fields

        self.ode_w[:] = self.old_ode_w



    def save_df_dt(self, dt, df_dt):

        """

        Store current time drivative of all state vars

        """

        self._df_dt_hist.append(df_dt)

        self._df_dt_hist_dt.append(dt)



    def get_df_dt_hist(self):

        """

        Get the past time derivatives of all state variables

        """

        return np.asarray(self._df_dt_hist_dt), np.asarray(self._df_dt_hist)



    def step(self, delta_t, method="rk1", n_inners=12, relax=0.9, eps=1e-12):

        """

        Step the system forward by delta_t.
@@ -381,10 +422,12 @@ class PipeNetwork1D(object): 
        self.old_ode_w = deepcopy(self.ode_w)

        if method == "rk1":

            df_dt, dw_dt = self.sweep_system()

            self.save_df_dt(delta_t, df_dt)

            self.fields.cc_fields += df_dt * delta_t

            self.ode_w += dw_dt * delta_t

        elif method == "rk2" or method == "midpoint":

            df_dt_1, dw_dt_1 = self.sweep_system()

            self.save_df_dt(delta_t, df_dt_1)

            self.fields.cc_fields += df_dt_1 * (delta_t / 2.)

            self.ode_w += dw_dt_1 * (delta_t / 2.)

            df_dt_2, dw_dt_2 = self.sweep_system()
@@ -393,6 +436,7 @@ class PipeNetwork1D(object): 
            self.ode_w += dw_dt_2 * delta_t

        elif method == "heun":

            df_dt_1, dw_dt_1 = self.sweep_system()

            self.save_df_dt(delta_t, df_dt_1)

            self.fields.cc_fields += df_dt_1 * delta_t

            self.ode_w += dw_dt_1 * delta_t

            df_dt_2, dw_dt_2 = self.sweep_system()
@@ -401,6 +445,7 @@ class PipeNetwork1D(object): 
            self.ode_w += (delta_t / 2.) * (dw_dt_1 + dw_dt_2)

        elif method == "rk4":

            k1_f, k1_w = self.sweep_system()

            self.save_df_dt(delta_t, k1_f)

            self.fields.cc_fields += k1_f * (delta_t / 2.)

            self.ode_w += k1_w * (delta_t / 2.)

            k2_f, k2_w = self.sweep_system()
@@ -421,6 +466,7 @@ class PipeNetwork1D(object): 
        elif method == "impl_trap":

            # Implicit trapazoid rule

            df_dt_1, dw_dt_1 = self.sweep_system()

            self.save_df_dt(delta_t, df_dt_1)

            k1_f = df_dt_1 * delta_t

            k2_f = k1_f

            k1_w = dw_dt_1 * delta_t
@@ -440,6 +486,7 @@ class PipeNetwork1D(object): 
        elif method == "impl_midpoint":

            # Implicit midpoint rule

            df_dt_1, dw_dt_1 = self.sweep_system()

            self.save_df_dt(delta_t, df_dt_1)

            k1_f = df_dt_1 * delta_t

            k1_w = dw_dt_1 * delta_t

            # Inner fixed point iterations

dthydro_cfd/fv_fields.py

+51 −9
Original line number Diff line number Diff line
@@ -35,7 +35,7 @@ def _face_avg( 


def mesh_source_term(

        mesh: Mesh1D, global_cc_field, source_fn,

        source_fn_kwargs={}, face_average=False, **kwargs):

        source_fn_kwargs={}, face_average=False, du_dt=None, **kwargs):

    """

    Computes the source term over the mesh.

    """
@@ -52,19 +52,26 @@ def mesh_source_term( 
            u_cell_avg = (l_val + r_val) / 2.0

        else:

            u_cell_avg = global_cc_field[idx]

        source_field[idx] = source_fn(u_cell_avg, cell, **source_fn_kwargs)

        cell_du_dt = None

        if du_dt is not None:

            cell_du_dt = du_dt[idx]

        source_field[idx] = source_fn(

                u_cell_avg, cell, du_dt=cell_du_dt, **source_fn_kwargs)

    return source_field





def point_source_term(

        mesh: Mesh1D, idx: int, global_cc_field, source_fn,

        source_fn_kwargs={}, face_average=False, **kwargs):

        source_fn_kwargs={}, face_average=False, du_dt=None, **kwargs):

    """

    Computes the source term on one cell in the mesh

    """

    method = kwargs.get("method", "avg")

    source_field = np.zeros(global_cc_field.shape)

    cell = mesh.cell_dict[idx]

    cell_du_dt = None

    if du_dt is not None:

        cell_du_dt = du_dt[idx]

    if face_average:

        # trapazoid rule

        l_val, r_val = _face_avg(mesh, idx, global_cc_field, method=method)
@@ -73,7 +80,8 @@ def point_source_term( 
        source_field[idx] = source_face_avg_2

    else:

        u_cell_avg = global_cc_field[idx]

        cell_source_avg = source_fn(u_cell_avg, cell, **source_fn_kwargs)

        cell_source_avg = source_fn(

                u_cell_avg, cell, du_dt=cell_du_dt, **source_fn_kwargs)

        source_field[idx] = cell_source_avg

    return source_field
@@ -551,19 +559,21 @@ class MeshField1D(object): 
        return fluxes



    def cell_sources(

            self, field_name_list, source_fn, source_fn_kwargs={}, **kwargs):

            self, field_name_list, source_fn, source_fn_kwargs={},

            src_idx=-1, du_dt_avg=None, **kwargs):

        """

        Evaluate a source term on each cell in the mesh

        Evaluate a source term on each cell in the mesh or at a specified

        source idx location.

        """

        src_idx = kwargs.get("src_idx", -1)

        combo_field = self.select_cc_fields(field_name_list)

        if src_idx > -1:

            cell_sources = point_source_term(

                    self.mesh, src_idx, combo_field, source_fn,

                    source_fn_kwargs)

                    source_fn_kwargs, du_dt=du_dt_avg)

        else:

            cell_sources = mesh_source_term(

                self.mesh, combo_field, source_fn, source_fn_kwargs)

                self.mesh, combo_field, source_fn,

                    source_fn_kwargs, du_dt=du_dt_avg)

        return cell_sources



    def mesh_grad(self, field_name_list, b_fn, b_fn_kwargs={}, **kwargs):
@@ -770,6 +780,38 @@ class SourceTerm(FieldTerm): 
        return x





class SourceTermDuDt(SourceTerm):

    def __init__(self, *args, **kwargs):

        super().__init__(*args, **kwargs)



    def apply(self,

              field: MeshField1D,

              df_dt_hist_dt: np.ndarray,

              df_dt_hist: np.ndarray,

              control_vars={}, *args, **kwargs):

        x = np.zeros(field.cc_fields.shape)

        # Do nothing

        if len(df_dt_hist_dt) < 2:

            return x

        # compute past time derivative averages

        wgt_0 = kwargs.get("wgt_0", 0.1)

        # up-weight more recent times in the avg

        wgts = np.linspace(wgt_0, 1., df_dt_hist_dt.size) * df_dt_hist_dt

        hist_sum = np.sum(wgts.reshape((-1, 1, 1)) *

                          np.abs(df_dt_hist), axis=0)

        du_dt_avg = hist_sum / np.sum(wgts)

        assert du_dt_avg.shape == x.shape



        # compute the source term

        select_field_idxs = field.get_select_field_idxs(self.var_names)

        sources = field.cell_sources(

            self.var_names, self.source_fn,

            du_dt_avg=du_dt_avg,

            source_fn_kwargs=self.source_fn_kwargs | control_vars)

        x[:, select_field_idxs] = self.scale * sources

        return x





class FluxTerm(FieldTerm):

    def __init__(self, term_name: str, var_names: List[str], flux_fn: callable,

                 flux_fn_jac=None, flux_fn_kwargs={}, *args, **kwargs):

dthydro_cfd/rom_cfd/pod_cfd_dash_app.py

+4 −1
Original line number Diff line number Diff line
@@ -17,7 +17,6 @@ Author: William Gurecky. wll@ornl.gov. June 2023 
import numpy as np

import h5py

import time

from dthydro_cfd.rom_cfd.pod_rom import PODInterpModel

# plotly imports

import plotly.graph_objects as go

import plotly.express as px
@@ -25,6 +24,10 @@ import plotly.express as px 
import dash

from dash import dcc, html

from dash.dependencies import Input, Output

try:

    from dthydro_cfd.rom_cfd.pod_rom import PODInterpModel

except ImportError:

    from pod_rom import PODInterpModel





# Connect to h5 file

dthydro_cfd/rom_cfd/pod_cfd_model.py

+50 −14
Original line number Diff line number Diff line
@@ -17,7 +17,10 @@ import time 
import os

from scipy.interpolate import interp1d, RBFInterpolator

import h5py

try:

    from dthydro_cfd.rom_cfd.pod_rom import PODInterpModel

except ImportError:

    from pod_rom import PODInterpModel





def construct_pod_model(snapshot_h5file, data_name="Pressure (Pa)", n_pod_modes=5):
@@ -49,20 +52,51 @@ def construct_pod_model(snapshot_h5file, data_name="Pressure (Pa)", n_pod_modes= 
    return pod_model, xyz_points





class PodCfdModel(object):

    """

    Creates a Reduced Order Model (ROM) of the 3D CFD data

    of the turbine chamber.

    """

    def __init__(self, snapshot_h5file: str, data_name: str="Pressure (Pa)", n_pod_modes: int=7):

        pod_model, xyz_points = construct_pod_model(snapshot_h5file, data_name, n_pod_modes)

        self._pod_model = pod_model

        self._xyz_points = xyz_points



    def predict(self, vane_angle: float, min_vane_angle: float=4., max_vane_angle: float=9.):

        """

        Evaluate the POD ROM.



        Args:

            vane_angle:  turbine guide vane angle in degrees

            min_vane_angle:  min vane angle (model only works

                             between min and max CFD vane angle data)

            min_vane_angle:  max vane angle (model only works

                             between min and max CFD vane angle data)

        """

        va = np.clip(vane_angle, min_vane_angle, max_vane_angle)

        return self._pod_model.predict(va)



    @property

    def xyz_points(self):

        return self._xyz_points





def test_pod_cfd_model(use_plotly=True,

                       snapshot_h5file="cfd_pod_snapshots.h5",

                       n_pod_modes=5):

                       n_pod_modes=7):

    # create the reduced order model(s)

    pod_model, xyz_points = construct_pod_model(snapshot_h5file, data_name="Pressure (Pa)", n_pod_modes=n_pod_modes)

    pod_model = PodCfdModel(

            snapshot_h5file, data_name="Pressure (Pa)",

            n_pod_modes=n_pod_modes)

    xyz_points = pod_model.xyz_points



    pod_model_vx, _ = construct_pod_model(

    pod_model_vx = PodCfdModel(

            snapshot_h5file, data_name="Velocity[i] (m/s)")

    pod_model_vy, _ = construct_pod_model(

    pod_model_vy = PodCfdModel(

            snapshot_h5file, data_name="Velocity[j] (m/s)")

    pod_model_vz, _ = construct_pod_model(

    pod_model_vz = PodCfdModel(

            snapshot_h5file, data_name="Velocity[k] (m/s)")





    # evaluate the pod model at an vane angle not in training set

    ts = time.time()

    current_vane_angle = 8.02
@@ -125,16 +159,18 @@ def test_pod_cfd_model(use_plotly=True, 
        print("Time to make plot: %0.4e (sec)" % (te-ts))





def dynamic_pod_model(snapshot_h5file="cfd_pod_snapshots.h5", n_pod_modes=7):



def test_dynamic_pod_model(snapshot_h5file="cfd_pod_snapshots.h5", n_pod_modes=7):

    # create the reduced order model(s)

    pod_model, xyz_points = construct_pod_model(snapshot_h5file, data_name="Pressure (Pa)", n_pod_modes=n_pod_modes)

    pod_model = PodCfdModel(

            snapshot_h5file, data_name="Pressure (Pa)",

            n_pod_modes=n_pod_modes)

    xyz_points = pod_model.xyz_points



    pod_model_vx, _ = construct_pod_model(

    pod_model_vx = PodCfdModel(

            snapshot_h5file, data_name="Velocity[i] (m/s)")

    pod_model_vy, _ = construct_pod_model(

    pod_model_vy = PodCfdModel(

            snapshot_h5file, data_name="Velocity[j] (m/s)")

    pod_model_vz, _ = construct_pod_model(

    pod_model_vz = PodCfdModel(

            snapshot_h5file, data_name="Velocity[k] (m/s)")



    # update vane angle and store data in tmp hdf5 file
@@ -212,4 +248,4 @@ def dynamic_pod_model(snapshot_h5file="cfd_pod_snapshots.h5", n_pod_modes=7): 


if __name__ == "__main__":

    # test_pod_cfd_model()

    dynamic_pod_model()
 
    test_dynamic_pod_model()

examples/alder_inelastic_penstock.py

+83 −5
Original line number Diff line number Diff line
@@ -3,11 +3,17 @@ 
"""

import numpy as np

import h5py

import os

import time

from copy import deepcopy

import matplotlib.pyplot as plt

from dthydro_cfd.inelastic import *

from dthydro_cfd.constants import *

from dthydro_cfd.turbine_models import *

# cfd pod rom model of turbine chamber imports

from dthydro_cfd.rom_cfd.pod_cfd_model import PodCfdModel

pwd = os.path.dirname(os.path.abspath(__file__))





# globals (only used for plotting)

global flow_torque
@@ -23,6 +29,25 @@ def alder_inelastic(t_end=20., h_in=60., h_out=0.0, pid_control=True, *args, **k 
    fig_dir = "./ex_inelastic_alder_figs/"

    os.makedirs(fig_dir, exist_ok=True)



    # Storage for plots

    result_dict = {}



    # create the POD ROM model if given

    snapshot_h5file = kwargs.get(

            "pod_rom_snapshots",

            pwd + "/../dthydro_cfd/rom_cfd/cfd_pod_snapshots.h5")

    if snapshot_h5file:

        pod_model_p = PodCfdModel(

                snapshot_h5file, data_name="Pressure (Pa)", n_pod_modes=7)

        pod_xyz_points = pod_model_p.xyz_points

        result_dict['3d_xyz_points'] = pod_xyz_points

        pod_model_vx = PodCfdModel(

                snapshot_h5file, data_name="Velocity[i] (m/s)", n_pod_modes=7)

        pod_model_vy = PodCfdModel(

                snapshot_h5file, data_name="Velocity[j] (m/s)", n_pod_modes=7)

        pod_model_vz = PodCfdModel(

                snapshot_h5file, data_name="Velocity[k] (m/s)", n_pod_modes=7)



    # define pipe segments

    pipe_e = 8e-4  # pipe roughness

    pipe_d = 3.0   # pipe diamter
@@ -216,14 +241,60 @@ def alder_inelastic(t_end=20., h_in=60., h_out=0.0, pid_control=True, *args, **k 
        vis_w_setpoint.append(deepcopy(w_setpoint))

        vis_t.append(t)

        vis_a.append(a_turb)

        tq = inelastic_alder_model.control_vars["load_torque"]

        # tq = inelastic_alder_model.control_vars["load_torque"]

        vis_pow_shaft.append(w * flow_torque)

        vis_pow_tot.append(simple_turbine_work_tot(inelastic_alder_model.q, a_turb, w, **turbine_params))

        pow_tot = simple_turbine_work_tot(inelastic_alder_model.q, a_turb, w, **turbine_params)

        vis_pow_tot.append(pow_tot)

        vis_flow_torque.append(deepcopy(flow_torque))

        vis_load_torque.append(load_t)

        p_s = inelastic_alder_model.calc_pressure_profile()

        vis_p_s.append(p_s)

        #vis_load_torque.append(lt_t)



        # Results for external plots

        result_dict['time'] = t

        # Flow in (kg/s)

        result_dict['flow_rate'] = inelastic_alder_model.q

        # 1d array of Pressures in (Pa)

        result_dict['penstock_pressures'] = p_s

        # Location from inlet (m) where pressure is calculated

        result_dict['penstock_pressure_points'] = inelastic_alder_model.pipe_endpoints_s

        # Rotation rate in rad/s

        result_dict['turbine_rotation_rate'] = w

        # From PID controller output

        result_dict['turbine_rotation_rate_setpoint'] = w_setpoint

        # Vane angle in rad

        result_dict['turbine_vane_angle'] = a_turb

        # Shaft power in Watts

        result_dict['turbine_shaft_power'] = w * flow_torque

        # Shaft power plus friction loss (for eff calc)

        result_dict['turbine_tot_power'] = pow_tot

        # Torque from water

        result_dict['flow_torqe'] = flow_torque

        # Torque from electric generator to produce emf

        result_dict['load_torqe'] = load_t



        # for 3D plotting of turbine chamber

        if snapshot_h5file and i%10 == 0:

            ts = time.time()

            rom_vane_angle = np.rad2deg(a_turb)

            # NOTE: Current POD ROM snapshot data only valid between

            # 4.0 and 9.0 degrees of vane opening!

            pred_p = pod_model_p.predict(rom_vane_angle, min_vane_angle=4., max_vane_angle=9.)

            pred_vx = pod_model_vx.predict(rom_vane_angle, min_vane_angle=4., max_vane_angle=9.)

            pred_vy = pod_model_vy.predict(rom_vane_angle, min_vane_angle=4., max_vane_angle=9.)

            pred_vz = pod_model_vz.predict(rom_vane_angle, min_vane_angle=4., max_vane_angle=9.)

            pred_vmag = np.sqrt(pred_vx**2 + pred_vy**2 + pred_vz**2)

            print("ROM Vane angle (deg): ", rom_vane_angle)

            print("ROM Max, Min 3D P (Pa): ", max(pred_p), min(pred_p))

            print("ROM Max, Min 3D V (m/s): ", max(pred_vmag), min(pred_vmag))

            te = time.time()

            print("Time to eval POD ROM: %0.4e (sec)" % (te-ts))

            result_dict['3d_pred_p'] = pred_p

            result_dict['3d_pred_vx'] = pred_vx

            result_dict['3d_pred_vy'] = pred_vy

            result_dict['3d_pred_vz'] = pred_vz

        yield result_dict



    vis_q = np.asarray(vis_q)

    vis_w = np.asarray(vis_w)

    vis_w_setpoint = np.asarray(vis_w_setpoint)
@@ -329,12 +400,19 @@ def alder_inelastic(t_end=20., h_in=60., h_out=0.0, pid_control=True, *args, **k 
    plt.close()







if __name__ == "__main__":

    import argparse

    parser = argparse.ArgumentParser()

    parser.add_argument("-t_f", type=float, help="Final time of the simulation.", default=50.0)

    parser.add_argument("-h_in", type=float, help="Penstock inlet depth", default=57.28)  # +6.87 = 64.15

    parser.add_argument("-h_out", type=float, help="tailwater depth", default=0.0)  # 6.87

    parser.add_argument("--pod_rom_snapshots", type=str, help="file name of CFD POD snapshots", default="")

    args = parser.parse_args()

    alder_inelastic(t_end=args.t_f, h_in=args.h_in, h_out=args.h_out)
 
    for res_dict in alder_inelastic(

            t_end=args.t_f,

            h_in=args.h_in,

            h_out=args.h_out,

            pod_rom_snapshots=args.pod_rom_snapshots):

        # print(res_dict)

        # Could do realtime additional plotting here

        pass