# Space-Time IsoGeometric Analysis for nonlinear transient heat transfer problems

Joaquín Cornejo Fuentes, Giancarlo Sangalli, Mattia Tani, Thomas Elguedj, David Dureisseix, Arnaud Duval

► To cite this version:

Joaquín Cornejo Fuentes, Giancarlo Sangalli, Mattia Tani, Thomas Elguedj, David Dureisseix, et al.. Space-Time IsoGeometric Analysis for nonlinear transient heat transfer problems. 2025. hal-05038710

HAL Id: hal-05038710

https://hal.science/hal-05038710v1

Preprint submitted on 17 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Space-Time IsoGeometric Analysis for nonlinear transient heat
transfer problems

Joaquín Cornejo Fuentes¹, Giancarlo Sangalli², Mattia Tani², Thomas Elguedj¹, *, David Dureisseix¹, and Arnaud Duval¹

¹Univ Lyon, INSA Lyon, CNRS, LaMCoS, UMR5259, Villeurbanne, 69621, France

²Università di Pavia, Dipartimento di Matematica, Pavia, 27100, Italy

*Corresponding author: thomas.elguedj@insa-lyon.fr

April 17, 2025

Abstract

This paper presents a study of Space-Time IsoGeometric Analysis for solving nonlinear tran-
sient heat transfer problems, and overviews some numerical techniques to improve computational
efficiency. To deal with nonlinear materials, common in industrial processes, we implement ro-
bust nonlinear solvers such as Newton's method and Picard's method. Among the optimiza-
tion strategies, we adopt Matrix-Free and Weighted Quadrature approaches to cope with the
high computational costs and memory requirements associated with Space-Time methods. In
addition, we implement a Fast Diagonalization-based preconditioner to efficiently improve the
performance of the iterative solver. Numerical tests demonstrate the superior performance of
Space-Time methods over traditional incremental approaches, particularly when optimized al-
gorithms are used. Our analysis highlights as well as the advantages of inexact nonlinear solvers
in preventing oversolving and accelerating convergence. This work provides an overview of the
development of efficient and scalable algorithms for nonlinear heat transfer problems.

***Keywords-*** IsoGeometric Analysis, Weighted Quadrature, Matrix-Free, Fast Diagonalization

# 1 Introduction

This work investigates the application of Space-Time IsoGeometric Analysis (ST-IGA) for solving non-linear transient heat transfer problems. Traditionally, such Partial Differential Equations are addressed using the standard Finite Element Method for spatial discretization, which leads to a system of Ordinary Differential Equations in time. The resulting Initial Value Problem is typically solved using Finite Difference schemes that advance incrementally in time. While straightforward to implement, this approach imposes stability constraints on certain schemes, thereby limiting the time step size and reducing computational efficiency. Furthermore, the separation of time and space discretization complicates the development of adaptive techniques and parallel computing strategies. The inherently sequential nature of these methods poses significant challenges in designing efficient algorithms suited for large-scale parallel computing environments.

IsoGeometric Analysis, first introduced by Hughes *et al.* [1], represents a significant advancement over standard FEM. IGA employs spline functions or their generalizations both to represent the computational domain and to approximate solutions to PDEs. This integration of geometric modeling and numerical simulation facilitates the transition from Computer Aided Design to analysis. Moreover, the high continuity of splines offers superior accuracy compared to traditional *C*⁰ piecewise polynomials as stated by Bressan and Sande [2]. A promising alternative to conventional time stepping schemes is the Space-Time (ST) approach, which treats time as an additional dimension and discretizes it similarly to the spatial domain. This unified discretization strategy offers several advantages, including high accuracy in both space and time,
unconditional stability, and the removal of time step size restrictions (see Langer and Steinbach [3]). These features also enable the derivation of a *posteriori* error estimates, making ST methods particularly well-suited for adaptive algorithm design [4] and efficient parallel computing strategies [5].

The concept of ST methods for FEM originated in the late 1960s with applications to elastodynamics as proposed by Fried [6] and Oden [7], [8], and was later extended to heat conduction by Bruch Jr and Zyvoloski [9]. Despite their potential, these methods were not widely adopted at the time due to the high computational cost associated with increasing the dimensionality of the problem. Recent advances in computational power, however, have revitalized interest in ST methods as practical tools for solving complex problems. In recent years, significant progress has been made in the mathematical analysis of Galerkin methods for ST parabolic equations, including the work on wavelet discretization [10] and standard FEM discretization [11]. One of the earliest applications of ST-IGA for parabolic problems is found in [12], where a stabilized variational formulation was introduced, resulting in an elliptic discrete bilinear form under a specific energy norm. More recent works have proposed ST discontinuous Galerkin methods with stabilized formulations, which still yield sequential schemes [13]–[16].

This paper focuses on the nonlinear heat equation and its ST continuous discretization, which enables smooth approximation across both space and time. The method traces its origin to the Least-Squares formulation by Montardini et al. [17], followed by the exploration of ST continuous Galerkin formulations for linear parabolic problems by Loli et al. [18]. In this work, we extend the approach to account for nonlinear materials, as commonly encountered in industrial processes. Mathematically, this entails a variational formulation in which the constitutive matrix depends on the unknown temperature field. Various strategies exist for handling such nonlinear systems, but we focus on two primary techniques: Newton's method and Picard's method. Newton's method is renowned for its quadratic convergence, making it highly efficient for solving nonlinear problems, provided a good initial guess is available. In contrast, Picard's method, while generally slower in convergence, is valued for its robustness and simplicity [19]. It is worth noting that the use of smooth splines in the time domain introduces certain challenges, particularly regarding the causality principle. While our benchmark examples do not exhibit numerical instabilities, potential issues such as unphysical behavior or spurious oscillations propagating in time under unstable conditions merit further investigation. Addressing these concerns through stable schemes, such as the approach proposed by Loli et al. [20], is left for future work.

Additionally, we address the computational challenges associated with high memory usage and solver performance. To mitigate memory constraints, we adopt the Weighted Quadrature approach introduced by Calabrò et al. [21] and Hiemstra et al. [22], which aims to minimize the number of quadrature points by defining specific quadrature rules for each basis function. To improve solver efficiency, we abandon the traditional matrix formation and storage, opting instead for a Matrix-Free approach, as proposed by Sangalli and Tani [23]. In this approach, only the necessary components for computing matrix-vector products are stored, which is the minimal requirement for iterative solvers. By combining the Matrix-Free and Weighted Quadrature techniques, we achieve optimal algorithms. This is further enhanced by exploiting the tensor-product structure of multivariate basis functions. Rather than using the standard element-wise assembly, we apply a row-wise assembly alongside a sum-factorization algorithm, as proposed by Antolin et al. [24]. Since iterative solvers are preferred, it becomes essential to use a well-suited preconditioner. As noted by Collier et al. [25], standard preconditioners commonly employed in FEM perform poorly when applied to IGA. In our framework, we adopt the preconditioner proposed by Sangalli and Tani [26], which takes advantage of the tensor structure of multivariate functions. This preconditioner is based on Fast Diagonalization, a direct method developed for solving elliptic PDEs in the Finite Difference Method Lynch et al. [27]. It is particularly efficient for IGA, as its setup and application have a computational cost comparable to that of a matrix-vector product [26]. This ensures both high computational efficiency and robustness with respect to the mesh discretization.

Numerical tests were performed using Fortran libraries integrated within a Python framework, as part of the open-source IGA code YETI [28]. The results include several application cases where we analyze the precision-to-CPU time ratio for both incremental and ST methods, demonstrating the superior performance of the latter, particularly when optimized algorithms are employed. Additionally, we evaluated the ratio between the CPU time required for the matrix-vector product and that needed for applying the preconditioner. Our findings show that as the polynomial degree increases and the mesh is refined, this ratio grows, primarily
due to modern hardware's ability to efficiently handle matrix-matrix operations and the optimization of memory usage. Furthermore, we extended our analysis by comparing the performance of nonlinear solvers under inexact strategies. As expected, these solvers experience a reduction in optimal convergence rates, but this trade-off effectively mitigates oversolving, resulting in faster overall solution times. This highlights the practical advantages of inexact strategies in balancing computational efficiency and accuracy in nonlinear systems.

The structure of the paper is as follows. Sec. 2 introduces the nonlinear transient heat transfer problem, along with its IsoGeometric discretization, and discuss key concepts related to nonlinear solvers. Sec. 3 provides a concise overview of the Weighted Quadrature and Matrix-Free approaches, while Sec. 4 details the adaptation of the Fast Diagonalization method as a preconditioner for our problem. Sec. 5 presents numerical experiments that validate the theoretical results, and Sec. 6 summarizes the main findings and conclusions.

# 2 Nonlinear transient heat equation

## 2.1 Strong form

We consider the space-time computational domain $\Pi = \Omega \times \mathcal{T}$, where $\Omega \subset \mathbb{R}^d$ represents the spatial domain in $d$-dimensions, and $\mathcal{T} = [0, T] \subset \mathbb{R}^+$ denotes the time interval. The boundary of $\Omega$, denoted by $\Gamma$, is considered to be piecewise smooth. The temperature $u$ will explicitly depend on the position $x \in \Omega$ and the time $t \in \mathcal{T}$. Furthermore, as in most thermomechanical processes, the model is considered to have anisotropic heterogeneous thermal properties: the volumetric heat capacity $\sigma(u, x)$ and the thermal conductivity $k(u, x)$ may explicitly depend on the temperature $u$ and the position $x$. We assume that $\sigma(u, x)$ and the terms $k_{ij}(u, x)$ are Lipschitz continuous with respect to $u$. We also consider that the strictly positive scalar $\sigma(u, x)$ and the symmetric positive definite $d \times d$ matrix $k(u, x)$ are uniformly elliptic and bounded, i.e.,

there exist independent scalar parameters $\mu_1, \mu_2, \sigma_0, k_0 > 0$ such that:

$$ \begin{aligned} \mu_1 \|w\|^2 &\le \sigma(u, x)w^2, \quad \|\sigma(u, x)w\| \le \sigma_0 \|w\|, \quad \forall x \in \Omega, \forall u, w \in \mathbb{R}; \\ \mu_2 \|w\|^2 &\le w \cdot k(u, x)w, \quad \|k(u, x)w\| \le k_0 \|w\|, \quad \forall x \in \Omega, \forall u \in \mathbb{R}, \forall w \in \mathbb{R}^d. \end{aligned} \quad (1) $$

Then, the strong form of the heat problem is given by:

$$ (S) = \begin{cases} \text{Find } u : \Pi \to \mathbb{R} \text{ such that:} \\ \sigma \partial_t u - \nabla \cdot (k \nabla u) = f, & \text{on } \Omega \times \mathcal{T}; \\ u = 0, & \text{on } \Gamma \times \mathcal{T}; \\ u = 0, & \text{on } \Omega \times \{0\}. \end{cases} \quad (2) $$

From here onward, the notation $\partial_t$ represents the derivative with respect to time $t$, while $\nabla$ refers to the spatial differential operator, containing the derivatives with respect to spatial coordinates $x_i$. Additionally, we consider a heat source function $f \in L^2(\mathcal{T}; H^{-1}(\Omega))$ independent of the temperature field $u$, homogeneous Dirichlet boundary conditions and a zero initial temperature across the spatial domain.

## 2.2 Weak form and Galerkin method

We introduce the Hilbert spaces:

$$ \begin{aligned} \mathcal{X} &= \{w \in L^2(\mathcal{T}; H_0^1(\Omega)) \cap H^1(\mathcal{T}; H^{-1}(\Omega)) \mid w(x, 0) = 0\}; \\ \mathcal{Y} &= \{w \in L^2(\mathcal{T}; H_0^1(\Omega))\}. \end{aligned} $$

Then, the variational formulation of the problem is given by:

$$ (W) = \begin{cases} \text{Find } u \in \mathcal{X}, \text{ such that } \forall w \in \mathcal{Y}: \\ \mathcal{A}(u; u, w) = \mathcal{F}(w); \end{cases} \quad (3) $$
where the operators $\mathcal{A}$ and $\mathcal{F}$ are defined as

$$ \mathcal{A}(s; u, w) = \int_{\Pi} \partial_t u \sigma(s, \mathbf{x}) w d\Pi + \int_{\Pi} \nabla u \cdot [\mathbf{k}(s, \mathbf{x}) \nabla w] d\Pi; $$

$$ \mathcal{F}(w) = \int_{\Pi} f(\mathbf{x}, t) w d\Pi. $$

And, after constructing the finite-dimensional approximations $\mathcal{X}^h \subset \mathcal{X}$ and $\mathcal{Y}^h \subset \mathcal{Y}$, the Galerkin formulation is given by:

$$ (G) = \begin{cases} \text{Find } u^h \in \mathcal{X}^h, \text{ such that } \forall w^h \in \mathcal{Y}^h : \\ \mathcal{A}(u^h; u^h, w^h) = \mathcal{F}(w^h). \end{cases} \qquad (4) $$

The well-posedness and stability of the Galerkin formulation described by Eq. (4) have been studied for the linear case, particularly in the framework of ST-FEM as proposed by Steinbach [11]. However, for the nonlinear case, a rigorous theoretical analysis—including a formal proof of well-posedness and a priori stability estimates—remains, to the best of our knowledge, an open research question. Some progress has been made in related nonlinear ST formulations; for instance, Zhao and Li [4] provide a detailed analysis of a ST continuous Galerkin method for the nonlinear Sobolev equation, proving existence and uniqueness of the numerical solution as well as a priori error estimates using Lagrange polynomials. Extending such results to the nonlinear heat equation within the ST-IGA framework presents significant challenges due to the complexities of IsoGeometric discretization and the intricate numerical treatment of nonlinear terms. Instead of providing a formal proof, we assume the existence of a mathematical solution and focus on developing computational strategies within the ST-IGA continuous Galerkin setting, including linearization techniques, iterative solvers, and efficient quadrature schemes. Our findings indicate that the proposed numerical techniques yield stable and accurate approximations, offering insights into the behavior of the nonlinear heat equation in an ST variational formulation. These results contribute to the advancement of ST-IGA methods for nonlinear problems and lay the groundwork for future theoretical investigations into their well-posedness and stability.

## 2.3 The matrix formulation

Before continuing, we introduce briefly B-spline functions. In the univariate case, let $p$ and $n$ be positive integers, then an open knot-vector in the domain $[0, 1]$ is defined as $\Xi = \{\xi_1, \dots, \xi_{n+p+1}\}$, with $0 = \xi_1 = \dots = \xi_{p+1} < \xi_{p+2} \le \dots \le \xi_n < \xi_{n+1} = \dots = \xi_{n+p+1} = 1$, where knots can be repeated up to $m$-times. Here, $m$ defines the continuity of the B-spline, denoted by $r$, with $m \le p$ and $r = p-m$. The intervals $[\xi_A, \xi_{A+1}]$, with $A = p+1, \dots, n$, are referred to as knot-spans, provided that $\xi_A \ne \xi_{A+1}$, ensuring a nonzero length. To draw a connection with FEM, we also use the term “elements” interchangeably. Moreover, we denote by $h$ the mesh-size, i.e., $h = \max\{|\xi_{A+1} - \xi_A|, A = 1, \dots, n+p\}$. Then, from the knot-vector $\Xi$, we denote by $\hat{b}_A$, $A = 1, \dots, n$, the univariate B-splines basis functions of degree $p$ and continuity $r$ in the space $S_r^p = \text{span}\{\hat{b}_{A, p} | A = 1, \dots, n\}$ (see Cottrell et al. [29]). In the multivariate case, more specifically in the space-time context, we consider functions that depend on $d$ spatial variables and the time variable. Thus, given positive integers $n_l, p_l$ for $l = 1, \dots, d$ and $n_t, p_t$, we introduce $d+1$ univariate knot-vectors $\Xi_l = \{\xi_{l, 1}, \dots, \xi_{l, n_l+p_l+1}\}$ for $l = 1, \dots, d$ and $\Xi_t = \{\xi_{t, 1}, \dots, \xi_{t, n_t+p_t+1}\}$. Then, let us also define $\mathbf{p}$ as the vector that contains the degree indices, i.e., $\mathbf{p} = (\mathbf{p}_s, p_t)$, where $\mathbf{p}_s = (p_1, \dots, p_d)$. Multivariate B-splines are defined as the tensor product of the univariate B-splines, i.e.,

$$ \hat{N}_{\mathbf{A}, \mathbf{p}}(\xi_s, \xi_t) = \hat{B}_{\mathbf{A}_s, \mathbf{p}_s}(\xi_s) \hat{b}_{A_t, p_t}(\xi_t), \quad \text{with} \quad \hat{B}_{\mathbf{A}_s, \mathbf{p}_s}(\xi_s) = \hat{b}_{A_1, p_1}(\xi_1) \dots \hat{b}_{A_d, p_d}(\xi_d), \qquad (5) $$

$\mathbf{A}_s = (A_1, \dots, A_d)$, $\mathbf{A} = (\mathbf{A}_s, A_t)$ and $\xi_s = (\xi_1, \dots, \xi_d)$. From now on, we assume to have the same polynomial degree and discretization in all spatial dimensions, i.e., with abuse of notations, we set $p_1 = \dots = p_d = p_s$ and $n_1 = \dots = n_d = n_s$. As a consequence, we have the same mesh-size in all spatial dimensions, i.e., $h_1 = \dots = h_d = h_s$. For sake of simplicity, we will drop the degree index of B-spline notation.

As shown in Fig. 1, in our IGA framework, we consider that the spatial domain $\Omega$ is a single-patch spline parameterized by a conforming mapping $F_s$ such that $\mathbf{x} = F_s(\xi_s) = \sum_{A_s \in \eta_s} C_{A_s} \hat{B}_{A_s}$, with $C_{A_s} \in \mathbb{R}^d$
and $\xi_s \in \hat{\Omega} = [0, 1]^d$. Here, the scalar index $A_s$, obtained from the vector index $\mathbf{A}_s$, is in the set of indices of spatial control points $\eta_s$. We consider as well that the time interval $\mathcal{T}$ is also described by a B-spline curve parameterized by the conforming mapping $F_t$ such that $t = F_t(\xi_t) = \sum_{A_t \in \eta_t} C_{A_t} \hat{b}_{A_t}$, with $C_{A_t} \in \mathbb{R}$ and $\xi_t \in \hat{\mathcal{T}} = [0, 1]$. Here, the scalar index $A_t$ is in the set of indices of temporal control points $\eta_t$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Space-Time IGA paraphernalia. Picture adapted from Langer et al. [12].</figcaption></figure>

In this framework, we consider that there exists a conforming parameterization that maps the parametric domain $\hat{\Pi} = [0, 1]^{d+1}$ into the computational domain $\Pi = \Omega \times \mathcal{T}$. In the parametric space, we define $\xi \in \hat{\Pi}$ such that $\xi = (\xi_s, \xi_t) \in \hat{\Omega} \times \hat{\mathcal{T}}$. Then, we assume that finite-dimensional solution field $u^h$ may be discretized such as

$$u^h(\mathbf{x}, t) = \sum_{A \in \eta} N_A(\mathbf{x}, t)u_A,$$

where $N_A(\mathbf{x}, t) = \hat{N}_A(F_s^{-1}(\mathbf{x}), F_t^{-1}(t))$ and $u_A \in \mathbb{R}$ are coefficients that could be prescribed in order to respect the initial and boundary conditions. In addition, we define $\eta$ as the set of indices of all the control points in space and time, obtained from the combination of $\eta_s$ and $\eta_t$. Similarly, $\eta_g$ denotes the set of indices of inactive control points (without prescribed values). Hereafter, we denote as $n_{cp} = \text{card}(\eta) = n_1 \dots n_d n_t = n_s^d n_t$ the total number of control points. Finally, the matrix formulation is given by:

