#  A Knowledge-driven Physics-Informed Neural Network model; Pyrolysis and
Ablation of Polymers

Aref Ghaderi  Ramin Akbari  Yang Chen  Roozbeh Dargazany111Corresponding
author. Department of Civil and Environmental Engineering, Michigan State
University<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In aerospace applications, multiple safety regulations were introduced to
address damage associated with pyrolysis. Predictive modeling of pyrolysis is
a challenging task since multiple thermo-chemo-mechanical laws need to be
concurrently solved at each time step. So far, classical modeling approaches
were mostly focused on defining the basic chemical processes (pyrolysis and
ignite) at micro-scale by decoupling them from thermal solution at the micro-
scale and then validating them using meso-scale experimental results. The
advent of Machine Learning (ML) and AI in recent years has provided an
opportunity to construct quick surrogate ML models to replace high fidelity
multi-physics models, which have a high computational cost and may not be
applicable for high nonlinear equations.

We developed a hybrid knowledge-infused Physics-informed neural network (PINN)
engine to simulate multiple stiff, and semi-stiff ODEs that govern Pyrolysis
and Ablation in elastomers. Our Engine is particularly developed to calculate
the char formation and degree of burning in the course of pyrolysis of
crosslinked polymeric systems. A multi-task learning approach is hired to
reduce training costs, while the search space is constrained using our
imprecise knowledge/expectations of the solution field. The proposed Hybrid-
PINN (H-PINN) solver was bench-marked against finite element high fidelity
solutions on different examples. We developed PINN architectures using
collocation training to forecast temperature distributions and the degree of
burning in the course of pyrolysis in multiple one- and two-dimensional
examples. By decoupling thermal and mechanical equations, we can predict the
loss of performance in the system by predicting the char formation pattern and
localized degree of burning at each continuum.

###### keywords:

Physics-Informed Neural Networks , Coupled PDEs , Deep Learning , Pyrolysis

††journal: Engineering Applications of Artificial Intelligence

##  1 Introduction

Due to their outstanding thermal adaptability and high thermal protection
efficiency, ablative thermal protection systems have been increasingly being
employed for aerospace vehicle and scramjet engine applications. Thermal
protection materials are faced with a significant problem due to the next
generation aerospace vehicles and engines’ growing structural efficiency
requirements. It is strategically important to use thermal protection coatings
to assure survivability of aircraft components from extreme heat. There are
two types of thermal barrier; (1) ceramic and (2) Polymer-matrix thermal
protection coatings.

Polymeric coatings used in supersonic aircraft/spacecrafts can be exposed to
+800oC temperatures or +1MW/m2 heat fluxes for few minutes. Among polymeric
coatings, silicone rubber has the heat resistance of inorganic polymer and
flexibility of organic polymer, due to the presence of Si-O-Si bond and side
organic group in backbone chains. Because of its exceptional heat tolerance,
rubber coatings as thermal barriers are becoming increasingly popular and thus
predicting the pyrolysis process of rubber has become the subject of interest
. While different aspects of pyrolysis, such as yhermal decomposition,
degradation, and evolution of the rubber matrix have been extensively studied
in the last decade, simulating pyrolysis remains to be extremely expensive due
to the multi-physics nature of the events. Here, our goal is to leverage a
physics-informed machine-learned engine to substitute the governing equations
and subsequently simulate the pyrolysis at a substantially lower cost.

Pyrolysis is a chemical decomposition process that initiates at elevated
temperatures. It highly differs from external processes like combustion and
hydrolysis which requires presence of additional reagents such as oxygen
($O_{2}$) or water (in hydrolysis). Pyrolysis results in solids (char),
condensable liquids (tar), and non-condensing/permanent gases
[natali2016ablation, chronopoulos2013thermal, liao2022hybrid].

To model pyrolysis in elastomeric materials, one need to describe the coupled
effects of thermal, chemical and physical processes [luo2007thermo,
luo2011thermo]. Table 1 summarizes the different thermal, chemical, and
physical processes that may be investigated in burning.

Table 1: A summary of processes that can be anticipated in burning analysis. | Heat conduction/temperature
---
[henderson1985model], [henderson1987mathematical], [florio1991study]
| Pyrolysis
---
[sullivan1992finite], [sullivan1990finite], [mouritz2009review]
| Volatile convective flow
---
[carpier2020tensile], [feih2012tensile], [shi2013high]
| Char formation/mass loss
---
[carpier2020tensile], [feih2012tensile], [shi2013high]
| Internal pressure
---
[sullivan1992finite], [sullivan1990finite], [mouritz2009review]
| Thermal expansion
---
[mouritz2009review], [feih2012tensile], [shi2013high]
| Thermal stresses
---
[kulkarni2018current], [li2020simulation], [vatikiotis1980heat]

Before getting into the details of the models, it is crucial to understand
what happens when the polymer is exposed to a high-temperature heat source.

Chemical Processes: Polymers can decompose thermally by oxidative reactions or
just by being exposed to heat. Several main groups of chemical processes play
significant roles in the heat decomposition of polymers, including (1) Chain-
stripping, in which atoms or groups not a part of the polymer chain (or
backbone) are cleaved; (2) cross-linking, in which bonds are formed between
polymer chains. (3) Random-chain scission, in which chain scissions occur at
seemingly random locations in the polymer chain; (4) End-chain scission, in
which individual monomer units are successively removed at the chain end; It
suffices to state here that heat decomposition of a polymer often involves
many reactions from each of these kinds. However, these broad categories offer
a conceptual framework that is helpful for comprehending and categorizing
polymer decomposing behavior (see Fig. 1).

Figure 1: General decomposition mechanisms.

Physical Processes: The nature of the material may have an impact on the many
physical processes that take place during heat decomposition. For instance,
simple phase transitions upon heating are not achievable for thermosetting
polymeric materials because they are infusible and insoluble after they have
been synthesized. On the other hand, if heating does not go over the minimum
thermal breakdown temperature, thermoplastics can be softened by heating
without the material undergoing irreversible alterations. This gives
thermoplastic materials a significant edge in terms of how simple it is to
mold or thermoform items. Carbonaceous chars are created during the thermal
decomposition of several materials, including cellulosic, thermosetting, and
thermoplastic ones. The ongoing thermal decomposition process will be
significantly impacted by the physical makeup of these chars. The pace of heat
decomposition of the remaining polymer is frequently determined by the
physical properties of the char. Even though char production is a chemical
process, its relevance is mostly because of its physical characteristics.

The polymer continues to decompose endothermically until the reaction zone
reaches the material’s back-face when the rest of the polymer is degraded to
volatiles and char.

Figure 2: Physical and chemical changes during thermal decomposition
[beyler2002thermal].

