Merge branch 'alder_demo' into 'main' (4e1c0697) · Commits · dt-hydro / dt-hydro-cfd

dthydro_cfd/common_models.py

+19 −5

Original line number	Diff line number	Diff line
		@@ -98,7 +98,7 @@ def reynolds(rho: float, v: float, L: float, mu: float):
		return rho * v * L / mu


		def intake_loss_coeff(open_fraction: float, base_loss=0.0001, loss_scale=0.2):
		def valve_loss_coeff(open_fraction: float, base_loss=0.0001, loss_scale=0.2, exponent=1.9):
		"""
		Compute loss coefficeint of intake restriction
		that restricts flow from the reservoir to the penstock
		@@ -110,7 +110,7 @@ def intake_loss_coeff(open_fraction: float, base_loss=0.0001, loss_scale=0.2):
		loss_scale: Scale loss coefficeint values.
		"""
		assert 0.0 <= open_fraction <= 1.0
		loss_coeff = 1.0 / open_fraction ** 1.3
		loss_coeff = 1.0 / open_fraction ** exponent
		loss_coeff *= loss_scale
		loss_coeff += base_loss
		# 2k method
		@@ -153,8 +153,22 @@ def valve_h_loss(q, u_valve, pipe_a, **kwargs):
		t: time in s
		pipe_a: pipe flow area
		"""
		v_flow = q / (RHO_H2O * pipe_a)
		k_l = intake_loss_coeff(u_valve, **kwargs)
		v_flow = q / (kwargs.get("rho", RHO_H2O) * pipe_a)
		k_l = valve_loss_coeff(u_valve, **kwargs)
		valve_h_l = minor_head_loss(k_l, v_flow)
		return valve_h_l


		def inout_h_loss(q, pipe_a, **kwargs):
		"""
		Computes head loss over inlets or outlets.

		Args:
		q: mass flow rate
		pipe_a: pipe flow area
		"""
		v_flow = q / (kwargs.get("rho", RHO_H2O) * pipe_a)
		k_l = kwargs.get("k_l", 0.1)
		valve_h_l = minor_head_loss(k_l, v_flow)
		return valve_h_l

		@@ -233,7 +247,7 @@ def bernoulli_penstock_flow_steady(
		head_loss_major = pipe_head_loss

		# Minor head loss is losses to to restrictions and bends
		k_l_intake = intake_loss_coeff(intake_f)
		k_l_intake = valve_loss_coeff(intake_f)
		intake_loss = minor_head_loss(k_l_intake, v)
		discharge_head_loss = minor_head_loss(discharge_loss_coeff(), v)
		head_loss_minor = intake_loss + turbine_head_loss + discharge_head_loss

dthydro_cfd/constants.py

0 → 100644

+16 −0

Original line number	Diff line number	Diff line
		# === Constants ===

		# Density of water in [kg/m^3]
		RHO_H2O = 1000.0

		# dynamic viscosity of water [Pa*s]
		MU_H2O = 8.9e-4

		# Grav constant [m/s^2]
		G = 9.804

		# Youngs modulus of pipe wall
		K_PIPE = 2e11

		# Youngs modulus of water
		K_WATER = 2e9

dthydro_cfd/control_models.py

0 → 100644

+56 −0

Original line number	Diff line number	Diff line
		"""
		Constants implementation of PID controller for shaft speed
		and excitation systems.
		"""
		import numpy as np
		from collections import deque


		class PidControlSpeed(object):
		"""
		PID controller for shaft speed

		Args:
		k_p: potential err gain
		k_i: integral err gain
		k_d: derivative err gain
		"""
		def __init__(self, k_p=1.0, k_i=1.0, k_d=0.01, max_hist=4):
		self._k_p = k_p
		self._k_i = k_i
		self._k_d = k_d
		self.v_k = 0.0
		self.delta_w = deque(maxlen=max_hist)

		def update(self, w, w_setpoint, min_out=0, max_out=3.14, dt=1., **kwargs):
		"""
		PID control of vane angle to maintain set shaft speed.

		Args:
		w: shaft speed rad/s
		w_setpoint: shaft speed setpoint rad/s
		min_out: min guide vane open angle rad
		max_out: max guide vane open angle rad
		dt: time step size
		"""
		k_p = self._k_p
		k_i = self._k_i
		k_d = self._k_d
		e_w = w_setpoint - w
		self.delta_w.append(e_w)
		# delta error
		potential_err = self.delta_w[-1]
		# integral error
		self.v_k = self.v_k + k_i * dt * self.delta_w[-1]
		self.v_k = np.clip(self.v_k, min_out, max_out)
		# derivative error
		div_err = 0.0
		if len(self.delta_w) > 2:
		div_err = (self.delta_w[-1] - self.delta_w[-2]) / dt
		# Apply gains
		update_va = k_p * potential_err + \
		self.v_k + \
		k_d * div_err
		# print("pid p_e: %0.5e, v_k: %0.5e, d_e: %0.5e" %
		# (potential_err, self.v_k, div_err))
		return update_va

dthydro_cfd/flux_limiters.py

+6 −0

Original line number	Diff line number	Diff line
		@@ -2,14 +2,17 @@
		Library of flux limiters.
		"""
		import numpy as np
		from numba import jit


		@jit(nopython=True)
		def limit_minmod(r):
		""" The minmod flux limiter function """
		lim = np.maximum(0.0, np.minimum(1.0, r))
		return lim


		@jit(nopython=True)
		def limit_superbee(r):
		""" The superbee flux limiter function """
		return np.maximum(
		@@ -19,17 +22,20 @@ def limit_superbee(r):
		)


		@jit(nopython=True)
		def limit_van_albada_1(r):
		""" The van albada limiter """
		return (r 2.0 + r) / (r 2.0 + 1)


		@jit(nopython=True)
		def limit_r(u_l, u_o, u_r):
		""" Helper function to compute r parameter used in flux limiters """
		r = (u_o - u_l) / (u_r - u_o + 1e-20)
		return np.clip(r, -1e16, 1e16)


		@jit(nopython=True)
		def slope_limit_minmod(a, b):
		""" Alternate slope limit minmod formulation """
		return 0.5 * (np.sign(a) + np.sign(b)) * np.min((np.abs(a), np.abs(b)))

