Source code for sharpy.solvers.lindynamicsim

import numpy as np
import os
import h5py as h5
from sharpy.utils.solver_interface import solver, BaseSolver, initialise_solver
import sharpy.utils.settings as settings
import sharpy.linear.src.libss as libss
import scipy.linalg as sclalg
import sharpy.utils.h5utils as h5utils
from sharpy.utils.datastructures import LinearTimeStepInfo
import sharpy.utils.cout_utils as cout
import time
import warnings

[docs]@solver class LinDynamicSim(BaseSolver): """Time-domain solution of Linear Time Invariant Systems Uses the derived linear time invariant systems and solves it in time domain. Requires a ``case_name.lininput.h5`` file in the case root folder that contains: * ``x0`` (optional): Initial state vector * ``input_vec``: Input vector ``(n_tsteps, n_inputs)``. Note: This solver is seldom used in SHARPy (its focus is on nonlinear time domain aeroelasticity) hence you may find this solver lacking in features. If you use it, you may need to make modifications. We would greatly appreciate that you contribute these modifications by means of a pull request! """ solver_id = 'LinDynamicSim' solver_classification = 'Coupled' settings_types = dict() settings_default = dict() settings_description = dict() settings_types['write_dat'] = 'list(str)' settings_default['write_dat'] = [] settings_description['write_dat'] = 'List of vectors to write: ``x``, ``y``, ``u`` and/or ``t``' settings_types['reference_velocity'] = 'float' settings_default['reference_velocity'] = 1. settings_description['reference_velocity'] = 'Velocity to scale the structural equations when using a non-dimensional system' settings_default['n_tsteps'] = 10 settings_types['n_tsteps'] = 'int' settings_description['n_tsteps'] = 'Number of time steps to run' settings_types['physical_time'] = 'float' settings_default['physical_time'] = 2. settings_description['physical_time'] = 'Time to run' settings_default['dt'] = 0.001 settings_types['dt'] = 'float' settings_description['dt'] = 'Time increment for the solution of systems without a specified dt' settings_types['postprocessors'] = 'list(str)' settings_default['postprocessors'] = list() settings_types['postprocessors_settings'] = 'dict' settings_default['postprocessors_settings'] = dict() settings_table = settings.SettingsTable() __doc__ += settings_table.generate(settings_types, settings_default, settings_description) def __init__(self): = None self.settings = dict() self.postprocessors = dict() self.with_postprocessors = False self.input_data_dict = dict() self.input_file_name = "" self.folder = None def initialise(self, data, custom_settings=None): = data if custom_settings: self.settings = custom_settings else: self.settings = data.settings[self.solver_id] settings.to_custom_types(self.settings, self.settings_types, self.settings_default) # Read initial state and input data and store in dictionary self.read_files() # Output folder self.folder = data.output_folder + '/lindynamicsim/' if not os.path.exists(self.folder): os.makedirs(self.folder) # initialise postprocessors self.postprocessors = dict() if len(self.settings['postprocessors']) > 0: self.with_postprocessors = True for postproc in self.settings['postprocessors']: self.postprocessors[postproc] = initialise_solver(postproc) self.postprocessors[postproc].initialise(, self.settings['postprocessors_settings'][postproc], caller=self) def run(self): ss = n_steps = self.settings['n_tsteps'] x0 = self.input_data_dict.get('x0', np.zeros(ss.states)) u = self.input_data_dict['u'] if len(x0) != ss.states: warnings.warn('Number of states in the initial state vector not equal to the number of states') x0 = np.zeros(ss.states) if u.shape[1] != ss.inputs: warnings.warn('Dimensions of the input vector not equal to the number of inputs') cout.cout_wrap('Number of inputs: %g' % ss.inputs, 3) cout.cout_wrap('Number of timesteps: %g' % n_steps, 3) cout.cout_wrap('Number of UVLM inputs: %g' %, 3) cout.cout_wrap('Number of beam inputs: %g' %, 3) breakpoint() try: dt = ss.dt except AttributeError: dt = self.settings['dt'] # Total time to run T = n_steps*dt u_ref = self.settings['reference_velocity'] # If the system is scaled: if u_ref != 1.: scaling_factors = dt_dimensional = scaling_factors['length'] / u_ref T_dimensional = n_steps * dt_dimensional T = T_dimensional / scaling_factors['time'] ss =['reference_velocity']) t_dom = np.linspace(0, T, n_steps) # Use the scipy linear solver sys = libss.ss_to_scipy(ss) cout.cout_wrap('Solving linear system using scipy...') t0 = time.time() # breakpoint() out = sys.output(u, t=t_dom, x0=x0) ts = time.time() - t0 cout.cout_wrap('\tSolved in %.2fs' % ts, 1) t_out = out[0] x_out = out[2] y_out = out[1] if self.settings['write_dat']: cout.cout_wrap('Writing linear simulation output .dat files to %s' % self.folder) if 'y' in self.settings['write_dat']: np.savetxt(self.folder + '/y_out.dat', y_out) cout.cout_wrap('Output vector written', 2) if 'x' in self.settings['write_dat']: np.savetxt(self.folder + '/x_out.dat', x_out) cout.cout_wrap('State vector written', 2) if 'u' in self.settings['write_dat']: np.savetxt(self.folder + '/u_out.dat', u) cout.cout_wrap('Input vector written', 2) if 't' in self.settings['write_dat']: np.savetxt(self.folder + '/t_out.dat', t_out) cout.cout_wrap('Time domain written', 2) cout.cout_wrap('Success', 1) # Pack state variables into linear timestep info cout.cout_wrap('Plotting results...') for n in range(len(t_out)-1): tstep = LinearTimeStepInfo() tstep.x = x_out[n, :] tstep.y = y_out[n, :] tstep.t = t_out[n] tstep.u = u[n, :] # TODO: option to save to h5 # Pack variables into respective aero or structural time step infos (with the + f0 from lin) # Need to obtain information from the variables in a similar fashion as done with the database # for the beam case aero_tstep, struct_tstep = state_to_timestep(, tstep.x, tstep.u, tstep.y) # run postprocessors if self.with_postprocessors: for postproc in self.postprocessors: = self.postprocessors[postproc].run(online=True) return def read_files(self): self.input_file_name =['SHARPy']['route'] + '/' +['SHARPy']['case'] + '.lininput.h5' # Check that the file exists try: h5utils.check_file_exists(self.input_file_name) # Read and store with h5.File(self.input_file_name, 'r') as input_file_handle: self.input_data_dict = h5utils.load_h5_in_dict(input_file_handle) except FileNotFoundError: pass
def state_to_timestep(data, x, u=None, y=None): """ Warnings: Under development Writes a state-space vector to SHARPy timesteps Args: data: x: u: y: Returns: """ if data.settings['LinearAssembler']['linear_system_settings']['beam_settings']['modal_projection'] and \ data.settings['LinearAssembler']['linear_system_settings']['beam_settings']['inout_coords'] == 'modes': modal = True else: modal = False # modal = True if data.linear.linear_system.uvlm.gust_assembler: start_x_aero = data.linear.linear_system.uvlm.gust_assembler.ss_gust.states else: start_x_aero = 0 x_aero = x[] x_struct = x[] # u_aero = TODO: external velocities phi = data.linear.linear_system.beam.sys.U Kas = data.linear.linear_system.couplings['Kas'] # Beam output y_beam = x_struct u_q = np.zeros( if u is not None: u_q += u[] u_q[:y_beam.shape[0]] += y_beam else: u_q[:y_beam.shape[0]] += y_beam if modal: # add eye matrix for extra inputs n_modes = phi.shape[1] n_inputs_aero_only = len(u_q) - 2*n_modes # Inputs to the UVLM other than structural inputs u_aero =, phi, np.eye(n_inputs_aero_only)).dot(u_q)) else: # if u_q.shape[0] != # u_aero_zero = data.linear.tsaero0 u_aero = # Unpack input zeta, zeta_dot, u_ext = data.linear.linear_system.uvlm.unpack_input_vector(u_aero) # Also add the beam forces. I have a feeling there is a minus there as well.... # Aero forces, gamma, gamma_dot, gamma_star = data.linear.linear_system.uvlm.unpack_ss_vector( data, x_n=x_aero, aero_tstep=data.linear.tsaero0, track_body=True) current_aero_tstep =[-1].copy() current_aero_tstep.forces = [forces[i_surf] + data.linear.tsaero0.forces[i_surf] for i_surf in range(len(gamma))] current_aero_tstep.gamma = [gamma[i_surf] + data.linear.tsaero0.gamma[i_surf] for i_surf in range(len(gamma))] current_aero_tstep.gamma_dot = [gamma_dot[i_surf] + data.linear.tsaero0.gamma_dot[i_surf] for i_surf in range(len(gamma))] current_aero_tstep.gamma_star = [gamma_star[i_surf] + data.linear.tsaero0.gamma_star[i_surf] for i_surf in range(len(gamma))] current_aero_tstep.zeta = zeta current_aero_tstep.zeta_dot = zeta_dot current_aero_tstep.u_ext = u_ext # aero_forces = beam_forces = data.linear.linear_system.couplings['Ksa'].dot(aero_forces) if u is not None: u_struct = u[] # y_struct = y[[sys_id].lsys['LinearBeam'].ss.outputs] # Reconstruct the state if modal if modal: phi = data.linear.linear_system.beam.sys.U x_s = sclalg.block_diag(phi, phi).dot(x_struct) else: x_s = x_struct y_s = beam_forces #+ # y_s =['LinearBeam'] current_struct_step = data.linear.linear_system.beam.unpack_ss_vector(x_s, y_s, data.linear.tsstruct0) # data.structure.timestep_info.append(current_struct_step) return current_aero_tstep, current_struct_step