We consider numerical methods for the efficient and stable solution of advanced, multidimensional partial integro-differential equations (PIDEs) and partial integro-differential complementarity problems (PIDCPs) arising in financial energy option valuation. Here the underlying uncertain factors, e.g. the electricity price, are modelled by exponential Lévy processes to account for the jumps that are often observed in the markets. These jumps give rise to the integralterm in the PIDEs and PIDCPs and is nonlocal. For the effective numerical solution, we investigate in this project operator splitting methods. In particular we study swing options which are a popular type of exotic energy options that have multiple exercise times. We develop novel, second-order operator splitting methods of the implicit-explicit (IMEX) and alternating direction implicit (ADI) kind. We analyse their fundamental properties of stability, consistency, monotonicity and convergence. The acquired theoretical results are validated by ample numerical experiments.