The physics of this thermo-chemical reaction is well understood, and it is
described by a collection of coupled nonlinear partial differential equations
(PDEs) that describe heat conduction and burning kinetics
[johnston1996process, boyard2016heat]. These PDEs (see section 2), however, do
not have a closed-form solution; computational approximation is necessary.

Traditional solvers vs. Data-driven methods. The finite element method (FEM)
is a common approach that (1) uses a set of basis functions within discretized
sub-domains of the geometry to approximate the solution to the PDE at each
time point and (2) iteratively simulates the evolution of the basis function
coefficients forward in time [zienkiewicz2005finite, boisse2005mesoscopic].
Comes with a restrictive computational cost in the case of repetitive
simulations, such as optimization, control, real-time monitoring,
probabilistic modeling, and uncertainty quantification. This is due to the
fact that forward simulations require solving large non-linear equations
repeatedly.

In general engineering applications, the advent of Machine Learning (ML) and
AI in recent years has provided an opportunity to construct quick surrogate ML
models to replace classical FEM. Classic neural networks, on the other hand,
map across finite-dimensional spaces and can thus only learn discretization-
specific solutions. This is frequently a constraint in actual applications,
necessitating the creation of mesh-invariant neural networks. The finite-
dimensional operators and Neural-FEM are two popular neural network-based
techniques for solving multiple PDEs [kochkov2021machine]

Finite-dimensional operators. A new line of research has suggested using
neural networks to learn mesh-free, infinite-dimensional operators
[hornik1989multilayer, khoo2021solving]. Only a forward pass of the network is
required to provide a solution for a new instance of the parameter, avoiding
the significant computational difficulties that plague Neural-FEM approaches.
The neural operator just requires data, not knowledge of the underlying PDE.
Due to the difficulty of evaluating integral operators, neural operators have
not provided effective numerical techniques that can match the success of
convolutional or recurrent neural networks in the finite-dimensional scenario.

Neural-FEM. The third method is known as the physics-informed neural network
(PINN). PINN differs from other machine learning paradigms widely employed for
mechanics and physics challenges in terms of how data is needed and employed.
Unlike supervised learning, which is often used for materials laws and
requires artificial intelligence to be trained with labels in order to
generate forecasts, the search for the solution in PINN does not require any
data other than the ones required to form the loss function
[bar2019unsupervised, karniadakis2021physics, chen2021physics,
haghighat2021physics].

PINN challenges. The training of PINN, however, is far from simple, especially
for nonlinear systems of equations. To construct the multi-layer perceptron,
non-linearities should be applied to each element of the output of the linear
transformation. This is unlike the finite element method, which is a more
entrenched framework with clear strategies and established mathematical
analysis that guarantees convergence and stability for both the solution and
weighting functions in predetermined finite-dimensional spaces. Furthermore,
for both forward and inverse problems, the physical constraints or controlling
equations could be expressed in several ways; for instance, the collocation-
based loss function, which evaluates the solution at specific collocation
points, or the energy-based method that can reduce the order of the
derivatives in governing equations despite requiring numerical integrations. A
large number of tunable hyperparameters, such as the configurations of the
neural network, the types of activation functions, and the neuron weight
initialization, as well as different techniques to impose boundary conditions
while providing significant flexibility, may bewilder researchers who are
unacquainted with neural networks [bhattacharya2020model, psaros2021meta,
fuhg2021model].

Goals. Modifying physics-informed neural network (PINN) to account for stiff
PDEs that govern Pyrolysis is the main goal of this article. Furthermore, we
proposed different strategies to reduce PINN training cost, to lower the
number of required trial-and-error tasks to tune the hyper-parameters, and to
increase the training process’s robustness. While our proposed strategies can
be used in other deep neural networks, the scope of this study is limited to
the collocation physics-informed neural network. Our current work develops a
hybrid PINN engine by
(i) employing non-dimensionalization and normalization of the physical
parameter to reduce complex equations,
(ii) infusing the loss function with our crude understanding of the solution
field or imprecise knowledge,
(iii) using an adaptive weighted-sum scheme to enhance optimization algorithms
in the context of multi-objective optimization,
(IV) exploiting different weight initialization approaches due to unbalanced
gradients

##  2 Exothermic Heat Transfer

When charring polymers are subjected to high temperatures, thermal energy is
transferred into the polymer via thermal convection. Pyrolytic gases and solid
residues will arise from the material’s decomposition. According to the degree
of pyrolysis, the polymer may be divided into three zones, as indicated in
Figure 2: The polymer decomposes into three zones: charring, pyrolysis, and
virgin material. [dennis2019combustion].

In the thermal study of heat transfer in polymers, three main types of thermal
energy transmission are often considered: conduction, convection, and
radiation. However, for the sake of simplicity, all mathematical models for
polymers cover the effect of heat conduction in the case of one-sided heating
only. The effect of external convection on heat transfer is rarely explored.
Similarly, heat radiation from a polymer is rarely taken into account. So, the
following PDE governs heat transfer in polymers, i.e., heat transfer with
internal heat generation is expressed as [bogetti1992process,
tadini2017thermal, niaki2021physics]

$\frac{\partial}{\partial t}\left(\rho C_{p}T\right)=\frac{\partial}{\partial
x}\left(k_{xx}\frac{\partial T}{\partial x}\right)+\frac{\partial}{\partial
y}\left(k_{yy}\frac{\partial T}{\partial y}\right)+\frac{\partial}{\partial
z}\left(k_{zz}\frac{\partial T}{\partial z}\right)+\dot{Q},$ (1)

where $T$ is the temperature, $C_{p}(T,\alpha)$, $k(T,\alpha)$, and
$\rho(T,\alpha)$ are the solid’s specific heat capacity, conductivity, and
density which are a function of temperature and degree of burning, and
$\dot{Q}$ is the rate of internal heat consumption. The change in thermal
energy per unit volume is represented on the left-hand side of the equation,
while the energy flux owing to conduction is represented on the right-hand
side. The through-thickness direction is defined by the x-direction, whereas
the planar directions are defined by the y- and z-directions.

#### Burning index

$\alpha\in(0,1)$ is introduced as a scalar unit that describes the pyrolysis
progress in each continuum. The relationship between $\dot{Q}$ and $\alpha$,
in particular, can be expressed as [mouritz2006modelling,
mouritz2004thermomechanical]

$\dot{Q}=-Q_{r}\frac{d\alpha}{dt},\qqref\frac{d\alpha}{dt}=g(\alpha,T)>0.$ (2)

where $Q_{r}$ is the heat of reaction generated per unit mass of polymer
during pyrolysis [henderson1987mathematical]. The evolving degree of burning,
$\alpha$ is usually governed by an ordinary differential equation(ODE) that
represents the rate of burning, $\frac{d\alpha}{dt}$, as a function of
temperature and degree of burning [henderson1985model, gibson2003modelling,
zobeiry2021physics].

