Lombardy Young NUmerical Analysts (LYNUM VI)


The 6th edition of the "LYNUM - Lombardy Young NUmerical analysts Meeting" workshop

is scheduled for June 14, 2024,

from 10:00 am to 5:00 pm

at Politecnico di Milano, Aula Consiglio, Building 14 "La Nave" (Via Bonardi 9),

The aim of the meeting is to bring together young researchers in Numerical Analysis and Scientific Computing from universities in Lombardy to exchange their work and foster new connections and potential future collaborations.

In this edition, we will have the opportunity to listen to the work of the following young researchers:

  • Achini Federico (Università Statale di Milano)
  • Bignardi Paolo (Università di Pavia)
  • Capuano Emilia (Politecnico di Milano)
  • Dimola Nunzio (Politecnico di Milano)
  • Ghiotto Massimiliano (Università di Pavia)
  • Haile Kirubell Biniamin (Università di Milano Bicocca)
  • Molinari Loris (Università Statale di Milano)
  • Mosconi Marialetizia (Università di Milano Bicocca)
  • Preda Silvia (Università dell’Insubria)

The registration is free of charge but mandatory: you can register here!

See you at LYNUM VI!

9:45 - 10:00: Introduction

10:00 - 10:30: "Numerical modelling of polymer mixing processes: towards multi-phase simulation", E. Capuano

10:30 - 11:00: "Virtual element methods for the obstacle problem", L. Molinari

11:00 - 11:30: Coffee Break

11:30 - 12:00: "Space-time coercive variational formulation and adaptive methods for wave problems", P. Bignardi

12:00 - 12:30: "Least-square finite elements for parabolic problems", M. Mosconi

12:30 - 14:00: Lunch

14:00 - 14:30: "High order semi-Lagrangian schemes and applications", S. Preda

14:30 - 15:00: "Operator Learning for multi-patch domains", M. Ghiotto

15:00 - 15:30: Coffee Break

15:30 - 16:00: "Deep-Learning preconditioners for mixed-dimensional partial differential equations on 3D-1D domains", N. Dimola

16:00 - 16:30: "A Péclet-robust discontinuous Galerkin method for nonlinear diffusion with advection", K. B. Haile

Error minimization of the Non-uniform Fast Fourier Transform

Federico Achini*

*Dipartimento di Matematica "Federigo Enriques" - Università Statale di Milano

The FFT allows to compute the Fourier Transform of a finite signal with only \(O(N \log(N))\) operations, compared to \(O(N^2)\) of the naïve approach. However it requires data to be equally spaced, and this is rarely the case in real-world applications such as MRI and inverse problems in tomography. To overcome this issue, in 1993 Dutt and Rokhlin developed a generalization of FFT to non-equispaced data (NUFFT) based on interpolation by suitable functions. Since then, several improvements have been made to reduce the error size by finding the best interpolants. In this exposition, I’ll show the most significant ones, why they are important and their limits.

Space–time coercive variational formulation and adaptive methods for wave problems

Paolo Bignardi*

*Dipartimento di Matematica "Felice Casorati" - Università di Pavia

Space–time methods discretise the full space-time domain of definition of initial-boundary value problems (IBVP), contrary to the classical method of lines and Rothe’s method. For the acoustic wave equation \(\frac1{c^2}\partial_{tt}^2u-\Delta u=f\) many non-conforming space–time formulations have been proposed (e.g. DG), but no stable variational formulation is available that also allows the use of general discrete spaces. Following a strategy used for Helmholtz problems, we propose the use of a Morawetz-multiplier to derive a general space–time variational formulation for the second-order time-domain wave equation that is continuous and coercive in a certain Hilbert space \(V\) whose norm is stronger than \(H^1(\Omega\times(0,T))\).

From Lax Milgram theorem it follows that any discrete formulation set in \(V_h\subset C^1(\Omega\times(0,T))\) is well posed and Céa’s lemma provides error bounds. 

