%~Mouliné par MaN_auto v.0.37.1 (7875ea12) 2025-11-03 10:07:29
\documentclass[CRMECA,Unicode,biblatex,published]{cedram}

\addbibresource{CRMECA_de-Chaisemartin_20241072.bib}

\usepackage{listings}
\usepackage{siunitx}
\usepackage{pgfplots}
\usepackage{subcaption}

\newcommand{\ArcNum}{ArcNum}

\newcommand{\Modules}{\emph{Modules}}
\newcommand{\Module}{\emph{Module}}
\newcommand{\Services}{\emph{Services}}
\newcommand{\Service}{\emph{Service}}
\newcommand{\Contributions}{\emph{Contributions}}
\newcommand{\Contribution}{\emph{Contribution}}

\newcommand{\perm}{\mathbb{K}}
\newcommand{\diff}{\mathbb{D}}
\newcommand{\diffT}{\mathbb{L}}

\newcommand{\VEC}[1]{\boldsymbol{#1}}
\newcommand{\vv}{\boldsymbol{v}}

\newcommand{\GRAD}{\nabla}

\newcommand{\emphb}[1]{#1}

\let\div\relax
\DeclareMathOperator{\div}{div}

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\usepackage{mathtools}
\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}

%%%---------------------------------------------------------------------------

%%% FIGURES

\makeatletter
\renewcommand\thesubfigure{(\alph{subfigure})}
\DeclareCaptionLabelFormat{mysublabelfmt}{(\alph{subfigure})}
\makeatother
\captionsetup[sub]{labelfont=bf,labelsep=space,labelformat=mysublabelfmt}

%%%---------------------------------------------------------------------------

%%% PACKAGE LISTINGS

\definecolor{lstbackgroundcolor}{RGB}{245,245,244}
\definecolor{XMLBlue} {RGB} {40,40,255}
\definecolor{XMLGreen} {RGB} {0,96,96}
\definecolor{XMLRed} {RGB} {128,0,0}
\definecolor{mpi}{RGB}{65,105,225}
\definecolor{linear}{RGB}{255,31,91}

\newlength{\MaxSizeOfLineNumbers}%
\settowidth{\MaxSizeOfLineNumbers}{99}%
\addtolength{\MaxSizeOfLineNumbers}{2.5ex}%

\lstdefinelanguage{customXML}
{ 
  xleftmargin=\MaxSizeOfLineNumbers,
  basicstyle=\footnotesize\ttfamily,
  captionpos=b,
  columns=fixed,       
  numbers=left,                                   
  frame=lines,                                         
  backgroundcolor=\color{lstbackgroundcolor},
  numberstyle=\tiny\color{darkgray},                        
  commentstyle=\it\color{XMLGreen},   
  stringstyle=\color{XMLBlue},
  identifierstyle=\color{XMLRed},
  showstringspaces=false,
  morestring=[b]",
  moredelim=[s][\bfseries\color{black}]{<}{\ },
  moredelim=[s][\bfseries\color{black}]{</}{>},
  moredelim=[l][\bfseries\color{black}]{/>},
  moredelim=[l][\bfseries\color{black}]{>},
  morecomment=[s]{<?}{?>},
  morecomment=[s]{<!--}{-->},
}

\definecolor{CPPBlue}   {RGB} {40,40,255}
\definecolor{CPPGreen}  {RGB} {0,96,96}
\definecolor{CPPRed}    {RGB} {128,0,0}
\definecolor{CPPViolet} {HTML} {7040A0}

\lstdefinelanguage{customCpp}
{
%	escapeinside={(*@}{@*)},
	xleftmargin=\MaxSizeOfLineNumbers,
    basicstyle=\footnotesize\ttfamily, 
    captionpos=b,
    columns=fixed,   
    numbers=left,                                   
    frame=lines,
    backgroundcolor=\color{lstbackgroundcolor},                   
    numberstyle=\tiny\color{darkgray},         
    commentstyle=\it\color{CPPGreen},              
    stringstyle=\rmfamily\slshape\color{CPPRed},
    keywordstyle=\color{CPPBlue},  
    showstringspaces=false,       
    mathescape=true,
    escapechar=\!,
    language=c++,                              
    morekeywords={alignas,continute,friend,register,true,alignof,decltype,goto,
    reinterpret_cast,try,asm,defult,if,return,typedef,auto,delete,inline,short,
    typeid,bool,do,int,signed,typename,break,double,long,sizeof,union,case,
    dynamic_cast,mutable,static,unsigned,catch,else,namespace,static_assert,using,
    char,enum,new,static_cast,virtual,char16_t,char32_t,explict,noexcept,struct,
    void,export,nullptr,switch,volatile,class,extern,operator,template,wchar_t,
    const,false,private,this,while,constexpr,float,protected,thread_local,
    const_cast,for,public,throw,std,
    % add RG
    ENUMERATE_,ENUMERATE_NODE,ENUMERATE_FACE,ENUMERATE_PHASE,ENUMERATE_SUBSYSTEM,ENUMERATE_SPECIES,ENUMERATE_CELL, RUNCOMMAND_ENUMERATE, RUNCOMMAND_LOOP,RUNCOMMAND_LOOP1, RUNCOMMAND_},
    emph=[1]{Arcane,map,set,multimap,multiset,unordered_map,unordered_set,
    unordered_multiset,unordered_multimap,vector,string,list,deque,
    array,stack,forwared_list,iostream,memory,shared_ptr,unique_ptr,
    random,bitset,ostream,istream,cout,cin,endl,move,default_random_engine,
    uniform_int_distribution,iterator,algorithm,functional,bing,numeric,
    % add RG
    ArcRes,ArcNum,Law,Accelerator,math,XPath,FluidDensityLaw,FluidDensityLawType},
    emphstyle=[1]\color{CPPViolet}, 
    emph=[2]{IMesh,VariableCellInt32,VariableNodeInt32,Node,Cell,Integer,VariableCellReal,Runner,VariableBuildInfo,Real,property,FluidDensityLawConfigService,FunctionManager,System,Signature,Function,FunctionEvaluator,IVariableMng,Vector,Matrix,Pressure,CapillaryPressure,FluidDensity,FaceGroup,TwoPointsStencil,contribution,values,Temperature,VolumeFraction,Density,HeatConductivity,Viscosity,VariableFaceReal,RelativePermeability,Space,MatrixDistribution,VectorDistribution,DirectMatrixBuilder},
    emphstyle=[2]\color{CPPGreen},
}

%%%---------------------------------------------------------------------------

%%% VERTICAL BAR

%%% ''Evaluated at'' macro
\newcommand\restr[3][]{{				
 \left.\kern-\nulldelimiterspace 
 #2 								
 \vphantom{\big\vert} 				
 \right\vert_{#3}^{#1}
 }}

%%%---------------------------------------------------------------------------


\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}}
\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel}


\title{ArcNum: an Arcane-based numerical framework used in porous media flow simulation applications}
\alttitle{ArcNum\tralicstex{}{~}: un framework numérique à base d'Arcane utilisé pour des simulations d'écoulements en milieu poreux}

\author[S. de Chaisemartin]{\firstname{Stéphane} \lastname{de Chaisemartin}\IsCorresp}
\address{IFPEN, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex, France}
\email{stephane.de-chaisemartin@ifpen.fr}

\author{\firstname{Sylvain} \lastname{Desroziers}}
\address[1]{IFPEN, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex, France}
\address{Michelin, 23 Place des Carmes Déchaux, 63040 Clermont Ferrand Cedex 9, France}

\author{\firstname{Guillaume} \lastname{Enchéry}}
\address[1]{IFPEN, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex, France}

\author{\firstname{Raphaël} \lastname{Gayno}}
\address[1]{IFPEN, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex, France}

\author{\firstname{Jean-Marc} \lastname{Gratien}}
\address[1]{IFPEN, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex, France}

\author{\firstname{Gilles} \lastname{Grospellier}}
\address{CEA, DAM, DIF, 91297 Arpajon, France}

\author{\firstname{Thomas} \lastname{Guignon}}
\address[1]{IFPEN, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex, France}

\author{\firstname{Pascal} \lastname{Havé}}
\address[1]{IFPEN, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex, France}
\address{NUANT, Baarerstrasse 20, 6300 Zug, Suisse}

\author{\firstname{Benoît} \lastname{Lelandais}}
\address[6]{CEA, DAM, DIF, 91297 Arpajon, France}

\author{\firstname{Alexandre} \lastname{l'Héritier}}
\address[6]{CEA, DAM, DIF, 91297 Arpajon, France}

\author{\firstname{Anthony} \lastname{Michel}}
\address[1]{IFPEN, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex, France}

\author{\firstname{Aboul Karim} \lastname{Mohamed El Maarouf}}
\address[1]{IFPEN, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex, France}

\author{\firstname{Valentin} \lastname{Postat}}
\address[6]{CEA, DAM, DIF, 91297 Arpajon, France}

\author{\firstname{Xavier} \lastname{Tunc}}
\address[1]{IFPEN, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex, France}

\author{\firstname{Soleiman} \lastname{Yousef}}
\address[1]{IFPEN, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex, France}

\keywords{\kwd{Mathematical software framework} \kwd{code generation} \kwd{automatic differentiation} \kwd{Arcane} \kwd{proxy-app} \kwd{multiphase porous media flow simulation} \kwd{massively parallel computing}}
\altkeywords{\kwd{Composants logiciels de simulation} \kwd{génération de code} \kwd{différenciation automatique} \kwd{Arcane} \kwd{démonstrateur} \kwd{simulation d'écoulement polyphasiques en milieu poreux} \kwd{calcul massivement parallèle}}

\editornote{Article submitted by invitation}
\alteditornote{Article soumis sur invitation}

\begin{abstract}
This article presents ArcNum, a framework built on the Arcane platform, designed to easily and efficiently develop and maintain the numerical core in finite-volume and finite-element applications. This framework first enables the automatic generation of code needed to instantiate and handle complex physical models from a textual description, through its component called GUMP. ArcNum also offers software components to create, register and evaluate a set of physical laws, requiring only the specification of their inputs, outputs, and corresponding mathematical formulations. Finally, the framework includes a component named Contribution, which combines law evaluation with automatic differentiation to assemble linear systems efficiently. The framework ArcNum has been used to develop an open source Arcane-based porous media flow simulation proxy-app, named ShArc. Single and two-phase porous media flow simulations performed with ShArc are presented to complete the framework description. In order to illustrate the ability to use ArcNum for High Performance Computing, massively parallel simulations conducted with ShArc are finally presented.
\end{abstract}

\begin{altabstract}
Cet article présente ArcNum, un framework conçu au-dessus de la plateforme Arcane pour faciliter le développement et la maintenance du coeur numérique des applications volumes finis ou éléments finis. Cette structure logicielle permet tout d'abord de générer le code nécessaire à l'instanciation et à la manipulation de modèles physiques complexes à partir d'une description textuelle, grâce à son composant appelé GUMP. ArcNum propose également des composants logiciels pour créer, enregistrer et évaluer un ensemble de lois physiques, ne nécessitant que la description de leurs entrées et sorties et la fourniture des formules mathématiques. Enfin, le framework fournit un composant appelé Contribution qui permet de combiner l'évaluation des lois et la différenciation automatique en vue d'assembler un système linéaire. Les composants d'ArcNum ont été utilisés pour construire un démonstrateur \emph{open source} de simulation d'écoulement en milieu poreux basé sur Arcane, nommé ShArc. Des simulations d'écoulements en milieu poreux monophasique et diphasique, réalisées avec ShArc, sont présentées pour compléter la description du framework. Afin d'illustrer la capacité d'utiliser ArcNum pour le calcul à haute performance, des simulations massivement parallèles réalisées avec ShArc sont finalement présentées.
\end{altabstract}

\shortrunauthors

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

\maketitle

\section{Introduction}


The open source \emphb{Arcane} platform for scientific computing~\cite{arcane09} is an advanced framework designed to facilitate the development and execution of complex scientific simulations applications. It provides researchers, scientists, and engineers with a comprehensive set of tools and libraries specifically tailored for high-performance computing (HPC) applications.

Arcane started being developed at the \emphb{CEA-DAM}, the Military Applications Division of the French Alternative Energies and Atomic Energy Commission, in the early 2000s, before being co-developed with \emphb{IFPEN}, French Institute for Research and Innovation in Energy, Mobility and Environment, from 2007.

During these twenty years of existence, Arcane has been the foundation of several industrial applications, from laser physics to CO\textsubscript{2} storage or geothermal energy simulations. Today eight industrial simulators are based on this platform, for nearly three million code lines from CEA-DAM and IFPEN.

Starting in the 2020s it has been decided to gradually move toward an open source strategy to foster the collaborative development between CEA-DAM and IFPEN, to make the platform available for the community and to foster potential new collaborations. Since the end of 2021, the platform is fully available on GitHub~\cite{arcaneframework-github}.

Arcane offers services dedicated to the application developers, handling for them all the low level computer science services needed by the application and giving them tools to dynamically build the application with an expandable and flexible option system.

Above the Arcane platform, dedicated to handle data and parallelism, numerical frameworks were developed to share the numerical tools (numerical schemes, (non)linear solvers...) between the applications. The works on linear solvers gave birth to \emphb{Alien}~\cite{alien22}, a solver interface library designed to provide an application with a unique entry point to the world of linear solver libraries. On top of these (non)linear systems to solve, it was necessary to supply the applications with a way to simplify the complex writing of the Jacobian assembly arising for example in porous media flow simulations. \emphb{\ArcNum{}} was developed with this objective.

In order to share the use of such Arcane-based numerical frameworks with the community, a numerical demonstrator, \emphb{ShArc}, has been released in open source on GitHub~\cite{sharc-github}. This Arcane-based proxy-application uses the frameworks Alien and \ArcNum{} to solve porous-media flow problems. ShArc is now a proxy-app of Exa-DI NumPEx project; in this context an installation guide is provided on GitLab~\cite{sharc-doc-gitlab}.

This paper is especially dedicated to the \ArcNum{} framework and is organized in the following way. As an introduction, we first recall the main tools provided by the Arcane framework and highlight recent evolutions of the platform. We then present the \ArcNum{} framework itself, detailing its three main components:
\begin{enumerate}[(i)]
\item a physical data model generator named GUMP;
\item a physical law framework;
\item a tool to assemble the Jacobian contributions, associated with automatic differentiation.
\end{enumerate}
Finally, numerical results obtained with the ShArc demonstrator, combining all these tools, are presented for two porous media flow simulation test cases. In addition, a performance study of ShArc on massively parallel configurations is reported.

\begin{figure}[h!]
\includegraphics[width=0.8\textwidth]{ArcaneLoop.PNG}
\caption{The different stages of the Arcane loop. The init, compute-loop and exit stages are the main stages used by the application developer to plug its \Modules{} entry-points.}\label{fig:arcaneloop}
\end{figure}

\section{The Arcane Framework}

Before diving into the \ArcNum{} tools, we present in this part the Arcane Framework, an open source modular platform for HPC, on which it is based.

\subsection{An open source modular structure}

The Arcane Framework is made of several components, available on a single GitHub repository~\cite{arcane-github}:
\begin{description}
\item[Arccon] made of CMake functions and package finders;
\item[Arccore] containing the base types and parallel tools of Arcane;
\item[Axlstar] gathering the code-generation tools;
\item[Alien] the linear algebra API (cf.\ Section~\ref{subsec:Alien});
\item[Arcane] the Arcane platform itself (cf.\ Section~\ref{subsec:Arcane}).
\end{description}
With this modular structure, both Arcane and Alien use the common base provided by Arccore, the set of generators provided by Axlstar and the common CMake tools included in Arccon.

The Arcane framework committee decided in 2021 to release the platform with an \emphb{open source} \emphb{Apache-2.0} license. All Arcane Framework components are now available on GitHub, see~\cite{arcaneframework-github}, through the \emph{ArcaneFramework} organization, gathering the Arcane developers at CEA-DAM and IFPEN.

\subsection{Arcane: a complete platform for scientific computing}\label{subsec:Arcane}

The C++ platform Arcane provides all the building blocks needed to develop a parallel application. In particular, it provides a mesh structure with easy-to-use parallelization tools, distributed variables and parallel IOs, to name just a few of its features. It is based on an efficient component architecture allowing to define simulation problems as plugins or \Modules{}, in the platform language, and to share functionalities through \Services{}, defined by an interface, i.e.\ a class defining a programming contract. We refer to~\cite{arcane09} for a more detailed description of the Arcane platform. We here only briefly describe the main functionalities and then present the platform recent developments.

The Arcane framework was originally designed to conduct simulations through the call of various entry-points, implemented by the developer in the \Modules{}. Within the \Modules{} entry-points, the application developer has access, amongst others, to the simulation mesh, natively distributed over subdomains, the simulation variables, or the simulated case options. The main program of the application only consists in calling the Arcane launcher. The application developer can then focus on developing \Modules{}, and possibly \Services{} to share treatment between \Modules{}. The \Modules{} entry-points are registered in a \emph{time-loop}, mainly composed of three stages: the init, compute-loop and exit stages, see Figure~\ref{fig:arcaneloop}. Both \Modules{} and \Services{} are made up of the following elements:
\begin{enumerate}[(i)]
\item a descriptor file, in xml format, used to defined the entry-points for a \Module{} or the interface of a \Service{}, and a list of options;
\item a C++ class, used to implement the \Module{} entry-points and the \Service{} interface.
\end{enumerate}

An Arcane mesh is made of the classical geometric items such as cells, faces, nodes and edges. Non-geometric items, dofs (degrees of freedom) or particles, can also be added and connected to the geometric items. The platform can read various formats (vtk, msh and med formats) and generate a few mesh types too (see Section~\ref{subsec:Meshes}). It enables the use of several mesh partitioners (Metis, Parmetis, PTScotch and Zoltan) in a transparent manner. The Arcane mesh is intrinsically unstructured and dynamic, which means that we can add new cells during the simulation and enrich the existing mesh. We can also create dynamically group of items to identify regions of the mesh. The mesh variables are then adapted automatically to take into account these changes in an efficient way. Let us add that a simulation may use several meshes if needed. Distributed variables can be defined on this/these mesh(es). A mesh variable is linked to a specific item family (a family of cells, faces, edges, nodes or dofs) and can have scalar or vector values of various types: Real (Arcane name for double precision floating points), Int16, Int32, Int64...

Arcane is designed to make parallel operations easier. By providing the concept of \emph{SubDomain} with its \emph{own} and \emph{ghost} mesh items, it allows the application developer to easily implement an application based on domain decomposition. Ghost item values can be easily synchronized and explicit communications between the subdomains are not necessary. Arcane also offers load balance functionalities. The communications are classically carried out through mpi. Threaded and hybrid mpi+thread implementations are also available, though not used currently by the main Arcane-based applications.

An Arcane code example is given in Listing~\ref{list:MeshVarEx} to illustrate how mesh, variables and parallel synchronization operations are used to compute a variable defined on the mesh nodes from a second variable defined on the mesh cells. Note that, in this example, the group \verb+mesh->allNodes()+ only contains the nodes of the current subdomain. A final synchronization step is therefore required to ensure the correctness of the computation, as some ghost node values may remain inconsistent.

\begin{center}
\begin{lstlisting}[xleftmargin=\MaxSizeOfLineNumbers,language=customCpp, basicstyle=\small\ttfamily,caption=Example of the computation of a node variable from a cell variable in parallel., label={list:MeshVarEx}]
Arcane::IMesh* mesh; // given by Arcane
Arcane::VariableCellInt32 cell_var {VariableBuildInfo{
                                        mesh,"MyCellVariable"}};
Arcane::VariableNodeInt32 node_var {VariableBuildInfo{
                                        mesh,"MyNodeVariable"}};
cell_var.fill(1);
node_var.fill(0);
ENUMERATE_(Node, inode, mesh->allNodes()) {
    for (Cell cell : inode->cells()) {
        node_var[inode] += cell_var[cell];
    }
}
node_var.synchronize();
\end{lstlisting}
\end{center}

This abstraction of mesh enumeration and parallel communications allows applications developers to concentrate on the coding of physics and numerical operators. The Arcane level of abstraction remains low with direct access to mesh variables. This allows to easily adapt existing numerical algorithms, and offers the possibility to define finite element or finite volume methods on the top of this Arcane API.

Thanks to this subdomain-based structure, simulators at CEA and IFPEN have been used in parallel simulations on meshes ranging from thousands of millions up to dozens of billion of cells on up to hundreds of thousands of cores, see Section~\ref{sec:GrandChallenge}.

\subsection{Alien: a generic extensible linear algebra framework}\label{subsec:Alien}

The Arcane framework now integrates a generic and extensible linear algebra framework called \emphb{Alien}~\cite{alien22}. This linear algebra API was initiated at IFPEN in the 2010s. It was designed to answer the problems arising when using different linear algebra packages in scientific software, and to easily switch between different sparse matrix storage formats. Alien provides the application developer with a unique API to handle matrices, vectors and to solve linear systems. An example of system assembling and resolution is given in Listing~\ref{list:AlienEx}. These Alien matrices and vectors are then converted at runtime into the representation of the chosen linear algebra package among PETSc, HYPRE, Trilinos, MTL or SuperLU. Alien also gives access to IFPEN internal solvers such as IFPSolver and MCG Solver~\cite{IfpenSolvers14}. It has been shown in~\cite{alien22} that these data structure conversions don't introduce a significant overhead in CPU-time.

\begin{center}
\begin{lstlisting}[language=customCpp,  caption={Handling distributed Alien matrices, vectors and solvers.}, label={list:AlienEx}]
auto pm = Env::parallelMng();
auto s = Space{10,"MySpace"};
auto md = MatrixDistribution{s,s,pm}; 
auto vd = VectorDistribution{s,pm};
auto A = Matrix{md}; // build a 10$\textcolor{CPPGreen}{\times}$10 square matrix
auto x = Vector{vd}; // build a 10 elements vector
auto b = Vector{vd}; // build a 10 elements vector
{
    // fill a tri-diagonal matrix
    auto builder = DirectMatrixBuilder(A,eResetValues);
    builder.reserve(30); // non-zero values
    for (Integer index = 0; index < 10; ++ index) {
        builder(index,index) = 2;
        if (index + 1 < 10) builder(index, index + 1) = -1;
        if (index - 1 >= 0) builder(index, index - 1) = -1;
    }
}
auto* solver = createSolver(/*...*/); // choose your favorite solver
auto status = solver->solve(A,x,b)    
\end{lstlisting}
\end{center}

\subsection{Handling an always wider variety of mesh}\label{subsec:Meshes}

While initially designed only for standard unstructured meshes, the Arcane framework has been gradually enriched with:
\begin{enumerate}[(i)]
\item degenerate hexahedral cells,
\item generated Cartesian or honeycomb meshes,
\item Adaptive Mesh Refinement (AMR), both for unstructured meshes and Cartesian meshes,
\item polyhedral meshes,
\item uniform mesh refinement, to generate a finer mesh.
\end{enumerate}

All these mesh types or features are available in parallel. Degenerate heaxahedral cells, see Figure~\ref{fig:DegHex}, have been added to handle complex soil geometries, see Figure~\ref{fig:GeoMesh}. Such cells can be refined by the Arcane AMR procedure.

\begin{figure}
\includegraphics[width=1.\textwidth]{GeoscienceMeshCells.PNG}
\caption{Degenerate hexahedral cells occurring in geological meshes.}\label{fig:DegHex}
\end{figure}

\begin{figure}
\includegraphics[width=0.6\textwidth]{MeshGeologyBW.jpg}
\caption{An example of geological mesh.}\label{fig:GeoMesh}
\end{figure}

A mesh can be loaded by the platform under various formats (e.g.\ vtk, gmsh, med) or directly generated from a few input data. For instance, 2-D or 3-D Cartesian meshes as well as 2-D or 3-D honeycomb meshes can be created by just specifying the sizes of the domain and the number of cells.

For Cartesian meshes, Arcane provides a dedicated API giving access to neighbour elements in each direction. Cartesian meshes also come along with a dedicated block-refined AMR procedure.

Polyhedral meshes can currently be defined through the vtk format. This enables, for instance, the use of Voronoi meshes, see Figure~\ref{fig:Voronoi}. To handle such meshes, Arcane relies on the internal library \emphb{Neo}. This library is designed to handle complex polyhedral meshes that can possibly evolve during the simulation.\footnote{Evolving meshes can, for instance, occur in geosciences when simulating the evolution of a sedimentary basin during geological times.} One of its interesting features is the way it asynchronously handles the mesh evolution by means of a dependency graph.

\begin{figure}[!h]
\includegraphics[width=0.6\textwidth]{mar_precond_lbfgs.jpg}
\caption{An example of Voronoi mesh colored by cell volumes.}\label{fig:Voronoi}
\end{figure}

\subsection{A recent API for computing on accelerators}

Recently, a new API has been introduced in Arcane to handle accelerators with the following characteristics:
\begin{enumerate}[(i)]
\item it supports Cuda, HIP and SYCL backends;
\item it allows to use dynamically different types of GPU and multi-threaded CPU;
\item it ensures that the mesh connectivities and variables are available on the accelerator.
\end{enumerate}

\begin{figure}
\begin{minipage}[t]{.48\textwidth}
\begin{lstlisting}[language=customCpp, basicstyle=\scriptsize\ttfamily,caption=Cell loop classical CPU version., label={list:MeshVarCPUEx}]
// My mesh cell variables
Arcane::VariableCellReal m_a;
Arcane::VariableCellReal m_b;
Arcane::VariableCellReal m_c;

void compute_cpu_only()
{
 // given by Arcane
 Arcane::IMesh* mesh; 
















 // classical enumerate cell
 ENUMERATE_
    (Cell,icell,mesh->allCells()) {	!\label{line_28}!
  Real a = m_a(icell);
  Real b = m_b(icell);
  m_c(icell) = math::sqrt(a*a+b*b);	!\label{line_31}!
 };
}
\end{lstlisting}
\end{minipage}
\hfill
\begin{minipage}[t]{.48\textwidth}
\begin{lstlisting}[language=customCpp, basicstyle=\scriptsize\ttfamily,caption=Cell loop accelerator version., label={list:MeshVarGPUEx}]
// My mesh cell variables
Arcane::VariableCellReal m_a;
Arcane::VariableCellReal m_b;
Arcane::VariableCellReal m_c;

void compute_cpu_or_gpu()
{
 // given by Arcane
 Arcane::IMesh* mesh; 

 // define runner and queue
 // executing resource
 Arcane::Accelerator::Runner runner;		!\label{line_13}!
 // execution flow
 auto queue = makeQueue(runner);			!\label{line_15}!
 // computing kernel
 auto command = makeCommand(queue);		!\label{line_17}!
 
 // take for the device algorithm
 // 'a' and 'b' as inputs
 auto in_a = viewIn(command,m_a);		!\label{line_21}!
 auto in_b = viewIn(command,m_b);
 // 'c' as output
 auto out_c = viewOut(command,m_c);		!\label{line_24}!

 // Launch compute on device
 command << RUNCOMMAND_ENUMERATE
    (Cell,icell,mesh->allCells()) {
  Real a = in_a(icell);
  Real b = in_n(icell)
  out_c(icell) = math::sqrt(a*a+b*b);
 };
}
\end{lstlisting}
\end{minipage}
\end{figure}

The API was inspired by the SYCL API from Khronos group~\cite{sycl}. As SYCL and Kokkos~\cite{CarterEdwards20143202}, it handles different architectures (NVIDIA, AMD, Intel). Its runtime-based approach allows for instance to address CPU (sequential or using multi-threading) or GPU or both with a same executable. The Arcane Accelerator API aims at providing the basic concepts from Arcane on GPU. The API allows to write classical mesh-based loops on GPU with few differences from the CPU version. In Listings~\ref{list:MeshVarCPUEx} and~\ref{list:MeshVarGPUEx} we compare a classical mesh compute loop on CPU and its implementation with accelerator API. The differences between both codes come only from a dedicated setup needed on GPU, but the core of the loop itself is identical, see lines~\ref{line_28}--\ref{line_31}. In the GPU version, it is necessary to define a runner to specify the resource used for execution (line~\ref{line_13}), a queue on this runner to manage the execution flow (line~\ref{line_15}) and a computing kernel as a command (line~\ref{line_17}). To manage variable access on specific accelerator devices, this API uses specific views of variables with access rules: In, Out or Inout, see lines~\ref{line_21}--\ref{line_24}.

This API manages Arcane variables on items (Cell, Node...), and gives access to the mesh connectivity. It also handles specific multi-dimensional array of dimension one to four. To submit command and process loop on target device we use specific loop macro based on C++ 11 lambda functions and specific loop macro, see Listing~\ref{list:RuncommandEx}.

\begin{figure}[h!]
\begin{lstlisting}[language=customCpp,   caption={Arcane API to execute parallel loops on a target device.}, label={list:RuncommandEx}]
 // pseudo-code describing the structure of the RUNCOMMAND_ API
 command<<RUNCOMMAND_"LoopType"(args lambda consistent with "LoopType")
  { 
    /* body lambda code to process on target device */ 
  };

 // two examples of the RUNCOMMAND_ API
 
 // loop over mesh items (Cells, Nodes...). 
 // example for Nodes:
 RUNCOMMAND_ENUMERATE(Node,inode,mesh()->allNodes()){};
  
 // loop over multi-dimensional array of size one to four. 
 // example for one dimension:
 RUNCOMMAND_LOOP1(iter, nb_value){};
  
\end{lstlisting}
\end{figure}

Finally, it also provides some high level functionalities like reductions and algorithms such as filtering, partitioning or sorting. It can also automatically profile the accelerator loops.

\subsection{From Arcane to a simulation application}\label{subsec:ArcApp}

To build a numerical application upon the Arcane framework, the application developer will create \Modules{} to implement the different physical models simulated. The entry-points of a \Module{} structure the various computing steps of the model and are dynamically plugged into the application time-loop through a scheduler defined, by means of an xml \texttt{.config} file, by the application end-user.

Numerical and cross-functional tools are usually implemented in Arcane \Services{}, so that they can be shared between \Modules{} or even applications. These \Services{} may handle numerical schemes, mesh reader/writer, etc. The library Alien, for instance, offers a collection of Arcane \Services{} plugging different solver libraries, gathered under a unique interface.

Examples showing how to build an application upon Arcane are available in the platform repository~\cite{arcane-github} (tutorial or samples directories) or in mini applications arcane-benchs repository~\cite{arcanebenchs-github}. Two proxy-apps are also available in the Arcane Framework page on GitHub:
\begin{enumerate}[(i)]
\item the finite element proxy-app \emphb{ArcaneFEM}~\cite{arcanefem-github};
\item the porous media flow proxy-app \emphb{ShArc}~\cite{sharc-github}, cf.\ Sections~\ref{sec:ShArcEx} and~\ref{sec:GrandChallenge}.
\end{enumerate}

\section{The ArcNum framework numerical tools}\label{sec:numTool}

In this section, we present the \ArcNum{} framework itself: a set of software components provided to facilitate the work of application developers. These components are located at the interface between Arcane's kernel and the simulation application. These tools are for example plugged into the proxy-app ShArc~\cite{sharc-github}, used in Sections~\ref{sec:ShArcEx} and~\ref{sec:GrandChallenge} of this article.

The components of the framework \ArcNum{}, namely \emphb{\emph{GUMP}}, \emphb{\emph{Law Framework}} and \emphb{\emph{AuDi}}, are respectively dedicated to physical model handling, to the registration and evaluation of physical laws, and to automatic differentiation. They provide the needed tools to build the numerical core of an application in a highly readable and long-term maintainable way.

\subsection{GUMP data model}\label{sec:gump}

For a given physical problem, GUMP (Generated Unified Manager Properties) allows to define the corresponding data model in an xml-style format file. This xml model has a very simple grammar defined in the \texttt{gump.xsd} scheme file. This grammar makes it possible to define a hierarchical structure of entities as shown in Listing~\ref{lis:gumpentities} and, thus, to define all the physical properties carried by each entity of the model as in Listing~\ref{lis:gumpproperties}. When the program is compiled, this file is parsed to generate C++ code that can be used by the application's developer. Even if most of the properties are defined in this xml model file, it is still possible to add additional properties within the application code.

\begin{center}
\begin{lstlisting}[language=customXML, caption= Example of an instance of a GUMP tree of entities., label = lis:gumpentities]
<gump>
  <model namespace="ArcRes">
    <entity name="System">
      <contains>
        <entity name="FluidSubSystem" unique="true" /> 
      </contains>
    </entity>
    <entity name="FluidSubSystem">
      <contains>
        <entity name="FluidPhase" />
      </contains>
    </entity>
    <entity name="FluidPhase">
      <contains>
        <entity name="Species" />        
      </contains>
    </entity>
    <entity name="Species">
      <contains>
        <entity name="Component" />
      </contains>
    </entity>
  </model>
</gump>
\end{lstlisting}
\end{center}

\begin{center}
\begin{lstlisting}[language=customXML,caption= Example of GUMP properties associated to an entity., label = lis:gumpproperties]
 <entity name="FluidPhase">
  <supports>
    <property name="Viscosity" dim="scalar" type="real" />
    <property name="Density" dim="scalar" type="real" />
    <property name="CapillaryPressure" dim="scalar" type="real" />
  </supports>
</entity>
\end{lstlisting}
\end{center}

This model is instantiated at runtime. In general, this action is carried out using a dedicated Arcane \Service{}, see Section~\ref{subsec:Arcane} and~\cite{arcane09} for the description of Arcane \Services{}. This \Service{} allows the user to specify its model in the simulation input file, as shown in Listing~\ref{lis:arcgumpmodel}.

\begin{center}
\begin{lstlisting}[language=customXML, caption=Example of an input simulation file instantiating a GUMP model., label = lis:arcgumpmodel]
<physical-model>
  <system name="UserSystem">
    <name>System</name>
    <fluid-system>
        <name>Fluid</name>
        <fluid-phase>
            <name>Water</name>
            <species>H2O</species>
        </fluid-phase>
        <fluid-phase>
          <name>Gas</name>
          <species>CO2</species>
        </fluid-phase>
    </fluid-system>
  </system>
</physical-model>
\end{lstlisting}
\end{center}

This approach aims at bringing flexibility to flow models, for instance. Each property of the system is referred by a specific name called ``xpath''. It consists in specifying the property name with the path to the entity carrying it as a prefix. As an example, the input data of an initial condition \Service{} is printed in Listing~\ref{lis:init}. The pressure property is uniquely designed through the path \texttt{[System]System::Pressure} and, here, its value is set to \texttt{1.38E+7}\unit{Pa}.

\begin{center}
\begin{lstlisting}[language=customXML,caption=Example of an input simulation file where xpath properties are used for initial conditions., label = lis:init]
<initial-condition name="Constant">
  <property>[System]System::Pressure</property>
  <condition>
    <value>1.38E+7</value>
  </condition>
</initial-condition>
\end{lstlisting}
\end{center}

The generated C++ API also allows the application's developer to iterate over the model entities and manipulate the associated properties. An example of user code, parsing the physical model defined in Listing~\ref{lis:init}, is given in Listing~\ref{lis:DataModelParsing}. This code snippet prints the different entities of the physical model and shows how to access to the system pressure property by using the xpath we mentioned earlier.

\begin{center}
\begin{lstlisting}[language=customCpp,   caption= C++ API to parse the physical system and access to the properties., label = lis:DataModelParsing]
ENUMERATE_SUBSYSTEM(isubsystem, system.subSystems()) {
    info() << *isubsystem; // prints: Fluid
    ENUMERATE_PHASE(iphase, isubsystem->phases()) {
      info() << *iphase;  // prints: Water Gas
      ENUMERATE_SPECIES(ispecies, iphase->species()) {
        info() << *ispecies; // prints: H20 CO2
      }
    }
}
// Access System Pressure
auto system_pressure = ArcRes::XPath::property(system,"System::Pressure");
\end{lstlisting}
\end{center}

By defining a path to access to the different physical entities, GUMP provides a way to define the input and output properties of the law framework, described in the following section.

\subsection{Law framework}\label{sec:laws}

When solving systems of partial differential equations, we frequently have to evaluate physical laws in groups of mesh items (cells, nodes...). For instance, in the context of porous-media flow simulations, these laws may correspond to petrophysical empirical relationships such as the relative permeabilities or the capillary pressures. The proposed law framework allows users to easily define new laws for their application.

In our framework, a law is simply a function which calculates output properties from input ones defined on groups of mesh items. It also evaluates the derivatives of the output properties with respect to the input ones.

In the law API, properties can be scalar, multi-scalar or vector-based. The access to a scalar law is simply done using its property name. It is for instance the case of the porosity. A multi-scalar law can be retrieved using several indexes of properties. Thus, the access to the molar fraction of a component is done using the name of the component and the name of the fluid phase.

The definition of a new law is carried out in three steps. First, the user defines the input and output properties of the law in an xml-style file (see Listing~\ref{lis:sig}). This xml file is parsed at compilation time to generate a C++ header in the build directory. This header contains two classes, the \texttt{lawName::signature} defined in the xml file and the \texttt{lawName::law} class managing the evaluations of the law on different supports: mesh, mesh groups, item vectors or scalars. Second, the user has to implement the function which will be used to evaluate the law on a single mesh item (see Listing~\ref{lis:eval}).

\begin{center}
\begin{lstlisting}[language=customXML, caption=Example of xml file defining the law signature., label = lis:sig]
<law name="FluidDensityLawType">
  <input name="temperature" type="real" dim="scalar" />
  <output name="fluidDensity" type="real" dim="scalar" />
</law>
\end{lstlisting}
\end{center}

\begin{center}
\begin{lstlisting}[language=customCpp,  caption=User evaluation method corresponding to the xml definition., label = lis:eval]
void FluidDensityLaw::eval(const Real T, Real& rho, Real& drho_dt)const
{
    rho = T + exp(T);
    drho_dt = 1 + exp(T);
}
\end{lstlisting}
\end{center}

At last, to be able to use a new law in an input \texttt{.arc} file, an arcane \Service{} has to be created. This \Service{} allows the user to choose the input and output properties at runtime. Listing~\ref{lis:arcfilelaw} shows an example of \texttt{.arc} input file and Listing~\ref{lis:intanciate} is the implementation of the corresponding service. In this implementation, let us note that the two classes, generated in the header file, are instantiated: the \texttt{lawName::law} object taking in arguments the law evaluation function \texttt{FluidDensityLaw::eval} and the \texttt{law::signature} instance. At this stage the law is configured and added to the \texttt{function\_mng} object (see last line of Listing~\ref{lis:intanciate}). This last \texttt{Law::FunctionManager} object is used as a ``law-database''. When the application is initialized, each \Service{}, which instantiates a law, adds it to the \texttt{Law::FunctionManager} object. The API also offers a \texttt{Law::FunctionEvaluator} object which allows one to select a subset of the available laws. This object is used to evaluate the selected laws on a specific support, at each time step, for instance. Let us note that it is only at this stage that the evaluation support is specified. An example of law evaluation is given in Listing~\ref{lis:law_eval}.

\begin{center}
\begin{lstlisting}[language=customXML, caption=Example of input simulation file configuring the law input and output., label = lis:arcfilelaw]
   <law name="FluidDensityLawConfig">
      <output>
        <fluid-density>[Phase]Water::Density</fluid_density>
      </output>
      <input>
        <temperature>[System]System::Pressure</temperature>
      </input>
      <parameters>
        <initial-temperature>0</initial-temperature>
      </parameters>
    </law>
\end{lstlisting}
\end{center}

\begin{center}
\begin{lstlisting}[language=customCpp, basicstyle=\scriptsize\ttfamily, caption=Example of C++ code to create a law and register it in the function manager., label = lis:intanciate]
#include "FluidDensityLawType_law.h" // generated from xml file 
#include "FluidDensityLaw.h" // user definition of the law (cf. Listing !\ref{lis:eval}!)
void
FluidDensityLawConfigService::
configure(Law::FunctionManager& function_mng, ArcRes::System& system)
{
    // Initialize arguments from .arc input file
    auto fluid_density = options()->output().fluidDensity();
    auto temperature = options()->input().temperature();

    // Initialize parameters from .arc input file
    auto initial_temperature = options()->parameters().initialTemperature();
    
    // Set law arguments
    FluidDensityLawType::Signature signature;
    signature.fluidDensity = ArcRes::XPath::property(system, fluid_density);
    signature.temperature  = ArcRes::XPath::property(system, temperature);

    // Set law parameters
    FluidDensityLaw law;
    law.setParameters(initial_temperature);

    // Create law
    auto function = std::make_shared<FluidDensityLawType::Function>(
        signature, law, &FluidDensityLaw::eval);

    // Add law to function manager
    function_mng << function;
}
\end{lstlisting}
\end{center}

\begin{center}
\begin{lstlisting}[language=customCpp, basicstyle=\scriptsize\ttfamily, caption=Example of C++ code to evaluate the registered laws on a given support., label = lis:law_eval]
// Prepare law evaluation
Law::FunctionManager function_mng; // Build and filled by the application
Law::FunctionEvaluator function_evaluator(function_mng);

IVariableMng variable_mng; // Given by the application
auto variable_accessor = variable_mng.variableAccessor<Cell>(); 

// Evaluate laws registered in function_mng
function_evaluator.evaluate(accessor, allCells(), Law::eWithDerivative);

// Access to law outputs
auto system_pressure = ArcRes::XPath::property(system,"System::Pressure");
auto computed_pressure = accessor.values(system_pressure);
auto computed_pressure_derivatives = accessor.derivatives(system_pressure);
\end{lstlisting}
\end{center}

\goodbreak

To conclude, the workflow to define and evaluate a physical law consists of four steps:
\begin{enumerate}[(i)]
\item definition of the law inputs and outputs in an xml-style file (Listing~\ref{lis:sig});
\item computation of its values and derivatives through a user-defined function (Listing~\ref{lis:eval});
\item definition of a configuration \Service{} to instantiate and register the law object (Listing~\ref{lis:intanciate});
\item call to the law evaluator on a given support (Listing~\ref{lis:law_eval}).
\end{enumerate}

At runtime, users only have to complete the \texttt{.arc} file to set the actual inputs, outputs and parameters of the law, see Listing~\ref{lis:arcfilelaw}. After evaluating the laws, their values and derivatives can be read for any item included in the support defined in the evaluator. This functionality is a key point of the API since it provides all necessary data to build the matrices for the linear system assembly. These data are then gathered in \Contribution{} objects which are described in the sequel.

\subsection{Contributions}\label{sec:contributions}

The resolution of nonlinear problems such as $\mathbf{F}(\mathbf{u})=0$ is often carried out thanks to a Newton's method. Each iteration $k$ requires the construction of a linearized problem involving the Jacobian matrix of $\mathbf{F}$ and the residual vector $\mathbf{F}(\mathbf{u}^k)$. The \emphb{\Contribution{}} API provides a convenient way to assemble the Jacobian and the residual terms, to compute automatically their derivatives to be stored in the matrix, and to add all these values in Alien systems.

\subsubsection{Contributions --- Accessors}\label{subsec:contrib_accessor}

For any GUMP property, the \Contribution{} framework allows to recover its value and derivatives with respect to the main unknowns of the problem. When computing fluxes on a face, the derivatives often depend on various unknowns located around this face. This set of local elements is called a stencil. In the case of a diffusion flux discretized with the finite-volume two-point flux, this stencil includes the cells that share this face. To make the assembly of such fluxes easier in our framework, \Contributions{} can be accessed according to the stencil of a given face. Let us give an example. In the Listing~\ref{lis:contrib_accessor}, lines~\ref{lis:contrib_accessor_7}--\ref{lis:contrib_accessor_12} illustrate how to get \Contributions{} related to three GUMP properties which are \texttt{ArcRes::Pressure}, \texttt{ArcRes::CapillaryPressure} and \texttt{ArcRes::FluidDensity}. These \Contributions{} depend on an unknown manager that specifies the variables defined in each cell. This manager here appears twice since it is used with a two-point stencil. This stencil is initialized for each face on line~\ref{lis:contrib_accessor_17}. The following lines~\ref{lis:contrib_accessor_18}--\ref{lis:contrib_accessor_25} show how to get the values for one stencil cell (back or front to the face).

\begin{center}
\begin{lstlisting}[language=customCpp,  caption=Example of C++ code to access to the \Contributions{}., label = lis:contrib_accessor]
_buildFluxInternal(ArcNum::Vector& residual, ArcNum::Matrix& jacobian)
{
    // unknown manager of the system (List of unknown properties)
    const auto& um = unknownsManager();

    // contribution wrapper by property
    auto P = Law::contribution<ArcRes::Pressure>(domain(),				!\label{lis:contrib_accessor_7}!
             functionMng(),um,um,system());
    auto Pc = Law::contribution<ArcRes::CapillaryPressure>(domain(),
             functionMng(),um,um,system().fluidSubSystem().phases());
    auto Rho = Law::contribution<ArcRes::FluidDensity>(domain(),
             functionMng(),um,um,system().fluidSubSystem().phases());		!\label{lis:contrib_accessor_12}!

    Arcane::FaceGroup inner_faces = mesh()->allCells().innerFaceGroup();
    // Arcane loop on faces
    ENUMERATE_FACE(iface,inner_faces) {
        ArcNum::TwoPointsStencil stencil(iface);							!\label{lis:contrib_accessor_17}!
        const Law::Cell& cell_k = stencil.back();							!\label{lis:contrib_accessor_18}!
        const Law::Cell& cell_l = stencil.front();
        // GUMP loop on fluid phases
        ENUMERATE_PHASE(iphase, fluid.phases()) {
            // Acces contribution by mesh item [cell_k] or [cell_l] 
            const auto Pk = P[cell_k];
            const auto Pl = P[cell_l];
            const auto Pck = Pc[iphase][cell_k];							!\label{lis:contrib_accessor_25}!
            // idem for Pcl, Rhok and Rhol
            ...
        }
    }
}
\end{lstlisting}
\end{center}

\subsubsection{Contributions --- Combination and automatic differentiation with AuDi}\label{subsec:contrib_linop_audi}

New \Contributions{} can be created as a result of operations, for operators such as \texttt{*, /, +, -}, between several existing \Contributions{}. An illustration is given in Listing~\ref{lis:contrib_linop_audi}. Here, the \Contribution{} \texttt{grad\_kl} is computed from the \Contributions{} defined in the Listing~\ref{lis:contrib_accessor}. These operations make use of an automatic differentiation library, implemented in our framework and called AuDi.

\begin{center}
\begin{lstlisting}[language=customCpp, caption=Example of operations on \Contributions{} computing values and derivatives., label = lis:contrib_linop_audi]
_buildFluxInternal(ArcNum::Vector& residual, ArcNum::Matrix& jacobian)
{
    ...
        ENUMERATE_PHASE(iphase, fluid.phases()) {
            ...
            const auto grad_kl = ((Rhok + Rhol) / 2) * 
            (Pk + Pck - Pl - Pcl]);
            ...
            const auto flux_kl = ...// computed from grad_kl
        }
    }
}
\end{lstlisting}
\end{center}

AuDi is a C++ library based on generic programming techniques. It uses expression templates~\cite{ExprTemplate95,CppTemplates02} and operator overloading. Operator overloading is a classical C++ technique that enables to define and implement operators on any object. It is applied in our case to implement the operators \texttt{+,-,/,*} on contributions. On the other hand, expression templates prevent the runtime from creating unnecessary variables or temporary objects and thus offer good CPU and memory performance in critical calculation sections such as flux evaluations.

\subsubsection{Contributions --- Alien's linear system assembly}\label{subsec:contributions_assembly}

Finally, adding \Contributions{} to the Jacobian matrix or to the residual, defined by means of Alien's matrix and vector (see Section~\ref{subsec:Alien}), can simply be done by using the operators \texttt{+=} or \texttt{-=}. Listing~\ref{lis:Jac} shows an example where a finite-volume flux contribution \texttt{flux\_lk} is inserted into the lines related to the equation \texttt{iequation} of both stencil cells. Let us notice that the insertion of the derivatives into the corresponding matrix columns is here done in a very compact way thanks to the use of the \texttt{stencil} object.

\begin{center}
\begin{lstlisting}[language=customCpp,  caption=Example of C++ code to assemble a Jacobian matrix using the \Contributions{}., label = lis:Jac]
_buildFluxInternal(ArcNum::Vector& residual, ArcNum::Matrix& jacobian)
{
    ...
        ENUMERATE_PHASE(iphase, fluid.phases()) {
            ...
            const Arcane::Integer iequation = iphase.index();
            if (cell_k.isOwn()) {
                residual[iequation][cell_k]          += flux_kl;
                jacobian[iequation][cell_k][stencil] += flux_kl;
            }
            if (cell_l.isOwn()) {
                residual[iequation][cell_l]          -= flux_kl;
                jacobian[iequation][cell_l][stencil] -= flux_kl;
            }
        }
    }
}
\end{lstlisting}
\end{center}

\section{First examples with ShArc, an Arcane-based demonstrator}\label{sec:ShArcEx}

In this section, we present numerical illustrations derived from two porous-media models implemented in ShArc: a single-phase thermohaline convection model and an isothermal two-phase flow model. We first introduce the common notations and subsequently provide detailed descriptions of both models.

\subsection{Notations}\label{sec:notation}

The physical quantities used in both examples are introduced in this section. In the sequel, we also indicate their units where $M$ stands for mass, $L$ for length, $T$ for time and $\Theta$ for temperature. The numerical values are given in the SI system. Vectors in $\mathbb{R}^3$ are denoted with bold letters and tensors in $\mathbb{R}^3\times\mathbb{R}^3$ with blackboard-bold style letters.

We denote by $\phi [-]$ the porosity of the porous medium $\Omega$, $\perm [L^2]$ its permeability tensor. For a fluid phase $p$, $P_p [ML^{-1}T^{-2}]$ denotes its pressure, $\mu_p [ML^{-1}T^{-1}]$ the viscosity and $\rho_p [ML^{-3}]$ the mass density. We denote by $\vv_p$ Darcy's velocity which is defined by
\begin{equation}\label{eq:darcy.vel}
\vv_p(P_p) = -\frac{\perm}{\mu_p}\parens{\GRAD P_p - \rho_p \VEC{g}}
\end{equation}
where $\VEC{g} = -g\GRAD z [LT^{-2}]$ denotes the gravity vector.

\subsection{Thermohaline convection}\label{sec:saltDome}

For the first model, we consider the transport of salt by water through the porous medium $\Omega$. This problem, known as thermohaline convection~\cite{DIERSCH1998401}, can be modelled by the following partial differential equations:
\begin{itemize}
\item the water mass balance,
\begin{equation}\label{eq:mass.bal}
\phi \partial_t \rho_{\mathrm{w}} + \div \parens{\rho_{\mathrm{w}} \vv_{\mathrm{w}} } = 0;
\end{equation}
\item the salt mass balance,
\begin{equation}\label{eq:saltMass.bal}
\phi \partial_t C + \div \parens{\VEC{F}_C } = 0,
\end{equation}
where $\VEC{F}_C=C\vv_{\mathrm{w}} - \phi \diff \GRAD C $, $\diff [L^2T^{-1}]$ is the porous-medium molecular diffusivity, $C [ML^{-3}]$ the mass concentration;
\item the heat balance,
\begin{equation}\label{eq:heat.bal}
\partial_t \parens[\big]{(\phi c_{\mathrm{w}} \rho_{\mathrm{w}}+ (1-\phi) \rho_\mathrm{s} c_\mathrm{s})\theta } +
\div \parens{\VEC{F}_\theta } = 0
\end{equation}
where $\theta [\Theta]$ is the temperature,
\begin{equation}\label{eq:heatFlux}
\VEC{F}_\theta=\rho_{\mathrm{w}} \theta \vv_{\mathrm{w}} - \parens[\big]{(1-\phi) \diffT_\mathrm{s}+\phi \diffT_{\mathrm{w}}} \GRAD \theta
\end{equation}
and, for $p \in \{\mathrm{w}, \mathrm{s}\}$ (the index $s$ standing for salt), $c_p$ is the heat capacity $[L^2T^{-2}\Theta^{-1}]$ and $ \diffT_\mathrm{p} [MLT^{-3}\Theta^{-1}]$ the heat conductivity.
\end{itemize}
The boundary of $\Omega$ is partitioned into $\partial \Omega=\Gamma_\mathrm{top} \cup \Gamma_\mathrm{bottom} \cup \Gamma_\mathrm{vertical}$ where $\Gamma_\mathrm{top}$ is the highest part of the border, $\Gamma_\mathrm{bottom}$ the lowest one and $\Gamma_\mathrm{vertical}$ the remaining one. We denote by $\VEC{n}$ the unit normal vector outwards to $\Omega$. On the boundaries, we fix the following conditions:
\begin{align}
	\hspace{3cm} \vv_w \cdot \VEC{n} & = 0			&& \text{on $\partial \Omega$}, \label{eq:neumannFlow} \hspace{3cm}
\\	\hspace{3cm} \chi & = \chi_\mathrm{p} 			&& \text{on $\Gamma_p$, $\chi \in \{C, \theta\}$, $p \in \{\mathrm{top}, \mathrm{bottom}\}$}, \label{eq:dirichlet} \hspace{3cm}
\\	\hspace{3cm} \VEC{F}_\chi \cdot \VEC{n} & = 0		&& \text{on $\Gamma_\mathrm{vertical}$, $\chi \in \{C, \theta\}$}. \label{eq:neumannOther} \hspace{3cm}
\end{align}
Initial conditions are also defined by zero for $C$ and $\theta$.

\subsubsection{Numerical setting}

We consider a thermohaline convection problem inspired from a test case proposed in~\cite{DIERSCH1998401}. Here the 2D space domain $\Omega= (0, 300)\times(0, 300)$ represents a homogeneous and isotropic aquifer where
\begin{itemize}
\item $\theta_\mathrm{bottom} = 200$, $\theta_\mathrm{top} = 0$, $C_\mathrm{bottom} = 10$, $C_\mathrm{top} = 0$;
\item $\phi=0.1$, $\perm=1e-13$, $\diff=1e-13$;
\item $\rho_\mathrm{s} = 2700$, $c_\mathrm{s}=1180$, $\diffT_\mathrm{s}=2$;
\item $c_{\mathrm{w}}=4200$, $\lambda_{\mathrm{w}}=0.65$;
\item $\rho_{\mathrm{w}}$ is given by
\[
\rho_{\mathrm{w}}= \rho_{\mathrm{w},0}\parens[\big]{1-\beta_\theta (\theta - \theta_0) + \beta_C(C-C_0)},
\]
with $\rho_{\mathrm{w},0} = 1000$, $\beta_\theta=3.37e-4$, $\beta_c = 6.38e-4$, $\theta_0 = 0$, $C_0=0$, {$\mu_{\mathrm{w}}=1e-3$}.
\end{itemize}

A cell-centered finite volume discretization~\cite{Eymard.Gallouet.ea:00} is applied on the system~\eqref{eq:darcy.vel}--\eqref{eq:neumannOther}, with an implicit Euler time stepping, a two-point flux approximation (TPFA) for the diffusive terms and upwinding for the convective ones. We perform the simulation on a grid discretized with $100 \times 100$ uniform cells. The initial time step is $\tau^0=\SI{8.64e5}{s}$. At each time step we apply Newton's method~\cite{kelley2003solving} to solve the corresponding non-linear system. The linear system obtained at each Newton iteration is solved using an LU solver. If a divergence of the Newton algorithm occurs we divide the time step by 2. Otherwise, the time step is increased by multiplying it by 1.5. Figure~\ref{fig:saltDome} depicts the evolution of temperature and salt concentration at different time steps.

\begin{figure}[h!]
\includegraphics[height=4.83cm]{Temp28.jpeg}
\includegraphics[height=4.83cm]{Temp34.jpeg}
\includegraphics[height=4.83cm]{Temp40.jpeg}
\\
\includegraphics[height=4.83cm]{conc28.jpeg}
\includegraphics[height=4.83cm]{conc34.jpeg}
\includegraphics[height=4.83cm]{conc40.jpeg}
\caption{Evolution of the temperature (top) and salt concentration (bottom) after 28 years (left), 34 years (middle), and 40 years (right).}\label{fig:saltDome}
\end{figure}

\subsubsection{Implementation within ShArc}\label{sec:salt-NT}

This section illustrates how the numerical tools presented in Section~\ref{sec:numTool} are used to implement in ShArc the resolution of this thermohaline convection system.

Using GUMP (see Section~\ref{sec:gump}), temperature and salt concentration are declared into the physical model by means of an xml file. An example of a section from this file is shown in Listing~\ref{lis:TempConc}.

\begin{center}
\begin{lstlisting}[language=customXML, caption=Definition of salt concentration and temperature GUMP properties., label = lis:TempConc]
    <property name="Temperature" dim="scalar" type="real" />
    <property name="Concentration" dim="scalar" type="real" />
\end{lstlisting}
\end{center}

The law framework API, described in Section~\ref{sec:laws}, allows to create the different laws related to the flow problem which is here considered. In our example, the fluid density $\rho_{\mathrm{w}}$ is defined by the code given in Listing~\ref{lis:EvalDensity}. This piece of code will be evaluated in each mesh cell.

\begin{center}
\begin{lstlisting}[language=customCpp, caption=C++ code computing fluid density law., label = lis:EvalDensity]
void FluidDensityLaw::eval(const Real T, const Real C, Real& rho, 
                            Real& drho_dt, Real& drho_dc ) const
{
    Real factor = (1.0 - m_betat*(T - m_t0) + m_betac*(C - m_c0));
    rho = m_rho0 * factor;
    drho_dt = -1.0 * m_rho0 * m_betat;
    drho_dc = m_rho0 * m_betac;
}
\end{lstlisting}
\end{center}

The parameters of the law $\rho_{\mathrm{w}}$ are then specified within the input file like in Listing~\ref{lis:fluidDens}. Note that the previous generated properties here appear as inputs of the law.

\begin{center}
\begin{lstlisting}[language=customXML,  caption={Input simulation file defining fluid density law inputs, output and parameters.}, label = lis:fluidDens]
    <law name="FluidDensityLawConfig">
      <output>
        <fluid-density>[Phase]Water::FluidDensity</fluid-density>
      </output>
      <input>
        <concentration>[System]System::Concentration</concentration>
        <temperature>[System]System::Temperature</temperature>
      </input>
      <parameters>
        <rho0>1000</rho0>
        <betat>3.37e-4</betat>
        <betac>6.38e-4</betac>
        <t0>0</t0>
        <c0>0</c0>
      </parameters>
    </law>
\end{lstlisting}
\end{center}

When assembling the Jacobian matrix, contributions (see Section~\ref{sec:contributions}) provide the local values and derivatives of the variables and laws according to the primary unknowns of the system (here temperature and concentration) within the scheme stencil. Thus for a face $\sigma$ between two cells $K$ and $L$, the contributions contain the derivatives with respect to the unknowns of these two cells. Examples on how to get such contributions are given in Listing~\ref{lis:c1}.

\begin{center}
\begin{lstlisting}[language=customCpp, basicstyle=\scriptsize\ttfamily, caption=C++ code retrieving the contributions for the evaluation of the thermal flux., label = lis:c1]
auto T = Law::contribution<ArcRes::Temperature>(domain(), ...,system());

auto phi = Law::values<ArcRes::VolumeFraction>(domain(),fluid);
auto mu = Law::contribution<ArcRes::Viscosity>(domain(), ...,fluid.phases());

auto rho = Law::contribution<ArcRes::Density>(domain(), ...,fluid.phases());
auto rhof = Law::contribution<FluidDensity>(domain(), ...,fluid.phases());

auto lambda_f = Law::contribution<HeatConductivity>( ..., fluid.phases());
auto lambda_s = Law::contribution<HeatConductivity>( ..., solid.phases());
\end{lstlisting}
\end{center}

Finally, Listing~\ref{lis:c2} shows how to compute the discrete version of the thermal flux~\eqref{eq:heatFlux} and add it to the residual vector and to the Jacobian matrix.

\begin{center}
\begin{lstlisting}[language=customCpp,  caption=C++ code building the upwind TPFA heat flux~\eqref{eq:heatFlux} using the contributions., label = lis:c2]
const auto darcy_kl = (permeability / mu[fluid.phase(0)][cell_k])
        * tau[iface] * (P[cell_k] - P[cell_l] + rho_kl * dgz_kl);

const auto cond_s = ( (1.0-phi[K]) * L_s[solid.phase(0)][cell_k] +
                           phi[K]  * L_f[fluid.phase(0)][cell_k] );

const auto flux_t_kl =  darcy_kl
        * (one(darcy_kl >= 0)*T[cell_k]*rho[fluid.phase(0)][cell_k]
         + one(darcy_kl <  0)*T[cell_l]*rho[fluid.phase(0)][cell_l])
         + cond_s * tau[iface] * (T[cell_k] - T[cell_l]);

residual[iequation][cell_k]          += flux_t_kl;
jacobian[iequation][cell_k][stencil] += flux_t_kl;
\end{lstlisting}
\end{center}

\subsection{Two-phase flow model}\label{sec:TPH}

We now consider an isothermal immiscible two-phase model of water and gas~\cite{chen2007reservoir}. Each phase has only one component: H\textsubscript{2}O and CO\textsubscript{2}. Phase changes of these components are not here taken into account. Consequently, we use the indices ``$\mathrm{w}$'' and ``$\mathrm{g}$'' to refer indifferently to the corresponding phase or component. The difference between both phase pressures is equal to the capillary pressure:
\begin{equation}\label{pressure}
P_\mathrm{g} - P_{\mathrm{w}} = P_{\mathrm{c}_{\mathrm{g},\mathrm{w}}},
\end{equation}
where $P_{\mathrm{c}_{\mathrm{g},\mathrm{w}}} [ML^{-1}T^{-2}]$ is a function of the gas saturation. Consequently, the model unknowns are the water pressure $P_{\mathrm{w}}$ and the saturations $S_p$, $p \in \{\mathrm{g}, \mathrm{w}\}$. Writing the conservation of the volume of each phase leads to the following system:
\begin{align}
	\phi \partial_t(\rho_{\mathrm{w}} S_{\mathrm{w}}) + \div ({\rho_{\mathrm{w}}} k_{\mathrm{r},{\mathrm{w}}} \vv_{\mathrm{w}}) & = q_{\mathrm{w}}, \label{eq:waterBal}
\\	\phi \partial_t(\rho_\mathrm{g} S_\mathrm{g}) + \div ({\rho_\mathrm{g}}k_{\mathrm{r},\mathrm{g}} \vv_\mathrm{g}) & = q_\mathrm{g}, \label{eq:gasBal}
\\	S_{\mathrm{w}} + S_\mathrm{g} & = 1, \label{eq:porousVolCons}
\end{align}
where $k_{\mathrm{r},p}(S_p) [-]$ stands for the relative permeability and $q_\mathrm{p} [ML^{-3}T^{-1}]$ are sink or source terms which are non-zero where well perforations are present. No-flow boundary conditions are used for both phases. Initial conditions are given for $S_\mathrm{g}$.

For further details about the two-phase flow model and the more general multiphase compositional flows in porous media we refer to~\cite{Acs, aziz2002, Coats, Young}. More numerical experiments on practical problems illustrated on other Arcane-based IFPEN simulators with advanced resolution methods (adaptive mesh refinement, adaptive solvers, multilevel refinement, mixed multiscale finite elements) on general meshes can be found in~\cite{refId0,Di_Pi_Voh_Yous_a_post_comp_14, DIPIETRO20142331, SY2016, MALLIK2020112367, VOHRALIK2018728, YOUSEF2025116162, WOS:000446551700001}.

\subsubsection{Numerical setting}

We consider the injection of $\mathrm{CO_2}$ into a porous rock initially saturated with water. The spatial domain $\Omega =(0,1010)\times(0,1010)$ is discretized with a grid composed of $101\times 101$ uniform cells. The process is organized following a five-spot pattern, see Figure~\ref{fiveSpot}, where one injection well is located in the middle of the domain and producers are in the four corners (these producers are here only used to create gradients of pressure through the domain for the purpose of our test). Peaceman well index models~\cite{chen2007reservoir} are used to compute the discrete values of $q_{\mathrm{w}}$ and $q_\mathrm{g}$ where these wells are perforated (well radius is set to $0.5$ and the skin factor is zero). This injection-production process is simulated during $t_\mathrm{F}=10$ years with an initial time step $\tau^0=\SI{8.64e3}{s}$, which is equal to $0.1$~days. If one time step iteration is completed, the time step is increased by a factor of $1.5$. We divide it by $2$ if a divergence of the Newton algorithm occurs.

The injection pressure is set constant in time and equal to $P_{\mathrm{inj}}=5.6e7$. In the same way, the producer pressure is fixed to $P_{\mathrm{pro}}=2.7e7$.

\begin{figure}[!h]
\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[height=5.1cm]{5spots.pdf}
\caption{Five-spot pattern}\label{fiveSpot}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
\centering
\includegraphics[height=5.5cm]{perm5spots.PNG}
\caption{Permeability}\label{A-perm}
\end{subfigure}
\caption{Test-case configuration of Section~\ref{sec:TPH}.}
\end{figure}

\goodbreak

The problem parameters are chosen as follows:
\begin{itemize}
\item $\phi=0.1$, $\perm$ is given by a geostatistical realization shown in Figure~\ref{A-perm};
\item $\rho_{\mathrm{w}} = 1025$, $\mu_{\mathrm{w}}=1e-3$;
\item $\rho_\mathrm{g} = 1.89$, $\mu_{\mathrm{g}} = 1.42e-5$;
\item a Brooks--Corey model~\cite{Bro_Cor_por_64} was used for the petrophysical properties:
\begin{itemize}
\item the relative permeability of a phase $p \in \{\mathrm{w},\mathrm{g}\}$ is given by
\begin{equation}\label{eq:BCkr}
k_{\mathrm{r},p}(S_p) =
\begin{cases}
	1 & \text{if $S_p \geq 1$},
\\	\parens[\Big]{\frac{S_p-S_p^{\mathrm{res}}}{1-S_p^{\mathrm{res}}}}^n & \text{if $S_p^{\mathrm{res}} < S_p <1$},
\\	0 &	\text{if $S_p \leq S_p^{\mathrm{res}}$},
\end{cases}
\end{equation}
where the residual saturations are respectively given by $S_{\mathrm{w}}^{\mathrm{res}}=0.2$ and $S_{\mathrm{g}}^{\mathrm{res}}=0.1$ and $n=1$;
\item $P_{c_{\mathrm{g}, \mathrm{w}}}$ is given by
\[
P_{c_{\mathrm{{g,w}}}}(S_{\mathrm{g}})= P_{e} \cdot (S_e)^m,\qquad S_e=1-\frac{S_{\mathrm{g}}-S_{\mathrm{g}}^{\mathrm{res}}}{1-S_{\mathrm{w}}^{\mathrm{res}}-S_{\mathrm{g}}^{\mathrm{res}}},
\]
with $P_{e}=8.73e5$, $m=-\frac{1}{2.89}$.
\end{itemize}
\end{itemize}
As for the previous test case, a two-point finite volume scheme with upwinding in space and the backward Euler scheme in time are applied to discretize the system~\eqref{eq:darcy.vel}, \eqref{eq:waterBal}--\eqref{eq:porousVolCons}. Newton's method~\cite{kelley2003solving} along with LU solver is used to solve the resulting nonlinear system.

Figure~\ref{fig:sat} shows the evolution of the water saturation at three different time steps and, in particular, the faster front displacement towards the upper left and lower right corners where larger permeability values are present, see Figure~\ref{A-perm} too.

\begin{figure}[!h]
\includegraphics[height=4.83cm]{water_sat_0,5Y.jpeg}
\includegraphics[height=4.83cm]{water_sat_1Y.jpeg}
\includegraphics[height=4.83cm]{water_sat_1,5Y.jpeg}
\caption{Water saturation at $0.5$ year (left), $1$ year (middle), and $1.5$ year (right).}\label{fig:sat}
\end{figure}

\subsubsection{Implementation within ShArc}

As for the previous example, this section gives an illustration of the use of the GUMP data model and of the law framework to define the capillary pressure: Listing~\ref{lis:Pc} shows the line defining this property in GUMP, Listing~\ref{lis:PcCode} the C++ code defining the law, and Listing~\ref{lis:Pc-input} the definition of the parameters of this law in the input file. At last, Listings~\ref{lis:varTP} and~\ref{lis:fluxTP} detail the construction of the \Contributions{} required for the phase fluxes, their evaluation and their storage into the matrix and into the right-hand side of the linear system.

\begin{center}
\begin{lstlisting}[language=customXML, caption=Definition of capillary pressure GUMP property., label = lis:Pc]
<property name="CapillaryPressure" dim="scalar" type="real" />
\end{lstlisting}
\end{center}

\begin{center}
\begin{lstlisting}[language=customCpp,  caption=C++ code computing the capillary pressure law, label = lis:PcCode]
void CapillaryPressureLaw::eval(const Real S, Real& Pc, Real& dPc_dS)const
{
    Real Se = Arcane::math::min(1 - (S - m_Sr)/(1 - m_Sr_ref -m_Sr), 1.0);
    Real dSe_dS = - 1.0/(1 - m_Sr_ref -m_Sr);
    Real alpha = - 1./m_lambda;
    Pc = m_Pe * pow(Se, alpha);
    dPc_dS = alpha * m_Pe * dSe_dS * pow(Se, alpha - 1);
}
\end{lstlisting}
\end{center}

\begin{center}
\begin{lstlisting}[language=customXML, basicstyle=\scriptsize\ttfamily, caption= Input simulation file defining the capillary pressure law parameters., label = lis:Pc-input]
    <law name="CapillaryPressureLawConfig">
      <output>
        <capillary-pressure>[Phase]Gas::CapillaryPressure</capillary-pressure>
      </output>
      <input>
        <saturation>[Phase]Gas::Saturation</saturation>
      </input>
      <parameters>
        <Pe>8.73e5</Pe>
        <Sr-ref>0.2</Sr-ref>
        <Sr>0.1</Sr>
        <lambda>2.89</lambda>
      </parameters>
    </law>
\end{lstlisting}
\end{center}

\begin{center}
\begin{lstlisting}[language=customCpp, basicstyle=\scriptsize\ttfamily, caption=C++ code retrieving the \Contributions{} for the evaluation of the phase fluxes., label = lis:varTP]
const Arcane::VariableFaceReal& T = m_transmissivities["T"];
auto P  = Law::contribution<ArcRes::Pressure>(...,system());
auto Pc = Law::contribution<ArcRes::CapillaryPressure>(...,fluid.phases());

auto mu  = Law::contribution<ArcRes::Viscosity>(...,fluid.phases());
auto rho = Law::contribution<ArcRes::Density>(...,fluid.phases());
auto kr  = Law::contribution<ArcRes::RelativePermeability>(...,fluid.phases());
\end{lstlisting}
\end{center}

\begin{center}
\begin{lstlisting}[language=customCpp, caption=C++ code building an upwind TPFA phase flux using the \Contributions{}., label = lis:fluxTP]
const auto grad_kl = T[iface] * (P[cell_k] + Pc[iphase][cell_k]  
                               - P[cell_l] - Pc[iphase][cell_l]);

const auto mobility_k =
        rho[iphase][cell_k] * kr[iphase][cell_k] / mu[iphase][cell_k];

const auto mobility_l =
        rho[iphase][cell_l] * kr[iphase][cell_l] / mu[iphase][cell_l];

const auto flux_kl = (audi::value(grad_kl)>=0) ? mobility_k*grad_kl 
                                               : mobility_l*grad_kl ;

residual[iequation][cell_k]          += flux_kl;
jacobian[iequation][cell_k][stencil] += flux_kl;
\end{lstlisting}
\end{center}

\section{Massively parallel simulations on a simplified two-phase flow model}\label{sec:GrandChallenge}

\subsection{Tests on \emph{Irene} supercomputer}

We now present the first scalability tests performed on the Rome partition of \emph{Irene} supercomputer~\cite{ireneConf}. The Rome partition contains 2292~dual-processor AMD Rome (Epyc) nodes at \qty{2.6}{GHz} with 64~cores per processor, for a total of \num[group-separator = {,}]{293376} computing cores. All nodes are connected through a HDR-100 Infiniband network.

The physical model used for these tests is similar to the one presented in Section~\ref{sec:TPH}. As a first attempt, it has been deliberately simplified in order to first test the code on larger grids than the ones usually used with other Arcane applications. This simplified model differs in terms of boundary conditions, rock and fluid properties. Namely, here:
\begin{itemize}
\item $\Omega=(0,L)\times(0,L)\times(0,H)$ with $L=80000$, $H=1000$, $t_\mathrm{F}=100$ days;
\item the relative permeability of a phase $p \in \{\mathrm{w},\mathrm{g}\}$ is given by~\eqref{eq:BCkr} with $n=2$ and $S_{\mathrm{w}}^{\mathrm{res}}=S_{\mathrm{g}}^{\mathrm{res}}=0$, $P_{c_{\mathrm{g}, \mathrm{w}}}=0$;
\item $\phi=0.2$, $\perm=e^{f(x,y,z)}$ with
\[
f(x, y, z)=-30+3*\sin(2*\pi*x/L)*\sin(2*\pi*y/L)+y/L;
\]
\item $\rho_{\mathrm{w}} = \rho_\mathrm{g}= 1000$, $\mu_{\mathrm{w}}=\mu_{\mathrm{g}}=10^{-3}$;
\item $q_{\mathrm{w}}=q_\mathrm{g}=0$, no-flow condition are used on $\partial \Omega$ except on $x=0$ and $x=L$ where $\restr{P}{x=0} = 8.10^7$, $\restr{S_{\mathrm{g}}}{x=0} = 1$, $\restr{P}{x=L} = 10^7$.
\end{itemize}
The computation domain $\Omega$ is discretized with a Cartesian grid of size $N^2 \times N_z$ with $N=2048$ and $N_z=100$. The time step evolution was chosen according to the evolution of Newton's iterations starting with an initial time step of one day with a maximal value equal to $10$ days. The linear systems were solved with a BiCGStab method preconditioned with CprAMG (Hypre Boomer AMG + block Jacobi ILU0). Note that the defined grid and model lead to $\approx\num{1.258}$ billion equations at each Newton iteration. Starting from 10~nodes, MPI parallel simulations were carried out up to 128~nodes (\num[group-separator = {,}]{16384}~cores) without changing the grid size (strong scalability study). Figure~\ref{fig:SharcIrenePerf} shows the scalability of ShArc's time loop from 10 up to 128~nodes (blue curve) and the theoretical linear performance compared to the 10-node run (red curve) which is the lowest number of nodes that can run this test case. These results were obtained as part of the Genci \emph{Grands Challenges 2024}~\cite{GenciGC2024}. We observe a super-linear scalability of ShArc performances up to 100 nodes due to better cache usage as per core problem size decreases.
On this range of nodes, this study has highlighted a parallel efficiency of ShArc's time loop at around~90\%.

Secondly, simulations have also been carried out at the limit of the Alien platform's capacity:\footnote{Due to the equation indexing with 32-bit integers.} a grid containing 700 million of elements or so was distributed on 256~nodes (\num[group-separator = {,}]{32768} cores and \num{2.1}~billion equations on the whole). These tests performed on Irene supercomputer have also provided an opportunity to compare the performances of Alien's MCGSolver in both MPI and MPI-OpenMP contexts. Further details are given in~\cite{eccomas}. The input data file of these case can be found on GitHub~\cite{SharcTestCases}.

\begin{figure}[h!]
\begin{tikzpicture}[scale=0.7]
\begin{axis}[
title={Grid size: 2048x2048x100, 10 to 128 \emph{Irene} Rome nodes}, title style={align=center}, xlabel={Number of nodes}, ylabel={1/time (sec$^{-1})$}, scaled y ticks=false, y label style={at={(axis description cs:-0.1,.5)}}, xmin=0, xmax=128, ymin=0, ymax=0.007, xtick={0,16,32,48,64,80,96,112,128}, ytick={0,0.001,0.002,0.003,0.004,0.005,0.006,0.007}, grid=major, grid style=dashed, width=\textwidth, height=0.75\textwidth, legend pos=north west, legend style={font=\small}, every axis label/.append style={font=\small},]

\addplot[mpi, mark=square, mark size=1pt] coordinates { (10,0.000520569) (12,0.000694172) (14,0.000739051) (16,0.000957237) (18,0.00103677) (20,0.001137322) (22,0.001330325) (24,0.001394132) (26,0.001446602) (28,0.00160264) (30,0.001656639) (32,0.00189344) (34,0.001752571) (36,0.002037467) (38,0.002128687) (40,0.002247195) (42,0.002319285) (44,0.002458567) (46,0.00269147) (48,0.002710234) (50,0.002583496) (52,0.003076761) (54,0.003097416) (56,0.003115016) (58,0.003333013) (60,0.00346053) (62,0.003564333) (64,0.003668136) (66,0.003166343) (68,0.003497409) (70,0.003828475) (72,0.003949405) (74,0.00415538) (76,0.003987342) (78,0.004470736) (80,0.004475903) (82,0.00448773) (84,0.004298782) (86,0.004596886) (88,0.004574919) (90,0.004552952) (92,0.005024257) (94,0.005059699) (96,0.005095141) (98,0.005194374) (100,0.005353258) (102,0.005210309) (104,0.00534274) (106,0.004987437) (108,0.005562099) (110,0.005349658) (112,0.005538382) (114,0.005727106) (116,0.005758742) (118,0.00579307) (120,0.006357355) (122,0.006169829) (124,0.005982304) (126,0.006198564) (128,0.006254068) };

\addplot[linear, mark size=4pt] coordinates { (10,0.000520569) (12,0.000624682) (14,0.000728796) (16,0.00083291) (18,0.000937023) (20,0.001041137) (22,0.001145251) (24,0.001249364) (26,0.001353478) (28,0.001457592) (30,0.001561706) (32,0.001665819) (34,0.001769933) (36,0.001874047) (38,0.00197816) (40,0.002082274) (42,0.002186388) (44,0.002290502) (46,0.002394615) (48,0.002498729) (50,0.002602843) (52,0.002706956) (54,0.00281107) (56,0.002915184) (58,0.003019298) (60,0.003123411) (62,0.003227525) (64,0.003331639) (66,0.003435752) (68,0.003539866) (70,0.00364398) (72,0.003748093) (74,0.003852207) (76,0.003956321) (78,0.004060435) (80,0.004164548) (82,0.004268662) (84,0.004372776) (86,0.004476889) (88,0.004581003) (90,0.004685117) (92,0.004789231) (94,0.004893344) (96,0.004997458) (98,0.005101572) (100,0.005205685) (102,0.005309799) (104,0.005413913) (106,0.005518027) (108,0.00562214) (110,0.005726254) (112,0.005830368) (114,0.005934481) (116,0.006038595) (118,0.006142709) (120,0.006246822) (122,0.006350936) (124,0.00645505) (126,0.006559164) (128,0.006663277) };
\end{axis}
\end{tikzpicture}
\caption{Scalability of ShArc's time loop from 10 up to 128 nodes (blue curve) and the theoretical linear performance compared to the 10-node run (red curve).}\label{fig:SharcIrenePerf}
\end{figure}

\subsection{Tests on IFPEN \emph{Ener440} supercomputer}

In this second series of parallel tests, we take the model a step further introducing a capillary pressure $P_c$, given in the description below. The tests are executed on IFPEN \emph{Ener440} supercomputer, composed of 240 computing nodes with dual Intel Skylake-X G-6140 processors clocked at \qty{2.3}{GHz}. Each node holds $2 \times 18$ computing cores so this supercomputer provides a total of 8640 computing cores. All nodes are connected with an Intel Omni-Path high speed network.

The physical model used for these tests is close to the model used with \emph{Irene} supercomputer, but with a capillary pressure $P_c$. It has the following characteristics:
\begin{itemize}
\item $\Omega=(0,L)\times(0,L)\times(0,H)$ with $L=20000$, $H=100$, $t_\mathrm{F}=1000$ days;
\item $\rho_{\mathrm{w}} = \rho_\mathrm{g}= 1000$, $\mu_{\mathrm{w}}=10^{-3}, \mu_{\mathrm{g}}=10^{-4}$;
\item $q_{\mathrm{w}}=q_\mathrm{g}=0$, no-flow condition are used on $\partial \Omega$ except on $x=0$ and $x=L$ where $\restr{P}{x=0} = 6.10^7$, $\restr{S_{\mathrm{g}}}{x=0} = 1$, $\restr{P}{x=L} = 10^7$;
\item $P_{c_{\mathrm{g}, \mathrm{w}}}$ is given by
\[
P_{c_{\mathrm{g}, \mathrm{w}}}(S_{\mathrm{g}})= P_{e} \cdot (S_e)^m,\qquad S_e=1-\frac{S_{\mathrm{g}}-S_{\mathrm{g}}^{\mathrm{res}}}{1-S_{\mathrm{w}}^{\mathrm{res}}-S_{\mathrm{g}}^{\mathrm{res}}},
\]
with $P_{e}=10^5$, $m=-1$, $S_{\mathrm{w}}^{\mathrm{res}}=\num{0.2}$ and $S_{\mathrm{g}}^{\mathrm{res}}=\num{0.1}$.
\end{itemize}
The domain $\Omega$ is discretized with a Cartesian grid of size $N^2 \times N_z$ with $N=2048$ and $N_z=10$.

\begin{figure}[h!]
\begin{tikzpicture}[scale=0.7]
\begin{axis}[
title={Grid size: 2048x2048x10, 3 to 20 \emph{Ener440} nodes}, title style={align=center}, xlabel={Number of nodes}, ylabel={1/time (sec$^{-1})$}, scaled y ticks=false, y label style={at={(axis description cs:-0.1,.5)}}, xmin=0, xmax=20, ymin=0, ymax=0.00022, xtick={0,2,4,6,8,10,12,14,16,18,20}, ytick={0,0.00004,0.00008,0.00012,0.00016,0.00020,0.00024}, grid=major, grid style=dashed, width=\textwidth, height=0.75\textwidth, legend pos=north west, legend style={font=\small}, every axis label/.append style={font=\small},]

\addplot[mpi, mark=square, mark size=1pt] coordinates { (3,3.04618967423238E-05) (4,4.08545236459209E-05) (5,5.1527562460397E-05) (6,6.05096370801474E-05) (7,7.27030170549861E-05) (8,8.03958666748392E-05) (9,9.1788057747564E-05) (10,0.00010030974416) (11,0.000111475498194) (12,0.000121735434427) (13,0.000131323382749) (14,0.000139919279648) (15,0.000151296428137) (16,0.000156526559514) (17,0.000165525864231) (18,0.000177670250293) (19,0.000185263496103) (20,0.00019258042166) };

\addplot[linear, mark size=4pt] coordinates { (3,3.04618967423238E-05) (4,4.06158623230984E-05) (5,5.0769827903873E-05) (6,6.09237934846476E-05) (7,7.10777590654222E-05) (8,8.12317246461967E-05) (9,9.13856902269713E-05) (10,0.000101539655808) (11,0.000111693621389) (12,0.000121847586969) (13,0.00013200155255) (14,0.000142155518131) (15,0.000152309483712) (16,0.000162463449292) (17,0.000172617414873) (18,0.000182771380454) (19,0.000192925346035) (20,0.000203079311615) };
\end{axis}
\end{tikzpicture}
\caption{Scalability of ShArc's time loop from 3 up to 20 nodes (blue curve) and the theoretical linear performance compared to the 3-node run (red curve).}\label{fig:SharcEnerPerf}
\end{figure}

The performance results in Figure~\ref{fig:SharcEnerPerf} are similar to those observed with \emph{Irene} system but at a lower scale: $N_z$ is reduced from 100 to 10 and tests are performed up to 20 nodes (720 cores). Relative efficiency, compared to the 3 nodes test, is $>94\%$. The test case is also available at~\cite{SharcTestCases}.

\section{Conclusion}

The \ArcNum{} framework presented in this article offers easy-to-use numerical tools significantly simplifying the writing of numerical methods within scientific computing codes. These numerical tools are built upon the Arcane platform, providing mesh management, data structures and parallelism. For linear system resolutions, they rely on the Alien library, giving access to a wide range of linear solvers. Used together, \ArcNum{}, Arcane, and Alien offer a complete toolbox for building numerical simulation applications. This framework has been used to build the open source proxy-app ShArc for Geosciences simulation. The two case studies presented in this publication and conducted with ShArc highlight the flexibility, development time savings and robustness of the \ArcNum{} framework. A scalability study shows that this gain is not at the expense of parallel performance.

The \ArcNum{} framework is also used within complex industrial applications at IFPEN. For example, it is the foundation of a compositional multiphase reactive transport simulator, coupled with mechanics. In this context, more advanced schemes, such as (non)linear multipoint finite volume and virtual element methods, have been implemented~\cite{WOS:000415304900005, encheryAgelas2023}. In this much more demanding context, the framework has been used with more complex data models, a great number of physical laws and varied assembly algorithms. It was therefore necessary to add new functionalities to ensure the needed expressiveness and maintain the performance requirement. More precisely, the following features were added to the framework for internal use at IFPEN:
\begin{itemize}
\item the ability to create a graph of laws to enable the composition of physical laws;
\item the management of dynamic multi-point stencils;
\item a sparse automatic differentiation.
\end{itemize}

The next step for this framework will concern the physical law package, where the usage of AI trained models and the ability to compute asynchronously on accelerators will be added.

\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.

\printbibliography

\end{document}