Skip to content

Update Celeris #396

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
fmckenna merged 2 commits into from
May 29, 2025
Merged

Update Celeris #396

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4bb2583
Update Celeris
JustinBonus May 28, 2025
95adc17
Fix Celeris sine_waves definition messing with things
JustinBonus May 28, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
106 changes: 66 additions & 40 deletions modules/createEVENT/Celeris/Celeris.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.'
)
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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],
Expand Down Expand Up @@ -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),
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -397,4 +424,3 @@ def main():

# write the event file
writeEVENT(forces, filenameEVENT, floorsCount=floorsCount)
# writeEVENT(forces, arguments.filenameEVENT)
31 changes: 26 additions & 5 deletions modules/createEVENT/Celeris/celeris/domain.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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:
Expand Down Expand Up @@ -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


Expand Down Expand Up @@ -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:
Expand Down
26 changes: 13 additions & 13 deletions modules/createEVENT/Celeris/celeris/solver.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -371,8 +371,8 @@
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]
Expand All @@ -381,10 +381,10 @@
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'])
Expand All @@ -401,27 +401,27 @@
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:
Expand Down Expand Up @@ -1431,7 +1431,7 @@
# HNear_vec neighbours [hN,hE,hS,hW]
self.Hnear[ i , j ] = [0.0, h_east, 0.0 , h_west]

# To avoid unnecesary computations

Check failure on line 1434 in modules/createEVENT/Celeris/celeris/solver.py

View workflow job for this annotation

GitHub Actions / spell-check 

unnecesary ==> unnecessary
h_cut = self.delta
if h_here <= h_cut:
if h_east <= h_cut and h_west <= h_cut:
Expand Down
Loading