diff --git a/modules/createEVENT/Celeris/Celeris.py b/modules/createEVENT/Celeris/Celeris.py index a555df1d2..ce5e2eaff 100644 --- a/modules/createEVENT/Celeris/Celeris.py +++ b/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,7 +73,7 @@ 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) @@ -81,7 +81,6 @@ def __init__(self, recorderID=-1): # noqa: N803 # 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) diff --git a/modules/createEVENT/Celeris/celeris/domain.py b/modules/createEVENT/Celeris/celeris/domain.py index 91ff8d7a1..51faa60b0 100644 --- a/modules/createEVENT/Celeris/celeris/domain.py +++ b/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: diff --git a/modules/createEVENT/Celeris/celeris/solver.py b/modules/createEVENT/Celeris/celeris/solver.py index 6eb388cf3..7f4187655 100644 --- a/modules/createEVENT/Celeris/celeris/solver.py +++ b/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: