Implement blowing slot boundary condition for arbitarary continuous x locations

python/config.py

@@ -1,6 +1,7 @@
 import json
 from typing import Any, Dict, Optional
 import numpy as np
+from enum import Enum, unique
 
 class Length():
     def __init__(self, lx: float, ly: float, lz: float):
@@ -109,13 +110,42 @@ def from_json(json_config: Dict[str, Any]):
 
         return Physics(mach, reynolds_friction, temp_ratio, visc_type, Tref, turb_inflow)
 
+@unique
+class JetMethod(Enum):
+    none = "None"
+    constant = "Constant"
+    sinusoidal = "Sinusoidal"
+    adaptive = "Adaptive"
+
+class Jet():
+    def __init__(self, jet_method: JetMethod, extra_json: Optional[Dict]):
+        self.jet_method = jet_method
+        self.extra_json = extra_json
+
+    @staticmethod
+    def from_json(json_config: Dict[str, Any]):
+        jet = json_config["blowing_bc"]
+
+        if jet == "None":
+            jet_method_str = "None"
+            extra_json = None
+        else:
+            jet_method_str = list(jet.keys())[0]
+            extra_json = jet[jet_method_str]
+
+        jet_method = JetMethod(jet_method_str)
+        print(jet_method)
+
+        return Jet(jet_method, extra_json)
+
 class Config():
-    def __init__(self, length: Length, grid: Grid, mpi: Mpi, temporal: Temporal, physics: Physics):
+    def __init__(self, length: Length, grid: Grid, mpi: Mpi, temporal: Temporal, physics: Physics, jet: Jet):
         self.length = length
         self.grid = grid
         self.mpi = mpi
         self.temporal = temporal
         self.physics = physics
+        self.jet = jet
 
     @staticmethod
     def from_json(json_config: Dict[str, Any]):
@@ -124,8 +154,9 @@ def from_json(json_config: Dict[str, Any]):
         mpi = Mpi.from_json(json_config)
         temporal = Temporal.from_json(json_config)
         physics = Physics.from_json(json_config)
+        jet = Jet.from_json(json_config)
 
-        return Config(length, grid, mpi, temporal, physics)
+        return Config(length, grid, mpi, temporal, physics, jet)
 
     def x_start(self) -> int:
         return self.grid.ng
@@ -160,3 +191,7 @@ def slice_flowfield_array(self, array: np.ndarray) -> np.ndarray:
                 self.y_start():self.y_end(), \
                 self.z_start():self.z_end()\
         ]
+
+    def local_to_global_x(self, x: int, rank: int):
+        previous_mpi_x = rank * self.nx_mpi();
+        return previous_mpi_x + x

python/jet_actuator.py