Because of the \(C^1\)-regularity requirement of the Galerkin discrete space, spline spaces are well suited to discretise this problem, as we show in numerical experiments that confirm the theoretical results.

The generality of the setting easily allows to derive an a-posteriori estimator for the \(V\)-norm of the error which predicts accurately the expected order of convergence and makes solving the IBVP computationally less expensive, allowing for space-time adaptive refinement.

Numerical modelling of polymer mixing processes: towards multi-phase simulation

Emilia Capuano*

*MOX, Dipartimento di Matematica - Politecnico di Milano

The numerical simulation of polymer mixing technologies is aimed at supporting the industrial reality for improving process efficiency and productivity. The mathematical and numerical modeling is continually improving in order to obtain increasingly reliable predictions of real industrial processes. The main challenges to be dealt with are the non-Newtonian rheologies of the materials involved, leading to highly nonlinear terms inside the Navier-Stokes equations, in addition to complex moving domain geometries consisting of rotating screws that mix and push forward the polymer compounds inside the machines.

This work falls within the framework of Immersed Boundary (IB) methods that handle complex solid boundaries inside the fluid domain, avoiding the need for a conformal grid with respect to the background fluid mesh.

Given a real mixing process, depending on the amount of polymer compound feeding the machine and on the screw velocity, the mixing device may be only partially filled with material. In order to account for this scenario, we focus on the integration of free-surface models [1] within the currently available immersed boundary method. Specifically, we rely on an Immersed Boundary - Volume Of Fluid (IB-VOF) model that is able to account for the presence of two materials inside a domain whose boundaries can be treated with the IB method.Appropriate wall boundary conditions are required to deal with contact lines, i.e. intersections of the biphase interface with solid walls. The Navier-slip boundary conditions [2] replace the standard no-slip, being able to impose a non-zero tangential velocity proportional to the normal stress at the boundary. Preliminary simulations and validations are carried out to test the implemented IB-VOF method combined with Navier-slip conditions, in comparison with either literature results or equivalent test cases resolved in the conformal framework, when possible.


[1]  S. Elgeti and H. Sauerland. "Deforming fluid domains within the finite element method Five mesh-based tracking methods in comparison". Archives of Computational Methods in Engineering 23, pp. 323-361, 2016.

[2]  L. L. Ferras, J. M. Nóbrega and F. T. Pinho. "Implementation of slip boundary conditions in the finite volume method: new techniques". International Journal for Numerical Methods in Fluids 72(7), pp. 724-747, 2013.

Deep-Learning preconditioners for mixed-dimensional partial differential equations on 3D-1D domains

Nunzio Dimola*, Nicola Rares Franco*, Paolo Zunino*

*MOX, Dipartimento di Matematica - Politecnico di Milano

Many physical problems that involve heterogeneous spatial scales, such as flow through fractured porous media, the study of fiber-reinforced materials, or the modeling of small circulation in living tissues, just to mention a few examples, can be described as mixed-dimensional problems. These are coupled partial differential equations (PDEs) defined in domains of heterogeneous dimensions that are embedded into each other. However, the algebraic structure of the resulting discretized system presents challenges in constructing efficient solution algorithms. In fact, the study of the interplay between the mathematical structure of the problem and solvers, as well as preconditioners for its discretization, is still in its infancy. The results presented in [4] for the solution of one-dimensional (1D) differential equations embedded in two-dimensional (2D), and more recently extended to the 3D-1D case in [3], have paved the way, but much more must be understood [2]. We address this challenge with a radically new approach originally proposed in [1]. We introduce a data-based technique to address the system resulting from the finite element discretisation of both 3D and 1D partial differential equations. In our method, we merge traditional iterative solvers with convolutional neural networks (CNNs) to create a preconditioned Krylov solver. The preconditioner employs the U-Net architecture, making use of CNN layers. Through our training approach, our CNN-based preconditioner is capable of generalizing across residuals coming from broad range of data that represent the functionality of 1D graphs within a 3D space. This work paves the way for the efficient solution of many-query problems where a mixed-dimensional problem is solved for many configurations of a microstructure shaped as a 1D/2D graph.


