\section{Introduction}
\label{sec:intro}

Texture synthesis, that is the process of synthesizing new image samples from a given examplar, has experienced a clear breakthrough with the work of Gatys et al. \cite{gatys2015texture}. Inheriting ideas from wavelet-based methods, it was first proposed in this work to synthesize textures by constraining second order statistics of the responses to classical classification neural networks. While this approach outperformed existing ones in term of visual fidelity, it relies on a relatively heavy optimization procedure that has to be carried for each new examplar. Therefore, methods have been proposed to perform synthesis using generative neural networks, such as networks learned for each new examplar \cite{ulyanovtexture} or GANs \cite{SGAN}. More recently, such approaches have been extended to perform the synthesis from a whole set of textures, by using adaptive normalizations of generative networks \cite{WCT} \cite{texton}.

The present work is in the continuation of these approaches and introduces an auto-encoder architecture enabling one to synthesize arbitrary textures from a compact latent representation. The autoencoder architecture avoids the training instability inherent to GAN-based approaches and is naturally suited to synthesis and editing tasks.
The encoding step yields a latent code adapted to an input texture, from which arbitrary samples can be synthesized. In order to do so, the latent space should be spatially agnostic, or in other words the network should treat in the same way different translations of a given input, which is made possible by carefully designing both the encoder and the generator.
Moreover, we wish the architecture to deal with both stochastic and structured, periodic-like textures. This is achieved thanks to the inclusion of sine waves in the design of the generator. Contrarily to what is done in the literature, we deal with these periodic components in a way that allows arbitrary rotations and scalings of the exemplar textures, a property that is clearly desirable in order to achieve a generic representation of textures. To the best of our knowledge, the resulting architecture is the first approach proposing a generic latent representation for textures, including an encoder to embed texture images, from which new samples can be synthesized or in which editing operations can be performed.

\section{Related Work}
\label{sec:related}
\subsection{Periodic Texture Synthesis}
Synthesizing periodic textures and textures with long range dependency is a challenging aspect of texture synthesis.
In order to enable the synthesis of multiple stationary periodic textures, the authors of \cite{PSGAN} improved the GAN-based synthesis method from \cite{SGAN} by providing additional long range periodic information to the generator. This concept has proven efficient and was reused in \cite{texton}, coupled with the structure of a StyleGAN2 \cite{StyleGAN2} network. Another approach is proposed in \cite{NeuralTextureGonthier} by imposing some constraints on the Fourier power spectrum.
\vspace{-1em}
\subsection{Universality and latent representations of textures}
Beyond the classical problem of exemplar-based texture synthesis, where the goal is to generate samples from a single input texture, a difficult problem is to develop a neural architecture that has the ability to synthesize arbitrary texture inputs, possibly not seen during the training phase. Incidentally, this also raises the question of how to efficiently encodes textures, e.g. through latent spaces. The authors of \cite{WCT} were the first to propose such an approach, leveraging the concept of adaptive instance normalization \cite{adain} to provide a flexible and universal style transfer and texture synthesis method, where each texture is represented by several thousands normalizing parameters related to VGG19 features. To allow for user control \cite{user-controllable} proposes a modification of PSGAN \cite{PSGAN}
by introducing a convolutional encoder. They focus on very small latent spaces
well suited to the representation of patches from a given image. Authors from \cite{neural} also use a convolutional encoder to generate textures, with exemplar-specific geometric transformations applied to the noise input of the generator.
Building on attention mechanisms, similarity maps or Fourier representations, some works \cite{Uatt}\cite{transposer} also tackle the problem of universal neural texture synthesis but suffer from verbatim copy of structures, such as those usually observed in patch-based approaches.
The method from \cite{texton} also has the ability to synthesize a variety of textures, but relies on a latent space which needs inverting the generator to represent given inputs.

\section{Method}
\label{sec:method}

\subsection{Architecture overview}
Our texture synthesis framework relies on an auto-encoder structure, that we call in short \textbf{TAE}. As explained in the introduction, we wish to develop a latent space that is spatially agnostic, deals with simple geometric transforms and with periodic textures.

The encoder $E$ is used to map an image to a latent texture representation space $\mathcal{W}$. From this representation, new samples with different spatial arrangements can be generated with the generator $G$. For this generator, we rely on the classical StyleGAN architecture \cite{StyleGAN2}. This architecture is augmented with spatially periodic information, namely sine maps as used in \cite{PSGAN}\cite{texton}. But unlike them, the periodic information is modulated using the latent representation and is thus adapted to each input sample. This adaptation of the periodic content is crucial and achieved by predicting scale and orientation parameters from the latent variable $w$.

The whole pipeline is mainly supervised by the texture reconstruction loss $\mathcal{L}_{style}$ \cite{gatys2015texture} between the input and its reconstruction by the texture auto-encoder.
Additionally, we introduce a network $T$ designed for direct texture sampling, enabling us to map a noise distribution onto our latent space $\mathcal{W}$ in a meaningful way. This allows the exploration of the latent space and possibly the creation of new textures. An overview of our method is shown in figure \ref{TAE}.

\vspace*{-.2cm} \begin{figure}[t!]
    \centering
      \includegraphics[angle=0,width=1.\linewidth]{figs/overview.pdf}
      \caption{Left: Texture Auto-encoder architecture. Right: Modified style block}\,

      \label{TAE}
      \vspace*{-.2cm}
\end{figure}

\subsection{Texture encoder}
In order to develop an encoder embedding texture images in the latent space of StyleGAN \cite{StyleGAN2}, we classically rely on the second order statistics of deep features, in the spirit of \cite{gatys2015texture} and the numerous works that have followed.  Given 5 different depths $l$, we retrieve features $F^l$ from the VGG19 network, each one having dimension $(C_l,H_l W_l)$, to extract their second order statistics. The main difficulty to build a texture encoder from these statistics lies in their high dimensionality. In order to extract a single low dimensional representation from these features, we first extract compact information from the VGG features at each scale (see Fig.\ref{quad}, right), which we then combine into a single compact vector $w=E(I)$ using an MLP.
At each scale, the second order features are extracted using the technique presented in \cite{second_order_cnn}. Indeed, for each depth $l$, instead of computing whole $(C,C)$ Gram matrices $G^l_{c,c'}=\sum_{k=1}^{H_lW_l} F^l_{c,k} F^l_{c',k}$ as in \cite{gatys} and then reduce their dimensionality, we directly compute from $G$ a representation vector $Q(G)$ of size $n_w$, the chosen dimension of our latent space. These vectors are defined as:
\begin{equation}
    Q(G)_i =m_i^T\cdot G\cdot m_i
\end{equation}
with $\{m_i\}_{1\leq i \leq n}$ a set of trainable vectors of size $(C,1)$. If we write $G=F^l \cdot F^{lT}$, the previous equation becomes:
\begin{equation}
    Q(G)_i =m_i^T\cdot F^l\cdot F^{lT}\cdot m_i =(m_i^T\cdot F^l)\cdot(m_i^T\cdot F^l)^T
\end{equation}
This allows the extraction of relevant second order features without having to compute the whole Gram matrix (as illustrated in Fig.\ref{quad}), reducing the complexity from $n(2HWC^2+4C) \approx 2nHWC^2$ to $n(2*CHW+2HW)  \approx 2nHWC$,
i.e. by a factor C, which goes up to 512 for the feature maps used in the latest depths of the VGG network.

\vspace*{-.2cm} \begin{figure}[]
    \centering
\includegraphics[width=.80\linewidth]{figs/encoder.pdf}
\caption{Left: Encoder, Right: Quadratic information extraction module}
      \label{quad}
      \vspace*{-.2cm} \end{figure} \vspace*{-.2cm}

This process results in a compact texture encoder, built only from second order statistics of some learned combination of features and carrying non localized  textural information. This non-localized nature is of course a desirable property in view of texture synthesis and will be shared by the generator.
\vspace{-0.7em}
\subsection{Texture generator}
We use a modified StyleGAN architecture as our generator G. To perform image synthesis, such an architecture uses a latent variable $w$ in a non-localized way (contrarily to DCGAN-like architectures): $w$ is used to predict and set the first and second moments of feature maps  during the multi-scale synthesis process (through Adaptive Instance Normalization \cite{adain}).

Aiming to remove any systematic spatial information in the generation process, we modify the StyleGAN architecture so that all spatial information is due to the realizations of noise maps at each step of the generation process: instead of learning a tensor of size (512*4*4), we discard spatial dimensions by only learning an input tensor (n*1*1), which we then expand by replication at the start of the generation process. Besides, convolutions use no padding to avoid border artifacts.

Thanks to these changes, our auto-encoder achieves full disentanglement of texture information and pattern localization.
Periodic information is also added at every level of the generator, as will be further described in sections \ref{sec:scale_indep}-\ref{sec:st}.
\vspace{-1em}
\subsection{Sine waves}
In order to synthesize periodic-like textures having long range dependency, we add sine waves to our generator, in a similar way to \cite{PSGAN}, \cite{texton}. Given a set of $n_{freq}$ frequencies $f'_i$  (detailed in \ref{sec:st}),
we build $S\in \mathbb{R}^{n_{freq}*H*W}$ a volume where each channel $S_i$ is a sine wave of frequency $f'_i$: $S=(sin(f'_i\cdot\boldsymbol{x}+\phi))_i,$ with $\boldsymbol{x} \in \mathbb{R}^{H*W}$  the spatial position, and $\phi$  a random phase.
In contrast to the works in \cite{PSGAN}, \cite{texton}, we modulate elements $S$ in the following ways:
\begin{itemize}
    \item we weight each channel $S_i$ of the sine maps $S$ with the coefficient $weight_{level}(f_i)$ defined in Section \ref{sec:scale_indep}. Its purpose is to inject the frequency $f_i$ in the right level of the generator given its magnitude.
    \item we use the latent variable $w$ to weight the use of every frequency according to the latent representation, using a fully connected layer $A \in \mathcal{M}_{n_{freq},n_w}$, in the spirit of \cite{StyleGAN2}, to project $w$ onto a weighting vector $A\cdot w$.