$$(M) = \left\{ \begin{array}{l} \text{Find } \mathbf{u} = [u_A]_{A \in \eta} \in \mathbb{R}^{n_{cp}}, \text{ such that:} \\ \sum_{B \in \eta} [\mathbf{A}(u^h)]_{AB} u_B = \mathbf{f}_A, \quad \forall A \in \eta - \eta_g; \\ \int_{\Pi} N_A \sigma(s, \mathbf{x}) \partial_t N_B d\Pi + \int_{\Pi} \nabla N_A \cdot [\mathbf{k}(s, \mathbf{x}) \nabla N_B] d\Pi; \\ \int_{\Pi} N_A f d\Pi. \end{array} \right. \quad (6)$$

where the operators are defined as follows:

$$[\mathbf{A}(s)]_{AB} = \int_{\Pi} N_A \sigma(s, \mathbf{x}) \partial_t N_B d\Pi + \int_{\Pi} \nabla N_A \cdot [\mathbf{k}(s, \mathbf{x}) \nabla N_B] d\Pi;$$

$$\mathbf{f}_A = \int_{\Pi} N_A f d\Pi.$$

## 2.4 Nonlinear methods

Among the various iterative methods for solving nonlinear problems, in this work, we focus exclusively on Newton's and Picard's methods. In particular, for the Newton's method [19], we define the terms of the nonlinear residual as follows:

$$r_A^{\text{NL}}(\mathbf{u}) = \mathbf{f}_A - \sum_{B \in \eta} [\mathbf{A}(u^h)]_{AB} u_B, \quad \forall A \in \eta - \eta_g, \quad (7)$$
recalling that the approximation solution is given by $u^h = \sum N_C u_C$, $\forall C \in \eta$. Starting with an initial guess $\mathbf{u}_{(0)}$, we construct iteratively a sequence of approximations $\mathbf{u}_{(k)}$ by solving a linearized problem at each step. Using a Taylor series expansion of Eq. (7) and omitting higher order terms, we obtain

$$ \mathcal{T}(u_{(k)}^h) \cdot (\mathbf{u}_{(k+1)} - \mathbf{u}_{(k)}) = \mathbf{r}_A^{\text{NL}}(\mathbf{u}_{(k)}). $$

Here, we consider that $u_{(k)}^h = \sum N_A[u_{(k)}]_A$, $\forall A \in \eta$, which allows us to express the (negative) tangent operator as $[\mathcal{T}(u^h)]_{AB} = \partial \mathbf{r}_A^{\text{NL}}(\mathbf{u})/\partial u_B$. With this formulation, the Newton update scheme follows as

$$ [u_{(k+1)}]_A = [u_{(k)}]_A + [\delta u]_A, \quad \text{where} \\ [\delta u]_A = \sum_{B \in \eta - \eta_g} [\mathcal{T}^{-1}(k)]_{AB} [\mathbf{r}_{(k)}^{\text{NL}}]_B \quad \text{for all } A \in \eta - \eta_g. \qquad (8) $$

For clarity, from now on, the subscript (k) indicates that a given quantity is evaluated at $u_{(k)}^h$, e.g., $\mathcal{T}_{(k)}$ denotes $\mathcal{T}(u_{(k)}^h)$. Moreover, the tangent matrix $\mathcal{T}$ is derived from Eqs. (6) and (7) as

$$ [\mathcal{T}(u^h)]_{AC} = [\mathbf{A}(u^h)]_{AC} + \sum_{B \in \eta} \left[ \frac{\partial \mathbf{A}(u^h)}{\partial u_C} \right]_{AB} u_B = [\mathbf{A}(u^h)]_{AC} + [\mathbf{B}(u^h)]_{AC} \quad \text{with } A, C \in \eta - \eta_g. $$

Here $[\mathbf{A}(u^h)]_{AC} = \int_{\Pi} N_A \sigma(u^h, \mathbf{x}) \partial_t N_C d\Pi + \int_{\Pi} \nabla N_A \cdot [\mathbf{k}(u^h, \mathbf{x}) \nabla N_C] d\Pi$ and

$$ \begin{aligned} [\mathbf{B}(u^h)]_{AC} &= \frac{\partial}{\partial u_C} \left[ \int_{\Pi} N_A \sigma(u^h, \mathbf{x}) \partial_t N_B d\Pi \right] u_B + \frac{\partial}{\partial u_C} \left[ \int_{\Pi} \nabla N_A \cdot (\mathbf{k}(u^h, \mathbf{x}) \nabla N_B) d\Pi \right] u_B \\ &= \left[ \int_{\Pi} N_A \frac{\partial \sigma}{\partial u^h} \frac{\partial u^h}{\partial u_C} \partial_t N_B d\Pi \right] u_B + \left[ \int_{\Pi} \nabla N_A \cdot \left( \frac{\partial \mathbf{k}}{\partial u^h} \frac{\partial u^h}{\partial u_C} \nabla N_B \right) d\Pi \right] u_B \\ &= \int_{\Pi} N_A \left( \frac{\partial \sigma}{\partial u^h} \partial_t u^h \right) N_C d\Pi + \int_{\Pi} \nabla N_A \cdot \left( \frac{\partial \mathbf{k}}{\partial u^h} \nabla u^h \right) N_C d\Pi. \end{aligned} $$

Consequently, we may write the various matrices in the parametric space as follows:

$$ \begin{aligned} [\mathbf{A}(\hat{u}^h)]_{AB} &= \int_{\hat{\Pi}} \hat{N}_A c(\hat{u}^h, \boldsymbol{\xi}_s) \hat{\partial}_\tau \hat{N}_B d\hat{\Pi} + \int_{\hat{\Pi}} \hat{\nabla} \hat{N}_A \cdot [D(\hat{u}^h, \boldsymbol{\xi}_s, \xi_t) \hat{\nabla} \hat{N}_B] d\hat{\Pi}, \\ \text{where } c(\hat{u}^h, \boldsymbol{\xi}_s) &= \hat{\sigma}(\hat{u}^h, \boldsymbol{\xi}_s) \det J_{F_s}(\boldsymbol{\xi}_s) \text{ and} \\ D(\hat{u}^h, \boldsymbol{\xi}_s, \xi_t) &= (J_{F_s}^{-1}(\boldsymbol{\xi}_s) \hat{\mathbf{k}}(\hat{u}^h, \boldsymbol{\xi}_s) J_{F_s}^{-\top}(\boldsymbol{\xi}_s)) \det J_{F_s}(\boldsymbol{\xi}_s) \det J_{F_t}(\xi_t); \end{aligned} \qquad (9a) $$

$$ \begin{aligned} [\mathbf{B}(\hat{u}^h)]_{AB} &= \int_{\hat{\Pi}} \hat{N}_A \tilde{c}(\hat{u}^h, \boldsymbol{\xi}_s) \hat{N}_B d\hat{\Pi} + \int_{\hat{\Pi}} [\hat{\nabla} \hat{N}_A \cdot \tilde{D}(\hat{u}^h, \boldsymbol{\xi}_s, \xi_t)] \hat{N}_B d\hat{\Pi}, \\ \text{where } \tilde{c}(\hat{u}^h, \boldsymbol{\xi}_s) &= \frac{\partial \hat{\sigma}}{\partial \hat{u}^h}(\hat{u}^h, \boldsymbol{\xi}_s) \hat{\partial}_\tau \hat{u}^h \det J_{F_s}(\boldsymbol{\xi}_s) \text{ and} \\ \tilde{D}(\hat{u}^h, \boldsymbol{\xi}_s, \xi_t) &= (J_{F_s}^{-1}(\boldsymbol{\xi}_s) \frac{\partial \hat{\mathbf{k}}}{\partial \hat{u}^h}(\hat{u}^h, \boldsymbol{\xi}_s) \hat{\nabla} \hat{u}^h J_{F_s}^{-1}(\boldsymbol{\xi}_s)) \det J_{F_s}(\boldsymbol{\xi}_s) \det J_{F_t}(\xi_t). \end{aligned} \qquad (9b) $$

In Eqs. (9a) and (9b), the notation $\hat{\partial}_\tau$ represents differentiation with respect to $\xi_t$, while $\hat{\nabla}$ contains the derivatives with respect to $\xi_i$, with $i = 1, \dots, d$. In general, the notation $\bullet$ indicates that the expression is evaluated in the parametric space rather than the physical space. Accordingly, we define $\hat{u}^h(\boldsymbol{\xi}_s, \xi_t) = u^h(F_s(\boldsymbol{\xi}_s), F_t(\xi_t))$, $\hat{\sigma}(\hat{s}, \boldsymbol{\xi}_s) = \sigma(s, F_s(\boldsymbol{\xi}_s))$, and $\hat{\mathbf{k}}(\hat{s}, \boldsymbol{\xi}_s) = \mathbf{k}(s, F_s(\boldsymbol{\xi}_s))$. Similarly, the Picard scheme as proposed by Reddy [19] can be outlined as

$$ \sum_{B \in \eta} [\mathbf{A}_{(k)}]_{AB} [u_{(k+1)}]_B = \mathbf{f}_A, \quad \text{with } A \in \eta - \eta_g. \qquad (10) $$
Here, according to the previously stated notation, $A_{(k)}$ denotes $A(u_{(k)}^t)$. Moreover, it is possible to rewrite Eq. (10) in order to obtain a Newton-like scheme where we neglect the contribution of the matrix B. In that case, in the nonlinear scheme Eq. (8), instead of considering a “full” tangent matrix $T = A + B$, we use a “partial” tangent matrix $T = A$. In Alg. 1, we explain in broad terms the different steps in order to solve the nonlinear problem. In this algorithm, the tangent matrix $T$ is defined as $T = A + \theta B$, where the parameter $\theta$ could be 0 (Picard’s scheme) or 1 (Newton’s scheme).

The nonlinear methods considered in this work have its advantages and disadvantages. On one hand, it can be shown theoretically that under suitable conditions the Newton scheme exhibits quadratic convergence (see Ortega and Rheinboldt [30]), while Picard converges only linearly. It is also important to remark that the global cost of both methods are quite different since Newton method needs to compute, in addition to matrix A, the matrix B. On the other hand, the Newton scheme is in practice quite sensitive to the initial solution estimate, i.e., it may diverge with a poor initial guess; while Picard scheme does not generally diverge. In order to exploit the best features of both methods, some authors propose to use the Picard scheme only for the first few iterations to improve the initial solution estimate, and then to switch to the Newton scheme for the last iterations to take advantage of the quadratic convergence rate. This approach is also known as Picard-Newton method. Finally, we propose the following stopping criterion for the nonlinear solver:

$$ \|\mathbf{r}_{(k)}^{\text{NL}}\|_2 \le \max \left( \epsilon_{\text{svg}}, \epsilon^{\text{NL}} \|\mathbf{r}_{(0)}^{\text{NL}}\|_2 \right), \quad \text{with } \epsilon_{\text{svg}} = 10^{-12} \text{ and } \epsilon^{\text{NL}} = 10^{-8}. \qquad (11) $$

<caption>**Algorithm 1: Nonlinear scheme**</caption>

**Input:** An initial guess $\mathbf{u}_{(0)}$ and $\theta \in \{0, 1\}$.
**Output:** A better guess of $\mathbf{u}$, written as $\mathbf{u}_{(k+1)}$.
1 **for** $k = 0, \dots, \text{maximum number of nonlinear iterations} \text{ do}$
2 &nbsp;&nbsp;&nbsp;Compute the residual $\mathbf{r}_{(k)}^{\text{NL}} = \mathbf{f} - A_{(k)}\mathbf{u}_{(k)}$.
3 &nbsp;&nbsp;&nbsp;**if** convergence **then** return
4 &nbsp;&nbsp;&nbsp;Solve the linearized system $T_{(k)}\delta\mathbf{u} = \mathbf{r}_{(k)}^{\text{NL}}$, where $T = A + \theta B$.
5 &nbsp;&nbsp;&nbsp;Update active nodes $\mathbf{u}_{(k+1)} = \mathbf{u}_{(k)} + \delta\mathbf{u}$.

In order to reduce the global computational cost of Alg. 1, in most cases, we would rather use an iterative method for solving the linear system in line 4. More specifically, we prefer Krylov methods (see [31] for further information) such as GMRES or BiCGstab, since $T_{(k)}$ is usually a huge nonsymmetric sparse matrix. These methods iteratively build a sequence $\delta\mathbf{u}_{(j)}$ that minimizes some norm of the linear (Krylov)-residual $\mathbf{r}_{(k,j)}^{\text{LK}} = \mathbf{r}_{(k)}^{\text{NL}} - T_{(k)}\delta\mathbf{u}_{(j)}$ over the affine space $\delta\mathbf{u}_{(0)} + \mathcal{K}_{(k,j)}$. Here $\mathcal{K}_{(k,j)}$, also known as Krylov-subspace, is a vector space defined as:

$$ \mathcal{K}_{(k,j)} = \operatorname{span} \left( \mathbf{r}_{(k,0)}^{\text{LK}}, \mathbf{r}_{(k,1)}^{\text{LK}}, \dots, \mathbf{r}_{(k,j-1)}^{\text{LK}} \right), \quad \text{with} \quad \mathbf{r}_{(k,l)}^{\text{LK}} = T_{(k)} \mathbf{r}_{(k,l-1)}^{\text{LK}}. $$

It can be shown theoretically that these vectors are linearly independent which is essential to prove the convergence of Krylov methods within a finite number of iterations [31]. The stopping criterion for the linear solver is written as:

$$ \|\mathbf{r}_{(k,j)}^{\text{LK}}\|_2 \le \epsilon^{\text{LK}} \|\mathbf{r}_{(k,0)}^{\text{LK}}\|_2. \qquad (12) $$

Usually, $\epsilon^{\text{LK}}$ is a fixed value of around $\epsilon_{\text{tight}}^{\text{LK}} = 10^{-8}$. Since it is a tight tolerance, we may say that the linear solver is “exact”, from an arithmetical point of view. In contrast, some authors propose an “inexact” approach which goal is to reduce drastically the cost of solving the linear system by limiting the number of iterations. As proposed by Dembo et al. [32], inexact methods can be described concisely as follows:

For some $\epsilon_{(k)}^{\text{LK}} \in (0, 1)$, find:

$$ \|\mathbf{r}_{(k)}^{\text{NL}} - T_{(k)}\delta\mathbf{u}\|_2 \le \epsilon_{(k)}^{\text{LK}} \|\mathbf{r}_{(k)}^{\text{NL}}\|_2. \qquad (13) $$

It must be noticed the condition in Eq. (13) is equivalent to Eq. (12). Thus, the key idea of “inexact” approach is to choose $\epsilon^{\text{LK}}$, also called “forcing term”, for every nonlinear iteration. In fact, the basic idea
is that when far from the solution, it is not necessary to compute an accurate linear solution, when getting closer to the solution, the linear solution must be precise in order to recover the convergence rate of the nonlinear method. Following the works of Birken [33] and Elman et al. [34], we set $\epsilon_{(k)}^{LK} = 0.1$ for the Picard method, as a relatively loose tolerance improves convergence. On the other hand, for the Newton method, following the work of Eisenstat and Walker [35] and Porcelli et al. [36], we propose that the forcing term follows the scheme:

$$ \text{For } k > 0, \gamma \in (0, 1], \zeta \in (1, 2], \text{ choose} \\ \epsilon_{(k)}^{LK} = \begin{cases} \epsilon_k, & \text{if } \gamma (\epsilon_{(k-1)}^{LK})^{\zeta} \le 0.1; \\ \max (\epsilon_k, \gamma (\epsilon_{(k-1)}^{LK})^{\zeta}), & \text{otherwise.} \end{cases} \qquad (14) $$

Here $\epsilon_k = \gamma (\|\mathbf{r}_{(k)}^{NL}\|_2 / \|\mathbf{r}_{(k-1)}^{NL}\|_2)^{\zeta}$. Under standard assumptions, the order of convergence of the nonlinear method is determined by the chosen exponent $\zeta$ (see [35][Theorem 2.3]). However, in practice, we set $\zeta < 2$ and use a loose initial threshold $\epsilon_{(0)}^{LK} = \epsilon_{loose}^{LK} \le 0.5$ to enhance performance. As with the nonlinear scheme, we add a safeguard to the stopping criterion of the linear solver, defined as:

$$ \|\mathbf{r}_{(k,j)}^{LK}\|_2 \le \max \left( \epsilon_{\text{svg}}, \max \left( \epsilon_{\text{tight}}^{LK}, \min \left( \epsilon_{\text{loose}}^{LK}, \epsilon_{(k)}^{LK} \right) \right) \|\mathbf{r}_{(k,0)}^{LK}\|_2 \right). \qquad (15) $$

As pointed out for many authors (see e.g. [37]–[39]), the choice of forcing term is very important since it not only determines the asymptotic speed of convergence to a solution but also affects the efficiency and robustness of the algorithm. When a “good” forcing term is chosen, we avoid oversolving, i.e., solving the linear system with unnecessary precision without decreasing the nonlinear residual in the same proportion of the linear one. Finally, taking into the way of computing the tangent matrix and the precision of the linear method, we consider the methods: Exact Newton, Exact Picard, Inexact Newton and Inexact Picard.

# 3 Matrix-Free approach

In our theoretical framework, we have to deal with matrix-vector products such as Au in line 2 or even Tδu within the linear solver in line 4 of Alg. 1. In both cases, the formation and storage of the matrices may be completely avoided by providing a routine that takes a vector v and returns the application Av or Tv. As a consequence, we may reduce drastically the computational resources in terms of CPU time and memory usage. Since, depending on the nonlinear scheme chosen, the tangent matrix could be A or A + B, we detail in the next paragraphs how to build the Matrix-Free algorithms for Av and Bv. However, before continuing, we introduce the mathematical tools that make it possible for such algorithms to be implemented.

## 3.1 Kronecker product and sum-factorization

Let $A \in \mathbb{R}^{n_1 \times n_2}$ and $B \in \mathbb{R}^{n_3 \times n_4}$. The Kronecker product between A and B is defined as

$$ A \otimes B = \begin{bmatrix} A_{11}B & \cdots & A_{1n_2}B \\ \vdots & \ddots & \vdots \\ A_{n_11}B & \cdots & A_{n_1n_2}B \end{bmatrix} \in \mathbb{R}^{n_1n_3 \times n_2n_4}. $$

The following properties of the Kronecker product are essential for understanding this work:

* it holds $(A \otimes B)^{\top} = A^{\top} \otimes B^{\top}$;

* if A, B, C and D are matrices of conforming order, then it holds $(A \otimes B) \cdot (C \otimes D) = (AC) \otimes (BD)$;

* if A and B are invertible, then $(A \otimes B)^{-1} = A^{-1} \otimes B^{-1}$.

Furthermore, we introduce the l-mode product $\times_l$ for $l \in \{1, \dots, d\}$, which represents the multiplication of a tensor $V \in \mathbb{R}^{n_1 \times \dots \times n_d}$ by a matrix $J \in \mathbb{R}^{m \times n_l}$. The elements of the resulting tensor are given by

$$ [V \times_l J]_{i_1 \dots i_d} = \sum_{j=1}^{n_l} J_{il} \cdot V_{i_1 \dots i_{l-1} j i_{l+1} \dots i_d}. \qquad (16) $$
Then, as introduced by Kolda and Bader [40], given $J_l \in \mathbb{R}^{m_l \times n_l}$, with $l = 1, \dots, d$, we may define the following property:

$$ (J_d \otimes \dots \otimes J_1) \cdot \mathrm{vec}(\mathcal{V}) = \mathrm{vec}(\mathcal{V} \times_1 J_1 \times_2 \dots \times_d J_d), \quad (17) $$

where the vectorization operator “vec” applied to a tensor stacks its entries into a column vector (unfolding operation), following a row-major ordering, i.e., $[\mathrm{vec}(\mathcal{V})]_i = [\mathcal{V}]_{i_1 \dots i_d}$, where $i$ is computed using the standard row-major index mapping. The great advantage of this property, also called “sum-factorization” in FEM context, is that the matrix $J_d \otimes \dots \otimes J_1$ does not need to be formed explicitly, hence the name “Matrix-Free”. In Alg. 2, we summarize the procedure to apply Eq. (17). Moreover, as mentioned by Sangalli and Tani [23], the complexity analysis reveals that the total number of floating-point operations (FLOPs), when $m_l = m$ and $n_l = n$ for $l = 1, \dots, d$, is $2mn(m^{d-1} + nm^{d-2} + \dots + n^{d-1})$ when the matrices are dense; and is bounded by $2 \max_l (nnz(J_l)) (m^{d-1} + nm^{d-2} + \dots + n^{d-1})$ when the matrices are sparse, where nnz($J_l$) denotes the number of nonzero terms of $J_l$.

### Algorithm 2: Sum-factorization technique

**Input:** $J_1, \dots, J_d$, where $J_l \in \mathbb{R}^{m_l \times n_l}$ and $v_{in} \in \mathbb{R}^{n_1 \dots n_d}$

**Output:** $v_{out} \in \mathbb{R}^{m_1 \dots m_d}$

1. Tensorize the vector $v_{in}$ such as $V_0 = \mathrm{vec}^{-1}(v_{in})$.

2. **for** $l = 1, \dots, d$ **do**

3. &nbsp;&nbsp;&nbsp;&nbsp;Initialize $V_l \in \mathbb{R}^{n_1 \times \dots \times n_{l-1} \times m \times n_{l+1} \times \dots \times n_d}$.

4. &nbsp;&nbsp;&nbsp;&nbsp;Compute the tensor-matrix product (or $l$-mode product) $V_l = V_{l-1} \times_l J_l$ using Eq. (16).

5. &nbsp;&nbsp;&nbsp;&nbsp;Clear or deallocate $V_{l-1}$.

6. Vectorize the tensor $V_d$ such as $v_{out} = \mathrm{vec}(V_d)$.

## 3.2 Quadrature techniques

To understand how to apply the “sum-factorization” technique to our problem, we must first focus on the way the different terms of the matrix are computed. In what follows we assume that the integrals in Eqs. (9a) and (9b) are estimated by means of a quadrature formula that exploits the tensor structure of basis functions. The main reason for this choice is that it is no longer necessary to use element-by-element loops, but instead use row-by-row loops that allow parallel processing. In our numerical benchmarking, for almost all examples we have used the standard Gauss-Legendre (GL) quadrature. However, some other approaches such as Gauss-Galerkin [41] or Gauss-Greville quadrature [42] are conceivable. Additionally, we have also included some results using Weighted Quadrature (WQ) approach [22]. Further information about these techniques will be provided in the next paragraphs. First, let us consider the univariate case where the integrals to be computed are of the kind

$$ \int_0^1 \hat{b}_{A_l}^{(\alpha_l)}(\xi_l)g(\xi_l)d\xi_l, $$

where $\alpha_l$ is a parameter that could be 1 or 0 depending on the function’s nature, i.e., $\bullet^{(0)} = \bullet$ and $\bullet^{(1)} = \partial\bullet/\partial\xi_l$. Here $\hat{b}_{A_l}$ is the univariate function in the $l$-dimension and $g(\xi_l)$ is a function that may depend, among others, on the geometry mapping and material properties. In the standard Gauss-Legendre technique, it is possible to estimate the integral by:

$$ \int_0^1 \hat{b}_{A_l}^{(\alpha_l)}(\xi_l)g(\xi_l)d\xi_l \approx \sum_{q_l \in Q_l} \omega_{l, q_l} \hat{b}_{A_l}^{(\alpha_l)}(\tilde{\xi}_{q_l})g(\tilde{\xi}_{q_l}). \quad (18) $$

Here, we define $Q_l$ as the set of indices of quadrature points following the $l$-dimension which are denoted as $\tilde{\xi}_{q_l}$. In addition, $\omega_{l, q_l}$ denotes the quadrature weight at the quadrature knot $\tilde{\xi}_{q_l}$. Both of them are usually found in lookup tables (see e.g. [43]). A well-known property of this technique is that it guarantees the exact result if the integrand is a polynomial of degree $(2p_l + 1)$ or less over the $l$-dimension and at least $(p_l + 1)$ quadrature points are used for each element. In the Weighted-Quadrature approach, unlike Eq. (18), the integral is estimated by:

$$ \int_0^1 \hat{b}_{A_l}^{(\alpha_l)}(\xi_l)g(\xi_l)d\xi_l \approx \sum_{q_l \in Q_l} \omega_{l, A_l, q_l} g(\tilde{\xi}_{q_l}). \quad (19) $$
In this case, the trial function $\hat{b}_{A_l}^{(\alpha_l)}$ is part of the integral measure $(\hat{b}_{A_l}^{(\alpha_l)}(\xi_l)d\xi_l)$. As a consequence, the price to pay is that the quadrature weights depend now on the basis function while the quadrature points does not. However, as mentioned by Hiemstra et al. [22], the great advantage with respect to standard Gauss-Legendre technique is that one may choose a convenient set of quadrature points for their problem. In most of the cases, the number of quadrature points will be independent of the degree and as few as possible and as long as the Schoenberg-Whitney interpolation theorem [44] is satisfied. Furthermore, since the quadrature weights are not known analytically, they are computed numerically as solution of the exactness conditions:

$$ \int_0^1 \hat{b}_{A_l}^{(\alpha_l)}(\xi_l) \bar{b}_{B_l}(\xi_l) d\xi_l = \sum_{q_l \in Q_l} \omega_{l, A_l, q_l}^{(\alpha_l)} \bar{b}_{B_l}(\tilde{\xi}_{q_l}), \qquad (20) $$

where $\bar{b}_{B_l}$ are the test functions in the space $S_{r_l-1}^{p_l}$. It is worth noticing that Eq. (20) is a linear system with respect to the weights, and a vast literature exists to resolve it efficiently and accurately. In our benchmark, the number of quadrature points per element (far from the boundaries) is considered to be 3. As proposed by Calabrò et al. [21], we may reformulate the strategy the weights are defined to further reduce the number of quadrature points per element to 2. However, it is beyond the scope of this paper, and we refer to [21] for more information on the latter approach.

Finally, in the multivariate case, since the B-spline functions are the tensor product of univariate B-spline functions, the multivariate quadrature rule may be constructed as the tensor product of univariate quadrature rules. For example, let us consider the integral

$$ \int_{\hat{\Pi}} \hat{N}_A^{(\alpha)}(\xi) g(\xi) d\hat{\Pi}, \text{ with } \alpha = (\alpha_1, \dots, \alpha_d, \alpha_t). $$

Then, using Weighted Quadrature approach, it may be estimated by

$$ \int_{\hat{\Pi}} \hat{N}_A^{(\alpha)}(\xi) g(\xi) d\hat{\Pi} \approx \sum_{q_l \in Q_l} \omega_{l, A_l, q_l}^{(\alpha_l)} \sum_{q_d \in Q_d} \omega_{d, A_d, q_d}^{(\alpha_d)} \cdots \sum_{q_1 \in Q_1} \omega_{1, A_1, q_1}^{(\alpha_1)} g(\tilde{\xi}_{q_1}, \dots, \tilde{\xi}_{q_d}, \tilde{\xi}_{q_t}). \qquad (21) $$

Actually, Eq. (21) may be interpreted as an expansion of Eq. (17) considering the tensor $\mathcal{G}$ with $\mathcal{G}_{q_1, \dots, q_d, q_t} = g(\tilde{\xi}_{q_1}, \dots, \tilde{\xi}_{q_d}, \tilde{\xi}_{q_t})$, and the matrices $W_l^{(\alpha_l)}$, for the $l$-dimension, which terms are defined by

$$ \left[ W_l^{(\alpha_l)} \right]_{A_l, q_l} = \omega_{l, A_l, q_l}^{(\alpha_l)}, \quad \text{for } A_l = 1, \dots, n_l. $$

The same remarks can be drawn when using Gauss-Legendre quadrature in Eq. (18). In that case, the terms of the matrices are given by

$$ \left[ W_l^{(\alpha_l)} \right]_{A_l, q_l} = \hat{b}_{A_l}^{(\alpha_l)}(\tilde{\xi}_{q_l}) \omega_{l, q_l}, \quad \text{for } A_l = 1, \dots, n_l. $$

We could even go on in this reflection and say that in Gauss-Legendre the matrix $W_l^{(\alpha_l)}$ is the resulting product of the sparse matrix $B_l^{(\alpha_l)}$ and the diagonal matrix $W_l = \text{diag}(\{\omega_{l, q_l}\})$. Since it will be used later in our discussion, we define the terms of $B_l^{(\alpha_l)}$ as

$$ \left[ B_l^{(\alpha_l)} \right]_{q_l, A_l} = \hat{b}_{A_l}^{(\alpha_l)}(\tilde{\xi}_{q_l}), \quad \text{for } A_l = 1, \dots, n_l. $$

Furthermore, it is important to note that, as a consequence of the tensor structure of the basis functions and the quadrature rules, the mesh grid of quadrature points in the multivariate case is constructed from the combination of the univariate quadrature points, i.e., $\tilde{\xi}_q = (\tilde{\xi}_{q_1}, \dots, \tilde{\xi}_{q_d}, \tilde{\xi}_{q_t})$. Hereafter the scalar index $q$ is obtained from the set of indices of quadrature points $\mathcal{Q} = \{(q_1, \dots, q_d, q_t) | q_t \in Q_t, q_l \in Q_l, \text{ for } l = 1, \dots, d\}$. In the particular case of Gauss-Legendre, the weight of a quadrature point is also computed as $\omega_q = \omega_{1, q_1} \dots \omega_{d, q_d} \omega_{t, q_t}$.
## 3.3 Matrix-Free algorithms in Space-Time parabolic problems

### 3.3.1 MF calculations for Av

For a given vector $v \in \mathbb{R}^{n_{cp}}$, we want to compute the output vector $Av \in \mathbb{R}^{n_{cp}}$. Then, as explained by Sangalli and Tani [23], using Eq. (9a) we may compute the A-th component of $Av$ by:

$$ 
\begin{aligned}
[A(\bullet)v]_A &= \sum_{B \in \eta} [A(\bullet)]_{AB} v_B \approx \sum_{B \in \eta} \sum_{q \in Q} \omega_q \left[ \hat{N}_A c(\bullet, \tilde{\xi}_q) \hat{\partial}_\tau \hat{N}_B + \hat{\nabla} \hat{N}_A \cdot \left( D(\bullet, \tilde{\xi}_q) \hat{\nabla} \hat{N}_B \right) \right] v_B \\
&\approx \sum_{q \in Q} \omega_q \left[ \hat{N}_A c(\bullet, \tilde{\xi}_q) \sum_{B \in \eta} \hat{\partial}_\tau \hat{N}_B v_B + \hat{\nabla} \hat{N}_A \cdot \left( D(\bullet, \tilde{\xi}_q) \sum_{B \in \eta} \hat{\nabla} \hat{N}_B v_B \right) \right].
\end{aligned} 
\quad (22) $$

To write Eq. (22) in a compact form, we introduce the diagonal matrices $C = \text{diag}(\{C(\bullet, \tilde{\xi}_q)\})$ and $D_{lm} = \text{diag}(\{D_{lm}(\bullet, \tilde{\xi}_q)\})$. Then, we write

$$ 
\begin{aligned}
A(\bullet)v &\approx (W_t^{(0)} \otimes W_d^{(0)} \otimes \dots \otimes W_1^{(0)}) C (B_t^{(1)} \otimes B_d^{(0)} \otimes \dots \otimes B_1^{(0)}) \cdot v \\
&\quad + \sum_{l=1}^{d} \sum_{m=1}^{d} (W_t^{(0)} \otimes W_d^{(\alpha_d)} \otimes \dots \otimes W_1^{(\alpha_1)}) D_{lm} (B_t^{(0)} \otimes B_d^{(\beta_d)} \otimes \dots \otimes B_1^{(\beta_1)}) \cdot v,
\end{aligned} 
$$

where $\alpha_k, \beta_k$ depend on the nested loop piloted by $l, m$. Finally, we compute $Av$ by means of Alg. 3.

<table><caption>**Algorithm 3:** Compute $v_{out} = Av_{in}$</caption>
  <thead>
    <tr>
      <th>Input:</th>
      <th>W<sub>l</sub><sup>(α<sub>l</sub>)</sup>, D<sub>lm</sub>, B<sub>l</sub><sup>(β<sub>l</sub>)</sup></th>
      <th>for l, m = 1, ..., d, t and α<sub>l</sub>, β<sub>l</sub> ∈ {0, 1}; and v<sub>in</sub></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>Output:</td>
      <td>v<sub>out</sub></td>
      <td></td>
    </tr>
    <tr>
      <td>1</td>
      <td>Compute v&#x0303; = (B<sub>t</sub><sup>(1)</sup> &otimes; B<sub>d</sub><sup>(0)</sup> &otimes; &hellip; &otimes; B<sub>1</sub><sup>(0)</sup>) &middot; v<sub>in</sub> using Alg. 2.</td>
      <td></td>
    </tr>
    <tr>
      <td>2</td>
      <td>Compute v&#x0303; = Cv&#x0303;.</td>
      <td></td>
    </tr>
    <tr>
      <td>3</td>
      <td>Compute v<sub>out</sub> = (W<sub>t</sub><sup>(0)</sup> &otimes; W<sub>d</sub><sup>(0)</sup> &otimes; &hellip; &otimes; W<sub>1</sub><sup>(0)</sup>) &middot; v&#x0303; using Alg. 2.</td>
      <td></td>
    </tr>
    <tr>
      <td>4</td>
      <td>for m = 1, ..., d do</td>
      <td></td>
    </tr>
    <tr>
      <td>5</td>
      <td>Set β<sub>k</sub> = 0 for k = 1, ..., d, and update β<sub>m</sub> = 1.</td>
      <td></td>
    </tr>
    <tr>
      <td>6</td>
      <td>Compute v&#x0303; = (B<sub>t</sub><sup>(0)</sup> &otimes; B<sub>d</sub><sup>(β<sub>d</sub>)</sup> &otimes; &hellip; &otimes; B<sub>1</sub><sup>(β<sub>1</sub>)</sup>) &middot; v<sub>in</sub> using Alg. 2.</td>
      <td></td>
    </tr>
    <tr>
      <td>7</td>
      <td>for l = 1, ..., d do</td>
      <td></td>
    </tr>
    <tr>
      <td>8</td>
      <td>Initialize α<sub>k</sub> = 0 for k = 1, ..., d, and update α<sub>l</sub> = 1.</td>
      <td></td>
    </tr>
    <tr>
      <td>9</td>
      <td>Compute v&#x0303; = D<sub>lm</sub>v&#x0303;.</td>
      <td></td>
    </tr>
    <tr>
      <td>10</td>
      <td>Compute v<sub>out</sub> = v<sub>out</sub> + (W<sub>t</sub><sup>(0)</sup> &otimes; W<sub>d</sub><sup>(α<sub>d</sub>)</sup> &otimes; &hellip; &otimes; W<sub>1</sub><sup>(α<sub>1</sub>)</sup>) &middot; v&#x0303; using Alg. 2.</td>
      <td></td>
    </tr>
  </tbody>
</table>

### 3.3.2 MF calculations for Bv

Applying the previous ideas and using Eq. (9b), we may compute the A-th component of $Bv$ by:

$$ [B(\bullet)v]_A \approx \sum_{q \in Q} \omega_q \left[ \hat{N}_A \tilde{c}(\bullet, \tilde{\xi}_q) + \hat{\nabla} \hat{N}_A \cdot \tilde{D}(\bullet, \tilde{\xi}_q) \right] \sum_{B \in \eta} \hat{N}_B v_B. \quad (23) $$

We also introduce the diagonal matrices $\tilde{C} = \text{diag}(\{\tilde{C}(\bullet, \tilde{\xi}_q)\})$ and $\tilde{D}_l = \text{diag}(\{\tilde{D}_l(\bullet, \tilde{\xi}_q)\})$. Then, we may write Eq. (23) in a compact form as

$$ 
\begin{aligned}
B(\bullet)v \approx & \left\{ (W_t^{(0)} \otimes W_d^{(0)} \otimes \dots \otimes W_1^{(0)}) \tilde{C} + \sum_{l=1}^{d} (W_t^{(0)} \otimes W_d^{(\alpha_d)} \otimes \dots \otimes W_1^{(\alpha_1)}) \tilde{D}_l \right\} \\
& \left( B_t^{(0)} \otimes B_d^{(0)} \otimes \dots \otimes B_1^{(0)} \right) \cdot v.
\end{aligned} 
$$

Here the indices $\alpha_k$ depend on the loop piloted by $l$. Finally, we compute $Bv$ by means of Alg. 4.
**Algorithm 4:** Compute $v_{out} = Bv_{in}$

Input: $W_l^{(\alpha_l)}$, $\tilde{D}_l$, $B_l^{(0)}$ for $l = 1, \dots, d$ and $\alpha_l \in \{0, 1\}$; $W_t^{(0)}$, $\tilde{C}$, $B_t^{(0)}$; and $v_{in}$

Output: $v_{out}$

1 Compute $\tilde{v} = (B_t^{(0)} \otimes B_d^{(0)} \otimes \dots \otimes B_1^{(0)}) \cdot v_{in}$ using Alg. 2.

2 Compute $\tilde{\tilde{v}} = \tilde{C}\tilde{v}$.

3 Compute $v_{out} = (W_t^{(0)} \otimes W_d^{(0)} \otimes \dots \otimes W_1^{(0)}) \cdot \tilde{\tilde{v}}$ using Alg. 2.

4 for $l = 1, \dots, d$ do

5 Initialize $\alpha_k = 0$ for $k = 1, \dots, d$, and update $\alpha_l = 1$.

6 Compute $\tilde{\tilde{v}} = \tilde{D}_l \tilde{v}$.

7 Compute $v_{out} = v_{out} + (W_t^{(0)} \otimes W_d^{(\alpha_d)} \otimes \dots \otimes W_1^{(\alpha_1)}) \cdot \tilde{\tilde{v}}$ using Alg. 2.

## 3.4 Computational cost

We briefly analyze the cost of Alg. 3 and 4. For the sake of simplicity, we assume that space and time are discretized using the same polynomial degree and mesh-size, i.e., $p_s = p_t =: p$, $h_s = h_t := h$ and $n_s = n_t :=: n$. Moreover, we use the same quadrature rule for all the dimensions of the computational domain. Therefore, let us assume that $W_l^{(\alpha_l)}$ and $W_t^{(\alpha_t)}$ are sparse matrices of size $n \times rn$ with $prn$ nonzero entries; $C$, $\tilde{C}$, $\tilde{D}_l$ and $D_{lm}$ are diagonal matrices of size $(rn)^{d+1} \times (rn)^{d+1}$; $B_l^{(\beta_l)}$ and $B_t^{(\beta_t)}$ are sparse matrices of size $rn \times n$ with $prn$ nonzero entries. Here, $r$ is the ratio between the number of quadrature points and control points following one dimension. According to Sangalli and Tani [23], when computing a matrix-vector product with the sparse matrix having a tensor structure, like in line 1 of Alg. 3, its computational cost is bounded by $2p(rn)^{d+1}$ FLOPs. Approximately the same number of operations is required for steps 3, 6 and 10 of Alg. 3. Lines 2 and 8 of Alg. 3 are diagonal scaling steps whose computational cost is $(rn)^{d+1}$ FLOPs for each of them. In addition, if $D$ is considered nonsymmetric, the cost of the two loops in Alg. 3 is roughly $d(d+1)/2$ times the cost of its three first lines, otherwise it could be reduced by half [23]. The total computational cost of $Av$ will roughly be $(d(d+1)/2 + 1)(4p+1)(rn)^{d+1}$ FLOPs. On the other hand, the cost of the first three lines of Alg. 4 should be bounded by $(2p+1)(rn)^{d+1}$. Additionally, the cost of the loop in Alg. 4 is roughly $d/2$ times the cost of its first three lines. The total computational cost of $Bv$ will roughly be $(d/2+1)(4p+1)(rn)^{d+1}$ FLOPs. Theoretically, the matrix-vector product $Av$ is approximately $d$-times the cost of computing $Bv$. Furthermore, it is important to note that the total cost of $Av$ and $Bv$ are in both cases $O(pr^{d+1}n_{cp})$ FLOPs. Regarding memory usage, it is proportional to the number of quadrature points, since just $O(r^{d+1}n_{cp})$ points are required. It is important to remark that, if we combine MF and WQ approaches, the requirements in CPU time and memory usage will be respectively $O(pn_{cp})$ and $O(n_{cp})$. However, if we rely on standard Gauss-Legendre quadrature, we would require $O(p^{d+2}n_{cp})$ FLOPs to compute the matrix-vector product. Note that this is even more costly than standard product which is around $O(p^{d+1}n_{cp})$ FLOPs.

## 4 Preconditioning based on Fast Diagonalization

In order to increase the convergence rate of the linear solver, it is necessary to find an efficient preconditioner for the matrix $T$. The preconditioner we propose in this work is based on Fast Diagonalization technique which was introduced by Sangalli and Tani [26]. This method takes advantage of the tensor structure of the matrices and allows to efficiently find their inverse. It is important to note that we have adopted the following theoretical developments because they can be more easily integrated into our IGA code. However, it should be recognized that even more parallel algorithms with lower asymptotic costs have been developed in the works of Montardini et al. [45].
## 4.1 FD with partial inclusion of geometry and material properties

We propose the following preconditioner:

$$[\mathbf{P}]_{AB} = \int_{\hat{\Pi}} \hat{N}_A \tilde{c}_\sigma \hat{\partial}_\tau \hat{N}_B d\Pi + \int_{\hat{\Pi}} \hat{\nabla} \hat{N}_A \tilde{C}_K \hat{\nabla} \hat{N}_B d\Pi, \qquad (24)$$

with $\tilde{C}_K = \text{diag}\{\tilde{c}_{K,1}, \dots, \tilde{c}_{K,d}\}$. Here, based on the work of Brubeck and Farrell [46], the scalar parameters are computed by:

$$\tilde{c}_{\sigma} = \int_{\hat{\Pi}} \hat{\sigma} \det J_{F_s} d\hat{\Pi} \quad \text{and} \quad \tilde{C}_{K} = \text{diag} \left( \int_{\hat{\Pi}} \left( J_{F}^{-1} \hat{k} J_{F}^{\top} \right) \det J_{F_s} \det J_{F_t} d\hat{\Pi} \right).$$

Next, by exploring the tensor product of the basis functions, the matrix P described in Eq. (24) can be expressed for a 3D geometry as follows:

$$\mathbf{P} = \tilde{c}_{\sigma} \mathbf{L}_{t} \otimes \mathbf{M}_{3} \otimes \mathbf{M}_{2} \otimes \mathbf{M}_{1} + \mathbf{M}_{t} \otimes [\tilde{c}_{K,1} \mathbf{M}_{3} \otimes \mathbf{M}_{2} \otimes \mathbf{K}_{1} + \tilde{c}_{K,2} \mathbf{M}_{3} \otimes \mathbf{K}_{2} \otimes \mathbf{M}_{1} + \tilde{c}_{K,3} \mathbf{K}_{3} \otimes \mathbf{M}_{2} \otimes \mathbf{M}_{1}];$$

where,

$$[\mathbf{L}_t]_{A_t B_t} = \int_0^1 \hat{b}_{A_t}(\xi_t) \hat{b}_{B_t}'(\xi_t) d\xi_t, \quad [\mathbf{M}_t]_{A_t B_t} = \int_0^1 \hat{b}_{A_t}(\xi_t) \hat{b}_{B_t}(\xi_t) d\xi_t,$$

$$[\mathbf{K}_t]_{A_t B_t} = \int_0^1 \hat{b}_{A_t}'(\xi_t) \hat{b}_{B_t}'(\xi_t) d\xi_t, \quad [\mathbf{M}_t]_{A_t B_t} = \int_0^1 \hat{b}_{A_t}(\xi_t) \hat{b}_{B_t}(\xi_t) d\xi, \text{ for } l = 1, \dots, d.$$

As proposed by Sangalli and Tani [26], we may solve the generalized eigenvalue problem for the pencils $(\mathbf{K}_l, \mathbf{M}_l)$ in order to express them like:

$$\mathbf{K}_l = \mathbf{U}_l \Lambda_l \mathbf{U}_l^\top \quad \text{and} \quad \mathbf{M}_l = \mathbf{U}_l \mathbf{U}_l^\top,$$

where $\mathbf{U}_l$ and $\Lambda_l$ are respectively the eigenvectors and eigenvalues following the $l$-dimension. Each matrix $\mathbf{U}_l$ is real $\mathbf{M}_l$-orthogonal, i.e., $\mathbf{U}_l^{-1} = \mathbf{U}_l^\top$. On the other hand, for the pencil $(\mathbf{L}_t, \mathbf{M}_t)$ we may consider the generalized Schur decomposition instead since $\mathbf{L}_t$ is not symmetric (nor skew-symmetric). Then,

$$\mathbf{L}_t = \mathbf{U}_t \mathbf{X}_t \mathbf{V}_t^\dagger \quad \text{and} \quad \mathbf{M}_t = \mathbf{U}_t \mathbf{Y}_t \mathbf{V}_t^\dagger,$$

where $\mathbf{X}_t$ and $\mathbf{Y}_t$ are complex upper triangular matrices and, $\mathbf{U}_t$ and $\mathbf{V}_t$ are unitary matrices, i.e., $\mathbf{U}_t^{-1} = \mathbf{U}_t^\dagger$ and $\mathbf{V}_t^{-1} = \mathbf{V}_t^\dagger$. Hereafter, $\bullet^\dagger$ is the Hermitian conjugate of $\bullet$. Thus, inverse of preconditioner is given by

$$\begin{gathered} \mathbf{P}^{-1} = (\mathbf{U}_t \otimes \mathbf{U}_d \otimes \dots \otimes \mathbf{U}_1) \mathbf{Q}^{-1} (\mathbf{V}_t^\dagger \otimes \mathbf{U}_d^\top \otimes \dots \otimes \mathbf{U}_1^\top), \quad \text{with} \\ \mathbf{Q} = \tilde{c}_\sigma \mathbf{X}_t \otimes \mathbb{I}_d \otimes \dots \otimes \mathbb{I}_1 + \mathbf{Y}_t \otimes \mathbf{\Lambda} \quad \text{and} \quad \mathbf{\Lambda} = \left( \sum_{l=1}^d \tilde{c}_{K,l} \mathbb{I}_d \otimes \dots \otimes \mathbf{\Lambda}_l \otimes \dots \otimes \mathbb{I}_1 \right). \end{gathered} \qquad (25)$$

In Eq. (25), the matrix $\mathbf{Q}$ is a main-diagonal block matrix, and its inverse $\mathbf{Q}^{-1}$ is composed of the inverse of its blocks. In that way, multithreading is conceivable. Then, applying $\mathbf{Q}^{-1}$ is the same as solving Eq. (26) for every block:

$$\mathbf{z}_{out} = (\tilde{c}_{\sigma} \mathbf{X}_{t} + \lambda \mathbf{Y}_{t})^{-1} \mathbf{z}_{in}, \qquad (26)$$

where $\mathbf{z}_{in}, \mathbf{z}_{out} \in \mathbb{R}^{n_t}$ and $\lambda$ is an eigenvalue from the diagonal matrix $\mathbf{\Lambda}$. In Alg. 5, we summarize the procedure of applying our preconditioner to a vector.

**Algorithm 5: Compute $\mathbf{v}_{out} = \mathbf{P}^{-1}\mathbf{v}_{in}$**

**Input:** The Golub2013 of pencils $(\mathbf{K}_l, \mathbf{M}_l)$ for $l = 1, \dots, d$, the Schur decomposition of pencil $(\mathbf{L}_t, \mathbf{M}_t)$, and $\mathbf{v}_{in}$

**Output:** $\mathbf{v}_{out}$

1. Compute $\tilde{\mathbf{z}} = (\mathbf{V}_t^\dagger \otimes \mathbf{U}_d^\top \otimes \dots \otimes \mathbf{U}_1^\top) \cdot \mathbf{v}_{in}$ using Alg. 2.

2. Solve each triangular block in $\tilde{\mathbf{z}} = \mathbf{Q}^{-1}\tilde{\mathbf{z}}$.

3. Compute $\mathbf{v}_{out} = (\mathbf{U}_t \otimes \mathbf{U}_d \otimes \dots \otimes \mathbf{U}_1) \cdot \tilde{\mathbf{z}}$ using Alg. 2.
## 4.2 Computational cost

We briefly analyze the cost of Alg. 5. For that purpose, we assume that all the matrices $K_l, M_l, L_t$ and $M_t$ are the same size $n \times n$. As explained by Sangalli and Tani [26], the cost of the Golub2013 and Schur decomposition depends on the algorithm used. We recommend Golub and Loan [47] for more information on the different existing methods. In our implementation, we have used the function `dsygvd` from the LAPACK library (see Anderson et al. [48]) which computes the eigenvalues and eigenvectors of a real generalized symmetric definite eigenproblem using a divide and conquer algorithm whose cost is bounded by $O(n^3)$ FLOPs. In addition, we have used the function `zgges`, also from the LAPACK library, which computes the Schur decomposition in $O(n^3)$ FLOPs. By this approach, the cost of input requirements is bounded by $O(n^3)$ FLOPs. Lines 1 and 3 involve a dense matrix-vector product with a matrix having a tensor structure. Using Alg. 2, each line is equivalent to performing $d+1$ products between dense matrices of size $n \times n$ and $n \times n^d$. Then, the cost of lines 1 and 3 is roughly $4(d+1)n^{d+2}$ FLOPs. Line 2 is the same as solving $n^d$ times the linear system $z_{out} = Q_t^{-1}z_{in}$, where $Q_t \in \text{span}(X_t, Y_t)$ is upper triangular since $X_t$ and $Y_t$ are also upper triangular. Then, the cost of line 2 is roughly $n^{d+2}$ FLOPs. The overall computational cost of Alg. 5 is therefore $(4d+5)n^{d+2}$ or $O(n^{d+2})$ FLOPs which does not depend on the B-spline degree. In the general case, even if we assume an equal parametrization in all dimensions, i.e., $p_s = p_t = p$, $h_s = h_t = h$, the size of the matrices involved in the preconditioner will still depend on the boundary conditions, more precisely, by the number of active control points. However, as a good estimate, we can approximate $n^{d+1} \approx n_{cp}$, where $n_{cp}$ is the total number of control points. Under this assumption, the final cost scales as $O(n_{cp}^{1+1/(d+1)})$.

# 5 Results

All the results presented here are performed using a laptop equipped with Intel i7-11850H processors running at 2.5 GHz and with 16 GB of RAM. The different tests were programmed in Python, version 3.8.10. However, the main algorithms presented throughout the paper, nowadays part of the open-source IGA code YETI [28], were programmed in Fortran90, compiled into a static library using a GNU Fortran compiler, version 9.4.0, and called in Python using `f2py` module. Moreover, in our implementation, we used functions from the LAPACK library, version 3.9.0. In our benchmark, parallelization has not been considered and all the computations are restricted to one computational thread.

## 5.1 Convergence analysis

Consider that we study a domain $\Omega = \{x = (x_1, x_2) | 0 \le x_1, 0 \le x_2, 0.25^2 \le x_1^2 + x_2^2 \le 1^2\}$ with material properties $\sigma = 1$ and

$$k(u) = (3 + 2 \tanh(u/50)) \begin{bmatrix} 1 & 0.5 \\ 0.5 & 2 \end{bmatrix},$$

during the time interval $T = [0, 1]$. The heat source imposed is computed analytically such as the manufactured solution $u$ is given by $u(x, t) = 50 \tanh(1 - x_1^2 - x_2^2) \sin(\pi(x_1^2 + x_2^2 - 0.25^2)) \sin(\pi x_1 x_2) g(t)$, with $g(t) = \sin(\pi t/2)(1 + 0.75 \cos(3\pi t/2))$. Moreover, we consider a zero temperature over the domain at $t = 0$ and homogeneous Dirichlet boundary conditions over $\Gamma \times T$. Figs. 2a and 2b show respectively the temperature field and the thermal conductivity over the considered geometry $\Omega$ at the end of the interval $T$. Since space and time variables can be separated in the exact solution field considered in this example, the temperature over $\Omega$ at a given time $t$ might be seen as a scaling of Fig. 2a. However, this is not longer the case for the thermal conductivity because of the nonlinear function in its definition.

Then, we solve the nonlinear problem using the Exact Picard method in a set of successively refined meshes, with $h_s = h_t = h$, for different degrees varying from 1 to 5, with $p_s = p_t = p$. In Fig. 3 we present $u - u^h$, which is the error between the exact and approximate solution, in the $L^2(\Pi)$-norm, i.e.,

$$\|u - u^h\|_{L^2(\Pi)}^2 = \int_{\Pi} (u - u^h)^2 d\Pi.$$

Moreover, we considered relevant to show the differences between the quadrature methods: standard Gauss-Legendre (ST-IGA-GL), Weighted Quadrature proposed by Calabrò et al. [21] (ST-IGA-WQ-1) and Weighted
Quadrature proposed by Hiemstra et al. [22] (ST-IGA-WQ-2). As mentioned in Sec. 3.2, ST-IGA-WQ-2 uses 3$^{d+1}$ quadrature points for every inner element, while ST-IGA-WQ-1, due to its definition, may use 2$^{d+1}$. As a reference, ST-IGA-GL uses $(p+1)^{d+1}$ quadrature points per element. The results for ST-IGA-GL confirm the theoretical estimate of the optimal convergence rate $O(h^{p+1})$. Moreover, for this example, no loss of accuracy is observed when using the Weighted Quadrature method ST-IGA-WQ-2. Unfortunately, this is not the case for the other method ST-IGA-WQ-1. The reason is probably the lack of sufficient quadrature points to capture the nonlinear behavior. In our example, since $d=2$, the total number of quadrature points used by ST-IGA-WQ-1 represents at most 30% of the total used in ST-IGA-WQ-2, and only 3.7% of the total used in ST-IGA-GL when $p=5$. Although, this feature makes it less accurate, it also makes it faster because of the reduction of the FLOPs as we show later on our discussion.

<figure><img src="image_3.png" /><figcaption>Figure 2: Graphic representation of the exact temperature and the thermal conductivity at $t=1$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: $L^2(\Pi)$-norm of $u - u^h$ using Gauss-Legendre and Weighted Quadrature methods.</figcaption></figure>
## 5.2 Comparison of incremental and Space-Time methods

The nonlinear heat transfer problem could be addressed using a conventional method for nonlinear parabolic PDEs: the predictor-corrector method (see e.g. [49]). The main idea is to use a time-integration scheme, such as the generalized-$\theta$ method, to obtain a first guess of the solution at a given time step (predictor phase) and to use until convergence a nonlinear solver, such as Picard's Method, to update this estimate (corrector phase). Henceforth, since this method solves step-by-step the whole ST problem, we will refer to it as “incremental” method. Its popularity, at least in engineering, is mainly due to its intuitive and low-cost algorithmic implementation if compared to ST approach. However, as we discuss below, its lack of higher continuity in the time dimension means that, to obtain sufficiently accurate results, it is necessary to use a tiny time step size.

<figure><img src="image_5.png" /><figcaption>Figure 4: Evolution of the relative error $u - u^h$ using Crank-Nicolson scheme and Space-Time approach for a given degree $p_t$ and mesh-size $h_t$.</figcaption></figure>

As an example, we consider the same nonlinear problem above to be solved by an incremental approach and, for comparative purposes, by its ST counterpart. The former, also denoted by INC-IGA-GL, uses a Crank-Nicolson scheme ($\theta = 0.5$) with a time step size of 1/16; while the latter, denoted by ST-IGA-GL, is obtained by fixing the time mesh-size to $h_t = 1/16$ and the time degree to $p_t = 1$. Moreover, to avoid any loss of accuracy, we consider that both methods use standard Gauss-Legendre quadrature. Then, in Fig. 4a, we show the $L^2(\Pi)$-norm of the error using different spatial mesh-size $h_s$ and different spatial degree $p_s$. As shown for both methods, despite the refinement, it is not possible to bring the error below a certain threshold. This behavior is simply explained by the fact that the error due to time discretization dominates over spatial discretization one.

It is important to note that, for this example, at least one order of magnitude of error could be gained by using ST method instead of incremental approach. The underperformance of the incremental method is due to a combination of at least two aspects: the error made in calculating the initial conditions and the error accumulated and propagated at each time step. As a consequence, the only way to lower the threshold is refining the time dimension: using smaller time step-size $h_t$ or using higher degree $p_t$. For example, in Fig. 4b, we use a time step size of $h_t = 1/32$ for both methods and $p_t = 2$ for ST approach. As expected, the threshold or the minimum error in the $L^2(\Pi)$-norm is lower for both cases and continues to be independent of the spatial discretization.
Another way to reach the same remarks is to calculate the minimum error that can be obtained for different time refinements with a given spatial discretization. In Fig. 5, considering a fine spatial discretization where $p_s = 8$ and $h_s = 1/64$, we show the evolution of the minimum error as function of the time degree $p_t$ and time mesh-size $h_t$ for incremental and ST methods. It is important to note that, for comparison purposes, instead of the number of elements, we consider the number of degrees of freedom in the time dimension, also known as the number of control points in time or $n_t$. Then, the results for INC-IGA-GL confirm the quadratic convergence rate as theoretically assumed for Crank-Nicolson, whereas for ST-IGA-GL, the optimal convergence rate $O(h_t^{p_t+1})$ is achieved. From Fig. 5, the most important finding is that an accurate solution could be obtained using ST approach, with a coarse time discretization $h_t$ but with a high degree $p_t$, rather than using an incremental approach with a lot of time steps. Or, from another point of view, the solution could be some orders of magnitude more accurate for the same number of degrees of freedom in time when using the ST approach instead of the incremental method.

<figure><img src="image_6.png" /><figcaption>Figure 5: Evolution of the relative error for Crank-Nicolson scheme and Space-Time methods with different degrees $p_t$ on time.</figcaption></figure>

A final question to consider is how expensive one method is relative to the others; and the answer is in the analysis of Fig. 6. As in Fig. 3, for a given mesh-size with $h_s = h_t = h$, we solve the nonlinear problem for different degrees considering $p_s = p_t = p$. However, in this case we present the evolution of the error with respect to the total running time instead of the number of elements. This procedure is performed for ST-IGA-GL and ST-IGA-WQ-2. Also for comparative reasons, despite the differences between the methods, we perform simulations using the incremental method, with a time step size of $h_t = h$, for a given spatial mesh-size $h_s = h$ and different spatial degrees $p_s = p$. In addition, we apply the incremental method, denoted by INC-IGA-WQ-2, with the improvement algorithms described in [50]. Briefly, the idea is to use Weighted Quadrature together with Matrix-Free and Fast Diagonalization algorithms.

From Figs. 6a and 6b, purely analyzing the most advantageous combination of precision and computational cost (in running time), the best method is ST-IGA-WQ-2. This remark derives from the simple fact that the running time for ST-IGA-WQ-2 is practically the same, at least the same order of magnitude, whatever the polynomial degree $p$. Therefore, accurate solutions could be obtained for almost the same price. As we will present latter, this is actually a consequence of the independence of the number of quadrature points with respect to the degree. To prove this point it is only necessary to look at the curve obtained using the standard Gauss-Legendre quadrature (ST-IGA-GL). It is worth mentioning that ST-IGA-GL uses 2 and
3 quadrature points per element when $p = 1$ and $p = 2$ respectively; and that is the reason why it has a better performance than ST-IGA-WQ-2 which uses always 3 quadrature points per element. Moreover, while it is true that INC-IGA-WQ-2 is slightly faster than ST-IGA-WQ-2, at least in the case presented, it is far less accurate. Of course, there are cases in which the lack of quadrature points penalizes the accuracy of the solution; but in any case the results presented in this section show that using the correct algorithms we can have a better performance than with conventional methods. In fact, an alternative approach to presenting the previous results in Fig. 6 is shown in Fig. 7 and Tables 1 and 2, which is inspired by the work published by Sangalli and Tani [23].

<figure><img src="image_7.png" /><figcaption>Figure 6: Precision vs. running time (in seconds) for Crank-Nicolson scheme and Space-Time methods.</figcaption></figure>

<table><caption>Table 1: Relative $L^2$ error for Matrix-Free Weighted Quadrature [22] across different polynomial degrees and mesh-sizes.</caption><thead><tr><th rowspan="2">Degree</th><th colspan="5">Number of elements</th></tr><tr><th>4</th><th>8</th><th>16</th><th>32</th><th>64</th></tr></thead><tbody><tr><td>1</td><td>1.15 &times; 10<sup>-1</sup></td><td>3.37 &times; 10<sup>-2</sup></td><td>8.49 &times; 10<sup>-3</sup></td><td>2.11 &times; 10<sup>-3</sup></td><td>5.28 &times; 10<sup>-4</sup></td></tr><tr><td>2</td><td>5.06 &times; 10<sup>-2</sup></td><td>3.75 &times; 10<sup>-3</sup></td><td>4.07 &times; 10<sup>-4</sup></td><td>4.90 &times; 10<sup>-5</sup></td><td>6.12 &times; 10<sup>-6</sup></td></tr><tr><td>3</td><td>1.06 &times; 10<sup>-2</sup></td><td>7.09 &times; 10<sup>-4</sup></td><td>3.11 &times; 10<sup>-5</sup></td><td>1.77 &times; 10<sup>-6</sup></td><td>1.10 &times; 10<sup>-7</sup></td></tr><tr><td>4</td><td>5.73 &times; 10<sup>-3</sup></td><td>1.19 &times; 10<sup>-4</sup></td><td>3.78 &times; 10<sup>-6</sup></td><td>1.12 &times; 10<sup>-7</sup></td><td>3.43 &times; 10<sup>-9</sup></td></tr><tr><td>5</td><td>8.48 &times; 10<sup>-4</sup></td><td>3.44 &times; 10<sup>-5</sup></td><td>3.41 &times; 10<sup>-7</sup></td><td>4.54 &times; 10<sup>-9</sup></td><td>6.88 &times; 10<sup>-11</sup></td></tr><tr><td>6</td><td>4.76 &times; 10<sup>-4</sup></td><td>5.82 &times; 10<sup>-6</sup></td><td>3.83 &times; 10<sup>-8</sup></td><td>2.36 &times; 10<sup>-10</sup></td><td>1.60 &times; 10<sup>-11</sup></td></tr><tr><td>7</td><td>1.54 &times; 10<sup>-4</sup></td><td>1.82 &times; 10<sup>-6</sup></td><td>4.46 &times; 10<sup>-9</sup></td><td>1.56 &times; 10<sup>-11</sup></td><td>5.84 &times; 10<sup>-14</sup></td></tr><tr><td>8</td><td>1.17 &times; 10<sup>-5</sup></td><td>3.74 &times; 10<sup>-7</sup></td><td>5.80 &times; 10<sup>-10</sup></td><td>7.52 &times; 10<sup>-13</sup></td><td>1.78 &times; 10<sup>-15</sup></td></tr></tbody></table>
<table><caption>Table 2: Total computational time (in seconds) for Matrix-Free Weighted Quadrature [22] across different polynomial degrees and mesh-sizes.</caption><thead><tr><th rowspan="2">Degree</th><th colspan="5">Number of elements</th></tr><tr><th>4</th><th>8</th><th>16</th><th>32</th><th>64</th></tr></thead><tbody><tr><td>1</td><td>1.41 × 10<sup>-1</sup></td><td>5.47 × 10<sup>-1</sup></td><td>4.23</td><td>2.11 × 10<sup>1</sup></td><td>5.06 × 10<sup>2</sup></td></tr><tr><td>2</td><td>1.56 × 10<sup>-1</sup></td><td>9.69 × 10<sup>-1</sup></td><td>5.02</td><td>4.64 × 10<sup>1</sup></td><td>5.43 × 10<sup>2</sup></td></tr><tr><td>3</td><td>2.50 × 10<sup>-1</sup></td><td>1.14</td><td>6.19</td><td>5.38 × 10<sup>1</sup></td><td>5.92 × 10<sup>2</sup></td></tr><tr><td>4</td><td>3.91 × 10<sup>-1</sup></td><td>1.30</td><td>7.78</td><td>6.45 × 10<sup>1</sup></td><td>6.57 × 10<sup>2</sup></td></tr><tr><td>5</td><td>6.25 × 10<sup>-1</sup></td><td>1.66</td><td>9.22</td><td>7.03 × 10<sup>1</sup></td><td>6.92 × 10<sup>2</sup></td></tr><tr><td>6</td><td>6.88 × 10<sup>-1</sup></td><td>2.13</td><td>1.10 × 10<sup>1</sup></td><td>7.83 × 10<sup>1</sup></td><td>7.88 × 10<sup>2</sup></td></tr><tr><td>7</td><td>8.28 × 10<sup>-1</sup></td><td>2.63</td><td>1.30 × 10<sup>1</sup></td><td>1.06 × 10<sup>2</sup></td><td>8.21 × 10<sup>2</sup></td></tr><tr><td>8</td><td>1.23</td><td>3.42</td><td>1.48 × 10<sup>1</sup></td><td>1.09 × 10<sup>2</sup></td><td>8.60 × 10<sup>2</sup></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 7: Precision vs. running time (in seconds) for Space-Time methods using different Weighted Quadrature [22]. Graphic representation of Tables 1 and 2.</figcaption></figure>

In Fig. 7, the relative error is compared against computational time for different mesh sizes and polynomial degrees using Weighted Quadrature. As previously demonstrated, the computational cost increases linearly with the number of elements and is almost independent of the polynomial degree, except for low-degree polynomials. For example, suppose we aim for a target accuracy around $10^{-3}$. In that case, we can either choose a polynomial degree $p_s = p_t = 6$ on a mesh with $4^3$ elements or $p_s = p_t = 1$ on a mesh with $64^3$ elements. The former approximation is achieved in less than one second, while the latter takes approximately 500 seconds on our workstations. If the desired accuracy is $10^{-8}$, we could use a polynomial degree $p_s = p_t = 5$ with $32^3$ elements; while, using an incremental method (equivalent to $p_t = 1$), the computational time would be impractical in real life, possibly around $10^8$ seconds, due to the severe penalty imposed by the order of convergence of the temporal scheme. This leads to the most important conclusion of our work: IGA k-refinement, combining Matrix-Free and Weighted Quadrature methods, is superior in every aspect, in terms of accuracy and computational speed, to low-order FEM.
## 5.3 Computational cost analysis

We now present a numerical analysis of the performance of Weighted Quadrature and Matrix-Free approaches described in Sec. 3. As previously mentioned, we employ the LAPACK functions `dsygvd` and `zgges` for Eigen and Schur decomposition, respectively. Dense-dense matrix products are performed using the Fortran90 intrinsic `matmul` function, while sparse-dense matrix products are computed with a custom algorithm using the CSR format. To measure execution time, we use the `process_time()` function from Python's `time` module.

For the results in Figs. 8a and 8b, we consider a 4D space-time domain, i.e. $\Omega \times \mathcal{T} \in \mathbb{R}^3 \times \mathbb{R}$, discretized such as $p_s = p_t = p$ and $h_s = h_t = h$. Then, we report in Fig. 8a the performance of Alg. 5. According to the results, the degree $p$ has no direct impact on the computational time required by Fast Diagonalization method. Moreover, the timings look impressive since, in the finest discretization, a problem with more than $128^4$ ($\approx 270$ million) elements is solved in less than 400 s. It is important to remark that within the running time shown in Fig. 8a, the time that it takes to perform Eigen and Schur decomposition is included; however, as mentioned in Sec. 4.2, these types of decomposition are actually offline computations that in a general procedure, such as Alg. 1, are performed right after defining the mesh discretization.

<figure><img src="image_9.png" /><figcaption>(a) Running time (in seconds) for Alg. 5 for different types of discretization.</figcaption></figure>

<caption>(b) Running time (in seconds) for Algs. 3 and 4 for Gauss-Legendre and Weighted Quadrature.</caption>

<caption>Figure 8: Computational cost for Fast Diagonalization and Matrix-Free methods.</caption>

For a given coarse mesh-size $h_s = h_t = 1/16$, we compared in Fig. 8b the running time that takes to compute a matrix-vector product using Algs. 3 and 4 considering Gauss-Legendre quadrature (MF-GL) and Weighted Quadrature (MF-WQ). As stated in Sec. 3.4, the cost of Matrix-Free algorithms combined with WQ is $O(n_{cp}p)$ while with GL is $O(n_{cp}p^{d+2})$. These theoretical results are reflected in the experimental data where it is shown that MF-WQ algorithms are faster than MF-GL ones, except when $p \le 2$. The latter case is explained because Weighted Quadrature, as defined, uses more quadrature points than Gauss-Legendre quadrature when using lower degrees. In our experiments, the computer ran out of memory for degrees $p \ge 7$ when using Gauss-Legendre quadrature, and that is why those results are not present in Fig. 8b. Furthermore, it is important to remark that, regardless of the quadrature rule, the time ratio between the matrix-vector products $Av$ and $Bv$ in our experiments is almost 2, which is slightly better than the theoretical ratio for a 3D geometry. As a result, one iteration of the linear solver incurs approximately 50% more computational cost when using Newton's linearization compared to Picard's method. This observation is crucial when deciding the accuracy of the solver, whether an exact or inexact approach is employed, as these options differ significantly in the associated computational expense. Finally, it is worth noticing that, although theoretically one would expect better times than with the Fast Diagonalization algorithm, this is not true due to the high
performance of modern computers in computing products of dense matrices and optimizing memory usage. In this case, the time of applying the preconditioner using a coarse mesh-size $h_s = h_t = 1/16$ is less than a 0.1 s which is at least 10 times lower than the time performing the matrix-vector products.

## 5.4 Comparison of nonlinear solvers

For the same nonlinear problem already described in Sec. 5.1, we now compare the performance of nonlinear solvers described in Sec. 2.4: Newton and Picard methods combined with exact or inexact linear solvers. It is worth recalling that the parameters used in the stopping criterion, as defined in Eq. (15), are $\epsilon_{svg} = 10^{-12}$, $\epsilon_{tight}^{LK} = 10^{-8}$ and $\epsilon_{loose}^{LK} = 0.5$. For "exact" methods, the forcing term $\epsilon_{(k)}^{LK}$ is consistently set to $\epsilon_{tight}^{LK}$. For inexact Picard, we use a fixed loose threshold of $\epsilon_{(k)}^{LK} = 0.1$. For inexact Newton, $\epsilon_{(k)}^{LK}$ is dynamically computed based on Eq. (14), assuming that $\epsilon_{(0)}^{LK} = \epsilon_{loose}^{LK}$, $\gamma = 1$. and $\zeta = 1.5$. A preliminary parametric study, inspired by a similar analysis conducted by Eisenstat and Walker [35][Tables 3.1 and 3.2], was performed to assess the influence of these parameters, leading to the selected values. Further details on this study can be found in Appendix A. Moreover, for the following results we have considered Gauss-Legendre quadrature to prevent any loose of accuracy.

<figure><img src="image_10.png" /><figcaption>(a) Relative residue w.r.t. nonlinear iterations.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Relative residue w.r.t. linear iterations.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) Relative error wr.r.t. nonlinear iterations.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(d) Relative error wr.r.t. linear iterations.</figcaption></figure>

