\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRMECA,Unicode,published]{cedram}
\newenvironment{noXML}{}{}%
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\sfrac#1#2{{#1}/{#2}}
\def\sofrac#1#2{({#1}/{#2})}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}\hline}
\def\tabnote#1{\vskip4pt\parbox{.94\linewidth}{#1}}
\def\botline{\\\hline}
\def\mathbi#1{\boldsymbol{#1}}
\def\bnabla{\bm{\nabla}}
\def\sfrac#1#2{{#1}/{#2}}
\def\stfrac#1#2{({#1})/{#2}}
\usepackage{etoolbox}
\usepackage{bm}
\usepackage{upgreek}
\newcounter{runlevel}
\let\MakeYrStrItalic\relax
\def\jobid{crmeca20240874}
\def\morerows#1#2{#2}
%\graphicspath{{/tmp/\jobid_figs/web/}} 
\graphicspath{{./figures/}} 
\def\refinput#1{}
\def\back#1{}
\csdef{Seqnsplit}{\\}
\DOI{10.5802/crmeca.281}
\datereceived{2024-10-23}
\daterevised{2024-12-18}
\dateaccepted{2024-12-18}
\ItHasTeXPublished
}
{\documentclass[crmeca]{article}
\usepackage{upgreek}
\usepackage[T1]{fontenc} 
\def\CDRdoi{10.5802/crmeca.281}
\def\0{\phantom{0}}
\def\newline{\break}
\def\sfrac#1#2{{#1}/{#2}}
\def\stfrac#1#2{({#1})/{#2}}
}
\makeatother

\usepackage{upgreek}

\dateposted{2025-01-09}
\begin{document}

\begin{noXML}

\CDRsetmeta{articletype}{research-article}

\title{An accurate and robust FFT-based solver for transient diffusion
in heterogeneous materials}