\end{itemize}
Finally, the addition of periodic content to a given feature map $F \in \mathbb{R}^{B*C*H*W}$ is performed by addition, after applying a 1 by 1 convolution filter to $S$ to reach the number of feature maps $C$ (see Fig.\ref{TAE}, right):
\begin{equation}
F \longleftarrow F+conv\Big(S *weight_{k_{level}} * A\cdot w)\Big)
\end{equation}
\vspace{-.8cm}
\subsection{Scale independent learnable frequencies}\label{sec:scale_indep}
In this section, we describe a novel way to generate periodic information independently from the architecture of the generator.
We call a level of the network the operations delimited by a change of resolution, in our case each level is a StyleGAN block, delimited by upsampling by a factor 2.

In the literature, periodic content is usually incorporated at specific levels of generators. Sine maps are either used at the lowest resolution level of the network in \cite{PSGAN} (not allowing for high-frequency content to be added), or at each level of the network as in \cite{texton} (leading to sudden discontinuity in the generation process due to the change of level). Instead of learning a different set of frequencies at each level, we choose to learn them independently of the architecture of the generator, and incorporate them in the generator at the adequate levels, depending on their magnitude. The idea is that a frequency $f$ has to be added at the level of the network where its magnitude (relative to the resolution of the level) is not too low nor too high.
The full procedure is illustrated in Fig.\ref{freqs}, left.
\begin{figure}[]
    \centering
      \includegraphics[width=1.1\linewidth]{figs/bias.pdf}
      \caption{\textbf{Left:} Use of a global frequency at different scales in the generator, depending on its magnitude. \textbf{Right:} Periodic information to incorporate to the network: learnable frequences are transformed using parameters depending on the input (through $w$).}
      \label{freqs}
      \vspace*{-.2cm} \end{figure}

In short, we learn a set of $n_{freq}$ frequencies $f_i$ independently of the architecture of the generator. Depending on the magnitude of a given frequency, we use it to create sine maps in the two levels of the generator where its magnitude is adequate (relative to the resolution of the level).
The use of a frequency $f$ is then weighted based on $|f|$: for each level number $k_{level}$ in the network , we define the weighting as:
\begin{equation}
weight_{k_{level}}(f)=\begin{cases} 1-\big|1-log_2(\frac{|f|}{f_0\cdot 2^{k_{level}}}) \big| \;\;\;\;\;\; \text{ if }\frac{f_0}{2}\leq \frac{|f|}{ 2^{k_{level}}} \leq  2f_0 \\
0 \text{ otherwise}
\end{cases}
\end{equation}
where $f_0$ corresponds to a reference magnitude.
This allows for a frequency's use level to smoothly move along the dyadic structure of the generator.

\subsection{Image specific scale and orientation estimation}\label{sec:st}
In this section, we describe a key component of our auto-encoder: each learned frequency $f_i$ is rotated and scaled accordingly to the input, both at training time and for inference. This allows to generate a new sample with appropriate periodic content, aligned with the input.

For each input image, scale and orientation parameters are learned from the corresponding latent variable $w$.
We directly infer these scale and rotation parameters  with an auxiliary network $Pred$: $Pred(w)=(s,\theta)$.
In order to define a scale and an orientation, one would usually need a reference, with scale $1$ and orientation $0$. We avoid this issue by forcing the network's predictions to be consistent with geometric transformations in the image space. This approach is reminiscent of self supervised learning (SSL) approaches.

Given an input image $I$, we geometrically transform $I$ with two distinct scaling and rotation parameters $\hat{s}_i$ and $\hat{\theta}_i$ and then encode these images to get latent codes:
$w_i=E(Scaling_{\hat{s}_i}(Rot_{\hat{\theta}_i}(I)))$. Finally, we predict transformation parameters from these latent codes: $s_i,\theta_i=Pred(w_i)$.
The corresponding loss $\mathcal{L}_{ssl}$ reads:
\begin{equation}
\begin{split}
    \mathcal{L}_{ssl}=1 - \Big|\frac{\pi}{2} -\big( (\hat{\theta}_1-\hat{\theta}_0) -(\theta_1-\theta_0)\big)mod_{\pi}\Big|
   +\big( \frac{\hat{s}_1}{\hat{s}_0}  - \frac{s_1}{s_0}\big)^2
\end{split}
\end{equation}

The first part of the loss takes into account the $\pi$-periodicity of the orientation of an image and is minimal when the difference in image rotation angles $\hat{\theta}_1-\hat{\theta}_0$ is equal to the difference of predicted angles $\theta_1-\theta_0 \;\; mod\; \pi$. Similarly, the second part of the equation is minimal when the ratio between estimated scales matches the ratio of the scales used in the geometrical augmentation. Eventually, these input-dependent parameters are used to scale and rotate the set of global frequencies $f_i$, effectively adapting the frequency patterns to each specific image, as illustrated in Fig.\ref{freqs}, right. Writing $s_{pred},\theta_{pred}=Pred(w)$, we get:
\begin{equation}
f'_i=Scaling_{s_{pred}}(Rot_{\theta_{pred}}(f_i))
\end{equation}

\vspace{-1em}
\subsection{Direct sampling}
Using our architecture, the natural way to synthesize a texture is to encode and decode an image into the latent space, performing synthesis by example. Now, it may be desirable to synthesize new textures without examplar input. In order to do so, we learn the distribution of the encoded training images $E(p_{data})$ in the $\mathcal{W}$ space. We choose an adversarial approach, where a network $T$ is trained to map a noise distribution $p_\mathcal{Z}$ to $E(p_{data})$. We train $T$ jointly with a discriminator $D_\mathcal{W}$, following the WGAN-GP framework \cite{WGANGP}, with a gradient penalty $\mathcal{L}_{GP}$:
\begin{equation}
    \mathcal{L}_{adv}(T)=\mathbb{E}_{z} \big[ D_\mathcal{W}(T(z)) \big]
\end{equation}
\begin{equation}
    \mathcal{L}_{adv}(D_\mathcal{W})=\mathbb{E}_{I\sim p_{data}} \big[ D_\mathcal{W}(E(I))\big]-\mathbb{E}_{z} \big[D_\mathcal{W}(T(z))\big]
    +\lambda_{GP} \mathcal{L}_{GP}
\end{equation}

\subsection{Losses}

We further constrain the texture synthesis process using the following losses. To compensate for low-frequency artifacts created by the optimization from \cite{gatys2015texture}, we add the spectral loss $\mathcal{L}_{Spe}$ from \cite{NeuralTextureGonthier} (equation 3.3) to force the spectrum of the generated image to match the input's spectrum.

Using sliced histogram matching \cite{sliced_color}, we also implement a color histogram loss:
\begin{equation}
\mathcal{L}_{Hist}(I_1,I_2)=\mathbb{E}_{x \in \mathbb{R}^3,\|v\|=1}\sum (sort(I_1 \cdot x^T)-sort(I_2\cdot x^T))^2
\end{equation}
Writing: $\mathcal{L}_{rec}(G,E)=\mathcal{L}_{style}(I,G(E(I)))$ our final optimization objective reads:
\begin{equation}
    \mathcal{L}=\mathcal{L}_{rec}+\mathcal{L}_{adv}+\mathcal{L}_{Hist}+\mathcal{L}_{Spe}
\end{equation}

\section{Experiments}

\label{sec:exp}
\subsection{Assessed methods, datasets and evaluation metrics}
Among the many texture synthesis algorithms, we choose to compare ourselves to these methods:
\begin{itemize}
    \item TextureCNN \cite{gatys2015texture}, as one of our goals is to shortcut this process with a feed-forward auto-encoder;
    \item PSGAN \cite{PSGAN}, as they introduced the idea of augmenting a generator with periodic content;
    \item Neural texture \cite{neural}, since they also use a texture encoder;
    \item Whitening Coloring Transform (WCT) \cite{WCT}, because this methods performs universal texture synthesis with no learning step.
\end{itemize}

We do not compare against \cite{user-controllable} (unavailable code) nor \cite{texton} as no encoder is included in their architecture. For all methods requiring a training (ours, PSGAN and neural texture), we choose a latent space dimension of 32.

PSGAN does not come with an encoder to project an image into a global latent space variable $z_g$. As a baseline, we perform inversion of PSGAN as in \cite{texture_mixer}: given an data sample, randomly initialize a latent variable $z_g$ and perform gradient descent in the latent space to minimize the texture loss $\mathcal{L}_{style}$ between the generated image $G(z_g)$ and the data sample. We keep the best result out of 10 random initializations. We call this resynthesis method PSGAN-GD.

We use the texture loss from \cite{gatys} as a distance to measure performance, along with SIFID (introduced in \cite{SinGAN}) and LPIPS \cite{LPIPS}. We use the dataset MacroTextures introduced in \cite{texton}.
We strongly augment each dataset with geometric transformations: rotations of any angles and scalings to train our method, Neural texture, and PSGAN.
We use the Adam optimizer with learning rate $1e^{-4}$ to train the whole pipeline end-to-end.
\vspace{-0.9em}
\subsection{Visual and quantitative results}
Synthesis results are simply obtained by first embedding an input exemplar into the latent space (using the encoder) and then generating a new sample from this latent code. Variability of the synthesized samples is classically obtained thanks to the noise maps in the generator. A comparison of syntheses obtained with the different methods is presented in Fig. \ref{comp_fig}. We observe that we indeed maintain appropriate long range dependencies, which is not the case in for TextureCNN nor WCT. We also get rid of low frequency artifacts mostly present in TextureCNN.

We show quantitative evaluations in Table 1. Our approach is very close to TextureCNN for the style metric, which leads to the best result as it directly optimizes $\mathcal{L}_{style}$ over the image.