[1] Y. Azulay and E. Treister. "Multigrid-augmented deep learning preconditioners for the Helmholtz equation". SIAM Journal on Scientific Computing 45(3), pp. S127–S151, 2023.

[2] N. Dimola, M. Kuchta, K.-A. Mardal and P. Zunino. "Robust preconditioning of mixed-dimensional PDEs on 3D-1D domains coupled with Lagrange multipliers". ArXiv, 2023.

[3] M. Kuchta, K.-A. Mardal and M. Mortensen. "Preconditioning trace coupled 3D-1D systems using fractional Laplacian". Numerical Methods for Partial Differential Equations 35, pp. 375–393, 2019.

[4] M. Kuchta, M. Nordaas, J. Verschaeve, M. Mortensen and K.-A. Mardal. "Preconditioners for saddle point systems with trace constraints coupling 2D and 1D domains". SIAM Journal on Scientific Computing 38(6), pp. B962–B987, 2016.

Operator learning for multi-patch domains

Massimiliano Ghiotto*

*Dipartimento di Matematica "Felice Casorati" - Università di Pavia

Traditional neural networks development has primarily focused on learning mappings between finite-dimensional Euclidean spaces. Recently, this has been generalized to neural operators that learn mappings between function spaces. For partial differential equations, neural operators directly learn the solution operator that maps from the functional space of parametric functions to the functional space of solutions. This approach have the potential to accelerate traditional numerical methods when a mathematical description is known. When no model is available, this data-driven methodology has the potential to approximate the underlying input-output map. 

Several neural operator architectures have emerged in recent years. Fourier Neural Operators (FNOs) and the broader class of Spectral Neural Opeartors (SNOs) have received considerable interest due to their state-of-the-art performance across numerous tasks. 

In this talk we introduce the neural operator framework, underlining that FNOs and SNOs have intrinsic geometric limitations to rectangular domains uniformly discretized. Subsequently, we propose a novel framework to extend neural operators to multi-patch domains, addressing this limitation.

A Péclet-robust discontinuous Galerkin method for nonlinear diffusion with advection

Lourenço Beirão da Veiga*, Daniele A. Di Pietro**, Kirubell B. Haile*

*Dipartimento di Matematica e Applicazioni - Università di Milano-Bicocca

**Département de Mathématiques - Université de Montpellier

Discontinuous Galerkin (DG) methods were introduced in the 70s [1] and they are nowadays widely regarded as the reference methods for advection-dominated problems. When a polynomial degree \(k\ge 1\) is used, classical error estimates for linear diffusion-advection(-reaction) problems show that the error contribution stemming from diffusive terms is \(O(h^k)\) (with \(h\) denoting the meshsize), while the one stemming from advective terms is \(O(h^{k+\frac12})\).

The present talk, based on [2], aims to show new Péclet-dependent error estimates for a problem with linear advection-reaction and nonlinear \(p\)-type diffusion, with Sobolev indices \(p\in (1, \infty)\). Convergence analyses for various DG schemes applied to pure $p$-type diffusion setting can be found, e.g., in [3, 4]. To the authors' knowledge, an investigation of model problems including both advection and nonlinear diffusion is missing in the literature of DG elements. Especially if one aims at developing sharp estimates which respect the local nature of diffusion and convection,  the interaction between the (linear) advection and the nonlinear diffusion cannot be accounted for through a simple combination of known techniques, as the estimate of each term becomes dependent on the local regime.