Assumptions: Energy transfer via convection is considered to be insignificant
in this study, and the volatile gases created by the pyrolysis reaction are
expected to be evacuated from the material promptly and hence have no effect
on the temperature. Temperature and the stage of the breakdown reaction
influence thermal conductivity, density, and heat capacity. However, because
the change in thermal conductivity with temperature cannot be calculated
theoretically, this term must be determined experimentally throughout the
temperature range of interest.

In the space-time domain $(\bm{\mathit{x}},t)$, this results in a coupled
system of differential equations for temperature and degree of burning. The
temperature $T(\bm{\mathit{x}},t)$ and degree of burning
$\alpha(\bm{\mathit{x}},t)$ in the polymeric material are predicted by solving
this system of differential equations with initial and boundary conditions.
This work aims to solve the system of differential equations for the polymeric
system burning shown in Fig. 3.

Figure 3: A schematic of polymeric materials during burning.

??this part seems to be out of place

The boundary conditions at both ends can be set to a specific temperature by
using the following equations

$\small\left.T\right|_{x=0,y,z}=T_{bx}(t),\qquad\left.T\right|_{x=L,y,z}=T_{tx}(t),\qquad\left.T\right|_{x,y=0,z}=T_{by}(t),\qquad\left.T\right|_{x,y=L,z}=T_{ty}(t),\qquad\left.T\right|_{x,y,z=0}=T_{bz}(t),\qquad\left.T\right|_{x,y,z=L}=T_{tz}(t),\qquad$
(3)

where the subscripts $b$ and $t$ represent the bottom and top of materials,
respectively. The initial conditions are written as follows

$\left.T\right|_{t=0}=T_{0}(x)\qquad\text{and}\qquad\left.\alpha\right|_{t=0}=\alpha_{0}(x),$
(4)

where $T_{0}$ is the system’s initial temperature, which is normally assumed
to be constant in the geometry, and $\alpha_{0}$ is the polymer’s initial
degree of burning, which is assumed to be zero for unburned materials in the
entire spatial domain.

##  3 Modeling by Physics-Informed Neural Network (PINN)

###  3.1 Neural Network Architecture

Neural networks are well-known for their ability to represent information.
Based on the universal approximation theorem, any continuous function can be
arbitrarily estimated by a multi-layer perceptron containing one hidden layer
and a finite number of neurons [yarotsky2017error, berman2019survey]. While
neural networks can compactly express very complicated functions, obtaining
the precise parameters (weights and biases) required to solve a particular PDE
can be challenging [bin2021pinneik].

The bulk of solutions has used feed-forward neural networks since Raissi et
al. [raissi2019physics, haghighat2021sciann] original vanilla PINN. Some
researchers, on the other hand, have tested with several types of neural
networks to investigate their effect on the overall PINN performance.

We start by building a simple $D$-layer multilayer feed-forward neural network
comprising an input layer, $D-1$ hidden layers, and an output layer. We
suppose that the $d^{th}$ hidden layer has $N_{d}$ neurons. The previous
layer’s post-activation output
$\bm{\mathit{x}}^{d-1}\in{\mathbb{R}}^{N_{d-1}}$ is then fed into the $d^{th}$
hidden layer, and the specific affine transformation is of the form

$\mathcal{H}_{d}\left(\mathbf{x}^{d-1}\right)\triangleq\mathbf{W}^{d}\mathbf{x}^{d-1}+\mathbf{b}^{d},$
(5)

where the network weight $\mathbf{W}^{d}\in{\mathbb{R}}^{N_{d}\times N_{d-1}}$
and the bias term $\mathbf{b}^{d}\in{\mathbb{R}}^{N_{d}}$ to be learned are
both initialized using unique procedures like Xavier or He initialization
[glorot2010understanding, he2015delving].

The nonlinear activation function $\sigma(.)$ is applied component-by-
component to the current layer’s affine output $\mathcal{H}_{d}$. Furthermore,
for some regression issues, this nonlinear activation is not employed in the
output layer. As a result, the neural network may be denoted as

$\mathcal{N}(\mathbf{x};\Theta)=\left(\mathcal{H}_{D}\circ\sigma\circ\mathcal{H}_{D-1}\circ\cdots\circ\sigma\circ\mathcal{H}_{1}\right)(\mathbf{x}),$
(6)

where $\circ$ denotes the composition operator,
$\Theta=\left\\{\mathbf{W}^{d},\mathbf{b}^{d}\right\\}_{d=1}^{D}\in\mathcal{P}$
denotes the learnable parameters to be optimized later in the network, and
$\mathcal{P}$ denotes the parameter space, and $\mathcal{N}$ and
$\mathbf{x}^{0}=\mathbf{x}$ denote the network’s output and input,
respectively.

#### Activation function

DNN training performance is influenced by the activation function. Activations
like as ReLU, Sigmoid, and Tanh are frequently utilized [sun2020surrogate,
ertuugrul2018novel]. Because the activation function in a PINN framework is
evaluated using the second-order derivative, it is critical to choose it
carefully. Because most activation functions (such as Sigmoid, Tanh, and
Swish) are nonlinear around $0$, it is preferable to pick a range of [$0,1$]
rather than a wider domain when rescaling the PDE to a dimensionless form.
Furthermore, smooth activation functions such as the sigmoid and hyperbolic
tangent can be used to ensure the regularity of PINNs, allowing for accurate
calculations of PINN generalization error[mishra2021estimates,
krishnapriyan2021characterizing].

###  3.2 Automatic Differentiation

In order to solve a PDE in PINNs, it is necessary to take derivatives of the
network’s output with respect to the input. The function $\bm{u}$ can be
differentiated since it is approximated by a NN with smooth activation
function $\hat{\bm{u}}_{\Theta}$. Calculating derivatives can be done in four
ways: manually, symbolically, numerically, or automatically. When applied to
complex functions, symbolic and numerical methods such as finite
differentiation perform poorly; on the other hand, automatic differentiation
(AD) conquers many limitations such as floating-point precision errors,
numerical differentiation, and memory-intensive symbolic approaches
[wang2022and, bolte2020mathematical].

#### Nondimensionalization & Normalization of Equations

Nondimensionalization $\bar{X}=\frac{X}{X_{0}}$ will be used to remove of
physical dimensions of the governing parameters $X$ to describe the physical
process as the deviation of the parameters from their reference point $X_{0}$
which is used for normalization. This procedure allows us to simplify the loss
function and prevent the effects of measured units in dimensional analysis.
Further scaling of the non-dimensionless parameters is also used to prevent
the effects of massive fluctuation of certain quantities that are less
important relative to some appropriate unit. These units refer to quantities
intrinsic to the physics of the process rather than measured units. Consider
the dimensionless variables (denoted by $\overline{\mathrm{o}}$) as follows:

$\bar{t}=\frac{t}{t^{*}},\quad\overline{{x}}=\frac{{x}}{x^{*}},\quad\overline{{y}}=\frac{{y}}{y^{*}},\quad\bar{T}=\frac{T}{T^{*}},\quad\bar{\alpha}=\frac{\alpha}{\alpha^{*}},$
(7)

$t^{*},x^{*},y^{*},T^{*},\alpha^{*}$ are dimensionless variables that written
as

$x^{*}={L_{x}},\qquad y^{*}={L_{y}},\qquad t^{*}=\frac{\rho
C_{p}{L_{x}^{2}L_{y}^{2}}}{k_{xx}L_{y}^{2}+k_{yy}L_{x}^{2}},\qquad
T^{*}={T_{t}-T_{0}},\qquad\alpha^{*}=1.$ (8)

As a result, the partial derivatives and the divergence operator are likewise
stated in this manner [haq2020heat].

$\frac{\partial\circ}{\partial
t}=\frac{1}{t^{*}}\frac{\partial\circ}{\partial\bar{t}},\quad{\nabla}\cdot\circ=\frac{1}{x^{*}}\overline{{\nabla}}\cdot\circ,\quad{\nabla}\cdot\circ=\frac{1}{y^{*}}\overline{{\nabla}}\cdot\circ.$
(9)

The degree of burning relations and the heat equation are written in a
dimensionless form as

$\frac{1}{t^{*}}\frac{\partial}{\partial\bar{t}}\left(\rho
C_{p}\bar{T}\right)=\frac{1}{x^{*}}\frac{\partial}{\partial\bar{x}}\left(k_{xx}\frac{1}{x^{*}}\frac{\partial\bar{T}}{\partial\bar{x}}\right)+\frac{1}{y^{*}}\frac{\partial}{\partial\bar{y}}\left(k_{yy}\frac{1}{y^{*}}\frac{\partial\bar{T}}{\partial\bar{y}}\right)+\frac{1}{z^{*}}\frac{\partial}{\partial\bar{z}}\left(k_{zz}\frac{1}{z^{*}}\frac{\partial\bar{T}}{\partial\bar{z}}\right)+\dot{Q},\quad\textit{while}\quad\dot{Q}=-Q_{r}\frac{1}{t^{*}}\frac{\partial\bar{\alpha}}{\partial\bar{t}}.$
(10)

###  3.3 Model Estimation by Learning

In the PINN methodology, network training takes place by minimizing the total
loss of the network parameters $\bm{\Theta}$,

$\bm{\Theta}^{*}=\underset{\bm{\Theta}\in\mathbb{R}^{D}}{\operatorname{argmin}}\mathcal{L}_{\mathcal{T}}(\mathbf{X};\bm{\Theta}).$
(11)

An error or loss function is defined in PINNs utilizing the network’s
processed outputs and derivatives based on the equations guiding the problem’s
physics. As a result, the network’s total loss, $\mathcal{L_{T}}$, is made up
of the sum of loss terms for the PDE ($\mathcal{L_{\mathcal{F}}}$) and the
initial and boundary conditions ($\mathcal{L_{\mathcal{B}}}$). Also, let us
assume that we have some imprecise knowledge which can provide a general idea
of the solution field. By updating the loss function, the imprecise knowledge
can be easily included in the total loss function to guide the training
process as

$\mathcal{L_{T}}(\Theta)=\mathcal{L_{\mathcal{F}}}(\Theta)+\mathcal{L_{\mathcal{B}}}(\Theta)+\sum_{i=1}^{N_{back}}\eta_{i}\mathcal{L}_{i},$
(12)

where $\mathcal{L}_{i}$ is the loss term for the $i_{th}$ background
knowledge, and $\eta_{i}$ is the control term for each knowledge’s
contribution. Since the imprecise knowledge is not exact, the minimizer is far
away from any true optimal point. As a result, we are interested in an
adaptive optimization problem in which imprecise knowledge is originally
integrated into a supervised learning assignment. Setting $\eta_{i}=0$ turns
off the associated objective once the solution is close enough to finish the
supervised learning work.. $\mathcal{L_{\mathcal{F}}}$ represents due to error
in satisfying the governing differential equations $\mathcal{F}$. It imposes
the differential equation $\mathcal{F}$ at collocation points over the domain,
$\Omega$, which can be chosen uniformly or unevenly. The other term,
$\mathcal{B}\left(\hat{\bm{u}}_{\Theta}\right)=\bm{g}$, represents due to the
error in satisfying of the boundary or initial conditions. A mean square error
formulation is used in a common implementation of the loss, where

$\mathcal{L}_{\mathcal{F}}(\Theta)={{MSE}}_{\mathcal{F}}=\frac{1}{N_{c}}\sum_{i=1}^{N_{c}}\|\mathcal{F}\left(\hat{\bm{u}}_{\Theta}\left(\bm{z}_{i}\right)\right)-\bm{f}\left(\bm{z}_{i}\right))\|^{2}=\frac{1}{N_{c}}\sum_{i=1}^{N_{c}}\|r_{\Theta}\left(\bm{u}_{i}\right)-\bm{r}_{i}\|^{2},$
(13)

and

$\mathcal{L}_{\mathcal{B}}(\Theta)=MSE_{\mathcal{B}}=\frac{1}{N_{B}}\sum_{i=1}^{N_{B}}\|\mathcal{B}\left(\hat{\bm{u}}_{\Theta}(\bm{z})\right)-\bm{g}\left(\bm{z}_{i}\right))\|^{2}.$
(14)

Figure 4: A schematic of PINN blocks. Differential equation residual (loss)
terms, as well as initial and boundary conditions, make up PINNs.

Optimization is used to minimize the loss function. Based on the literature,
optimization of loss functions is performed using minibatch sampling with the
Adam and LBFGS method, which is a quasi-Newton optimization procedure in most
of the PINN literature. The gradient-based optimizer will almost likely become
stuck in one of the local minima for the loss function [pang2019fpinns,
maclaurin2015gradient]. Because stochastic gradient descent (SGD) struggles
with random collocation points, especially in 3D setups, the Adam approach,
which combines adaptive learning rate and momentum approaches, is used to
speed up convergence [zhu2021machine, bock2019proof].

Note that because it may incorporate numerous different objectives for
separate tasks, the PINN problem might be regarded as a multi-objective multi-
task problem. This view of modeling has several advantages. It lowers the
order of regularity, which is advantageous for problems involving
discontinuities. Additionally, because the order of partial derivatives
increases exponentially in terms of time complexity in the automated
differentiation technique, it may lower the total computing cost
[zhu2021local]. Multiple objectives used to assess the correctness of the
solution may result in gradients that are conflicting. As a result of the
opposing gradient, a compromised solution that is not always Pareto optimal
may emerge. We use the gradient surgery technique [yu2020gradient] for the
PINN framework to address the conflicting gradient and the challenge of
balancing different objectives. We have the following update rule in the most
common form of GD-based methods,