Our method comes second best, notably beating WCT by a relatively big margin. Additionally, our method is orders of magnitude faster. We also obtain the best performance for the SIFID and LPIPS metrics, in particular outperforming the original TextureCNN approach. The results of PSGAN-GD are poor, which is essentially due to the instability of the inversion procedure in the latent space.

\newcommand{9}{9}
\newcommand{16}{16}
\newcommand{15}{15}
\newcommand{29}{29}
\newcommand{35}{35}
\newcommand{6}{6}
\newcommand{12}{12}

\vspace*{-.2cm} \begin{figure}[h!]
\centering
\begin{tabular}{cccccc}
\setlength{\tabcolsep}{.1pt}
\renewcommand{\arraystretch}{.1pt}
\centering
GT & TAE (ours) & WCT & TextureCNN & Neural Texture  \\

\includegraphics[width=.2\linewidth]{figs/test/GT/9.png}&
\includegraphics[width=.2\linewidth]{figs/test/TAE/9.png}&
\includegraphics[width=.2\linewidth]{figs/test/WCT/9.png}&
\includegraphics[width=.2\linewidth]{figs/test/gatys/9.png}&
\includegraphics[width=.2\linewidth]{figs/test/neural/9.png} \\

\includegraphics[width=.2\linewidth]{figs/test/GT/16.png}&
\includegraphics[width=.2\linewidth]{figs/test/TAE/16.png}&
\includegraphics[width=.2\linewidth]{figs/test/WCT/16.png}&
\includegraphics[width=.2\linewidth]{figs/test/gatys/16.png}&
\includegraphics[width=.2\linewidth]{figs/test/neural/16.png} \\

\includegraphics[width=.2\linewidth]{figs/test/GT/15.png}&
\includegraphics[width=.2\linewidth]{figs/test/TAE/15.png}&
\includegraphics[width=.2\linewidth]{figs/test/WCT/15.png}&
\includegraphics[width=.2\linewidth]{figs/test/gatys/15.png}&
\includegraphics[width=.2\linewidth]{figs/test/neural/15.png} \\

\includegraphics[width=.2\linewidth]{figs/test/GT/35.png}&
\includegraphics[width=.2\linewidth]{figs/test/TAE/35.png}&
\includegraphics[width=.2\linewidth]{figs/test/WCT/35.png}&
\includegraphics[width=.2\linewidth]{figs/test/gatys/35.png}&
\includegraphics[width=.2\linewidth]{figs/test/neural/35.png} \\

\end{tabular}
\caption{Comparison of syntheses obtained from 5 different methods.}
\label{comp_fig}
      \vspace*{-.2cm} \end{figure} \vspace*{-.2cm}

\begin{table}[h!]
\caption{Evaluation of the methods using $\mathcal{L}_{style}$, SIFID, LPIPS and runtime in seconds (all on one NVIDA RTX6000 GPU) on the MacroTextures dataset}
\centering
\scalebox{1}{
 \begin{tabular}{|l|c|c|c|c|}
 \hline
 & \;\;\;\;$\mathcal{L}_{style}$ ($\cdot 10^{3}$) \;\;\;\;&\;\; SIFID ($\cdot 10^{-6}$) \;\;&\;\;\;\;\; LPIPS \;\;\;\;\;&\;\;\;\; runtime\;\;\;\;\\
  \hline
  TextureCNN \cite{gatys2015texture}& $ \bold{3.8 \pm 16.8}  $ & $  14 \pm 14  $ & $  0.51 \pm 0.09 $ & $ 60$ \\ \hline

  TAE (ours)  & $ 5.2 \pm 11.7 $ & $ \bold{9 \pm  11}$ & $ \bold{0.35 \pm 0.09} $ & $ \bold{0.0072} $\\ \hline

  Neural texture \cite{neural}& $ 135 \pm 468 $ & $ 178 \pm 190 $ & $ 0.56 \pm 0.09 $ &  $0.085$  \\ \hline
  WCT \cite{WCT} & $ 29.8 \pm 132.9 $ & $ 10 \pm 18 $ & $ 0.44 \pm 0.07 $ & $ 3.5 $ \\ \hline

  PSGAN-GD \cite{PSGAN}& $ 208 \pm 611 $ & $ 381\pm 607 $ & $ 0.64 \pm 0.11 $ & $ 24$\\ \hline
 \end{tabular}}
 \label{table1}
 \vspace*{-.4cm}
\end{table}

\subsection{Geometric completeness}
In order to illustrate the ability of the method to deal with simple geometric transforms, we perform auto-encoding on various scalings and rotations of the same input image, thus showing successful adaptation to the input in Fig.\ref{augmentation}.

\begin{figure}[h!]
    \centering
      \includegraphics[width=.90\linewidth]{figs/augmentation_3.png}
      \caption{Top row: texture image geometrically augmented, Bottom row: feed forward texture resynthesis $G(E(I))$}
      \label{augmentation}       \vspace*{-.2cm} \end{figure}
      \begin{figure}[h!]
    \centering
    \includegraphics[width=.45\linewidth]{figs/fig_capture.png}
\caption{Illustration of the difficulty of PSGAN to deal with arbitrary rotations \newline \textbf{First row}: Rotations of input texture $R_\theta(I)$, $I=G(z_g)$ - \textbf{Second row}: PSGAN samples where $z_{g,opt}$ is obtained by gradient descent to match the texture $R_\theta(I)$ - \textbf{Third row}: Directly rotating the frequencies without changing $z_{g}$}
    \label{PSGAN_rot}
      \vspace*{-.2cm} \end{figure}

As a comparison, we show that PSGAN and thus the methods based upon it (e.g. \cite{user-controllable}) lack this ability.  To exhibit that PSGAN indeed does not have the capacity of handling different rotations of the same textures, we start by sampling a global variable $z_g$, generating an image $I=G(z_g)$.

Then, we create different rotations of $I$: $R_\theta(I)$ as shown in the first row of Fig.\ref{PSGAN_rot}, and try to reproduce them with the generator $G$. Indeed, taking $I=G(z_g)$ avoids inversion problems. On the contrary, as observed in Fig. 6a) of \cite{user-controllable}, the latent space of PSGAN is not smooth, thus the inversion of this network is
unstable.

The second row shows $G(z_{g,opt})$, where $z_{g,opt}$ is obtained by minimizing  $\mathcal{L}_{style}(R_\theta(I),G(z_{g,opt}))$ gradient descent. We notice that the periodic structure is not rotated, while smaller details are correctly oriented. This is indeed not a satisfactory result. In the PSGAN architecture, the frequencies of sine maps are inferred from $z_g$. We also tried to directly rotate these frequencies with the same angle $\theta$ as we rotated the image $I$ with, thus aligning the periodic content with the image we want to reconstruct, $R_\theta(I)$. This operation yields poor pattern fidelity albeit somewhat aligned, as can be seen on the third row.

These two tests prove that, although trained on a dataset augmented with rotations, locally in the latent space of PSGAN, each texture is highly linked to an orientation. There is no way to change the orientation without harming texture fidelity.
\captionsetup[subfloat]{labelformat=empty}

\vspace{-1em}
\subsection{Spatial interpolation}
\begin{figure}[h!]
    \centering
    \includegraphics[width=.46\linewidth]{figs/inference/nuancier.png}\hfill
      \includegraphics[width=.18\linewidth]{figs/interpolation_0.png}\includegraphics[width=.36\linewidth]{figs/interpolation_zoom_0.png}
      \caption{\textbf{Left:} Texture Palette. Four input textures of various sizes  and aspect ratio are displayed in the corners. In the center, an image shows a field obtained by performing spatial interpolation between the corner textures. On this field, the red dots are positions from which the images in the middle of each side are synthesized.
      \newline \textbf{Right:} Interpolation between 4 textures without and with expansion by a factor 2}
      \label{interpolation}
      \vspace*{-.2cm} \end{figure} \vspace*{-.2cm} \
Our approach relies on a texture latent space and is therefore naturally adapted to texture interpolation, for instance by simply linearly interpolating latent variables. Spatial interpolation (building an image across which a texture is progressively varying), on the other hand, is more challenging. Indeed, AdaIN layers normalize features based on statistics computed across the whole image, therefore globally defining new features. We solve this problem by computing statistics locally. Additionally, the mean and variance imposed after normalization are computed from a texture representation $w$ varying in the spatial domain, to allow a smooth transition of the content. An asset of this method is that the periodic content has a natural variation in space. We can also use it to create interpolated texture from a visual palette (Fig.\ref{interpolation}, left) We show an example of texture interpolation, and interpolation combined with expansion in Fig.\ref{interpolation}, right.

\section{Conclusion}
\label{sec:conclusion}
We introduced a feed forward auto-encoder network having the ability to learn a set of textures, yielding fast exemplar-based synthesis results.
Thanks to adaptive periodic content, our method allows texture encoding and reconstruction of any orientation and scale. Its disentanglement of texture characteristics and spatial distribution allow for common texture manipulations such as interpolation and expansion. The proposed approach outperforms closely related methods for usual texture fidelity metrics.

\section{}


\section{Introduction}

\begin{figure*}[htbp]
\begin{center}
\centerline{\includegraphics[width=0.95\linewidth]{Figure1-Schematic-Updated-Revised.png}}
\caption{PINNs learning the solution of 2D N-S equations in a complex geometry (wavy channel flow problem, $Re=100$). Our LSA-PINN method can learn accurate solution with fewer (50\% less) training samples and faster speed (50\% less training iterations).}
\label{fig:fast-pinn-advantage}
\end{center}
\vspace*{-5mm}
\end{figure*}

Physics-informed neural networks (PINNs) have emerged as a promising method for learning the solution of dynamical system from the governing physics~\cite{raissi2019physics}. PINNs have recently been studied for a wide range of physical phenomena and applications across science and engineering domains — electromagnetic, fluid dynamics, heat transfer, etc~\cite{karniadakis2021physics,cuomo2022scientific}. The distinctive feature of PINNs is the use of governing physics law, typically in the form of partial differential equations (PDEs), as the learning objective. This physics-informed learning constrains the PINN from violating the underlying physics at all training points sampled from the problem domain.