The discretization of the nonlinear diffusion term is based on the full gradient including jump liftings and interior-penalty stabilization while, for the advective contribution, we consider a strengthened version of the classical upwind scheme. The peculiarity of our error estimates is that they track the dependence of the local contributions to the error on local Péclet numbers. In the linear case, corresponding to \(p = 2\), local Péclet numbers can be computed based on the sole knowledge of the problem data and the mesh, making it possible to identify a priori advection- and diffusion-dominated elements/faces. We emphasize that our results hold for general polygonal and polyhedral meshes, which we believe is an important asset. The present contribution furthermore sets the stage for future publications developing pressure robust and advection-robust finite elements for time-dependent Navier–Stokes type equations modeling incompressible fluid flows with non-Newtonian rheology.

In the present talk, after presenting the model and the numerical method, we will outline the theoretical results. Finally, a set of numerical tests supporting the theory will be shown.


[1] G. A. Baker. "Finite element methods for elliptic equations using nonconforming elements". Mathematics of Computation 31(137), pp. 45-49, 1977.

[2] L. Beirão da Veiga, D. A. Di Pietro and K. B. Haile. "A Péclet-robust discontinuous Galerkin method for nonlinear diffusion with advection". Mathematical Models and Methods in Applied Sciences, in press, 2024.

[3] E. Burman and A. Ern. "Discontinuous Galerkin approximation with discrete variational principle for the nonlinear Laplacian".  Comptes Rendus de l'Académie des Sciences Paris 346(17-18), pp. 1013–1016, 2008.

[4] T. Malkmus, M. Ruzicka, S. Eckstein and I. Toulopoulos. "Generalizations of SIP methods to systems with p-structure". IMA Journal of Numerical Analysis 38 (3), pp. 1420–1451, 2018.

Virtual element methods for the obstacle problem

Loris Molinari*

*Dipartimento di Matematica "Federigo Enriques" - Università Statale di Milano

Virtual Element Methods (VEM) are a recent family of numerical methods widely employed today for approximating partial differential equations. This class of methods naturally adapts to arbitrary polygonal decomposition of the domain, due to the choice of suitable discrete spaces that are no longer restricted to polynomials. Moreover, through the use of virtual functions and appropriate degrees of freedom, VEM can achieve high-order of accuracy without explicitly integrating non-polynomial functions.

In our work, we apply the lowest-degree virtual element method to address a specific constrained elliptic differential problem, commonly referred to as the obstacle problem. This problem, in its variational formulation, represents an interesting example of variational inequality - an important class of nonlinear problems frequently encountered in various fields of applied sciences.

In this talk, we will focus on both the primal and mixed variational formulations of the obstacle problem. We will begin by briefly presenting the mathematical framework of VEM. Subsequently, we will provide an explicit VEM discretization of the variational formulations, by introducing the employed discrete spaces and forms. In addition, for each formulation we will theoretically analyze  the convergence of the method and we will finally report a few computational results, which will be able to numerically validate the presented a priori error estimates. 


[1] L. Beirão Da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L.D. Marini and A.Russo. "Basic principles of virtual element methods". Mathematical Models and Methods in Applied Sciences 23, pp. 199-214, 2013.

[2] J. Haslinger, I. Hlaváček, J. Nečas. "Numerical methods for unilateral problems in solid mechanics". Handbook of Numerical Analysis 4, pp. 313-485, 1996.

[3] T. Gustafsson, R. Stenberg and J. Videman. "Mixed and stabilized finite element methods for the obstacle problem". SIAM Journal on Numerical Analysis 55(6), pp. 2718–2744, 2017.

Least-squares finite elements for parabolic problems

Marialetizia Mosconi*

*Dipartimento di Matematica e Applicazioni - Università di Milano-Bicocca

Compared to standard primal formulations, first-order system least-squares (FOSLS) finite element methods (FEM) have the advantage of leading to linear systems with a symmetric and positive definite coefficient matrix; moreover, they allow for the automatic design of error estimators that are reliable and efficient, with constants that are independent of the discretization parameters. On the other hand, they require more degrees of freedom, since an extra flux variable is employed, and may lead to linear systems that are more ill-conditioned than those in the corresponding primal formulation.