dthydro_cfd/fv_domain.py

+99 −14

Original line number	Diff line number	Diff line
		"""
		Defines domain and PDE
		Defines domain and time integration
		"""
		import numpy as np
		from typing import List, Union, Dict, Tuple
		@@ -12,7 +12,7 @@ from six import iteritems
		from scipy.optimize import fsolve, bisect, newton
		import matplotlib.pyplot as plt
		from scipy.integrate import odeint, solve_ivp
		from dthydro_cfd.common_models import reynolds, moody_colebrook, intake_loss_coeff, minor_head_loss
		from dthydro_cfd.common_models import reynolds, moody_colebrook, valve_loss_coeff, minor_head_loss
		from dthydro_cfd.fv_meshes import *
		from dthydro_cfd.fv_fields import *
		from dthydro_cfd.fv_domain import *
		@@ -63,6 +63,10 @@ class PipeSegment1D(object):
		self.fields = MeshField1D(self.mesh)
		self.scalars = {}

		@property
		def ds(self):
		return np.sqrt(np.sum((self.c_in - self.c_out) ** 2.0))

		@property
		def cell_idxs(self):
		return self.mesh.cell_idxs
		@@ -157,6 +161,19 @@ class PipeNetwork1D(object):
		term_name, u_var_name, flux_fn, flux_fn_jac=flux_fn_jac,
		flux_fn_kwargs=flux_fn_kwargs, scale=scale, coeff_fn=coeff_fn, **kwargs)

		def add_flux_body_term(
		self, term_name: str, u_var_name: List[str],
		flux_fn: callable, b_fn: callable,
		flux_fn_kwargs={}, b_fn_kwargs={}, scale=1.0):
		"""
		This term arises from applying the divergence theorem to

		\int B(U) \frac{d F(U)}{ds} dV = \int BF \cdot \hat n dS - \int F(\nabla B) dV
		"""
		self.terms_dict[term_name] = VolumeGradTerm(
		term_name, u_var_name, flux_fn, b_fn,
		flux_fn_kwargs, b_fn_kwargs, scale)

		def add_vol_source_term(
		self, term_name: str, u_var_name: List[str], source_fn: callable,
		source_fn_kwargs={}, scale=1.0):
		@@ -288,6 +305,7 @@ class PipeNetwork1D(object):
		print("WARNING: Reduce time step size or increase cell width.")
		return cr

		#@profile
		def sweep_system(self, fields=None):
		"""
		Computes the fluxes and vol source terms to obtain
		@@ -313,6 +331,8 @@ class PipeNetwork1D(object):
		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)
		# Add volume terms from b cdot dF(u)/ds divergence theorem
		pass
		# eval coupled ODEs
		for i, (ode_name, ode_term) in enumerate(iteritems(self.ode_dict)):
		dw_dt_s, s_ode = ode_term.apply(
		@@ -331,9 +351,9 @@ class PipeNetwork1D(object):
		adaptive time step methods.
		"""
		self.fields.cc_fields[:] = self.old_cc_fields
		self.ode_w = self.old_ode_w
		self.ode_w[:] = self.old_ode_w

		def step(self, delta_t, method="rk1"):
		def step(self, delta_t, method="rk1", n_inners=12, relax=0.9, eps=1e-12):
		"""
		Step the system forward by delta_t.
		Solves the following
		@@ -349,7 +369,7 @@ class PipeNetwork1D(object):
		---
		U_{t_1} \approx U_{t_0} + \delta t f_t(u_0, t_0)

		RK2
		RK2: midpoint rule
		---
		U_{t_{1/2}} \approx U_{t_0} + \delta t f_t(u_0, t_0) / 2
		U_{t_1} \approx U_{t_0} + \delta t f_t(u_{t_{1/2}}, t_0)
		@@ -360,18 +380,83 @@ class PipeNetwork1D(object):
		self.old_cc_fields[:] = self.fields.cc_fields
		self.old_ode_w = deepcopy(self.ode_w)
		if method == "rk1":
		f_t, dw_dt = self.sweep_system()
		self.fields.cc_fields += f_t * delta_t
		df_dt, dw_dt = self.sweep_system()
		self.fields.cc_fields += df_dt * delta_t
		self.ode_w += dw_dt * delta_t
		elif method == "rk2":
		f_t_1, dw_dt_1 = self.sweep_system()
		self.fields.cc_fields += f_t_1 * (delta_t / 2.)
		elif method == "rk2" or method == "midpoint":
		df_dt_1, dw_dt_1 = self.sweep_system()
		self.fields.cc_fields += df_dt_1 * (delta_t / 2.)
		self.ode_w += dw_dt_1 * (delta_t / 2.)
		f_t_2, dw_dt_2 = self.sweep_system()
		self.fields.cc_fields[:] = self.old_cc_fields
		self.ode_w[:] = self.old_ode_w
		self.fields.cc_fields += f_t_2 * delta_t
		df_dt_2, dw_dt_2 = self.sweep_system()
		self.rewind_step()
		self.fields.cc_fields += df_dt_2 * delta_t
		self.ode_w += dw_dt_2 * delta_t
		elif method == "heun":
		df_dt_1, dw_dt_1 = self.sweep_system()
		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()
		self.rewind_step()
		self.fields.cc_fields += (delta_t / 2.) * (df_dt_1 + df_dt_2)
		self.ode_w += (delta_t / 2.) * (dw_dt_1 + dw_dt_2)
		elif method == "rk4":
		k1_f, k1_w = self.sweep_system()
		self.fields.cc_fields += k1_f * (delta_t / 2.)
		self.ode_w += k1_w * (delta_t / 2.)
		k2_f, k2_w = self.sweep_system()
		#
		self.rewind_step()
		self.fields.cc_fields += k2_f * (delta_t / 2.)
		self.ode_w += k2_w * (delta_t / 2.)
		k3_f, k3_w = self.sweep_system()
		#
		self.rewind_step()
		self.fields.cc_fields += k3_f * delta_t
		self.ode_w += k3_w * delta_t
		k4_f, k4_w = self.sweep_system()
		#
		self.rewind_step()
		self.fields.cc_fields += (delta_t / 6.) * (k1_f + 2.k2_f + 2.k3_f + k4_f)
		self.ode_w += (delta_t / 6.) * (k1_w + 2.k2_w + 2.k3_w + k4_w)
		elif method == "impl_trap":
		# Implicit trapazoid rule
		df_dt_1, dw_dt_1 = self.sweep_system()
		k1_f = df_dt_1 * delta_t
		k2_f = k1_f
		k1_w = dw_dt_1 * delta_t
		k2_w = k1_w
		# Inner fixed point iterations
		for i in range(n_inners):
		self.rewind_step()
		self.fields.cc_fields += 0.5 * (k1_f + k2_f)
		self.ode_w += 0.5 * (k1_w + k2_w)
		df_dt_i, dw_dt_i = self.sweep_system()
		k_norm = np.linalg.norm(delta_t * df_dt_i - k2_f) / np.linalg.norm(k2_f)
		print("Implicit step %d coeff L2 norm=%0.5e " % (i, k_norm))
		k2_f = relax * (delta_t * df_dt_i) + (1.-relax) * k2_f
		k2_w = relax * (delta_t * dw_dt_i) + (1.-relax) * k2_w
		if k_norm < eps:
		break
		elif method == "impl_midpoint":
		# Implicit midpoint rule
		df_dt_1, dw_dt_1 = self.sweep_system()
		k1_f = df_dt_1 * delta_t
		k1_w = dw_dt_1 * delta_t
		# Inner fixed point iterations
		for i in range(n_inners):
		self.rewind_step()
		self.fields.cc_fields += 0.5 * k1_f
		self.ode_w += 0.5 * k1_w
		df_dt_i, dw_dt_i = self.sweep_system()
		k_norm = np.linalg.norm(delta_t * df_dt_i - k1_f) / np.linalg.norm(k1_f)
		print("Implicit step %d coeff L2 norm=%0.5e " % (i, k_norm))
		k1_f = relax * (delta_t * df_dt_i) + (1.-relax) * k1_f
		k1_w = relax * (delta_t * dw_dt_i) + (1.-relax) * k1_w
		if k_norm < eps:
		break
		self.rewind_step()
		self.fields.cc_fields += k1_f
		self.ode_w += k1_w
		else:
		raise NotImplementedError
		self.t += delta_t

Admin message