$\Theta_{n}=\Theta_{n-1}-\alpha_{n}\nabla_{\Theta}{\mathcal{L}_{\mathcal{T}}}_{n-1},$
(15)

where $\Theta_{n}$ is the updated unknown vector at iteration $n^{th}$,
$\nabla_{\Theta}{\mathcal{L}_{\mathcal{T}}}_{n-1}$ is the gradient of the
total loss function with respect to the unknown vector at iteration $n-1$, and
$\alpha_{n}$ is the learning rate that regulates the solution’s stability
during iterations.

Let’s call each objective’s gradient vector
$\bm{\mathit{g}}_{i}=\nabla_{\Theta}{\mathcal{L}_{i}}(\Theta)$ in gradient
surgery, and the total gradient is
$\bm{\mathit{g}}=\nabla_{\Theta}{\mathcal{L}_{T}}(\Theta)=\sum_{i=1}^{N_{obj}}\bm{\mathit{g}}_{i}$.
You can find the concept in Fig. 5. For further information, we direct
interested readers to [yu2020gradient].

Algorithm 1 Algorithm of gradient surgery

1:$\mathbf{\Theta}$ $\leftarrow$ Initialize random weights and biases.

2:$n,i$ $\leftarrow 1$

3:for $n$ to $N_{itr}$ do

4: for $i$ to $N_{obj}$ do

5: $\mathcal{L}_{i}(\mathbf{\Theta_{n}})$,
$\bm{\mathit{g_{i}}}=\nabla\mathcal{L}_{i}$ $\leftarrow$ compute objectives
and gradients.

6: if $\bm{\mathit{g_{i}}}.\bm{\mathit{g_{j}}}<0$ then

7:
$\bm{\mathit{g_{i}}}=\bm{\mathit{g_{i}}}-\frac{\bm{\mathit{g_{i}}}.\bm{\mathit{g_{j}}}}{||\bm{\mathit{g_{i}}}||^{2}}$
$\leftarrow$ modified gradient.

8: end if

9: $\bm{\mathit{g}}=\sum_{i=1}^{N_{obj}}\bm{\mathit{g}}_{i}$ $\leftarrow$
Aggregate gradient vectors.

10: end for

11: $\mathbf{\Theta}_{n+1}$ $\leftarrow$
$GD(\mathbf{\Theta}_{n},\bm{\mathit{g}})$

12:end for

Figure 5: Gradient surgery idea [yu2020gradient, bahmani2021training]

#### Weight initialization

Because training deep models is such a challenging operation, the choice of
initialization has a significant impact on most algorithms. The beginning
point can decide whether or not the method converges at all, with certain
initial points being so insecure that the algorithm runs into arithmetic
issues and fails. Over the previous decade, more customized techniques have
been the de facto norm, as they may result in a slightly more successful
optimization process. Pang et al. [pang2019fpinns] describe a method for
picking the most appropriate one.

For initialization, Glorot and Bengio [glorot2010understanding] advocated
using a correctly scaled uniform distribution. This is referred to as the
”Xavier” initialization. The assumption that the activations are linear is
used to derive it. This approach is computed as a random number with a uniform
probability distribution ($\mathcal{U}$) as
”$weight=\mathcal{U}[-\frac{1}{\sqrt{n}},\frac{1}{\sqrt{n}}]$”, where $n$ is
the number of node inputs. When used to initialize networks that employ the
rectified linear (ReLU) activation function, the ”Xavier” weight
initialization was discovered to have issues.

The ”He” initialization method is the current standard for initializing the
weights of neural network layers and nodes that employ the rectified linear
(ReLU) activation function [he2015delving]. The he initialization technique is
computed as a random number with a Gaussian probability distribution
($\mathcal{G}$) as $weight=\mathcal{G}(0,\sqrt{\frac{2}{n}})$, where $n$ is
the number of node inputs.

#### Adaptive weights

Training multi-objective total loss functions pose challenges for optimization
techniques, as others have pointed out [raissi2019physics]. The unique
solution to the governing equations is obtained once initial and boundary
conditions are imposed strongly, especially in the case of boundary value
problems. It is common to use a weighted-sum scheme to enhance optimization
algorithms in the context of multi-objective optimization
[marler2010weighted]. Wang et al. [wang2020understanding] showed that
unbalanced gradients cause inaccuracy in optimization and developed an
adaptive loss weight approach where gradients of individual terms in the loss
function are normalized in order to decrease the stiffness of the gradient
dynamics in PINNs. As a result, the total loss term $\mathcal{L_{Tot}}$ is
changed to

$\mathcal{L_{Tot}}=\omega_{\mathcal{F}}\mathcal{L_{\mathcal{F}}}+\omega_{\mathcal{B}}\mathcal{L_{\mathcal{B}}},$
(16)

where $\omega$ represents the weight loss associated with each loss term.
Based on the parameter values $(\Theta)^{e}$ at epoch $e$, the updated scaling
weight $\omega^{e+1}$ for each loss term is calculated as

$\omega^{e+1}=\beta\omega^{e}+(1-\beta)\hat{\omega}^{e+1},$ (17)

where $\beta=0.9$ based on [wang2020understanding] and

$\hat{\omega}_{i}^{e+1}=\frac{1}{\omega_{i}^{e}}\frac{\max\left(\left|\nabla\mathcal{L}_{T}\left(\bm{\Theta}^{e}\right)\right|\right)}{\operatorname{mean}\left(\left|\nabla\mathcal{L}_{i}\left(\bm{\Theta}^{e}\right)\right|\right)}$
(18)

Such an adaptive strategy for loss weights improves the system’s robustness
and improves the model’s prediction accuracy. ?? it is not clear which method
has been used for adaptive weight optimization

Algorithm 2 PINN Algorithm for Burning Problem

1:$\textbf{X},\textbf{t}$ $\leftarrow$ Sample uniformly spatial and temporal
domains.

2:$\mathbf{\Theta_{n}}$ $\leftarrow$ Initialize randomly using Xavier scheme.

3:$n$ $\leftarrow 1$

4:while $err>TOL$ do

5: $\mathbf{\Theta_{n}}$ $\leftarrow$ Optimize $\mathcal{L}_{T}$ for
$\mathbf{\Theta}$.

6: $T_{n},\alpha_{n}$ $\leftarrow$ Evaluate temperature and degree of burning
using $\mathbf{\Theta}$.

7: $err$ $\leftarrow$
$\left\|\mathbf{{\Theta}_{n}}-\mathbf{{\Theta}_{n-1}}\right\|/\left\|\mathbf{{\Theta}_{n}}\right\|$
with $\left\|o\right\|$ as the $L^{2}$ norm.

8: $n$ $\leftarrow$ $n+1$

9:end while

##  4 Case Study 1: 1D Burning Problem of Polymer

