Update Celeris for hydrostatic force coupling with OpenSees

Author: JustinBonus

modules/createEVENT/Celeris/Celeris.py

@@ -60,6 +60,24 @@
     print()  # noqa: T201
 
 
+# Check if bresenham is installed before importing bresenham
+try:
+    import bresenham
+except ImportError:
+    print('Bresenham is not installed. Please install it using "pip install bresenham".')
+    # noqa: T201
+    print()
+    subprocess.run([sys.executable, '-m', 'pip', 'install', 'bresenham'], check=False)  # noqa: S603
+    try:
+        import bresenham
+    except ImportError:
+        print('Bresenham installation failed. Please install it manually.')
+        # noqa: T201
+        sys.exit(1)
+    print('Bresenham is installed successfully.')
+    # noqa: T201
+    print()
+
 class FloorForces:  # noqa: D101
     def __init__(self, recorderID=-1):  # noqa: N803
         if recorderID < 0:
@@ -96,7 +114,7 @@ def __init__(self, recorderID=-1):  # noqa: N803
                             continue
                         # Assume there is no header in the file
                         # Assume recorder IDs are sequential, starting from 1
-                        if (j + 1) == recorderID:
+                        if (j) == recorderID:
                             # Strip away leading / trailing white-space,
                             # Delimit by regex to capture " ", \s, "  ", tabs, etc.
                             # Each value should be a number, rep. the force on recorder j at a time-step i
