diff --git a/cadet/cadet.py b/cadet/cadet.py index a705ebe..b2ec514 100644 --- a/cadet/cadet.py +++ b/cadet/cadet.py @@ -237,11 +237,18 @@ def load_results(self): def run(self, timeout=None, check=None): data = self.cadet_runner.run(simulation=self.root.input, filename=self.filename, timeout=timeout, check=check) + # TODO: Why is this commented out? # self.return_information = data return data - def run_load(self, timeout=None, check=None, clear=True): - data = self.cadet_runner.run(simulation=self.root.input, filename=self.filename, timeout=timeout, check=check) + def run_load(self, timeout = None, check=None, clear=True): + data = self.cadet_runner.run( + simulation=self.root.input, + filename=self.filename, + timeout=timeout, + check=check + ) + # TODO: Why is this commented out? # self.return_information = data self.load_results() if clear: diff --git a/cadet/cadet_dll.py b/cadet/cadet_dll.py index 88920c5..ae1c36c 100644 --- a/cadet/cadet_dll.py +++ b/cadet/cadet_dll.py @@ -1,85 +1,126 @@ - -from cmath import sin import ctypes import numpy import addict import cadet.cadet_dll_parameterprovider as cadet_dll_parameterprovider + def log_handler(file, func, line, level, level_name, message): log_print('{} ({}:{:d}) {}'.format(level_name.decode('utf-8') , func.decode('utf-8') , line, message.decode('utf-8') )) -c_cadet_result = ctypes.c_int - -array_double = ctypes.POINTER(ctypes.POINTER(ctypes.c_double)) - - -point_int = ctypes.POINTER(ctypes.c_int) - -def null(*args): +def _no_log_output(*args): pass if 0: log_print = print else: - log_print = null + log_print = _no_log_output + +# Some common types +CadetDriver = ctypes.c_void_p +array_double = ctypes.POINTER(ctypes.POINTER(ctypes.c_double)) +point_int = ctypes.POINTER(ctypes.c_int) + +# Values of cdtResult +# TODO: Convert to lookup table to improve error messages below. +c_cadet_result = ctypes.c_int +_CDT_OK = 0 +_CDT_ERROR = -1 +_CDT_ERROR_INVALID_INPUTS = -2 +_CDT_DATA_NOT_STORED = -3 class CADETAPIV010000_DATA(): - _data_ = {} - _data_['createDriver'] = ('drv',) - _data_['deleteDriver'] = (None, 'drv') - _data_['runSimulation'] = ('return', 'drv', 'parameterProvider') - _data_['getNParTypes'] = ('return', 'drv', 'unitOpId', 'nParTypes') - _data_['getSolutionInlet'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nPort', 'nComp') - _data_['getSolutionOutlet'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nPort', 'nComp') - _data_['getSolutionBulk'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nComp') - _data_['getSolutionParticle'] = ('return', 'drv', 'unitOpId', 'parType', 'time', 'data', 'nTime', 'nParShells', 'nAxialCells', 'nRadialCells', 'nComp') - _data_['getSolutionSolid'] = ('return', 'drv', 'unitOpId', 'parType', 'time', 'data', 'nTime', 'nParShells', 'nAxialCells', 'nRadialCells', 'nBound') - _data_['getSolutionFlux'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nComp') - _data_['getSolutionVolume'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime') - _data_['getSolutionDerivativeInlet'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nPort', 'nComp') - _data_['getSolutionDerivativeOutlet'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nPort', 'nComp') - _data_['getSolutionDerivativeBulk'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nComp') - _data_['getSolutionDerivativeParticle'] = ('return', 'drv', 'unitOpId', 'parType', 'time', 'data', 'nTime', 'nParShells', 'nAxialCells', 'nRadialCells', 'nComp') - _data_['getSolutionDerivativeSolid'] = ('return', 'drv', 'unitOpId', 'parType', 'time', 'data', 'nTime', 'nParShells', 'nAxialCells', 'nRadialCells', 'nBound') - _data_['getSolutionDerivativeFlux'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nComp') - _data_['getSolutionDerivativeVolume'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime') - _data_['getSensitivityInlet'] = ('return', 'drv', 'unitOpId', 'idx', 'time', 'data', 'nTime', 'nPort', 'nComp') - _data_['getSensitivityOutlet'] = ('return', 'drv', 'unitOpId', 'idx', 'time', 'data', 'nTime', 'nPort', 'nComp') - _data_['getSensitivityBulk'] = ('return', 'drv', 'unitOpId', 'idx', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nComp') - _data_['getSensitivityParticle'] = ('return', 'drv', 'unitOpId', 'idx', 'parType', 'time', 'data', 'nTime', 'nParShells', 'nAxialCells', 'nRadialCells', 'nComp') - _data_['getSensitivitySolid'] = ('return', 'drv', 'unitOpId', 'idx', 'parType', 'time', 'data', 'nTime', 'nParShells', 'nAxialCells', 'nRadialCells', 'nBound') - _data_['getSensitivityFlux'] = ('return', 'drv', 'unitOpId', 'idx', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nComp') - _data_['getSensitivityVolume'] = ('return', 'drv', 'unitOpId', 'idx', 'time', 'data', 'nTime') - _data_['getSensitivityDerivativeInlet'] = ('return', 'drv', 'unitOpId', 'idx', 'time', 'data', 'nTime', 'nPort', 'nComp') - _data_['getSensitivityDerivativeOutlet'] = ('return', 'drv', 'unitOpId', 'idx', 'time', 'data', 'nTime', 'nPort', 'nComp') - _data_['getSensitivityDerivativeBulk'] = ('return', 'drv', 'unitOpId', 'idx', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nComp') - _data_['getSensitivityDerivativeParticle'] = ('return', 'drv', 'unitOpId', 'idx', 'parType', 'time', 'data', 'nTime', 'nParShells', 'nAxialCells', 'nRadialCells', 'nComp') - _data_['getSensitivityDerivativeSolid'] = ('return', 'drv', 'unitOpId', 'idx', 'parType', 'time', 'data', 'nTime', 'nParShells', 'nAxialCells', 'nRadialCells', 'nBound') - _data_['getSensitivityDerivativeFlux'] = ('return', 'drv', 'unitOpId', 'idx', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nComp') - _data_['getSensitivityDerivativeVolume'] = ('return', 'drv', 'unitOpId', 'idx', 'time', 'data', 'nTime') + """ + Definition of CADET-C-API v1.0 + + signatures : dict with signatures of exported API functions. (See CADET/include/cadet/cadet.h) + lookup_prototype : ctypes for common parameters + lookup_output_argument_type : ctypes for API output parameters (e.g., double* time or double** data) + + """ + + # Order is important, has to match the cdtAPIv010000 struct of the C-API + signatures = {} + signatures['createDriver'] = ('drv',) + signatures['deleteDriver'] = (None, 'drv') + signatures['runSimulation'] = ('return', 'drv', 'parameterProvider') + + signatures['getNumParTypes'] = ('return', 'drv', 'unitOpId', 'nParTypes') + signatures['getNumSensitivities'] = ('return', 'drv', 'nSens') + + signatures['getSolutionInlet'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nPort', 'nComp') + signatures['getSolutionOutlet'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nPort', 'nComp') + signatures['getSolutionBulk'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nComp') + signatures['getSolutionParticle'] = ('return', 'drv', 'unitOpId', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nComp') + signatures['getSolutionSolid'] = ('return', 'drv', 'unitOpId', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nBound') + signatures['getSolutionFlux'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParTypes', 'nComp') + signatures['getSolutionVolume'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime') + + signatures['getSolutionDerivativeInlet'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nPort', 'nComp') + signatures['getSolutionDerivativeOutlet'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nPort', 'nComp') + signatures['getSolutionDerivativeBulk'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nComp') + signatures['getSolutionDerivativeParticle'] = ('return', 'drv', 'unitOpId', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nComp') + signatures['getSolutionDerivativeSolid'] = ('return', 'drv', 'unitOpId', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nBound') + signatures['getSolutionDerivativeFlux'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParTypes', 'nComp') + signatures['getSolutionDerivativeVolume'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime') + + signatures['getSensitivityInlet'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'time', 'data', 'nTime', 'nPort', 'nComp') + signatures['getSensitivityOutlet'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'time', 'data', 'nTime', 'nPort', 'nComp') + signatures['getSensitivityBulk'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nComp') + signatures['getSensitivityParticle'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nComp') + signatures['getSensitivitySolid'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nBound') + signatures['getSensitivityFlux'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParTypes', 'nComp') + signatures['getSensitivityVolume'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'time', 'data', 'nTime') + + signatures['getSensitivityDerivativeInlet'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'time', 'data', 'nTime', 'nPort', 'nComp') + signatures['getSensitivityDerivativeOutlet'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'time', 'data', 'nTime', 'nPort', 'nComp') + signatures['getSensitivityDerivativeBulk'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nComp') + signatures['getSensitivityDerivativeParticle'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nComp') + signatures['getSensitivityDerivativeSolid'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nBound') + signatures['getSensitivityDerivativeFlux'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParTypes', 'nComp') + signatures['getSensitivityDerivativeVolume'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'time', 'data', 'nTime') + + signatures['getLastState'] = ('return', 'drv', 'state', 'nStates') + signatures['getLastStateTimeDerivative'] = ('return', 'drv', 'state', 'nStates') + signatures['getLastUnitState'] = ('return', 'drv', 'unitOpId', 'state', 'nStates') + signatures['getLastUnitStateTimeDerivative'] = ('return', 'drv', 'unitOpId', 'state', 'nStates') + signatures['getLastSensitivityState'] = ('return', 'drv', 'sensIdx', 'state', 'nStates') + signatures['getLastSensitivityStateTimeDerivative'] = ('return', 'drv', 'sensIdx', 'state', 'nStates') + signatures['getLastSensitivityUnitState'] = ('return', 'drv', 'sensIdx', 'unitOpId', 'state', 'nStates') + signatures['getLastSensitivityUnitStateTimeDerivative'] = ('return', 'drv', 'sensIdx', 'unitOpId', 'state', 'nStates') + + signatures['getPrimaryCoordinates'] = ('return', 'drv', 'unitOpId', 'coords', 'nCoords') + signatures['getSecondaryCoordinates'] = ('return', 'drv', 'unitOpId', 'coords', 'nCoords') + signatures['getParticleCoordinates'] = ('return', 'drv', 'unitOpId', 'parType', 'coords', 'nCoords') + + signatures['getSolutionTimes'] = ('return', 'drv', 'time', 'nTime') lookup_prototype = { 'return': c_cadet_result, - 'drv': ctypes.c_void_p, + 'drv': CadetDriver, 'unitOpId': ctypes.c_int, - 'idx': ctypes.c_int, + 'sensIdx': ctypes.c_int, 'parType': ctypes.c_int, 'time': array_double, 'data': array_double, + 'coords': array_double, 'nTime': point_int, + 'nCoords': point_int, 'nPort': point_int, 'nAxialCells': point_int, 'nRadialCells': point_int, 'nParTypes': point_int, + 'nSens': point_int, 'nParShells': point_int, 'nComp': point_int, 'nBound': point_int, + 'state': array_double, + 'nStates': point_int, None: None, - 'parameterProvider': ctypes.POINTER(cadet_dll_parameterprovider.PARAMETERPROVIDER) + 'parameterProvider': ctypes.POINTER(cadet_dll_parameterprovider.PARAMETERPROVIDER), } - lookup_call = { + lookup_output_argument_type = { 'time': ctypes.POINTER(ctypes.c_double), 'nTime': ctypes.c_int, 'data': ctypes.POINTER(ctypes.c_double), @@ -87,16 +128,22 @@ class CADETAPIV010000_DATA(): 'nAxialCells': ctypes.c_int, 'nRadialCells': ctypes.c_int, 'nParTypes': ctypes.c_int, + 'nSens': ctypes.c_int, 'nParShells': ctypes.c_int, 'nComp': ctypes.c_int, 'nBound': ctypes.c_int, + 'state': ctypes.POINTER(ctypes.c_double), + 'nStates': ctypes.c_int, + 'coords': ctypes.POINTER(ctypes.c_double), + 'nCoords': ctypes.c_int, } -def setup_api(): +def _setup_api(): + """list: Tuples with function names and ctype functions""" _fields_ = [] - for key, value in CADETAPIV010000_DATA._data_.items(): + for key, value in CADETAPIV010000_DATA.signatures.items(): args = tuple(CADETAPIV010000_DATA.lookup_prototype[key] for key in value) _fields_.append( (key, ctypes.CFUNCTYPE(*args)) ) @@ -104,141 +151,478 @@ def setup_api(): class CADETAPIV010000(ctypes.Structure): - _fields_ = setup_api() + """Mimic cdtAPIv010000 struct of CADET C-API in ctypes""" + _fields_ = _setup_api() -def null(obj): - "do nothing" - return obj class SimulationResult: - def __init__(self, api, driver): + def __init__(self, api: CADETAPIV010000, driver: CadetDriver): self._api = api self._driver = driver - def load_data(self, unit, get_solution, get_solution_str, idx=None, parType=None, own_data=True): - vars = {} - wrappers = {} - for key in CADETAPIV010000_DATA._data_[get_solution_str]: + def _load_data( + self, + get_solution_str: str, + unitOpId: int | None = None, + sensIdx: int = None, + parType: int = None): + + get_solution = getattr(self._api, get_solution_str) + + # Collect actual values + call_args = [] + call_outputs = {} + + # Construct API call function arguments + for key in CADETAPIV010000_DATA.signatures[get_solution_str]: if key == 'return': + # Skip, this is the actual return value of the API function continue elif key == 'drv': - vars['drv'] = self._driver - wrappers[key] = null - elif key == 'unitOpId': - vars['unitOpId'] = unit - wrappers[key] = null - elif key == 'idx': - vars['idx'] = idx - wrappers[key] = null + call_args.append(self._driver) + elif key == 'unitOpId' and unitOpId is not None: + call_args.append(unitOpId) + elif key == 'sensIdx': + call_args.append(sensIdx) elif key == 'parType': - vars['parType'] = parType - wrappers[key] = null + call_args.append(parType) else: - vars[key] = CADETAPIV010000_DATA.lookup_call[key]() - wrappers[key] = ctypes.byref + _obj = CADETAPIV010000_DATA.lookup_output_argument_type[key]() + call_outputs[key] = _obj + call_args.append(ctypes.byref(_obj)) + + result = get_solution(*call_args) + + if result == _CDT_DATA_NOT_STORED: + # Call successful, but data is not available + return + elif result == _CDT_ERROR_INVALID_INPUTS: + raise ValueError("Error reading data: Invalid call arguments") + elif result != _CDT_OK: + # Something else failed + raise Exception("Error reading data.") + + return call_outputs + + def _process_array( + self, + call_outputs, + data_key: str, + len_key: str, + own_data: bool = True): + + array_length = call_outputs[len_key].value + if array_length == 0: + return + + data = numpy.ctypeslib.as_array(call_outputs[data_key], shape=(array_length, )) + if own_data: + return data.copy() + return data - result = get_solution(*tuple(wrappers[var_key](var_value) for var_key, var_value in vars.items())) + def _process_data( + self, + call_outputs, + own_data: bool = True): shape = [] dims = [] - dimensions = ['nTime', 'nPort', 'nParShells', 'nAxialCells', 'nRadialCells', 'nComp', 'nBound'] + + # Ordering of multi-dimensional arrays, all possible dimensions: + # bulk: 'nTime', ('nAxialCells',) ('nRadialCells' / 'nPorts',) 'nComp' + # particle_liquid: 'nTime', ('nParTypes',) ('nAxialCells',) ('nRadialCells' / 'nPorts',) ('nParShells',) 'nComp' + # particle_solid: 'nTime', ('nParTypes',) ('nAxialCells',) ('nRadialCells' / 'nPorts',) ('nParShells',) 'nComp', 'nBound' + # flux: 'nTime', ('nParTypes',) ('nAxialCells',) ('nRadialCells' / 'nPorts',) 'nComp' + dimensions = [ + 'nTime', + 'nParTypes', + 'nAxialCells', + 'nPort', + 'nRadialCells', + 'nParShells', + 'nComp', + 'nBound', + ] for dim in dimensions: - if dim in vars and vars[dim].value: - shape.append(vars[dim].value) + if dim in call_outputs and call_outputs[dim].value: + shape.append(call_outputs[dim].value) dims.append(dim) - data = numpy.ctypeslib.as_array(vars['data'], shape=shape) - time = numpy.ctypeslib.as_array(vars['time'], shape=(vars['nTime'].value, )) + if len(shape) == 0: + return - if own_data: - return time.copy(), data.copy(), dims - else: - return time, data, dims + if 'data' in call_outputs: + if 'nParShells' in dims: + nParShells = call_outputs['nParShells'].value + if nParShells == 1: + shape.pop(dims.index('nParShells')) + + data = numpy.ctypeslib.as_array(call_outputs['data'], shape=shape) + if own_data: + data = data.copy() + + if 'time' in call_outputs: + n_time = call_outputs['nTime'].value + time = numpy.ctypeslib.as_array(call_outputs['time'], shape=(n_time, )) + if own_data: + time = time.copy() + + return time, data, dims + + def _load_and_process(self, *args, own_data=True, **kwargs): + call_outputs = self._load_data(*args, **kwargs) + + if call_outputs is None: + return + + processed_results = self._process_data(call_outputs, own_data) + return processed_results + + def _load_and_process_array(self, data_key: str, len_key: str, *args, own_data=True, **kwargs): + call_outputs = self._load_data(*args, **kwargs) + + if call_outputs is None: + return None + + return self._process_array(call_outputs, data_key, len_key, own_data) + + def npartypes(self, unitOpId: int): + call_outputs = self._load_data( + 'getNumParTypes', + unitOpId=unitOpId, + ) + + return int(call_outputs['nParTypes'].value) + + def nsensitivities(self): + call_outputs = self._load_data( + 'getNumSensitivities', + ) + + return int(call_outputs['nSens'].value) + + def solution_inlet(self, unitOpId: int, own_data=True): + return self._load_and_process( + 'getSolutionInlet', + unitOpId=unitOpId, + own_data=own_data, + ) + + def solution_outlet(self, unitOpId: int, own_data=True): + return self._load_and_process( + 'getSolutionOutlet', + unitOpId=unitOpId, + own_data=own_data, + ) + + def solution_bulk(self, unitOpId: int, own_data=True): + return self._load_and_process( + 'getSolutionBulk', + unitOpId=unitOpId, + own_data=own_data, + ) + + def solution_particle(self, unitOpId: int, parType: int, own_data=True): + return self._load_and_process( + 'getSolutionParticle', + unitOpId=unitOpId, + parType=parType, + own_data=own_data, + ) + + def solution_solid(self, unitOpId: int, parType: int, own_data=True): + return self._load_and_process( + 'getSolutionSolid', + unitOpId=unitOpId, + parType=parType, + own_data=own_data, + ) + + def solution_flux(self, unitOpId: int, own_data=True): + return self._load_and_process( + 'getSolutionFlux', + unitOpId=unitOpId, + own_data=own_data, + ) + + def solution_volume(self, unitOpId: int, own_data=True): + return self._load_and_process( + 'getSolutionVolume', + unitOpId=unitOpId, + own_data=own_data, + ) + + def soldot_inlet(self, unitOpId: int, own_data=True): + return self._load_and_process( + 'getSolutionDerivativeInlet', + unitOpId=unitOpId, + own_data=own_data, + ) + + def soldot_outlet(self, unitOpId: int, own_data=True): + return self._load_and_process( + 'getSolutionDerivativeOutlet', + unitOpId=unitOpId, + own_data=own_data, + ) + + def soldot_bulk(self, unitOpId: int, own_data=True): + return self._load_and_process( + 'getSolutionDerivativeBulk', + unitOpId=unitOpId, + own_data=own_data, + ) + + def soldot_particle(self, unitOpId: int, parType: int, own_data=True): + return self._load_and_process( + 'getSolutionDerivativeParticle', + unitOpId=unitOpId, + parType=parType, + own_data=own_data, + ) - def inlet(self, unit, own_data=True): - return self.load_data(unit, self._api.getSolutionInlet, 'getSolutionInlet', own_data=own_data) + def soldot_solid(self, unitOpId: int, parType: int, own_data=True): + return self._load_and_process( + 'getSolutionDerivativeSolid', + unitOpId=unitOpId, + parType=parType, + own_data=own_data, + ) - def outlet(self, unit, own_data=True): - return self.load_data(unit, self._api.getSolutionOutlet, 'getSolutionOutlet', own_data=own_data) + def soldot_flux(self, unitOpId: int, own_data=True): + return self._load_and_process( + 'getSolutionDerivativeFlux', + unitOpId=unitOpId, + own_data=own_data, + ) - def bulk(self, unit, own_data=True): - return self.load_data(unit, self._api.getSolutionBulk, 'getSolutionBulk', own_data=own_data) + def soldot_volume(self, unitOpId: int, own_data=True): + return self._load_and_process( + 'getSolutionDerivativeVolume', + unitOpId=unitOpId, + own_data=own_data, + ) - def particle(self, unit, parType, own_data=True): - return self.load_data(unit, self._api.getSolutionBulk, 'getSolutionBulk', own_data=own_data) + def sens_inlet(self, unitOpId: int, sensIdx: int, own_data=True): + return self._load_and_process( + 'getSensitivityInlet', + unitOpId=unitOpId, + sensIdx=sensIdx, + own_data=own_data, + ) - def solid(self, unit, parType, own_data=True): - return self.load_data(unit, self._api.getSolutionSolid, 'getSolutionSolid', own_data=own_data) + def sens_outlet(self, unitOpId: int, sensIdx: int, own_data=True): + return self._load_and_process( + 'getSensitivityOutlet', + unitOpId=unitOpId, + sensIdx=sensIdx, + own_data=own_data, + ) - def flux(self, unit, own_data=True): - return self.load_data(unit, self._api.getSolutionFlux, 'getSolutionFlux', own_data=own_data) + def sens_bulk(self, unitOpId: int, sensIdx: int, own_data=True): + return self._load_and_process( + 'getSensitivityBulk', + unitOpId=unitOpId, + sensIdx=sensIdx, + own_data=own_data, + ) - def volume(self, unit, own_data=True): - return self.load_data(unit, self._api.getSolutionVolume, 'getSolutionVolume', own_data=own_data) + def sens_particle(self, unitOpId: int, sensIdx: int, parType: int, own_data=True): + return self._load_and_process( + 'getSensitivityParticle', + unitOpId=unitOpId, + sensIdx=sensIdx, + parType=parType, + own_data=own_data, + ) - def derivativeInlet(self, unit, own_data=True): - return self.load_data(unit, self._api.getSolutionDerivativeInlet, 'getSolutionDerivativeInlet', own_data=own_data) + def sens_solid(self, unitOpId: int, sensIdx: int, parType: int, own_data=True): + return self._load_and_process( + 'getSensitivitySolid', + unitOpId=unitOpId, + sensIdx=sensIdx, + parType=parType, + own_data=own_data, + ) - def derivativeOutlet(self, unit, own_data=True): - return self.load_data(unit, self._api.getSolutionDerivativeOutlet, 'getSolutionDerivativeOutlet', own_data=own_data) + def sens_flux(self, unitOpId: int, sensIdx: int, own_data=True): + return self._load_and_process( + 'getSensitivityFlux', + unitOpId=unitOpId, + sensIdx=sensIdx, + own_data=own_data, + ) - def derivativeBulk(self, unit, own_data=True): - return self.load_data(unit, self._api.getSolutionDerivativeBulk, 'getSolutionDerivativeBulk', own_data=own_data) + def sens_volume(self, unitOpId: int, sensIdx: int, own_data=True): + return self._load_and_process( + 'getSensitivityVolume', + unitOpId=unitOpId, + sensIdx=sensIdx, + own_data=own_data, + ) - def derivativeParticle(self, unit, parType, own_data=True): - return self.load_data(unit, self._api.getSolutionDerivativeParticle, 'getSolutionDerivativeParticle', own_data=own_data) + def sensdot_inlet(self, unitOpId: int, sensIdx: int, own_data=True): + return self._load_and_process( + 'getSensitivityDerivativeInlet', + unitOpId=unitOpId, + sensIdx=sensIdx, + own_data=own_data, + ) - def derivativeSolid(self, unit, parType, own_data=True): - return self.load_data(unit, self._api.getSolutionDerivativeSolid, 'getSolutionDerivativeSolid', own_data=own_data) + def sensdot_outlet(self, unitOpId: int, sensIdx: int, own_data=True): + return self._load_and_process( + 'getSensitivityDerivativeOutlet', + unitOpId, + sensIdx=sensIdx, + own_data=own_data, + ) - def derivativeFlux(self, unit, own_data=True): - return self.load_data(unit, self._api.getSolutionDerivativeFlux, 'getSolutionDerivativeFlux', own_data=own_data) + def sensdot_bulk(self, unitOpId: int, sensIdx: int, own_data=True): + return self._load_and_process( + 'getSensitivityDerivativeBulk', + unitOpId=unitOpId, + sensIdx=sensIdx, + own_data=own_data, + ) - def derivativeVolume(self, unit, own_data=True): - return self.load_data(unit, self._api.getSolutionDerivativeVolume, 'getSolutionDerivativeVolume', own_data=own_data) + def sensdot_particle(self, unitOpId: int, sensIdx: int, parType: int, own_data=True): + return self._load_and_process( + 'getSensitivityDerivativeParticle', + unitOpId=unitOpId, + sensIdx=sensIdx, + parType=parType, + own_data=own_data, + ) - def sensitivityInlet(self, unit, idx, own_data=True): - return self.load_data(unit, self._api.getSensitivityInlet, 'getSensitivityInlet', idx=idx, own_data=own_data) + def sensdot_solid(self, unitOpId: int, sensIdx: int, parType: int, own_data=True): + return self._load_and_process( + 'getSensitivityDerivativeSolid', + unitOpId=unitOpId, + sensIdx=sensIdx, + parType=parType, + own_data=own_data, + ) - def sensitivityOutlet(self, unit, idx, own_data=True): - return self.load_data(unit, self._api.getSensitivityOutlet, 'getSensitivityOutlet', idx=idx, own_data=own_data) + def sensdot_flux(self, unitOpId: int, sensIdx: int, own_data=True): + return self._load_and_process( + 'getSensitivityDerivativeFlux', + unitOpId=unitOpId, + sensIdx=sensIdx, + own_data=own_data, + ) - def sensitivityBulk(self, unit, idx, own_data=True): - return self.load_data(unit, self._api.getSensitivityBulk, 'getSensitivityBulk', idx=idx, own_data=own_data) + def sensdot_volume(self, unitOpId: int, sensIdx: int, own_data=True): + return self._load_and_process( + 'getSensitivityDerivativeVolume', + unitOpId=unitOpId, + sensIdx=sensIdx, + own_data=own_data, + ) - def sensitivityParticle(self, unit, idx, parType, own_data=True): - return self.load_data(unit, self._api.getSensitivityParticle, 'getSensitivityParticle', idx=idx, parType=parType, own_data=own_data) + def last_state_y(self, own_data=True): + return self._load_and_process_array( + 'state', + 'nStates', + 'getLastState', + own_data=own_data, + ) - def sensitivitySolid(self, unit, idx, parType, own_data=True): - return self.load_data(unit, self._api.getSensitivitySolid, 'getSensitivitySolid', idx=idx, parType=parType, own_data=own_data) + def last_state_ydot(self, own_data=True): + return self._load_and_process_array( + 'state', + 'nStates', + 'getLastStateTimeDerivative', + own_data=own_data, + ) - def sensitivityFlux(self, unit, idx, own_data=True): - return self.load_data(unit, self._api.getSensitivityFlux, 'getSensitivityFlux', idx=idx, own_data=own_data) + def last_state_y_unit(self, unitOpId: int, own_data=True): + return self._load_and_process_array( + 'state', + 'nStates', + 'getLastUnitState', + unitOpId=unitOpId, + own_data=own_data, + ) + + def last_state_ydot_unit(self, unitOpId: int, own_data=True): + return self._load_and_process_array( + 'state', + 'nStates', + 'getLastUnitStateTimeDerivative', + unitOpId=unitOpId, + own_data=own_data, + ) - def sensitivityVolume(self, unit, idx, own_data=True): - return self.load_data(unit, self._api.getSensitivityVolume, 'getSensitivityVolume', idx=idx, own_data=own_data) + def last_state_sens(self, sensIdx: int, own_data=True): + return self._load_and_process_array( + 'state', + 'nStates', + 'getLastSensitivityState', + sensIdx=sensIdx, + own_data=own_data, + ) - def sensitivityDerivativeInlet(self, unit, idx, own_data=True): - return self.load_data(unit, self._api.getSensitivityDerivativeInlet, 'getSensitivityDerivativeInlet', idx=idx, own_data=own_data) + def last_state_sensdot(self, sensIdx: int, own_data=True): + return self._load_and_process_array( + 'state', + 'nStates', + 'getLastSensitivityStateTimeDerivative', + sensIdx=sensIdx, + own_data=own_data, + ) - def sensitivityDerivativeOutlet(self, unit, idx, own_data=True): - return self.load_data(unit, self._api.getSensitivityDerivativeOutlet, 'getSensitivityDerivativeOutlet', idx=idx, own_data=own_data) + def last_state_sens_unit(self, unitOpId: int, sensIdx: int, own_data=True): + return self._load_and_process( + 'getLastSensitivityUnitState', + unitOpId=unitOpId, + sensIdx=sensIdx, + own_data=own_data, + ) - def sensitivityDerivativeBulk(self, unit, idx, own_data=True): - return self.load_data(unit, self._api.getSensitivityDerivativeBulk, 'getSensitivityDerivativeBulk', idx=idx, own_data=own_data) + def last_state_sensdot_unit(self, unitOpId: int, sensIdx: int, own_data=True): + return self._load_and_process( + 'getLastSensitivityUnitStateTimeDerivative', + unitOpId=unitOpId, + sensIdx=sensIdx, + own_data=own_data, + ) - def sensitivityDerivativeParticle(self, unit, idx, parType, own_data=True): - return self.load_data(unit, self._api.getSensitivityDerivativeParticle, 'getSensitivityDerivativeParticle', idx=idx, parType=parType, own_data=own_data) + def primary_coordinates(self, unitOpId: int, own_data=True): + return self._load_and_process_array( + 'coords', + 'nCoords', + 'getPrimaryCoordinates', + unitOpId=unitOpId, + own_data=own_data, + ) - def sensitivityDerivativeSolid(self, unit, idx, parType, own_data=True): - return self.load_data(unit, self._api.getSensitivityDerivativeSolid, 'getSensitivityDerivativeSolid', idx=idx, parType=parType, own_data=own_data) + def secondary_coordinates(self, unitOpId: int, own_data=True): + return self._load_and_process_array( + 'coords', + 'nCoords', + 'getSecondaryCoordinates', + unitOpId=unitOpId, + own_data=own_data, + ) - def sensitivityDerivativeFlux(self, unit, idx, own_data=True): - return self.load_data(unit, self._api.getSensitivityDerivativeFlux, 'getSensitivityDerivativeFlux', idx=idx, own_data=own_data) + def particle_coordinates(self, unitOpId: int, parType: int, own_data=True): + return self._load_and_process_array( + 'coords', + 'nCoords', + 'getParticleCoordinates', + unitOpId=unitOpId, + parType=parType, + own_data=own_data, + ) - def sensitivityDerivativeVolume(self, unit, idx, own_data=True): - return self.load_data(unit, self._api.getSensitivityDerivativeVolume, 'getSensitivityDerivativeVolume', idx=idx, own_data=own_data) + def solution_times(self, own_data=True): + return self._load_and_process_array( + 'time', + 'nTime', + 'getSolutionTimes', + own_data=own_data + ) class CadetDLL: @@ -295,6 +679,7 @@ def __init__(self, dll_path): cdtGetAPIv010000(ctypes.byref(self._api)) self._driver = self._api.createDriver() + self.res = None def clear(self): if hasattr(self, "res"): @@ -306,205 +691,314 @@ def __del__(self): log_print('deleteDriver()') self._api.deleteDriver(self._driver) - - def run(self, filename = None, simulation=None, timeout = None, check=None): + def run(self, filename=None, simulation=None, timeout=None, check=None): pp = cadet_dll_parameterprovider.PARAMETERPROVIDER(simulation) self._api.runSimulation(self._driver, ctypes.byref(pp)) self.res = SimulationResult(self._api, self._driver) + + # TODO: Return if simulation was successful or crashed return self.res - def load_solution(self, sim, solution_fun, solution_str): - # - [ ] Split Ports (incl `SINGLE_AS_MULTI_PORT`) - # - [ ] Split Partype (Particle + Solid) - # - [ ] Coordinates? - # - [ ] Sensitivities / IDs - # - [ ] LAST_STATE_Y / LAST_STATE_YDOT - # - [ ] LAST_STATE_SENSY_XXX / LAST_STATE_SENSYDOT_XXX + def load_results(self, sim): + if self.res is None: + return + + self.load_solution_times(sim) + self.load_coordinates(sim) + self.load_solution(sim) + self.load_sensitivity(sim) + self.load_state(sim) + + def load_solution_times(self, sim): + """Load solution_times from simulation results.""" + write_solution_times = sim.root.input['return'].get('write_solution_times', 0) + if write_solution_times: + sim.root.output.solution.solution_times = self.res.solution_times() + + def load_coordinates(self, sim): + """Load coordinates data from simulation results.""" + coordinates = addict.Dict() + # TODO: Use n_units from API? + for unit in range(sim.root.input.model.nunits): + unit_index = self._get_index_string('unit', unit) + write_coordinates = sim.root.input['return'][unit_index].get('write_coordinates', 0) + if write_coordinates: + pc = self.res.primary_coordinates(unit) + if pc is not None: + coordinates[unit_index]['axial_coordinates'] = pc + + sc = self.res.secondary_coordinates(unit) + if sc is not None: + coordinates[unit_index]['radial_coordinates'] = sc + + num_par_types = self.res.npartypes(unit) + for pt in range(num_par_types): + par_coords = self.res.particle_coordinates(unit, pt) + if par_coords is not None: + par_idx = self._get_index_string('particle_coordinates', pt) + coordinates[unit_index][par_idx] = par_coords + + if len(coordinates) > 0: + sim.root.output.coordinates = coordinates + + def load_solution(self, sim): + """Load solution data from simulation results.""" solution = addict.Dict() - if self.res is not None: - for key, value in sim.root.input['return'].items(): - if key.startswith('unit'): - if value[f'write_{solution_str}']: - unit = int(key[-3:]) - t, out, dims = solution_fun(unit) + # TODO: Use n_units from API? + for unit in range(sim.root.input.model.nunits): + unit_index = self._get_index_string('unit', unit) + unit_solution = addict.Dict() + + unit_solution.update(self._load_solution_io(sim, unit, 'solution_inlet')) + unit_solution.update(self._load_solution_io(sim, unit, 'solution_outlet')) + unit_solution.update(self._load_solution_trivial(sim, unit, 'solution_bulk')) + unit_solution.update(self._load_solution_particle(sim, unit, 'solution_particle')) + unit_solution.update(self._load_solution_particle(sim, unit, 'solution_solid')) + unit_solution.update(self._load_solution_trivial(sim, unit, 'solution_flux')) + unit_solution.update(self._load_solution_trivial(sim, unit, 'solution_volume')) + + unit_solution.update(self._load_solution_io(sim, unit, 'soldot_inlet')) + unit_solution.update(self._load_solution_io(sim, unit, 'soldot_outlet')) + unit_solution.update(self._load_solution_trivial(sim, unit, 'soldot_bulk')) + unit_solution.update(self._load_solution_particle(sim, unit, 'soldot_particle')) + unit_solution.update(self._load_solution_particle(sim, unit, 'soldot_solid')) + unit_solution.update(self._load_solution_trivial(sim, unit, 'soldot_flux')) + unit_solution.update(self._load_solution_trivial(sim, unit, 'soldot_volume')) + + if len(unit_solution) > 1: + solution[unit_index].update(unit_solution) + + if len(unit_solution) > 1: + sim.root.output.solution.update(solution) - if not len(solution.solution_times): - solution.solution_times = t - - solution[key][solution_str] = out - return solution - - def load_solution_io(self, sim, solution_fun, solution_str): - solution = addict.Dict() - if self.res is not None: - for key,value in sim.root.input['return'].items(): - if key.startswith('unit'): - if value[f'write_{solution_str}']: - unit = int(key[-3:]) - t, out, dims = solution_fun(unit) - - if not len(solution.solution_times): - solution.solution_times = t - - split_components_data = value.get('split_components_data', 1) - split_ports_data = value.get('split_ports_data', 1) - single_as_multi_port = value.get('single_as_multi_port', 0) - - nComp = dims.index('nComp') - try: - nPorts = dims.index('nPorts') - except ValueError: - nPorts = None - - if split_components_data: - if split_ports_data: - if nPorts is None: - if single_as_multi_port: - for comp in range(out.shape[nComp]): - comp_out = numpy.squeeze(out[..., comp]) - solution[key][f'{solution_str}_port_000_comp_{comp:03d}'] = comp_out - else: - for comp in range(out.shape[nComp]): - comp_out = numpy.squeeze(out[..., comp]) - solution[key][f'{solution_str}_comp_{comp:03d}'] = comp_out - else: - for port in range(out.shape[nPorts]): - for comp in range(out.shape[nComp]): - comp_out = numpy.squeeze(out[..., port, comp]) - solution[key][f'{solution_str}_port_{port:03d}_comp_{comp:03d}'] = comp_out - else: - for comp in range(out.shape[nComp]): - comp_out = numpy.squeeze(out[...,comp]) - solution[key][f'{solution_str}_comp_{comp:03d}'] = comp_out - else: - if split_ports_data: - if nPorts is None: - if single_as_multi_port: - solution[key][f'{solution_str}_port_000'] = out - else: - solution[key][solution_str] = out - else: - for port in range(out.shape[nPorts]): - port_out = numpy.squeeze(out[..., port, :]) - solution[key][f'{solution_str}_port_{port:03d}'] = port_out - else: - solution[key][solution_str] = out return solution - def load_inlet(self, sim): - return self.load_solution_io(sim, self.res.inlet, 'solution_inlet') - - def load_outlet(self, sim): - return self.load_solution_io(sim, self.res.outlet, 'solution_outlet') - - def load_bulk(self, sim): - return self.load_solution(sim, self.res.bulk, 'solution_bulk') - - def load_particle(self, sim): - return self.load_solution(sim, self.res.particle, 'solution_particle') - - def load_solid(self, sim): - return self.load_solution(sim, self.res.solid, 'solution_solid') - - def load_flux(self, sim): - return self.load_solution(sim, self.res.flux, 'solution_flux') - - def load_flux(self, sim): - return self.load_solution(sim, self.res.flux, 'solution_flux') - - def load_volume(self, sim): - return self.load_solution(sim, self.res.volume, 'solution_volume') - - def load_derivative_inlet(self, sim): - return self.load_solution_io(sim, self.res.derivativeInlet, 'soldot_inlet') - - def load_derivative_outlet(self, sim): - return self.load_solution_io(sim, self.res.derivativeOutlet, 'soldot_outlet') - - def load_derivative_bulk(self, sim): - return self.load_solution(sim, self.res.derivativeBulk, 'soldot_bulk') - - def load_derivative_particle(self, sim): - return self.load_solution(sim, self.res.derivativeParticle, 'soldot_particle') - - def load_derivative_solid(self, sim): - return self.load_solution(sim, self.res.derivativeSolid, 'soldot_solid') + def load_sensitivity(self, sim): + """Load sensitivity data from simulation results.""" + sensitivity = addict.Dict() + nsens = sim.root.input.sensitivity.get('nsens', 0) + + for sens in range(nsens): + sens_index = self._get_index_string('param', sens) + + for unit in range(sim.root.input.model.nunits): + unit_sensitivity = addict.Dict() + unit_index = self._get_index_string('unit', unit) + + unit_sensitivity.update( + self._load_solution_io(sim, unit, 'sens_inlet', sens) + ) + unit_sensitivity.update( + self._load_solution_io(sim, unit, 'sens_outlet', sens) + ) + unit_sensitivity.update( + self._load_solution_trivial(sim, unit, 'sens_bulk', sens) + ) + unit_sensitivity.update( + self._load_solution_particle(sim, unit, 'sens_particle', sens) + ) + unit_sensitivity.update( + self._load_solution_particle(sim, unit, 'sens_solid', sens) + ) + unit_sensitivity.update( + self._load_solution_trivial(sim, unit, 'sens_flux', sens) + ) + unit_sensitivity.update( + self._load_solution_trivial(sim, unit, 'sens_volume', sens) + ) + + unit_sensitivity.update( + self._load_solution_io(sim, unit, 'sensdot_inlet', sens) + ) + unit_sensitivity.update( + self._load_solution_io(sim, unit, 'sensdot_outlet', sens) + ) + unit_sensitivity.update( + self._load_solution_trivial(sim, unit, 'sensdot_bulk', sens) + ) + unit_sensitivity.update( + self._load_solution_particle(sim, unit, 'sensdot_particle', sens) + ) + unit_sensitivity.update( + self._load_solution_particle(sim, unit, 'sensdot_solid', sens) + ) + unit_sensitivity.update( + self._load_solution_trivial(sim, unit, 'sensdot_flux', sens) + ) + unit_sensitivity.update( + self._load_solution_trivial(sim, unit, 'sensdot_volume', sens) + ) + + if len(unit_sensitivity) > 0: + sensitivity[sens_index][unit_index] = unit_sensitivity + + if len(sensitivity) > 0: + sim.root.output.sensitivity = sensitivity + + def load_state(self, sim): + """Load last state from simulation results.""" + # System state + write_solution_last = sim.root.input['return'].get('write_solution_last', 0) + if write_solution_last: + sim.root.output['last_state_y'] = self.res.last_state_y() + sim.root.output['last_state_ydot'] = self.res.last_state_ydot() + + # System sensitivities + write_sens_last = sim.root.input['return'].get('write_sens_last', 0) + if write_sens_last: + for idx in range(self.res.nsensitivities()): + idx_str_y = self._get_index_string('last_state_sensy', idx) + sim.root.output[idx_str_y] = self.res.last_state_sens(idx) + idx_str_ydot = self._get_index_string('last_state_sensydot', idx) + sim.root.output[idx_str_ydot] = self.res.last_state_sensdot(idx) + + + # TODO: Use n_units from API? + solution = sim.root.output.solution + for unit in range(sim.root.input.model.nunits): + unit_index = self._get_index_string('unit', unit) + write_solution_last = sim.root.input['return'][unit_index].get('write_solution_last_unit', 0) + if write_solution_last: + solution_last_unit = self.res.last_state_y_unit(unit) + solution[unit_index]['last_state_y'] = solution_last_unit + soldot_last_unit = self.res.last_state_ydot_unit(unit) + solution[unit_index]['last_state_ydot'] = soldot_last_unit + + def _checks_if_write_is_true(func): + """Decorator to check if unit operation solution should be written out.""" + def wrapper(self, sim, unitOpId, solution_str, *args, **kwargs): + unit_index = self._get_index_string('unit', unitOpId) + + write = sim.root.input['return'][unit_index].get(f'write_{solution_str}', 0) + + if not write: + return {} + + solution = func(self, sim, unitOpId, solution_str, *args, **kwargs) + + if solution is None: + return {} + + return solution + + return wrapper + + def _loads_data(func): + """Decorator to load data from simulation results before processing further.""" + def wrapper(self, sim, unitOpId, solution_str, sensIdx=None, *args, **kwargs): + solution_fun = getattr(self.res, solution_str) + if sensIdx is None: + data = solution_fun(unitOpId) + else: + data = solution_fun(unitOpId, sensIdx) - def load_derivative_flux(self, sim): - return self.load_solution(sim, self.res.derivativeFlux, 'soldot_flux') + if data is None: + return - def load_derivative_volume(self, sim): - return self.load_solution(sim, self.res.derivativeVolume, 'soldot_volume') + solution = func(self, sim, data, unitOpId, solution_str, *args, **kwargs) - def load_sensitivity_inlet(self, sim): - return self.load_solution_io(sim, self.res.sensitivityInlet, 'sens_inlet') + return solution - def load_sensitivity_outlet(self, sim): - return self.load_solution_io(sim, self.res.sensitivityOutlet, 'sens_outlet') + return wrapper - def load_sensitivity_bulk(self, sim): - return self.load_solution(sim, self.res.sensitivityBulk, 'sens_bulk') + @_checks_if_write_is_true + @_loads_data + def _load_solution_trivial(self, sim, data, unitOpId, solution_str, sensIdx=None): + solution = addict.Dict() + _, out, _ = data + solution[solution_str] = out - def load_sensitivity_particle(self, sim): - return self.load_solution(sim, self.res.sensitivityParticle, 'sens_particle') + return solution - def load_sensitivity_solid(self, sim): - return self.load_solution(sim, self.res.sensitivitySolid, 'sens_solid') + @_checks_if_write_is_true + @_loads_data + def _load_solution_io(self, sim, data, unitOpId, solution_str, sensIdx=None): + solution = addict.Dict() + _, out, dims = data + + split_components_data = sim.root.input['return'].get('split_components_data', 1) + split_ports_data = sim.root.input['return'].get('split_ports_data', 1) + single_as_multi_port = sim.root.input['return'].get('single_as_multi_port', 0) + + nComp_idx = dims.index('nComp') + nComp = out.shape[nComp_idx] + try: + nPort_idx = dims.index('nPort') + nPort = out.shape[nPort_idx] + except ValueError: + nPort_idx = None + nPort = 1 + + if split_components_data: + if split_ports_data: + if nPort == 1: + if single_as_multi_port: + for comp in range(nComp): + comp_out = numpy.squeeze(out[..., comp]) + solution[f'{solution_str}_port_000_comp_{comp:03d}'] = comp_out + else: + for comp in range(nComp): + comp_out = numpy.squeeze(out[..., comp]) + solution[f'{solution_str}_comp_{comp:03d}'] = comp_out + else: + for port in range(nPort): + for comp in range(nComp): + comp_out = numpy.squeeze(out[..., port, comp]) + solution[f'{solution_str}_port_{port:03d}_comp_{comp:03d}'] = comp_out + else: + for comp in range(nComp): + comp_out = numpy.squeeze(out[..., comp]) + solution[f'{solution_str}_comp_{comp:03d}'] = comp_out + else: + if split_ports_data: + if nPort == 1: + if single_as_multi_port: + solution[f'{solution_str}_port_000'] = out + else: + solution[solution_str] = numpy.squeeze(out[..., 0, :]) + else: + for port in range(nPort): + port_out = numpy.squeeze(out[..., port, :]) + solution[f'{solution_str}_port_{port:03d}'] = port_out + else: + if nPort == 1 and nPort_idx is not None: + solution[solution_str] = numpy.squeeze(out[..., 0, :]) + else: + solution[solution_str] = out - def load_sensitivity_flux(self, sim): - return self.load_solution(sim, self.res.sensitivityFlux, 'sens_flux') + return solution - def load_sensitivity_volume(self, sim): - return self.load_solution(sim, self.res.sensitivityVolume, 'sens_volume') + @_checks_if_write_is_true + def _load_solution_particle(self, sim, unitOpId, solution_str, sensIdx=None): + solution = addict.Dict() + solution_fun = getattr(self.res, solution_str) - def load_sensitivity_derivative_inlet(self, sim): - return self.load_solution_io(sim, self.res.sensitivityDerivativeInlet, 'sensdot_inlet') + npartype = self.res.npartypes(unitOpId) - def load_sensitivity_derivative_outlet(self, sim): - return self.load_solution_io(sim, self.res.sensitivityDerivativeOutlet, 'sensdot_outlet') + for partype in range(npartype): + if sensIdx is None: + data = solution_fun(unitOpId, partype) + else: + data = solution_fun(unitOpId, sensIdx, partype) - def load_sensitivity_derivative_bulk(self, sim): - return self.load_solution(sim, self.res.sensitivityDerivativeBulk, 'sensdot_bulk') + if data is None: + continue - def load_sensitivity_derivative_particle(self, sim): - return self.load_solution(sim, self.res.sensitivityDerivativeParticle, 'sensdot_particle') + _, out, dims = data - def load_sensitivity_derivative_solid(self, sim): - return self.load_solution(sim, self.res.sensitivityDerivativeSolid, 'sensdot_solid') + if npartype == 1: + solution[solution_str] = out + else: + par_index = self._get_index_string('partype', partype) + solution[f'{solution_str}_{par_index}'] = out - def load_sensitivity_derivative_flux(self, sim): - return self.load_solution(sim, self.res.sensitivityDerivativeFlux, 'sensdot_flux') + if len(solution) == 0: + return - def load_sensitivity_derivative_volume(self, sim): - return self.load_solution(sim, self.res.sensitivityDerivativeVolume, 'sensdot_volume') + return solution - def load_results(self, sim): - sim.root.output.solution.update(self.load_inlet(sim)) - sim.root.output.solution.update(self.load_outlet(sim)) - sim.root.output.solution.update(self.load_bulk(sim)) - sim.root.output.solution.update(self.load_particle(sim)) - sim.root.output.solution.update(self.load_solid(sim)) - sim.root.output.solution.update(self.load_flux(sim)) - sim.root.output.solution.update(self.load_volume(sim)) - sim.root.output.solution.update(self.load_derivative_inlet(sim)) - sim.root.output.solution.update(self.load_derivative_outlet(sim)) - sim.root.output.solution.update(self.load_derivative_bulk(sim)) - sim.root.output.solution.update(self.load_derivative_particle(sim)) - sim.root.output.solution.update(self.load_derivative_solid(sim)) - sim.root.output.solution.update(self.load_derivative_flux(sim)) - sim.root.output.solution.update(self.load_derivative_volume(sim)) - #sim.root.output.solution.update(self.load_sensitivity_inlet(sim)) - #sim.root.output.solution.update(self.load_sensitivity_outlet(sim)) - #sim.root.output.solution.update(self.load_sensitivity_bulk(sim)) - #sim.root.output.solution.update(self.load_sensitivity_particle(sim)) - #sim.root.output.solution.update(self.load_sensitivity_solid(sim)) - #sim.root.output.solution.update(self.load_sensitivity_flux(sim)) - #sim.root.output.solution.update(self.load_sensitivity_volume(sim)) - #sim.root.output.solution.update(self.load_sensitivity_derivative_inlet(sim)) - #sim.root.output.solution.update(self.load_sensitivity_derivative_outlet(sim)) - #sim.root.output.solution.update(self.load_sensitivity_derivative_bulk(sim)) - #sim.root.output.solution.update(self.load_sensitivity_derivative_particle(sim)) - #sim.root.output.solution.update(self.load_sensitivity_derivative_solid(sim)) - #sim.root.output.solution.update(self.load_sensitivity_derivative_flux(sim)) - #sim.root.output.solution.update(self.load_sensitivity_derivative_volume(sim)) + @staticmethod + def _get_index_string(prefix, index): + """Helper method to get string indices (e.g. (unit, 0) -> 'unit_000').""" + return f'{prefix}_{index:03d}' diff --git a/cadet/cadet_dll_utils.py b/cadet/cadet_dll_utils.py index d4b4045..b2c8a99 100644 --- a/cadet/cadet_dll_utils.py +++ b/cadet/cadet_dll_utils.py @@ -120,7 +120,7 @@ def param_provider_get_int_array(reader, name, n_elem, val): o = c[n] if isinstance(o, list): o = numpy.ascontiguousarray(o) - if (not isinstance(o, numpy.ndarray)) or (o.dtype != numpy.int) or (not o.flags.c_contiguous): + if (not isinstance(o, numpy.ndarray)) or (o.dtype != int) or (not o.flags.c_contiguous): return -1 n_elem[0] = ctypes.c_int(o.size) @@ -189,25 +189,20 @@ def param_provider_get_bool_array_item(reader, name, index, val): def param_provider_get_string_array_item(reader, name, index, val): - n = name.decode('utf-8') - c = reader.current() + name = name.decode('utf-8') + current_reader = reader.current() - if n in c: - o = c[n] - if index == 0: - if hasattr(o, 'encode'): - bytes_val = o.encode('utf-8') - elif hasattr(o, 'decode'): - bytes_val = o + if name in current_reader: + str_value = current_reader[name] + if isinstance(str_value, bytes): + bytes_val = str_value else: - if hasattr(o[index], 'encode'): - bytes_val = o[index].encode('utf-8') - elif hasattr(o[index], 'decode'): - bytes_val = o[index] + bytes_val = str_value[index] reader.buffer = bytes_val val[0] = ctypes.cast(reader.buffer, ctypes.c_char_p) - log_print('GET array [string] ({}) {}: {}'.format(index, n, reader.buffer.decode('utf-8'))) + log_print(f"GET array [string] ({index}) {name}: {reader.buffer.decode('utf-8')}") + return 0 return -1 @@ -264,7 +259,7 @@ def param_provider_num_elements(reader, name): def param_provider_push_scope(reader, name): n = name.decode('utf-8') - + if reader.push_scope(n): return 0 else: @@ -273,4 +268,4 @@ def param_provider_push_scope(reader, name): def param_provider_pop_scope(reader): reader.pop_scope() - return 0 \ No newline at end of file + return 0 diff --git a/tests/__init__.py b/tests/__init__.py index e69de29..633f866 100644 --- a/tests/__init__.py +++ b/tests/__init__.py @@ -0,0 +1,2 @@ +# -*- coding: utf-8 -*- + diff --git a/tests/test_dll.py b/tests/test_dll.py new file mode 100644 index 0000000..871173f --- /dev/null +++ b/tests/test_dll.py @@ -0,0 +1,935 @@ +import copy +import os +from pathlib import Path +import platform +import subprocess +import sys + +from addict import Dict +import pytest + +from cadet import Cadet + + +# %% Utility methods + +# TODO: Remove once #14 is merged +cadet_root = Path('/home/jo/code/CADET/install/capi/') + +def setup_model( + cadet_root, + use_dll=True, + model='GENERAL_RATE_MODEL', + n_partypes=1, + include_sensitivity=False, + file_name='LWE.h5', + ): + """ + Set up and initialize a CADET model template. + + This function prepares a CADET model template by invoking the `createLWE` executable + with specified parameters. It supports the configuration of the model type, number + of particle types, inclusion of sensitivity analysis, and the name of the output + file. Depending on the operating system, it adjusts the executable name accordingly. + After creating the model, it initializes a Cadet instance with the specified or + default CADET binary and the created model file. + + Parameters + ---------- + cadet_root : str or Path + The root directory where the CADET software is located. + use_dll : bool, optional + If True, use the in-memory interface for CADET. Otherwise, use the CLI. + The default is True. + model : str, optional + The model type to set up. The default is 'GENERAL_RATE_MODEL'. + n_partypes : int, optional + The number of particle types. The default is 1. + include_sensitivity : bool, optional + If True, included parameter sensitivities in template. The default is False. + file_name : str, optional + The name of the file to which the CADET model is written. + The default is 'LWE.h5'. + + Returns + ------- + Cadet + An initialized Cadet instance with the model loaded. + + Raises + ------ + Exception + If the creation of the test simulation encounters problems, + detailed in the subprocess's stdout and stderr. + FileNotFoundError + If the CADET executable or DLL file cannot be found at the specified paths. + + Notes + ----- + The function assumes the presence of `createLWE` executable within the `bin` + directory of the `cadet_root` path. The sensitivity analysis, if included, is + configured for column porosity. + + See Also + -------- + Cadet : The class representing a CADET simulation model. + + Examples + -------- + >>> cadet_model = setup_model( + '/path/to/cadet', + use_dll=False, + model='GENERAL_RATE_MODEL', + n_partypes=2, + include_sensitivity=True, + file_name='my_model.h5' + ) + This example sets up a GENERAL_RATE_MODEL with 2 particle types, includes + sensitivity analysis, and writes the model to 'my_model.h5', using the command-line + interface. + """ + executable = 'createLWE' + if platform.system() == 'Windows': + executable += '.exe' + + create_lwe_path = Path(cadet_root) / 'bin' / executable + + args = [ + create_lwe_path.as_posix(), + f'--out {file_name}', + f'--unit {model}', + f'--parTypes {n_partypes}', + ] + + if include_sensitivity: + args.append('--sens COL_POROSITY/-1/-1/-1/-1/-1/-1/0') + + ret = subprocess.run( + args, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + cwd='./' + ) + + if ret.returncode != 0: + if ret.stdout: + print('Output', ret.stdout.decode('utf-8')) + if ret.stderr: + print('Errors', ret.stderr.decode('utf-8')) + raise Exception( + "Failure: Creation of test simulation ran into problems" + ) + + model = Cadet() + # TODO: This should be simplified once #14 is merged. + if use_dll: + cadet_path = cadet_root / 'lib' / 'libcadet.so' + if not cadet_path.exists(): + cadet_path = cadet_root / 'lib' / 'libcadet_d.so' + else: + cadet_path = cadet_root / 'bin' / 'cadet-cli' + + if not cadet_path.exists(): + raise FileNotFoundError("Could not find CADET at: {cadet_root}") + + model.cadet_path = cadet_path + model.filename = file_name + model.load() + + return model + + +def setup_solution_recorder( + model, + split_components=0, + split_ports=0, + single_as_multi_port=0, + ): + """ + Configure the solution recorder for the simulation. + + This function adjusts the model's settings to specify what simulation data should + be recorded, including solutions at various points (inlet, outlet, bulk, etc.), + sensitivities, and their derivatives. It allows for the configuration of how + components and ports are treated in the output data, potentially splitting them for + individual analysis or aggregating them for a more holistic view. + + Parameters + ---------- + model : Cadet + The model instance to be configured for solution recording. + split_components : int, optional + If 1, split component data in the output. The default is 0. + split_ports : int, optional + If 1, split port data in the output. The default is 0. + single_as_multi_port : int, optional + If 1, treat single ports as multiple ports in the output. The default is 0. + + Examples + -------- + >>> model = Cadet() + >>> setup_solution_recorder(model, split_components=1, split_ports=1, single_as_multi_port=1) + This example demonstrates configuring a Cadet model instance for detailed solution + recording, with component and port data split, and single ports treated as multiple + ports. + + """ + + model.root.input['return'].write_solution_times = 1 + model.root.input['return'].write_solution_last = 1 + model.root.input['return'].write_sens_last = 1 + + model.root.input['return'].split_components_data = split_components + model.root.input['return'].split_ports_data = split_ports + model.root.input['return'].single_as_multi_port = single_as_multi_port + + model.root.input['return'].unit_000.write_coordinates = 1 + + model.root.input['return'].unit_000.write_solution_inlet = 1 + model.root.input['return'].unit_000.write_solution_outlet = 1 + model.root.input['return'].unit_000.write_solution_bulk = 1 + model.root.input['return'].unit_000.write_solution_particle = 1 + model.root.input['return'].unit_000.write_solution_solid = 1 + model.root.input['return'].unit_000.write_solution_flux = 1 + model.root.input['return'].unit_000.write_solution_volume = 1 + + model.root.input['return'].unit_000.write_soldot_inlet = 1 + model.root.input['return'].unit_000.write_soldot_outlet = 1 + model.root.input['return'].unit_000.write_soldot_bulk = 1 + model.root.input['return'].unit_000.write_soldot_particle = 1 + model.root.input['return'].unit_000.write_soldot_solid = 1 + model.root.input['return'].unit_000.write_soldot_flux = 1 + model.root.input['return'].unit_000.write_soldot_volume = 1 + + model.root.input['return'].unit_000.write_sens_inlet = 1 + model.root.input['return'].unit_000.write_sens_outlet = 1 + model.root.input['return'].unit_000.write_sens_bulk = 1 + model.root.input['return'].unit_000.write_sens_particle = 1 + model.root.input['return'].unit_000.write_sens_solid = 1 + model.root.input['return'].unit_000.write_sens_flux = 1 + model.root.input['return'].unit_000.write_sens_volume = 1 + + model.root.input['return'].unit_000.write_sensdot_inlet = 1 + model.root.input['return'].unit_000.write_sensdot_outlet = 1 + model.root.input['return'].unit_000.write_sensdot_bulk = 1 + model.root.input['return'].unit_000.write_sensdot_particle = 1 + model.root.input['return'].unit_000.write_sensdot_solid = 1 + model.root.input['return'].unit_000.write_sensdot_flux = 1 + model.root.input['return'].unit_000.write_sensdot_volume = 1 + + model.root.input['return'].unit_000.write_solution_last_unit = 1 + model.root.input['return'].unit_000.write_soldot_last_unit = 1 + + for unit in range(model.root.input.model['nunits']): + model.root.input['return']['unit_{0:03d}'.format(unit)] = model.root.input['return'].unit_000 + + if model.filename is not None: + model.save() + + +def run_simulation_with_options(use_dll, model_options, solution_recorder_options): + """Run a simulation with specified options for the model and solution recorder. + + Initializes and configures a simulation model with given options, sets up the + solution recording parameters, and executes the simulation. This function leverages + `setup_model` to create and initialize the model and `setup_solution_recorder` to + configure how the simulation results should be recorded based on the specified + options. + + Parameters + ---------- + use_dll : bool, optional + If True, use the in-memory interface for CADET. Otherwise, use the CLI. + The default is True. + model_options : dict + A dictionary of options to pass to `setup_model` for initializing the model. + Keys should match the parameter names of `setup_model`, excluding `use_dll`. + solution_recorder_options : dict + A dictionary of options to pass to `setup_solution_recorder` for configuring the + solution recorder. Keys should match the parameter names of + `setup_solution_recorder`. + + Returns + ------- + Cadet + An instance of the Cadet class with the model simulated and loaded. + + Examples + -------- + >>> use_dll = True + >>> model_options = { + ... 'model': 'GENERAL_RATE_MODEL', + ... 'n_partypes': 2, + ... 'include_sensitivity': True, + ... 'file_name': 'model_output.h5' + ... } + >>> solution_recorder_options = { + ... 'split_components': 1, + ... 'split_ports': 1, + ... 'single_as_multi_port': True + ... } + >>> model = run_simulation_with_options(use_dll, model_options, solution_recorder_options) + This example configures and runs a GENERAL_RATE_MODEL with sensitivity analysis + and two particle types, records the solution with specific options, and loads the + simulation results for further analysis. + """ + model = setup_model(cadet_root, use_dll, **model_options) + setup_solution_recorder(model, **solution_recorder_options) + + model.run_load() + + return model + + +# %% Model templates + +cstr_template = { + 'model': 'CSTR', + 'n_partypes': 1, + 'include_sensitivity': False, +} + +lrm_template = { + 'model': 'LUMPED_RATE_MODEL_WITHOUT_PORES', + 'n_partypes': 1, + 'include_sensitivity': False, +} + +lrmp_template = { + 'model': 'LUMPED_RATE_MODEL_WITH_PORES', + 'n_partypes': 1, + 'include_sensitivity': False, +} + +grm_template = { + 'model': 'GENERAL_RATE_MODEL', + 'n_partypes': 1, + 'include_sensitivity': False, +} + +grm_template_sens = { + 'model': 'GENERAL_RATE_MODEL', + 'n_partypes': 1, + 'include_sensitivity': True, +} + +grm_template_partypes = { + 'model': 'GENERAL_RATE_MODEL', + 'n_partypes': 2, + 'include_sensitivity': False, +} + +_2dgrm_template = { + 'model': 'GENERAL_RATE_MODEL_2D', + 'n_partypes': 1, + 'include_sensitivity': False, +} + + +# %% Solution recorder templates + +no_split_options = { + 'split_components': 0, + 'split_ports': 0, + 'single_as_multi_port': 0, +} + +split_ports_options = { + 'split_components': 0, + 'split_ports': 1, + 'single_as_multi_port': 0, +} + +split_all_options = { + 'split_components': 1, + 'split_ports': 1, + 'single_as_multi_port': 1, +} + +# %% Test cases + +class Case(): + def __init__(self, name, model_options, solution_recorder_options, expected_results): + self.name = name + self.model_options = model_options + self.solution_recorder_options = solution_recorder_options + self.expected_results = expected_results + + def __str__(self): + return self.name + + def __repr__(self): + return \ + f"Case('{self.name}', {self.model_options}, " \ + f"{self.solution_recorder_options}, {self.expected_results})" + +# %% CSTR + +cstr = Case( + name='cstr', + model_options=cstr_template, + solution_recorder_options=no_split_options, + expected_results={ + 'last_state_y': (21,), + 'last_state_ydot': (21,), + 'coordinates_unit_000': { + 'axial_coordinates': (1,), + 'particle_coordinates_000': (1,), + }, + 'solution_times': (1501,), + 'solution_unit_000': { + 'solution_inlet': (1501, 4), + 'solution_outlet': (1501, 4), + 'solution_bulk': (1501, 4), + 'solution_solid': (1501, 4), + 'solution_volume': (1501, 1), + 'soldot_inlet': (1501, 4), + 'soldot_outlet': (1501, 4), + 'soldot_bulk': (1501, 4), + 'soldot_solid': (1501, 4), + 'soldot_volume': (1501, 1), + 'last_state_y': (13,), + 'last_state_ydot': (13,), + }, + 'solution_unit_001': { + 'solution_inlet': (1501, 4), + 'solution_outlet': (1501, 4), + 'soldot_inlet': (1501, 4), + 'soldot_outlet': (1501, 4), + 'last_state_y': (4,), + 'last_state_ydot': (4,), + }, + }, +) + +# %% LRM + +lrm = Case( + name='lrm', + model_options=lrm_template, + solution_recorder_options=no_split_options, + expected_results={ + 'last_state_y': (92,), + 'last_state_ydot': (92,), + 'coordinates_unit_000': { + 'axial_coordinates': (10,), + }, + 'solution_times': (1501,), + 'solution_unit_000': { + 'solution_inlet': (1501, 4), + 'solution_outlet': (1501, 4), + 'solution_bulk': (1501, 10, 4), + 'solution_solid': (1501, 10, 4), + 'soldot_inlet': (1501, 4), + 'soldot_outlet': (1501, 4), + 'soldot_bulk': (1501, 10, 4), + 'soldot_solid': (1501, 10, 4), + 'last_state_y': (84,), + 'last_state_ydot': (84,), + }, + 'solution_unit_001': { + 'solution_inlet': (1501, 4), + 'solution_outlet': (1501, 4), + 'soldot_inlet': (1501, 4), + 'soldot_outlet': (1501, 4), + 'last_state_y': (4,), + 'last_state_ydot': (4,), + }, + }, +) + +# %% LRMP + +lrmp = Case( + name='lrmp', + model_options=lrmp_template, + solution_recorder_options=no_split_options, + expected_results={ + 'last_state_y': (172,), + 'last_state_ydot': (172,), + 'coordinates_unit_000': { + 'axial_coordinates': (10,), + 'particle_coordinates_000': (1,), + }, + 'solution_times': (1501,), + 'solution_unit_000': { + 'last_state_y': (164,), + 'last_state_ydot': (164,), + 'soldot_bulk': (1501, 10, 4), + 'soldot_flux': (1501, 1, 10, 4), + 'soldot_inlet': (1501, 4), + 'soldot_outlet': (1501, 4), + 'soldot_particle': (1501, 10, 4), + 'soldot_solid': (1501, 10, 4), + 'solution_bulk': (1501, 10, 4), + 'solution_flux': (1501, 1, 10, 4), + 'solution_inlet': (1501, 4), + 'solution_outlet': (1501, 4), + 'solution_particle': (1501, 10, 4), + 'solution_solid': (1501, 10, 4), + }, + 'solution_unit_001': { + 'last_state_y': (4,), + 'last_state_ydot': (4,), + 'soldot_inlet': (1501, 4), + 'soldot_outlet': (1501, 4), + 'solution_inlet': (1501, 4), + 'solution_outlet': (1501, 4), + }, + }, +) + +# %% GRM + +grm = Case( + name='grm', + model_options=grm_template, + solution_recorder_options=no_split_options, + expected_results={ + 'last_state_y': (412,), + 'last_state_ydot': (412,), + 'coordinates_unit_000': { + 'axial_coordinates': (10,), + 'particle_coordinates_000': (4,), + }, + 'solution_times': (1501,), + 'solution_unit_000': { + 'last_state_y': (404,), + 'last_state_ydot': (404,), + 'soldot_bulk': (1501, 10, 4), + 'soldot_flux': (1501, 1, 10, 4), + 'soldot_inlet': (1501, 4), + 'soldot_outlet': (1501, 4), + 'soldot_particle': (1501, 10, 4, 4), + 'soldot_solid': (1501, 10, 4, 4), + 'solution_bulk': (1501, 10, 4), + 'solution_flux': (1501, 1, 10, 4), + 'solution_inlet': (1501, 4), + 'solution_outlet': (1501, 4), + 'solution_particle': (1501, 10, 4, 4), + 'solution_solid': (1501, 10, 4, 4), + }, + 'solution_unit_001': { + 'last_state_y': (4,), + 'last_state_ydot': (4,), + 'soldot_inlet': (1501, 4), + 'soldot_outlet': (1501, 4), + 'solution_inlet': (1501, 4), + 'solution_outlet': (1501, 4), + }, + }, +) + +# %% GRM Split + +grm_split = Case( + name='grm_split', + model_options=grm_template, + solution_recorder_options=split_all_options, + expected_results={ + 'last_state_y': (412,), + 'last_state_ydot': (412,), + 'coordinates_unit_000': { + 'axial_coordinates': (10,), + 'particle_coordinates_000': (4,), + }, + 'solution_times': (1501,), + 'solution_unit_000': { + 'last_state_y': (404,), + 'last_state_ydot': (404,), + 'soldot_bulk': (1501, 10, 4), + 'soldot_flux': (1501, 1, 10, 4), + 'soldot_inlet_port_000_comp_000': (1501,), + 'soldot_inlet_port_000_comp_001': (1501,), + 'soldot_inlet_port_000_comp_002': (1501,), + 'soldot_inlet_port_000_comp_003': (1501,), + 'soldot_outlet_port_000_comp_000': (1501,), + 'soldot_outlet_port_000_comp_001': (1501,), + 'soldot_outlet_port_000_comp_002': (1501,), + 'soldot_outlet_port_000_comp_003': (1501,), + 'soldot_particle': (1501, 10, 4, 4), + 'soldot_solid': (1501, 10, 4, 4), + 'solution_bulk': (1501, 10, 4), + 'solution_flux': (1501, 1, 10, 4), + 'solution_inlet_port_000_comp_000': (1501,), + 'solution_inlet_port_000_comp_001': (1501,), + 'solution_inlet_port_000_comp_002': (1501,), + 'solution_inlet_port_000_comp_003': (1501,), + 'solution_outlet_port_000_comp_000': (1501,), + 'solution_outlet_port_000_comp_001': (1501,), + 'solution_outlet_port_000_comp_002': (1501,), + 'solution_outlet_port_000_comp_003': (1501,), + 'solution_particle': (1501, 10, 4, 4), + 'solution_solid': (1501, 10, 4, 4), + }, + 'solution_unit_001': { + 'last_state_y': (4,), + 'last_state_ydot': (4,), + 'soldot_inlet_port_000_comp_000': (1501,), + 'soldot_inlet_port_000_comp_001': (1501,), + 'soldot_inlet_port_000_comp_002': (1501,), + 'soldot_inlet_port_000_comp_003': (1501,), + 'soldot_outlet_port_000_comp_000': (1501,), + 'soldot_outlet_port_000_comp_001': (1501,), + 'soldot_outlet_port_000_comp_002': (1501,), + 'soldot_outlet_port_000_comp_003': (1501,), + 'solution_inlet_port_000_comp_000': (1501,), + 'solution_inlet_port_000_comp_001': (1501,), + 'solution_inlet_port_000_comp_002': (1501,), + 'solution_inlet_port_000_comp_003': (1501,), + 'solution_outlet_port_000_comp_000': (1501,), + 'solution_outlet_port_000_comp_001': (1501,), + 'solution_outlet_port_000_comp_002': (1501,), + 'solution_outlet_port_000_comp_003': (1501,), + }, + }, +) + +# %% GRM Sens + +grm_sens = Case( + name='grm_sens', + model_options=grm_template_sens, + solution_recorder_options=no_split_options, + expected_results={ + 'last_state_y': (412,), + 'last_state_ydot': (412,), + 'coordinates_unit_000': { + 'axial_coordinates': (10,), + 'particle_coordinates_000': (4,), + }, + 'solution_times': (1501,), + 'solution_unit_000': { + 'last_state_y': (404,), + 'last_state_ydot': (404,), + 'soldot_bulk': (1501, 10, 4), + 'soldot_flux': (1501, 1, 10, 4), + 'soldot_inlet': (1501, 4), + 'soldot_outlet': (1501, 4), + 'soldot_particle': (1501, 10, 4, 4), + 'soldot_solid': (1501, 10, 4, 4), + 'solution_bulk': (1501, 10, 4), + 'solution_flux': (1501, 1, 10, 4), + 'solution_inlet': (1501, 4), + 'solution_outlet': (1501, 4), + 'solution_particle': (1501, 10, 4, 4), + 'solution_solid': (1501, 10, 4, 4), + }, + 'solution_unit_001': { + 'last_state_y': (4,), + 'last_state_ydot': (4,), + 'soldot_inlet': (1501, 4), + 'soldot_outlet': (1501, 4), + 'solution_inlet': (1501, 4), + 'solution_outlet': (1501, 4), + }, + 'sens_param_000_unit_000': { + 'sensdot_bulk': (1501, 10, 4), + 'sensdot_flux': (1501, 1, 10, 4), + 'sensdot_inlet': (1501, 4), + 'sensdot_outlet': (1501, 4), + 'sensdot_particle': (1501, 10, 4, 4), + 'sensdot_solid': (1501, 10, 4, 4), + 'sens_bulk': (1501, 10, 4), + 'sens_flux': (1501, 1, 10, 4), + 'sens_inlet': (1501, 4), + 'sens_outlet': (1501, 4), + 'sens_particle': (1501, 10, 4, 4), + 'sens_solid': (1501, 10, 4, 4), + }, + 'sens_param_000_unit_001': { + 'sensdot_inlet': (1501, 4), + 'sensdot_outlet': (1501, 4), + 'sens_inlet': (1501, 4), + 'sens_outlet': (1501, 4), + }, + }, +) + +# %% GRM ParTypes + +grm_par_types = Case( + name='grm_par_types', + model_options=grm_template_partypes, + solution_recorder_options=no_split_options, + expected_results={ + 'last_state_y': (772,), + 'last_state_ydot': (772,), + 'coordinates_unit_000': { + 'axial_coordinates': (10,), + 'particle_coordinates_000': (4,), + 'particle_coordinates_001': (4,), + }, + 'solution_times': (1501,), + 'solution_unit_000': { + 'last_state_y': (764,), + 'last_state_ydot': (764,), + 'soldot_bulk': (1501, 10, 4), + 'soldot_flux': (1501, 2, 10, 4), + 'soldot_inlet': (1501, 4), + 'soldot_outlet': (1501, 4), + 'soldot_particle_partype_000': (1501, 10, 4, 4), + 'soldot_particle_partype_001': (1501, 10, 4, 4), + 'soldot_solid_partype_000': (1501, 10, 4, 4), + 'soldot_solid_partype_001': (1501, 10, 4, 4), + 'solution_bulk': (1501, 10, 4), + 'solution_flux': (1501, 2, 10, 4), + 'solution_inlet': (1501, 4), + 'solution_outlet': (1501, 4), + 'solution_particle_partype_000': (1501, 10, 4, 4), + 'solution_particle_partype_001': (1501, 10, 4, 4), + 'solution_solid_partype_000': (1501, 10, 4, 4), + 'solution_solid_partype_001': (1501, 10, 4, 4), + }, + 'solution_unit_001': { + 'last_state_y': (4,), + 'last_state_ydot': (4,), + 'soldot_inlet': (1501, 4), + 'soldot_outlet': (1501, 4), + 'solution_inlet': (1501, 4), + 'solution_outlet': (1501, 4), + }, + }, +) + +# %% 2D GRM +_2dgrm = Case( + name='_2dgrm', + model_options=_2dgrm_template, + solution_recorder_options=no_split_options, + expected_results={ + 'last_state_y': (1228,), + 'last_state_ydot': (1228,), + 'coordinates_unit_000': { + 'axial_coordinates': (10,), + 'particle_coordinates_000': (4,), + 'radial_coordinates': (3,), + }, + 'solution_times': (1501,), + 'solution_unit_000': { + 'last_state_y': (1212,), + 'last_state_ydot': (1212,), + 'soldot_bulk': (1501, 10, 3, 4), + 'soldot_flux': (1501, 1, 10, 3, 4), + 'soldot_inlet': (1501, 3, 4), + 'soldot_outlet': (1501, 3, 4), + 'soldot_particle': (1501, 10, 3, 4, 4), + 'soldot_solid': (1501, 10, 3, 4, 4), + 'solution_bulk': (1501, 10, 3, 4), + 'solution_flux': (1501, 1, 10, 3, 4), + 'solution_inlet': (1501, 3, 4), + 'solution_outlet': (1501, 3, 4), + 'solution_particle': (1501, 10, 3, 4, 4), + 'solution_solid': (1501, 10, 3, 4, 4), + }, + 'solution_unit_001': { + 'last_state_y': (4,), + 'last_state_ydot': (4,), + 'soldot_inlet': (1501, 4), + 'soldot_outlet': (1501, 4), + 'solution_inlet': (1501, 4), + 'solution_outlet': (1501, 4), + }, + }, +) + +# %% 2D GRM Split Ports (single_as_multi_port=False) + +_2dgrm_split_ports = Case( + name='_2dgrm_split_ports', + model_options=_2dgrm_template, + solution_recorder_options=split_ports_options, + expected_results={ + 'last_state_y': (1228,), + 'last_state_ydot': (1228,), + 'coordinates_unit_000': { + 'axial_coordinates': (10,), + 'particle_coordinates_000': (4,), + 'radial_coordinates': (3,), + }, + 'solution_times': (1501,), + 'solution_unit_000': { + 'last_state_y': (1212,), + 'last_state_ydot': (1212,), + 'soldot_bulk': (1501, 10, 3, 4), + 'soldot_flux': (1501, 1, 10, 3, 4), + 'soldot_inlet_port_000': (1501, 4), + 'soldot_inlet_port_001': (1501, 4), + 'soldot_inlet_port_002': (1501, 4), + 'soldot_outlet_port_000': (1501, 4), + 'soldot_outlet_port_001': (1501, 4), + 'soldot_outlet_port_002': (1501, 4), + 'soldot_particle': (1501, 10, 3, 4, 4), + 'soldot_solid': (1501, 10, 3, 4, 4), + 'solution_bulk': (1501, 10, 3, 4), + 'solution_flux': (1501, 1, 10, 3, 4), + 'solution_inlet_port_000': (1501, 4), + 'solution_inlet_port_001': (1501, 4), + 'solution_inlet_port_002': (1501, 4), + 'solution_outlet_port_000': (1501, 4), + 'solution_outlet_port_001': (1501, 4), + 'solution_outlet_port_002': (1501, 4), + 'solution_particle': (1501, 10, 3, 4, 4), + 'solution_solid': (1501, 10, 3, 4, 4), + }, + 'solution_unit_001': { + 'last_state_y': (4,), + 'last_state_ydot': (4,), + 'soldot_inlet': (1501, 4), + 'soldot_outlet': (1501, 4), + 'solution_inlet': (1501, 4), + 'solution_outlet': (1501, 4), + }, + }, +) + +# %% 2D GRM Split All + +_2dgrm_split_all = Case( + name='_2dgrm_split_all', + model_options=_2dgrm_template, + solution_recorder_options=split_all_options, + expected_results={ + 'last_state_y': (1228,), + 'last_state_ydot': (1228,), + 'coordinates_unit_000': { + 'axial_coordinates': (10,), + 'particle_coordinates_000': (4,), + 'radial_coordinates': (3,), + }, + 'solution_times': (1501,), + 'solution_unit_000': { + 'last_state_y': (1212,), + 'last_state_ydot': (1212,), + 'soldot_bulk': (1501, 10, 3, 4), + 'soldot_flux': (1501, 1, 10, 3, 4), + 'soldot_inlet_port_000_comp_000': (1501,), + 'soldot_inlet_port_000_comp_001': (1501,), + 'soldot_inlet_port_000_comp_002': (1501,), + 'soldot_inlet_port_000_comp_003': (1501,), + 'soldot_inlet_port_001_comp_000': (1501,), + 'soldot_inlet_port_001_comp_001': (1501,), + 'soldot_inlet_port_001_comp_002': (1501,), + 'soldot_inlet_port_001_comp_003': (1501,), + 'soldot_inlet_port_002_comp_000': (1501,), + 'soldot_inlet_port_002_comp_001': (1501,), + 'soldot_inlet_port_002_comp_002': (1501,), + 'soldot_inlet_port_002_comp_003': (1501,), + 'soldot_outlet_port_000_comp_000': (1501,), + 'soldot_outlet_port_000_comp_001': (1501,), + 'soldot_outlet_port_000_comp_002': (1501,), + 'soldot_outlet_port_000_comp_003': (1501,), + 'soldot_outlet_port_001_comp_000': (1501,), + 'soldot_outlet_port_001_comp_001': (1501,), + 'soldot_outlet_port_001_comp_002': (1501,), + 'soldot_outlet_port_001_comp_003': (1501,), + 'soldot_outlet_port_002_comp_000': (1501,), + 'soldot_outlet_port_002_comp_001': (1501,), + 'soldot_outlet_port_002_comp_002': (1501,), + 'soldot_outlet_port_002_comp_003': (1501,), + 'soldot_particle': (1501, 10, 3, 4, 4), + 'soldot_solid': (1501, 10, 3, 4, 4), + 'solution_bulk': (1501, 10, 3, 4), + 'solution_flux': (1501, 1, 10, 3, 4), + 'solution_inlet_port_000_comp_000': (1501,), + 'solution_inlet_port_000_comp_001': (1501,), + 'solution_inlet_port_000_comp_002': (1501,), + 'solution_inlet_port_000_comp_003': (1501,), + 'solution_inlet_port_001_comp_000': (1501,), + 'solution_inlet_port_001_comp_001': (1501,), + 'solution_inlet_port_001_comp_002': (1501,), + 'solution_inlet_port_001_comp_003': (1501,), + 'solution_inlet_port_002_comp_000': (1501,), + 'solution_inlet_port_002_comp_001': (1501,), + 'solution_inlet_port_002_comp_002': (1501,), + 'solution_inlet_port_002_comp_003': (1501,), + 'solution_outlet_port_000_comp_000': (1501,), + 'solution_outlet_port_000_comp_001': (1501,), + 'solution_outlet_port_000_comp_002': (1501,), + 'solution_outlet_port_000_comp_003': (1501,), + 'solution_outlet_port_001_comp_000': (1501,), + 'solution_outlet_port_001_comp_001': (1501,), + 'solution_outlet_port_001_comp_002': (1501,), + 'solution_outlet_port_001_comp_003': (1501,), + 'solution_outlet_port_002_comp_000': (1501,), + 'solution_outlet_port_002_comp_001': (1501,), + 'solution_outlet_port_002_comp_002': (1501,), + 'solution_outlet_port_002_comp_003': (1501,), + 'solution_particle': (1501, 10, 3, 4, 4), + 'solution_solid': (1501, 10, 3, 4, 4), + }, + 'solution_unit_001': { + 'last_state_y': (4,), + 'last_state_ydot': (4,), + 'soldot_inlet_port_000_comp_000': (1501,), + 'soldot_inlet_port_000_comp_001': (1501,), + 'soldot_inlet_port_000_comp_002': (1501,), + 'soldot_inlet_port_000_comp_003': (1501,), + 'soldot_outlet_port_000_comp_000': (1501,), + 'soldot_outlet_port_000_comp_001': (1501,), + 'soldot_outlet_port_000_comp_002': (1501,), + 'soldot_outlet_port_000_comp_003': (1501,), + 'solution_inlet_port_000_comp_000': (1501,), + 'solution_inlet_port_000_comp_001': (1501,), + 'solution_inlet_port_000_comp_002': (1501,), + 'solution_inlet_port_000_comp_003': (1501,), + 'solution_outlet_port_000_comp_000': (1501,), + 'solution_outlet_port_000_comp_001': (1501,), + 'solution_outlet_port_000_comp_002': (1501,), + 'solution_outlet_port_000_comp_003': (1501,), + }, + }, +) + +# %% Actual tests + +use_dll = [False, True] +test_cases = [ + cstr, + lrm, + lrmp, + grm, + grm_split, + grm_sens, + grm_par_types, + _2dgrm, + _2dgrm_split_ports, + _2dgrm_split_all +] + +@pytest.mark.parametrize("use_dll", use_dll) +@pytest.mark.parametrize("test_case", test_cases) +def test_simulator_options(use_dll, test_case): + model_options = test_case.model_options + solution_recorder_options = test_case.solution_recorder_options + expected_results = test_case.expected_results + + model = run_simulation_with_options( + use_dll, model_options, solution_recorder_options + ) + + assert model.root.output.last_state_y.shape == expected_results['last_state_y'] + assert model.root.output.last_state_ydot.shape == expected_results['last_state_ydot'] + + for key, value in expected_results['coordinates_unit_000'].items(): + assert model.root.output.coordinates.unit_000[key].shape == value + + assert model.root.output.solution.solution_times.shape == expected_results['solution_times'] + + for key, value in expected_results['solution_unit_000'].items(): + assert model.root.output.solution.unit_000[key].shape == value + + for key, value in expected_results['solution_unit_001'].items(): + assert model.root.output.solution.unit_001[key].shape == value + + if model_options['include_sensitivity']: + for key, value in expected_results['sens_param_000_unit_000'].items(): + assert model.root.output.sensitivity.param_000.unit_000[key].shape == value + + if model_options['include_sensitivity']: + for key, value in expected_results['sens_param_000_unit_001'].items(): + assert model.root.output.sensitivity.param_000.unit_001[key].shape == value + + +if __name__ == "__main__": + pytest.main(["test_dll.py"])