Existing PINNs evaluate the PDE constraints in their training loss by either automatic differentiation (AD) or numerical differentiation (ND)-type method. \textit{AD-loss} has become a popular choice because it can exactly compute the derivative terms at any sample location directly from the input and the network weights. While both methods have their pros and cons, \textit{ND-loss} can be flexibly implemented across many different neural network (NN) architectures, including both multi-layer perceptrons (MLPs) and convolutional neural networks (CNNs), because they do not require the NN to retain differentiability, unlike AD. Recent studies ~\cite{gao2021phygeonet,fang2021high,chiu2022can} have suggested that ND-type methods and especially \textit{coupled-automatic-numerical
differentiation (CAN)-loss} ~\cite{chiu2022can} can more robustly and efficiently produce accurate solutions with fewer training samples, whereas conventional \textit{AD-loss} is prone to failure during training. This is because ND-type methods approximate high order derivatives using \emph{PINN output from neighbouring samples}, hence, they can effectively connect sparse samples into piecewise regions via these local approximations, thereby facilitating fast physics-informed learning across the entire domain with sparser samples.

When dealing with irregular geometries, it becomes difficult for existing PINNs to perfectly connect training samples in the domain’s interior to the boundary. Failing to do so can cause undesirable training failure as the PINN starts to overfit at the near boundary region. Since many PDEs of practical interest are boundary-value problems, it is desirable to have the PINNs model the correct boundary behaviors. While adding dense samples to better refine the piecewise local regions near the boundary may improve accuracy, the extent to which sampling needs to be increased is empirical, and a denser sampling strategy may adversely impede training convergence. Moreover, it may not be feasible to locally refine the near boundary sampling density for PINNs with CNN-architecture.

Hence, we propose a loss function formulation in this work that helps provide linear approximation to the local gradient behaviour at the boundary, thereby restoring connectivity between domain boundary and near boundary interior samples. Such \textit{local structure approximation} also provides certain structural bias during PINN training to prevent inaccurate representation near the boundaries. This new boundary connectivity (\textit{BCXN})-loss function is key to a novel class of LSA-PINN method which can more efficiently learn the solution to PDEs with sparser training samples, regardless of domain geometry; see example in Fig.~\ref{fig:fast-pinn-advantage}. In addition, this method can be jointly implemented with other PINN advances in optimization such as loss balancing, domain decomposition and adaptive sampling.

In the rest of this work, we present comprehensive experiments to demonstrate i) the flexible implementation of LSA-PINN method for both MLP and CNN architectures; and ii) the effectiveness of LSA-PINN for learning multiple complex fluid dynamical systems, including for meta-modelling. Compared to conventional PINNs with the \textit{AD-} and \textit{ND-loss}, our LSA-PINNs with linear \textit{BCXN-loss} are shown to be capable of tackling challenging PDE problems while using sparser training samples, hence expanding the exciting potential of PINNs for learning complex, real-world dynamics.

\section{Related Work}

\textbf{Efficient sampling in PINNs.} The theoretical limit of physics-informed loss learning in relation to training samples has been provided by prior studies ~\cite{lu2021machine,mishra2022estimates}. To improve PINN training speed for practical applications ~\cite{markidis2021old}, several studies have focused on efficient sampling strategies such as importance sampling, adaptive sampling, and sequential sampling to reduce the amount of training samples required during PINN training ~\cite{anitescu2019artificial,mcclenny2020self,wight2020solving,nabian2021efficient,lu2021deepxde,lye2021iterative,daw2022rethinking,mattey2022novel, wu2023comprehensive}. Domain decomposition and parallelization strategies have also been explored to speed up training ~\cite{jagtap2020conservative,jagtap2021extended,shukla2021parallel,dong2021local,li2019d3m}. Our LSA-PINN differs from these works in that we make physics-informed learning more robust in the sparse sample regime via a newly-proposed \textit{BCXN-loss}.

\textbf{CNN architecture and numerical differentiation (ND)-type loss for PINNs.} CNN-based formulation allows us to design and train larger, more powerful networks, hence it has potential to be scalable for more complex, large-scale PDE problem~\cite{wandel2020learning,gao2021phygeonet,wandel2021teaching,ranade2021discretizationnet,wandel2022spline,ren2022phycrnet}. However, cumbersome coordinate transformations may need to be performed to better handle irregular domains. CNN-architecture PINNs usually utilize ND methods for computing the PDE loss on their input grid, which is much cheaper to compute as compared to AD-based loss ~\cite{gao2021phygeonet}. Hybrid frameworks that couple both AD- and ND-type loss have also been proposed for MLP architecture PINNs to unify the advantages of both methods ~\cite{jagtap2020conservative,mitusch2021hybrid,fang2021high,chiu2022can}. Our method further augments the efficiency and applicability of ND-type loss on irregular geometries.

\textbf{Enforcing BC constraint in PINN loss.} It is important for the PINNs to prioritize learning the correct boundary behaviors~\cite{shin2020convergence, wang2022respecting}. The BCs are usually enforced into the PINN loss as a soft constraint via penalty method. Various strategies have been explored in this context to dynamically calibrate the relative important between PDE and BC constraints during PINN training~\cite{elhamod2020cophy,bischof2021multi,jin2021nsfnets,thanasutives2021adversarial,maddu2022inverse,de2022mopinns,van2022optimally,xiang2022self,huanguniversal, wang2021understanding}. There are also other approaches to bypass the loss balancing issue. For example, one can devise an ansatz function such that the BCs are exactly satisfied by construction or implicitly formulate the BC constraints into the PDE loss~\cite{lagaris1998artificial,lagaris2000neural,mcfall2009artificial,berg2018unified,nabian2019deep,karumuri2020simulator,wang2020mesh}, leaving only single loss term from PDE residuals to be optimized. Another example is the use of augmented Lagrangian method to impose the BC constraints into PINN loss~\cite{lu2021physics}. Our method implicitly incorporates BCs into the PDE loss via a different strategy.

\section{Preliminary}

\subsection{Governing physics - incompressible Navier-Stokes (N-S) equations}

Fundamentally, PINNs have been shown to be applicable to many different physical systems. The present study focuses on learning fluid dynamics with PINNs. We consider fluid problems where the governing physics are the steady-state, incompressible N-S equations derived from the conservation of mass and momentum:
\begin{subequations}
\label{eq:ns_eqn}
\begin{align}
&\nabla\cdot\vec{u} = 0 \\
&\left(\vec{u}\cdot\nabla\right)\vec{u} = Re^{-1}\Delta\vec{u} - \nabla p
\end{align}
\end{subequations}
In the above PDEs, the primitive variables $\vec{u}$ and $p$ are velocity vector and pressure while Reynolds number ($Re$) represents the ratio of inertial to viscous forces. $\vec{u}$ consists of 2 components $\left(u,v\right)$ for a 2D case. While this work uses Cartesian coordinates, this formulation is extendable to other coordinate systems.

\subsection{MLP- and CNN-architecture PINNs}

We use a fully connected DNN architecture (e.g., MLP, CNN) to represent the solution of the dynamical process \textbf{\textit{U}} $= [\vec{u},p]^{T}$. For MLPs, the input $\vec{x}$ is a point coordinate in spatial domain. For CNNs, the input \textbf{\textit{X}} is a tensor with a fixed shape, comprising the entire (discretized) spatial domain. The accuracy of the PINN outputs \textbf{\textit{U}}($\vec{x}$;\textbf{\textit{w}}) given input $\vec{x}$ is determined by the network parameters \textbf{\textit{w}}, which are optimized w.r.t. the PINN loss function during training. The PINN loss function is defined as the composition of a PDE loss component $\left(L_{PDE}\right)$ and a BC loss component $\left(L_{BC}\right)$:
\begin{subequations}
\label{eq:loss_fn}
\begin{align}
L_{PINN} &= \lambda_{PDE} L_{PDE} + \lambda_{BC} L_{BC}  \\
\begin{split}
L_{PDE} &= \big\|\nabla\cdot\vec{u}(\vec{x};\textbf{\textit{w}})\big\|_{\Omega}^2 + \big\|\left(\vec{u}(\vec{x};\textbf{\textit{w}})\cdot\nabla\right)\vec{u}(\vec{x};\textbf{\textit{w}}) \\
& - Re^{-1}\Delta\vec{u}(\vec{x};\textbf{\textit{w}}) + \nabla p(\vec{x};\textbf{\textit{w}}) \big\|_{\Omega}^2 \end{split} \\
L_{BC} &= \big\|B[u(\vec{x};\textbf{\textit{w}})] - \vec{u}(\vec{x}) \big\|_{\partial\Omega}^2
\end{align}
\end{subequations}
The PDE loss penalizes deviation from governing N-S equations for the PINN output \textbf{\textit{U}}($\vec{x}$;\textbf{\textit{w}}) over the fluid domain $\vec{x}\in\Omega$, whereas the BC loss penalizes deviation from the desired Dirichlet boundary condition $\vec{u}_{BC}(\vec{x})$ which constrains the velocity field at the domain boundary $\vec{x}\in\partial\Omega$. The relative weights $\lambda_{PDE}$ and $\lambda_{BC}$ in Eq.(\ref{eq:loss_fn})a control the trade-off between different loss components.

The PDE and BC loss components are defined over a continuous domain, but for practical reasons, we evaluate the discretized PINN loss from a finite set of $n$ samples $D=\{( \vec{x}_{i} )\}_{i=1}^{n}$ during training. The CNN-architecture PINNs naturally acquire an input with equidistantly spaced grid where all the samples are coming from. For the MLP-architecture PINNs, it is also convenient to design training samples based on equidistantly spaced grid, then $m (< n)$ samples can be randomly drawn to compute the PINN loss during each stochastic gradient descent (SGD) mini-batch training iteration.

Although it is ideal for a PINN to satisfy the PDE constraint on a very dense set of samples, i.e., $\Delta\vec{x} \rightarrow0$ where $\Delta\vec{x}$ is the distance between two adjacent sample points, training such PINN may not be computationally feasible even with state-of-the-art SGD variants. For a more complex physics problem, training on dense samples is correspondingly more computationally expensive, and is also more likely to become unsuccessful due to adverse impact by some local region with very complex gradient behaviours, e.g., abnormally high gradient~\cite{fuks2020limitations,michoski2020solving,ramabathiran2021spinn,ji2021stiff,lucor2021physics,huanguniversal}. Given the same problem, training PINNs on sparser samples means fewer restrictions (the PDE constraint that needs to be satisfied) and an increased likelihood of avoiding problematic gradient points, thereby accelerating training. The \emph{relaxation of PDE constraint via reducing the sample density makes learning physically meaningful solution feasible}, although it may introduce some approximation error to the solution. Note that there is still an open question as to how best balance sample density and solution accuracy. Methods that can consistently produce more accurate solutions with sparse training samples like LSA-PINN are hugely advantageous for learning complex, real-world dynamics.

\subsection{PINN learning challenge on sparse samples}

When samples are sparse relative to the local complexity, a highly flexible, over-parameterized PINN can become “overfit”, leading to incorrect solution even as the PDE constraint is fulfilled at all sample points (i.e., the convergence is good). For \textit{AD-loss} PINN, this phenomenon can occur at any location within the computation domain. On the contrary, \textit{ND-loss} are more robust to sample density, because ND-type methods approximate the derivative terms by using \emph{output from neighbouring samples}, e.g., a finite difference-type stencil, which connects sparse samples into piecewise continuous regions to facilitate fast training across the entire problem domain. The use of lower-order \textit{local structure approximation} for the derivative terms also provides certain structural bias for PINNs to reduce overfitting. For example, abnormally high gradients are regularized so that it is easier for PINNs to learn the solution, with the accuracy being dependent on the approximation employed.

However, \textit{ND-loss} can still exhibit training failure if the near-boundary training samples do not perfectly connect to the domain boundary through their stencils (as illustrated in Fig.~\ref{fig:stencil-schematic}). Their outputs essentially act as free parameters during optimization for satisfying the PDE constraint at near-boundary samples, and the PINN training tends to “overfit” the near-boundary region to more easily jointly satisfy the PDE and BC constraints, leading to incorrect solution. \emph{This becomes a major issue for ND-loss when dealing with irregular geometries}, because it is almost impossible to design sample locations such that their stencils perfectly connect the domain boundary to other training samples.

One potential solution is to increase the sampling density near the boundary, i.e., $\Delta\vec{x} \rightarrow0$, to avoid such over-fitting. This however increases the complexity of the optimization problem and can cause adverse effects on training. It is also very challenging to locally refine the sampling density for PINNs with CNN-architecture. To resolve this issue without increasing the sample density, we propose \textit{BCXN-loss} to impose \textit{local structure approximation} to govern the near boundary gradient behaviour and restore connectivity between domain boundary and near-boundary samples.

\begin{figure}[!htbp]
\begin{center}
\mbox{\subfigure[\scriptsize Stencil points connect to boundary]{\includegraphics[height=2cm,keepaspectratio]{f03_stencil_connect.png}}
\hspace*{5mm}
      \subfigure[\scriptsize Stencil points fall outside domain]{\includegraphics[height=2cm,keepaspectratio]{f03_stencil_disconnect.png}}
\hspace*{5mm}
      \subfigure[\scriptsize Extrapolation along normal direction]{\includegraphics[height=2cm,keepaspectratio]{f04_BIPINN_2_ES_BIG.png}}}
\caption{(a)-(b) Schematic of a near-boundary training sample (red circle) and its stencil points (square) typically used for evaluating PDE constraint in \textit{ND-loss}. BC information is passed to stencils in fluid domain (blue square) (a) successfully in regular domain; and (b) unsuccessfully in an irregular domain due to under-defined stencils (unfilled square) outside the domain. (c) Schematic of related definitions under the present LSA-PINN framework. \textit{BCXN-loss} enforces a constraint across 3 points: external stencil point A, boundary point P, and mirror point Q inside fluid domain.}
\label{fig:stencil-schematic}
\end{center}
\vspace*{-5mm}
\end{figure}

\section{Method}
\label{methods}
\subsection{Enforcing linear constraint at near-boundary samples}
This section outlines the proposed LSA-PINN method. Here we denote the external (out-of-domain) stencil of a near-boundary sample as ES point(s). As per Fig.~\ref{fig:stencil-schematic}c, the boundary condition $\vec{u}_{BC}$ and field value $\vec{u}_{MI}$ are defined at the boundary point P and a chosen mirror point Q inside the domain, while the field value $\vec{u}_{ES}$ at ES point A is the value to be determined.

By employing Taylor series expansion on points A and Q relative to point P along $\overline{AQ}$ with respect to local coordinate $n$ (normal direction of surface), the following can be derived:
\begin{subequations}
\label{eq:bc_relation}
\begin{align}
&\vec{u}_{ES} = \vec{u}_{BC} + \overline{AP}\frac{\partial\vec{u}_{BC}}{\partial n} + \frac{\overline{AP}^2}{2}\frac{\partial^{2}\vec{u}_{BC}}{\partial n^{2}} + O(\overline{AP}^3) \\
&\vec{u}_{MI} = \vec{u}_{BC} - \overline{PQ}\frac{\partial\vec{u}_{BC}}{\partial n} + \frac{\overline{PQ}^2}{2}\frac{\partial^{2}\vec{u}_{BC}}{\partial n^{2}} + O(\overline{PQ}^3)
\end{align}
\end{subequations}
where $\overline{AP}$ and $\overline{PQ}$ are the distances between A and P and between P and Q, and $\overline{AQ} = \overline{AP} + \overline{PQ}$. The field value $\vec{u}_{ES}$ at point A can then be derived as:
\begin{equation}
\label{eq:bc_eqn}
\vec{u}_{ES} = \vec{u}_{MI} + (\vec{u}_{BC} - \vec{u}_{MI})\times\frac{\overline{AQ}}{\overline{PQ}} + O(\overline{AQ}^2)
\end{equation}
From the above equation, it can be seen that $\vec{u}_{ES}$ is linearly constrained by $\vec{u}_{BC}$ and $\vec{u}_{MI}$. While the mirror point Q can be conveniently chosen at the center of stencil, i.e., the near-boundary sample location, we choose mirror point Q along the normal direction of boundary P (ref. to Fig.~\ref{fig:stencil-schematic}c) as it is more general and performs better in our experiments. Hence, $\overline{AP} = \overline{PQ}$, and Eq. (\ref{eq:bc_eqn}) becomes:
\begin{equation}
\label{eq:bc_normal}
\vec{u}_{ES} = 2\vec{u}_{BC} - \vec{u}_{MI}
\end{equation}
A similar process can be applied for the derivation of corresponding linear constraints for Neumann- and Robin-type boundary conditions.

\subsection{Boundary connectivity (\textit{BCXN})-loss with direct forcing} \label{df-bcxn}

We can \emph{directly apply Eq. (\ref{eq:bc_normal}) to compute the field value $\vec{u}_{ES}$ for any external stencil point} during the evaluation of PDE constraint on near-boundary samples with ND-type schemes. This direct forcing approach has an association to the direct forcing immersed boundary methods in numerical computing. We can then use the following \textit{BCXN-loss}:
\begin{equation}
\label{eq:bcxn_loss}
L_{LSA-PINN} = L_{PDE(BCXN)}
\end{equation}
as our LSA-PINN's PDE loss term whereby the PDE constraints on near-boundary samples are modulated by the local linear constraint. In such implementation, the BCs are implicitly infused into the training loss, obviating the BC loss term.

\subsection{Procedure to evaluate the field value $\vec{u}_{MI}$ at mirror point} \label{mirror-point}
Note that PINN methods in general require a routine to sample from inside the geometry of interest as well as the boundary. As compared to conventional PINNs, LSA-PINN involves additional steps. We first describe the steps to compute the mirror point location (along normal direction):
\begin{itemize}
\item \textbf{Determine whether a stencil is external.} In this study, we first construct a level set function $\phi$ of the geometry of interest ~\cite{osher2003geometric}, where the shortest distance $\overline{AP}$ between the stencil point A and boundary point P can be found. From the sign of level set function, we can determine whether a stencil is external.
\item \textbf{Compute mirror point location for external stencil.} Once $\overline{AP}$ is determined, we compute the location of mirror point Q inside the fluid domain, with the distance $\overline{AP} = \overline{PQ}$.
\end{itemize}
In practice, given a fixed set of training samples, all the external stencil points and their mirror locations can be pre-computed. On the other hand, the field value $\vec{u}_{MI}$ at mirror points depends on the LSA-PINN’s output and are evaluated inside a training iteration:
\begin{itemize}
\item \textbf{MLP-architecture.} The $\vec{u}_{MI}$ value can be directly obtained by evaluating $\vec{u}_{MI}(\vec{x};\textbf{\textit{w}})$ at the mirror point.
\item \textbf{CNN-architecture.} When the location of mirror point does not coincide with the CNN grid, the following inverse-distance-weighted interpolation function can be utilized to obtain the $\vec{u}_{MI}$ value at mirror point ~\cite{chiu2021development}.
\end{itemize}
\begin{subequations}
\label{eq:interpolations}
\begin{align}
&u(\Tilde{x},\Tilde{y}) = \Sigma u(x,y) \phi '\Bigl(\frac{x-\Tilde{x}}{h}\Bigr) \phi '\Bigl(\frac{y-\Tilde{y}}{h}\Bigr) \\
&\phi '\Bigl(\frac{x-\Tilde{x}}{h}\Bigr) = \phi\Bigl(\frac{x-\Tilde{x}}{h}\Bigr) / \Sigma\phi \Bigl(\frac{x-\Tilde{x}}{h}\Bigr) \\
&\phi '\Bigl(\frac{y-\Tilde{y}}{h}\Bigr) = \phi\Bigl(\frac{y-\Tilde{y}}{h}\Bigr) / \Sigma\phi \Bigl(\frac{y-\Tilde{y}}{h}\Bigr) \\
&\phi(r) =
\begin{cases}
    \frac{1}{8}\Bigl(3 - 2|r| + \sqrt{1 + 4|r| + 4r^{2}}\Bigr), & |r| \leq 1 \\
    \frac{1}{8}\Bigl(5 - 2|r| - \sqrt{-7 + 12|r| - 4r^{2}}\Bigr), & 1 < |r| \leq 2 \\
    0,              & |r| > 2