@@ -173,16 +191,14 @@ def addFloorForceToEvent(  # noqa: N802
     Add force (one component) time series and pattern in the event file
     Use of Wind is just a placeholder for now, since its more developed than Hydro
     """  # noqa: D205, D400
-    seriesName = '1'  # noqa: N806
-    patternName = '1'  # noqa: N806
     seriesName = 'WindForceSeries_' + str(floor) + direction  # noqa: N806
     patternName = 'WindForcePattern_' + str(floor) + direction  # noqa: N806
 
     pattern = {
         'name': patternName,
         'timeSeries': seriesName,
         'numSteps': len(force.X),
-        'dT': 0.01,
+        'dT': dt,
         'type': 'WindFloorLoad',
         'floor': str(floor),
         'story': str(floor),
@@ -196,8 +212,8 @@ def addFloorForceToEvent(  # noqa: N802
         'dof': directionToDof(direction),
         'floor': str(floor),
         'story': str(floor),
-        'dT': 0.01,
-        'dt': 0.01,
+        'dT': dt,
+        'dt': dt,
         'numSteps': len(force.X),
         'data': force.X,
     }
@@ -242,8 +258,8 @@ def writeEVENT(forces, eventFilePath='EVENT.json', floorsCount=1):  # noqa: N802
         'timeSeries': timeSeriesArray,
         'pressure': pressure,
         'numSteps': len(forces[0].X),
-        'dT': 1.0,
-        'dt': 1.0,
+        'dT': dt,
+        'dt': dt,
         'units': {'force': 'Newton', 'length': 'Meter', 'time': 'Sec'},
     }
 
@@ -341,13 +357,19 @@ def main():
         configFilename = event['configFile']  # noqa: N816
         bathymetryFilename = event['bathymetryFile']  # noqa: N816
         waveFilename = event['waveFile']  # noqa: N816
+        # Check for event['config']['dt'] in the event file
+        # and set dt to the value in the event file
+        # if not found, set dt to 0.01
+        dt = event['config']['dt'] if 'dt' in event['config'] else 0.01
 
     print('Running Celeris with script:', scriptName)  # noqa: T201
     print('Running Celeris with directory:', caseDirectory)  # noqa: T201
     print('Running Celeris with config file:', configFilename)  # noqa: T201
     print('Running Celeris with bathymetry:', bathymetryFilename)  # noqa: T201
     print('Running Celeris with waves:', waveFilename)  # noqa: T201
 
+    floorsCount = 1  # noqa: N816
+
     if arguments.getRV == True:  # noqa: E712
         print('RVs requested')  # noqa: T201
         # Read the number of floors
@@ -413,7 +435,6 @@ def main():
         )
 
         forces = []
-        floorsCount = 1  # noqa: N816
         for i in range(floorsCount):
             forces.append(FloorForces(recorderID=(i + 1)))
 

modules/createEVENT/Celeris/celeris/runner.py

@@ -5,7 +5,6 @@
 from celeris.solver import *  # noqa: F403
 from celeris.utils import *  # noqa: F403
 from taichi import tools  # For saving images
-
 # Done at top of translation unit for now, clean-up later - JB
 base_frame_dir = './plots/'  # ! Assumes we are in CelerisAi/ directory
 os.makedirs(  # noqa: PTH103
@@ -72,8 +71,8 @@ def Evolve_0(self):  # noqa: N802, D102
 
     def Evolve_Steps(self, step=0):  # noqa: C901, N802, D102
         i = step
-
-        self.solver.Pass1()
+        self.solver.update_step()
+        self.solver.Pass1(int(i))
 
         if self.solver.useSedTransModel:
             self.solver.Pass1_SedTrans()
@@ -125,7 +124,7 @@ def Evolve_Steps(self, step=0):  # noqa: C901, N802, D102
                     src=self.solver.NewState_Sed, dst=self.solver.State_Sed
                 )
 
-            self.solver.Pass1()
+            self.solver.Pass1(int(i))
 
             if self.solver.useSedTransModel:
                 self.solver.Pass1_SedTrans()
@@ -545,6 +544,8 @@ def Evolve_Display(  # noqa: C901, N802, D102
 
         start_time = time.time()
 
+        self.solver.overwrite_force()
+
         while window.running:
             # self.paint()
             # self.paint_new()
@@ -569,6 +570,8 @@ def Evolve_Display(  # noqa: C901, N802, D102
                 )  # using the Taichi tensors to render the image
             self.Evolve_Steps(i)
 
+            window.line(self.solver.force_sensor_begin_scaled, self.solver.force_sensor_end_scaled, radius=1, color=0x39FF14)
+            
             self.buffer_step = 1
             self.render_step = 10
             self.image_step = 100
@@ -635,7 +638,11 @@ def Evolve_Display(  # noqa: C901, N802, D102
                 # Improve the performance.The visualization is done only every 5 timesteps
                 window.show()
 
-            if i > self.maxsteps:
+            self.output_forces = True
+            if (self.output_forces):
+                self.solver.write_hydrostatic_force()
+                
+            if i > self.maxsteps:          
                 if self.saveimg:
                     print('Creating GIF...')  # noqa: T201
                     if frame_paths:  # Check if there are frames to create a GIF

modules/createEVENT/Celeris/celeris/solver.py

@@ -3,6 +3,8 @@
 import matplotlib.colors as clr
 import matplotlib.pyplot as plt
 import numpy as np
+from bresenham import bresenham
+import sys
 import taichi as ti
 import taichi.math as tm
 from celeris.utils import (
@@ -233,7 +235,35 @@ def __init__(  # noqa: C901, PLR0913, PLR0915
         self.dzdt_F_coef = dzdt_F_coef
         self.dzdt_I_coef = dzdt_I_coef
 
+        self.force_sensor_begin_scaled = [0.0, 0.0]
+        self.force_sensor_end_scaled = [1.0, 1.0]
+        self.force_sensor_begin_pixel = [0.0, 0.0]
+        self.force_sensor_end_pixel = [1.0, 1.0]
+        self.hydrostatic_forces = ti.field(dtype=ti.f32, shape=(self.maxsteps))
+        self.hydrostatic_forces_from_numpy = ti.field(dtype=ti.f32, shape=(self.maxsteps))
+        self.hydrostatic_forces.fill(-1.0)
+        self.hydrostatic_forces_numpy = np.zeros((self.maxsteps))
+        self.hydrostatic_forces_list = []
+        self.hydrostatic_force = ti.field(dtype=ti.f32, shape=(1))
+        # self.hydrostatic_forces_numpy = np.zeros((1, self.maxsteps))
+        self.step = 0
         if self.bc.celeris == True:  # noqa: E712
+            if checjson('force_sensor_begin', self.bc.configfile) == 1:
+                self.force_sensor_begin = self.bc.configfile['force_sensor_begin']
+                self.force_sensor_begin_scaled[0] = self.force_sensor_begin[0] / self.bc.configfile['WIDTH']
+                self.force_sensor_begin_scaled[1] = self.force_sensor_begin[1] / self.bc.configfile['HEIGHT']
+            else:
+                self.force_sensor_begin = [0.0, 0.0]
+                self.force_sensor_begin_scaled = [0.0, 0.0]
+                
+            if checjson('force_sensor_end', self.bc.configfile) == 1:
+                self.force_sensor_end = self.bc.configfile['force_sensor_end']
+                self.force_sensor_end_scaled[0] = self.force_sensor_end[0] / self.bc.configfile['WIDTH']
+                self.force_sensor_end_scaled[1] = self.force_sensor_end[1] / self.bc.configfile['HEIGHT']
+            else:
+                self.force_sensor_end = [self.bc.configfile['WIDTH'], self.bc.configfile['HEIGHT']]
+                self.force_sensor_end_scaled = [1.0, 1.0]
+                
             if checjson('whiteWaterDecayRate', self.bc.configfile) == 1:
                 self.whiteWaterDecayRate = float(
                     self.bc.configfile['whiteWaterDecayRate']
@@ -969,10 +999,23 @@ def tridiag_coeffs_Y(self):  # noqa: N802, D102
                 self.coefMaty[i, j] = ti.Vector([a, b, c, 0.0], self.precision)
 
     @ti.kernel
-    def Pass1(self):  # noqa: N802, D102
+    def Pass1(self, step: int):  # noqa: N802, D102
         # PASS 0 and Pass1 - edge value construction
         # using Generalized minmod limiter
         zro = ti.Vector([0.0, 0.0, 0.0, 0.0], self.precision)
+
+        # Bresenham's line algorithm
+        # bresenham_list = list(
+        #     bresenham(
+        #         self.force_sensor_begin[0],
+        #         self.force_sensor_begin[1],
+        #         self.force_sensor_end[0],
+        #         self.force_sensor_end[1],
+        #     )
+        # )
+        # print("Bresenham List: ", bresenham_list)
+        hydrostatic_force = 0.0
+        
         for i, j in self.State:
             # Compute the coordinates of the neighbors
             rightIdx = ti.min(i + 1, self.nx - 1)  # noqa: N806
@@ -1128,6 +1171,41 @@ def Pass1(self):  # noqa: N802, D102
                 [maxInundatedDepth, maxVelocityU, maxVelocityV, maxVelocityC],
                 self.precision,
             )
+                    
+            if (j >= self.force_sensor_begin[1] and j <= self.force_sensor_end[1]) and (i >= self.force_sensor_begin[0] and i <= self.force_sensor_end[0]):
+                # print("Calculating hydrostatic force at: ", i, j)
+                # Calculate the hydrostatic force
+                # Assuming a simple hydrostatic force calculation
+                # hydrostatic_force = 1/2 density * g * w * h
+                # where density is assumed to be 1 for simplicity
+                # and h is the water height at (i, j)
+                hsqr = float(self.H[i, j][0]) * float(self.H[i, j][0])
+                hydrostatic_force += 1000 * 0.5 * 1000.0 * float(self.g) * float(hsqr) * float(self.dy)
+
+        self.hydrostatic_forces[int(0)] = float(hydrostatic_force)
+        self.hydrostatic_forces[int(step)] = float(hydrostatic_force)
+
+
+    def overwrite_force(self):
+        # open forces.evt to overwrite
+        print("Overwriting forces.evt")
+        with open("forces.evt", "w") as f:
+            max_steps = 10000 # self.hydrostatic_forces.shape[0]
+            for i in range(max_steps):
+                if (i < max_steps - 1):
+                    f.write(f"{0.0} ")
+                else:
+                    f.write(f"{0.0}\n")
+                
+
+    def write_hydrostatic_force(self):
+        # open forces.evt to append        
+        with open("forces.evt", "a") as f:
+            f.write(f"{self.hydrostatic_forces[int(0)]} ")
+
+    def update_step(self):
+        self.step += 1
+
 
     @ti.kernel
     def Pass1_SedTrans(self):  # noqa: N802, D102