To the best of our knowledge, error estimators for primal formulations of space-time (XT) finite elements for parabolic problems are not analyzed in the literature and cannot be straightforwardly carried over from their elliptic counterparts. For this reason, a FOSLS XT FEM is analyzed in [1], where a priori and a posteriori error estimates are discussed, and a reliable and efficient error estimator is constructed; special simplicial meshes and Lagrangian elements are employed both for the primal and flux variables, which yields a loss of convergence for singular functions.

Improved convergence rates are derived in [2]; this is accomplished using prismatic meshes, tensor product elements in space and time, and Raviart-Thomas elements for the spatial flux. The global linear system is solved monolithically.

In this talk, we discuss the pros and cons of the approaches in [1, 2] and investigate whether there is room for proposing further improved formulations.


[1] T. Führer and M. Karkulik. "Space-time least-squares finite elements for parabolic equations". Computers & Mathematics with Applicaitions 92, pp. 27–36, 2021.

[2] G. Gantner and R. Stevenson. "Improved rates for a space-time FOSLS of parabolic PDEs". Numerische Mathematik 156(1), pp. 133–157, 2024.

High order semi-Lagrangian schemes and applications

Silvia Preda*, Elisabetta Carlini**, Roberto Ferretti***, Matteo Semplice*

*Dipartimento di Scienza e Alta Tecnologia - Università degli Studi dell'Insubria

**Dipartimento di Matematica - Università La Sapienza

***Dipartimento di Matematica e Fisica - Università degli Studi Roma Tre

We examine a high order numerical scheme for time-dependent first order Hamilton-Jacobi-Bellman equations. In particular we propose to combine a semi-Lagrangian scheme with a Central Weighted Non-Oscillatory (CWENO) reconstruction. Convergency is studied in the case of state- and time-independent Hamiltonians and numerical simulations are presented in oneand two space dimensions, also in more general cases, demonstrating superior performance in terms of CPU time gain by the Central version of the scheme with respect to the one coupled with traditional WENO reconstructions. This type of numerical scheme is then generalized to the case of second order evolutive Hamilton–Jacobi equation which arises in the level set formulation of mean curvature motion and, in particular, in the case of a curvature-related level set model first proposed by Zhao et al., to reconstruct unknown surfaces from an unorganized set of points. The problem follows a variational and partial differential equation (PDE) formulation with a curvature constraint that minimizes the surface area weighted by its distance from the point cloud. Level set methods are used in this framework to track the evolution of an initial surface and to find an implicit representation of the final shape. Among all the possible representations, we compute the signed distance function at least in the vicinity of the reconstructed surface. Regarding the numerical approximation of the solution, the use of the semi-Lagrangian scheme coupled with a local interpolator allows us to save computational costs, while also enhancing the parallel implementation of the algorithm. Numerical tests in two and three dimensions are presented to evaluate the quality of the approximated solution and the efficiency of the algorithm in terms of CPU times.


[1] E. Carlini, R. Ferretti, S. Preda and M. Semplice. "A CWENO large time-step scheme for Hamilton–Jacobi equations". Submitted, 2024.

[2] S. Preda and M. Semplice. "Surface reconstruction from point cloud using a semi-Lagrangian scheme with local interpolator". Submitted, 2024.

[3] M. Falcone and R. Ferretti. "Consistency of a large time-step scheme for mean curvature motion". Numerical Mathematics and Advanced Applications, pp. 495–502, 2003.

[4] H.-K. Zhao, S. Osher, B. Merriman and M. Kang. "Implicit and nonparametric shape reconstruction from unorganized data using a variational level set method". Computer Vision and Image Understanding 80 (3), pp. 295–314, 2000.


  • Corti Mattia (
  • Temellini Erika (




Registration Form

deadline 10-06-2024