\end{cases}
\end{align}
\end{subequations}
$(\Tilde{x},\Tilde{y})$ is the location of mirror point to be interpolated, while $(x,y)$ is the location of training samples inside the fluid domain. In this study, $h$ is chosen as $\Delta x$, which leads to a $2\Delta x$ interpolation radius. Once $\vec{u}_{MI}$ and $\vec{u}_{BC}$ are evaluated, we can compute $\vec{u}_{ES}$ based on Eq. (\ref{eq:bc_normal}) during training.

\section{Results}

\begin{figure*}[htbp]
\begin{center}
\centerline{\includegraphics[width=0.95\linewidth]{Figure2-Wavy-Flow-Revised.pdf}}
\caption{Comparison of $u$-, $v$-velocity and pressure $p$ contour between the ground truth and PINN solutions (median MSE from 10 runs), for the wavy channel flow problem, $Re=100$.}
\label{fig:wavy-channel-show}
\end{center}
\vspace*{-5mm}
\end{figure*}

In this work, we study the performance of the proposed LSA-PINN for learning multiple complex fluid dynamical systems, using both MLP and CNN architectures. We further demonstrate that LSA-PINN is broadly applicable to different ND-type training losses, including ND-type (\textit{BCXN-ND}) and CAN-type (\textit{BCXN-CAN}) schemes.

The LSA-PINN (\textit{BCXN-ND-}, \textit{BCXN-CAN-loss}) is compared with baseline PINN (\textit{ND-}, \textit{CAN-}, \textit{AD-loss}) models, based on mean squared error (MSE) of the respective PINN solutions' velocity vector $\vec{u}$ relative to ground truth. To facilitate comparison, identical network architectures, initialization distribution, and training settings are employed. Relevant settings are summarized in Table~\ref{tab01} and Table~\ref{tab02} in Appendix \ref{table-settings}. For each experiment, we perform and present results based on 10 independent runs with different initialization. The ground truth solutions are obtained by an in-house numerical solver based on the improved divergence-free condition compensated (IDFC) method~\cite{chiu2018improved}. The 3 fluid dynamics test cases are:
\begin{itemize}
    \item \textbf{2D semi-circle lid-driven cavity flow, $Re=1000$.} The fluid flow inside the semi-circle cavity is driven by a lid velocity $u_{lid}=1,v_{lid}=0$ at the top wall. No-slip condition $(u=v=0)$ is applied to the other wall. There is a primary eddy near the cavity's center.
    \item \textbf{2D lid-driven cavity flow, $Re=1000$.} The fluid flow inside a $1 \times 1$ unit square cavity is driven by a lid velocity $u_{lid}=1,v_{lid}=0$ at the top wall. No-slip condition is applied to the other walls. The lid-driven cavity flow has been widely chosen as a benchmark case for many numerical methods due to the complex physics encapsulated within. At $Re=1000$, there is a primary eddy near the center of the cavity, with secondary and tertiary eddies at the bottom-right and bottom-left regions.
    \item \textbf{2D wavy channel flow, $Re=100$.} The fluid flow passing through a long wavy channel is studied. The inlet profile at left boundary is defined as $u(0,y)=-\frac{3}{2} y^{2}+\frac{3}{2},v=0$. A non-slip condition is applied to the top and bottom walls, while outlet boundary conditions $(\frac{\partial u}{\partial x}=\frac{\partial v}{\partial x}=p=0)$ are applied to the right boundary.

\end{itemize}

\subsection{MLP-architecture LSA-PINN \label{mlp-pinn}}

In the MLP-architecture experiment, we evaluate the LSA-PINN's performance and compare them to other \textit{AD-}, \textit{CAN-}, and \textit{ND-loss} baseline methods previously published for the MLP architecture. All models were trained on a relatively sparse set of training sample (4,230, 10,400, and 15,200 points for semi-circle lid-driven cavity, lid-driven cavity, and wavy channel flow problems, respectively) for 300,000 iterations using Adam algorithm. The average MSE and relative L2 error across multiple independent runs are summarized in Table \ref{tab: mlp-error-all}.

\begin{table}[htbp]
\caption{Model Performance Comparisons for MLP-architecture PINNs}
\begin{center}
\begin{adjustbox}{width=.48\textwidth}
\begin{tabular}{cccc}
\hline
\textbf{Problem} & \textbf{Model} & \textbf{MSE} & \textbf{Rel. Error} \\
\hline
\multirow{6}{*}{\parbox{3cm}{\centering Semi-circle lid-driven cavity flow, $Re = 1000$}} & {\textit{BCXN-CAN}} & {$6.21 \times 10^{-5}$} & {$3.55 \times 10^{-2}$} \\
\cline{2-4}
{} & {\textit{BCXN-ND}} & {$1.61 \times 10^{-4}$} & {$5.95 \times 10^{-2}$} \\
\cline{2-4}
{} & {\textit{CAN}} & {$2.82 \times 10^{-2}$} & {$7.62 \times 10^{-1}$} \\
\cline{2-4}
{} & {\textit{ND}} & {$3.91 \times 10^{-2}$} & {$8.56 \times 10^{-1}$} \\
\cline{2-4}
{} & {\textit{AD (SIREN)}} & {$1.74 \times 10^{-2}$} & {$5.95 \times 10^{-1}$} \\
\cline{2-4}
{} & {\textit{AD} (Vanilla)} & {$5.46 \times 10^{-2}$} & {$9.83 \times 10^{-1}$} \\
\hline
\multirow{3}{*}{\parbox{3cm}{\centering Lid-driven cavity flow, $Re = 1000$}} & {\textit{BCXN-CAN}} & {$5.10 \times 10^{-5}$} & {$3.38 \times 10^{-2}$} \\
\cline{2-4}
{} & {\textit{CAN}} & {$4.23 \times 10^{-2}$} & {$9.77 \times 10^{-1}$} \\
\cline{2-4}
{} & {\textit{AD (SIREN)}} & {$3.47 \times 10^{-2}$} & {$8.87 \times 10^{-1}$} \\
\hline
\multirow{3}{*}{\parbox{3cm}{\centering Wavy channel flow, $Re = 100$}} & {\textit{BCXN-CAN}} & {$1.34 \times 10^{-4}$} & {$7.28 \times 10^{-2}$} \\
\cline{2-4}
{} & {\textit{CAN}} & {$5.96 \times 10^{-1}$} & {$2.55 \times 10^{0}$} \\
\cline{2-4}
{} & {\textit{AD (SIREN)}} & {$3.52 \times 10^{-2}$} & {$7.19 \times 10^{-1}$} \\
\hline
\end{tabular}
\end{adjustbox}
\begin{enumerate}[(I)]

\item
Implementation for \textit{CAN-}PINN models based on previous work~\cite{chiu2022can}.
\item
Implementation for \textit{SIREN-}PINN models based on previous work~\cite{sitzmann2020implicit,wang2021eigenvector,wong2022learning}.
\item
Implementation for vanilla PINN models based on previous work~\cite{raissi2019physics}.
\end{enumerate}

\label{tab: mlp-error-all}
\end{center}
\end{table}

The results indicate that our present LSA-PINN method can successfully learn a good solution with low MSE across a diverse set of problems with irregular domains and complex physics. All 3 test cases show noticeable improvement in solution accuracy ($>$2 orders of magnitude lower $\vec{u}$ MSE) with the LSA-PINN method. The baseline PINNs fail to produce a reasonable solution, i.e., they cannot properly learn the correct flow with current sampling density and training iteration, see example in Fig.~\ref{fig:wavy-channel-show}.

\textbf{Case study of information propagation.} Information propagation from the boundary to the inner domain has been suggested to be critical for successful PINN modelling ~\cite{wang2022respecting, daw2022rethinking}. Hence, we investigate the evolution of the PDE residual across multiple training runs for the semi-circle lid-driven cavity with and without \textit{BCXN-loss} to better understand the impact of enhanced connectivity and \textit{local structure approximation} imposed by the \textit{BCXN-loss}. Illustrative contour plots of the PDE residuals and velocity at various training iterations are provided in Fig. \ref{fig:BCXN-error-prop}.

\begin{figure}[H]
\begin{center}
\centerline{\includegraphics[width=0.5\textwidth,keepaspectratio]{Figure3-PDE-Error-Propagation-Revised}}
\vspace*{-3mm}
\caption{Contour plots of PDE error residual and velocity for (a) LSA-PINN (\textit{BCXN-ND}) and (b) PINN (\textit{ND}), for the semi-circle lid-driven cavity flow problem, $Re=1000$.}
\label{fig:BCXN-error-prop}
\end{center}
\vspace*{-6mm}
\end{figure}