To evaluate the performance of the proposed PINN, the heat transfer PDEs in a
polymeric material are trained in Python (V3.6.8), using Tensorflow and Keras
libraries (V2.10). We can define the PINN solution approach for the 1D
pyrolysis problem now that we have covered the method of Physics-Informed
Neural Networks (PINNs) and heat equations.

$\frac{1}{t^{*}}\frac{\partial}{\partial\bar{t}}\left(\rho
C_{p}\bar{T}\right)=\frac{1}{x^{*}}\frac{\partial}{\partial\bar{x}}\left(k_{xx}\frac{1}{x^{*}}\frac{\partial\bar{T}}{\partial\bar{x}}\right)+\dot{Q},\quad\textit{while}\quad\dot{Q}=-Q_{r}\frac{1}{t^{*}}\frac{\partial\bar{\alpha}}{\partial\bar{t}}.$
(19)

$T$ and $\alpha$ are the unknown solution variables in the case of 1D burning.
For the dimensionless version of these variables, neural networks are,

$\bar{T}:(\bar{x},\bar{t})\mapsto\mathcal{N}_{\bar{T}}\left(\bar{x},\bar{t};\bm{\Theta}_{\bar{T}}\right),\qquad\textit{and}\qquad\bar{\alpha}:(\bar{x},\bar{t})\mapsto\mathcal{N}_{\bar{\alpha}}\left(\bar{x},\bar{t};\bm{\Theta}_{\bar{\alpha}}\right),$
(20)

where $\bm{\Theta}_{o}$ denotes that these networks have their own set of
parameters. After that, the total coupled loss function is written as,

$\mathcal{L}_{Tot}=\omega_{1}\mathcal{L}_{{T}}+\omega_{2}\mathcal{L}_{\alpha}+\omega_{3}\mathcal{L}_{\alpha_{0}}+\omega_{4}\mathcal{L}_{T_{0}}+\omega_{5}\mathcal{L}_{T_{BC1}}+\omega_{6}\mathcal{L}_{T_{BC2}}+\sum_{i=1}^{N_{back}}\eta_{i}\mathcal{L}_{i},$
(21)

while

$\displaystyle\mathcal{L}_{T}$
$\displaystyle=\frac{1}{N_{c}}\sum_{i=1}^{N_{c}}\|\frac{1}{t^{*}}\frac{\partial}{\partial\bar{t}}\left(\rho
C_{p}\bar{T}\right)-\frac{1}{x^{*}}\frac{\partial}{\partial\bar{x}}\left(k_{xx}\frac{1}{x^{*}}\frac{\partial\bar{T}}{\partial\bar{x}}\right)+Q_{r}\frac{1}{t^{*}}\frac{\partial\bar{\alpha}}{\partial\bar{t}}\|^{2},$
(22) $\displaystyle\mathcal{L}_{\alpha}$
$\displaystyle=\frac{1}{N_{c}}\sum_{i=1}^{N_{c}}\|\frac{1}{t^{*}}\frac{\partial\bar{\alpha}}{\partial\bar{t}}-g(\bar{\alpha},\bar{T})\|^{2},\qquad\textit{loss
of burning degree}$ $\displaystyle\mathcal{L}_{T_{0}}$
$\displaystyle=\frac{1}{N_{B}}\sum_{i=1}^{N_{B}}\|\bar{T}-\bar{T}_{0}(x)\|^{2},\qquad\textit{loss
of initial condition of temperature}$ $\displaystyle\mathcal{L}_{T_{bc_{1}}}$
$\displaystyle=\frac{1}{N_{B}}\sum_{i=1}^{N_{B}}\|\bar{T}-\bar{T}_{b}(t)\|^{2},\qquad\textit{loss
of boundary condition at bottom}$ $\displaystyle\mathcal{L}_{T_{bc_{2}}}$
$\displaystyle=\frac{1}{N_{B}}\sum_{i=1}^{N_{B}}\|\bar{T}-\bar{T}_{t}(t)\|^{2},\qquad\textit{loss
of boundary condition at top}$ $\displaystyle\mathcal{L}_{\alpha_{0}}$
$\displaystyle=\frac{1}{N_{B}}\sum_{i=1}^{N_{B}}\|\bar{\alpha}-\bar{\alpha}_{0}(x)\|^{2}\qquad\textit{loss
of initial condition of burning degree},$

#### Training and Hyperparameter Searches

It is hardly surprising that network size, architecture, and optimizer
hyperparameters like learning rate can significantly impact PINN solution
quality. We used the hyperbolic-tangent and sigmoid activation functions to
create neural networks with four hidden layers and 20 neurons in each layer
for all of the examples in this paper. We used the Adam optimizer with an
initial learning rate of $2*10^{-3}$ and an exponential learning decay to
$10^{-5}$ at the end of training. We have also explored a range of
hyperparameters, but we have not noticed any substantial improvements in terms
of improving the pyrolysis problem studied here. Fig. 6 shows the training
history for various hidden layers and neurons.

Figure 6: Losses for a variety of hidden layer (left) and neuron numbers
(right).

#### Comparison with Numerical Results

The reference numerical model results will be compared to the above-mentioned
best-performing PINN model outputs in this part. Fig. 7 shows the PINN,
numerical solution, and errors. As you can see, the PINN solution has a high
level of agreement with the expected solution when compared to the findings of
numerical analysis. The upper face is given a temperature of $T_{t}=700$,
while the bottom face is given a temperature of
$T_{t}=700sin(\frac{\pi}{20}t)$. $T_{0}=700x^{3}$ is the initial temperature.
??please discuss how you calculated the numerical, is it matlab solver

Figure 7: PINN and numerical predictions for full cycle of burning for
temperature (left), and degree of burning (right). ?? I would suggest make t
like a table, one row temp, other one alpha, colomns can be PINN, numerical,
and error

#### Role of Weight Initialization

Because training deep models is such a tough operation, the choice of
initialization has a significant impact on most algorithms. So, we trained the
model with Xavier and He initialization for tanh and ReLU activation functions
hereunder. Fig. 8 shows the results for Xavier and He initialization.

Figure 8: Evolution of loss term for Xavier weight initialization and He
weight initialization.

Training with and without Adaptive Weights. The performance of the suggested
PINN with and without loss weights is shown in Fig. 9 to emphasize the need to
use loss weights.

Figure 9: Evolution of loss term for optimization with and without loss
weights.

Training with and without Background Knowledge. To show the ability of
background knowledge, we compare two case studies with and without background
knowledge. We consider that the top of the geometry of polymer has a higher
temperature and degree of burning always compared to the bottom of that. We
considered it as the background knowledge and added it to the total loss. Fig.
10 shows the result for comparing of this factor.

Figure 10: Evolution of loss term for optimization with and without background
knowledge.

##  5 Case Study 2: 2D Burning Problem of Matrix Cube

For the second example, we modeled the pyrolysis of ethylene propylene diene
monomer (EPDM) composites at the high temperature and pressure of solid rocket
motors (SRMs), which results in a char layer with a non-uniform distribution
pattern. Interestingly, the compact structure of the char can increase the
ablation resistance of the composites. To better understand the ablation
mechanism and inform the design of EPDM composites for SRMs, it is crucial to
understand and describe the char formation pattern throughout the life of the
material and also understand the material behaviour at the interface. We used
our H-PINN to predict the spatial evolution of chemical processes (pyrolysis
and charring) before validating them against mesoscale finite element results.
It is of great theoretical and practical significance to study the formation
mechanism of the compact structure of the char layer for further studying the
ablation mechanism and guiding the development of EPDM composites. To date,
the char formation pattern in polymeric materials cannot be predicted with
reduced order models, such as machine-learned approaches. In this work, the
hybrid PINN for 2D pyrolysis of EPDM composites are designed from the basic 2D
heat equation, namely

$\frac{1}{t^{*}}\frac{\partial}{\partial\bar{t}}\left(\rho
C_{p}\bar{T}\right)=\frac{1}{x^{*}}\frac{\partial}{\partial\bar{x}}\left(k_{xx}\frac{1}{x^{*}}\frac{\partial\bar{T}}{\partial\bar{x}}\right)+\frac{1}{y^{*}}\frac{\partial}{\partial\bar{y}}\left(k_{yy}\frac{1}{y^{*}}\frac{\partial\bar{T}}{\partial\bar{y}}\right)+\dot{Q},\quad\textit{while}\quad\dot{Q}=-Q_{r}\frac{1}{t^{*}}\frac{\partial\bar{\alpha}}{\partial\bar{t}}.$
(23)

Figure 11: SEM image of vertical-section of the char layer (left) and Micro-CT
images of horizontal-section of the char layer (right) [xi2019redeposition].

The following steps were required to expand the PINN to 2D: adding one input
(y dimension) and providing a pre-layer for the y input, adjusting error
calculation to include a term for the second derivative of the PINN prediction
with respect to y, and modifying the error calculation for the heat transfer
boundary conditions. The total loss function’s many terms, as well as their
complexity, present huge complications in network training. The only way to
get meaningful results from the Adam optimizer was to utilize scored-adaptive
weights; however, we have found that the process has a high computational cost
and lacks robustness. We offer a unique strategy, inspired by
[pirk2017interactive], to overcome the limitations of standard network
training methods for multiphysics issues.

Figure 12: A mesh is used to discretize the material (left), the pyrolyzing
front proceeds towards the last spot of burning as pyrolysis progresses
(right).

The pyrolyzing front is shown as a two-dimensional triangle element. It means
that we reduce the order of 2D problems by decomposing them into many 1D
problems and training them as a collection of learning. In this concept, we
need to assume the last point of burning in the material. Next, we decompose
2D geometry into many 1D elements, which the last point of burning is shared
among all 1D elements. Fig. 13 and Fig. 14 show the PINN solution for
temperature and degree of burning for different time steps, respectively. The
domain is subjected to a constant temperature $300+700e^{-15x^{2}}$ at the top
and bottom faces and a constant temperature $300+600e^{-15y^{2}}$ at the right
and left faces.

Figure 13: Temperature predictions using PINN for various time steps in
burning $t=0.01,0.1,0.5,2,4,7$.

Figure 14: Degree of burning predictions using PINN for various time steps in
burning $t=0.01,0.1,0.5,2,4,7$.

##  6 Case Study 3: Shen’s Model for Burning

The second category of burning, which considers that the specific heat remains
constant, with only thermal conductivity and material density being
temperature dependent, is presented in this section [shen2007modeling]. The
heat transfer PDE in a polymeric material is presented as

$\frac{\partial}{\partial
t}\left(\rho_{s}C_{s}T+\rho_{c}C_{c}T\right)=\frac{\partial}{\partial
x}\left(k_{xx}\frac{\partial T}{\partial
x}\right)+\dot{Q},\quad\text{while}\quad\dot{Q}=K\rho_{s}(C_{s}-C_{c})(T).$
(24)

Figure 15: Schematic of 1D burning and char formation.

Results are presented for boundary and initial conditions such that $L=1mm$,
$T_{0}=700x^{3}$, $T_{b}=0$ and $T_{t}=700$ for duration of $10$ seconds. The
training is conducted on $500$ uniformly distributed points in $(x,t)$.

Thermal characteristic values of polymers

The thermal conductivity and density of the polymer, as well as char, have
temperature-dependent properties, as shown in Fig. 16. The thermal
conductivity of the material is shown in Fig. 16.a. When the temperature rises
above $250F^{\circ}$, the sample’s exposure time to heat flux increases, and
the material’s thermal conductivity rapidly decreases as a result of the
temperature rise. The material’s thermal conductivity before this temperature
is around $10W/mK$. The thermal conductivity of the solid stabilizes at
$0.3W/mK$ at temperatures up to $550F^{\circ}$. Beyond $250$ degrees
Fahrenheit, the density of the polymer rapidly declines until it reaches zero
at $550$ degrees Fahrenheit. The density of the char, on the other hand,
rapidly increases to $800kg/m3$ at temperature $700$ (see Fig. 16.b).

Figure 16: Thermo-physical properties of polymer

Comparison between numerical and predicted temperature profiles

To assess the accuracy of the PINN in terms of predicting $T$, results are
compared to numerical predictions from MATLAB. Fig. 17 shows the temperature
prediction for the case study. The obtained relative error of the predicted
temperature when comparing PINN and numerical methods is below $1\%$.

Figure 17: PINN and numerical temperature predictions for full cycle of
burning.

The temperature evolution at the left side of the polymeric material, middle
of the polymeric material, and right side of the polymer are shown for the
full burning cycle in Fig. 18. Also, The mean square error of the loss
function is reported in Fig. 19.

Figure 18: Temperature history for (x = 0.2), middle of polymer (x = 0.5), and
for boundary (x = 0.8).

Figure 19: Evolution of loss term in the training of the network.

It is somewhat remarkable that a neural network with hundreds of parameters
cannot approximate the analytical solution of the 1D hyperbolic PDE with a
non-convex flux function with any degree of accuracy. This is particularly
striking because, in theory, there should exist a network that can offer a
near approximation of the continuous solution for any arbitrarily chosen PDE.
However, this is not what is seen. As a result, this leads us to conclude that
the issue is not with the solution in nature but with the methodology we
employ to obtain it, i.e., the optimization procedure or the loss function.

##  7 Discussion and Concluding Remarks

In this paper, we studied a comprehensive study on Physics-Informed Neural
Networks (PINNs) for the forward solution of pyrolysis problems by making the
training more straightforward. We explored how Physics-Informed Neural
Networks (PINNs) can be employed to solve pyrolysis problems in the forward
phase. Our research is the first to explore fully coupled temperature-degree-
of-burning relationships in pyrolysis problems. We presented a dimensionless
version of these relations that leads to the optimizer’s stable and convergent
behavior.