<caption>Figure 9: Comparison of nonlinear solvers for a quart-of-annulus geometry using degree $p = 3$ and mesh-size $h = 1/16$.</caption>
In Fig. 9a, we show the evolution of the relative norm of the nonlinear residue over the nonlinear iterations. As expected from a theoretical point of view, “exact” methods preserve optimal convergence rate, whereas “inexact” methods do not. To be precise, Inexact Picard recovers optimal linear convergence after 3 (nonlinear) iterations; while Inexact Newton shows first linear and then superlinear convergence, unfortunately not becoming quadratic. For this specific nonlinear problem, the choice of “inexact” methods has a moderate impact on the number of nonlinear iterations for the same order of magnitude of the nonlinear residue. In the case of using a direct solver for the linearized problem, the “inexact” methods would increase the computational cost compared to their respective “exact” counterparts. However, since linear iterative solvers are preferred, we have to take into account the total number of iterations in our analysis. In Fig. 9b, instead of the number of nonlinear iterations as in Fig. 9a, we present the nonlinear residue with respect to the total number of linear iterations accumulated. Here, the motivation behind the use of “inexact” methods can be seen in more detail. Both “inexact” Picard and Newton converge after about 20 linear iterations, which represents 20% and 50% of what would be obtained using their respective “exact” counterparts. These results support what is mentioned in Sec. 2.4: at least in the first nonlinear iterations, it is not necessary to compute a precise solution to the linearized problem.