As described in the previous sections, the \textit{ND} methods typically fail to converge to a good model even after 300,000 training iterations but this can be improved via the addition of \textit{BCXN-loss} during training. The velocity contours clearly illustrate the improved convergence. In addition, there is a striking difference in the distribution of PDE residuals across the domain during the early stages of training. In a typical \textit{ND}-PINN, the PDE residuals rapidly reduce to zero, especially in the parts of the domain far away from the boundary condition. This also corresponds to the \textit{ND-loss} PINN velocity profile converging to a wrong solution, and is indicative of the information from the boundary condition being lost within the domain. On the other hand, the PDE residuals within the domain persist for a much longer period when the \textit{BCXN-loss} is included, which is indicative of information from the top boundary being persistently transmitted to the inner domain until the LSA-PINN model converges. This is also consistent with the physical dynamics of the system, as we note that the lid-driven cavity flow, as the name suggests, is heavily dependent on the top boundary's motion to drive the dynamical flow within the entire domain.

\textbf{Case study of samples with different boundary connectivity.} The lid-driven cavity test case has a regular domain, but we chose Cell-centered (whereby near-boundary points need support points outside the domain) sample locations to show that the proposed method can be useful for this kind of scenarios. To further illustrate the impact of having near-boundary stencil points outside of the domain, we also investigate the impact of Cell-centered sampling points against Node-centered (whereby near-boundary points use support points that are coincident with the domain boundary). These two sets of sampling points are tested on a lid-driven cavity flow at $Re = 1000$ with different batch sizes and point sampling density for the \textit{BCXN-CAN}, \textit{CAN} and \textit{AD} models, and are trained for $300,000$ iterations each. The MSE and L2 relative errors are collated across multiple independent runs and tabulated in Table \ref{tab: ldc-sampling-effect}. A qualitatively good solution for this problem corresponds to an MSE of $1 \times 10^{-3}$ and below.

\begin{table}[htbp]
\caption{Model Performance for Different Sampling Densities. [C] and [N] refer to cell- and node-centered sample arrangements.}
\begin{center}
\begin{adjustbox}{width=.49\textwidth}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multirow{2}{*}{\textbf{Batch Size}} & \multirow{2}{*}{\textbf{Points}} &
\multirow{2}{*}{\textbf{Model}} & \multicolumn{2}{|c|}{\textbf{$50 \times 50$}} &
\multicolumn{2}{|c|}{\textbf{$100 \times 100$}} \\
\cline{4-7}
{} & {} & {} & {MSE} & {Rel. Error} & {MSE} & {Rel. Error} \\
\hline
\multirow{4}{*}{$500$} & \multirow{2}{*}{[C]} & {\textit{BCXN-CAN}} & {\scriptsize$2.63 \times 10^{-4}$} & {\scriptsize$7.75 \times 10^{-2}$} & {\scriptsize$3.06 \times 10^{-3}$} & {\scriptsize$2.63 \times 10^{-1}$} \\
\cline{3-7}
{} & {} & {\textit{CAN \& AD}} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} \\
\cline{2-7}
{} & \multirow{2}{*}{[N]} & {\textit{CAN}} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} \\
\cline{3-7}
{} & {} & {\textit{AD}} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} \\
\hline
\multirow{4}{*}{$1000$} & \multirow{2}{*}{[C]} & {\textit{BCXN-CAN}} & {\scriptsize$2.23 \times 10^{-4}$} & {\scriptsize$7.13 \times 10^{-2}$} & {\scriptsize$3.97 \times 10^{-4}$} & {\scriptsize$9.47 \times 10^{-2}$} \\
\cline{3-7}
{} & {} & {\textit{CAN \& AD}} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} \\
\cline{2-7}
{} & \multirow{2}{*}{[N]} & {\textit{CAN}} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} \\
\cline{3-7}
{} & {} & {\textit{AD}} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} \\
\hline
\multirow{4}{*}{$2000$} & \multirow{2}{*}{[C]} & {\textit{BCXN-CAN}} & {\scriptsize$2.36 \times 10^{-4}$} & {\scriptsize$7.34 \times 10^{-2}$} & {\footnotesize$7.14 \times 10^{-5}$} & {\footnotesize$4.02 \times 10^{-2}$} \\
\cline{3-7}
{} & {} & {\textit{CAN \& AD}} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} \\
\cline{2-7}
{} & \multirow{2}{*}{[N]} & {\textit{CAN}} & {\scriptsize$7.41 \times 10^{-4}$} & {\scriptsize$1.32 \times 10^{-1}$} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} \\
\cline{3-7}
{} & {} & {\textit{AD}} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} \\
\hline
\multirow{4}{*}{$4000$} & \multirow{2}{*}{[C]} & {\textit{BCXN-CAN}} & {\cellcolor{lightgray}} & {\cellcolor{lightgray}} & {\footnotesize$2.38 \times 10^{-5}$} & {\footnotesize$2.31 \times 10^{-2}$} \\
\cline{3-7}
{} & {} & {\textit{CAN \& AD}} & {\cellcolor{lightgray}} & {\cellcolor{lightgray}} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} \\
\cline{2-7}
{} & \multirow{2}{*}{[N]} & {\textit{CAN}} & {\cellcolor{lightgray}} & {\cellcolor{lightgray}} & {\scriptsize$2.94 \times 10^{-4}$} & {\scriptsize$8.24 \times 10^{-2}$} \\
\cline{3-7}
{} & {} & {\textit{AD}} & {\cellcolor{lightgray}} & {\cellcolor{lightgray}} & {\tiny$>1 \times 10^{-2}$} & {\tiny$>5 \times 10^{-1}$} \\
\hline
\end{tabular}
\label{tab: ldc-sampling-effect}
\end{adjustbox}
\end{center}
\end{table}

From this experiment, we notice that the Node-centered \textit{CAN-loss} PINN model has improved performance relative to the Cell-centered version, for a larger batch size and training samples. This further emphasizes our hypothesis that accurate modelling of the boundary conditions is an important factor in obtaining good performance. As this is a square domain, the use of a sampling point distribution that has perfect connectivity with the boundary of the domain naturally leads to improved performance.

In addition, the inclusion of our proposed \textit{BCXN-loss} on a Cell-centered LSA-PINN can improve model performance to achieve a lower error than the properly aligned baseline PINN models, especially for sparse training samples setting. This suggests that the use of LSA-PINN can afford flexibility in the choice of sampling strategies, whereas conventional sampling strategies may require sampling that maintains connectivity with the boundaries.

\textbf{Case study with different sampling density.} From the above experiment, we notice that a 4x increase in sampling density (from $50 \times 50$ to $100 \times 100$) can improve model performance. Hence, we further study the lid-driven cavity flow problem to understand the trade-off between convergence in accuracy and sampling density. We train on a denser (i.e., $200 \times 200$) set of samples until the accuracy of the \textit{CAN} and \textit{AD} models start to improve. Note that denser sampling requires larger training iteration and batch size, and the tuning of these training hyper-parameters can be very time consuming. Although all PINN models can eventually achieve improved accuracy with denser sampling, Fig.~\ref{fig:convergence-ldc} indicates that it is much harder and slower to learn from a dense set of samples, thereby demonstrating the benefits of LSA-PINN. Encouragingly, LSA-PINN learns more accurate solutions ($>$1 order of magnitude lower $\vec{u}$ MSE) than the baseline \textit{CAN-} and \textit{AD-loss} PINNs with less ($<$50\%) training iterations, while using fewer ($<$50\%) training samples for the lid-driven cavity flow problem.

\begin{figure}
\centering
\subfigure{\includegraphics[width=0.5\textwidth,keepaspectratio]{Figure4-Convergence-Revised.png}}
\vspace*{-7mm}
\caption{PINNs learning the solution of lid-driven cavity flow under (a) sparse and (b) dense training scenarios respectively. The convergence trends of $\vec{u}$ MSE are plotted. Bold lines indicate the median convergence path, and the shaded areas indicate the inter-percentile range from $10^{th}-90^{th}$ percentiles.}
\label{fig:convergence-ldc}
\vspace*{-5mm}
\end{figure}

\subsection{CNN-architecture LSA-PINNs} \label{cnn-fast}
The \textit{BCXN-loss} method does not pose any requirement on the differentiable property of the networks, hence it can be freely implemented in any neural network architecture, including CNN. Hence, we compare the LSA-PINN \textit{(BCXN-ND)} to baseline PINN \textit{(ND)} method previously published for the CNN architecture and summarize their average MSE and relative L2 error across multiple independent runs in Table \ref{tab: cnn-error-all}. The U-Net architecture and training settings used are summarized in Table~\ref{tab02}. All PINN models for this section are trained on a single discretized 2D input which consists of 3 channels, namely the x-, y-coordinate, and an indicator to differentiate fluid and non-fluid domain. They were trained for 200,000 iterations using Adam algorithm on a relatively sparse set of training samples, i.e., the semi-circle and square lid-driven cavity flow problems use 3,930 and 10,000 interior points which belong to the fluid domain respectively.

\begin{table}[htbp]
\caption{Model Performance Comparisons for CNN-architecture PINNs}
\begin{adjustbox}{width=.48\textwidth}
\begin{tabular}{|c|c|c|c|}
\hline
\textbf{Problem} & \textbf{Model} & \textbf{MSE} & \textbf{Rel. Error} \\
\hline
\multirow{2}{*}{\parbox{3cm}{\centering Semi-circle lid-driven cavity flow, $Re = 1000$}} & {\textit{BCXN-ND}} & {$1.14 \times 10^{-4}$} & {$4.92 \times 10^{-2}$} \\
\cline{2-4}
{} & {\textit{ND}} & {$3.96 \times 10^{-2}$} & {$8.55 \times 10^{-1}$} \\
\hline
\multirow{2}{*}{\parbox{3cm}{\centering Lid-driven cavity flow, $Re = 1000$}} & {\textit{BCXN-ND}} & {$2.52 \times 10^{-5}$} & {$4.46 \times 10^{-2}$} \\
\cline{2-4}
{} & {\textit{ND}} & {$2.38 \times 10^{-2}$} & {$1.00 \times 10^{0}$} \\
\hline
\multirow{2}{*}{\parbox{4cm}{\centering Transient semi-circle lid-driven cavity flow, $Re = 500$}} & {\textit{BCXN-ND}} & {$3.95 \times 10^{-4}$} & {$1.29 \times 10^{-1}$} \\
\cline{2-4}
{} & {\textit{ND}} & {$9.66 \times 10^{-3}$} & {$6.81 \times 10^{-1}$} \\
\hline
\end{tabular}
\label{tab: cnn-error-all}
\end{adjustbox}
\end{table}