While the achieved results are close to our expectations, it should be noted
that training PINNs is time-consuming. We relate the training challenge to the
multi-objective optimization issue and the application of a first-order
optimization algorithm, as reported by others. Given the difficulties
encountered and overcome in this work for the forward problem, the next step
is to use PINNs to inverse burning situations.

A multi-task learning approach has emerged in which a NN must fit observed
data while decreasing a PDE residual. This article introduces single- and
separated-network PINN architectures to forecast temperature distributions and
the degree of burning of a pyrolysis problem in a one-dimensional (1D) and
two-dimensional (2D) rectangular domain for the first time to model through
collocation training.

The complex and non-convex multi-objective loss function presents substantial
obstacles for forward problems in training PINNs. We discovered that adding
several differential relations to the loss function causes an unstable
optimization issue, which may lead to convergence to the trivial null solution
or significant deviation of the solution. The dimensionless form of the
coupled governing equations that we find most beneficial to the optimizer is
used to address this problem. The numerical results are compared with results
obtained from PINN to show the performance of the solution. Our research is
the first to explore fully coupled temperature-degree-of-burning relationships
in pyrolysis problems. Unlike classical numerical methods, the proposed PINN
does not depend on domain discretization. In addition to these
characteristics, the proposed PINN achieves good accuracy in predicting
solution variables, which makes it a candidate to be utilized for surrogate
modeling of pyrolysis problems. While the achieved results are close to our
expectations, it should be noted that training PINNs is time-consuming.

##  Appendix A General PDEs

In the case of a function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ for
which we want to calculate the Jacobian $J$, AD can calculate the numerical
derivative in either forward or reverse mode once all operations’ graphs that
make up the mathematical expression are calculated. As a result, PINNs can
solve differential equations in their most general form, such as

$\displaystyle\mathcal{F}(\bm{u}(\bm{z});\gamma)$
$\displaystyle=\bm{f}(\bm{z})$ $\displaystyle\bm{z}\text{ in }\Omega,$ (25)
$\displaystyle\mathcal{B}(\bm{u}(\bm{z}))$ $\displaystyle=\bm{g}(\bm{z})$
$\displaystyle\bm{z}\text{ in }\partial\Omega,$

defined with the boundary $\partial\Omega$ on the domain $\Omega$.
$\bm{z}=\left[x_{1},\ldots,x_{n};t\right]$ denotes the space-time coordinate
vector, $\bm{u}$ denotes the unknown solution, $\gamma$ are the physics
parameters, $\bm{f}$ denotes the function identifying the problem’s data, and
$\mathcal{F}$ denotes the non linear differential operator. Finally, because
the initial condition is a sort of Dirichlet boundary condition in the spatio-
temporal domain, $\mathcal{B}$ can be used to signify arbitrary initial or
boundary conditions relating to the problem, and $\bm{g}$ can be used to
denote the boundary function.

In essence, AD includes a PDE into the loss of the neural network, where the
differential equation residual is

$r_{\Theta}(\bm{z}):=\mathcal{F}\left(\hat{\bm{u}}_{\Theta}(\bm{z});\gamma\right)-\bm{f}.$
(26)

The following is the result of the original formulation of the aforementioned
equation.

$r_{\Theta}(\bm{x},t):=\frac{\partial}{\partial
t}\hat{u}_{\Theta}(\bm{x},t)+\mathcal{F}_{\bm{x}}\hat{u}_{\Theta}(\bm{x},t).$
(27)

In the deep learning framework, the chain rule is employed to calculate
hierarchical derivatives from the output layer to the input layer with regard
to input spatio-temporal coordinates using automatic differentiation for
numerous layers.

##  Appendix B Mapping

We’ll discuss how to map a circular region to a square region and back in this
part. There are countless variations on how to carry out this mapping. We
focus on mappings that have attractive closed-form invertible equations. In
terms of mathematics, we are looking for functions $f$ that transform each
point $(u,v)$ in the circular disc to a point $(x,y)$ in the square region and
vice versa. To put it another way, we’re trying to come up with equations for
$f$ such that $(u,v)=f(x,y)$ and $(x,y)=f^{-1}(u,v)$. Additionally, we will
subject the mapping to the following three constraints:

  * 1.

$(0,0)=f(0,0)$, meaning that the center of the circle and the center of the
square, respectively.

  * 2.

The extreme points of the forms all along x-axis match up, or $(\pm 1,0)=f(\pm
1,0)$.

  * 3.

The extreme points of the shapes on the y-axis correspond because $(0,\pm
1)=f(0,\pm 1)$.

One of the easiest methods for mapping a circular disc to a square area is to
linearly extend the circle to its rim. Depending on which side of the square
the stretching takes place, the equations for extending from rim to rim must
be divided into four separate situations. We may further reduce these
equations by applying the signum function, and a technique developed by Dave
Cline [shirley1997low] to get:

$u=\left\\{\begin{array}[]{ll}\operatorname{sgn}(x)\frac{x^{2}}{\sqrt{x^{2}+y^{2}}}&\text{
when }x^{2}\geq y^{2}\\\
\operatorname{sgn}(y)\frac{xy}{\sqrt{x^{2}+y^{2}}}&\text{ when
}x^{2}<y^{2}\end{array}\quad
v=\begin{cases}\operatorname{sgn}(x)\frac{xy}{\sqrt{x^{2}+y^{2}}}&\text{ when
}x^{2}\geq y^{2}\\\
\operatorname{sgn}(y)\frac{y^{2}}{\sqrt{x^{2}+y^{2}}}&\text{ when
}x^{2}<y^{2}\end{cases}\right.$ (28)

Also, Manuel Fernandez Guasti developed an algebraic equation in 1992 to
depict a form that is between a circle and a square [guasti1992analytic]. The
parameter $s$ in his equation may be utilized to nicely mix the circle and the
square. We can create a method to seamlessly transfer a circular disc to a
square area using the Fernandez Guasti squircle. The key concept is to map
each circular contour on the inside of the disc to a square circular contour.
This may be combined with the radial restriction, and the resulting equations
for the mapping can be obtained.

Figure 20: Various type of mapping.

On the other hand, a regular rectangular grid is transformed into a regular
curvilinear grid made up of elliptical arcs via elliptical grid mapping. The
inverse of Nowell’s elliptical grid mapping is now obtained. The three steps
of the derivation can be summed up as follows: Step 1: Introduce the
trigonometric variables $\alpha$ and $\beta$ and convert circular coordinates
$(u,v)$ to polar coordinates. Step2: Obtain an equation for $(x,y)$ in terms
of $\alpha$ and $\beta$ in step two. Find an equation for $(x,y)$ in terms of
$u$ and $v$ in step three.