To complete our analysis, we show in Figs. 9c and 9d the evolution of the error $u - u^h$ in the relative $L^2(\Pi)$-norm in function of the nonlinear and linear iterations, respectively. As shown, after some nonlinear iterations, the error no longer decreases even if the nonlinear residue does as in Figs. 9a and 9b. The reason for this interesting observation is that after a few iterations, in the sense of the $L^2$-norm, the increment $\delta u$ becomes negligible with respect to the error caused by the discretization. It shows that, in order to avoid oversolving in the nonlinear loop, the stopping criterion Eq. (11) based on the Euclidean norm is not well-adapted for nonlinear systems derived from the Galerkin method. Specifically, for such systems, a more appropriate choice would be the $L^2$- and $H^1$-norms, as applied by Kamdem and Sayah [51]. Alternatively, Arioli and Loghin [37] suggest that more complex norms, based on dual spaces of the function spaces used in the formulation, may provide more accurate and reliable stopping criteria for nonlinear problems. However, this is beyond the scope of this paper and will be left for our future work.

## 5.5 A more complex example

Consider that we study a domain $\Omega = \{x = (x_1, x_2) | 0 \le x_1 \le 0.4, 0.75x_1 - 0.5 \le x_2 \le -0.75x_1 + 0.5\}$ with material properties $\sigma = 1$ and

$$k(u) = (3 + 2 \tanh(u/50) + \cos(3\pi u/50))) \begin{bmatrix} 1 & 0.5 \\ 0.5 & 2 \end{bmatrix}$$

during the time interval $T = [0, 1]$. The heat source imposed is computed analytically such as the manufactured solution is $u(x, t) = \sin(6\pi x_2) \sin(6\pi(x_2 + 0.75x_1 - 0.5)(-x_2 + 0.75x_1 - 0.5)) \sin(10\pi x_1)g(t)$, where $g(t)$ is the same function as previously defined. We consider as well a zero temperature over the domain at $t = 0$ and homogeneous Dirichlet boundary conditions over $\Gamma \times T$. In Fig. 10, the temperature field and thermal conductivity across the geometry $\Omega$ are shown at the end of the time interval $T$.

<figure><img src="image_14.png" /><figcaption>Figure 10: Graphic representation of the exact temperature and the thermal conductivity at $t = 1$.</figcaption></figure>

Compared to Fig. 2, this case deals with a pathological problem that has greater heterogeneity and
stronger property gradients. As a result, a decrease in performance is expected, meaning that a minimum number of elements is required to achieve the asymptotic convergence rate. In Fig. 11 we present the evolution of the error, in the sense of the $L^2(\Pi)$-norm, in a set of successively refined meshes, with $h_s = h_t = h$, for different degrees varying from 1 to 5, with $p_s = p_t = p$. To get a clear idea of the convergence problem, we considered relevant to show the differences between the quadrature methods: standard Gauss-Legendre (ST-IGA-GL) in Fig. 11a and Weighted Quadrature (ST-IGA-WQ-2) in Fig. 11b. In both cases, we observe that coarse meshes fail to achieve acceptable precision, unlike the results in Fig. 3 for the first example. This outcome aligns with our expectations. Moreover, above a certain level of refinement, the results obtained using Gauss-Legendre quadrature show the optimal convergence rate. In contrast, Weighted Quadrature does not reach this rate. Instead, it experiences a significant loss of accuracy and probably deviates from the optimal convergence rate. This discrepancy is probably due to an insufficient number of quadrature points, inadequate to effectively capture the heterogeneity of the material.

