Repository: cadet/CADET-Python (language: python)

Issue:
Add autodetect cadetpath functionality from CADET-Process to CADET-Python
This PR adds the CADET autodetection functionality from CADET-Process into CADET-Python


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/pipeline.yml b/.github/workflows/pipeline.yml
index c19d50a..14033df 100644
--- a/.github/workflows/pipeline.yml
+++ b/.github/workflows/pipeline.yml
@@ -63,4 +63,4 @@ jobs:
       run: |
         pip install .[testing]
         mamba install cadet -c conda-forge
-        pytest tests --rootdir=tests -m "not slow"
+        pytest tests --rootdir=tests -m "not slow and not local"
diff --git a/CADET-Python.pyproj b/CADET-Python.pyproj
deleted file mode 100644
index 04931fc..0000000
--- a/CADET-Python.pyproj
+++ /dev/null
@@ -1,104 +0,0 @@
-﻿<?xml version="1.0" encoding="utf-8"?>
-<Project ToolsVersion="4.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003" DefaultTargets="Build">
-  <PropertyGroup>
-    <Configuration Condition=" '$(Configuration)' == '' ">Debug</Configuration>
-    <SchemaVersion>2.0</SchemaVersion>
-    <ProjectGuid>{12356d41-d8b4-45c3-8413-01e847dcfac2}</ProjectGuid>
-    <ProjectHome />
-    <StartupFile>examples\dll_test.py</StartupFile>
-    <SearchPath />
-    <WorkingDirectory>.</WorkingDirectory>
-    <OutputPath>.</OutputPath>
-    <ProjectTypeGuids>{888888a0-9f3d-457c-b088-3a5042f75d52}</ProjectTypeGuids>
-    <LaunchProvider>Standard Python launcher</LaunchProvider>
-    <InterpreterId>CondaEnv|CondaEnv|cadet_devel</InterpreterId>
-    <IsWindowsApplication>False</IsWindowsApplication>
-  </PropertyGroup>
-  <PropertyGroup Condition="'$(Configuration)' == 'Debug'" />
-  <PropertyGroup Condition="'$(Configuration)' == 'Release'" />
-  <PropertyGroup>
-    <VisualStudioVersion Condition=" '$(VisualStudioVersion)' == '' ">10.0</VisualStudioVersion>
-  </PropertyGroup>
-  <ItemGroup>
-    <Compile Include="cadet\cadet.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="cadet\cadet_dll.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="cadet\LICENSE">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="cadet\__init__.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="examples\cadet_json.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="examples\cadet_library.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="examples\Cadet_Python.py" />
-    <Compile Include="examples\Cadet_Python_Scoop.py" />
-    <Compile Include="examples\common.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="examples\CSTR.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="examples\dll_test.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="examples\fix_h5.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="examples\HICWang.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="examples\library_test.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="examples\MSSMA.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="examples\MSSMA2.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="examples\MSSMA_2comp.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="examples\Pilot_Bypass_noBT_noC.py" />
-    <Compile Include="examples\ron.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="examples\SMB.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="examples\Test.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="examples\test_lwe.py">
-      <SubType>Code</SubType>
-    </Compile>
-    <Compile Include="setup.py">
-      <SubType>Code</SubType>
-    </Compile>
-  </ItemGroup>
-  <ItemGroup>
-    <InterpreterReference Include="CondaEnv|CondaEnv|cadet_devel" />
-    <InterpreterReference Include="Global|ContinuumAnalytics|Anaconda37-64" />
-  </ItemGroup>
-  <ItemGroup>
-    <Folder Include="cadet\" />
-    <Folder Include="examples\" />
-  </ItemGroup>
-  <ItemGroup>
-    <Content Include="README.md">
-      <SubType>Code</SubType>
-    </Content>
-    <Content Include="setup.cfg">
-      <SubType>Code</SubType>
-    </Content>
-  </ItemGroup>
-  <Import Project="$(MSBuildExtensionsPath32)\Microsoft\VisualStudio\v$(VisualStudioVersion)\Python Tools\Microsoft.PythonTools.targets" />
-</Project>
\ No newline at end of file
diff --git a/CADET-Python.sln b/CADET-Python.sln
deleted file mode 100644
index ad69a48..0000000
--- a/CADET-Python.sln
+++ /dev/null
@@ -1,20 +0,0 @@
-﻿
-Microsoft Visual Studio Solution File, Format Version 12.00
-# Visual Studio 15
-VisualStudioVersion = 15.0.26608.5
-MinimumVisualStudioVersion = 10.0.40219.1
-Project("{888888A0-9F3D-457C-B088-3A5042F75D52}") = "CADET-Python", "CADET-Python.pyproj", "{12356D41-D8B4-45C3-8413-01E847DCFAC2}"
-EndProject
-Global
-	GlobalSection(SolutionConfigurationPlatforms) = preSolution
-		Debug|Any CPU = Debug|Any CPU
-		Release|Any CPU = Release|Any CPU
-	EndGlobalSection
-	GlobalSection(ProjectConfigurationPlatforms) = postSolution
-		{12356D41-D8B4-45C3-8413-01E847DCFAC2}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
-		{12356D41-D8B4-45C3-8413-01E847DCFAC2}.Release|Any CPU.ActiveCfg = Release|Any CPU
-	EndGlobalSection
-	GlobalSection(SolutionProperties) = preSolution
-		HideSolutionNode = FALSE
-	EndGlobalSection
-EndGlobal
diff --git a/cadet/LICENSE b/LICENSE
similarity index 96%
rename from cadet/LICENSE
rename to LICENSE
index 6394452..b5cea56 100644
--- a/cadet/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright (c) 2019, William Heymann
+Copyright (c) 2019 - 2024, The CADET-Python Authors
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
diff --git a/cadet/__init__.py b/cadet/__init__.py
index 8978e6f..509a277 100644
--- a/cadet/__init__.py
+++ b/cadet/__init__.py
@@ -1,6 +1,6 @@
-name = 'CADET-python'
+name = 'CADET-Python'
 
 __version__ = "0.14.1"
 
-from .cadet import H5
+from .h5 import H5
 from .cadet import Cadet
diff --git a/cadet/cadet.py b/cadet/cadet.py
index d63d4e4..dcfa6e9 100644
--- a/cadet/cadet.py
+++ b/cadet/cadet.py
@@ -1,284 +1,261 @@
-import copy
-import json
+import os
+from pathlib import Path
 import platform
-import pprint
 import shutil
 import subprocess
+from typing import Optional
 import warnings
-from pathlib import Path
 
-with warnings.catch_warnings():
-    warnings.filterwarnings("ignore", category=FutureWarning)
-    import h5py
-import numpy
-from addict import Dict
+from cadet.h5 import H5
+from cadet.runner import CadetRunnerBase, CadetCLIRunner, ReturnInformation
+from cadet.cadet_dll import CadetDLLRunner
 
-import filelock
-import contextlib
 
-from cadet.cadet_dll import CadetDLL
+def is_dll(path: os.PathLike) -> bool:
+    """
+    Determine if the given path points to a shared library.
 
+    Parameters
+    ----------
+    path : os.PathLike
+        Path to the file.
 
-class H5():
-    pp = pprint.PrettyPrinter(indent=4)
+    Returns
+    -------
+    bool
+        True if the file has a shared library extension (.so, .dll), False otherwise.
+    """
+    suffix = Path(path).suffix
+    return suffix in {'.so', '.dll'}
 
-    def transform(self, x):
-        return x
 
-    def inverse_transform(self, x):
-        return x
+def install_path_to_cadet_paths(
+        install_path: Optional[os.PathLike],
+        ) -> tuple[Optional[Path], Optional[Path], Optional[Path], Optional[Path]]:
+    """
+    Get the correct paths (root_path, cadet_cli_path, cadet_dll_path, cadet_create_lwe_path)
+     from the installation path of CADET. This is extracted into a function to be used
+     in both the meta class and the cadet class itself.
+
+    Parameters
+    ----------
+    install_path : Optional[os.PathLike]
+        Path to the root of the CADET installation or the executable file 'cadet-cli'.
+        If a file path is provided, the root directory will be inferred.
+
+    Returns
+    -------
+    tuple[Optional[Path], Optional[Path], Optional[Path], Optional[Path]]
+        Tuple with CADET installation paths
+        (root_path, cadet_cli_path, cadet_dll_path, cadet_create_lwe_path)
+    """
+    if install_path is None:
+        return None, None, None, None
 
-    def __init__(self, *data):
-        self.root = Dict()
-        self.filename = None
-        for i in data:
-            self.root.update(copy.deepcopy(i))
+    install_path = Path(install_path)
 
-    def load(self, paths=None, update=False, lock=False):
-        if self.filename is not None:
+    if install_path.is_file():
+        cadet_root = install_path.parent.parent
+        warnings.warn(
+            "The specified install_path is not the root of the CADET installation. "
+            "It has been inferred from the file path."
+        )
+    else:
+        cadet_root = install_path
 
-            if lock is True:
-                lock_file = filelock.FileLock(self.filename + '.lock')
-            else:
-                lock_file = contextlib.nullcontext()
+    root_path = cadet_root
 
-            with lock_file:
-                with h5py.File(self.filename, 'r') as h5file:
-                    data = Dict(recursively_load(h5file, '/', self.inverse_transform, paths))
-                    if update:
-                        self.root.update(data)
-                    else:
-                        self.root = data
-        else:
-            print('Filename must be set before load can be used')
-
-    def save(self, lock=False):
-        if self.filename is not None:
-            if lock is True:
-                lock_file = filelock.FileLock(self.filename + '.lock')
-            else:
-                lock_file = contextlib.nullcontext()
-
-            with lock_file:
-                with h5py.File(self.filename, 'w') as h5file:
-                    recursively_save(h5file, '/', self.root, self.transform)
-        else:
-            print("Filename must be set before save can be used")
-
-    def save_json(self, filename):
-        with Path(filename).open("w") as fp:
-            data = convert_from_numpy(self.root, self.transform)
-            json.dump(data, fp, indent=4, sort_keys=True)
-
-    def load_json(self, filename, update=False):
-        with Path(filename).open("r") as fp:
-            data = json.load(fp)
-            data = recursively_load_dict(data, self.inverse_transform)
-            if update:
-                self.root.update(data)
-            else:
-                self.root = data
-
-    def append(self, lock=False):
-        "This can only be used to write new keys to the system, this is faster than having to read the data before writing it"
-        if self.filename is not None:
-            if lock is True:
-                lock_file = filelock.FileLock(self.filename + '.lock')
-            else:
-                lock_file = contextlib.nullcontext()
-
-            with lock_file:
-                with h5py.File(self.filename, 'a') as h5file:
-                    recursively_save(h5file, '/', self.root, self.transform)
-        else:
-            print("Filename must be set before save can be used")
-
-    def __str__(self):
-        temp = []
-        temp.append('Filename = %s' % self.filename)
-        temp.append(self.pp.pformat(self.root))
-        return '\n'.join(temp)
-
-    def update(self, merge):
-        self.root.update(merge.root)
-
-    def __getitem__(self, key):
-        key = key.lower()
-        obj = self.root
-        for i in key.split('/'):
-            if i:
-                obj = obj[i]
-        return obj
-
-    def __setitem__(self, key, value):
-        key = key.lower()
-        obj = self.root
-        parts = key.split('/')
-        for i in parts[:-1]:
-            if i:
-                obj = obj[i]
-        obj[parts[-1]] = value
-
-
-def is_dll(value):
-    suffix = Path(value).suffix
-    return suffix in {'.so', '.dll'}
+    cli_executable = 'cadet-cli'
+    lwe_executable = 'createLWE'
 
+    if platform.system() == 'Windows':
+        cli_executable += '.exe'
+        lwe_executable += '.exe'
 
-class CadetMeta(type):
-    """
-    A meta class for the CADET interface. This allows calls to Cadet.cadet_path = "..." to set
-    the cadet_path for all subsequent Cadet() instances.
-    """
-    _cadet_runner_class = None
-    _is_file_class = None
+    cadet_cli_path = cadet_root / 'bin' / cli_executable
+    if cadet_cli_path.is_file():
+        cadet_cli_path = cadet_cli_path
+    else:
+        raise FileNotFoundError(
+            "CADET could not be found. Please check the path."
+        )
 
-    @property
-    def is_file(cls):
-        return bool(cls._is_file_class)
+    cadet_create_lwe_path = cadet_root / 'bin' / lwe_executable
+    if cadet_create_lwe_path.is_file():
+        cadet_create_lwe_path = cadet_create_lwe_path.as_posix()
 
-    @property
-    def cadet_path(cls):
-        if cls._cadet_runner_class is not None:
-            return cls._cadet_runner_class.cadet_path
+    if platform.system() == 'Windows':
+        dll_path = cadet_root / 'bin' / 'cadet.dll'
+        dll_debug_path = cadet_root / 'bin' / 'cadet_d.dll'
+    else:
+        dll_path = cadet_root / 'lib' / 'libcadet.so'
+        dll_debug_path = cadet_root / 'lib' / 'libcadet_d.so'
 
-    @cadet_path.setter
-    def cadet_path(cls, value):
-        if cls._cadet_runner_class is not None and cls._cadet_runner_class.cadet_path != value:
-            del cls._cadet_runner_class
+    # Look for debug dll if dll is not found.
+    if not dll_path.is_file() and dll_debug_path.is_file():
+        dll_path = dll_debug_path
 
-        if is_dll(value):
-            cls._cadet_runner_class = CadetDLL(value)
-            cls._is_file_class = False
-        else:
-            cls._cadet_runner_class = CadetFile(value)
-            cls._is_file_class = True
+    if dll_path.is_file():
+        cadet_dll_path = dll_path
+    else:
+        cadet_dll_path = None
 
-    @cadet_path.deleter
-    def cadet_path(cls):
-        del cls._cadet_runner_class
+    return root_path, cadet_cli_path, cadet_dll_path, cadet_create_lwe_path
 
 
-class Cadet(H5, metaclass=CadetMeta):
-    # cadet_path must be set in order for simulations to run
-    def __init__(self, *data):
-        super().__init__(*data)
-        self._cadet_runner = None
+class CadetMeta(type):
+    """
+    Meta class for the CADET interface.
 
-        self._is_file = None  # Is CLI or DLL
-        # self.cadet_path  # from Bill, declared in meta class -> path to CLI-file or DLL-file
+    This meta class allows setting the `cadet_path` attribute for all instances of the
+    `Cadet` class.
+    """
 
-        self.install_path = None  # from Jo -> root of the CADET installation.
+    @property
+    def cadet_path(cls) -> Optional[Path]:
+        """
+        Get the current CADET path.
 
-        self.cadet_cli_path = None
-        self.cadet_dll_path = None
-        self.cadet_create_lwe_path = None
+        Returns
+        -------
+        Optional[Path]
+            The current CADET path if set, otherwise None.
+        """
+        if cls.use_dll and cls.cadet_dll_path is not None:
+            return cls.cadet_dll_path
+        elif cls.cadet_cli_path is not None:
+            return cls.cadet_cli_path
+        else:
+            return None
 
-        self.return_information = None
+    @cadet_path.setter
+    def cadet_path(cls, cadet_path: os.PathLike) -> None:
+        """
+        Set the CADET path and initialize the appropriate runner.
 
-    @property
-    def is_file(self):
-        if self._is_file is not None:
-            return bool(self._is_file)
-        if self._is_file_class is not None:
-            return bool(self._is_file_class)
+        Parameters
+        ----------
+        cadet_path : os.PathLike
+            Path to the CADET executable or library.
 
-    @property
-    def cadet_runner(self):
-        if self._cadet_runner is not None:
-            return self._cadet_runner
-        if hasattr(self, "_cadet_runner_class") and self._cadet_runner_class is not None:
-            return self._cadet_runner_class
-        self.autodetect_cadet()
-        return self._cadet_runner
+        Notes
+        -----
+        If the path is a DLL, a `CadetDLLRunner` runner is used.
+        Otherwise, a `CadetFileRunner` runner is used.
+        """
+        cadet_path = Path(cadet_path)
 
-    @property
-    def cadet_path(self):
-        runner = self.cadet_runner
-        if runner is not None:
-            return runner.cadet_path
+        warnings.warn(
+            "Support for setting cadet.cadet_path will be removed in a future version. "
+            "TODO: What alternative should be used?",
+            DeprecationWarning
+        )
 
-    @cadet_path.setter
-    def cadet_path(self, value):
-        if self._cadet_runner is not None and self._cadet_runner.cadet_path != value:
-            del self._cadet_runner
+        cls.use_dll = cadet_path.suffix in [".dll", ".so"]
 
-        if is_dll(value):
-            self._cadet_runner = CadetDLL(value)
-            self._is_file = False
-        else:
-            self._cadet_runner = CadetFile(value)
-            self._is_file = True
+        _install_path, cadet_cli_path, cadet_dll_path, cadet_create_lwe_path = install_path_to_cadet_paths(cadet_path)
 
-    @cadet_path.deleter
-    def cadet_path(self):
-        del self._cadet_runner
+        cls._install_path = _install_path
+        cls.cadet_create_lwe_path = cadet_create_lwe_path
+        cls.cadet_cli_path = cadet_cli_path
+        cls.cadet_dll_path = cadet_dll_path
 
-    def autodetect_cadet(self):
+
+class Cadet(H5, metaclass=CadetMeta):
+    """
+    CADET interface class.
+
+    This class manages the CADET runner, whether it's based on a CLI executable or a DLL,
+    and provides methods for running simulations and loading results.
+
+    Attributes
+    ----------
+    install_path : Optional[Path]
+        The root directory of the CADET installation.
+    cadet_cli_path : Optional[Path]
+        Path to the 'cadet-cli' executable.
+    cadet_dll_path : Optional[Path]
+        Path to the 'cadet.dll' or equivalent shared library.
+    cadet_create_lwe_path : Optional[Path]
+        Path to the 'createLWE' executable.
+    return_information : Optional[dict]
+        Stores the information returned after a simulation run.
+    """
+
+    def __init__(self, install_path: Optional[Path] = None, use_dll: bool = False, *data):
         """
-        Autodetect installation CADET based on operating system and API usage.
+        Initialize a new instance of the Cadet class.
+        Priority order of install_paths is:
+        1. install_path set in __init__ args
+        2. install_path set in CadetMeta
+        3. auto-detected install_path
 
-        Returns
-        -------
-        cadet_root : Path
-            Installation path of the CADET program.
+        Parameters
+        ----------
+        *data : tuple
+            Additional data to be passed to the H5 base class initialization.
         """
-        executable = 'cadet-cli'
-        if platform.system() == 'Windows':
-            executable += '.exe'
+        super().__init__(*data)
 
-        # Searching for the executable in system path
-        path = shutil.which(executable)
+        self.cadet_create_lwe_path: Optional[Path] = None
+        self.return_information: Optional[dict] = None
 
-        if path is None:
-            raise FileNotFoundError(
-                "Could not autodetect CADET installation. Please provide path."
+        # Regardless of settings in the Meta Class, if we get an install_path, we respect the install_path
+        if install_path is not None:
+            self.use_dll = use_dll
+            self.install_path = install_path  # This will set _cadet_dll_runner and _cadet_cli_runner
+            return
+
+        # If _cadet_cli_runner_class has been set in the Meta Class, use them, else instantiate Nones
+        if hasattr(self, "cadet_cli_path") and self.cadet_cli_path is not None:
+            self._cadet_cli_runner: Optional[CadetCLIRunner] = CadetCLIRunner(
+                self.cadet_cli_path
             )
+        else:
+            self._cadet_cli_runner: Optional[CadetCLIRunner] = None
+            self.use_dll = use_dll
 
-        cli_path = Path(path)
+        if hasattr(self, "cadet_dll_path") and self.cadet_dll_path is not None:
+            try:
+                self._cadet_dll_runner: Optional[CadetDLLRunner] = CadetDLLRunner(
+                    self.cadet_dll_path
+                )
+            except ValueError:
+                self.cadet_dll_path = None
+                self._cadet_dll_runner: Optional[CadetCLIRunner] = None
+                self.use_dll = False
+        else:
+            self._cadet_dll_runner: Optional[CadetCLIRunner] = None
+            self.use_dll = use_dll
 
-        cadet_root = None
-        if cli_path is not None:
-            cadet_root = cli_path.parent.parent
-            self.install_path = cadet_root
+        # If any runner has been set in the Meta Class, don't auto-detect Cadet, just return
+        if self._cadet_cli_runner is not None or self._cadet_dll_runner is not None:
+            return
 
-        return cadet_root
+        # If neither Meta Class nor install_path are given: auto-detect Cadet
+        self.install_path = self.autodetect_cadet()
 
     @property
-    def install_path(self):
-        """str: Path to the installation of CADET.
-
-        This can either be the root directory of the installation or the path to the
-        executable file 'cadet-cli'. If a file path is provided, the root directory will
-        be inferred.
-
-        Raises
-        ------
-        FileNotFoundError
-            If CADET cannot be found at the specified path.
-
-        Warnings
-        --------
-        If the specified install_path is not the root of the CADET installation, it will
-        be inferred from the file path.
+    def install_path(self) -> Optional[Path]:
+        """
+        Path to the installation of CADET.
 
-        See Also
-        --------
-        check_cadet
+        Returns
+        -------
+        Optional[Path]
+            The root directory of the CADET installation or the path to 'cadet-cli'.
         """
         return self._install_path
 
     @install_path.setter
-    def install_path(self, install_path):
+    def install_path(self, install_path: Optional[os.PathLike]) -> None:
         """
         Set the installation path of CADET.
 
         Parameters
         ----------
-        install_path : str or Path
-            Path to the root of the CADET installation.
-            It should either be the root directory of the installation or the path
-            to the executable file 'cadet-cli'.
+        install_path : Optional[os.PathLike]
+            Path to the root of the CADET installation or the executable file 'cadet-cli'.
             If a file path is provided, the root directory will be inferred.
         """
         if install_path is None:
@@ -288,205 +265,263 @@ def install_path(self, install_path):
             self.cadet_create_lwe_path = None
             return
 
-        install_path = Path(install_path)
+        root_path, cadet_cli_path, cadet_dll_path, create_lwe_path = install_path_to_cadet_paths(install_path)
 
-        if install_path.is_file():
-            cadet_root = install_path.parent.parent
-            warnings.warn(
-                "The specified install_path is not the root of the CADET installation. "
-                "It has been inferred from the file path."
-            )
-        else:
-            cadet_root = install_path
+        self._install_path = root_path
+        self.cadet_create_lwe_path = create_lwe_path
+
+        if cadet_cli_path is not None:
+            self._cadet_cli_runner = CadetCLIRunner(cadet_cli_path)
+            self.cadet_cli_path = cadet_cli_path
+
+        self.cadet_dll_path = cadet_dll_path
+        if cadet_dll_path is not None:
+            try:
+                self._cadet_dll_runner = CadetDLLRunner(cadet_dll_path)
+            except ValueError:
+                pass
+
+    @property
+    def cadet_path(self) -> Optional[Path]:
+        """
+        Get the path to the current CADET executable or library.
 
-        self._install_path = cadet_root
+        Returns
+        -------
+        Path
+            The path to the current CADET executable or library if set, otherwise None.
+        """
+        runner = self.cadet_runner
+        if runner is not None:
+            return runner.cadet_path
+        return None
 
-        cli_executable = 'cadet-cli'
-        lwe_executable = 'createLWE'
+    @cadet_path.setter
+    def cadet_path(self, cadet_path: os.PathLike) -> None:
+        """
+        Set the CADET path and initialize the appropriate runner.
 
+        Parameters
+        ----------
+        cadet_path : os.PathLike
+            Path to the CADET executable or library.
+
+        Notes
+        -----
+        If the path is a DLL, a `CadetDLLRunner` runner is used.
+        Otherwise, a `CadetFileRunner` runner is used.
+        """
+        cadet_path = Path(cadet_path)
+
+        warnings.warn(
+            "Deprecation warning: Support for setting cadet.cadet_path will be removed in a future version.",
+            FutureWarning
+        )
+        self.install_path = cadet_path
+
+    @staticmethod
+    def autodetect_cadet() -> Optional[Path]:
+        """
+        Autodetect the CADET installation path.
+
+        Returns
+        -------
+        Optional[Path]
+            The root directory of the CADET installation.
+
+        Raises
+        ------
+        FileNotFoundError
+            If CADET cannot be found in the system path.
+        """
+        executable = 'cadet-cli'
         if platform.system() == 'Windows':
-            cli_executable += '.exe'
-            lwe_executable += '.exe'
+            executable += '.exe'
 
-        cadet_cli_path = cadet_root / 'bin' / cli_executable
-        if cadet_cli_path.is_file():
-            self.cadet_cli_path = cadet_cli_path
-            self.cadet_path = cadet_cli_path
-        else:
+        path = shutil.which(executable)
+
+        if path is None:
             raise FileNotFoundError(
-                "CADET could not be found. Please check the path"
+                "Could not autodetect CADET installation. Please provide path."
             )
 
-        cadet_create_lwe_path = cadet_root / 'bin' / lwe_executable
-        if cadet_create_lwe_path.is_file():
-            self.cadet_create_lwe_path = cadet_create_lwe_path.as_posix()
+        cli_path = Path(path)
+        cadet_root = cli_path.parent.parent if cli_path else None
 
-        if platform.system() == 'Windows':
-            dll_path = cadet_root / 'bin' / 'cadet.dll'
-            dll_debug_path = cadet_root / 'bin' / 'cadet_d.dll'
-        else:
-            dll_path = cadet_root / 'lib' / 'lib_cadet.so'
-            dll_debug_path = cadet_root / 'lib' / 'lib_cadet_d.so'
+        return cadet_root
 
-        # Look for debug dll if dll is not found.
-        if not dll_path.is_file() and dll_debug_path.is_file():
-            dll_path = dll_debug_path
+    @property
+    def cadet_runner(self) -> CadetRunnerBase:
+        """
+        Get the current CADET runner instance.
 
-        # Look for debug dll if dll is not found.
-        if dll_path.is_file():
-            self.cadet_dll_path = dll_path.as_posix()
+        Returns
+        -------
+        Optional[CadetRunnerBase]
+            The current runner instance, either a DLL or file-based runner.
+        """
 
-    def transform(self, x):
-        return str.upper(x)
+        if self.use_dll and self.found_dll:
+            return self._cadet_dll_runner
 
-    def inverse_transform(self, x):
-        return str.lower(x)
+        if self.use_dll and not self.found_dll:
+            raise ValueError("Set Cadet to use_dll but no dll interface found.")
 
-    def load_results(self):
-        runner = self.cadet_runner
-        if runner is not None:
-            runner.load_results(self)
+        return self._cadet_cli_runner
 
-    def run(self, timeout=None, check=None):
-        data = self.cadet_runner.run(simulation=self.root.input, filename=self.filename, timeout=timeout, check=check)
-        # self.return_information = data
-        return data
+    def create_lwe(self, file_path=None):
+        """Create basic LWE example and loads the configuration into self.
 
-    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)
-        # self.return_information = data
-        self.load_results()
-        if clear:
-            self.clear()
-        return data
+        Parameters
+        ----------
+        file_path : Path, optional
+            Path to store HDF5 file. If None, temporary file will be created and
+            deleted after simulation.
 
-    def clear(self):
-        runner = self.cadet_runner
-        if runner is not None:
-            runner.clear()
+        """
+        file_path_input = file_path
+        if file_path is None:
+            file_name = "LWE.h5"
+            cwd = os.getcwd()
+            file_path = Path(cwd) / file_name
+        else:
+            file_path = Path(file_path).absolute()
+            file_name = file_path.name
+            cwd = file_path.parent.as_posix()
+
+        ret = subprocess.run(
+            [self.cadet_create_lwe_path, '-o', file_name],
+            stdout=subprocess.PIPE,
+            stderr=subprocess.PIPE,
+            cwd=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 RuntimeError(
+                "Failure: Creation of test simulation ran into problems"
+            )
 
+        self.filename = file_path
 
-class CadetFile:
+        self.load()
 
-    def __init__(self, cadet_path):
-        self.cadet_path = cadet_path
+        if file_path_input is None:
+            os.remove(file_path)
 
-    def run(self, filename=None, simulation=None, timeout=None, check=None):
-        if filename is not None:
-            data = subprocess.run([self.cadet_path, filename], timeout=timeout, check=check, capture_output=True)
-            return data
-        else:
-            print("Filename must be set before run can be used")
+        return self
+
+    @property
+    def found_dll(self):
+        """
+        Check if a DLL interface was found.
 
-    def clear(self):
-        pass
+        Returns
+        -------
+        bool
+            True if a cadet DLL interface was found.
+            False otherwise.
+        """
+        return self.cadet_dll_path is not None
 
-    def load_results(self, sim):
-        sim.load(paths=["/meta", "/output"], update=True)
+    def transform(self, x: str) -> str:
+        """
+        Transform the input string to uppercase.
 
+        Parameters
+        ----------
+        x : str
+            Input string.
 
-def convert_from_numpy(data, func):
-    ans = Dict()
-    for key_original, item in data.items():
-        key = func(key_original)
-        if isinstance(item, numpy.ndarray):
-            item = item.tolist()
+        Returns
+        -------
+        str
+            Transformed string in uppercase.
+        """
+        return str.upper(x)
 
-        if isinstance(item, numpy.generic):
-            item = item.item()
+    def inverse_transform(self, x: str) -> str:
+        """
+        Transform the input string to lowercase.
 
-        if isinstance(item, bytes):
-            item = item.decode('ascii')
+        Parameters
+        ----------
+        x : str
+            Input string.
 
-        if isinstance(item, Dict):
-            ans[key_original] = convert_from_numpy(item, func)
-        else:
-            ans[key] = item
-    return ans
+        Returns
+        -------
+        str
+            Transformed string in lowercase.
+        """
+        return str.lower(x)
 
+    def run(
+            self,
+            timeout: Optional[int] = None,
+    ) -> ReturnInformation:
+        """
+        Run the CADET simulation.
 
-def recursively_load_dict(data, func):
-    ans = Dict()
-    for key_original, item in data.items():
-        key = func(key_original)
-        if isinstance(item, dict):
-            ans[key] = recursively_load_dict(item, func)
-        else:
-            ans[key] = item
-    return ans
+        Parameters
+        ----------
+        timeout : Optional[int]
+            Maximum time allowed for the simulation to run, in seconds.
 
+        Returns
+        -------
+        ReturnInformation
+            Information about the simulation run.
+        """
+        return_information = self.cadet_runner.run(
+            self,
+            timeout=timeout,
+        )
+
+        return return_information
+
+    def run_load(
+            self,
+            timeout: Optional[int] = None,
+            clear: bool = True
+    ) -> ReturnInformation:
+        """
+        Run the CADET simulation and load the results.
 
-def set_path(obj, path, value):
-    "paths need to be broken up so that subobjects are correctly made"
-    path = path.split('/')
-    path = [i for i in path if i]
+        Parameters
+        ----------
+        timeout : Optional[int]
+            Maximum time allowed for the simulation to run, in seconds.
+        clear : bool
+            If True, clear previous results after loading new ones.
 
-    temp = obj
-    for part in path[:-1]:
-        temp = temp[part]
+        Returns
+        -------
+        ReturnInformation
+            Information about the simulation run.
+        """
+        return_information = self.run(timeout)
+        self.load_results()
 
-    temp[path[-1]] = value
+        if clear:
+            self.clear()
+        return return_information
 
+    def load_results(self) -> None:
+        """Load the results of the last simulation run into the current instance."""
+        runner = self.cadet_runner
+        if runner is not None:
+            runner.load_results(self)
 
-def recursively_load(h5file, path, func, paths):
-    ans = Dict()
-    if paths is not None:
-        for path in paths:
-            item = h5file.get(path, None)
-            if item is not None:
-                if isinstance(item, h5py._hl.dataset.Dataset):
-                    set_path(ans, path, item[()])
-                elif isinstance(item, h5py._hl.group.Group):
-                    set_path(ans, path, recursively_load(h5file, path + '/', func, None))
-    else:
-        for key_original in h5file[path].keys():
-            key = func(key_original)
-            local_path = path + key
-            item = h5file[path][key_original]
-            if isinstance(item, h5py._hl.dataset.Dataset):
-                ans[key] = item[()]
-            elif isinstance(item, h5py._hl.group.Group):
-                ans[key] = recursively_load(h5file, local_path + '/', func, None)
-    return ans
-
-
-def recursively_save(h5file, path, dic, func):
-    if not isinstance(path, str):
-        raise ValueError("path must be a string")
-    if not isinstance(h5file, h5py._hl.files.File):
-        raise ValueError("must be an open h5py file")
-
-    # argument type checking
-    if not isinstance(dic, dict):
-        raise ValueError("must provide a dictionary")
-
-    # save items to the hdf5 file
-    for key, item in dic.items():
-        key = str(key)
-        value = None
-
-        if not isinstance(key, str):
-            raise ValueError("dict keys must be strings to save to hdf5")
-
-        if isinstance(item, dict):
-            recursively_save(h5file, path + key + '/', item, func)
-
-        # handle int, float, string and ndarray of int32, int64, float64
-        elif isinstance(item, str):
-            value = numpy.array(item.encode('ascii'))
-        elif isinstance(item, list) and all(isinstance(i, str) for i in item):
-            value = numpy.array([i.encode('ascii') for i in item])
-        else:
-            try:
-                value = numpy.array(item)
-            except TypeError:
-                raise ValueError('Cannot save %s/%s key with %s type.' % (path, func(key), type(item)))
+    def clear(self) -> None:
+        """Clear the loaded results from the current instance."""
+        runner = self.cadet_runner
+        if runner is not None:
+            runner.clear()
 
-        if value is not None:
-            try:
-                h5file[path + func(key)] = value
-            except OSError as e:
-                if str(e) == 'Unable to create link (name already exists)':
-                    raise KeyError(f'Name conflict with upper and lower case entries for key "{path}{key}".')
-                else:
-                    raise
+    def __del__(self):
+        self.clear()
+        del self._cadet_dll_runner
+        del self._cadet_cli_runner
diff --git a/cadet/cadet_dll.py b/cadet/cadet_dll.py
index 88920c5..487853b 100644
--- a/cadet/cadet_dll.py
+++ b/cadet/cadet_dll.py
@@ -1,510 +1,2134 @@
-
-from cmath import sin
 import ctypes
-import numpy
+import io
+import os
+from pathlib import Path
+from typing import Any, Optional
+
 import addict
-import cadet.cadet_dll_parameterprovider as cadet_dll_parameterprovider
+import numpy
 
-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') ))
+from cadet.runner import CadetRunnerBase
+import cadet.cadet_dll_parameterprovider as cadet_dll_parameterprovider
 
-c_cadet_result = ctypes.c_int
 
+# Common types for ctypes function signatures
+CadetDriver = ctypes.c_void_p
 array_double = ctypes.POINTER(ctypes.POINTER(ctypes.c_double))
-
-
+point_bool = ctypes.POINTER(ctypes.c_bool)
 point_int = ctypes.POINTER(ctypes.c_int)
+point_double = ctypes.POINTER(ctypes.c_double)
 
-def null(*args):
-    pass
-
-if 0:
-    log_print = print
-else:
-    log_print = null
-
-
-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')
-
+# Values of cdtResult
+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:
+    """
+    Definition of CADET-C-API v1.0 function signatures and type mappings.
+
+    Attributes
+    ----------
+    signatures : dict
+        Signatures of exported API functions.
+    lookup_prototype : dict
+        Mapping of common ctypes parameters.
+    lookup_output_argument_type : dict
+        Mapping of ctypes output parameters for the API.
+    """
+
+    # API function signatures
+    # Note, order is important, it has to match the cdtAPIv010000 struct of the C-API
+    signatures = {}
+
+    signatures['getFileFormat'] = ('return', 'fileFormat')
+
+    signatures['createDriver'] = ('drv',)
+    signatures['deleteDriver'] = (None, 'drv')
+    signatures['runSimulation'] = ('return', 'drv', 'parameterProvider')
+
+    signatures['getNumUnitOp'] = ('return', 'drv', 'nUnits')
+    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', 'keepAxialSingletonDimension')
+    signatures['getSolutionParticle'] = ('return', 'drv', 'unitOpId', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nComp', 'keepAxialSingletonDimension', 'keepParticleSingletonDimension')
+    signatures['getSolutionSolid'] = ('return', 'drv', 'unitOpId', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nBound', 'keepAxialSingletonDimension', 'keepParticleSingletonDimension')
+    signatures['getSolutionFlux'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParTypes', 'nComp', 'keepAxialSingletonDimension')
+    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', 'keepAxialSingletonDimension')
+    signatures['getSolutionDerivativeParticle'] = ('return', 'drv', 'unitOpId', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nComp', 'keepAxialSingletonDimension', 'keepParticleSingletonDimension')
+    signatures['getSolutionDerivativeSolid'] = ('return', 'drv', 'unitOpId', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nBound', 'keepAxialSingletonDimension', 'keepParticleSingletonDimension')
+    signatures['getSolutionDerivativeFlux'] = ('return', 'drv', 'unitOpId', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParTypes', 'nComp', 'keepAxialSingletonDimension')
+    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', 'keepAxialSingletonDimension')
+    signatures['getSensitivityParticle'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nComp', 'keepAxialSingletonDimension', 'keepParticleSingletonDimension')
+    signatures['getSensitivitySolid'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nBound', 'keepAxialSingletonDimension', 'keepParticleSingletonDimension')
+    signatures['getSensitivityFlux'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParTypes', 'nComp', 'keepAxialSingletonDimension')
+    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', 'keepAxialSingletonDimension')
+    signatures['getSensitivityDerivativeParticle'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nComp', 'keepAxialSingletonDimension', 'keepParticleSingletonDimension')
+    signatures['getSensitivityDerivativeSolid'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'parType', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParShells', 'nBound', 'keepAxialSingletonDimension', 'keepParticleSingletonDimension')
+    signatures['getSensitivityDerivativeFlux'] = ('return', 'drv', 'unitOpId', 'sensIdx', 'time', 'data', 'nTime', 'nAxialCells', 'nRadialCells', 'nParTypes', 'nComp', 'keepAxialSingletonDimension')
+    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')
+
+    signatures['getTimeSim'] = ('return', 'drv', 'timeSim')
+
+    # Mappings for common ctypes parameters
     lookup_prototype = {
+        None: None,
         'return': c_cadet_result,
-        'drv': ctypes.c_void_p,
-        'unitOpId': ctypes.c_int,
-        'idx': ctypes.c_int,
-        'parType': ctypes.c_int,
-        'time': array_double,
-        'data': array_double,
+        'fileFormat': point_int,
+        'drv': CadetDriver,
+        'parameterProvider': ctypes.POINTER(cadet_dll_parameterprovider.PARAMETERPROVIDER),
+        'nUnits': point_int,
         '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,
-        None: None,
-        'parameterProvider': ctypes.POINTER(cadet_dll_parameterprovider.PARAMETERPROVIDER)
+        'nStates': point_int,
+        'unitOpId': ctypes.c_int,
+        'sensIdx': ctypes.c_int,
+        'parType': ctypes.c_int,
+        'time': array_double,
+        'data': array_double,
+        'state': array_double,
+        'coords': array_double,
+        'keepAxialSingletonDimension': point_bool,
+        'keepParticleSingletonDimension': point_bool,
+        'timeSim': point_double,
     }
 
-    lookup_call = {
-        'time': ctypes.POINTER(ctypes.c_double),
+    lookup_output_argument_type = {
+        'fileFormat': ctypes.c_int,
+        'nUnits': ctypes.c_int,
         'nTime': ctypes.c_int,
-        'data': ctypes.POINTER(ctypes.c_double),
+        'nCoords': ctypes.c_int,
         'nPort': ctypes.c_int,
         '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,
+        'nStates': ctypes.c_int,
+        'time': ctypes.POINTER(ctypes.c_double),
+        'data': ctypes.POINTER(ctypes.c_double),
+        'state': ctypes.POINTER(ctypes.c_double),
+        'coords': ctypes.POINTER(ctypes.c_double),
+        'keepAxialSingletonDimension': ctypes.c_bool,
+        'keepParticleSingletonDimension': ctypes.c_bool,
+        'timeSim': ctypes.c_double,
     }
 
 
-def setup_api():
-    _fields_ = []
+def _setup_api() -> list[tuple[str, ctypes.CFUNCTYPE]]:
+    """
+    Set up the API function prototypes for CADETAPIV010000.
 
-    for key, value in CADETAPIV010000_DATA._data_.items():
+    Returns
+    -------
+    list of tuple
+        List of function names and corresponding ctypes function prototypes.
+    """
+    _fields_ = []
+    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)) )
+        _fields_.append((key, ctypes.CFUNCTYPE(*args)))
 
     return _fields_
 
 
 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:
+    """
+    Handles reading results from a CADET simulation.
+
+    Parameters
+    ----------
+    api : CADETAPIV010000
+        The CADET API instance.
+    driver : CadetDriver
+        The driver handle used for the simulation.
 
-    def __init__(self, api, driver):
+    """
+
+    def __init__(self, api: CADETAPIV010000, driver: CadetDriver) -> None:
         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: Optional[int] = None,
+            sensIdx: Optional[int] = None,
+            parType: Optional[int] = None
+            ) -> Optional[dict[str, Any]]:
+        """
+        Load data from the CADET API.
+
+        Parameters
+        ----------
+        get_solution_str : str
+            The name of the API function to call.
+        unitOpId : Optional[int], default=None
+            The unit operation ID.
+        sensIdx : Optional[int], default=None
+            The sensitivity index.
+        parType : Optional[int], default=None
+            The particle type index.
+
+        Returns
+        -------
+        Optional[dict[str, Any]]
+            The output arguments mapped to their values, or None if data is not available.
+        """
+        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
-
-        result = get_solution(*tuple(wrappers[var_key](var_value) for var_key, var_value in vars.items()))
-
+                _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: dict[str, Any],
+            data_key: str,
+            len_key: str,
+            own_data: bool = True,
+            ) -> Optional[numpy.ndarray]:
+        """
+        Process array data from the API.
+
+        Parameters
+        ----------
+        call_outputs : dict
+            The output arguments returned from the API call.
+        data_key : str
+            The key for the data array in the output arguments.
+        len_key : str
+            The key for the length of the data array.
+        own_data : bool, default=True
+            Whether to create a copy of the data.
+
+        Returns
+        -------
+        Optional[numpy.ndarray]
+            The processed data array or None if the length is zero.
+        """
+        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
+
+    def _process_data(
+            self,
+            call_outputs: dict[str, Any],
+            own_data: bool = True,
+            ) -> Optional[tuple[numpy.ndarray, numpy.ndarray, list[str]]]:
+        """
+        Process multi-dimensional data from the API.
+
+        Parameters
+        ----------
+        call_outputs : dict
+            The output arguments returned from the API call.
+        own_data : bool, default=True
+            Whether to create a copy of the data.
+
+        Returns
+        -------
+        Optional[tuple[numpy.ndarray, numpy.ndarray, list[str]]]
+            A tuple of time, data, and dimensions, or None if no data is available.
+        """
         shape = []
         dims = []
-        dimensions = ['nTime', 'nPort', 'nParShells', 'nAxialCells', 'nRadialCells', 'nComp', 'nBound']
+        drop_indices = []
+
+        # 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 'data' in call_outputs:
+            if 'nAxialCells' in dims:
+                nAxialCells = call_outputs['nAxialCells'].value
+                keep_axial_singleton_dim = call_outputs['keepAxialSingletonDimension'].value
+                if nAxialCells == 1 and not keep_axial_singleton_dim:
+                    drop_indices.append(dims.index('nAxialCells'))
+
+            if 'nParShells' in dims:
+                nParShells = call_outputs['nParShells'].value
+                keep_particle_singleton_dim = call_outputs['keepParticleSingletonDimension'].value
+                if nParShells == 1 and not keep_particle_singleton_dim:
+                    drop_indices.append(dims.index('nParShells'))
+
+            shape = numpy.delete(shape, drop_indices)
+
+            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: Any,
+            own_data: bool = True,
+            **kwargs: Any
+            ) -> Optional[tuple[numpy.ndarray, numpy.ndarray, list[str]]]:
+        """
+        Load and process data from the API in a single step.
+
+        Parameters
+        ----------
+        own_data : bool, default=True
+            Whether to create a copy of the data.
+
+        Returns
+        -------
+        Optional[tuple[numpy.ndarray, numpy.ndarray, list[str]]]
+            A tuple of time, data, and dimensions, or None if no data is available.
+        """
+        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: Any,
+            own_data: bool = True,
+            **kwargs: Any
+            ) -> Optional[numpy.ndarray]:
+        """
+        Load and process a 1D array from the API.
+
+        Parameters
+        ----------
+        data_key : str
+            The key for the data array in the output arguments.
+        len_key : str
+            The key for the length of the data array.
+        own_data : bool, default=True
+            Whether to create a copy of the data.
+
+        Returns
+        -------
+        Optional[numpy.ndarray]
+            The processed data array or None if no data is available.
+        """
+        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 file_format(
+            self,
+            own_data: bool = True,
+            ) -> float:
+        """
+        Load the file format used for the simulation config.
+
+        Parameters
+        ----------
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        int
+            The file format.
+        """
+        call_outputs = self._load_data('getFileFormat')
+        return int(call_outputs['fileFormat'].value)
+
+    def nunits(self) -> int:
+        """
+        Get the number of unit operations in the system.
+
+        Returns
+        -------
+        int
+            The number of unit operations in the system.
+        """
+        call_outputs = self._load_data('getNumUnitOp')
+        return int(call_outputs['nUnits'].value)
+
+    def npartypes(self, unitOpId: int) -> int:
+        """
+        Get the number of particle types for a given unit operation.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+
+        Returns
+        -------
+        int
+            The number of particle types.
+        """
+        call_outputs = self._load_data('getNumParTypes', unitOpId=unitOpId)
+        return int(call_outputs['nParTypes'].value)
+
+    def nsensitivities(self) -> int:
+        """
+        Get the number of sensitivities defined for the simulation.
+
+        Returns
+        -------
+        int
+            The number of sensitivities.
+        """
+        call_outputs = self._load_data('getNumSensitivities')
+        return int(call_outputs['nSens'].value)
+
+    def solution_inlet(
+            self,
+            unitOpId: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the inlet solution for a given unit operation.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the solution.
+        """
+        return self._load_and_process(
+            'getSolutionInlet', unitOpId=unitOpId, own_data=own_data
+        )
 