\alttitle{Un solveur FFT pr\'{e}cis et robuste pour la diffusion
transitoire dans les mat\'{e}riaux h\'{e}t\'{e}rog\`{e}nes}

\author{\firstname{Lionel} \lastname{G\'{e}l\'{e}bart}\CDRorcid{0000-0002-1387-1978}}
\address{Universit\'{e} Paris Saclay, CEA, SRMA, 91191, Gif/Yvette, France}
\email{lionel.gelebart@cea.fr}

\begin{abstract}
The purpose of the present letter is to propose an efficient, accurate
and robust FFT-based solver for transient diffusion in heterogeneous
materials with ``realistic'' BC, taking advantage of two recent
advances in terms of boundary conditions and finite difference schemes
to overcome their actual limitations (periodic BC and spurious
oscillations). It is an essential step towards couplings between
mechanics and other physics (such as the diffusion of species) through
FFT-based solvers. Discrete Trigonometric Transform are used to
implement non-periodic boundary conditions, and a finite difference
(FD) scheme recently proposed by Finel is advantageously compared to
the common hexahedral FD scheme. ``Accurate'' refers to two
properties: accurate in term of locality with a small size of Finite
Difference pencil to capture fluctuations around material
heterogeneities, and accurate in term of precision with the absence of
spurious spatial oscillations (at least in the reported cases with
well-separated inclusions). The ``robustness'' is here associated to
the stability of the solver, especially associated to the implicit time
integration method. The description of the method focuses on thermal
diffusion but applies to any other similar diffusion process (with the
same type of parabolic equation). As a by-product, the FD scheme
proposed by Finel is introduced in a more general (for mixing finite
different schemes) and simple way (no introduction of FCC subgrids),
extending its usage to any type of grid parity (not only even grids).
\end{abstract}

\begin{altabstract}
Cette note propose un solveur FFT efficace, pr\'{e}cis et robuste pour
la diffusion transitoire dans des mat\'{e}riaux h\'{e}t\'{e}rog\`{e}nes
avec des conditions aux limites (CL) dites
\og r\'{e}alistes\fg{}. Cette approche
exploite deux avanc\'{e}es r\'{e}centes : le traitement des conditions
aux limites non p\'{e}riodiques et l'utilisation d'un nouveau
sch\'{e}ma de diff\'{e}rences finies limitant les oscillations
parasites. Il s'agit d'une \'{e}tape essentielle vers le couplage entre
la m\'{e}canique et d'autres ph\'{e}nom\`{e}nes physiques (comme la
diffusion d'esp\`{e}ces) \`{a} l'aide de solveurs FFT. Les
Transform\'{e}es Trigonom\'{e}triques Discr\`{e}tes sont utilis\'{e}es
pour impl\'{e}menter des conditions aux limites non p\'{e}riodiques,
tandis qu'un sch\'{e}ma de diff\'{e}rences finies (DF) r\'{e}cemment
propos\'{e} par Finel est compar\'{e} avantageusement au sch\'{e}ma DF
classique hexa\'{e}drique.

Le terme \og pr\'{e}cis \fg{} se
r\'{e}f\`{e}re \`{a} deux aspects : (1) la pr\'{e}cision locale,
gr\^{a}ce \`{a} un sch\'{e}ma de DF de petite taille permettant de
capturer les fluctuations autour des h\'{e}t\'{e}rog\'{e}n\'{e}it\'{e}s
du mat\'{e}riau, et (2) la pr\'{e}cision en termes de fid\'{e}lit\'{e},
avec l'absence d'oscillations spatiales parasites (au moins dans les
cas \'{e}tudi\'{e}s avec des inclusions bien s\'{e}par\'{e}es). La 
\og robustesse \fg{} fait r\'{e}f\'{e}rence
\`{a} la stabilit\'{e} du solveur, notamment gr\^{a}ce une m\'{e}thode
d'int\'{e}gration temporelle implicite. Bien que la description se
concentre sur la diffusion thermique, la m\'{e}thode s'applique \`{a}
tout processus de diffusion similaire (r\'{e}gie par une \'{e}quation
parabolique du m\^{e}me type).

En compl\'{e}ment, le sch\'{e}ma DF propos\'{e} par Finel est
pr\'{e}sent\'{e} de mani\`{e}re plus g\'{e}n\'{e}rale (pour
m\'{e}langer diff\'{e}rents sch\'{e}mas de diff\'{e}rences finies) et
simplifi\'{e}e (sans introduction de sous-grilles FCC), \'{e}largissant
\'{e}galement son usage \`{a} tout type de parit\'{e} de grille (et pas
seulement aux grilles paires).
\clearpage
\end{altabstract}

\keywords{\kwd{FFT}
\kwd{Heterogeneous materials}
\kwd{Transient diffusion} 
\kwd{Finite difference}
\kwd{Tetrahedra}}


\altkeywords{\kwd{FFT}
\kwd{Mat\'{e}riaux h\'{e}t\'{e}rog\`{e}nes}
\kwd{Diffusion transitoire}
\kwd{Diff\'{e}rences finies}
\kwd{T\'{e}tra\`{e}dres}}


\maketitle

\end{noXML}

\section{Introduction}

Among many other applications, the development of an FFT-based solver
for transient diffusion in heterogeneous material is essential to study
the diffusion of species coming from the surface of the material and
their coupling with mechanics. For such coupled simulations, the
concentration of atoms induces a deformation of the crystal taken into
account through an eigenstrain (stress-free strain) directly
proportional to the concentration of atoms \cite{11,12}. Therefore, an
accurate description of the concentration field is required. Besides,
the boundary conditions generally used in this context are Dirichlet
boundary conditions (applied concentration at the surface of the
material, e.g.\ \cite{22,13}) or Neumann BC (applied flux, e.g.\ 
\cite{10}). Hence, for such simulations, the finite element method is
an appropriate numerical tool (e.g.\ \cite{22,13,10}). However, the
method is quite heavy in terms of memory or computation time,
especially when considering 3D complex microstructures. FFT-based
solvers for heterogeneous materials, initially proposed in the context
of static mechanics by Moulinec and Suquet \cite{14}, provide an
efficient alternative to standard FE solvers. The method introduces a
homogeneous material on which the solution of an auxiliary problem is
directly obtained by using Discrete Fourier Transforms (named FFT, for
Fast Fourier Transform, in the following) through the application of
the discrete Green operator. This efficient resolution of the auxiliary
problem is used as a step of a fixed-point algorithm \cite{14}.
Described by similar partial differential equations, the method applies
similarly to stationary diffusion problems on heterogeneous materials.
Besides their efficiency, these solvers suffer from two important
drawbacks in the present context. First, using Discrete Fourier
Transforms, the boundary conditions (BC) are implicitly periodic. To
circumvent this limitation, various teams have recently investigated
the question for diffusion \cite{9,15,17} and mechanics \cite{7,18,16},
taking advantage of Discrete Trigonometric Transforms (DTT), using both
Finite Difference or Galerkin based spatial approximations. Second, the
original method \cite{14} suffers from spurious oscillations attributed
to the Gibbs phenomenon. Using discrete Green operators built on Finite
Difference schemes drastically reduces these artefacts. Among them, the
most popular scheme, initially proposed by Willot and referred as the
rotated scheme \cite{23}, is equivalent to the use of hexaedral linear
finite elements with reduced integration (with a single Gauss point),
noted HEX8R \cite{19}. However, if the HEX8R is insensitive to Gibbs
phenomenon, it comes with another kind of oscillations known as
hourglass phenomenon \cite{2}. Various techniques propose to stabilize
hourglass \cite{2}, even in the context of FFT-based solvers \cite{21},
but they introduce an additional numerical parameter to be adjusted.
Very recently, Finel proposed another type of Finite Difference scheme
and its FFT-based implementation for static mechanics, completely
removing the spurious oscillations observed with HEX8R \cite{5},
without any additional numerical\unskip\break parameter.


The purpose of the present letter is to propose an efficient, accurate
and robust FFT-based solver for transient diffusion in heterogeneous
materials with ``realistic'' BC, taking advantage of two recent
advances in terms of boundary conditions and finite difference schemes
to overcome actual limitations (periodic BC and spurious oscillations).
``Accurate'' refers to two properties: accurate in term of locality
with small sizes of Finite Difference pencils to capture fluctuations
around material heterogeneities, and accurate in term of precision with
the absence of spurious spatial oscillations. The ``robustness'' is
here associated to the stability of the solver, especially associated
to the implicit time integration method.\vadjust{\pagebreak} The description of the method
focuses on thermal diffusion but applies to any other similar diffusion
process (with the same type of parabolic equation). As a by-product,
the FD scheme proposed by Finel \cite{5} is
introduced in a more general (for mixing finite different schemes) and
simple way (no introduction of FCC subgrids), extending its usage to
any type of grid parity (not only even grids).


\section{Method}\label{sec2}

\subsection{Implicit FFT-based resolution with non-periodic boundary conditions}\label{sec2.1}

The volume equations governing the transient thermal diffusion problem reads:
{\begin{eqnarray}\label{eqn1}
\begin{cases} 
\displaystyle
\rho c_{p}\frac{\partial T}{\partial t}+\mathrm{div}(\mathbi{q})=0
\vspace*{4pt}\\
\displaystyle
\mathbi{q}=-\underline{K}\cdot\bnabla \mathbi{T} 
\end{cases} 
\end{eqnarray}}\unskip
with $\rho$, $c_{p}$ and $\underline{K}$, the density, specific
heat capacity and thermal conductivity; material parameters that can be
heterogeneous over the domain $\Omega$, a parallelepiped domain of
size $L\times L\times L$~$\mathrm{m}^{3}$. $T$ and
$\mathbi{q}$ {are respectively the temperature and the heat flux. A
set of ``realistic'' initial and boundary conditions is proposed. The
initial temperature} $T_{0}$ is assumed homogeneous except on the
face $x=0$, where $T(x=0)=T^{\mathrm{left}}$. The boundary
conditions are Dirichlet BC on faces $x=0$ and $x=L$ with
respectively $T(x=0)=T^{\mathrm{left}}$ and
$T(x=L)=T_{0}$, for any time. The four other faces are
isolated (i.e.\ Neumann BC with $\mathbi{q}\cdot\mathbi{n}=0$,
$\mathbi{n}$ being the outer normal vector to the face).


Equation (\ref{eqn2}), with $\mathrm{d}T$ as unknown field, discretizes Equation
(\ref{eqn1}) in space and time, with a $\theta$-method in time, purely
implicit when $\theta =1$ (Euler implicit) and purely explicit when
$\theta =0$ (Euler explicit). The case $\theta =0.5$ corresponds to
the implicit Cranck--Nicholson method. The method is implicit as soon
as $\theta > 0$. The equation reads:
{\begin{eqnarray}\label{eqn2}
\rho c_{p}\frac{\mathrm{d}T}{\mathrm{d}t}-\mathrm{div}_{\mathrm{D}}
(\underline{K}\cdot\bnabla _{\mathbf{D}}(T^{t}+\theta \,\mathrm{d}T))=0,
\end{eqnarray}}\unskip
with $T^{t}$ the temperature field, known at the beginning of the time
step $(t)$, $T^{t}+\mathrm{d}T$, the temperature at the end of the time
step ($t+\mathrm{d}t$), and $\mathrm{d}T$ the unknown field. The
subscript ($\mathrm{D}$)
corresponds to a finite differences (FD) approximation for operators
$\mathrm{div}$ and $\bnabla$. The domain $\Omega $ is regularly
discretized in small elements (or voxels). For the two FD schemes
considered below, the temperature is defined at voxels' corners, as
well as $\mathrm{div}_{\mathrm{D}}(\mathbi{q})$, and the heat flux
$\mathbi{q}$ is defined at voxels' centres, as $\bnabla
_{\mathbf{D}}\mathbi{T}$. The grid of centres has a size $n^{3}$
and the grid of corners has a size $(n+1)^{3}$.



Introducing homogeneous quantities $(\rho c_{p})_{0}$ and
$k_{0}$ allows for defining a fix-point algorithm on $\mathrm{d}T$:
{\begin{eqnarray}\label{eqn3}
\frac{(\rho c_{p})_{0}}{\mathrm{d}t}\mathrm{d}T-\theta {k}_{0}\Delta 
_{\mathrm{D}}\,\mathrm{d}T=
\frac{(\rho c_{p})_{0}-\rho c_{p}}{\mathrm{d}t}\mathrm{d}T+\mathrm{div}_{\mathrm{D}}
(\underline{K}\cdot\bnabla _{\mathbf{D}}(T^{t}+\theta \,\mathrm{d}T)-\theta k_{0}
\bnabla _{\mathbf{D}}(\mathrm{d}T))
\end{eqnarray}}\unskip

Actually, for a given field $\mathrm{d}T$, the right hand-side (rhs) term
can be evaluated in real space using the FD scheme ($\mathrm{D}$).
Noting $p$ the rhs term, and assuming periodic BC for $\mathrm{d}T$
(non-periodic BC are considered later), the left hand-side (lhs)
term is inverted using Discrete Fourier Transforms to obtain
$\mathrm{d}T$  as follows:
{\begin{eqnarray}\label{eqn4}
\widehat{\mathrm{d}T}=\frac{\hat{p}}{\frac{(\rho c_{p})_{0}}{\mathrm{d}t}+
\theta k_{0}\Vert\boldsymbol{\xi }_{\mathrm{D}}\Vert^{2}}=
\widehat{G_{0}}(\boldsymbol{\xi })\hat{p}
\end{eqnarray}}\unskip
where $\boldsymbol{\xi }_{\mathrm{D}}$ are ``modified''
wave-vectors that are functions of the classical wave-vectors
$\boldsymbol{\xi}$. The expression of $\boldsymbol{\xi
}_{\mathrm{D}}$ depends on the chosen FD scheme ($\mathrm{D}$)
(various expressions can be found for different FD schemes in~\cite{20}).

In the case of non-periodic boundary conditions, expression
(\ref{eqn4}) still holds but with Discrete Fourier Transform performed
on 2 (or 4) times extended signals, evaluated by Discrete Trigonometric
Transforms (DTTs). The correspondence between the DTTs and the types of
BC are reported in Appendix~\ref{appA}. For a detailed description,
readers are invited to consult various recent papers considering the
extension of FFT-based solvers to non-periodic BC \cite{9,15,17} and
especially \cite{9} written for various FD schemes, including the HEX8R
scheme used below, and for any type of BC.


Regarding the BC applied in the present case, $\mathrm{d}T$ is assigned to null
Dirichlet BC (on faces $x=0$ and $x=L$), and null Neumann BC on the
other faces. As explained in \cite{9}, the null Dirichlet BC on faces
$x=0$ and $x=L$ are enforced during derivation in real space by
assuming an anti-symmetric field $\mathrm{d}T$ for the points located outside
the domain (i.e.~$x< 0$ and $x> L${). The equivalent in Fourier space
is achieved by using the Discrete Sine Transform of type I (DST1). For
the four other faces, the null flux is enforced during derivation by
assuming a symmetric field} $\mathrm{d}T$ {on the points located outside the
domain (i.e.~$y< 0$ and $y> L$, $z< 0$ and $z> L$). The equivalent in
Fourier space is achieved using the Discrete Cosine Transform of type I
(DCT1). Note that for the two FD schemes considered below (HEX8R and
TETRA2), the gradient terms $\bnabla _{\mathbf{D}}(\cdot)$ of 
Equation~(\ref{eqn3}) are evaluated at voxels' centres from the temperature
field defined at voxels' nodes. Hence, their evaluation with these FD
schemes do not rely on quantities located outside the domain. On the
contrary, the divergence term $\mathrm{div}_{\mathrm{D}}(\cdot)$ of
Equation (\ref{eqn3}) is evaluated (with HEX8R or TETRA2 scheme) at
voxels' nodes from the quantity ($\underline{K}\cdot\bnabla
_{\mathbf{D}}(T^{t}+\theta \,\mathrm{d}T)-\theta k_{0}\bnabla 
_{\mathbf{D}}(\mathrm{d}T)$)
defined at voxel's centres. Hence, the evaluation of the
$\mathrm{div}_{\mathrm{D}}(\cdot)$ term relies on a quantity defined
outside the domain using appropriate symmetries. For the quantity
$\underline{K}\cdot\bnabla _{\mathbf{D}}(T^{t}+\theta \,\mathrm{d}T)-\theta
k_{0}\bnabla _{\mathbf{D}}(\mathrm{d}T)$, the symmetries are those of $\bnabla
_{\mathbf{D}}(\mathrm{d}T)$ (for each component of the gradient, the derivation
in a given direction inverses the type of symmetry in that direction
(symmetry $\leftrightarrow$ anti-symmetry) and let unchanged the type of symmetry
in the two other directions, see~\cite{9}). 


\subsection{Focus on two finite difference schemes}

The presentation above introduces the $\theta$-method in time without
any precision on the FD scheme used in space. Focus is put on small
size Finite Difference pencils: the widely adopted HEX8R scheme
\cite{23,19} and the recently proposed tetrahedral-based scheme
\cite{5}, denoted TETRA2 in the following (the ``2'' standing for the
two tetrahedras used to define the scheme).


For the HEX8R scheme, the gradient is evaluated at a voxel's centre
from the 8 values at voxel's corners and the divergence at a
voxel's corner from the 8 values defined at the 8 neighbouring
voxel's centres. In both cases, the finite difference derivation can
be deduced from the approximated volume averages given below (where
$\Omega_{\mathrm{D}}=\Omega _{\mathrm{H}}$ is an hexahedra depicted on
Figure~\ref{fig1}, whose vertices are whether voxels'
corners, for the temperature gradient, or the voxels' centres for the
flux divergence): 
{\begin{eqnarray}\label{eqn5}
\begin{cases} 
\displaystyle
\bnabla \mathbi{T}_{\mathrm{D}}=\frac{1}{V_{\mathrm{D}}}
\iiint_{\Omega _{\mathrm{D}}}\bnabla \mathbi{T}\,\mathrm{d}v=\frac{1}{V_{\mathrm{D}}}
\iint _{\partial {\Omega _{\mathrm{D}}}}T\mathbi{n}\,\mathrm{d}s 
\vspace*{4pt}\\
\displaystyle
\mathrm{div}_{\mathrm{D}}(\mathbi{q})=\frac{1}{V_{\mathrm{D}}}
\iiint_{\Omega _{\mathrm{D}}}\mathrm{div}(\mathbi{q})\,\mathrm{d}v=
\frac{1}{V_{\mathrm{D}}}\iint _{\partial {\Omega _{\mathrm{D}}}}\mathbi{q}\cdot\mathbi{n}
\,\mathrm{d}s
\end{cases} 
\end{eqnarray}}\unskip


\begin{figure}
\vspace*{-3pt}
\includegraphics{fig01}
\vspace*{-4pt}
\caption{\label{fig1}Definition of the volumes used for the different Finite
Difference schemes (see Equation~(\ref{eqn5})).}
\vspace*{-3pt}
\end{figure}

The approximation of integrals, on each face of $\partial \Omega
_{\mathrm{D}}$, by a discrete sum on nodal quantities allows defining the
complete FD scheme. Note that in that case, the same finite difference
scheme is used for gradient and divergence, the only difference comes
from the location of the derived quantity (whether a grid of corners,
for the temperature, or a grid of centres, for the flux). The usage of
the HEX8R scheme in FFT-based methods is well documented and the
expression of the associated modified wave-vectors, necessary for the
FFT-based implementation, can be found in~\cite{23,20}. 

The TETRA2 scheme corresponds to a very recent proposition by Finel
\cite{5}, in the context of numerical periodic homogenization for
mechanics. Introduced in a rather complex way, a simplified description
is proposed below. The main idea is that the temperature gradient is
evaluated on two grids, both located at voxels' centres (same
location as with the HEX8R scheme). Referring to a finite element
perspective, they can be regarded as two Integration Point (IP) grids
(instead of a single one for HEX8R). For each IP grid, the temperature
gradient is evaluated with one of the two tetrahedral FD schemes, T1
and T2 (see Figure~\ref{fig1}). As the two tetrahedras are defined
within an hexaedra, this finite difference scheme has the same support
as HEX8R (i.e.\ the same points are used).


In practice, the two temperature gradient fields, $\bnabla
_{\mathbf{T1}}\mathbi{T}$ and $\bnabla
_{\mathbf{T2}}\mathbi{T}$, are evaluated using a
tetrahedral Finite difference (using Equation (\ref{eqn5}) with $\Omega
_{\mathrm{D}}=\Omega _{\mathrm{T1}}$ or $\Omega_{\mathrm{D}}=\Omega 
_{\mathrm{T2}}$, the two tetrahedras
defined on Figure~\ref{fig1}). Taking into account the two associated
fluxes $\mathbi{q}_{\mathbf{T1}}$ and
$\mathbi{q}_{\mathbf{T2}}$ (with
$\mathbi{q}_{\mathbf{Ti}}=-\underline{\mathbi{K}}\cdot\bnabla
_{\mathbf{Ti}}\mathbi{T}$), the discrete flux divergence is
defined from a mixture of the two fluxes,
$\mathbi{q}_{\mathbf{T1}}$ and
$\mathbi{q}_{\mathbf{T2}}$, as follows:
{\begin{eqnarray}\label{eqn6}
\mathrm{div}_{\mathrm{T}}(\mathbi{q})=
\tfrac{1}{2}(\mathrm{div}_{\mathrm{T1}}
(\mathbi{q}_{\mathbf{T2}})+
\mathrm{div}_{\mathrm{T2}}(\mathbi{q}_{\mathbf{T1}})),
\end{eqnarray}}\unskip
so that, the implicit Equation (\ref{eqn2}) becomes:
{\begin{eqnarray}\label{eqn7}
\rho c_{p}\frac{\mathrm{d}T}{\mathrm{d}t}-\frac{1}{2}(\mathrm{div}_{\mathrm{T1}}
(\underline{K}\cdot\bnabla _{\mathbf{T2}}
(T^{t}+\mathrm{d}T))+\mathrm{div}_{\mathrm{T2}}(\underline{K}\cdot
\bnabla _{\mathbf{T1}}(T^{t}+\mathrm{d}T)))=0,
\end{eqnarray}}\unskip
and finally, the fix-point algorithm corresponding to Equation
(\ref{eqn3}) reads:
{\begin{eqnarray}
&&
\frac{(\rho c_{p})_{0}}{\mathrm{d}t}\mathrm{d}T-\theta \mathrm{k}_{0}
\Delta _{\mathrm{T1}\mathrm{T2}}\,\mathrm{d}T = 
\frac{(\rho c_{p})_{0}-\rho c_{p}}{\mathrm{d}t}\,\mathrm{d}T \nonumber\\&&\quad
 +  \frac{1}{2}(\mathrm{div}_{\mathrm{T2}}(\underline{K}\cdot
\bnabla _{\mathbf{T1}}(T^{t}+\mathrm{d}T)-\theta k_{0}
\bnabla _{\mathbf{T1}}(\mathrm{d}T))
 + \, \mathrm{div}_{\mathrm{T1}}(\underline{K}\cdot
\bnabla _{\mathbf{T2}}(T^{t}+\mathrm{d}T)-\theta k_{0}
\bnabla _{\mathbf{T2}}(\mathrm{d}T))) 
\label{eqn8}
\end{eqnarray}}\unskip

As previously, evaluating in real space the r.h.s term of Equation
(\ref{eqn8}), denoted $p$, from a given $\mathrm{d}T$ field, is the first step
of the fix-point algorithm, and the second step consists of inverting
the lhs term in Fourier space to determine a new $\mathrm{d}T$, used to
evaluate a new rhs term and so on. The notation $\Delta
_{\mathrm{T1}\mathrm{T2}}$ stands for
$\sfrac{1}{2}(\mathrm{div}_{\mathrm{T1}}(\bnabla
_{\mathbf{T2}}(\cdot))+\mathrm{div}_{\mathrm{T2}} (\bnabla
_{\mathbf{T1}}(\cdot)))$. In Fourier space, the divergence
and gradient operators reads:
{\begin{eqnarray}
\begin{array}{c}
\displaystyle
\widehat{\bnabla _{\mathbf{Ti}}\mathbi{T}}=\mathrm{i}
\hat{\mathrm{T}}\boldsymbol{\xi }_{\mathbf{Ti}}=\mathrm{i}
\hat{\mathrm{T}}\mathbi{d}_{\mathbf{Ti}} \exp 
\left(\frac{\mathrm{i}(\xi _{x}+\xi _{y}+\xi _{z})}{2}\right)
\vspace*{6pt}\\
\displaystyle
\widehat{\mathrm{div}_{\mathrm{Ti}}(\mathbi{q})} =\mathrm{i}
\hat{\mathbi{q}}\cdot\boldsymbol{\xi }_{\mathbf{Ti}}=\mathrm{i}
\hat{\mathbi{q}}\cdot\mathbi{d}_{\mathbf{Ti}}\exp 
\left(-\frac{\mathrm{i}(\xi _{x}+\xi _{y}+\xi _{z})}{2}\right)
\end{array}
\label{eqn9}
\end{eqnarray}}\unskip



The definition of $\mathbi{d}_{\mathbf{Ti}}$ can be found in
\cite{5} and it must be emphasized that
$\mathbi{d}_{\mathbf{T1}}$ and
$\mathbi{d}_{\mathbf{T2}}$ are complex conjugate
$(\mathbi{d}_{\mathbf{T1}}=
\overline{\mathbi{d}_{\mathbf{T2}}})$. Hence, the
application of the operator} $\Delta _{\mathrm{T1}\mathrm{T2}}$  in
Fourier space corresponds to a multiplication by the term
$-\mathbi{d}_{\mathbf{T1}}\cdot
\overline{\mathbi{d}_{\mathbf{T1}}}=-
\Vert\mathbi{d}_{\mathbf{T1}}\Vert^2
(=-\Vert\mathbi{d}_{\mathbf{T2}}
\Vert^2)$. Finally, the lhs term of
Equation (\ref{eqn8}) can be easily inverted, replacing
$\Vert\boldsymbol{\xi}_{\mathrm{D}}\Vert^{2}$ by
$\Vert\mathbi{d}_{\mathbf{T1}}\Vert^2$ 
in Equation (\ref{eqn4}).


\begin{proof}[Remarks]
This presentation is rather different from the presentation
proposed by Finel \cite{5}. It is believed
simplified (without the introduction of FCC subgrids) and offers a
quite general and versatile framework for mixing finite differences
operators. In addition, this description (as well as its
implementation) is independent of the grid size parity (odd or even
number of voxels), whereas Finel limits its usage to even grids.
\let\qed\relax
\end{proof} 


\subsection{Additional details}

To reduce the number of implicit iterations, the initialization of the
fix point algorithm, at each time step, consists of a linear
extrapolation from the solution fields $\mathrm{d}T_{n-2}$ and $\mathrm{d}T_{n-1}$
evaluated at the last two times steps.


In the present implementation, the fix-point algorithm can be
accelerated by an Anderson acceleration technique \cite{1} as proposed
in \cite{4} for accelerating FFT-based solvers in static mechanics.
However, the benefit of the acceleration was not so clear in the case
presented in Section~\ref{sec4}, and numbers of iterations in
Table~\ref{tab1} are given for the classical fix-point algorithm.

\begin{table}[h!]
\caption{Number of iterations as a function of the step size and the
explicit or implicit method (up to the time of interest
$t=300\mathrm{d}t_{\mathrm{CFL}}^{*}$)\label{tab1}}
\begin{tabular}{cccc}
\thead
& Number of steps & HEX8R & TETRA2 \\
\endthead 
Explicit $\mathrm{d}t=\mathrm{d}t_{\mathrm{CFL}}^{*}/2$\vspace*{4pt} & 600 & 600$^*$ & 600$^*$ \\
Implicit $\mathrm{d}t=\mathrm{d}t_{\mathrm{CFL}}^{*}/2$\vspace*{4pt} & 600 & 2876 & 3257 \\
Implicit $\mathrm{d}t=\mathrm{d}t_{\mathrm{CFL}}^{*}$\vspace*{4pt} & 300 & 2406 & 2712 \\
Implicit $\mathrm{d}t=2\mathrm{d}t_{\mathrm{CFL}}^{*}$ & 150 & 2246 & 2584 
\botline
\end{tabular}
\tabnote{*~For the explicit method, each
iteration is a direct evaluation (no iterative back and forth in
Fourier space).}
\end{table}

The choice made for the reference thermal conductivity comes from the
proposition of Moulinec for isotropic materials \cite{14}:
$k_{0}=\sfrac{1}{2}(\min (k)+\max (k))$. If materials are anisotropic,
an ``isotropized'' thermal conductivity can be take into account.
Note that the choice of $k_{0}$ do not modify the solution but the
convergence rate of the algorithm. 


For the sake of brevity in Section~\ref{sec4}, only the conductivity is
considered heterogeneous and ($\rho c_{p}$) is assumed homogeneous.
Note that simulations performed with heterogeneous ($\rho c_{p}$)
demonstrate similar results regarding spurious oscillations. 


Finally, the relative error between two successive iterates is used as
convergence criterion with a low threshold $\varepsilon _{0}$ to
ensure the quality of the results:
{\begin{eqnarray}\label{eqn10}
\frac{\Vert \mathrm{d}T^{k}-\mathrm{d}T^{k-1}\Vert_{2}}{\Vert
\mathrm{d}T^{k}\Vert_{2}}< \varepsilon _{0}=10^{-6}
\end{eqnarray}}\unskip


\section{Validation on a homogeneous case}

The algorithm and its implementation is partly validated on the
homogeneous case for which the analytical solution is known.
Considering the initial and boundary conditions described in
Section~\ref{sec2.1}, the 3D problem reduces to a 1D problem whose
solution reads:
{\begin{eqnarray}\label{eqn11}
T(x,t)=\left(1-\frac{x}{L}\right)T^{\mathrm{left}}+
\sum _{k\in \mathbb{N}^{*}}-\frac{2T^{\mathrm{left}}}{k\uppi }\exp 
\left(-\left(\frac{k\uppi }{L}\right)^{2}\frac{K}{\rho c_{p}}t\right)
\sin \left(\frac{k\uppi }{L}x\right)
\end{eqnarray}}\unskip

The parameters used for the validation are $L=1~\mathrm{m}$,
$K=10~\mathrm{W}{\cdot}\mathrm{m}^{-1}{\cdot}\mathrm{K}^{-1}$, $\rho
=1000~\mathrm{kg}{\cdot}\mathrm{m}^{-3}$,
$c_{p}=1000~\mathrm{J}{\cdot}\mathrm{kg}^{-1}{\cdot}\mathrm{K}^{-1}$,
$T^{\mathrm{left}}=100~\mathrm{K}$ and $T^{0}=0~\mathrm{K}$. For the
numerical simulation, the 3D cubic domain is discretized in $32^{3}$
voxels (even if the solution is known invariant in directions $y$ and
$z$). 


\subsection{Stability of the methods}

The numerical experiments performed with the HEX8R and TETRA2 schemes
combined with the explicit Euler method ($\theta =0$) were stable as
soon as $\mathrm{d}t\leq \sfrac{\mathrm{d}t_{\mathrm{CFL}}}{2}$, with $\mathrm{d}t_{\mathrm{CFL}}=\sfrac{\rho
c_{p}\mathrm{d}x^{2}}{K}$ and $\mathrm{d}x$ the voxel dimension, and systematically
stable with the implicit Euler method ($\theta =1$). These
observations are consistent with mathematical results obtained with an
equivalent 1D algorithm, a $\theta$-method in time and centred
second derivative in space (i.e.~$\sfrac{\mathrm{d}^{2}T}{\mathrm{d}x^{2}}\sim
\stfrac{T_{n+1}+T_{n-1}-2T_{n}}{\mathrm{d}x^{2}}$). Actually, considering the
loading and the homogeneous material, the solution is uniaxial and, in
that context, both 3D schemes (HEX8R and TETRA2) turns out equivalent
to this standard and well-known 1D algorithm that is conditionally
stable for $\theta =0$ (stable if $\mathrm{d}t\leq \sfrac{\mathrm{d}t_{\mathrm{CFL}}}{2}$), and
unconditionally stable if $\theta \geq 1/2$.



\subsection{Temperature profiles comparison}


The times of interest for that comparison are $2\mathrm{d}t_{\mathrm{CFL}}$,
$6\mathrm{d}t_{\mathrm{CFL}}$ and $18\mathrm{d}t_{\mathrm{CFL}}$, and the considered time
steps are $\sfrac{\mathrm{d}t_{\mathrm{CFL}}}{2}$, the limit of stability for the explicit
method, and $\mathrm{d}t_{\mathrm{CFL}}$, only accessible to the implicit
method. Figure~\ref{fig2}(left) plots the temperature profiles
comparing analytical and numerical profiles. As expected from the
remark in the previous subsection, HEX8R and TETRA2 schemes provide
exactly the same results (``$+$'' and ``$\times$'' symbols are
superimposed), in good agreement with the analytical profiles. The
normalized error between analytical and numerical results,
$(T_{\mathrm{an.}}-T_{\mathrm{num.}})/T^{\mathrm{left}}$, is plotted on
Figure~\ref{fig2}(right). For time $t=2\mathrm{d}t_{\mathrm{CFL}}$, results for
the explicit method ($\theta =0$), used with $\mathrm{d}t=\sfrac{\mathrm{d}t_{\mathrm{CFL}}}{2}$
(its limit of stability), exhibits an error of ${\sim}$6\% that reduces
to ${\sim}$3\% with the implicit method ($\theta =1$). With a larger
time of interest, $t=6\mathrm{d}t_{\mathrm{CFL}}$ (dashed lines), the error
decreases for both methods but the implicit scheme remains more precise
(still with $\mathrm{d}t=\sfrac{\mathrm{d}t_{\mathrm{CFL}}}{2}$). Using a twice-larger step size
(i.e.\ $\mathrm{d}t=\mathrm{d}t_{\mathrm{CFL}}$) the explicit method provides a similar
precision w.r.t.\ the explicit method (with a twice-smaller step). As
the purpose is more the validation of the solver and its implementation
than a detailed comparison between the two methods (explicit vs
implicit), widely explored in handbooks for this simple 1D case, the
next step focuses on the heterogeneous case, for which HEX8R and TETRA2
schemes are no longer equivalent.

\begin{figure}
\includegraphics{fig02}
\caption{\label{fig2}Transient diffusion in a homogeneous material. 
(left) Numerical temperature profiles compared to analytical solution
for  $t=2\mathrm{d}t_{\mathrm{CFL}}$, $6\mathrm{d}t_{\mathrm{CFL}}$ and
$18\mathrm{d}t_{\mathrm{CFL}}$. (right) Normalized error
$(T_{\mathrm{an.}}-T_{\mathrm{num.}})/T^{\mathrm{left}}$, for time
$t=2\mathrm{d}t_{\mathrm{CFL}}$ with continuous lines, and
$t=6\mathrm{d}t_{\mathrm{CFL}}$  with dashed lines. In the legend,
$\mathrm{d}t=1/2$ stands for $\mathrm{d}t=\mathrm{d}t_{\mathrm{CFL}}/2$, and
$\mathrm{d}t=1$ for $\mathrm{d}t=\mathrm{d}t_{\mathrm{CFL}}$.}
\end{figure}


\section{Application to the heterogeneous case}\label{sec4}

A simple case of application is proposed with a single spherical
inclusion of radius 0.3~m at the position (0.15~m, 0.4~m, 0.6~m) within
a cubic domain of size $1\times 1\times 1~\mathrm{m}^{3}$, discretized
in $32 \times 32 \times 32$ voxels. The mass density and specific heat
capacity are assumed homogeneous (same values as previously), while the
thermal conductivity is set to
$10^{4}~\mathrm{W}{\cdot}\mathrm{m}^{-1}{\cdot}\mathrm{K}^{-1}$. Compared to the
conductivity of the matrix,
$10~\mathrm{W}{\cdot}\mathrm{m}^{-1}{\cdot}\mathrm{K}^{-1}$, the contrast is quite
high (1000). As spurious oscillations generally correlate with high
contrasts, this case should allow for discriminating the two FD
schemes, which were equivalent in the homogeneous case used in the
previous section.


\subsection{Stability of the method}

The stability condition for the explicit method ($\theta =0$) has been
explored empirically on this example: whether combined with the HEX8R
or the TETRA2 scheme, the method is stable as soon as the step size
$\mathrm{d}t\leq \sfrac{{\mathrm{d}t}_{\mathrm{CFL}}^{*}}{2}$, with
$\mathrm{d}t_{\mathrm{CFL}}^{*}=\sfrac{\rho c_{p}\mathrm{d}x^{2}}{\max (K)}$, and unstable
otherwise. This condition also proved valid for simulations with a low
conductivity inclusion. Note that this definition of
$\mathrm{d}t_{\mathrm{CFL}}^{*}$ allows for using the same stability condition
than the homogeneous case with $\mathrm{d}t_{\mathrm{CFL}}$. For heterogeneous
$\rho c_{p}$ the condition proved valid for
$\mathrm{d}t_{\mathrm{CFL}}^{*}=\sfrac{\min (\rho
c_{p})\mathrm{d}x^{2}}{K}$. For both $K$ and
$\rho c_{p}$ heterogeneous, a similar definition for 
$\mathrm{d}t_{\mathrm{CFL}}^{*}$, if exists, is still under investigation.
Finally, all the tests performed with the implicit method ($\theta =1$)
remained stable.


\subsection{Temperature fields comparison}

For that comparison, the time of interest is
$t=300\mathrm{d}t_{\mathrm{CFL}}^{*}$, chosen to easily discriminate the two
schemes. Figure~\ref{fig3} displays the different temperature fields
for the two schemes HEX8R (up) and TETRA2 (bottom), obtained with the
explicit ($\theta =0$) or implicit ($\theta =1$) method (respectively
left and right side). Explicit simulations have been performed using
$\mathrm{d}t=\sfrac{{\mathrm{d}t}_{\mathrm{CFL}}^{*}}{2}$ and implicit simulations with
$\mathrm{d}t=2{\mathrm{d}t}_{\mathrm{CFL}}^{*}$. The HEX8R scheme induces important
spurious oscillations, whether using an implicit or explicit method.
Importantly, an overshoot is observed on the colorbar: the temperature
field can be locally higher than the applied temperature
$T^{\mathrm{left}}=100~\mathrm{K}$. On the other hand, the TETRA2
scheme used with the implicit or explicit method, do not exhibit any
sign of spurious oscillations or overshoot: the field is smooth within
each material and the transition is sharp at the interface between the
two materials. Finally, the quasi-identical fields obtained with the
explicit method and the implicit method with 4 times larger step
demonstrates the interest of the implicit method in that case
(especially for quite large times of interest).


\begin{figure}
\includegraphics{fig03}
\caption{\label{fig3}Temperature field distribution with the HEX8R
scheme (up) and the TETRA2 scheme (bottom), using the explicit method
($\theta=0$, $\mathrm{d}t=\mathrm{d}t_{\mathrm{CFL}}^{*}/2$, (left)),
or the implicit method  ($\theta=1$,
$\mathrm{d}t=2\mathrm{d}t_{\mathrm{CFL}}^{*}$, (right)). Note the
oscillations and overshoot (values higher than 100~K on the colorbar).}
\end{figure}


As a conclusion regarding the quality of the solution, using the TETRA2
scheme appears a relevant choice to avoid any artefact such as spurious
oscillations and overshoot. This is a crucial point especially when
considering the solution field as an input for another coupled physics
(mechanics or phase transition, for example). 


\subsection{Performance}

The purpose of this subsection is not to provide a deep analysis but
rather a few elements of discussion on the choice of the step size, or
between explicit or implicit method. Table~\ref{tab1} gathers the
cumulated number of iterations up to the time of interest
$t=300\mathrm{d}t_{\mathrm{CFL}}^{*}$ with different time
steps for the implicit method (all of them generate almost identical
fields, as observed~in~Figure~\ref{fig3}). Considering the choice
between HEX8R and TETRA2, the number of iterations is a bit higher for
TETRA2. In addition, each iteration of the TETRA2 scheme is a bit more
complex as the flux is evaluated on two grids of centres (see
Section~\ref{sec2}). However, the price to pay, regarding the
important enhancement of the solution quality, is rather low.
Considering the choice of the step size for implicit algorithms: it is
worth noting that while the step size is multiplied by 2 the number of
iterations is only divided by a coefficient between 1.02 and 1.2.
Actually, when increasing the step size, the number of fix-point
iterations, required to solve each step, increases (especially at the
beginning of the simulation). Finally, considering the choice between
explicit and implicit algorithm, the difference in terms of iterations
is rather important. In addition, an iteration of the explicit
algorithm has a much lower cost as it consists of a direct evaluation
of $\mathrm{div}(\mathbi{q})$ in real space and does not require
any back and force in Fourier space. From that point of view, and for
that given example (material, loading and time of interest), the
explicit method appears much more performant than the implicit scheme. 

\section{Example on a larger unit-cell}\label{sec5}

The present section is just an example showing that, despite of a
non-optimized \textit{python} implementation for this proof of concept,
the method applies straightforwardly to a larger microstructure that
consists of 125 randomly distributed spheres with a volume fraction of
22\%. The size of the unit-cell is $1 \times 1 \times
1~\mathrm{m}^{3}$, discretized in $128 \times 128 \times 128$ voxels.
Microstructure is generated with periodicity condition and periodic
conditions are considered in the transverse directions (wrt to the
temperature gradient direction) and the same conditions are applied for
the upper and lower face (100~K and 0~K homogeneous applied
temperature). Initial temperature is 0~K. The mass density and specific
heat capacity are assumed homogeneous (same values as previously),
while the thermal conductivity is set to
$10^{4}~\mathrm{W}{\cdot}\mathrm{m}^{-1}{\cdot}\mathrm{K}^{-1}$ for the inclusion and
$100~\mathrm{W}{\cdot}\mathrm{m}^{-1}{\cdot}\mathrm{K}^{-1}$ for the matrix. The
temperature field displayed on Figure~\ref{fig4} is evaluated after
$5000\mathrm{d}t_{\mathrm{CFL}}^{*}$. A second simulation considers
insulating inclusions with a thermal conductivity of
$1~\mathrm{W}{\cdot}\mathrm{m}^{-1}{\cdot}\mathrm{K}^{-1}$. The temperature field is
evaluated at the same time (corresponding to
$50\mathrm{d}t_{\mathrm{CFL}}^{*}$ in that case). On this more
realistic application, the temperatures fields are still very smooth
and comparison clearly exhibits the influence of the inclusions,
whether conducting or insulating, on the propagation of the temperature
with the composite material.




\begin{figure}[h!]
\includegraphics{fig04}
\caption{\label{fig4}Example exhibiting the influence of conducting
(left) or insulating (right) spherical inclusions within a matrix on
the temperature field after the same propagation time (the simulated
unit-cell is shown at the centre).} 
\end{figure}


Important remark: it must be mentioned that, in this example, the
inclusions are well separated from each others. Various tests, not
reported in this short note, have revealed the existence of spurious
oscillations in cases where inclusions are ``touching'' each other(s).
If the amount of oscillation is still less than HEX8R scheme, this
point must be kept in mind and needs further investigations. 


\section{Conclusion and prospects}

As a conclusion of this letter, it is demonstrated that combining
recent developments using DTTs, to apply non-periodic boundary
conditions, with a recent Finite Difference scheme (noted TETRA2),
allows for taking advantage of the efficiency of FFT-based solvers for
simulating transient diffusion in heterogeneous materials with
realistic boundary conditions. As demonstrated, the method is robust in
time (using an implicit solver), accurate in space as the solution
fields do not exhibit any spurious oscillations (at least if the
material consists of a matrix with well separated inclusions), and
accurate in term of ``locality'' as the TETRA2 finite difference
scheme only considers the 8 points of an hexahedra to approximate the
divergence or the gradient. In that sense, it is as local as the HEX8R
scheme but do not suffer from spurious oscillations originating from
the hourglass phenomenon. To the author's opinion, the presentation
of the TETRA2 scheme, proposed by Finel in his original paper \cite{5},
has been simplified, generalized (towards other mixing of finite
differences) and a bit extended (its usage is not limited to even
grids). From the numerical experiments in Section~\ref{sec4}, using the
TETRA2 scheme is a bit less efficient than the HEX8R scheme, however
the quality enhancement is so important that theTETRA2 choice should be
favoured (author's opinion). Regarding the choice of explicit or
implicit method, and step size for the implicit method, the question
has been addressed for a specific case of application (homogeneous
$\rho c_{p}$, high conductivity contrast, large time of interest) and
the explicit method proved more efficient than the implicit one.
However, a strategy mixing Euler explicit and implicit together with an
adaptive step size will probably be an efficient way to take benefit
from the advantages of the two methods.


As mentioned in Section~\ref{sec5}, for matrix inclusion
microstructures with inclusions touching each other, spurious
oscillations appear on the temperature field with the TETRA2 scheme. If
less important than with the HEX8R scheme, this point needs further
investigations. The use of composite voxels is generally an efficient
way to smooth oscillations is such situations~\cite{3}.


As a future prospect, the present FFT-based solver ($\theta$-method in
time and TETRA2 scheme in space) can be considered as a reliable and
essential component for coupling transient diffusion with mechanics.
Important applications in the context of materials science consider the
diffusion of species strongly coupled with mechanics through a
stress-free strain directly proportional to the concentration of
diffusing atoms (e.g.\ \cite{22,10,13}). Such couplings, generally
performed with standard Finite Element methods, could now take benefit
from the present efficient, robust and accurate FFT-based solver
adapted to heterogeneous materials and realistic boundary\break conditions. 


Finally, an implementation of this framework is in progress in the
FFT-based code AMITEX~\cite{8} in order to take advantage of massively
parallel computers and deal with large scale\break simulations.



\section*{Declaration of interests}
The authors do not work for, advise, own shares in, or receive funds
from any organization that could benefit from this article, and have
declared no affiliations other than their research organizations.

\back{}

\appendix{}

\def\appendixlabel{Appendix A}

\section{DTTs and boundary conditions}\label{appA}
Table~\ref{tab2} gathers the different discrete symmetry types with the
corresponding types of DTT (according to the names defined in the
\textit{FFTW} library \cite{6}) and the associated types of boundary
conditions.

For the definition of the symmetry TXTY: 
\begin{itemize}
\item T has the value W, if the symmetry axis are located on the last
points, or H, if the axis is located a half element just after the last
points. In the present paper, the temperature being located at
voxels' corners, only W types of symmetry are considered,
\item X has the value S or A, whether the signal is Symmetric or
Anti-symmetric on the left point, 
\item Y has the value S or A, whether the signal is Symmetric or
Anti-symmetric on the right point. 
\end{itemize}
Finally, an Anti-symmetry is associated to a Dirichlet BC and a
Symmetry to a Neumann BC.


\begin{table}[h!]
\caption{DTTs and boundary conditions\label{tab2}}
\begin{tabular}{ccc}
\thead
Discrete symmetry type (TXTY) & 
Discrete Trigo. Transf. \cite{6} & Boundary condition (left/right) \\
\endthead
WSWS & DCT1 & Neumann/Neumann \\ 
HSHS & DCT2 &  \\ 
WSWA & DCT3 & Neumann/Dirichlet \\ 
HSHA & DCT4 &  \\ 
WAWA & DST1 & Dirichlet/Dirichlet \\ 
HAHA & DST2 &  \\ 
WAWS & DST3 & Dirichlet/Neumann \\ 
HAHS & DST4 &   
\botline
\end{tabular}
\vspace*{-12pt}
\end{table}

\bibliographystyle{crunsrt}
\bibliography{crmeca20240874}
\refinput{crmeca20240874-reference.tex}

\end{document}