<figure><img src="image_15.png" /><figcaption>Figure 11: $L^2(\Pi)$-norm of $u - u^h$ using Gauss-Legendre and Weighted Quadrature methods.</figcaption></figure>

Similar observations can be made when analyzing the performance of the nonlinear solvers discussed in Sec. 2.4. To ensure accuracy, we used Gauss-Legendre quadrature in the following results, avoiding potential loss of accuracy from Weighted Quadrature methods. In Fig. 12a, we present the evolution of the relative norm of the nonlinear residual over the nonlinear iterations. Unlike the previous example, Exact Newton behaves almost identically to Inexact Newton. Initially, Inexact Newton converges at the same rate as Inexact Picard because their forcing terms are similar. However, after the third iteration, Inexact Newton recovers its optimal convergence, matching the asymptotic behavior of Exact Newton. For this particular problem, using “inexact” methods does not affect the total number of Newton iterations needed to achieve a similar reduction in the nonlinear residual. However, for the Picard method, the exact version slightly reduces the number of nonlinear iterations. As mentioned above, to make a fair comparison, we must focus on the total number of matrix-vector products to achieve convergence. Therefore, in Fig. 12b, we examine the nonlinear residual relative to the total number of accumulated linear iterations, instead of the number of nonlinear iterations as shown in Fig. 12a. This highlights the advantage of “inexact” methods: Inexact Picard and Newton converge faster than their “exact” counterparts. Finally, Figs. 12c and 12d shows the evolution of the error against the nonlinear and linear iterations, respectively. After a certain number of nonlinear iterations, the error stops decreasing, even though the nonlinear residual continues to drop, as seen in Figs. 12a and 12b. This indicates that using a stopping criterion based solely on the Euclidean norm is not appropriate for such problems. A more suitable stopping criterion is needed.
<figure><img src="image_16.png" /><figcaption>(a) Relative residue w.r.t. nonlinear iterations.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Relative residue w.r.t. linear iterations.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(c) Relative error wr.r.t. nonlinear iterations.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(d) Relative error w.r.t. linear iterations.</figcaption></figure>

Figure 12: Comparison of nonlinear solvers for a trapezoidal geometry using degree $p = 3$ and mesh-size $h = 1/16$.

# 6 Conclusion

In this work, we have introduced the essential components required to solve nonlinear transient heat transfer problems using Space-Time IsoGeometric Analysis (ST-IGA) with high-order and high-continuous splines. Our primary goal has been to tackle the significant computational challenge associated with increasing the problem dimensionality. To optimize computational cost in terms of both CPU time and memory usage, we employed innovative methods such as the Matrix-Free Weighted Quadrature (MF-WQ) approach for efficient matrix-vector operations and the Fast Diagonalization Technique to accelerate the convergence of the linear solver. We show in our theoretical framework that the resulting preconditioner is straightforward to implement in an IGA code. In our numerical benchmark, we addressed a nonlinear parabolic problem with smooth solution. The results consistently show that Space-Time $k$-refinement, when combined with the MF-WQ approach, significantly outperforms traditional low-order methods, even if it means solving a problem of higher dimension. Additionally, we explored further computational improvements by incorporating adaptive resolution strategies, such as the “inexact” approach. This method avoids the oversolving phenomenon, enhancing efficiency during the early iterations of the nonlinear solver. The results show that implementing
this idea into our algorithms can reduce the final computational cost by up to 50%. Together, these techniques culminate in a robust and efficient framework able to solve complex heat transfer problems with reduced computational overhead, achieving the main objective of this study.

Further research will focus on the extension of these techniques to more complex problems such as welding simulations involving moving heat sources. Such scenarios may introduce numerical instabilities that the current formulation might not adequately address. Therefore, additional stabilization matrices might need to be implemented to ensure a stable solution. However, many stabilization techniques in the literature rely on strongly nonlinear matrices that necessitate explicit computation, potentially conflicting with our Matrix-Free approach. Investigating the impact of these matrices on the computational cost and overall performance of the algorithm is an interesting avenue for future exploration.

## Acknowledgement

This work was partially supported by the ANR METALIC project, grant ANR-21-CE10-0003; and by the joint laboratory for the numerical simulation of manufacturing and repair processes between LaMCoS and the companies EDF and Framatome. J. Cornejo-Fuentes was also supported by a doctoral scholarship from INSA-Lyon. These supports are gratefully acknowledged.

## A Impact of forcing term parameters on convergence

To investigate the influence of different parameters on the number of nonlinear and linear iterations, we implemented an experimental design based on an exponential distribution. This approach allowed us to efficiently explore parameter dependencies. We then performed an analysis of variance (ANOVA) to assess the significance of each parameter, finding that while all parameters affect the number of Newton's and Krylov's iterations, the two most influential ones are $\gamma$ and $\zeta$. Based on this finding, we selected these parameters for visualization in Figs. 13 and 14. Figs. 13a and 13b show respectively the number of nonlinear and linear iterations as a function of $\gamma$ and $\zeta$, with a the parameter $\epsilon_{(0)}^{LK}$ fixed at 1/2. Similarly, Figs. 14a and 14b present the same analysis but with the parameter $\epsilon_{(0)}^{LK}$ set to 1/32. It is important to note that the results presented here guided the selection of parameter values for the simulations described in Sec. 5.4.

<figure><img src="image_20.png" /><figcaption>Figure 13: Number of nonlinear and linear iterations for different values of $\gamma$ and $\zeta$ with a fixed value of $\epsilon_{(0)}^{LK} = 1/2$.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 14: Number of nonlinear and linear iterations for different values of $\gamma$ and $\zeta$ with a fixed value of $\epsilon_{(0)}^{LK} = 1/32$.</figcaption></figure>

