Merge pull request #396 from JustinBonus/master

Update Celeris

Author: fmckenna

modules/createEVENT/Celeris/Celeris.py

@@ -61,7 +61,7 @@
 
 class FloorForces:  # noqa: D101
     def __init__(self, recorderID=-1):  # noqa: N803
-        if recorderID < 0:
+        if recorderID < 0 or recorderID is None:  # noqa: E701
             print(  # noqa: T201
                 'No recorder ID, or a negative ID, provided, defaulting to 0 for all forces.'
             )
@@ -73,15 +73,14 @@ def __init__(self, recorderID=-1):  # noqa: N803
             self.Y = []
             self.Z = []
             # prepend zeros to the list to account for the timeSeries transient analysis req in OpenSees
-            prependZero = False  # noqa: N806
+            prependZero = True  # noqa: N806
             if prependZero:
                 self.X.append(0.0)
                 self.Y.append(0.0)
                 self.Z.append(0.0)
 
             # Read in forces.[out or evt] file and add to EVENT.json
             # now using intermediary forces.evt for output of preceding Python calcs,
-            # prevents confusion with forces.out made by FEM tab
             if os.path.exists('forces.evt'):  # noqa: PTH110
                 with open('forces.evt') as file:  # noqa: PTH123
                     print('Reading forces from forces.evt to EVENT.json')  # noqa: T201
@@ -93,17 +92,12 @@ def __init__(self, recorderID=-1):  # noqa: N803
                         if not strip_line:
                             print('Empty line found in forces.evt... skip')  # noqa: T201
                             continue
-                        # Assume there is no header in the file
-                        # Assume recorder IDs are sequential, starting from 1
+                        # Assume recorder IDs are sequential
+                        # Often, the first ID (0) is set to 0.0 for all time-steps, as it usually maps to a rigid node at the structures base
                         if (j) == recorderID:
-                            # Strip away leading / trailing white-space,
-                            # Delimit by regex to capture " ", \s, "  ", tabs, etc.
+                            # 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
-                            # clean_line = re.split() # default is '\s+', which is any whitespace
                             clean_line = re.split(r';\s|;|,\s|,|\s+', strip_line)
-                            # clean_line = re.split(r';|,\s', strip_line)
-                            # clean_line = re.split("\s+", strip_line)
-
                             for k in range(len(clean_line)):
                                 self.X.append(float(clean_line[k]))
                                 self.Y.append(0.0)
@@ -121,23 +115,18 @@ def __init__(self, recorderID=-1):  # noqa: N803
                     self.Y = [0.0]
                     self.Z = [0.0]
                 else:
-                    # for a timeSeries with N elements, we append an element at N+1 to represent the max force of the series
-                    self.X.append(max(self.X))
-                    self.Y.append(max(self.Y))
-                    self.Z.append(max(self.Z))
-
                     print(  # noqa: T201
-                        'Length: ',
+                        'Force time-series length: ',
                         len(self.X),
-                        ', Max force: ',
+                        ', Max force (X,Y,Z): ',
                         max(self.X),
                         max(self.Y),
                         max(self.Z),
-                        ', Min force: ',
+                        ', Min force (X,Y,Z): ',
                         min(self.X),
                         min(self.Y),
                         min(self.Z),
-                        ', Last force: ',
+                        ', Last force (X,Y,Z): ',
                         self.X[-1],
                         self.Y[-1],
                         self.Z[-1],
@@ -180,6 +169,7 @@ def addFloorForceToEvent(  # noqa: N802
         'timeSeries': seriesName,
         'numSteps': len(force.X),
         'dT': dt,
+        'dt': dt,
         'type': 'WindFloorLoad',
         'floor': str(floor),
         'story': str(floor),
@@ -262,26 +252,55 @@ def GetFloorsCount(BIMFilePath):  # noqa: N802, N803, D103
     return int(bim['GeneralInformation']['stories'])
 
 
-def GetCelerisScript():  # noqa: N802, N803, D103
-    # filePath = BIMFilePath  # noqa: N806
-    # with open(filePath, encoding='utf-8') as file:  # noqa: PTH123
-    #     evt = json.load(file)
-    # file.close  # noqa: B018
-
-    # fileNameKey = 'simulationScript'  # noqa: N806
-    # filePathKey = fileNameKey + 'Path'  # noqa: N806
-
-    # for event in evt['Events']:
-    #     fileName = event[fileNameKey]  # noqa: N806
-    #     filePath = event[filePathKey]  # noqa: N806
-    #     return os.path.join(filePath, fileName)  # noqa: PTH118
-
+def GetCelerisScript(defaultScriptName = 'setrun.py'):  # noqa: N802, N803, D103
     defaultScriptPath = f'{os.path.realpath(os.path.dirname(__file__))}'  # noqa: N806, PTH120
-    defaultScriptName = 'setrun.py'  # noqa: N806
     defaultScriptPath = os.path.join(defaultScriptPath, defaultScriptName)  # noqa: PTH118
     return defaultScriptPath  # noqa: PTH118
 
 
+def dt(event=None):  # noqa: D102
+    """
+    Computes the time step based on the Courant criterion:
+        dt = Courant * dx / sqrt(g * maxdepth)
+
+    Returns:
+        float: The computed time step.
+    """
+    if 'Courant_num' in event['config']:
+        # Check is string
+        if isinstance(event['config']['Courant_num'], str):
+            Courant = 0.1
+        else:
+            Courant = event['config']['Courant_num']
+    else:
+        Courant = 0.1
+        
+    if 'base_depth' in event['config']:
+        # Check is string
+        if isinstance(event['config']['base_depth'], str):
+            maxdepth = 1.0
+        else:
+            maxdepth = event['config']['base_depth']
+    elif 'seaLevel' in event['config']:
+        # Check is string
+        if isinstance(event['config']['seaLevel'], str):
+            maxdepth = 1.0
+        else:
+            maxdepth = event['config']['seaLevel']
+    else:
+        maxdepth = 1.0
+    
+    if (maxdepth <= 0.0):
+        print('Warning: maxdepth is less than or equal to 0.0, setting it to 1.0. This will affect the time-step and simulation stability')  # noqa: T201
+        maxdepth = 1.0
+    dx = abs(event['config']['dx']) if 'dx' in event['config'] else 1.0
+    gravity = abs(event['config']['g']) if 'g' in event['config'] else 9.80665  # Acceleration due to gravity in m/s^2
+    dt = Courant * dx / np.sqrt(gravity * maxdepth)
+    if dt <= 0.0:
+        print('Warning: Computed dt is less than or equal to 0.0, setting it to 1.0e-3')  # noqa: T201
+        dt = 1.0e-3
+    return dt  # noqa: N806
+
 def main():
     """
     Entry point to generate event file using Celeris.
@@ -335,10 +354,18 @@ 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
+        # Determine dt for the force time series
+        # Try to compute using Courant_num (CFL), otherwise look for dt in the config
+        if 'Courant_num' in event['config']:
+            print('Courant_num found in event file. Compute dt.')  # noqa: T201
+            dt = dt(event=event)
+        else:
+            print('Courant_num not found in event file. Use provided dt')  # noqa: T201
+            if 'dt' in event['config']:
+                dt = event['config']['dt']  # noqa: N816
+            else:
+                print('dt not found in event file. Use default dt')  # noqa: T201
+                dt = 1e-3  # noqa: N806
 
     print('Running Celeris with script:', scriptName)  # noqa: T201
     print('Running Celeris with directory:', caseDirectory)  # noqa: T201
@@ -397,4 +424,3 @@ def main():
 
         # write the event file
         writeEVENT(forces, filenameEVENT, floorsCount=floorsCount)
-        # writeEVENT(forces, arguments.filenameEVENT)

modules/createEVENT/Celeris/celeris/domain.py

@@ -225,7 +225,7 @@ def __init__(  # noqa: C901
         self.WaveType = WaveType  # -1 to read from a file
         self.sine_wave = sine_wave
         self.data = np.zeros((2, 4))
-        self.N_data = 2
+        self.N_data = 1
         self.W_data = None
         self.init_eta = init_eta
         self.amplitude = Amplitude
@@ -242,6 +242,12 @@ def __init__(  # noqa: C901
                 self.amplitude = float(self.configfile['amplitude'])
             if checjson('period', self.configfile) == 1:
                 self.period = float(self.configfile['period'])
+            if checjson('sine_wave', self.configfile) == 1:
+                self.sine_wave = [float(self.configfile['sine_wave'][0]),
+                                 float(self.configfile['sine_wave'][1]),
+                                 float(self.configfile['sine_wave'][2]),
+                                 float(self.configfile['sine_wave'][3])]
+                self.N_data = 1 # Just one wave
             if checjson('BoundaryWidth', self.configfile) == 1:
                 self.BoundaryWidth = int(self.configfile['BoundaryWidth'])
             else:
@@ -343,10 +349,17 @@ def get_data(self):  # noqa: D102
         Returns:
             ti.field: A Taichi field of shape `(N_data, 4)` containing wave parameters.
         """
+        data = None  # noqa: N806
         if self.WaveType == -1:
             self.load_data()
-        data = ti.field(self.precision, shape=(self.N_data, 4))
-        data.from_numpy(self.data)
+            data = ti.field(self.precision, shape=(self.N_data, 4))
+            data.from_numpy(self.data)
+        elif self.WaveType == 1 or self.WaveType == 2:
+            data = ti.field(self.precision, shape=(self.N_data, 4))
+            data.from_numpy(np.array([self.sine_wave]))
+        else:
+            data = ti.field(self.precision, shape=(self.N_data, 4))
+            data.from_numpy(np.array([[0, 0, 0, 0]]))
         return data
 
 
@@ -478,11 +491,19 @@ def __init__(  # noqa: C901, PLR0913
             if checjson('WIDTH', self.configfile) == 1:
                 self.Nx = int(self.configfile['WIDTH'])
             else:
-                self.Nx = Nx
+                if (self.topodata != None) and (self.topodata.datatype == 'celeris' or self.topodata.datatype == 'txt'):
+                    self.Nx = self.topodata.z().shape[0]
+                else:
+                    self.Nx = Nx
+                self.configfile["WIDTH"] = self.Nx    
             if checjson('HEIGHT', self.configfile) == 1:
                 self.Ny = int(self.configfile['HEIGHT'])
             else:
-                self.Ny = Ny
+                if (self.topodata != None) and (self.topodata.datatype == 'celeris' or self.topodata.datatype == 'txt'):
+                    self.Ny = self.topodata.z().shape[1]
+                else:
+                    self.Ny = Ny
+                self.configfile["HEIGHT"] = self.Ny
             if checjson('dx', self.configfile) == 1:
                 self.dx = float(self.configfile['dx'])
             else:

modules/createEVENT/Celeris/celeris/solver.py

@@ -371,8 +371,8 @@ def __init__(  # noqa: C901, PLR0913, PLR0915
                 self.force_sensor_begin = self.bc.configfile['force_sensor_begin']
                 self.force_sensor_begin[0] = int(self.force_sensor_begin[0] / self.dx)
                 self.force_sensor_begin[1] = int(self.force_sensor_begin[1] / self.dy)
-                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']
+                self.force_sensor_begin_scaled[0] = self.force_sensor_begin[0] / self.domain.Nx
+                self.force_sensor_begin_scaled[1] = self.force_sensor_begin[1] / self.domain.Ny
             else:
                 self.force_sensor_begin = [0.0, 0.0]
                 self.force_sensor_begin_scaled = [0.0, 0.0]
@@ -381,10 +381,10 @@ def __init__(  # noqa: C901, PLR0913, PLR0915
                 self.force_sensor_end = self.bc.configfile['force_sensor_end']
                 self.force_sensor_end[0] = int(self.force_sensor_end[0] / self.dx)
                 self.force_sensor_end[1] = int(self.force_sensor_end[1] / self.dy)
-                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']
+                self.force_sensor_end_scaled[0] = self.force_sensor_end[0] / self.domain.Nx
+                self.force_sensor_end_scaled[1] = self.force_sensor_end[1] / self.domain.Ny
             else:
-                self.force_sensor_end = [self.bc.configfile['WIDTH'], self.bc.configfile['HEIGHT']]
+                self.force_sensor_end = [self.domain.Nx, self.domain.Ny]
                 self.force_sensor_end_scaled = [1.0, 1.0]
             if checjson('NumberOfTimeSeries', self.bc.configfile) == 1:
                 self.num_wave_gauges = int(self.bc.configfile['NumberOfTimeSeries'])
@@ -401,27 +401,27 @@ def __init__(  # noqa: C901, PLR0913, PLR0915
                     locationOfTimeSeries[1] = self.bc.configfile['locationOfTimeSeries'][ts]["yts"]
                     self.wave_gauge_pixels[ts,0] = int(locationOfTimeSeries[0] / self.dx)
                     self.wave_gauge_pixels[ts,1] = int(locationOfTimeSeries[1] / self.dy)
-                    self.wave_gauge_scaled[ts,0] = self.wave_gauge_pixels[ts,0] / self.bc.configfile['WIDTH']
-                    self.wave_gauge_scaled[ts,1] = self.wave_gauge_pixels[ts,1] / self.bc.configfile['HEIGHT']
+                    self.wave_gauge_scaled[ts,0] = self.wave_gauge_pixels[ts,0] / self.domain.Nx
+                    self.wave_gauge_scaled[ts,1] = self.wave_gauge_pixels[ts,1] / self.domain.Ny
                     print("Wave gauge location: ", locationOfTimeSeries[0], locationOfTimeSeries[1])
                     print("Wave gauge pixel location: ", self.wave_gauge_pixels[ts,0], self.wave_gauge_pixels[ts,1])
                     print("Wave gauge scaled location: ", self.wave_gauge_scaled[ts,0], self.wave_gauge_scaled[ts,1])
 
                     self.velocimeter_pixels[ts,0] = int(locationOfTimeSeries[0] / self.dx)
                     self.velocimeter_pixels[ts,1] = int(locationOfTimeSeries[1] / self.dy)
-                    self.velocimeter_scaled[ts,0] = self.velocimeter_pixels[ts,0] / self.bc.configfile['WIDTH']
-                    self.velocimeter_scaled[ts,1] = self.velocimeter_pixels[ts,1] / self.bc.configfile['HEIGHT']
+                    self.velocimeter_scaled[ts,0] = self.velocimeter_pixels[ts,0] / self.domain.Nx
+                    self.velocimeter_scaled[ts,1] = self.velocimeter_pixels[ts,1] / self.domain.Ny
 
             else:
                 self.wave_gauge_pixels[0,0] = int(1)
                 self.wave_gauge_pixels[0,1] = int(1)
-                self.wave_gauge_scaled[0,0] = self.wave_gauge_pixels[0,0] / self.bc.configfile['WIDTH']
-                self.wave_gauge_scaled[0,1] = self.wave_gauge_pixels[0,1] / self.bc.configfile['HEIGHT']
+                self.wave_gauge_scaled[0,0] = self.wave_gauge_pixels[0,0] / self.domain.Nx
+                self.wave_gauge_scaled[0,1] = self.wave_gauge_pixels[0,1] / self.domain.Ny
 
                 self.velocimeter_pixels[0,0] = int(1)
                 self.velocimeter_pixels[0,1] = int(1)
-                self.velocimeter_scaled[0,0] = self.velocimeter_pixels[0,0] / self.bc.configfile['WIDTH']
-                self.velocimeter_scaled[0,1] = self.velocimeter_pixels[0,1] / self.bc.configfile['HEIGHT']
+                self.velocimeter_scaled[0,0] = self.velocimeter_pixels[0,0] / self.domain.Nx
+                self.velocimeter_scaled[0,1] = self.velocimeter_pixels[0,1] / self.domain.Ny
 
 
             if checjson('whiteWaterDecayRate', self.bc.configfile) == 1: