Magmatic systems and volcanoes are found throughout the entire world, yet the processes responsible for the ascent of magmatic melt, as well as the structure of these systems, are still poorly constrained due to a lack of direct observations. All hypotheses in this field are derived from indirect observations (geophysical surveys, fieldwork and petrological studies on fossil exhumed systems or at the surface of currently active volcanic systems) and robust physical and geochemical models need to be developed to validate, quantify or refine these ideas. We address these challenges in this thesis in two parts. First, we study the melt depletion of a magmatic reservoir connected to the surface by a weak conduit. Using a 2D model, we simplify the rheology and structure of the system to a few parameters and determine how sensitive the velocities within the conduit are to variations in these parameters. Several modes of transport are identified and translated into analytical scaling laws. We then apply these scaling laws to the 2021 eruption of the Cumbre Vieja volcano to constrain the structure of the magmatic system located beneath La Palma in the Canary Islands. Second, we develop tools to model magmatic systems in more detail and in a self-consistent manner across the lithosphere. One of these tools is an extension of the numerical continuum approach commonly used in geodynamics to model both shear and tensile plastic failure. This enables us to include dykes, the main form of melt transport in the elasto-plastic upper crust, in our models. The constitutive equations presented include compressible visco-elasto-plasticity with viscoplastic regularization and non-linear rheologies. A new yield function, adapted to work reliably in tensile conditions without introducing unphysical stress states into the model, is also presented. Another tool described in this thesis is MAGEMin, a Gibbs energy minimizer applied to igneous systems. MAGEMin is an efficient and highly scalable minimization package that opens up new possibilities for petrological applications as well as for use in conjunction with thermomechanical models. It uses a combination of linear programming, extended Partitioning Gibbs Energy and gradient-based minimization. The implementation of the thermodynamic dataset (Holland et al., 2018) was benchmarked against THERMOCALC.