Like the MLP experiment, LSA-PINN for CNNs show a noticeable improvement in the solution accuracy (2-3 orders of magnitude lower $\vec{u}$ MSE) over the baseline \textit{ND-loss} PINN. The results indicate that LSA-PINN can effectively learn the solution of fluid dynamic problems with U-Net architecture and that the \textit{BCXN-loss} can benefit CNN-architecture PINNs by remedying training issues stemming from a discretized input being unable to exactly resolve complex geometries.

\textbf{Transient 2D semi-circle lid-driven cavity flow, $Re=500$.} In addition to the steady-state semi-circle lid-driven cavity flow, we further evaluate the improvement obtained from application of CNN-architecture LSA-PINN to the modelling of a transient version of this problem. In this instance, the problem's initial condition corresponds to that of a fluid domain at rest (at t = 0). The PINN model is then used to model the evolution of the flow within a semi-circle lid-driven cavity from $t = 0$ to $t = 2$, $\Delta{t} = 0.1$. Model settings are provided in Table~\ref{tab02}. The transient PINN models have additional $t$ channel in their input.

The performance of LSA-PINN \textit{(BCXN-ND)} and baseline PINN \textit{(ND)} are compared in Table \ref{tab: cnn-error-all}. The average errors for all time points between $t = 0$ and $t = 2.0$ yield MSEs of 3.95e-4 and 9.66e-3 with and without the inclusion of \textit{BCXN-loss} respectively. Hence, the inclusion of \textit{BCXN-loss} in this transient problem improves MSE by $>1$ order of magnitude. Velocity contour plots of the \textit{BCXN-ND} and \textit{ND} models are provided in Fig. \ref{fig:transient-contour} to illustrate the differences in convergence.

\begin{figure}[htbp]
\begin{center}
\centerline{\includegraphics[width=0.5\textwidth,keepaspectratio]{Figure5-Transient-Revised.pdf}}
\caption{Velocity contours of the transient simulation at $t = 0.5$, $t = 1.0$ and $t=2.0$ for the ground truth and PINN solutions.}
\label{fig:transient-contour}
\end{center}
\vspace*{-6mm}
\end{figure}

\textbf{Performance with increasing problem complexity.} To better understand the improvement in accuracy obtainable by the \textit{BCXN-loss}, we also investigate the utility of LSA-PINN for modelling more complex problems more efficiently. Hence, we train both LSA-PINN and baseline PINN for the semi-circle lid-driven cavity flow across a range of $Re$. As $Re$ increases, the complexity of the developed flow typically requires more sophisticated algorithms and training samples for better capture of the non-linear effects.

The average MSE and relative L2 error obtained across multiple runs for various $Re$ are presented in Table \ref{tab: semi-ldc-table}. We clearly see that the performance difference is increasingly pronounced for more complicated (higher $Re$) problems, thus demonstrating that strategies like LSA-PINN methods are increasingly important for more complex problems.

\begin{table}[htbp]
\caption{Model Performance for the 2D Semi-circular Lid-Driven Cavity Problem}
\begin{adjustbox}{width=.48\textwidth}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multirow{2}{*}{\textbf{Metric}} & \multirow{2}{*}{\textbf{Model}} & \multicolumn{5}{|c|}{\textbf{Reynolds Number (\textit{Re})}} \\
\cline{3-7}
{} & {} & {$50$} & {$100$} & {$200$} & {$400$} & {$1000$} \\
\hline
\multirow{2}{*}{MSE} & {\textit{BCXN-ND}} & {$2.51 \times 10^{-5}$} & {$2.76 \times 10^{-5}$} & {$3.56 \times 10^{-5}$} & {$4.58 \times 10^{-5}$} & {$1.14 \times 10^{-4}$}\\
\cline{2-7}
{} & {\textit{ND}} & {$2.98 \times 10^{-4}$} & {$2.79 \times 10^{-4}$} & {$3.74 \times 10^{-4}$} & {$2.33 \times 10^{-3}$} & {$3.96 \times 10^{-2}$} \\
\hline
\multirow{2}{*}{Rel. Error} & {\textit{BCXN-ND}} & {$0.0258$} & {$0.0269$} & {$0.0305$} & {$0.0335$} & {$0.0492$} \\
\cline{2-7}
{} & {\textit{ND}} & {$0.0823$} & {$0.0732$} & {$0.0780$} & {$0.211$} & {$0.855$} \\
\hline
\end{tabular}
\label{tab: semi-ldc-table}
\end{adjustbox}
\end{table}

\textbf{Meta-PINN model.} The CNN-architecture PINN models represent the entire spatial domain as a single input. We can train the model on multiple samples, each representing a different scenario such as changes in geometry or physical behaviour. Such a meta-PINN model can be used to predict new scenarios. Hence, we further study the ability of LSA-PINN to learn the solutions to multiple flow scenarios concurrently (spanning $Re=50-1000$) on semi-circle lid-driven cavity geometry. The model architecture and training setting is given in Appendix \ref{table-settings} Table~\ref{tab02}. The LSA-PINN \textit{(BCXN-ND)} model is used to predict $Re=300, 500, 600, 700, 900$ scenarios which are not used during training. The results are summarized in Fig.~\ref{fig:meta-ldc}. The present LSA-PINN method successfully learns a good solution with low $\vec{u}$ MSE (i.e., $<$1e-4) for most of the $Re$ scenarios, and achieves good prediction accuracy.

\begin{figure}[H]
\centering
\subfigure{\includegraphics[width=\linewidth,keepaspectratio]{Figure6-Meta-Model-Revised.png}}
\vspace*{-6mm}
\caption{LSA-PINN learning the solutions of multiple 2D semi-circle lid-driven cavity flow scenarios concurrently, using CNN architecture. The $\vec{u}$ MSE between ground truth and LSA-PINN solution for train and test $Re$ are plotted. Inset are $u$- and $v$-velocity contours of the ground truth and LSA-PINN solutions for $Re=600$ and $Re=900$.}
\label{fig:meta-ldc}
\end{figure}

\subsection{Impact of \textit{BCXN}-loss on training time} \label{training-time}
We further compare the time taken to train a model with and without the \textit{BCXN-loss} on the MLP- and CNN-type models as applied to the semi-circle lid-driven cavity problem at $Re = 1000$. The average MSE and relative L2 error obtained are averaged across 10 independent runs, and the minimum and maximum wall-clock time are recorded. All experiments are run with standard Tensorflow and Keras implementation on a 20-core workstation with an Intel Xeon Gold 6248 processor and an Nvidia RTX 2080 Ti GPU.

Quantitative comparisons for the MLP-architecture and CNN-architecture models are presented in Table \ref{tab: mlp-error-time}. For fair comparison within architecture types, all MLP models were run for an identical number of training iterations (300,000) with the same batch size (1,000) and model size while all CNN models were similarly run for an identical number of training iterations (200,000) and with the same model size.

\begin{table}[htbp]
\caption{Model Performance for different PINN architectures}
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline
\textbf{Model} & \textbf{Method} & \textbf{MSE} & \textbf{Rel. Error} & \textbf{Time (min)} \\
\hline
{MLP} & {\textit{BCXN-CAN}} & {$6.21 \times 10^{-5}$} & {$3.55 \times 10^{-2}$} & {$55-57$} \\
\hline
{MLP} & {\textit{BCXN-ND}} & {$1.61 \times 10^{-4}$} & {$5.95 \times 10^{-2}$} & {$23-24$} \\
\hline
{MLP} & {\textit{CAN}} & {$2.82 \times 10^{-2}$} & {$7.62 \times 10^{-1}$} & {$40-46$} \\
\hline
{MLP} & {\textit{ND}} & {$3.91 \times 10^{-2}$} & {$8.56 \times 10^{-1}$} & {$20-22$} \\
\hline
{MLP} & {\textit{AD (SIREN)}} & {$1.74 \times 10^{-2}$} & {$5.95 \times 10^{-1}$} & {$58-61$} \\
\hline
{MLP} & {\textit{AD} (Vanilla)} & {$5.46 \times 10^{-2}$} & {$9.83 \times 10^{-1}$} & {$50-55$} \\
\hline
{CNN} & {\textit{BCXN}} & {$1.14 \times 10^{-4}$} & {$4.92 \times 10^{-2}$} & {$126-136$} \\
\hline
{CNN} & {\textit{ND}} & {$3.96 \times 10^{-2}$} & {$8.55 \times 10^{-1}$} & {$121-137$} \\
\hline
\end{tabular}
\label{tab: mlp-error-time}
\end{center}
\end{table}

From this study, we note that the inclusion of \textit{BCXN-loss} across all the cases only slightly increases the amount of time required for model training. This increase was on the order of less than 5 \% for the CNN-type architecture and the ND-MLP-type model. Critically, while the time taken may have increased slightly, the MSE typically decreases more than 1 order of magnitude. In many instances, this trade-off between model training time and performance with the inclusion of a \textit{BCXN-loss} is probably preferable.

\section{Conclusion}
We present a LSA-PINN method which can more efficiently learn the solution to PDEs for complex geometry with sparse training samples. This is accomplished by enforcing a linear constraint during the training implicitly. While linear constraints are applied in this work due to truncation of the Taylor series expansion in the derivation, it is worth studying if different, higher order approximations can achieve better performance in future work. Overall, our comprehensive experimental studies demonstrate practical advantages of LSA-PINN on diverse test problems in fluid dynamics, improving the accuracy of solutions by orders of magnitude even while requiring much less training iterations and samples.

\bibliographystyle{IEEEtran}