-        if own_data:
-            return time.copy(), data.copy(), dims
-        else:
-            return time, data, dims
+    def solution_outlet(
+            self,
+            unitOpId: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the outlet solution for a given unit operation.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the solution.
+        """
+        return self._load_and_process(
+            'getSolutionOutlet', unitOpId=unitOpId, 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 solution_bulk(
+            self,
+            unitOpId: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the bulk solution for a given unit operation.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the solution.
+        """
+        return self._load_and_process(
+            'getSolutionBulk', unitOpId=unitOpId, 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 solution_particle(
+            self,
+            unitOpId: int,
+            parType: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the particle-phase solution for a given unit operation and particle type.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        parType : int
+            The particle type index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the solution.
+        """
+        return self._load_and_process(
+            'getSolutionParticle', unitOpId=unitOpId, parType=parType, 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 solution_solid(
+            self,
+            unitOpId: int,
+            parType: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the solid-phase solution for a given unit operation and particle type.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        parType : int
+            The particle type index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the solution.
+        """
+        return self._load_and_process(
+            'getSolutionSolid', unitOpId=unitOpId, parType=parType, 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 solution_flux(
+            self,
+            unitOpId: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the flux solution for a given unit operation.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the solution.
+        """
+        return self._load_and_process(
+            'getSolutionFlux', unitOpId=unitOpId, 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 solution_volume(
+            self,
+            unitOpId: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the volume solution for a given unit operation.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the solution.
+        """
+        return self._load_and_process(
+            'getSolutionVolume', unitOpId=unitOpId, 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 soldot_inlet(
+            self,
+            unitOpId: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the time derivative of the inlet solution for a given unit operation.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the solution derivative.
+        """
+        return self._load_and_process(
+            'getSolutionDerivativeInlet', unitOpId=unitOpId, 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 soldot_outlet(
+            self,
+            unitOpId: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the time derivative of the outlet solution for a given unit operation.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the solution derivative.
+        """
+        return self._load_and_process(
+            'getSolutionDerivativeOutlet', unitOpId=unitOpId, 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 soldot_bulk(
+            self,
+            unitOpId: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the time derivative of the bulk solution for a given unit operation.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the solution derivative.
+        """
+        return self._load_and_process(
+            'getSolutionDerivativeBulk', unitOpId=unitOpId, 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 soldot_particle(
+            self,
+            unitOpId: int,
+            parType: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the time derivative of the particle-phase solution for a given unit operation and particle type.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        parType : int
+            The particle type index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the solution derivative.
+        """
+        return self._load_and_process(
+            'getSolutionDerivativeParticle',
+            unitOpId=unitOpId,
+            parType=parType,
+            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 soldot_solid(
+            self,
+            unitOpId: int,
+            parType: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the time derivative of the solid-phase solution for a given unit operation and particle type.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        parType : int
+            The particle type index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the solution derivative.
+        """
+        return self._load_and_process(
+            'getSolutionDerivativeSolid',
+            unitOpId=unitOpId,
+            parType=parType,
+            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 soldot_flux(
+            self,
+            unitOpId: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the time derivative of the flux solution for a given unit operation.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the solution derivative.
+        """
+        return self._load_and_process(
+            'getSolutionDerivativeFlux', unitOpId=unitOpId, 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 soldot_volume(
+            self,
+            unitOpId: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the time derivative of the volume solution for a given unit operation.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the solution derivative.
+        """
+        return self._load_and_process(
+            'getSolutionDerivativeVolume', unitOpId=unitOpId, 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 sens_inlet(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the sensitivity data for the inlet of a given unit operation and sensitivity index.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity data.
+        """
+        return self._load_and_process(
+            'getSensitivityInlet', 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 sens_outlet(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load sensitivity data for the outlet of a given unit operation and sensitivity index.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity data.
+        """
+        return self._load_and_process(
+            'getSensitivityOutlet',
+            unitOpId=unitOpId,
+            sensIdx=sensIdx,
+            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 sens_bulk(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load sensitivity data for the bulk of a given unit operation and sensitivity index.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity data.
+        """
+        return self._load_and_process(
+            'getSensitivityBulk', unitOpId=unitOpId, sensIdx=sensIdx, 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 sens_particle(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            parType: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load sensitivity data for the particle phase of a given unit operation, sensitivity index, and particle type.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        parType : int
+            The particle type index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity data.
+        """
+        return self._load_and_process(
+            'getSensitivityParticle',
+            unitOpId=unitOpId,
+            sensIdx=sensIdx,
+            parType=parType,
+            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 sens_solid(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            parType: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load sensitivity data for the solid phase of a given unit operation, sensitivity index, and particle type.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        parType : int
+            The particle type index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity data.
+        """
+        return self._load_and_process(
+            'getSensitivitySolid',
+            unitOpId=unitOpId,
+            sensIdx=sensIdx,
+            parType=parType,
+            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 sens_flux(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load sensitivity data for the flux of a given unit operation and sensitivity index.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity data.
+        """
+        return self._load_and_process(
+            'getSensitivityFlux', unitOpId=unitOpId, sensIdx=sensIdx, 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 sens_volume(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load sensitivity data for the volume of a given unit operation and sensitivity index.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity data.
+        """
+        return self._load_and_process(
+            'getSensitivityVolume',
+            unitOpId=unitOpId,
+            sensIdx=sensIdx,
+            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 sensdot_inlet(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load sensitivity derivative data for the inlet of a given unit operation and sensitivity index.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity derivative data.
+        """
+        return self._load_and_process(
+            'getSensitivityDerivativeInlet',
+            unitOpId=unitOpId,
+            sensIdx=sensIdx,
+            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 sensdot_outlet(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load sensitivity derivative data for the outlet of a given unit operation and sensitivity index.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity derivative data.
+        """
+        return self._load_and_process(
+            'getSensitivityDerivativeOutlet',
+            unitOpId,
+            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 sensdot_bulk(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load sensitivity derivative data for the bulk of a given unit operation and sensitivity index.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity derivative data.
+        """
+        return self._load_and_process(
+            'getSensitivityDerivativeBulk',
+            unitOpId=unitOpId,
+            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 sensdot_particle(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            parType: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load sensitivity derivative data for the particle phase of a given unit operation, sensitivity index, and particle type.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        parType : int
+            The particle type index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity derivative data.
+        """
+        return self._load_and_process(
+            'getSensitivityDerivativeParticle',
+            unitOpId=unitOpId,
+            sensIdx=sensIdx,
+            parType=parType,
+            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 sensdot_solid(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            parType: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load sensitivity derivative data for the solid phase of a given unit operation, sensitivity index, and particle type.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        parType : int
+            The particle type index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity derivative data.
+        """
+        return self._load_and_process(
+            'getSensitivityDerivativeSolid',
+            unitOpId=unitOpId,
+            sensIdx=sensIdx,
+            parType=parType,
+            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 sensdot_flux(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load sensitivity derivative data for the flux of a given unit operation and sensitivity index.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity derivative data.
+        """
+        return self._load_and_process(
+            'getSensitivityDerivativeFlux',
+            unitOpId=unitOpId,
+            sensIdx=sensIdx,
+            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 sensdot_volume(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load sensitivity derivative data for the volume of a given unit operation and sensitivity index.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity derivative data.
+        """
+        return self._load_and_process(
+            'getSensitivityDerivativeVolume',
+            unitOpId=unitOpId,
+            sensIdx=sensIdx,
+            own_data=own_data
+        )
+
+    def last_state_y(
+            self,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the last state of the system.
+
+        Parameters
+        ----------
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        numpy.ndarray
+            The state data.
+        """
+        return self._load_and_process_array(
+            'state', 'nStates', 'getLastState', own_data=own_data
+        )
+
+    def last_state_ydot(
+            self,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the time derivative of the last state of the system.
+
+        Parameters
+        ----------
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        numpy.ndarray
+            The time derivative data.
+        """
+        return self._load_and_process_array(
+            'state', 'nStates', 'getLastStateTimeDerivative', own_data=own_data
+        )
+
+    def last_state_y_unit(
+            self,
+            unitOpId: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the last state of a given unit operation.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        numpy.ndarray
+            The state data.
+        """
+        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: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the time derivative of the last state of a given unit operation.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        numpy.ndarray
+            The time derivative data.
+        """
+        return self._load_and_process_array(
+            'state',
+            'nStates',
+            'getLastUnitStateTimeDerivative',
+            unitOpId=unitOpId,
+            own_data=own_data
+        )
+
+    def last_state_sens(
+            self,
+            sensIdx: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the last sensitivity state for a given sensitivity index.
+
+        Parameters
+        ----------
+        sensIdx : int
+            The sensitivity index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        numpy.ndarray
+            The sensitivity state data.
+        """
+        return self._load_and_process_array(
+            'state',
+            'nStates',
+            'getLastSensitivityState',
+            sensIdx=sensIdx,
+            own_data=own_data
+        )
+
+    def last_state_sensdot(
+            self,
+            sensIdx: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the time derivative of the last sensitivity state for a given sensitivity index.
+
+        Parameters
+        ----------
+        sensIdx : int
+            The sensitivity index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        numpy.ndarray
+            The sensitivity state time derivative data.
+        """
+        return self._load_and_process_array(
+            'state',
+            'nStates',
+            'getLastSensitivityStateTimeDerivative',
+            sensIdx=sensIdx,
+            own_data=own_data
+        )
+
+    def last_state_sens_unit(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the last sensitivity state for a given unit operation and sensitivity index.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity state.
+        """
+        return self._load_and_process(
+            'getLastSensitivityUnitState',
+            unitOpId=unitOpId,
+            sensIdx=sensIdx,
+            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 last_state_sensdot_unit(
+            self,
+            unitOpId: int,
+            sensIdx: int,
+            own_data: bool = True,
+            ) -> tuple[numpy.ndarray, numpy.ndarray, list[str]]:
+        """
+        Load the time derivative of the last sensitivity state for a given unit operation and sensitivity index.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        sensIdx : int
+            The sensitivity index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        tuple[numpy.ndarray, numpy.ndarray, list[str]]
+            The time, data, and dimensions of the sensitivity state time derivative.
+        """
+        return self._load_and_process(
+            'getLastSensitivityUnitStateTimeDerivative',
+            unitOpId=unitOpId,
+            sensIdx=sensIdx,
+            own_data=own_data
+        )
+
+    def primary_coordinates(
+            self,
+            unitOpId: int,
+            own_data: bool = True,
+            ) -> numpy.ndarray:
+        """
+        Load the primary coordinates for a given unit operation.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        numpy.ndarray
+            The primary coordinates.
+        """
+        return self._load_and_process_array(
+            'coords',
+            'nCoords',
+            'getPrimaryCoordinates',
+            unitOpId=unitOpId,
+            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 secondary_coordinates(
+            self,
+            unitOpId: int,
+            own_data: bool = True,
+            ) -> numpy.ndarray:
+        """
+        Load the secondary coordinates for a given unit operation.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        numpy.ndarray
+            The secondary coordinates.
+        """
+        return self._load_and_process_array(
+            'coords',
+            'nCoords',
+            'getSecondaryCoordinates',
+            unitOpId=unitOpId,
+            own_data=own_data
+        )
 
+    def particle_coordinates(
+            self,
+            unitOpId: int,
+            parType: int,
+            own_data: bool = True,
+            ) -> numpy.ndarray:
+        """
+        Load the particle coordinates for a given unit operation and particle type.
+
+        Parameters
+        ----------
+        unitOpId : int
+            The unit operation ID.
+        parType : int
+            The particle type index.
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        numpy.ndarray
+            The particle coordinates.
+        """
+        return self._load_and_process_array(
+            'coords',
+            'nCoords',
+            'getParticleCoordinates',
+            unitOpId=unitOpId,
+            parType=parType, own_data=own_data
+        )
 
-class CadetDLL:
+    def solution_times(
+            self,
+            own_data: bool = True,
+            ) -> numpy.ndarray:
+        """
+        Load the solution times from the simulation.
+
+        Parameters
+        ----------
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        numpy.ndarray
+            The solution times.
+        """
+        return self._load_and_process_array(
+            'time', 'nTime', 'getSolutionTimes', own_data=own_data
+        )
 
-    def __init__(self, dll_path):
-        self.cadet_path = dll_path
-        self._lib = ctypes.cdll.LoadLibrary(dll_path)
+    def time_sim(
+            self,
+            own_data: bool = True,
+            ) -> float:
+        """
+        Load the simulation run time.
+
+        Parameters
+        ----------
+        own_data : bool, optional
+            Whether to copy the data (default is True).
+
+        Returns
+        -------
+        float
+            The simulation run time.
+        """
+        call_outputs = self._load_data('getTimeSim')
+        return float(call_outputs['timeSim'].value)
+
+
+class CadetDLLRunner(CadetRunnerBase):
+    """
+    Runner for CADET simulations using a DLL-based interface.
+
+    This class loads and interacts with the CADET DLL, providing methods for running
+    simulations and loading results.
+    """
+
+    def __init__(self, dll_path: os.PathLike | str) -> None:
+        """
+        Initialize the CadetDLLRunner.
+
+        Parameters
+        ----------
+        dll_path : os.PathLike or str
+            Path to the CADET DLL.
+        """
+        self._cadet_path = Path(dll_path)
+        self._lib = ctypes.cdll.LoadLibrary(self._cadet_path.as_posix())
 
         # Query meta information
         cdtGetLibraryVersion = self._lib.cdtGetLibraryVersion
         cdtGetLibraryVersion.restype = ctypes.c_char_p
-        self.cadet_version = cdtGetLibraryVersion().decode('utf-8')
+        self._cadet_version = cdtGetLibraryVersion().decode('utf-8')
 
         cdtGetLibraryCommitHash = self._lib.cdtGetLibraryCommitHash
         cdtGetLibraryCommitHash.restype = ctypes.c_char_p
-        self.cadet_commit_hash = cdtGetLibraryCommitHash().decode('utf-8')
+        self._cadet_commit_hash = cdtGetLibraryCommitHash().decode('utf-8')
 
         cdtGetLibraryBranchRefspec = self._lib.cdtGetLibraryBranchRefspec
         cdtGetLibraryBranchRefspec.restype = ctypes.c_char_p
-        self.cadet_branch = cdtGetLibraryBranchRefspec().decode('utf-8')
+        self._cadet_branch = cdtGetLibraryBranchRefspec().decode('utf-8')
 
         cdtGetLibraryBuildType = self._lib.cdtGetLibraryBuildType
         cdtGetLibraryBuildType.restype = ctypes.c_char_p
-        self.cadet_build_type = cdtGetLibraryBuildType().decode('utf-8')
+        self._cadet_build_type = cdtGetLibraryBuildType().decode('utf-8')
 
-        # Enable logging
-        LOG_HANDLER_CLBK = ctypes.CFUNCTYPE(
+        # Define the log handler callback type
+        self.LOG_HANDLER_CLBK = ctypes.CFUNCTYPE(
             None,
-            ctypes.c_char_p,
-            ctypes.c_char_p,
-            ctypes.c_uint,
-            ctypes.c_int,
-            ctypes.c_char_p,
-            ctypes.c_char_p
-        )
-        set_log_handler = self._lib.cdtSetLogReceiver
-        set_log_handler.argtypes = [LOG_HANDLER_CLBK]
-        set_log_handler.restype = None
-        # Keep reference alive by assigning it to this Python object
-        self._log_handler = LOG_HANDLER_CLBK(log_handler)
-        set_log_handler(self._log_handler)
+            ctypes.c_char_p,  # file
+            ctypes.c_char_p,  # func
+            ctypes.c_uint,    # line
+            ctypes.c_int,     # level
+            ctypes.c_char_p,  # level_name
+            ctypes.c_char_p   # message
+        )
+
+        # Set up the C functions for setting the log handler and log level
+        self._set_log_handler = self._lib.cdtSetLogReceiver
+        self._set_log_handler.argtypes = [self.LOG_HANDLER_CLBK]
+        self._set_log_handler.restype = None
 
         self._set_log_level = self._lib.cdtSetLogLevel
         self._set_log_level.argtypes = [ctypes.c_int]
         self._set_log_level.restype = None
-        self._set_log_level(2)
 
-        # Query API
-        cdtGetAPIv010000 = self._lib.cdtGetAPIv010000
-        cdtGetAPIv010000.argtypes = [ctypes.POINTER(CADETAPIV010000)]
-        cdtGetAPIv010000.restype = c_cadet_result
+        self._default_log_level = 2
 
-        self._api = CADETAPIV010000()
-        cdtGetAPIv010000(ctypes.byref(self._api))
+        # Query API
+        try:
+            cdtGetLatestCAPIVersion = self._lib.cdtGetLatestCAPIVersion#
+        except AttributeError:
+            raise ValueError(
+                "CADET-Python does not support CADET-CAPI at all."
+            )
+        cdtGetLatestCAPIVersion.restype = ctypes.c_char_p
+        self._cadet_capi_version = cdtGetLatestCAPIVersion().decode('utf-8')
+
+        # Check which C-API is provided by CADET (given the current install path)
+        if self._cadet_capi_version == "1.0.0":
+            cdtGetAPIv010000 = self._lib.cdtGetAPIv010000
+            cdtGetAPIv010000.argtypes = [ctypes.POINTER(CADETAPIV010000)]
+            cdtGetAPIv010000.restype = c_cadet_result
+            self._api = CADETAPIV010000()
+            cdtGetAPIv010000(ctypes.byref(self._api))
+        else:
+            raise ValueError(
+                "CADET-Python does not support CADET-CAPI version "
+                f"({self._cadet_capi_version})."
+            )
 
         self._driver = self._api.createDriver()
+        self.res: Optional[SimulationResult] = None
+
+    def clear(self) -> None:
+        """
+        Clear the current simulation state.
 
-    def clear(self):
+        This method deletes the current simulation results and resets the driver.
+        """
         if hasattr(self, "res"):
             del self.res
+
         self._api.deleteDriver(self._driver)
         self._driver = self._api.createDriver()
 
-    def __del__(self):
-        log_print('deleteDriver()')
+    def __del__(self) -> None:
+        """
+        Clean up the CADET driver on object deletion.
+        """
         self._api.deleteDriver(self._driver)
 
-
-    def run(self, filename = None, simulation=None, timeout = None, check=None):
+    def setup_log_buffer(self, log_level: int = None) -> io.StringIO:
+        """
+        Set up a new log buffer for capturing log output with the specified log level.
+
+        Parameters
+        ----------
+        log_level : int, optional
+            The desired log level. Defaults to the instance's default log level.
+
+        Returns
+        -------
+        log_buffer : io.StringIO
+            A new log buffer for capturing log output.
+        """
+        if log_level is None:
+            log_level = self._default_log_level
+
+        # Create a new log buffer
+        log_buffer = io.StringIO()
+
+        def log_handler(file, func, line, level, level_name, message):
+            """Logging callback function."""
+            msg = f"{level_name.decode('utf-8')} ({func.decode('utf-8')}:{line}) {message.decode('utf-8')}"
+            log_buffer.write(msg + "\n")  # Write to the instance's buffer
+
+        # Set up the logging callback
+        self._log_handler = self.LOG_HANDLER_CLBK(log_handler)
+        self._set_log_handler(self._log_handler)
+
+        # Set the log level
+        self._set_log_level(log_level)
+
+        # Return the log buffer so it can be accessed after the run
+        return log_buffer
+
+    def run(
+            self,
+            simulation: Optional["Cadet"] = None,
+            timeout: Optional[int] = None,
+            ) -> None:
+        """
+        Run a CADET simulation using the DLL interface.
+
+        Parameters
+        ----------
+        simulation : Optional[Cadet]
+            Simulation object containing input data.
+        timeout : Optional[int]
+            Maximum time allowed for the simulation to run, in seconds.
+
+        Raises
+        ------
+        RuntimeError
+            If the simulation process returns a non-zero exit code.
+        """
         pp = cadet_dll_parameterprovider.PARAMETERPROVIDER(simulation)
 
-        self._api.runSimulation(self._driver, ctypes.byref(pp))
+        log_buffer = self.setup_log_buffer()
+
+        returncode = self._api.runSimulation(self._driver, ctypes.byref(pp))
+
+        if returncode != 0:
+            error_message = log_buffer.getvalue()
+            raise RuntimeError(f"Simulation failed with error: {error_message}")
+
         self.res = SimulationResult(self._api, self._driver)
-        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
+
+        return
+
+    def load_results(self, sim: "Cadet") -> None:
+        """
+        Load the simulation results into the provided simulation object.
+
+        Parameters
+        ----------
+        sim : Cadet
+            The simulation object where results will be loaded.
+        """
+        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)
+
+        self.load_meta(sim)
+
+    def load_solution_times(self, sim: "Cadet") -> None:
+        """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: "Cadet") -> None:
+        """Load coordinates data from simulation results."""
+        coordinates = addict.Dict()
+        for unit in range(self.res.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: "Cadet") -> addict.Dict:
+        """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)
+        for unit in range(self.res.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(solution) > 1:
+            sim.root.output.solution.update(solution)
 
-                        if not len(solution.solution_times):
-                            solution.solution_times = t
+        return solution
 
-                        solution[key][solution_str] = out
+    def load_sensitivity(self, sim: "Cadet") -> None:
+        """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(self.res.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: "Cadet") -> None:
+        """Load last state from simulation results."""
+        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()
+
+        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)
+
+        solution = sim.root.output.solution
+        for unit in range(self.res.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
+
+    @staticmethod
+    def _get_index_string(prefix: str, index: int) -> str:
+        """Get a formatted string index (e.g., ('unit', 0) -> 'unit_000')."""
+        return f'{prefix}_{index:03d}'
+
+    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)
+            if data is None:
+                return
+            solution = func(self, sim, data, unitOpId, solution_str, *args, **kwargs)
+            return solution
+        return wrapper
+
+    @_checks_if_write_is_true
+    @_loads_data
+    def _load_solution_trivial(
+            self,
+            sim: "Cadet",
+            data: tuple[numpy.ndarray, numpy.ndarray, list[int]],
+            unitOpId: int,
+            solution_str: str,
+            sensIdx: Optional[int] = None
+            ) -> addict.Dict:
+        """Load trivial solution data from simulation results."""
+        solution = addict.Dict()
+        _, out, _ = data
+        solution[solution_str] = out
         return solution
 
-    def load_solution_io(self, sim, solution_fun, solution_str):
+    @_checks_if_write_is_true
+    @_loads_data
+    def _load_solution_io(
+            self,
+            sim: "Cadet",
+            data: tuple[numpy.ndarray, numpy.ndarray, list[int]],
+            unitOpId: int,
+            solution_str: str,
+            sensIdx: Optional[int] = None
+            ) -> addict.Dict:
+        """Load IO-related 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)
-
-                        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
+        _, 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
+
+        return solution
+
+    @_checks_if_write_is_true
+    def _load_solution_particle(
+            self,
+            sim: "Cadet",
+            unitOpId: int,
+            solution_str: str,
+            sensIdx: Optional[int] = None
+            ) -> addict.Dict:
+        """Load particle-related solution data from simulation results."""
+        solution = addict.Dict()
+        solution_fun = getattr(self.res, solution_str)
+
+        npartype = self.res.npartypes(unitOpId)
+
+        for partype in range(npartype):
+            if sensIdx is None:
+                data = solution_fun(unitOpId, partype)
+            else:
+                data = solution_fun(unitOpId, sensIdx, partype)
+
+            if data is None:
+                continue
+
+            _, out, dims = data
+
+            if npartype == 1:
+                solution[solution_str] = out
+            else:
+                par_index = self._get_index_string('partype', partype)
+                solution[f'{solution_str}_{par_index}'] = out
+
+        if len(solution) == 0:
+            return addict.Dict()
+
         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_derivative_flux(self, sim):
-        return self.load_solution(sim, self.res.derivativeFlux, 'soldot_flux')
-
-    def load_derivative_volume(self, sim):
-        return self.load_solution(sim, self.res.derivativeVolume, 'soldot_volume')
-
-    def load_sensitivity_inlet(self, sim):
-        return self.load_solution_io(sim, self.res.sensitivityInlet, 'sens_inlet')
-
-    def load_sensitivity_outlet(self, sim):
-        return self.load_solution_io(sim, self.res.sensitivityOutlet, 'sens_outlet')
-
-    def load_sensitivity_bulk(self, sim):
-        return self.load_solution(sim, self.res.sensitivityBulk, 'sens_bulk')
-
-    def load_sensitivity_particle(self, sim):
-        return self.load_solution(sim, self.res.sensitivityParticle, 'sens_particle')
-
-    def load_sensitivity_solid(self, sim):
-        return self.load_solution(sim, self.res.sensitivitySolid, 'sens_solid')
-
-    def load_sensitivity_flux(self, sim):
-        return self.load_solution(sim, self.res.sensitivityFlux, 'sens_flux')
-
-    def load_sensitivity_volume(self, sim):
-        return self.load_solution(sim, self.res.sensitivityVolume, 'sens_volume')
-
-    def load_sensitivity_derivative_inlet(self, sim):
-        return self.load_solution_io(sim, self.res.sensitivityDerivativeInlet, 'sensdot_inlet')
-
-    def load_sensitivity_derivative_outlet(self, sim):
-        return self.load_solution_io(sim, self.res.sensitivityDerivativeOutlet, 'sensdot_outlet')
-
-    def load_sensitivity_derivative_bulk(self, sim):
-        return self.load_solution(sim, self.res.sensitivityDerivativeBulk, 'sensdot_bulk')
-
-    def load_sensitivity_derivative_particle(self, sim):
-        return self.load_solution(sim, self.res.sensitivityDerivativeParticle, 'sensdot_particle')
-
-    def load_sensitivity_derivative_solid(self, sim):
-        return self.load_solution(sim, self.res.sensitivityDerivativeSolid, 'sensdot_solid')
-
-    def load_sensitivity_derivative_flux(self, sim):
-        return self.load_solution(sim, self.res.sensitivityDerivativeFlux, 'sensdot_flux')
-
-    def load_sensitivity_derivative_volume(self, sim):
-        return self.load_solution(sim, self.res.sensitivityDerivativeVolume, 'sensdot_volume')
-
-    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}'
+
+    def load_meta(self, sim: "Cadet") -> None:
+        """Load meta information about CADET."""
+
+        sim.root.meta.cadet_branch = self.cadet_branch
+        sim.root.meta.cadet_commit = self.cadet_commit_hash
+        sim.root.meta.cadet_version = self.cadet_version
+        sim.root.meta.file_format = self.res.file_format()
+        sim.root.meta.time_sim = self.res.time_sim()
+
+    @property
+    def cadet_version(self) -> str:
+        return self._cadet_version
+
+    @property
+    def cadet_branch(self) -> str:
+        return self._cadet_branch
+
+    @property
+    def cadet_build_type(self) -> str:
+        return self._cadet_build_type
+
+    @property
+    def cadet_commit_hash(self) -> str:
+        return self._cadet_commit_hash
+
+    @property
+    def cadet_path(self) -> os.PathLike:
+        return self._cadet_path
diff --git a/cadet/cadet_dll_parameterprovider.py b/cadet/cadet_dll_parameterprovider.py
index 73e03d4..9b95153 100644
--- a/cadet/cadet_dll_parameterprovider.py
+++ b/cadet/cadet_dll_parameterprovider.py
@@ -1,47 +1,88 @@
-
 import ctypes
 import cadet.cadet_dll_utils as utils
 import addict
+from typing import Any, Dict, Optional, Union
 
-c_cadet_result = ctypes.c_int
 
+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 null(*args: Any) -> None:
     pass
 
+
 if 0:
     log_print = print
 else:
     log_print = null
 
+
 class NestedDictReader:
+    """
+    Utility class to read and navigate through nested dictionaries.
+    """
 
-    def __init__(self, data):
+    def __init__(self, data: Dict[str, Any]) -> None:
         self._root = data
-        self._cursor = []
-        self._cursor.append(data)
-        self.buffer = None
-
-    def push_scope(self, scope):
+        self._cursor = [data]
+        self.buffer: Optional[Any] = None
+
+    def push_scope(self, scope: str) -> bool:
+        """
+        Enter a nested scope within the dictionary.
+
+        Parameters
+        ----------
+        scope : str
+            The key representing the nested scope.
+
+        Returns
+        -------
+        bool
+            True if the scope exists and was entered, otherwise False.
+        """
         if scope in self._cursor[-1]:
-            log_print('Entering scope {}'.format(scope))
+            log_print(f'Entering scope {scope}')
             self._cursor.append(self._cursor[-1][scope])
             return True
-
         return False
 
-    def pop_scope(self):
-        self._cursor.pop()
-        log_print('Exiting scope')
-
-    def current(self):
+    def pop_scope(self) -> None:
+        """
+        Exit the current scope.
+        """
+        if len(self._cursor) > 1:
+            self._cursor.pop()
+            log_print('Exiting scope')
+
+    def current(self) -> Any:
+        """
+        Get the current scope data.
+
+        Returns
+        -------
+        Any
+            The current data under the scope.
+        """
         return self._cursor[-1]
 
-def recursively_convert_dict(data):
+
+def recursively_convert_dict(data: Dict[str, Any]) -> addict.Dict:
+    """
+    Recursively convert dictionary keys to uppercase while preserving nested structure.
+
+    Parameters
+    ----------
+    data : dict
+        The dictionary to convert.
+
+    Returns
+    -------
+    addict.Dict
+        A new dictionary with all keys converted to uppercase.
+    """
     ans = addict.Dict()
     for key_original, item in data.items():
         if isinstance(item, dict):
@@ -53,20 +94,23 @@ def recursively_convert_dict(data):
 
 
 class PARAMETERPROVIDER(ctypes.Structure):
-    """Implemented the CADET Parameter Provider interface which allows 
-    querying python for the necessary parameters for a CADET simulation
-    
-    _fields_ is a list of function names and signatures exposed by the
-    capi Parameter Provider interface
-    https://docs.python.org/3/library/ctypes.html
-    https://github.com/modsim/CADET/blob/master/src/libcadet/api/CAPIv1.cpp"""
+    """
+    Implement the CADET Parameter Provider interface, allowing querying Python for parameters.
 
-    def __init__(self, simulation):
-        sim_input = recursively_convert_dict(simulation)
+    This class exposes various function pointers as fields in a ctypes structure
+    to be used with CADET's C-API.
 
+    Parameters
+    ----------
+    simulation : Cadet
+        The simulation object containing the input data.
+    """
+
+    def __init__(self, simulation: "Cadet") -> None:
+        sim_input = recursively_convert_dict(simulation.root.input)
         self.userData = NestedDictReader(sim_input)
 
-        #figure out how to add this to the class
+        # Assign function pointers at instance level
         self.getDouble = self._fields_[1][1](utils.param_provider_get_double)
         self.getInt = self._fields_[2][1](utils.param_provider_get_int)
         self.getBool = self._fields_[3][1](utils.param_provider_get_bool)
@@ -112,25 +156,3 @@ def __init__(self, simulation):
         ('pushScope', ctypes.CFUNCTYPE(c_cadet_result, ctypes.py_object, ctypes.c_char_p)),
         ('popScope', ctypes.CFUNCTYPE(c_cadet_result, ctypes.py_object)),
     ]
-
-    #these don't work right now but work when places on an instance
-    #getDouble = _fields_[1][1](utils.param_provider_get_double)
-    #getInt = _fields_[2][1](utils.param_provider_get_int)
-    #getBool = _fields_[3][1](utils.param_provider_get_bool)
-    #getString = _fields_[4][1](utils.param_provider_get_string)
-
-    #getDoubleArray = _fields_[5][1](utils.param_provider_get_double_array)
-    #getIntArray = _fields_[6][1](utils.param_provider_get_int_array)
-    #getBoolArray = ctypes.cast(None, _fields_[7][1])
-    #getStringArray = ctypes.cast(None, _fields_[8][1])
-
-    #getDoubleArrayItem = _fields_[9][1](utils.param_provider_get_double_array_item)
-    #getIntArrayItem = _fields_[10][1](utils.param_provider_get_int_array_item)
-    #getBoolArrayItem = _fields_[11][1](utils.param_provider_get_bool_array_item)
-    #getStringArrayItem = _fields_[12][1](utils.param_provider_get_string_array_item)
-
-    #exists = _fields_[13][1](utils.param_provider_exists)
-    #isArray = _fields_[14][1](utils.param_provider_is_array)
-    #numElements = _fields_[15][1](utils.param_provider_num_elements)
-    #pushScope = _fields_[16][1](utils.param_provider_push_scope)
-    #popScope = _fields_[17][1](utils.param_provider_pop_scope)
\ No newline at end of file
diff --git a/cadet/cadet_dll_utils.py b/cadet/cadet_dll_utils.py
index d4b4045..daa1933 100644
--- a/cadet/cadet_dll_utils.py
+++ b/cadet/cadet_dll_utils.py
@@ -1,16 +1,36 @@
 import ctypes
-import numpy
+import numpy as np
+from typing import Any, Optional
 
-def null(*args):
-    pass
-
-if 0:
-    log_print = print
-else:
-    log_print = null
 
+def null(*args: Any) -> None:
+    """Do nothing (used as a placeholder function)."""
+    pass
 
-def param_provider_get_double(reader, name, val):
+log_print = print if 0 else null
+
+def param_provider_get_double(
+        reader: Any,
+        name: ctypes.c_char_p,
+        val: ctypes.POINTER(ctypes.c_double)
+        ) -> int:
+    """
+    Retrieve a double value from the reader based on the provided name.
+
+    Parameters
+    ----------
+    reader : Any
+        The reader object containing the current data scope.
+    name : ctypes.c_char_p
+        The name of the parameter to retrieve.
+    val : ctypes.POINTER(ctypes.c_double)
+        A pointer to store the retrieved double value.
+
+    Returns
+    -------
+    int
+        0 if the value was found and retrieved successfully, -1 otherwise.
+    """
     n = name.decode('utf-8')
     c = reader.current()
 
@@ -22,14 +42,34 @@ def param_provider_get_double(reader, name, val):
             float_val = float(o[0])
 
         val[0] = ctypes.c_double(float_val)
-
-        log_print('GET scalar [double] {}: {}'.format(n, float(val[0])))
+        log_print(f"GET scalar [double] {n}: {float(val[0])}")
         return 0
 
     return -1
 
 
-def param_provider_get_int(reader, name, val):
+def param_provider_get_int(
+        reader: Any,
+        name: ctypes.c_char_p,
+        val: ctypes.POINTER(ctypes.c_int)
+        ) -> int:
+    """
+    Retrieve an integer value from the reader based on the provided name.
+
+    Parameters
+    ----------
+    reader : Any
+        The reader object containing the current data scope.
+    name : ctypes.c_char_p
+        The name of the parameter to retrieve.
+    val : ctypes.POINTER(ctypes.c_int)
+        A pointer to store the retrieved integer value.
+
+    Returns
+    -------
+    int
+        0 if the value was found and retrieved successfully, -1 otherwise.
+    """
     n = name.decode('utf-8')
     c = reader.current()
 
@@ -41,41 +81,79 @@ def param_provider_get_int(reader, name, val):
             int_val = int(o[0])
 
         val[0] = ctypes.c_int(int_val)
-
-        log_print('GET scalar [int] {}: {}'.format(n, int(val[0])))
+        log_print(f"GET scalar [int] {n}: {int(val[0])}")
         return 0
 
     return -1
 
 
-def param_provider_get_bool(reader, name, val):
+def param_provider_get_bool(
+        reader: Any,
+        name: ctypes.c_char_p,
+        val: ctypes.POINTER(ctypes.c_uint8)
+        ) -> int:
+    """
+    Retrieve a boolean value from the reader based on the provided name.
+
+    Parameters
+    ----------
+    reader : Any
+        The reader object containing the current data scope.
+    name : ctypes.c_char_p
+        The name of the parameter to retrieve.
+    val : ctypes.POINTER(ctypes.c_uint8)
+        A pointer to store the retrieved boolean value.
+
+    Returns
+    -------
+    int
+        0 if the value was found and retrieved successfully, -1 otherwise.
+    """
     n = name.decode('utf-8')
     c = reader.current()
 
     if n in c:
         o = c[n]
-
         try:
             int_val = int(o)
         except TypeError:
             int_val = int(o[0])
 
         val[0] = ctypes.c_uint8(int_val)
-
-        log_print('GET scalar [bool] {}: {}'.format(n, bool(val[0])))
+        log_print(f"GET scalar [bool] {n}: {bool(val[0])}")
         return 0
 
     return -1
 
 
-def param_provider_get_string(reader, name, val):
+def param_provider_get_string(
+        reader: Any,
+        name: ctypes.c_char_p,
+        val: ctypes.POINTER(ctypes.c_char_p)
+        ) -> int:
+    """
+    Retrieve a string value from the reader based on the provided name.
+
+    Parameters
+    ----------
+    reader : Any
+        The reader object containing the current data scope.
+    name : ctypes.c_char_p
+        The name of the parameter to retrieve.
+    val : ctypes.POINTER(ctypes.c_char_p)
+        A pointer to store the retrieved string value.
+
+    Returns
+    -------
+    int
+        0 if the value was found and retrieved successfully, -1 otherwise.
+    """
     n = name.decode('utf-8')
     c = reader.current()
 
     if n in c:
         o = c[n]
 
-        #we have one of array of strings, array of bytestrings, bytestring or string or something convertable to one of these
         if hasattr(o, 'encode'):
             bytes_val = o.encode('utf-8')
         elif hasattr(o, 'decode'):
@@ -87,51 +165,121 @@ def param_provider_get_string(reader, name, val):
 
         reader.buffer = bytes_val
         val[0] = ctypes.cast(reader.buffer, ctypes.c_char_p)
-
         return 0
 
     return -1
 
 
-def param_provider_get_double_array(reader, name, n_elem, val):
+def param_provider_get_double_array(
+        reader: Any,
+        name: ctypes.c_char_p,
+        n_elem: ctypes.POINTER(ctypes.c_int),
+        val: ctypes.POINTER(ctypes.POINTER(ctypes.c_double))
+        ) -> int:
+    """
+    Retrieve a double array from the reader based on the provided name.
+
+    Parameters
+    ----------
+    reader : Any
+        The reader object containing the current data scope.
+    name : ctypes.c_char_p
+        The name of the parameter to retrieve.
+    n_elem : ctypes.POINTER(ctypes.c_int)
+        A pointer to store the number of elements in the array.
+    val : ctypes.POINTER(ctypes.POINTER(ctypes.c_double))
+        A pointer to store the retrieved array.
+
+    Returns
+    -------
+    int
+        0 if the array was found and retrieved successfully, -1 otherwise.
+    """
     n = name.decode('utf-8')
     c = reader.current()
 
-    if (n in c):
+    if n in c:
         o = c[n]
         if isinstance(o, list):
-            o = numpy.ascontiguousarray(o)
-        if (not isinstance(o, numpy.ndarray)) or (o.dtype != numpy.double) or (not o.flags.c_contiguous):
+            o = np.ascontiguousarray(o)
+        if not isinstance(o, np.ndarray) or o.dtype != np.double or not o.flags.c_contiguous:
             return -1
 
         n_elem[0] = ctypes.c_int(o.size)
         val[0] = o.ctypes.data_as(ctypes.POINTER(ctypes.c_double))
-        log_print('GET array [double] {}: {}'.format(n, o))
+        log_print(f"GET array [double] {n}: {o}")
         return 0
 
     return -1
 
 
-def param_provider_get_int_array(reader, name, n_elem, val):
+def param_provider_get_int_array(
+        reader: Any,
+        name: ctypes.c_char_p,
+        n_elem: ctypes.POINTER(ctypes.c_int),
+        val: ctypes.POINTER(ctypes.POINTER(ctypes.c_int))
+        ) -> int:
+    """
+    Retrieve an integer array from the reader based on the provided name.
+
+    Parameters
+    ----------
+    reader : Any
+        The reader object containing the current data scope.
+    name : ctypes.c_char_p
+        The name of the parameter to retrieve.
+    n_elem : ctypes.POINTER(ctypes.c_int)
+        A pointer to store the number of elements in the array.
+    val : ctypes.POINTER(ctypes.POINTER(ctypes.c_int))
+        A pointer to store the retrieved array.
+
+    Returns
+    -------
+    int
+        0 if the array was found and retrieved successfully, -1 otherwise.
+    """
     n = name.decode('utf-8')
     c = reader.current()
 
-    if (n in c):
+    if n in c:
         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):
+            o = np.ascontiguousarray(o)
+        if not isinstance(o, np.ndarray) or o.dtype != int or not o.flags.c_contiguous:
             return -1
 
         n_elem[0] = ctypes.c_int(o.size)
         val[0] = o.ctypes.data_as(ctypes.POINTER(ctypes.c_int))
-        log_print('GET array [int] {}: {}'.format(n, o))
+        log_print(f"GET array [int] {n}: {o}")
         return 0
 
     return -1
 
 
-def param_provider_get_double_array_item(reader, name, index, val):
+def param_provider_get_double_array_item(
+        reader: Any,
+        name: ctypes.c_char_p,
+        index: int, val: ctypes.POINTER(ctypes.c_double)
+        ) -> int:
+    """
+    Retrieve an item from a double array in the reader based on the provided name and index.
+
+    Parameters
+    ----------
+    reader : Any
+        The reader object containing the current data scope.
+    name : ctypes.c_char_p
+        The name of the parameter to retrieve.
+    index : int
+        The index of the array item to retrieve.
+    val : ctypes.POINTER(ctypes.c_double)
+        A pointer to store the retrieved double value.
+
+    Returns
+    -------
+    int
+        0 if the value was found and retrieved successfully, -1 otherwise.
+    """
     n = name.decode('utf-8')
     c = reader.current()
 
@@ -144,13 +292,36 @@ def param_provider_get_double_array_item(reader, name, index, val):
             float_val = float(o[index])
 
         val[0] = ctypes.c_double(float_val)
-        log_print('GET array [double] ({}) {}: {}'.format(index, n, val[0]))
+        log_print(f"GET array [double] ({index}) {n}: {val[0]}")
         return 0
 
     return -1
 
 
-def param_provider_get_int_array_item(reader, name, index, val):
+def param_provider_get_int_array_item(
+        reader: Any,
+        name: ctypes.c_char_p,
+        index: int, val: ctypes.POINTER(ctypes.c_int)
+        ) -> int:
+    """
+    Retrieve an item from an integer array in the reader based on the provided name and index.
+
+    Parameters
+    ----------
+    reader : Any
+        The reader object containing the current data scope.
+    name : ctypes.c_char_p
+        The name of the parameter to retrieve.
+    index : int
+        The index of the array item to retrieve.
+    val : ctypes.POINTER(ctypes.c_int)
+        A pointer to store the retrieved integer value.
+
+    Returns
+    -------
+    int
+        0 if the value was found and retrieved successfully, -1 otherwise.
+    """
     n = name.decode('utf-8')
     c = reader.current()
 
@@ -163,13 +334,36 @@ def param_provider_get_int_array_item(reader, name, index, val):
             int_val = int(o[index])
 
         val[0] = ctypes.c_int(int_val)
-        log_print('GET array [int] ({}) {}: {}'.format(index, n, val[0]))
+        log_print(f"GET array [int] ({index}) {n}: {val[0]}")
         return 0
 
     return -1
 
 
-def param_provider_get_bool_array_item(reader, name, index, val):
+def param_provider_get_bool_array_item(
+        reader: Any,
+        name: ctypes.c_char_p,
+        index: int, val: ctypes.POINTER(ctypes.c_uint8)
+        ) -> int:
+    """
+    Retrieve an item from a boolean array in the reader based on the provided name and index.
+
+    Parameters
+    ----------
+    reader : Any
+        The reader object containing the current data scope.
+    name : ctypes.c_char_p
+        The name of the parameter to retrieve.
+    index : int
+        The index of the array item to retrieve.
+    val : ctypes.POINTER(ctypes.c_uint8)
+        A pointer to store the retrieved boolean value.
+
+    Returns
+    -------
+    int
+        0 if the value was found and retrieved successfully, -1 otherwise.
+    """
     n = name.decode('utf-8')
     c = reader.current()
 
@@ -182,50 +376,103 @@ def param_provider_get_bool_array_item(reader, name, index, val):
             int_val = int(o[index])
 
         val[0] = ctypes.c_uint8(int_val)
-        log_print('GET array [bool] ({}) {}: {}'.format(index, n, bool(val[0])))
+        log_print(f"GET array [bool] ({index}) {n}: {bool(val[0])}")
         return 0
 
     return -1
 
 
-def param_provider_get_string_array_item(reader, name, index, val):
-    n = name.decode('utf-8')
-    c = 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
+def param_provider_get_string_array_item(
+        reader: Any,
+        name: ctypes.c_char_p,
+        index: int, val: ctypes.POINTER(ctypes.c_char_p)
+        ) -> int:
+    """
+    Retrieve an item from a string array in the reader based on the provided name and index.
+
+    Parameters
+    ----------
+    reader : Any
+        The reader object containing the current data scope.
+    name : ctypes.c_char_p
+        The name of the parameter to retrieve.
+    index : int
+        The index of the array item to retrieve.
+    val : ctypes.POINTER(ctypes.c_char_p)
+        A pointer to store the retrieved string value.
+
+    Returns
+    -------
+    int
+        0 if the value was found and retrieved successfully, -1 otherwise.
+    """
+    name_str = name.decode('utf-8')
+    current_reader = reader.current()
+
+    if name_str in current_reader:
+        str_value = current_reader[name_str]
+        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_str}: {reader.buffer.decode('utf-8')}")
         return 0
 
     return -1
 
 
-def param_provider_exists(reader, name):
+def param_provider_exists(
+        reader: Any,
+        name: ctypes.c_char_p
+        ) -> int:
+    """
+    Check if a given parameter name exists in the reader.
+
+    Parameters
+    ----------
+    reader : Any
+        The reader object containing the current data scope.
+    name : ctypes.c_char_p
+        The name of the parameter to check.
+
+    Returns
+    -------
+    int
+        1 if the name exists, 0 otherwise.
+    """
     n = name.decode('utf-8')
     c = reader.current()
 
-    log_print('EXISTS {}: {}'.format(n, n in c))
-
-    if n in c:
-        return 1
-
-    return 0
-
-
-def param_provider_is_array(reader, name, res):
+    log_print(f"EXISTS {n}: {n in c}")
+
+    return 1 if n in c else 0
+
+
+def param_provider_is_array(
+        reader: Any,
+        name: ctypes.c_char_p,
+        res: ctypes.POINTER(ctypes.c_uint8)
+        ) -> int:
+    """
+    Check if a given parameter is an array.
+
+    Parameters
+    ----------
+    reader : Any
+        The reader object containing the current data scope.
+    name : ctypes.c_char_p
+        The name of the parameter to check.
+    res : ctypes.POINTER(ctypes.c_uint8)
+        A pointer to store the result (1 if the parameter is an array, 0 otherwise).
+
+    Returns
+    -------
+    int
+        0 if the check was successful, -1 if the parameter does not exist.
+    """
     n = name.decode('utf-8')
     c = reader.current()
 
@@ -233,18 +480,31 @@ def param_provider_is_array(reader, name, res):
         return -1
 
     o = c[n]
-    res[0] = ctypes.c_uint8(0)
-    if isinstance(o, list):
-        res[0] = ctypes.c_uint8(1)
-    elif isinstance(o, numpy.ndarray):
-        res[0] = ctypes.c_uint8(1)
-
-    log_print('ISARRAY {}: {}'.format(n, bool(res[0])))
+    res[0] = ctypes.c_uint8(1 if isinstance(o, (list, np.ndarray)) else 0)
+    log_print(f"ISARRAY {n}: {bool(res[0])}")
 
     return 0
 
 
-def param_provider_num_elements(reader, name):
+def param_provider_num_elements(
+        reader: Any,
+        name: ctypes.c_char_p
+        ) -> int:
+    """
+    Get the number of elements in a given parameter if it is an array.
+
+    Parameters
+    ----------
+    reader : Any
+        The reader object containing the current data scope.
+    name : ctypes.c_char_p
+        The name of the parameter to check.
+
+    Returns
+    -------
+    int
+        The number of elements if the parameter is an array, 1 otherwise.
+    """
     n = name.decode('utf-8')
     c = reader.current()
 
@@ -253,24 +513,56 @@ def param_provider_num_elements(reader, name):
 
     o = c[n]
     if isinstance(o, list):
-        log_print('NUMELEMENTS {}: {}'.format(n, len(o)))
+        log_print(f"NUMELEMENTS {n}: {len(o)}")
         return len(o)
-    elif isinstance(o, numpy.ndarray):
-        log_print('NUMELEMENTS {}: {}'.format(n, o.size))
+    elif isinstance(o, np.ndarray):
+        log_print(f"NUMELEMENTS {n}: {o.size}")
         return o.size
 
-    log_print('NUMELEMENTS {}: {}'.format(n, 1))
+    log_print(f"NUMELEMENTS {n}: 1")
     return 1
 
-def param_provider_push_scope(reader, name):
-	n = name.decode('utf-8')
-	
-	if reader.push_scope(n):
-		return 0
-	else:
-		return -1
 
+def param_provider_push_scope(
+        reader: Any,
+        name: ctypes.c_char_p
+        ) -> int:
+    """
+    Push a new scope in the reader based on the provided name.
+
+    Parameters
+    ----------
+    reader : Any
+        The reader object containing the current data scope.
+    name : ctypes.c_char_p
+        The name of the scope to push.
+
+    Returns
+    -------
+    int
+        0 if the scope was successfully pushed, -1 otherwise.
+    """
+    n = name.decode('utf-8')
 
-def param_provider_pop_scope(reader):
-	reader.pop_scope()
-	return 0
\ No newline at end of file
+    if reader.push_scope(n):
+        return 0
+    else:
+        return -1
+
+
+def param_provider_pop_scope(reader: Any) -> int:
+    """
+    Pop the current scope from the reader.
+
+    Parameters
+    ----------
+    reader : Any
+        The reader object containing the current data scope.
+
+    Returns
+    -------
+    int
+        0 if the scope was successfully popped.
+    """
+    reader.pop_scope()
+    return 0
diff --git a/cadet/h5.py b/cadet/h5.py
new file mode 100644
index 0000000..a4cfa88
--- /dev/null
+++ b/cadet/h5.py
@@ -0,0 +1,432 @@
+from typing import Optional, Union, Any, List
+import copy
+import json
+import pprint
+import warnings
+from pathlib import Path
+
+from addict import Dict
+with warnings.catch_warnings():
+    warnings.filterwarnings("ignore", category=FutureWarning)
+    import h5py
+import numpy
+
+import filelock
+import contextlib
+
+
+class H5:
+    """
+    A class for handling hierarchical HDF5 data structures and JSON representations.
+
+    Attributes
+    ----------
+    root : Dict
+        The root data structure holding the HDF5/JSON data.
+    filename : Optional[str]
+        Path to the HDF5 file.
+
+    Methods
+    -------
+    transform(x: Any) -> Any
+        Applies a transformation to the data before saving.
+    inverse_transform(x: Any) -> Any
+        Applies an inverse transformation to the data after loading.
+    load(paths: Optional[List[str]] = None, update: bool = False, lock: bool = False) -> None
+        Loads data from the specified HDF5 file.
+    save(lock: bool = False) -> None
+        Saves the current data to the specified HDF5 file.
+    save_json(filename: Union[str, Path]) -> None
+        Saves the current data to a JSON file.
+    load_json(filename: Union[str, Path], update: bool = False) -> None
+        Loads data from a JSON file.
+    append(lock: bool = False) -> None
+        Appends new keys to the HDF5 file without reading existing data.
+    update(other: "H5") -> None
+        Merges another H5 object's data with the current one.
+    """
+
+    pp = pprint.PrettyPrinter(indent=4)
+
+    def transform(self, x: Any) -> Any:
+        """
+        Transform the data before saving.
+
+        Parameters
+        ----------
+        x : Any
+            Data to be transformed.
+
+        Returns
+        -------
+        Any
+            Transformed data.
+        """
+        return x
+
+    def inverse_transform(self, x: Any) -> Any:
+        """
+        Apply an inverse transformation to the data after loading.
+
+        Parameters
+        ----------
+        x : Any
+            Data to be transformed back.
+
+        Returns
+        -------
+        Any
+            Inversely transformed data.
+        """
+        return x
+
+    def __init__(self, *data: Any):
+        """
+        Initialize an H5 object with optional data.
+
+        Parameters
+        ----------
+        data : Any
+            Optional initial data to populate the object.
+        """
+        self.root = Dict()
+        self.filename: Optional[str] = None
+        for i in data:
+            self.root.update(copy.deepcopy(i))
+
+    def load(self, paths: Optional[List[str]] = None, update: bool = False, lock: bool = False) -> None:
+        """
+        Load data from the specified HDF5 file.
+
+        Parameters
+        ----------
+        paths : Optional[List[str]], optional
+            Specific paths to load within the HDF5 file.
+        update : bool, optional
+            If True, updates the existing data with the loaded data.
+        lock : bool, optional
+            If True, uses a file lock while loading.
+        """
+        if self.filename is not None:
+            lock_file = filelock.FileLock(self.filename + '.lock') if lock else contextlib.nullcontext()
+
+            with lock_file:
+                with h5py.File(self.filename, 'r') as h5file:
+                    data = Dict(recursively_load(h5file, '/', self.inverse_transform, paths))
+                    if update:
+                        self.root.update(data)
+                    else:
+                        self.root = data
+        else:
+            print('Filename must be set before load can be used')
+
+    def save(self, lock: bool = False) -> None:
+        """
+        Save the current data to the specified HDF5 file.
+
+        Parameters
+        ----------
+        lock : bool, optional
+            If True, uses a file lock while saving.
+        """
+        if self.filename is not None:
+            lock_file = filelock.FileLock(self.filename + '.lock') if lock else contextlib.nullcontext()
+
+            with lock_file:
+                with h5py.File(self.filename, 'w') as h5file:
+                    recursively_save(h5file, '/', self.root, self.transform)
+        else:
+            print("Filename must be set before save can be used")
+
+    def save_json(self, filename: str | Path) -> None:
+        """
+        Save the current data to a JSON file.
+
+        Parameters
+        ----------
+        filename : str | Path
+            Path to the JSON file.
+        """
+        with Path(filename).open("w") as fp:
+            data = convert_from_numpy(self.root, self.transform)
+            json.dump(data, fp, indent=4, sort_keys=True)
+
+    def load_json(self, filename: str | Path, update: bool = False) -> None:
+        """
+        Load data from a JSON file.
+
+        Parameters
+        ----------
+        filename : str | Path
+            Path to the JSON file.
+        update : bool, optional
+            If True, updates the existing data with the loaded data.
+        """
+        with Path(filename).open("r") as fp:
+            data = json.load(fp)
+            data = recursively_load_dict(data, self.inverse_transform)
+            if update:
+                self.root.update(data)
+            else:
+                self.root = data
+
+    def append(self, lock: bool = False) -> None:
+        """
+        Append new keys to the HDF5 file without reading existing data.
+
+        Parameters
+        ----------
+        lock : bool, optional
+            If True, uses a file lock while appending.
+        """
+        if self.filename is not None:
+            lock_file = filelock.FileLock(self.filename + '.lock') if lock else contextlib.nullcontext()
+
+            with lock_file:
+                with h5py.File(self.filename, 'a') as h5file:
+                    recursively_save(h5file, '/', self.root, self.transform)
+        else:
+            print("Filename must be set before save can be used")
+
+    def __str__(self) -> str:
+        """
+        Return a string representation of the object.
+
+        Returns
+        -------
+        str
+            String representation of the filename and root data.
+        """
+        temp = []
+        temp.append(f'Filename = {self.filename}')
+        temp.append(self.pp.pformat(self.root))
+        return '\n'.join(temp)
+
+    def update(self, other: "H5") -> None:
+        """
+        Merge another H5 object's data with the current one.
+
+        Parameters
+        ----------
+        other : H5
+            Another H5 object whose data will be merged.
+        """
+        self.root.update(other.root)
+
+    def __getitem__(self, key: str) -> Any:
+        """
+        Access data by key.
+
+        Parameters
+        ----------
+        key : str
+            Key for accessing nested data.
+
+        Returns
+        -------
+        Any
+            Retrieved data.
+        """
+        key = key.lower()
+        obj = self.root
+        for i in key.split('/'):
+            if i:
+                obj = obj[i]
+        return obj
+
+    def __setitem__(self, key: str, value: Any) -> None:
+        """
+        Set data by key.
+
+        Parameters
+        ----------
+        key : str
+            Key for accessing nested data.
+        value : Any
+            Value to set for the given key.
+        """
+        key = key.lower()
+        obj = self.root
+        parts = key.split('/')
+        for i in parts[:-1]:
+            if i:
+                obj = obj[i]
+        obj[parts[-1]] = value
+
+
+def convert_from_numpy(data: Dict, func: callable) -> Dict:
+    """
+    Convert numpy objects to native Python types.
+
+    Parameters
+    ----------
+    data : Dict
+        Input dictionary potentially containing numpy objects.
+    func : callable
+        Transformation function for dictionary keys.
+
+    Returns
+    -------
+    Dict
+        Dictionary with converted data.
+    """
+    ans = Dict()
+    for key_original, item in data.items():
+        key = func(key_original)
+        if isinstance(item, numpy.ndarray):
+            item = item.tolist()
+
+        if isinstance(item, numpy.generic):
+            item = item.item()
+
+        if isinstance(item, bytes):
+            item = item.decode('ascii')
+
+        if isinstance(item, Dict):
+            ans[key_original] = convert_from_numpy(item, func)
+        else:
+            ans[key] = item
+    return ans
+
+
+def recursively_load_dict(data: dict, func: callable) -> Dict:
+    """
+    Recursively load data from a dictionary.
+
+    Parameters
+    ----------
+    data : dict
+        Input dictionary to load.
+    func : callable
+        Transformation function for dictionary keys.
+
+    Returns
+    -------
+    Dict
+        Dictionary with loaded data.
+    """
+    ans = Dict()
+    for key_original, item in data.items():
+        key = func(key_original)
+        if isinstance(item, dict):
+            ans[key] = recursively_load_dict(item, func)
+        else:
+            ans[key] = item
+    return ans
+
+
+def set_path(obj: Dict, path: str, value: Any) -> None:
+    """
+    Set a value within a nested dictionary given a path.
+
+    Parameters
+    ----------
+    obj : Dict
+        Dictionary to set the value in.
+    path : str
+        Slash-separated path indicating where to set the value.
+    value : Any
+        Value to be set at the specified path.
+    """
+    path_parts = [i for i in path.split('/') if i]
+    temp = obj
+    for part in path_parts[:-1]:
+        temp = temp[part]
+    temp[path_parts[-1]] = value
+
+
+def recursively_load(h5file: h5py.File, path: str, func: callable, paths: Optional[List[str]]) -> Dict:
+    """
+    Recursively load data from an HDF5 file.
+
+    Parameters
+    ----------
+    h5file : h5py.File
+        The HDF5 file to load data from.
+    path : str
+        Path within the HDF5 file.
+    func : callable
+        Transformation function for dictionary keys.
+    paths : Optional[List[str]]
+        Specific paths to load, or None to load everything.
+
+    Returns
+    -------
+    Dict
+        Loaded data.
+    """
+    ans = Dict()
+    if paths is not None:
+        for path in paths:
+            item = h5file.get(path, None)
+            if item is not None:
+                if isinstance(item, h5py._hl.dataset.Dataset):
+                    set_path(ans, path, item[()])
+                elif isinstance(item, h5py._hl.group.Group):
+                    set_path(ans, path, recursively_load(h5file, path + '/', func, None))
+    else:
+        for key_original in h5file[path].keys():
+            key = func(key_original)
+            local_path = path + key
+            item = h5file[path][key_original]
+            if isinstance(item, h5py._hl.dataset.Dataset):
+                ans[key] = item[()]
+            elif isinstance(item, h5py._hl.group.Group):
+                ans[key] = recursively_load(h5file, local_path + '/', func, None)
+    return ans
+
+
+def recursively_save(h5file: h5py.File, path: str, dic: Dict, func: callable) -> None:
+    """
+    Recursively save data to an HDF5 file.
+
+    Parameters
+    ----------
+    h5file : h5py.File
+        The HDF5 file to save data to.
+    path : str
+        Path within the HDF5 file.
+    dic : Dict
+        Dictionary of data to save.
+    func : callable
+        Transformation function for dictionary keys.
+
+    Raises
+    ------
+    ValueError
+        If path or h5file types are invalid, or if the dictionary contains unsupported data types.
+    """
+    if not isinstance(path, str):
+        raise ValueError("path must be a string")
+    if not isinstance(h5file, h5py._hl.files.File):
+        raise ValueError("must be an open h5py file")
+    if not isinstance(dic, dict):
+        raise ValueError("must provide a dictionary")
+
+    for key, item in dic.items():
+        key = str(key)
+        value = None
+
+        if not isinstance(key, str):
+            raise ValueError("dict keys must be strings to save to hdf5")
+
+        if isinstance(item, dict):
+            recursively_save(h5file, path + key + '/', item, func)
+        elif isinstance(item, str):
+            value = numpy.array(item.encode('ascii'))
+        elif isinstance(item, list) and all(isinstance(i, str) for i in item):
+            value = numpy.array([i.encode('ascii') for i in item])
+        else:
+            try:
+                value = numpy.array(item)
+            except TypeError:
+                raise ValueError(f'Cannot save {path}/{func(key)} key with {type(item)} type.')
+
+        if value is not None:
+            try:
+                h5file[path + func(key)] = value
+            except OSError as e:
+                if str(e) == 'Unable to create link (name already exists)':
+                    raise KeyError(f'Name conflict with upper and lower case entries for key "{path}{key}".')
+                else:
+                    raise
diff --git a/cadet/runner.py b/cadet/runner.py
new file mode 100644
index 0000000..fff6304
--- /dev/null
+++ b/cadet/runner.py
@@ -0,0 +1,257 @@
+import os
+import pathlib
+import re
+import subprocess
+from abc import ABC, abstractmethod
+from dataclasses import dataclass
+from pathlib import Path
+from typing import Optional
+
+
+@dataclass
+class ReturnInformation:
+    """
+    Class to store information about a CADET run return status.
+
+    Parameters
+    ----------
+    return_code : int
+        An integer representing the return code. 0 indicates success, non-zero values indicate errors.
+    error_message : str
+        A string containing the error message if an error occurred. Empty if no error.
+    log : str
+        A string containing log information.
+    """
+    return_code: int
+    error_message: str
+    log: str
+
+
+class CadetRunnerBase(ABC):
+    """
+    Abstract base class for CADET runners.
+
+    Subclasses must implement the `run`, `clear`, and `load_results` methods.
+    """
+
+    @abstractmethod
+    def run(
+            self,
+            simulation: "Cadet",
+            timeout: Optional[int] = None,
+    ) -> ReturnInformation:
+        """
+        Run a CADET simulation.
+
+        Parameters
+        ----------
+        simulation : Cadet
+            The simulation object.
+        timeout : Optional[int]
+            Maximum time allowed for the simulation to run, in seconds.
+
+        Returns
+        -------
+        ReturnInformation
+            Information about the simulation run.
+        """
+        pass
+
+    @abstractmethod
+    def clear(self) -> None:
+        """
+        Clear the simulation data.
+        """
+        pass
+
+    @abstractmethod
+    def load_results(self, sim: "Cadet") -> None:
+        """
+        Load the results of the simulation into the provided object.
+
+        Parameters
+        ----------
+        sim : Cadet
+            The simulation object where results will be loaded.
+        """
+        pass
+
+    @property
+    @abstractmethod
+    def cadet_version(self) -> str:
+        pass
+
+    @property
+    @abstractmethod
+    def cadet_branch(self) -> str:
+        pass
+
+    @property
+    @abstractmethod
+    def cadet_build_type(self) -> str:
+        pass
+
+    @property
+    @abstractmethod
+    def cadet_commit_hash(self) -> str:
+        pass
+
+    @property
+    @abstractmethod
+    def cadet_path(self) -> Optional[os.PathLike]:
+        pass
+
+
+class CadetCLIRunner(CadetRunnerBase):
+    """
+    File-based CADET runner.
+
+    This class runs CADET simulations using a command-line interface (CLI) executable.
+    """
+
+    def __init__(self, cadet_path: str | os.PathLike) -> None:
+        """
+        Initialize the CadetFileRunner.
+
+        Parameters
+        ----------
+        cadet_path : os.PathLike
+            Path to the CADET CLI executable.
+        """
+        cadet_path = Path(cadet_path)
+
+        self._cadet_path = cadet_path
+        self._get_cadet_version()
+
+    def run(
+            self,
+            simulation: "Cadet",
+            timeout: Optional[int] = None,
+    ) -> ReturnInformation:
+        """
+        Run a CADET simulation using the CLI executable.
+
+        Parameters
+        ----------
+        simulation : Cadet
+            Not used in this runner.
+        timeout : Optional[int]
+            Maximum time allowed for the simulation to run, in seconds.
+
+        Raises
+        ------
+        RuntimeError
+            If the simulation process returns a non-zero exit code.
+
+        Returns
+        -------
+        ReturnInformation
+            Information about the simulation run.
+        """
+        if simulation.filename is None:
+            raise ValueError("Filename must be set before run can be used")
+
+        data = subprocess.run(
+            [self.cadet_path, str(simulation.filename)],
+            timeout=timeout,
+            capture_output=True
+        )
+
+        return_info = ReturnInformation(
+            return_code=data.returncode,
+            error_message=data.stderr.decode('utf-8'),
+            log=data.stdout.decode('utf-8')
+        )
+
+        return return_info
+
+    def clear(self) -> None:
+        """
+        Clear the simulation data.
+
+        This method can be extended if any cleanup is required.
+        """
+        pass
+
+    def load_results(self, sim: "Cadet") -> None:
+        """
+        Load the results of the simulation into the provided object.
+
+        Parameters
+        ----------
+        sim : Cadet
+            The simulation object where results will be loaded.
+        """
+        sim.load(paths=["/meta", "/output"], update=True)
+
+    def _get_cadet_version(self) -> dict:
+        """
+        Get version and branch name of the currently instanced CADET build.
+        Returns
+        -------
+        dict
+            Dictionary containing: cadet_version as x.x.x, cadet_branch, cadet_build_type, cadet_commit_hash
+        Raises
+        ------
+        ValueError
+            If version and branch name cannot be found in the output string.
+        RuntimeError
+            If any unhandled event during running the subprocess occurs.
+        """
+        try:
+            result = subprocess.run(
+                [self.cadet_path, '--version'],
+                check=True,
+                stdout=subprocess.PIPE,
+                stderr=subprocess.PIPE,
+                text=True
+            )
+            version_output = result.stdout.strip()
+
+            version_match = re.search(
+                r'cadet-cli version ([\d.]+) \((.*) branch\)\n',
+                version_output
+            )
+
+            commit_hash_match = re.search(
+                "Built from commit (.*)\n",
+                version_output
+            )
+
+            build_variant_match = re.search(
+                "Build variant (.*)\n",
+                version_output
+            )
+
+            if version_match:
+                self._cadet_version = version_match.group(1)
+                self._cadet_branch = version_match.group(2)
+                self._cadet_commit_hash = commit_hash_match.group(1)
+                if build_variant_match:
+                    self._cadet_build_type = build_variant_match.group(1)
+                else:
+                    self._cadet_build_type = None
+            else:
+                raise ValueError("CADET version or branch name missing from output.")
+        except subprocess.CalledProcessError as e:
+            raise RuntimeError(f"Command execution failed: {e}")
+
+    @property
+    def cadet_version(self) -> str:
+        return self._cadet_version
+
+    @property
+    def cadet_branch(self) -> str:
+        return self._cadet_branch
+
+    @property
+    def cadet_build_type(self) -> str:
+        return self._cadet_build_type
+
+    @property
+    def cadet_commit_hash(self) -> str:
+        return self._cadet_commit_hash
+
+    @property
+    def cadet_path(self) -> os.PathLike:
+        return self._cadet_path
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..1a1e629
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,57 @@
+[build-system]
+requires = [
+    "setuptools>=69",
+    "wheel",
+]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "CADET-Python"
+dynamic = ["version"]
+authors = [
+    { name = "William Heymann", email = "w.heymann@fz-juelich.de" },
+    { name = "Samuel Leweke", email = "s.leweke@fz-juelich.de" },
+    { name = "Johannes Schmölder", email = "j.schmoelder@fz-juelich.de" },
+    { name = "Ronald Jäpel", email = "r.jaepel@fz-juelich.de" },
+]
+description = "CADET-Python is a Python interface to the CADET-Core simulator"
+readme = "README.md"
+requires-python = ">=3.10"
+keywords = ["process modeling", "process optimization", "chromatography"]
+license = { text = "BSD-3-Clause" }
+classifiers = [
+    "Programming Language :: Python :: 3",
+    "Operating System :: OS Independent",
+    "License :: OSI Approved :: BSD License",
+    "Intended Audience :: Science/Research",
+]
+dependencies = [
+    "addict",
+    "numpy",
+    "h5py",
+    "filelock",
+]
+
+[project.optional-dependencies]
+testing = ["pytest"]
+
+[project.urls]
+"homepage" = "https://github.com/cadet/CADET-Python"
+"Bug Tracker" = "https://github.com/cadet/CADET-Python/issues"
+
+[tool.setuptools.dynamic]
+version = { attr = "cadet.__version__" }
+
+[tool.setuptools.packages.find]
+include = ["cadet*"]
+
+[tool.ruff]
+src = ["cadet"]
+line-length = 88
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]
+markers = [
+    "slow: marks tests as slow (deselect with '-m \"not slow\"')",
+    "local: marks tests as only useful on local installs (deselect with '-m \"not local\"')",
+]
\ No newline at end of file
diff --git a/setup.cfg b/setup.cfg
deleted file mode 100644
index 23d72ac..0000000
--- a/setup.cfg
+++ /dev/null
@@ -1,7 +0,0 @@
-# Inside of setup.cfg
-[metadata]
-description-file = README.md
-
-[tool:pytest]
-markers =
-    slow: marks tests as slow (deselect with '-m "not slow"')
diff --git a/setup.py b/setup.py
deleted file mode 100644
index 54cf7e2..0000000
--- a/setup.py
+++ /dev/null
@@ -1,35 +0,0 @@
-import setuptools
-
-with open("README.md", "r") as fh:
-    long_description = fh.read()
-    
-setuptools.setup(
-    name="CADET-Python",
-    version="0.14.1",
-    author="William Heymann",
-    author_email="w.heymann@fz-juelich.de",
-    description="CADET is a python interface to the CADET chromatography simulator",
-    long_description=long_description,
-    long_description_content_type="text/markdown",
-    url="https://github.com/modsim/CADET-Python",
-    packages=setuptools.find_packages(),
-    install_requires=[
-          'addict',
-          'numpy',
-          'h5py',
-          'filelock'
-      ],
-    classifiers=[
-        "Programming Language :: Python :: 3",
-        "License :: OSI Approved :: BSD License",
-        "Operating System :: OS Independent",
-    ],
-    extras_require={
-        "testing": [
-            "pytest",
-            "matplotlib",
-            "pandas",
-        ],
-    },
-    python_requires='>=3.7',
-) 


Test patch:
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_MSSMA_2comp.py b/tests/test_MSSMA_2comp.py
deleted file mode 100644
index 06f6713..0000000
--- a/tests/test_MSSMA_2comp.py
+++ /dev/null
@@ -1,201 +0,0 @@
-import os
-
-import matplotlib.pyplot as plt
-import numpy
-import pandas
-import pytest
-
-from cadet import Cadet
-from tests import common
-
-
-def gen_fraction_times(start, stop, bins):
-    return numpy.linspace(start, stop, bins + 1, endpoint=True)
-
-
-def gen_fractions(fraction_times, sim):
-    nComp = sim.root.input.model.unit_000.ncomp
-
-    df = pandas.DataFrame()
-
-    times = sim.root.output.solution.solution_times
-
-    for idx, (start, stop) in enumerate(zip(fraction_times[:-1], fraction_times[1:])):
-        selected = (times >= start) & (times <= stop)
-        temp = {'Start': start, 'Stop': stop}
-        for comp in range(nComp):
-            local_times = times[selected]
-            local_values = sim.root.output.solution.unit_001["solution_outlet_comp_%03d" % comp][selected]
-
-            temp[str(comp)] = numpy.trapz(local_values, local_times) / (stop - start)
-        df = df.append(temp, ignore_index=True)
-    return df
-
-
-def main():
-    if not os.path.exists("tmp"):
-        os.makedirs("tmp")
-    simulation = Cadet(common.common.root)
-    simulation.filename = "tmp/MSSMA_2comp.h5"
-    createSimulation(simulation)
-    simulation.save()
-    data = simulation.run()
-    print(data)
-    simulation.load()
-    plotSimulation(simulation)
-    return simulation
-
-
-def createSimulation(simulation):
-    root = simulation.root
-
-    root.input.model.nunits = 3
-
-    root.input.model.connections.nswitches = 1
-    root.input.model.connections.switch_000.section = 0
-    root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.0,
-                                                           1, 2, -1, -1, 1.0]
-    root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
-    root.input.model.unit_000.unit_type = 'INLET'
-    root.input.model.unit_000.ncomp = 3
-
-    root.input.model.unit_000.sec_000.const_coeff = [100.0, 0.15, 0.01]
-    root.input.model.unit_000.sec_000.lin_coeff = [0.0, 0.0, 0.0]
-    root.input.model.unit_000.sec_000.quad_coeff = [0.0, 0.0, 0.0]
-    root.input.model.unit_000.sec_000.cube_coeff = [0.0, 0.0, 0.0]
-
-    root.input.model.unit_000.sec_001.const_coeff = [75, 0.0, 0.0]
-    root.input.model.unit_000.sec_001.lin_coeff = [0.0, 0.0, 0.0]
-    root.input.model.unit_000.sec_001.quad_coeff = [0.0, 0.0, 0.0]
-    root.input.model.unit_000.sec_001.cube_coeff = [0.0, 0.0, 0.0]
-
-    root.input.model.unit_000.sec_002.const_coeff = [75, 0.0, 0.0]
-    root.input.model.unit_000.sec_002.lin_coeff = [0.05, 0.0, 0.0]
-    root.input.model.unit_000.sec_002.quad_coeff = [0.0, 0.0, 0.0]
-    root.input.model.unit_000.sec_002.cube_coeff = [0.0, 0.0, 0.0]
-
-    root.input.model.unit_001.adsorption_model = 'MULTISTATE_STERIC_MASS_ACTION'
-    root.input.model.unit_001.col_dispersion = 1.5E-7
-    root.input.model.unit_001.col_length = 0.25
-    root.input.model.unit_001.col_porosity = 0.3
-    root.input.model.unit_001.film_diffusion = [2.14E-4, 2.1e-5, 2.1e-5]
-    root.input.model.unit_001.init_c = [75, 0.0, 0.0]
-    root.input.model.unit_001.init_q = [225, 0.0, 0.0, 0.0, 0.0]
-    root.input.model.unit_001.ncomp = 3
-    root.input.model.unit_001.par_diffusion = [4.08E-10, 9.0E-12, 9.0e-12]
-    root.input.model.unit_001.par_porosity = 0.4
-    root.input.model.unit_001.par_radius = 3.25E-5
-    root.input.model.unit_001.par_surfdiffusion = [0.0, 0.0, 0.0, 0.0, 0.0]
-    root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL'
-    root.input.model.unit_001.velocity = 0.001
-
-    root.input.model.unit_001.adsorption.is_kinetic = 1
-    root.input.model.unit_001.adsorption.mssma_ka = [0.0, 1E11, 8E6, 1E11, 8E6]
-    root.input.model.unit_001.adsorption.mssma_kd = [0.0, 6E11, 2E16, 6E11, 2E16]
-    root.input.model.unit_001.adsorption.mssma_lambda = 225
-    root.input.model.unit_001.adsorption.mssma_nu = [0.0, 10, 25, 20, 50]
-    root.input.model.unit_001.adsorption.mssma_sigma = [0.0, 48, 66, 48 * 2, 66 * 2]
-    root.input.model.unit_001.adsorption.mssma_refq = 225
-    root.input.model.unit_001.adsorption.mssma_refc0 = 520.0
-    root.input.model.unit_001.adsorption.mssma_rates = [0.0, 0.0, 1e20, 10, 0.0, 0.0, 1e20, 10, 0.0]
-
-    root.input.model.unit_001.discretization.nbound = [1, 2, 2]
-    root.input.model.unit_001.discretization.ncol = 50
-    root.input.model.unit_001.discretization.npar = 5
-
-    root.input.model.unit_002.ncomp = 3
-    root.input.model.unit_002.unit_type = 'OUTLET'
-
-    root.input.solver.user_solution_times = numpy.linspace(0, 15000, 1000)
-    root.input.solver.sections.nsec = 3
-    root.input.solver.sections.section_continuity = [0, 0]
-    root.input.solver.sections.section_times = [0.0, 4500, 6100, 15000]
-
-    root.input.solver.consistent_init_mode = 3
-
-
-def plotSimulation(simulation):
-    f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8])
-    plotInlet(ax1, simulation)
-    plotOutlet(ax2, simulation)
-    f.tight_layout()
-    plt.savefig("tmp/MSSMA.png")
-
-
-def plotInlet(axis, simulation):
-    solution_times = simulation.root.output.solution.solution_times
-
-    inlet_salt = simulation.root.output.solution.unit_000.solution_inlet_comp_000
-    inlet_p1 = simulation.root.output.solution.unit_000.solution_inlet_comp_001
-    inlet_p2 = simulation.root.output.solution.unit_000.solution_inlet_comp_002
-
-    axis.set_title("Inlet")
-    axis.plot(solution_times, inlet_salt, 'b-', label="Salt")
-    axis.set_xlabel('time (s)')
-
-    # Make the y-axis label, ticks and tick labels match the line color.
-    axis.set_ylabel('mMol Salt', color='b')
-    axis.tick_params('y', colors='b')
-
-    axis2 = axis.twinx()
-    axis2.plot(solution_times, inlet_p1, 'r-', label="P1")
-    axis2.plot(solution_times, inlet_p2, 'g-', label="P2")
-    axis2.set_ylabel('mMol Protein', color='r')
-    axis2.tick_params('y', colors='r')
-
-    lines, labels = axis.get_legend_handles_labels()
-    lines2, labels2 = axis2.get_legend_handles_labels()
-    axis2.legend(lines + lines2, labels + labels2, loc=0)
-
-
-def plotOutlet(axis, simulation):
-    solution_times = simulation.root.output.solution.solution_times
-
-    outlet_salt = simulation.root.output.solution.unit_002.solution_outlet_comp_000
-    outlet_p1 = simulation.root.output.solution.unit_002.solution_outlet_comp_001
-    outlet_p2 = simulation.root.output.solution.unit_002.solution_outlet_comp_002
-
-    data = numpy.vstack([solution_times, outlet_p1]).transpose()
-    numpy.savetxt('tmp/comp2_1.csv', data, delimiter=',')
-
-    data = numpy.vstack([solution_times, outlet_p2]).transpose()
-    numpy.savetxt('tmp/comp2_2.csv', data, delimiter=',')
-
-    data = numpy.vstack([solution_times, outlet_p1 + outlet_p2]).transpose()
-    numpy.savetxt('tmp/comp2_comb.csv', data, delimiter=',')
-
-    fraction_times = gen_fraction_times(6000, 14000, 10)
-    df = gen_fractions(fraction_times, simulation)
-    df.to_csv("tmp/comp2_fraction.csv", columns=('Start', 'Stop', '1'), index=False)
-
-    axis.set_title("Output")
-    axis.plot(solution_times, outlet_salt, 'b-', label="Salt")
-    axis.set_xlabel('time (s)')
-
-    # Make the y-axis label, ticks and tick labels match the line color.
-    axis.set_ylabel('mMol Salt', color='b')
-    axis.tick_params('y', colors='b')
-
-    axis2 = axis.twinx()
-    axis2.plot(solution_times, outlet_p1, 'r-', label="P1")
-    axis2.plot(solution_times, outlet_p2, 'g-', label="P2")
-    axis2.set_ylabel('mMol Protein', color='r')
-    axis2.tick_params('y', colors='r')
-
-    lines, labels = axis.get_legend_handles_labels()
-    lines2, labels2 = axis2.get_legend_handles_labels()
-    axis2.legend(lines + lines2, labels + labels2, loc=0)
-
-
-@pytest.mark.slow
-def test_MSSMA_2comp():
-    sim = main()
-    assert isinstance(sim.root.output.solution.unit_002.solution_outlet_comp_001, numpy.ndarray)
-    assert isinstance(sim.root.output.solution.unit_002.solution_outlet_comp_000, numpy.ndarray)
-
-
-if __name__ == "__main__":
-    import sys
-
-    print(sys.version)
-    main()
diff --git a/tests/test_SMB.py b/tests/test_SMB.py
deleted file mode 100644
index 07a7e22..0000000
--- a/tests/test_SMB.py
+++ /dev/null
@@ -1,258 +0,0 @@
-# Basic 4-zone SMB setup
-
-import math
-import os.path
-
-import matplotlib.pyplot as plt
-import numpy
-import pytest
-
-from cadet import Cadet
-
-# number of columns in a cycle
-cycle_size = 8
-
-# number of cycles
-cycles = 4
-
-# number of times flows have to be expanded for a 4-zone model
-repeat_size = int(cycle_size / 4)
-
-
-def gen_connections(units, cycle_size, size, step, flows, flows_static):
-    temp = []
-    connections = list(zip(units, numpy.roll(units, -1), flows))
-    io = numpy.roll(units, step)[[0, size * 2, size - 1, size * 3 - 1]]
-    ios = list(zip([0, 1, 2, 3], io))
-
-    for connection in connections:
-        temp.append([connection[0], connection[1], -1, -1, connection[2]])
-    # inputs
-    idx = 0
-    for io in ios[:2]:
-        temp.append([io[0], io[1], -1, -1, flows_static[idx]])
-        idx += 1
-    # outputs
-    for io in ios[2:]:
-        temp.append([io[1], io[0], -1, -1, flows_static[idx]])
-        idx += 1
-    return numpy.array(temp)
-
-
-def expand_flow(seq, inlet1, inlet2, number):
-    "expand the flows for smb, this is more complex since we need link values"
-    temp = []
-    temp.extend([seq[3] + inlet1] * (number - 1))
-    temp.append(seq[0])
-
-    temp.extend([seq[0]] * (number - 1))
-    temp.append(seq[1])
-
-    temp.extend([seq[1] + inlet2] * (number - 1))
-    temp.append(seq[2])
-
-    temp.extend([seq[2]] * (number - 1))
-    temp.append(seq[3])
-
-    return temp
-
-
-def main():
-    if not os.path.exists("tmp"):
-        os.makedirs("tmp")
-    smb = Cadet()
-    smb.filename = 'tmp/SMB.h5'
-    createSimulation(smb)
-    print("Simulated Created")
-    smb.save()
-    smb.run()
-    smb.load()
-    print("Simulation Run")
-    plotSimulation(smb)
-    return smb
-
-
-def createSimulation(simulation):
-    simulation.root.input.model.nunits = 4 + cycle_size
-
-    simulation.root.input.model.solver.gs_type = 1
-    simulation.root.input.model.solver.max_krylov = 0
-    simulation.root.input.model.solver.max_restarts = 0
-    simulation.root.input.model.solver.schur_safety = 1e-8
-
-    # setup connections
-    simulation.root.input.model.connections.nswitches = cycle_size
-
-    units = range(4, 4 + cycle_size)
-
-    flows = expand_flow([7.66E-07, 7.66E-07, 8.08E-07, 8.08E-07], 0.98e-7, 1.96e-07, repeat_size)
-    flows_static = numpy.array([0.98e-7, 1.96e-7, 1.4e-7, 1.54e-7])
-
-    for i in range(cycle_size):
-        simulation.root.input.model.connections["switch_%03d" % i].section = i
-        simulation.root.input.model.connections["switch_%03d" % i].connections = gen_connections(units, cycle_size,
-                                                                                                 repeat_size, -i,
-                                                                                                 numpy.array(list(
-                                                                                                     numpy.roll(flows,
-                                                                                                                i))),
-                                                                                                 flows_static)
-
-    # setup inlets
-    simulation.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
-    simulation.root.input.model.unit_000.ncomp = 2
-    simulation.root.input.model.unit_000.unit_type = 'INLET'
-
-    for i in range(cycle_size):
-        # section
-        simulation.root.input.model.unit_000["sec_%03d" % i].const_coeff = [0.55 / 180.16, 0.55 / 180.16]
-        simulation.root.input.model.unit_000["sec_%03d" % i].lin_coeff = [0.0, 0.0]
-        simulation.root.input.model.unit_000["sec_%03d" % i].quad_coeff = [0.0, 0.0]
-        simulation.root.input.model.unit_000["sec_%03d" % i].cube_coeff = [0.0, 0.0]
-
-    simulation.root.input.model.unit_001.inlet_type = 'PIECEWISE_CUBIC_POLY'
-    simulation.root.input.model.unit_001.ncomp = 2
-    simulation.root.input.model.unit_001.unit_type = 'INLET'
-
-    for i in range(cycle_size):
-        # section
-        simulation.root.input.model.unit_001["sec_%03d" % i].const_coeff = [0.0, 0.0]
-        simulation.root.input.model.unit_001["sec_%03d" % i].lin_coeff = [0.0, 0.0]
-        simulation.root.input.model.unit_001["sec_%03d" % i].quad_coeff = [0.0, 0.0]
-        simulation.root.input.model.unit_001["sec_%03d" % i].cube_coeff = [0.0, 0.0]
-
-    # create columns
-    for unit in range(4, 4 + cycle_size):
-        simulation.root.input.model["unit_%03d" % unit].unit_type = 'GENERAL_RATE_MODEL'
-
-        col = simulation.root.input.model["unit_%03d" % unit]
-
-        col.ncomp = 2
-        col.cross_section_area = math.pi * (0.02 ** 2) / 4.0
-        col.col_dispersion = 3.8148e-20
-        col.col_length = 0.25 / repeat_size
-        col.col_porosity = 0.83
-        col.init_c = [0.0, 0.0]
-        col.init_q = [0.0, 0.0]
-        col.film_diffusion = [100.0, 100.0]
-        col.par_diffusion = [1.6e4, 1.6e4]
-        col.par_radius = 0.0005
-        col.par_porosity = 0.000001
-        col.par_surfdiffusion = [0.0, 0.0]
-        col.adsorption_model = 'LINEAR'
-        col.velocity = 1.0
-
-        col.adsorption.is_kinetic = 0
-        col.adsorption.lin_ka = [5.72, 7.7]
-        col.adsorption.lin_kd = [1.0, 1.0]
-
-        col.discretization.gs_type = 1
-        col.discretization.max_krylov = 0
-        col.discretization.max_restarts = 0
-        col.discretization.nbound = [1, 1]
-        col.discretization.ncol = 40
-        col.discretization.npar = 1
-        col.discretization.par_disc_type = 'EQUIDISTANT_PAR'
-        col.discretization.par_disc_vector = [0.0, 0.5, 1.0]
-        col.discretization.reconstruction = 'WENO'
-        col.discretization.schur_safety = 1e-8
-        col.discretization.use_analytic_jacobian = 1
-
-        col.discretization.weno.boundary_model = 0
-        col.discretization.weno.weno_eps = 1e-12
-        col.discretization.weno.weno_order = 3
-
-    # create outlets
-    simulation.root.input.model.unit_002.ncomp = 2
-    simulation.root.input.model.unit_002.unit_type = 'OUTLET'
-
-    simulation.root.input.model.unit_003.ncomp = 2
-    simulation.root.input.model.unit_003.unit_type = 'OUTLET'
-
-    # create output information
-
-    simulation.root.input['return'].write_solution_times = 1
-
-    ret = simulation.root.input['return']
-
-    ret.write_solution_last = 1
-
-    ret.unit_002.write_sens_column = 0
-    ret.unit_002.write_sens_column_inlet = 0
-    ret.unit_002.write_sens_column_outlet = 0
-    ret.unit_002.write_sens_flux = 0
-    ret.unit_002.write_sens_particle = 0
-
-    ret.unit_002.write_solution_column = 0
-    ret.unit_002.write_solution_column_inlet = 1
-    ret.unit_002.write_solution_column_outlet = 1
-    ret.unit_002.write_solution_flux = 0
-    ret.unit_002.write_solution_particle = 0
-
-    ret.unit_003.write_sens_column = 0
-    ret.unit_003.write_sens_column_inlet = 0
-    ret.unit_003.write_sens_column_outlet = 0
-    ret.unit_003.write_sens_flux = 0
-    ret.unit_003.write_sens_particle = 0
-
-    ret.unit_003.write_solution_column = 0
-    ret.unit_003.write_solution_column_inlet = 1
-    ret.unit_003.write_solution_column_outlet = 1
-    ret.unit_003.write_solution_flux = 0
-    ret.unit_003.write_solution_particle = 0
-
-    simulation.root.input.solver.nthreads = 0
-    simulation.root.input.solver.user_solution_times = numpy.linspace(0, cycles * 180 * 4, 1000 * cycle_size * cycles)
-    simulation.root.input.solver.sections.nsec = cycle_size * cycles
-    simulation.root.input.solver.sections.section_continuity = [0] * (cycle_size * cycles - 1)
-    simulation.root.input.solver.sections.section_times = [float(i) * 180 * 4.0 / cycle_size for i in
-                                                           range(cycle_size * cycles + 1)]
-
-    simulation.root.input.solver.time_integrator.abstol = 1e-10
-    simulation.root.input.solver.time_integrator.algtol = 1e-10
-    simulation.root.input.solver.time_integrator.init_step_size = 1e-14
-    simulation.root.input.solver.time_integrator.max_steps = 1e6
-    simulation.root.input.solver.time_integrator.reltol = 1e-6
-
-
-def plotSimulation(simulation):
-    solution_times = simulation.root.output.solution.solution_times
-
-    e_0 = simulation.root.output.solution.unit_002.solution_outlet_comp_000
-    e_1 = simulation.root.output.solution.unit_002.solution_outlet_comp_001
-
-    r_0 = simulation.root.output.solution.unit_003.solution_outlet_comp_000
-    r_1 = simulation.root.output.solution.unit_003.solution_outlet_comp_001
-
-    fig = plt.figure(figsize=[10, 2 * 10])
-
-    graph = fig.add_subplot(2, 1, 1)
-    graph.set_title("Extract")
-
-    graph.plot(solution_times, e_0, 'r', label='1')
-    graph.plot(solution_times, e_1, 'g', label='2')
-    graph.legend()
-
-    graph = fig.add_subplot(2, 1, 2)
-    graph.set_title("Raffinate")
-
-    graph.plot(solution_times, r_0, 'r', label='1')
-    graph.plot(solution_times, r_1, 'g', label='2')
-    graph.legend()
-    plt.savefig("tmp/SMB.png")
-
-
-@pytest.mark.slow
-def test_SMB():
-    from datetime import datetime
-    start = datetime.now()
-    sim = main()
-    print(datetime.now() - start)
-    assert isinstance(sim.root.output.solution.unit_003.solution_outlet_comp_001, numpy.ndarray)
-    assert isinstance(sim.root.output.solution.unit_003.solution_outlet_comp_000, numpy.ndarray)
-
-
-if __name__ == "__main__":
-    import sys
-
-    print(sys.version)
-    main()
diff --git a/tests/test_dll.py b/tests/test_dll.py
new file mode 100644
index 0000000..fddc2bf
--- /dev/null
+++ b/tests/test_dll.py
@@ -0,0 +1,938 @@
+from pathlib import Path
+import subprocess
+import pytest
+
+from cadet import Cadet
+
+
+# %% Utility methods
+
+# Use this to specify custom cadet_roots if you require it.
+cadet_root = None
+
+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.
+    """
+
+    cadet_model = Cadet(install_path=cadet_root, use_dll=True)
+
+    args = [
+        cadet_model.cadet_create_lwe_path,
+        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"
+        )
+
+    cadet_model.filename = file_name
+    cadet_model.load()
+
+    return cadet_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,),
+            'soldot_inlet': (1501, 4),
+            'soldot_outlet': (1501, 4),
+            'soldot_bulk': (1501, 4),
+            'soldot_solid': (1501, 4),
+            'soldot_volume': (1501, ),
+            '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
+
+
+@pytest.mark.parametrize("use_dll", use_dll)
+@pytest.mark.parametrize("test_case", [grm])
+def test_meta(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
+    )
+
+    meta_information = {
+        'cadet_branch': str,
+        'cadet_commit': str,
+        'cadet_version': str,
+        'file_format': int,
+        'time_sim': float,
+    }
+
+    assert model.root.meta.keys() == meta_information.keys()
+
+    for meta_key, meta_type in meta_information.items():
+        assert isinstance(model.root.meta[meta_key], meta_type)
+
+
+if __name__ == "__main__":
+    pytest.main(["test_dll.py"])
diff --git a/tests/test_install_path_settings.py b/tests/test_install_path_settings.py
new file mode 100644
index 0000000..99db92a
--- /dev/null
+++ b/tests/test_install_path_settings.py
@@ -0,0 +1,69 @@
+import os
+from pathlib import Path
+
+import pytest
+import re
+
+from cadet import Cadet
+
+""" These tests require two distinct CADET installations to compare between and should not run in the CI"""
+
+
+# Full path to cadet.dll or cadet.so, that is different from the system/conda cadet
+full_path_dll = Path("path/to/cadet")
+
+install_path_conda = Cadet.autodetect_cadet()
+
+
+def test_autodetection():
+    sim = Cadet()
+    assert sim.install_path == install_path_conda
+    assert sim.cadet_dll_path.parent.parent == install_path_conda
+    assert sim.cadet_cli_path.parent.parent == install_path_conda
+    assert sim.cadet_runner.cadet_path.suffix not in [".dll", ".so"]
+
+@pytest.mark.local
+def test_install_path():
+    sim = Cadet(install_path=full_path_dll, use_dll=True)
+    assert sim.cadet_dll_path == full_path_dll
+    assert sim.cadet_runner.cadet_path.suffix in [".dll", ".so"]
+
+    sim = Cadet()
+    sim.install_path = full_path_dll.parent.parent
+    sim.use_dll = True
+    assert sim.cadet_dll_path == full_path_dll
+    assert sim.cadet_runner.cadet_path.suffix in [".dll", ".so"]
+
+    sim = Cadet()
+    with pytest.deprecated_call():
+        sim.cadet_path = full_path_dll.parent.parent
+
+    sim.use_dll = True
+    assert sim.cadet_dll_path == full_path_dll
+    assert sim.cadet_runner.cadet_path.suffix in [".dll", ".so"]
+
+
+@pytest.mark.local
+def test_dll_runner_attrs():
+    cadet = Cadet(full_path_dll.parent.parent)
+    cadet_runner = cadet._cadet_dll_runner
+    assert re.match("\d\.\d\.\d", cadet_runner.cadet_version)
+    assert isinstance(cadet_runner.cadet_branch, str)
+    assert isinstance(cadet_runner.cadet_build_type, str | None)
+    assert isinstance(cadet_runner.cadet_commit_hash, str)
+    assert isinstance(cadet_runner.cadet_path, str | os.PathLike)
+
+
+@pytest.mark.local
+def test_cli_runner_attrs():
+    cadet = Cadet(full_path_dll.parent.parent)
+    cadet_runner = cadet._cadet_cli_runner
+    assert re.match("\d\.\d\.\d", cadet_runner.cadet_version)
+    assert isinstance(cadet_runner.cadet_branch, str)
+    assert isinstance(cadet_runner.cadet_build_type, str | None)
+    assert isinstance(cadet_runner.cadet_commit_hash, str)
+    assert isinstance(cadet_runner.cadet_path, str | os.PathLike)
+
+
+if __name__ == '__main__':
+    pytest.main(["test_install_path_settings.py"])
diff --git a/tests/test_lwe.py b/tests/test_lwe.py
deleted file mode 100644
index bd0acad..0000000
--- a/tests/test_lwe.py
+++ /dev/null
@@ -1,165 +0,0 @@
-# Normally the main function is placed at the bottom of the file, but I have placed it at the top so that
-# This interface is more like a tutorial and can be read from the top down and any given function
-# can be drilled into to learn more about it.
-import os
-
-import matplotlib.pyplot as plt
-import numpy
-
-from cadet import Cadet
-from tests import common
-
-
-# Helper functions that make it easier to set the values in the HDF5 file
-# In the CADET pdf file each value in the hdf5 file has a type. The functions
-# below match those types.
-
-def main():
-    if not os.path.exists("tmp"):
-        os.makedirs("tmp")
-    simulation = Cadet(common.common.root)
-    simulation.filename = "tmp/LWE.h5"
-    createSimulation(simulation)
-    simulation.save()
-    simulation.run()
-    simulation.load()
-
-    plotSimulation(simulation)
-    return simulation
-
-
-def createSimulation(simulation):
-    root = simulation.root
-
-    root.input.model.nunits = 3
-    root.input.model.unit_001.discretization.use_analytic_jacobian = 1
-
-    root.input.model.connections.nswitches = 1
-    root.input.model.connections.switch_000.section = 0
-    root.input.model.connections.switch_000.connections = [0, 1, -1, -1, 1.0,
-                                                           1, 2, -1, -1, 1.0]
-    root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
-    root.input.model.unit_000.unit_type = 'INLET'
-    root.input.model.unit_000.ncomp = 2
-
-    root.input.model.unit_000.sec_000.const_coeff = [50.0, 1.0]
-    root.input.model.unit_000.sec_000.lin_coeff = [0.0, 0.0]
-    root.input.model.unit_000.sec_000.quad_coeff = [0.0, 0.0]
-    root.input.model.unit_000.sec_000.cube_coeff = [0.0, 0.0]
-
-    root.input.model.unit_000.sec_001.const_coeff = [50.0, 0.0]
-    root.input.model.unit_000.sec_001.lin_coeff = [0.0, 0.0]
-    root.input.model.unit_000.sec_001.quad_coeff = [0.0, 0.0]
-    root.input.model.unit_000.sec_001.cube_coeff = [0.0, 0.0]
-
-    root.input.model.unit_000.sec_002.const_coeff = [100.0, 0.0]
-    root.input.model.unit_000.sec_002.lin_coeff = [0.2, 0.0]
-    root.input.model.unit_000.sec_002.quad_coeff = [0.0, 0.0]
-    root.input.model.unit_000.sec_002.cube_coeff = [0.0, 0.0]
-
-    root.input.model.unit_001.adsorption_model = 'STERIC_MASS_ACTION'
-    root.input.model.unit_001.col_dispersion = 5.75e-8
-    root.input.model.unit_001.col_length = 0.014
-    root.input.model.unit_001.col_porosity = 0.37
-    root.input.model.unit_001.film_diffusion = [6.9e-6, 6.9e-6]
-    root.input.model.unit_001.init_c = [50.0, 0.0]
-    root.input.model.unit_001.init_q = [1200.0, 0.0]
-    root.input.model.unit_001.ncomp = 2
-    root.input.model.unit_001.par_diffusion = [7e-10, 6.07e-11]
-    root.input.model.unit_001.par_porosity = 0.75
-    root.input.model.unit_001.par_radius = 4.5e-5
-    root.input.model.unit_001.par_surfdiffusion = [0.0, 0.0]
-    root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL'
-
-    # root.input.model.unit_001.velocity = 1
-    # root.input.model.unit_001.cross_section_area = 4700.352526439483
-    root.input.model.unit_001.velocity = 5.75e-4
-
-    root.input.model.unit_001.adsorption.is_kinetic = 0
-    root.input.model.unit_001.adsorption.sma_ka = [0.0, 35.5]
-    root.input.model.unit_001.adsorption.sma_kd = [0.0, 1000.0]
-    root.input.model.unit_001.adsorption.sma_lambda = 1200.0
-    root.input.model.unit_001.adsorption.sma_nu = [0.0, 4.7]
-    root.input.model.unit_001.adsorption.sma_sigma = [0.0, 11.83]
-
-    root.input.model.unit_001.discretization.nbound = [1, 1]
-    root.input.model.unit_001.discretization.ncol = 100
-    root.input.model.unit_001.discretization.npar = 5
-
-    root.input.model.unit_002.ncomp = 2
-    root.input.model.unit_002.unit_type = 'OUTLET'
-
-    root.input.solver.user_solution_times = numpy.linspace(0, 1500, 1500)
-    root.input.solver.sections.nsec = 3
-    root.input.solver.sections.section_continuity = [0, 0]
-    root.input.solver.sections.section_times = [0.0, 10.0, 90.0, 1500.0]
-
-
-def plotSimulation(simulation):
-    f, (ax1, ax2) = plt.subplots(1, 2, figsize=[16, 8])
-    plotInlet(ax1, simulation)
-    plotOutlet(ax2, simulation)
-    f.tight_layout()
-    plt.savefig("tmp/lwe.png")
-
-
-def plotInlet(axis, simulation):
-    solution_times = simulation.root.output.solution.solution_times
-
-    inlet_salt = simulation.root.output.solution.unit_000.solution_inlet_comp_000
-    inlet_p1 = simulation.root.output.solution.unit_000.solution_inlet_comp_001
-    # inlet_p2 = simulation.root.output.solution.unit_000.solution_inlet_comp_002
-    # inlet_p3 = simulation.root.output.solution.unit_000.solution_inlet_comp_003
-
-    axis.set_title("Inlet")
-    axis.plot(solution_times, inlet_salt, 'b-', label="Salt")
-    axis.set_xlabel('time (s)')
-
-    # Make the y-axis label, ticks and tick labels match the line color.
-    axis.set_ylabel('mMol Salt', color='b')
-    axis.tick_params('y', colors='b')
-
-    axis2 = axis.twinx()
-    axis2.plot(solution_times, inlet_p1, 'r-', label="P1")
-    # axis2.plot(solution_times, inlet_p2, 'g-', label="P2")
-    # axis2.plot(solution_times, inlet_p3, 'k-', label="P3")
-    axis2.set_ylabel('mMol Protein', color='r')
-    axis2.tick_params('y', colors='r')
-
-    lines, labels = axis.get_legend_handles_labels()
-    lines2, labels2 = axis2.get_legend_handles_labels()
-    axis2.legend(lines + lines2, labels + labels2, loc=0)
-
-
-def plotOutlet(axis, simulation):
-    solution_times = simulation.root.output.solution.solution_times
-
-    outlet_salt = simulation.root.output.solution.unit_002.solution_outlet_comp_000
-    outlet_p1 = simulation.root.output.solution.unit_002.solution_outlet_comp_001
-    # outlet_p2 = simulation.root.output.solution.unit_002.solution_outlet_comp_002
-    # outlet_p3 = simulation.root.output.solution.unit_002.solution_outlet_comp_003
-
-    axis.set_title("Output")
-    axis.plot(solution_times, outlet_salt, 'b-', label="Salt")
-    axis.set_xlabel('time (s)')
-
-    # Make the y-axis label, ticks and tick labels match the line color.
-    axis.set_ylabel('mMol Salt', color='b')
-    axis.tick_params('y', colors='b')
-
-    axis2 = axis.twinx()
-    axis2.plot(solution_times, outlet_p1, 'r-', label="P1")
-    # axis2.plot(solution_times, outlet_p2, 'g-', label="P2")
-    # axis2.plot(solution_times, outlet_p3, 'k-', label="P3")
-    axis2.set_ylabel('mMol Protein', color='r')
-    axis2.tick_params('y', colors='r')
-
-    lines, labels = axis.get_legend_handles_labels()
-    lines2, labels2 = axis2.get_legend_handles_labels()
-    axis2.legend(lines + lines2, labels + labels2, loc=0)
-
-
-def test_lwe():
-    sim = main()
-    assert isinstance(sim.root.output.solution.unit_002.solution_outlet_comp_001, numpy.ndarray)
-    assert isinstance(sim.root.output.solution.unit_002.solution_outlet_comp_000, numpy.ndarray)
diff --git a/tests/test_meta_class_install_path.py b/tests/test_meta_class_install_path.py
new file mode 100644
index 0000000..e46e40e
--- /dev/null
+++ b/tests/test_meta_class_install_path.py
@@ -0,0 +1,31 @@
+from pathlib import Path
+import pytest
+
+from cadet import Cadet
+
+
+""" These tests require two distinct CADET installations to compare between and should not run in the CI"""
+
+# Full path to cadet.dll or cadet.so, that is different from the system/conda cadet
+full_path_dll = Path("path/to/cadet")
+
+install_path_conda = Cadet.autodetect_cadet()
+
+
+@pytest.mark.local
+def test_meta_class():
+    Cadet.cadet_path = full_path_dll
+    assert Cadet.use_dll
+
+    # With a path set in the meta class, the sim instance should not autodetect and use the meta class cadet path
+    sim = Cadet()
+    assert sim.use_dll
+    assert sim.install_path == full_path_dll.parent.parent
+    assert sim.cadet_dll_path == full_path_dll
+    assert sim.cadet_cli_path.parent == full_path_dll.parent
+
+    # With an install path given, the sim instance should use the given install path
+    sim = Cadet(install_path=install_path_conda)
+    assert sim.install_path == install_path_conda
+    assert sim.cadet_dll_path.parent.parent == install_path_conda
+    assert sim.cadet_cli_path.parent.parent == install_path_conda