python jet actuator is more configurable, ghost node updates on bound…

…ary condition

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,19 +110,32 @@ def from_json(json_config: Dict[str, Any]):
 
         return Physics(mach, reynolds_friction, temp_ratio, visc_type, Tref, turb_inflow)
 
-class Jet():
-    def __init__(self, slot_start: int, slot_end: int):
-        self.slot_start = slot_start
-        self.slot_end = slot_end
+@unique
+class JetMethod(Enum):
+    none = "None"
+    constant = "Constant"
+    adaptive = "Adaptive"
 
-        self.has_slot = slot_start != -1
+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]):
-        slot_start = json_config["slot_start"]
-        slot_end = json_config["slot_end"]
+        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]
 
-        return Jet(slot_start, slot_end)
+        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, jet: Jet):
@@ -177,15 +191,6 @@ def slice_flowfield_array(self, array: np.ndarray) -> np.ndarray:
                 self.z_start():self.z_end()\
         ]
 
-    def global_slot_length(self):
-        # these values will be 1-indexed. For example: slot_start = 1, slot_end = 25 has
-        # 25 locations on the x axis that are blowing. However, end-start = 25 - 1 = 24 locations
-        # so we add an additional +1 to account for this
-        return self.jet.slot_end - self.jet.slot_start + 1
-
-    def local_slot_length(self):
-        return self.jet.slot_end - self.jet.slot_start + 1
-
     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

@@ -1,6 +1,9 @@
 import numpy as np
-from config import Config
+from config import Config, JetMethod, Jet
 import libstreams as streams
+from abc import ABC, abstractmethod
+from typing import Optional, Dict
+import utils
 
 # 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
@@ -32,20 +35,20 @@ def poly(self, amplitude: float) -> Polynomial:
         return Polynomial(a, b, c)
 
 class JetActuator():
-    def __init__(self, rank: int, config: Config):
-        # if x_start_slot == -1 then we do not have a matrix
-        # allocated on this MPI process
+    def __init__(self, rank: int, config: Config, slot_start: int, slot_end: int):
 
         self.rank = rank
         self.config = config
 
-        vertex_x = (config.jet.slot_start +  config.jet.slot_end) / 2
-        self.factory = PolynomialFactory(vertex_x, config.jet.slot_start, config.jet.slot_end)
+        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:
@@ -80,3 +83,47 @@ def set_amplitude(self, amplitude: float):
         streams.wrap_copy_blowing_bc_to_gpu()
         return None
 
+class AbstractActuator(ABC):
+    @abstractmethod
+    def step_actuator(self):
+        pass
+
+class NoActuation(AbstractActuator):
+    def __init__(self):
+        utils.hprint("skipping initialization of jet actuator")
+        pass
+
+    def step_actuator(self):
+       pass
+
+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)
+
+    def step_actuator(self):
+        self.actuator.set_amplitude(self.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.adaptive:
+        exit()
+    else:
+        exit()

python/main.py

@@ -78,6 +78,8 @@
 
 # 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)
@@ -104,12 +106,12 @@
 # Main solver loop, we start time stepping until we are done
 #
 
-actuator = jet_actuator.JetActuator(rank, config)
+actuator = jet_actuator.init_actuator(rank, config)
 
 time = 0
 for i in range(config.temporal.num_iter):
 
-    actuator.set_amplitude(1.0)
+    actuator.step_actuator()
 
     streams.wrap_step_solver()
 

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/bcblow.F90

@@ -51,6 +51,28 @@ subroutine bcblow(ilat)
                     ! w() is indexed differently on CPU and GPU
                     rho = w_gpu(i,j,k,1)
                     w_gpu(i,j,k,3) = rho * jet_velocity
+
+                    ! update the values at the ghost nodes. This is pretty similar to 
+                    ! the ghost nodes of bcwall, but instead we edit the y-velocity a little differently
+                    do l=1,ng
+                        rho  = w_gpu(i,1+l,k,1)
+                        uu   = w_gpu(i,1+l,k,2)/w_gpu(i,1+l,k,1)
+                        vv   = w_gpu(i,1+l,k,3)/w_gpu(i,1+l,k,1)
+                        ww   = w_gpu(i,1+l,k,4)/w_gpu(i,1+l,k,1)
+                        rhoe = w_gpu(i,1+l,k,5)
+
+                        qq   = 0.5_mykind*(uu*uu+vv*vv+ww*ww)
+                        pp   = gm1*(rhoe-rho*qq)
+                        tt   = pp/rho
+                        tt   = 2._mykind*twall-tt
+                        rho  = pp/tt
+
+                        w_gpu(i,1-l,k,1) =  rho
+                        w_gpu(i,1-l,k,2) = -rho*uu
+                        w_gpu(i,1-l,k,3) = -rho*(2.*jet_velocity  - vv)
+                        w_gpu(i,1-l,k,4) = -rho*ww
+                        w_gpu(i,1-l,k,5) = pp*gm+qq*rho
+                    enddo
 #else
                     jet_velocity = blowing_bc_slot_velocity(i-x_start_slot+1,k)
 

src/readinp.F90

@@ -90,9 +90,19 @@ subroutine readinp
  case (1) ! BL
   ibc(1) = ibc_inflow
   ibc(2) = 4 ! Extrapolation + non reflecting treatment
-  ibc(3) = 8 ! Wall + reflecting treatment
+  !ibc(3) = 8 
   ibc(4) = 4
   ibcnr(1) = 1
+
+  if (force_sbli_blowing_bc == 1) then
+      ! use blowing boundary condition
+      ibc(3) = blowing_sbli_boundary_condition
+  else
+      ! use default solver BC
+      ! Wall + reflecting treatment
+      ibc(3) = 8
+  endif
+
  case (2) ! SBLI
   ibc(1) = ibc_inflow
   ibc(2) = 4