Comparing these results, it is not easy to find a clear or systematic trend. However, we can draw some qualitative conclusions. In the Newton solver, by choosing parameters following the trend $\gamma \to 0$, $\zeta \to 2$, $\epsilon_{(0)}^{LK} \to 0$ (towards the upper-right corner of Figs. 14a and 14b), the forcing term becomes very tight, similar to the case of an exact Newton method. As a result, we recover the ideal convergence behavior, though the price to pay is oversolving the system. This explains why the number of linear iterations increases (in some cases, a factor of 2 can even be observed between the most favorable and least favorable cases).

For the analysis of the number of iterations in the Krylov solver, another phenomenon must also be considered: in the extreme case where the evolution of the forcing term is very slow, i.e., when $\gamma \to 1$, $\zeta \to 1$, $\epsilon_{(0)}^{LK} \to 1/2$ (towards the lower-left corner of Figs. 13a and 13b). In this case, the Krylov solver performs very few iterations (potentially less than 5 per nonlinear iteration), which leads to insufficient progress towards the solution. This results in an increase in the number of Newton iterations, losing the superlinear (or quadratic) convergence, and slightly increasing the number of linear iterations.

It is worth mentioning that for values of $\gamma$, $\zeta$, $\epsilon_{(0)}^{LK}$ that are not within the 'extreme' cases, their influence on the number of iterations (especially the linear iterations, which are of particular interest due to their computational cost) is not entirely clear. It is important to note that this analysis depends on the problem, as different parameter triplets may be optimal for different configurations of materials, geometries, mesh discretizations, etc. Overall, our results suggest a complex interplay between these parameters, without an immediately identifiable pattern. A more detailed analysis would be required to fully characterize these relationships. However, such an investigation falls beyond the scope of this work and is left for future research.

## References

[1] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs, “Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement,” *Computer Methods in Applied Mechanics and Engineering*, vol. 194, no. 39, pp. 4135–4195, 2005. DOI: 10.1016/j.cma.2004.10.008.
[2] A. Bressan and E. Sande, "Approximation in FEM, DG and IGA: A theoretical comparison," *Numerische Mathematik*, vol. 143, no. 4, pp. 923-942, 2019. DOI: 10.1007/s00211-019-01063-5.

[3] U. Langer and O. Steinbach, *Space-Time Methods: Applications to Partial Differential Equations*. Berlin, Boston: De Gruyter, 2019. DOI: 10.1515/9783110548488.

[4] Z. Zhao and H. Li, "Numerical analysis of a space-time continuous Galerkin method for the nonlinear Sobolev equation," *Journal of Mathematical Analysis and Applications*, vol. 531, no. 2, Part 2, p. 127-145, 2024. DOI: 10.1016/j.jmaa.2023.127845.

[5] Y. Maday and E. M. Rønquist, "Parallelization in time through tensor-product space-time solvers," *Comptes Rendus Mathematique*, vol. 346, no. 1, pp. 113-118, 2008. DOI: 10.1016/j.crma.2007.09.012.

[6] I. Fried, "Finite-element analysis of time-dependent phenomena," *AIAA Journal*, vol. 7, no. 6, pp. 1170-1173, 1969. DOI: 10.2514/3.5299.

[7] J. T. Oden, "A general theory of finite elements. I. Topological considerations," *International Journal for Numerical Methods in Engineering*, vol. 1, no. 2, pp. 205-221, 1969. DOI: 10.1002/nme.1620010209.

[8] J. T. Oden, "A general theory of finite elements. II. Applications," *International Journal for Numerical Methods in Engineering*, vol. 1, no. 3, pp. 247-259, 1969. DOI: 10.1002/nme.1620010304.

[9] J. C. Bruch Jr and G. Zyvoloski, "Transient two-dimensional heat conduction problems solved by the finite element method," *International Journal for Numerical Methods in Engineering*, vol. 8, no. 3, pp. 481-494, 1974. DOI: 10.1002/nme.1620080304.

[10] C. Schwab and R. Stevenson, "Space-time adaptive wavelet methods for parabolic evolution problems," *Mathematics of Computation*, vol. 78, no. 267, pp. 1293-1318, 2009. DOI: 10.1090/S0025-5718-08-02205-9.

[11] O. Steinbach, "Space-time finite element methods for parabolic problems," *Computational Methods in Applied Mathematics*, vol. 15, no. 4, pp. 551-566, 2015. DOI: doi:10.1515/cmam-2015-0026.

[12] U. Langer, S. E. Moore, and M. Neumüller, "Space-time isogeometric analysis of parabolic evolution problems," *Computer Methods in Applied Mechanics and Engineering*, vol. 306, pp. 342-363, 2016. DOI: 10.1016/j.cma.2016.03.042.

[13] K. Takizawa, T. E. Tezduyar, Y. Otoguro, T. Terahara, T. Kuraishi, and H. Hattori, "Turbocharger flow computations with the space-time isogeometric analysis (ST-IGA)," *Computers & Fluids*, vol. 142, pp. 15-20, 2017. DOI: 10.1016/j.compfluid.2016.02.021.

[14] K. Takizawa, T. E. Tezduyar, and T. Terahara, "Ram-air parachute structural and fluid mechanics computations with the space-time isogeometric analysis (ST-IGA)," *Computers & Fluids*, vol. 141, pp. 191-200, 2016, Advances in Fluid-Structure Interaction. DOI: 10.1016/j.compfluid.2016.05.027.

[15] C. Saadé, S. Lejeunes, D. Eyheramendy, and R. Saad, "Space-time isogeometric analysis for linear and non-linear elastodynamics," *Computers & Structures*, vol. 254, p. 106-109, 2021. DOI: 10.1016/j.compstruc.2021.106594.

[16] F. K. Feutang, S. Lejeunes, and D. Eyheramendy, "A space-time approach for the simulation of brittle fracture with phase-field models in elastodynamics," *Computers & Structures*, vol. 307, p. 107-116, 2025. DOI: 10.1016/j.compstruc.2024.107616.

[17] M. Montardini, M. Negri, G. Sangalli, and M. Tani, "Space-time least-squares isogeometric method and efficient solver for parabolic problems," *Mathematics of Computation*, vol. 89, no. 323, pp. 1193-1227, 2020. DOI: 10.1090/mcom/3471.
[18] G. Loli, M. Montardini, G. Sangalli, and M. Tani, "An efficient solver for space-time isogeometric galerkin methods for parabolic problems," *Computers & Mathematics with Applications*, vol. 80, no. 11, pp. 2586-2603, 2020. DOI: 10.1016/j.camwa.2020.09.014.

[19] J. N. Reddy, *An Introduction to Nonlinear Finite Element Analysis*. Oxford University Press, Mar. 2004. DOI: 10.1093/acprof:oso/9780198525295.001.0001.

[20] G. Loli, G. Sangalli, and P. Tesini, "High-order spline upwind for space-time isogeometric analysis," *Computer Methods in Applied Mechanics and Engineering*, vol. 417, p. 116408, 2023. DOI: 10.1016/j.cma.2023.116408.

[21] F. Calabrò, G. Sangalli, and M. Tani, "Fast formation of isogeometric Galerkin matrices by weighted quadrature," *Computer Methods in Applied Mechanics and Engineering*, vol. 316, pp. 606-622, 2017. DOI: 10.1016/j.cma.2016.09.013.

[22] R. Hiemstra, G. Sangalli, M. Tani, F. Calabrò, and T. J. R. Hughes, "Fast formation and assembly of finite element matrices with application to isogeometric linear elasticity," *Computer Methods in Applied Mechanics and Engineering*, vol. 355, pp. 234-260, 2019. DOI: 10.1016/j.cma.2019.06.020.

[23] G. Sangalli and M. Tani, "Matrix-free weighted quadrature for a computationally efficient isogeometric k-method," *Computer Methods in Applied Mechanics and Engineering*, vol. 338, pp. 117-133, 2018. DOI: 10.1016/j.cma.2018.04.029.

[24] P. Antolin, A. Buffa, F. Calabrò, M. Martinelli, and G. Sangalli, "Efficient matrix computation for tensor-product isogeometric analysis: The use of sum factorization," *Computer Methods in Applied Mechanics and Engineering*, vol. 285, pp. 817-828, 2015. DOI: 10.1016/j.cma.2014.12.013.

[25] N. Collier, L. Dalcin, D. Pardo, and V. M. Calo, "The cost of continuity: Performance of iterative solvers on isogeometric finite elements," *SIAM Journal on Scientific Computing*, vol. 35, no. 2, A767-A784, 2013. DOI: 10.1137/120881038.

[26] G. Sangalli and M. Tani, "Isogeometric preconditioners based on fast solvers for the Sylvester equation," *SIAM Journal on Scientific Computing*, vol. 38, no. 6, A3644-A3671, 2016. DOI: 10.1137/16M1062788.

[27] R. E. Lynch, J. R. Rice, and D. H. Thomas, "Tensor product analysis of partial difference equations," *Bulletin of the American Mathematical Society*, vol. 70, no. 3, pp. 378-384, 1964.

[28] A. Duval, T. Hirschler, J. Cornejo Fuentes, M. Guerder, and T. Elguedj, *YETI: Yet another iga code*, https://lamcosplm.insa-lyon.fr/projects/yeti, 2023.

[29] J. A. Cottrell, T. J. R. Hughes, and Y. Bazilevs, *Isogeometric Analysis: Toward Integration of CAD and FEA*. John Wiley & Sons, 2009. DOI: 10.1002/9780470749081.

[30] J. M. Ortega and W. C. Rheinboldt, *Iterative Solution of Nonlinear Equations in Several Variables*. Academic Press, 1970. DOI: 10.1137/1.9780898719468.

[31] R. Barrett, M. Berry, T. F. Chan, et al., *Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods*. SIAM, 1994. DOI: 10.1137/1.9781611971538.

[32] R. S. Dembo, S. C. Eisenstat, and T. Steihaug, "Inexact Newton methods," *SIAM Journal on Numerical Analysis*, vol. 19, no. 2, pp. 400-408, 1982. DOI: 10.1137/0719025.

[33] P. Birken, "Termination criteria for inexact fixed-point schemes," *Numerical Linear Algebra with Applications*, vol. 22, no. 4, pp. 702-716, 2015. DOI: 10.1002/nla.1982.

[34] H. Elman, D. Silvester, and A. Wathen, *Finite Elements and Fast Iterative Solvers: with Applications in Incompressible Fluid Dynamics*. Oxford University Press, Jun. 2014. DOI: 10.1093/acprof:oso/9780199678792.001.0001.

[35] S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an inexact Newton method," *SIAM Journal on Scientific Computing*, vol. 17, no. 1, pp. 16-32, 1996. DOI: 10.1137/0917003.
[36] M. Porcelli, V. Simoncini, and M. Tani, "Preconditioning of active-set Newton methods for PDE-constrained optimal control problems," *SIAM Journal on Scientific Computing*, vol. 37, no. 5, S472-S502, 2015. DOI: 10.1137/140975711.

[37] M. Arioli and D. Loghin, "Stopping criteria for mixed finite element problems," *Electronic Transactions on Numerical Analysis*, vol. 29, pp. 178-192, 2008.

[38] H.-B. An, Z.-Y. Mo, and X.-P. Liu, "A choice of forcing terms in inexact Newton method," *Journal of Computational and Applied Mathematics*, vol. 200, no. 1, pp. 47-60, 2007. DOI: 10.1016/j.cam.2005.12.030.

[39] N. Tardieu and E. Cheignon, "A Newton-Krylov method for solid mechanics," *European Journal of Computational Mechanics*, vol. 21, no. 3-6, pp. 374-384, 2012. DOI: 10.1080/17797179.2012.721501.

[40] T. G. Kolda and B. W. Bader, "Tensor decompositions and applications," *SIAM Review*, vol. 51, no. 3, pp. 455-500, 2009. DOI: 10.1137/07070111X.

[41] M. Bartoň and V. M. Calo, "Gauss-Galerkin quadrature rules for quadratic and cubic spline spaces and their application to isogeometric analysis," *Computer-Aided Design*, vol. 82, pp. 57-67, 2017. DOI: 10.1016/j.cad.2016.07.003.

[42] Z. Zou, T. J. R. Hughes, M. A. Scott, D. Miao, and R. A. Sauer, "Efficient and robust quadratures for isogeometric analysis: Reduced Gauss and Gauss-Greville rules," *Computer Methods in Applied Mechanics and Engineering*, vol. 392, p. 114722, 2022. DOI: 10.1016/j.cma.2022.114722.

[43] G. H. Golub and J. H. Welsch, "Calculation of Gauss quadrature rules," *Mathematics of computation*, vol. 23, no. 106, pp. 221-230, 1969. DOI: 10.2307/2004418.

[44] P. W. Smith, "A practical guide to splines (Carl de Boor)," *SIAM Review*, vol. 22, no. 4, p. 520, 1980. DOI: 10.1137/1022106.

[45] M. Montardini, G. Sangalli, and M. Tani, "Robust isogeometric preconditioners for the Stokes system based on the fast diagonalization method," *Computer Methods in Applied Mechanics and Engineering*, vol. 338, pp. 162-185, 2018. DOI: 10.1016/j.cma.2018.04.017.

[46] P. D. Brubeck and P. E. Farrell, "A scalable and robust vertex-star relaxation for high-order FEM," *SIAM Journal on Scientific Computing*, vol. 44, no. 5, A2991-A3017, 2022. DOI: 10.1137/21M1444187.

[47] G. H. Golub and C. F. van Loan, *Matrix Computations* (Johns Hopkins Studies in the Mathematical Sciences). Johns Hopkins University Press, 2013.

[48] E. Anderson, Z. Bai, C. Bischof, et al., *LAPACK Users' Guide*. Society for Industrial and Applied Mathematics, 1999.

[49] J. M. Winget and T. J. R. Hughes, "Solution algorithms for nonlinear transient heat conduction analysis employing element-by-element iterative strategies," *Computer Methods in Applied Mechanics and Engineering*, vol. 52, no. 1, pp. 711-815, 1985. DOI: 10.1016/0045-7825(85)90015-5.

[50] J. Cornejo Fuentes, T. Elguedj, D. Dureisseix, and A. Duval, "A cheap preconditioner based on fast diagonalization method for matrix-free weighted-quadrature isogeometric analysis applied to nonlinear transient heat transfer problems," *Computer Methods in Applied Mechanics and Engineering*, vol. 414, p. 116157, 2023. DOI: 10.1016/j.cma.2023.116157.

[51] J. D. Kamdem and T. Sayah, "Stopping criteria for nonlinear variational problems: Iterative approach," *Numerical Algorithms*, Nov. 2024. DOI: 10.1007/s11075-024-01969-8.