« prev ^ index » next coverage.py v7.13.4, created at 2026-02-27 18:22 +0000
1# -*- coding: utf-8 -*-
2"""
3Synthetic "hardware" that allows the responses to be simulated by integrating
4modal equations of motion nearly real-time.
6Rattlesnake Vibration Control Software
7Copyright (C) 2021 National Technology & Engineering Solutions of Sandia, LLC
8(NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
9Government retains certain rights in this software.
11This program is free software: you can redistribute it and/or modify
12it under the terms of the GNU General Public License as published by
13the Free Software Foundation, either version 3 of the License, or
14(at your option) any later version.
16This program is distributed in the hope that it will be useful,
17but WITHOUT ANY WARRANTY; without even the implied warranty of
18MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19GNU General Public License for more details.
21You should have received a copy of the GNU General Public License
22along with this program. If not, see <https://www.gnu.org/licenses/>.
23"""
25import multiprocessing as mp
26import time
27from typing import List
29import netCDF4
30import numpy as np
31import scipy.signal as signal
33from .abstract_hardware import HardwareAcquisition, HardwareOutput
34from .utilities import Channel, DataAcquisitionParameters, flush_queue
36DEBUG = False
38NOISE_LEVEL = 0.00000001
40if DEBUG:
41 import glob
42 import os
44 TEST_FILES_OUTPUT_NAMES = r"debug\output_data_{:}.npz"
45 TEST_FILES_ACQUISITION_NAMES = r"debug\acquire_data_{:}.npz"
46 TEST_FILES_SYSTEM_NAME = r"debug\system.npz"
47 # Delete existing test files
48 for str_format in [TEST_FILES_ACQUISITION_NAMES, TEST_FILES_OUTPUT_NAMES]:
49 files = glob.glob(str_format.format("*"))
50 for file in files:
51 try:
52 os.remove(file)
53 except (FileNotFoundError, PermissionError):
54 pass
57class ExodusAcquisition(HardwareAcquisition):
58 """Class defining the interface between the controller and synthetic acquisition
60 This class defines the interfaces between the controller and the
61 data acquisition portion of the hardware. In this case, the hardware is
62 actually simulated by integrating the modal equations of motion of a finite
63 element model containing a modal solution. It is run by the Acquisition
64 process, and must define how to get data from the test hardware into the
65 controller."""
67 def __init__(self, exodus_file: str, queue: mp.queues.Queue):
68 """Loads in the Exodus file and sets initial parameters to null values
71 Parameters
72 ----------
73 exodus_file : str :
74 Path to the Exodus file.
75 queue : mp.queues.Queue
76 A queue that passes input data from the ExodusOutput class to the
77 Exodus input class. Normally, this data transfer would occur through
78 the physical test object: the exciters would excite the test object
79 with the specified excitation and the Acquisition would record the
80 responses to that excitation. In the synthetic case, we need to
81 pass the output data to the acquisition which does the integration.
83 """
84 self.exo = Exodus(exodus_file)
85 self.phi = None
86 self.phi_full = None
87 self.response_channels: np.ndarray
88 self.response_channels = None
89 self.system = None
90 self.times = None
91 self.state = None
92 self.frame_time = None
93 self.queue = queue
94 self.force_buffer = None
95 self.integration_oversample = None
96 self.acquisition_delay = None
97 self.damping = None
99 def set_up_data_acquisition_parameters_and_channels(
100 self, test_data: DataAcquisitionParameters, channel_data: List[Channel]
101 ):
102 """
103 Initialize the hardware and set up channels and sampling properties
105 The function must create channels on the hardware corresponding to
106 the channels in the test. It must also set the sampling rates.
108 Parameters
109 ----------
110 test_data : DataAcquisitionParameters :
111 A container containing the data acquisition parameters for the
112 controller set by the user.
113 channel_data : List[Channel] :
114 A list of ``Channel`` objects defining the channels in the test
116 Returns
117 -------
118 None.
120 """
121 self.create_response_channels(channel_data)
122 self.set_parameters(test_data)
124 def create_response_channels(self, channel_data: List[Channel]):
125 """Method to set up response channels
127 This function takes channels from the supplied list of channels and
128 extracts the mode shape coefficients corresponding to those channels.
130 Parameters
131 ----------
132 channel_data : List[Channel] :
133 A list of ``Channel`` objects defining the channels in the test
135 """
136 # print('{:} Channels'.format(len(channel_data)))
137 displacements = self.exo.get_displacements()
138 node_numbers = self.exo.get_node_num_map()
139 try:
140 self.damping = float(channel_data[0].comment)
141 print(f"{self.damping} Damping")
142 except ValueError:
143 self.damping = 0.01
144 self.response_channels = np.array(
145 [
146 channel.feedback_device is None or channel.feedback_device == ""
147 for channel in channel_data
148 ],
149 dtype="bool",
150 )
151 self.phi_full = np.array(
152 [self._create_channel(channel, displacements, node_numbers) for channel in channel_data]
153 )
154 # Need to add a signal buffer in case the write size is not equal to
155 # the read size
156 self.force_buffer = np.zeros((0, np.sum(~self.response_channels)))
158 def set_parameters(self, test_data: DataAcquisitionParameters):
159 """Method to set up sampling rate and other test parameters
161 For the synthetic case, we will set up the integration parameters using
162 the sample rates provided.
164 Parameters
165 ----------
166 test_data : DataAcquisitionParameters :
167 A container containing the data acquisition parameters for the
168 controller set by the user.
170 """
171 # Get the number of modes that we will keep (bandwidth*1.5)
172 frequencies = self.exo.get_times()
173 frequencies[frequencies < 0] = 0 # Eliminate any negative frequencies
174 keep_modes = frequencies < test_data.nyquist_frequency * 1.5
175 # Oversampling due to integration
176 self.integration_oversample = test_data.output_oversample
177 self.phi = self.phi_full[..., keep_modes]
178 nmodes = self.phi.shape[-1]
179 frequencies = frequencies[keep_modes]
180 m = np.eye(nmodes)
181 k = np.diag((frequencies * 2 * np.pi) ** 2)
182 c = np.diag((2 * 2 * np.pi * frequencies * self.damping))
183 a, b, c, d = mck_to_state_space(m, c, k)
184 self.system = signal.StateSpace(a, b, c, d)
185 # Need to get one more sample than you would think because lsim doesn't bridge the gap
186 # between integrations
187 self.times = np.arange(test_data.samples_per_read * self.integration_oversample + 1) / (
188 test_data.sample_rate * self.integration_oversample
189 )
190 self.frame_time = test_data.samples_per_read / test_data.sample_rate
191 self.state = np.zeros(nmodes * 2)
192 self.acquisition_delay = test_data.samples_per_write / test_data.output_oversample
193 if DEBUG:
194 np.savez(TEST_FILES_SYSTEM_NAME, a=a, b=b, c=c, d=d, times=self.times)
196 def start(self):
197 """Method to start acquiring data.
199 For the synthetic case, it simply initializes the state of the system to zero"""
200 self.state[:] = 0
202 def get_acquisition_delay(self) -> int:
203 """
204 Get the number of samples between output and acquisition.
206 This function returns the number of samples that need to be read to
207 ensure that the last output is read by the acquisition. If there is
208 buffering in the output, this delay should be adjusted accordingly.
210 Returns
211 -------
212 int
213 Number of samples between when a dataset is written to the output
214 and when it has finished playing.
216 """
217 return self.acquisition_delay
219 def read(self):
220 """Method to read a frame of data from the hardware
222 This function gets the force from the output queue and adds it to the
223 buffer of time signals that represents the force. It then integrates
224 a frame of time and sends it to the acquisition.
226 Returns
227 -------
228 read_data :
229 2D Data read from the controller with shape ``n_channels`` x
230 ``n_samples``
232 """
233 start_time = time.time()
234 while self.force_buffer.shape[0] < self.times.size:
235 try:
236 forces = self.queue.get(timeout=self.frame_time)
237 except mp.queues.Empty:
238 # If we don't get an output in time, this likely means output has stopped so
239 # just put zeros.
240 forces = np.zeros((self.force_buffer.shape[-1], self.times.size))
241 self.force_buffer = np.concatenate((self.force_buffer, forces.T), axis=0)
243 # Now extract a force that is the correct size
244 this_force = self.force_buffer[: self.times.size]
245 # And leave the rest for next time
246 # Note we have to keep the last force sample still on the
247 # buffer because it will be the next force sample we use
248 self.force_buffer = self.force_buffer[self.times.size - 1 :]
250 # print('Got Force')
251 # print('this_force shape: {:}'.format(this_force.shape))
253 modal_forces = np.einsum("ij,ki->kj", self.phi[~self.response_channels], this_force)
254 # print('modal_forces shape: {:}'.format(modal_forces.shape))
256 # print('Integrating...')
257 # print('system: {:}'.format(self.system))
258 # print('forces: {:}\n{:}'.format(modal_forces.shape,modal_forces))
259 # print('times: {:}\n{:}'.format(self.times.shape,self.times))
260 # print('state: {:}\n{:}'.format(self.state.shape,self.state))
262 times_out, sys_out, _ = signal.lsim(self.system, modal_forces, self.times, self.state)
263 # print('output: {:}\n{:}'.format(sys_out.shape,sys_out))
265 sys_accels = sys_out[:, 2 * self.phi.shape[-1] : 3 * self.phi.shape[-1]]
267 accelerations = self.phi[self.response_channels] @ sys_accels.T
269 self.state[:] = sys_out[-1, 0 : 2 * self.phi.shape[-1]]
271 # Now we need to combine accelerations with forces in the same way
272 output = np.zeros((len(self.response_channels), len(self.times))) # n channels x n times
273 output[self.response_channels] = accelerations
274 output[~self.response_channels] = this_force.T
276 # print('output: {:}\n{:}'.format(output.shape,output))
278 integration_time = time.time() - start_time
279 remaining_time = self.frame_time - integration_time
280 if remaining_time > 0.0:
281 time.sleep(remaining_time)
283 if DEBUG:
284 # Find current files
285 num_files = len(glob.glob(TEST_FILES_ACQUISITION_NAMES.format("*")))
286 np.savez(
287 TEST_FILES_ACQUISITION_NAMES.format(num_files),
288 full_output=output,
289 integration_oversample=self.integration_oversample,
290 modal_forces=modal_forces,
291 forces=this_force,
292 integration_output=sys_out,
293 times=times_out,
294 state=self.state,
295 response_channels=self.response_channels,
296 )
297 # We don't want to return the last sample because it
298 # will be the initial state for the next sample
299 return output[..., : -1 : self.integration_oversample] + NOISE_LEVEL * np.random.randn(
300 *output[..., : -1 : self.integration_oversample].shape
301 )
303 def read_remaining(self):
304 """Method to read the rest of the data on the acquisition
306 This function simply returns one sample of zeros.
308 Returns
309 -------
310 read_data :
311 2D Data read from the controller with shape ``n_channels`` x
312 ``n_samples``
313 """
314 return np.zeros((len(self.response_channels), 1))
316 def stop(self):
317 """Method to stop the acquisition.
319 This simply sets the state to zero."""
320 self.state[:] = 0
322 def close(self):
323 """Method to close down the hardware
325 This simply closes the Exodus file."""
326 self.exo.close()
328 def _create_channel(self, channel: Channel, displacement, node_numbers):
329 """Helper function to create a channel from the Exodus file.
331 This function parses the channel information and then extracts the
332 mode shape row corresponding to that channel.
334 Parameters
335 ----------
336 channel: Channel :
337 ``Channel`` object specifying the channel information
339 displacement :
340 The 3D mode shape matrix consisting of direction x mode x node
341 dimensions.
343 node_numbers :
344 The node numbers in the finite element model, either increasing
345 from one to the largest node or from the node num map.
347 Returns
348 -------
349 phi_row :
350 A row of the mode shape matrix corresponding to the channel
352 """
353 node_number = int(channel.node_number)
354 node_index = np.where(node_numbers == node_number)[0][0]
355 if channel.node_direction.lower().replace(" ", "") in ["x+", "+x"]:
356 direction = np.array([1, 0, 0])
357 elif channel.node_direction.lower().replace(" ", "") in ["x-", "-x"]:
358 direction = np.array([-1, 0, 0])
359 elif channel.node_direction.lower().replace(" ", "") in ["y+", "+y"]:
360 direction = np.array([0, 1, 0])
361 elif channel.node_direction.lower().replace(" ", "") in ["y-", "-y"]:
362 direction = np.array([0, -1, 0])
363 elif channel.node_direction.lower().replace(" ", "") in ["z+", "+z"]:
364 direction = np.array([0, 0, 1])
365 elif channel.node_direction.lower().replace(" ", "") in ["z-", "-z"]:
366 direction = np.array([0, 0, -1])
367 else:
368 direction = np.array([float(val) for val in channel.node_direction.split(",")])
369 direction /= np.linalg.norm(direction)
370 phi_row = np.einsum("i,ij", direction, displacement[..., node_index])
371 return phi_row
374class ExodusOutput(HardwareOutput):
375 """Class defining the interface between the controller and synthetic output
377 Note that the only thing that this class does is pass data to the acquisition
378 hardware task which actually performs the integration. Therefore, many of
379 the functions here are actually empty."""
381 def __init__(self, queue: mp.queues.Queue):
382 """
383 Initializes the hardware by simply storing the data passing queue.
385 Parameters
386 ----------
387 queue : mp.queues.Queue
388 Queue used to pass data from output to acquisition for integration.
389 See ``ExodusAcquisition.__init__``
391 """
392 self.queue = queue
394 def set_up_data_output_parameters_and_channels(
395 self, test_data: DataAcquisitionParameters, channel_data: List[Channel]
396 ):
397 """
398 Initialize the hardware and set up sources and sampling properties
400 This does nothing for the synthetic hardware
402 Parameters
403 ----------
404 test_data : DataAcquisitionParameters :
405 A container containing the data acquisition parameters for the
406 controller set by the user.
407 channel_data : List[Channel] :
408 A list of ``Channel`` objects defining the channels in the test
410 Returns
411 -------
412 None.
414 """
416 def start(self):
417 """Method to start acquiring data
419 Does nothing for synthetic hardware."""
421 def write(self, data: np.ndarray):
422 """Method to write a frame of data
424 For the synthetic excitation, this simply puts the data into the data-
425 passing queue.
427 Parameters
428 ----------
429 data : np.ndarray
430 Data to write to the output.
432 """
433 if DEBUG:
434 # Find current files
435 num_files = len(glob.glob(TEST_FILES_OUTPUT_NAMES.format("*")))
436 np.savez(TEST_FILES_OUTPUT_NAMES.format(num_files), forces=data)
437 self.queue.put(data)
439 def stop(self):
440 """Method to stop the acquisition
442 Does nothing for synthetic hardware."""
443 flush_queue(self.queue)
445 def close(self):
446 """Method to close down the hardware
448 Does nothing for synthetic hardware."""
450 def ready_for_new_output(self):
451 """Signals that the hardware is ready for new output
453 Returns ``True`` if the data-passing queue is empty.
454 """
455 return self.queue.empty()
458class ExodusError(Exception):
459 """An exception to specify an error has occured in the Exodus reader"""
462class Exodus:
463 """A class with limited Exodus file read capabilities"""
465 def __init__(self, filename):
466 """
467 Read in an Exodus file given by ``filename``
469 Parameters
470 ----------
471 filename : str
472 Path to the Exodus file to read.
474 """
475 self.filename = filename
476 self._ncdf_handle = netCDF4.Dataset(filename, "r") # pylint: disable=no-member
478 @property
479 def num_dimensions(self):
480 """Property to get the number of dimensions in an Exodus file"""
481 return self._ncdf_handle.dimensions["num_dim"].size
483 @property
484 def num_nodes(self):
485 """Property to get the number of nodes in an Exodus file"""
486 return self._ncdf_handle.dimensions["num_nodes"].size
488 @property
489 def num_times(self):
490 """Property to get the number of time steps in an Exodus file"""
491 return self._ncdf_handle.dimensions["time_step"].size
493 def get_node_num_map(self):
494 """Retrieve the list of local node IDs from the Exodus file.
496 Returns
497 -------
498 node_num_map : np.array
499 Returns a 1D array with size num_nodes, denoting the node number
500 for the node in each index
502 Notes
503 -----
504 If there is no node_num_map in the Exodus file, this function simply
505 returns an array from 1 to self.num_nodes
506 """
507 if "node_num_map" in self._ncdf_handle.variables:
508 return self._ncdf_handle.variables["node_num_map"][:]
509 else:
510 return np.ma.MaskedArray(np.arange(self.num_nodes) + 1)
512 def get_displacements(self, displacement_name="Disp", capital_coordinates=True):
513 """Get the displacements from the Exodus file.
515 Parameters
516 ----------
517 displacement_name :
518 The prefix of the displacement variable (Default value = 'Disp')
519 capital_coordinates :
520 Whether or not the X,Y,Z values appended to the displacement name
521 are capitalized. (Default value = True)
523 Returns
524 -------
525 displacement_array : np.ndarray
526 Returns a 3D array with size num_dimensions, num_times,num_nodes
527 containing the displacements in each direction at each mode shape
528 at each node.
529 """
530 return np.array(
531 [
532 self.get_node_variable_values(
533 displacement_name + (val.upper() if capital_coordinates else val.lower())
534 )
535 for val in "xyz"[: self.num_dimensions]
536 ]
537 )
539 def get_coords(self):
540 """Retrieve the coordinates of the nodes in the Exodus file.
542 Returns
543 -------
544 coordinate_array : np.ndarray
545 Returns a 2D array with size num_dimensions, num_nodes
546 containing the position of each node.
547 """
548 # TODO Add error checking
549 coord_names = ("coordx", "coordy", "coordz")[: self.num_dimensions]
550 raw_list = [self._ncdf_handle.variables[name][:] for name in coord_names]
551 coords = np.array(raw_list)
552 return coords
554 def get_node_variable_names(self):
555 """Gets a tuple of nodal variable names from the Exodus file
557 Returns
558 -------
559 node_var_names : tuple(str)
560 A tuple of strings corresponding to the variable names in the
561 model."""
562 try:
563 raw_records = self._ncdf_handle.variables["name_nod_var"]
564 except KeyError as e:
565 raise ExodusError("Node Variable Names are not defined!") from e
566 node_var_names = tuple(
567 "".join(
568 value.decode() for value in line if not isinstance(value, np.ma.core.MaskedConstant)
569 )
570 for line in raw_records
571 )
572 return node_var_names
574 def get_node_variable_values(self, name_or_index, step=None):
575 """Gets the node variable values for the specified timestep
577 Parameters
578 ----------
579 name_or_index : str or int
580 Name or Index of the nodal variable that is desired. If
581 type(name_or_index) == str, then it is assumed to be the name. If
582 type(name_or_index) == int, then it is assumed to be the index.
583 step : int
584 Time step at which to recover the nodal variable (Default value = None)
586 Returns
587 -------
588 node_var_values : np.ndarray
589 A 1D or 2D numpy array of variable values depending on whether a
590 time step is specified or not.
592 """
593 if isinstance(name_or_index, (int, np.integer)):
594 index = name_or_index
595 elif isinstance(name_or_index, (str, np.character)):
596 try:
597 index = self.get_node_variable_names().index(name_or_index)
598 except ValueError as e:
599 raise ExodusError(
600 f"Name {name_or_index} not found in self.get_node_variable_names(). "
601 f"Options are {self.get_node_variable_names()}"
602 ) from e
603 else:
604 raise ValueError("name_or_index must be integer or string")
605 vals_nod_var_name = f"vals_nod_var{index + 1:d}"
606 if step is not None:
607 if step >= self.num_times:
608 raise ExodusError("Invalid Time Step")
609 return self._ncdf_handle.variables[vals_nod_var_name][step, :]
610 return self._ncdf_handle.variables[vals_nod_var_name][:]
612 def get_times(self):
613 """Gets the time values from the Exodus file
615 Returns
616 -------
617 times : np.ndarray
618 A vector of time values in the Exodus file. These may be
619 frequencies if the Exodus file contains a modal or frequency-based
620 solution.
621 """
622 return self._ncdf_handle.variables["time_whole"][:]
624 def close(self):
625 """Closes the Exodus file"""
626 self._ncdf_handle.close()
629def mck_to_state_space(m, c, k):
630 """Creates a state-space representation from a mass, stiffness, and damping matrix.
632 Parameters
633 ----------
634 m : np.ndarray
635 A square array defining the system mass matrix.
636 c : np.ndarray
637 A square array defining the system damping matrix
638 k : np.ndarray
639 A square array defining the system stiffness matrix
641 Returns
642 -------
643 a : np.ndarray
644 The state space A matrix
645 b : np.ndarray
646 The state space B matrix
647 c : np.ndarray
648 The state space C matrix
649 d : np.ndarray
650 The state space D matrix
651 """
653 ndofs = m.shape[0]
655 # a = [[ 0, I],
656 # [-m^-1*k,-m^-1*c]]
658 a_state = np.block(
659 [
660 [np.zeros((ndofs, ndofs)), np.eye(ndofs)],
661 [-np.linalg.solve(m, k), -np.linalg.solve(m, c)],
662 ]
663 )
665 # b = [[ 0, m^-1]]
667 b_state = np.block([[np.zeros((ndofs, ndofs))], [np.linalg.inv(m)]])
669 # c = [[ I, 0], # Displacements
670 # [ 0, I], # Velocities
671 # [-m^-1*k,-m^-1*c], # Accelerations
672 # [ 0, 0]] # Forces
674 c_state = np.block(
675 [
676 [np.eye(ndofs), np.zeros((ndofs, ndofs))],
677 [np.zeros((ndofs, ndofs)), np.eye(ndofs)],
678 [-np.linalg.solve(m, k), -np.linalg.solve(m, c)],
679 [np.zeros((ndofs, ndofs)), np.zeros((ndofs, ndofs))],
680 ]
681 )
683 # d = [[ 0], # Displacements
684 # [ 0], # Velocities
685 # [ m^-1], # Accelerations
686 # [ I]] # Forces
688 d_state = np.block(
689 [
690 [np.zeros((ndofs, ndofs))],
691 [np.zeros((ndofs, ndofs))],
692 [np.linalg.inv(m)],
693 [np.eye(ndofs)],
694 ]
695 )
697 return a_state, b_state, c_state, d_state
700# def integrate_MCK(m,c,k,times,forces,x0=None):
701# """Integrates a system defined by a Mass, Stiffness, and Damping Matrix
703# Parameters
704# ----------
705# m : np.ndarray
706# A square array defining the system mass matrix.
707# c : np.ndarray
708# A square array defining the system damping matrix
709# k : np.ndarray
710# A square array defining the system stiffness matrix
711# times :
712# A vector of times that will be passed to the ``scipy.signal.lsim``
713# function.
714# forces :
715# A vector of forces that will be passed to the ``scipy.signal.lsim``
716# function
717# x0 :
718# The initial state of the system that will be used by the ``scipy.signal.lsim``
719# function (Default value = None, meaning the system starts at rest)
721# Returns
722# -------
724# """
725# a,b,c,d = mck_to_state_space(m,c,k)
726# linear_system = signal.StateSpace(a,b,c,d)
728# times_out,sys_out,x_out = signal.lsim(linear_system,forces,times,x0)
730# sys_disps = sys_out[:,0*m.shape[0]:1*m.shape[0]]
731# sys_vels = sys_out[:,1*m.shape[0]:2*m.shape[0]]
732# sys_accels = sys_out[:,2*m.shape[0]:3*m.shape[0]]
733# sys_forces = sys_out[:,3*m.shape[0]:4*m.shape[0]]
735# return sys_disps,sys_vels,sys_accels,sys_forces
737# def integrate_modes(natural_frequencies,damping_ratios,modal_forces,times,q0,qd0):
738# """Integrate a set of natural frequencies, damping ratios, and modal forces
740# Parameters
741# ----------
742# natural_frequencies : np.ndarray
743# A 1D array of natural frequencies for each mode
744# damping_ratios : np.ndarray
745# A 1D array of critical damping ratios for each mode
746# modal_forces : np.ndarray
747# A set of modal forces that will be passed to the integrate_MCK function
748# times : np.ndarray
749# A set of times that will be passed to the integrate_MCK function
750# q0 : np.ndarray
751# Initial values of the modal coefficients
752# qd0 : np.ndarray
753# Initial velocities of the modal coefficients
756# Returns
757# -------
759# """
760# nmodes = natural_frequencies.size
761# m = np.eye(nmodes)
762# k = np.diag((natural_frequencies*2*np.pi)**2)
763# c = np.diag((2*2*np.pi*natural_frequencies*damping_ratios))
764# q = integrate_MCK(m,c,k,times,modal_forces,np.concatenate((q0,qd0)))[:2]
765# return q