@@ -0,0 +1,171 @@
+import numpy as np
+from config import Config, JetMethod, Jet
+import libstreams as streams
+from abc import ABC, abstractmethod
+from typing import Optional, Dict
+import utils
+import math
+
+# the equation of the polynomial for the jet actuation in coordinates local to the jet
+# actuator. This means that the jet actuator starts at x = 0 and ends at x = slot_end
+# 
+# this must be recomputed at every amplitude change
+class Polynomial():
+    def __init__(self, a: float, b: float, c: float):
+        self.a = a
+        self.b = b
+        self.c = c
+
+    def evaluate(self, x_idx: int) -> float:
+        return self.a * x_idx**2  + \
+                self.b * x_idx + \
+                self.c
+
+# helper class to recompute the polynomial of the jet actuator
+class PolynomialFactory():
+    def __init__(self, vertex_x: float, slot_start: float, slot_end: float):
+        self.slot_start = slot_start
+        self.slot_end = slot_end
+        self.vertex_x  = vertex_x
+
+    def poly(self, amplitude: float) -> Polynomial:
+        # see streams section of lab documentation for the derivation of this
+        #
+        # in essence, it is solving for the coefficients a,b,c of the polynomial
+        # y = ax^2 + bx +c 
+        # using the fact that
+        # y(jet start) = 0
+        # y(jet end) = 0
+        # y((jet start + jet_end) / 2) = amplitude
+        a = amplitude/(self.vertex_x**2 - self.vertex_x*self.slot_end- (self.vertex_x - self.slot_end)*self.slot_start)
+        b = -(amplitude*self.slot_end+ amplitude*self.slot_start)/(self.vertex_x**2 - self.vertex_x*self.slot_end- (self.vertex_x - self.slot_end)*self.slot_start)
+        c = amplitude*self.slot_end*self.slot_start/(self.vertex_x**2 - self.vertex_x*self.slot_end- (self.vertex_x - self.slot_end)*self.slot_start)
+
+        return Polynomial(a, b, c)
+
+class JetActuator():
+    def __init__(self, rank: int, config: Config, slot_start: int, slot_end: int):
+
+        self.rank = rank
+        self.config = config
+
+        vertex_x = (slot_start +  slot_end) / 2
+        self.factory = PolynomialFactory(vertex_x, slot_start, slot_end)
+
+        self.local_slot_start_x = streams.mod_streams.x_start_slot
+        self.local_slot_nx = streams.mod_streams.nx_slot
+        self.local_slot_nz = streams.mod_streams.nz_slot
+
+        # if x_start_slot == -1 then we do not have a matrix
+        # allocated on this MPI process
+        self.has_slot = streams.mod_streams.x_start_slot != -1
+
+        if self.has_slot:
+            #print(self.local_slot_nx, self.local_slot_nz, self.local_slot_start_x)
+            self.bc_velocity = streams.mod_streams.blowing_bc_slot_velocity[:, :]
+
+    def set_amplitude(self, amplitude: float):
+        # WARNING: copying to GPU and copying to CPU must happen on ALL mpi procs
+        # if we have a matrix, we can adjust the velocity of the blowing actuator
+        streams.wrap_copy_blowing_bc_to_cpu()
+
+        if not self.has_slot:
+            streams.wrap_copy_blowing_bc_to_gpu()
+            return None
+
+        # calculate the equation of the polynomial that the velocity of the jet
+        # actuator will have with this amplitude
+        poly = self.factory.poly(amplitude)
+
+        for idx in range(self.local_slot_nx):
+            local_x = self.local_slot_start_x + idx
+            global_x = self.config.local_to_global_x(local_x, self.rank)
+
+            velo =  poly.evaluate(global_x)
+            #print(f"velocity at idx {idx} / local x {local_x} / global x {global_x} ::: {velo}")
+
+            self.bc_velocity[idx, 0:self.local_slot_nz] = velo
+
+        #print(f"array shape for BC", self.bc_velocity.shape)
+
+        # finally, copy everything back to the GPU
+        streams.wrap_copy_blowing_bc_to_gpu()
+        return None
+
+class AbstractActuator(ABC):
+    @abstractmethod
+    # returns the amplitude of the jet that was used
+    def step_actuator(self, time: float) -> float:
+        pass
+
+class NoActuation(AbstractActuator):
+    def __init__(self):
+        utils.hprint("skipping initialization of jet actuator")
+        pass
+
+    # returns the amplitude of the jet that was used
+    def step_actuator(self, _:float) -> float:
+        return 0.
+
+class ConstantActuator(AbstractActuator):
+    def __init__(self, amplitude: float, slot_start: int, slot_end: int, rank: int, config: Config):
+        utils.hprint("initializing a constant velocity actuator")
+
+        self.slot_start = slot_start
+        self.slot_end = slot_end
+        self.amplitude = amplitude
+
+        self.actuator = JetActuator(rank, config, slot_start, slot_end)
+
+    # returns the amplitude of the jet that was used
+    def step_actuator(self, _: float) -> float:
+        self.actuator.set_amplitude(self.amplitude)
+        return self.amplitude
+
+class SinusoidalActuator(AbstractActuator):
+    def __init__(self, amplitude: float, slot_start: int, slot_end: int, rank: int, config: Config, angular_frequency:float ):
+        utils.hprint("initializing a constant velocity actuator")
+
+        self.slot_start = slot_start
+        self.slot_end = slot_end
+        self.amplitude = amplitude
+
+        self.actuator = JetActuator(rank, config, slot_start, slot_end)
+        self.angular_frequency = angular_frequency
+
+    # returns the amplitude of the jet that was used
+    def step_actuator(self, time: float) -> float:
+        adjusted_amplitude = math.sin(self.angular_frequency * time)
+
+        self.actuator.set_amplitude(adjusted_amplitude)
+
+        return adjusted_amplitude
+
+def init_actuator(rank: int, config: Config) -> AbstractActuator:
+    jet_config = config.jet
+
+    if jet_config.jet_method == JetMethod.none:
+        return NoActuation()
+    elif jet_config.jet_method == JetMethod.constant:
+        print(jet_config.extra_json)
+        # these should be guaranteed to exist in the additional json information
+        # so we can essentially ignore the errors that we have here
+        slot_start = jet_config.extra_json["slot_start"]
+        slot_end = jet_config.extra_json["slot_end"]
+        amplitude = jet_config.extra_json["amplitude"]
+
+        return ConstantActuator(amplitude, slot_start, slot_end, rank, config);
+    elif jet_config.jet_method == JetMethod.sinusoidal:
+        print(jet_config.extra_json)
+        # these should be guaranteed to exist in the additional json information
+        # so we can essentially ignore the errors that we have here
+        slot_start = jet_config.extra_json["slot_start"]
+        slot_end = jet_config.extra_json["slot_end"]
+        amplitude = jet_config.extra_json["amplitude"]
+        angular_frequency = jet_config.extra_json["angular_frequency"]
+
+        return SinusoidalActuator(amplitude, slot_start, slot_end, rank, config, angular_frequency);
+    elif jet_config.jet_method == JetMethod.adaptive:
+        exit()
+    else:
+        exit()

python/main.py

@@ -28,6 +28,7 @@
 import numpy as np
 from config import Config
 import utils
+import jet_actuator
 
 #
 # Load in config, initialize
@@ -43,7 +44,10 @@
 span_average = np.zeros([5, config.nx_mpi(), config.ny_mpi()], dtype=np.float64)
 temp_field = np.zeros((config.nx_mpi(), config.ny_mpi(), config.nz_mpi()), dtype=np.float64)
 dt_array = np.zeros(1)
+amplitude_array = np.zeros(1)
 time_array = np.zeros(1)
+dissipation_rate_array = np.zeros(1)
+energy_array = np.zeros(1)
 
 #
 # execute streams setup routines
@@ -73,14 +77,19 @@
 velocity_dset = io_utils.VectorField3D(flowfields, [5, *grid_shape], flowfield_writes, "velocity", rank)
 flowfield_time_dset = io_utils.Scalar1D(flowfields, [1], flowfield_writes, "time", rank)
 
-# span average files 
+# span average files
 numwrites = int(math.ceil(config.temporal.num_iter / config.temporal.span_average_io_steps))
+
+# this is rho, u, v, w, E (already normalized from the rho u, rho v... values from streams)
 span_average_dset = io_utils.VectorFieldXY2D(span_averages, [5, *span_average_shape], numwrites, "span_average", rank)
 shear_stress_dset = io_utils.ScalarFieldX1D(span_averages, [config.grid.nx], numwrites, "shear_stress", rank)
 span_average_time_dset = io_utils.Scalar0D(span_averages, [1], numwrites, "time", rank)
+dissipation_rate_dset = io_utils.Scalar0D(span_averages, [1], numwrites, "dissipation_rate", rank)
+energy_dset = io_utils.Scalar0D(span_averages, [1], numwrites, "energy", rank)
 
 # trajectories files
 dt_dset = io_utils.Scalar0D(trajectories, [1], config.temporal.num_iter, "dt", rank)
+amplitude_dset = io_utils.Scalar0D(trajectories, [1], config.temporal.num_iter, "jet_amplitude", rank)
 
 # mesh datasets
 x_mesh_dset = io_utils.Scalar1DX(mesh_h5, [config.grid.nx], 1, "x_grid", rank)
@@ -99,8 +108,13 @@
 # Main solver loop, we start time stepping until we are done
 #
 
+actuator = jet_actuator.init_actuator(rank, config)
+
 time = 0
 for i in range(config.temporal.num_iter):
+
+    amplitude = actuator.step_actuator(time)
+
     streams.wrap_step_solver()
 
     time += streams.mod_streams.dtglobal
@@ -121,10 +135,27 @@
         # write the time at which this data was collected
         span_average_time_dset.write_array(time_array)
 
+        # calculate dissipation rate on GPU and store the result
+        streams.wrap_dissipation_calculation()
+        dissipation_rate_array[:] = streams.mod_streams.dissipation_rate
+        dissipation_rate_dset.write_array(dissipation_rate_array)
+        utils.hprint(f"dissipation is {dissipation_rate_array[0]}")
+
+        # calculate energy on GPU and store the result
+        streams.wrap_energy_calculation()
+        energy_array[:] = streams.mod_streams.energy
+        energy_dset.write_array(energy_array)
+
+        utils.hprint(f"energy is {energy_array[0]}")
+
     # save dt information for every step
     dt_array[:] = streams.mod_streams.dtglobal
     dt_dset.write_array(dt_array)
 
+    # save amplitude at every step
+    amplitude_array[:] = amplitude
+    amplitude_dset.write_array(amplitude_array)
+
     if not (config.temporal.full_flowfield_io_steps is None):
         if (i % config.temporal.full_flowfield_io_steps) == 0:
             utils.hprint("writing flowfield")

python/utils.py

@@ -12,9 +12,12 @@ def hprint(*args):
 def calculate_span_averages(config: Config, span_average: np.ndarray, temp_field: np.ndarray, streams_data: np.ndarray):
     # for rhou, rhov, rhow, and rhoE
     for data_idx in [1,2,3,4]:
+        # divide rho * VALUE by rho so we exclusively have VALUE
         temp_field[:] = np.divide(streams_data[data_idx, :, :, :], streams_data[0, :, :, :], out = temp_field[:])
 
+        # sum the VALUE across the z direction
         span_average[data_idx, :, :] = np.sum(temp_field, axis=2, out=span_average[data_idx, :, :])
+        # divide the sum by the number of points in the z direction to compute the average
         span_average[data_idx, :, :] /= config.grid.nz
 
     span_average[0, :, :] = np.sum(streams_data[0, :, :, :], axis=2, out=span_average[0, :, :])

src/Makefile

@@ -54,7 +54,7 @@ OBJ_FILES = alloc.o bcdf.o bcextr.o bcfree.o bc.o bcrecyc.o bcrelax.o bcshk.o bc
     setup.o solver.o startmpi.o stats.o step.o target_reystress.o tbforce.o updateghost.o utility.o utyibm.o \
     visflx.o writedf.o writefield.o writefieldvtk.o writegridplot3d.o writerst.o \
     writestatbl.o writestatchann.o writestat.o writestatzbl.o write_wallpressure.o visflx_stag.o \
-	write_probe_data.o bcblow.o
+	write_probe_data.o bcblow.o dissipation.o
 
 PY_API_OBJS = min_api.o
 
@@ -67,7 +67,7 @@ STATIC_FILES = alloc.F90 bcdf.F90 bcextr.F90 bcfree.F90 bc.F90 bcrecyc.F90 bcrel
     setup.F90 solver.F90 startmpi.F90 stats.F90 step.F90 target_reystress.F90 tbforce.F90 updateghost.F90 utility.F90 utyibm.F90 \
     visflx.F90 writedf.F90 writefield.F90 writefieldvtk.F90 writegridplot3d.F90 writerst.F90 \
     writestatbl.F90 writestatchann.F90 writestat.F90 writestatzbl.F90 write_wallpressure.F90 visflx_stag.F90 \
-	write_probe_data.F90 bcblow.F90 
+	write_probe_data.F90 bcblow.F90 dissipation.F90
 
 STATIC_MODS = mod_streams.F90 mod_sys.F90 euler.F90
 

src/alloc.F90

@@ -215,6 +215,11 @@ subroutine allocate_vars()
  allocate(wrecycav_gpu(ng,ny,nv))
 
  allocate(tauw_x(1:nx))
+
+ allocate(fdm_y_stencil_gpu(5, ny))
+ allocate(fdm_y_stencil(5, ny))
+ allocate(fdm_individual_stencil(0:1, 0:4, 0:4))
+ allocate(fdm_grid_points(0:4))
 !
 endsubroutine allocate_vars