From 07fe0d5599a51a3cc766a2c921e3daf538a33c73 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 24 Jun 2024 17:48:59 -0600 Subject: [PATCH] Remove now no longer used directory. --- .../Non-linear_heat_conduction.tex | 511 ------------------ .../Documentation/refined.sty | 279 ---------- 2 files changed, 790 deletions(-) delete mode 100755 nonlinear-heat_transfer_with_AD_NOX/Documentation/Non-linear_heat_conduction.tex delete mode 100755 nonlinear-heat_transfer_with_AD_NOX/Documentation/refined.sty diff --git a/nonlinear-heat_transfer_with_AD_NOX/Documentation/Non-linear_heat_conduction.tex b/nonlinear-heat_transfer_with_AD_NOX/Documentation/Non-linear_heat_conduction.tex deleted file mode 100755 index f1b6292..0000000 --- a/nonlinear-heat_transfer_with_AD_NOX/Documentation/Non-linear_heat_conduction.tex +++ /dev/null @@ -1,511 +0,0 @@ -\documentclass[preprint,12pt]{elsarticle} -\usepackage{refined} -\usepackage{amsmath} -\usepackage{listings} -\usepackage{fancyvrb} -\usepackage[makeroom]{cancel} -\newcommand{\te}[1]{\lstinline{#1}} - -\usepackage{xcolor} - -\definecolor{codegreen}{rgb}{0,0.6,0} -\definecolor{codegray}{rgb}{0.5,0.5,0.5} -\definecolor{codepurple}{rgb}{0.58,0,0.82} -\definecolor{backcolour}{rgb}{0.95,0.95,0.96} - -\lstdefinestyle{mystyle}{ - backgroundcolor=\color{backcolour}, - commentstyle=\color{codegreen}, - keywordstyle=\color{magenta}, - %numberstyle=\tiny\color{codegray}, - stringstyle=\color{codepurple}, - basicstyle=\ttfamily\small, - breakatwhitespace=false, - breaklines=true, - captionpos=b, - keepspaces=true, - numbersep=7pt, - showspaces=false, - showstringspaces=false, - showtabs=false, - tabsize=1, - classoffset=0, - xleftmargin=0pt, - frameround=tttt, - language=C++, - columns=flexible -} - -\lstset{style=mystyle,emph={ADNumberType},emphstyle={\color{red}}} -\journal{Elsevier} -\numberwithin{equation}{section} -\title{Nonlinear heat transfer problem} -\author{Narasimhan Swaminathan} -\affiliation{Department of Mechanical Engieeering, Indian Institute of Technology Madras, Chennai, 600036} - -\begin{document} - \begin{abstract} - In this example, we solve a simple transient nonlinear heat transfer equation. The nonlinearity is due to the temperature dependence of the thermal conductivity. Two main aspects covered by this example are (a) it develops the residual and the jacobian using automatic differentiation and (b) solves the nonlinear equations using TRILINOS NOX. The actual code contains the comments which will explain how these aspects are executed. Here, we give the full derivation and set up the equations. We also provide explanations to some of the functions important for this applicaiton. - \end{abstract} - \maketitle -\section{Strong form} -\subsection{Governing equations} -We solve now a non-linear differential equation governing the heat transfer in a 2D domain using a combination of Automatic Differentiation and TRILINOS NOX in this example. The link for the original problem can be found in \href{https://www.mathworks.com/help/pde/ug/heat-transfer-problem-with-temperature-dependent-properties.html}{Link}. We shall consider the finite element model of the following differential equation -\begin{equation} - \rho C_p \pd{T(x,y)}{t} = \nabla \cdot k \nabla T + f \label{eq:sf1} -\end{equation} -subjected to appropriate boundary conditions. For the sake of convenience, we will write the \eref{eq:sf1} as -\begin{equation} - \rho C_p \dot{T} = \pl k T_{,i} \pr_{,i} + f \label{eq:sf2} -\end{equation} -The nonlinearity arises because the thermal conductivity is a function of the temperature $T$. We will write down the functional form of this dependence later. -\subsection{Boundary conditions and domain} -Please see the \href{https://www.mathworks.com/help/pde/ug/heat-transfer-problem-with-temperature-dependent-properties.html}{Link} for the actual domain over which the differential equation in solved. The \te{.geo} and \te{.msh} file for the domain can be found in the \te{mesh} folder. The boundary conditions are as follows: Left end (domain id 1) is fixed at 100$^\circ$C and the right end (domain id 3) has a outgoing flux of magnitude 10. All other material parameters are defined in file \te{source/nonlinear_heat_cons_des.cc} -\section{Weak form} -The weak form of the \eref{eq:sf2} over a typical finite element $\Omega_e$ with boundary $\Gamma_e$ is written as -\begin{equation} - \int_{\Omega_e} w \pl\rho C_p \dot{T}- \pl k T_{,i} \pr_{,i} - f \pr dx dy = 0\label{eq:wf1} -\end{equation} -which becomes -\begin{equation} - \int_{\Omega_e} w \rho C_p \dot{T}- \pl w \pl k T_{,i} \pr_{,i} + w f \pr dx dy = 0 \label{eq:wf2} -\end{equation} -integrating by parts we get -\begin{equation} - \int_{\Omega_e} w \rho C_p \dot{T} dx dy - \int_{\Omega_e} \pl -w_{,i} \pl k T_{,i} \pr + w f \pr dx dy -\int_{\Gamma_e} w kT_{,i}n_{i}ds^e = 0 \label{eq:wf3} -\end{equation} -Clearly, $k T_{,i}n_i$ is nothing but the heat flux along the boundary $\Gamma_e$, which we write as $q_{\Gamma^e}$. $ds^e$ is an infinitesimal element along the boundary. Hence, we have -\begin{equation} - \int_{\Omega_e} w \rho C_p \dot{T} dx dy - \int_{\Omega_e} \pl -w_{,i} \pl k T_{,i} \pr + w f \pr dx dy -\int_{\Gamma_e} w q_{\Gamma^e}ds^e = 0 \label{eq:wf4} -\end{equation} -We now consider the approximation for $T = \psi_J T_J$, where summation is implied over the repeated indices (Summation over all nodes in the finite element). $\psi_J$ is the shape function. We hence get the finite element model taking ($w=\psi_I$) as -\begin{equation} -\dot{T}_J \int_{\Omega_e}\rho C_p \psi_I \psi_J dx dy + \int_{\Omega_e} \pl \psi_{I,i} \pl k \psi_{J,i}T_J \pr - \psi_I f \pr dx dy -\int_{\Gamma_e} \psi_I q_{\Gamma^e}ds^e = 0 \label{eq:wf5} -\end{equation} -we note the $k$ is now also a function of $T$, hence we actually have -\begin{equation} - \dot{T}_J \int_{\Omega_e}\rho C_p \psi_I \psi_J dx dy + \int_{\Omega_e} \pl \psi_{I,i} \pl k (\psi_PT_P) \psi_{J,i}T_J \pr - \psi_I f \pr dx dy -\int_{\Gamma_e} \psi_I q_{\Gamma^e}ds^e = 0 \label{eq:wf6} -\end{equation} -calling $\int_{\Gamma_e} \psi_I q_{\Gamma^e}ds^e$ as $Q_I$, we have -\begin{equation} - \dot{T}_J \int_{\Omega_e}\rho C_p \psi_I \psi_J dx dy + \int_{\Omega_e} \pl \psi_{I,i} \pl k (\psi_PT_P) \psi_{J,i}T_J \pr - \psi_I f \pr dx dy -Q_I = 0 \label{eq:wf7} -\end{equation} -\begin{equation} - \dot{T}_J \int_{\Omega_e}\rho C_p \psi_I \psi_J dx dy + \int_{\Omega_e} \pl \psi_{I,i} \pl k (\psi_PT_P) \psi_{J,i}T_J \pr - \psi_I f \pr dx dy -Q_I =0 \label{eq:wf8} -\end{equation} -We write this in matrix form as -\begin{equation} - M_{IJ}\dot{T}_J + L_I - Q_I=0 \label{eq:matform} -\end{equation} -\begin{equation} - \underbrace{\pl \int_{\Omega_e}\rho C_p \psi_I \psi_J dx dy \pr}_{M_{IJ}} \dot{T}_J + \underbrace{\int_{\Omega_e} \pl \psi_{I,i} \pl k (\psi_PT_P) \psi_{J,i}T_J \pr - \psi_I f \pr dx dy}_{L_I} -Q_I = 0 \label{eq:matform2} -\end{equation} -Since, this a time dependent problem, let us write the \eref{eq:matform} at times $s$ and $s+1$ as follows. This will allow us write the time marching schemes for the problem. -\begin{eqnarray} - M_{IJ}^s\dot{T}_J^{s} + L_I^{s} - Q_I^{s}=0 \label{eq:matofrm2} \\ - M_{IJ}^{s+1}\dot{T}_J^{s+1} + L_I^{s+1} - Q_I^{s+1}=0 \label{eq:matofrm3} -\end{eqnarray} -We will consider the case where $M_{IJ}$ and $Q_I$ are both independent of time. Hence, we have -\begin{eqnarray} - M_{IJ}\dot{T}_J^{s} + L_I^{s} - Q_I=0 \label{eq:matofrm4} \\ - M_{IJ}\dot{T}_J^{s+1} + L_I^{s+1} - Q_I=0 \label{eq:matofrm5} -\end{eqnarray} -\section{Time marching} -Consider the following approximation for the time derivative -\begin{equation} - \Delta t \pl \pl1-\alpha \pr\dot{T}_J^{s} + \alpha \dot{T}_J^{s+1} \pr = T_J^{s+1} - T_J^{s}\label{scheme} -\end{equation} -where $\Delta t$ is the time step and $\alpha$ is $0<\alpha<1$. Multiplying \eref{scheme} with $M_{IJ}$ we get, -\begin{equation} - \Delta t \pl1-\alpha \pr M_{IJ}\dot{T}_J^{s} + \Delta t \alpha M_{IJ}\dot{T}_J^{s+1} = M_{IJ}T_J^{s+1} - M_{IJ}T_J^{s}\label{scheme1} -\end{equation} -combining terms properly, -\begin{equation} - \Delta t \alpha M_{IJ}\dot{T}_J^{s+1} - M_{IJ}T_J^{s+1} = - M_{IJ}T_J^{s} - \Delta t \pl1-\alpha \pr M_{IJ}\dot{T}_J^{s} \label{scheme2} -\end{equation} -We now use \eref{eq:matofrm4} and \eref{eq:matofrm5} to get, -\begin{equation} - \Delta t \alpha \pl -L_I^{s+1}+Q_I \pr - M_{IJ}T_J^{s+1} = - M_{IJ}T_J^{s} - \Delta t \pl1-\alpha \pr \pl -L_I^s+Q_I\pr \label{scheme3} -\end{equation} -expanding -\begin{equation} - - \alpha \Delta t L_I^{s+1} -M_{IJ}T_{J}^{s+1} = -M_{IJ}T_{J}^s + \Delta t \pl 1-\alpha \pr L_I^s-\Delta t \alpha Q_I -\Delta t \pl1-\alpha \pr Q_I -\end{equation} -which gives -\begin{equation} - \alpha \Delta t L_I^{s+1} + M_{IJ}T_{J}^{s+1} = M_{IJ}T_{J}^s - \Delta t \pl 1-\alpha \pr L_I^s+\Delta t \alpha Q_I +\Delta t \pl1-\alpha \pr Q_I -\end{equation} -giving -\begin{equation} - M_{IJ}T_{J}^{s+1} + \alpha \Delta t L_I^{s+1} = M_{IJ}T_{J}^s - \Delta t \pl 1-\alpha \pr L_I^s+ \Delta t Q_I -\end{equation} -Now we have -\begin{equation} - \boxed{ - M_{IJ}T_{J}^{s+1} + \alpha \Delta t L_I^{s+1} - M_{IJ}T_{J}^s + \Delta t \pl 1-\alpha \pr L_I^s - \Delta t Q_I - } - \label{eq:finaleq0} -\end{equation} -We want to find $T_K$s so that the above \eref{eq:finaleq0} is satisfied for every time step $s+1$, in essence we want to find out $T_K^{\ts}$. Consequently, for values of $T_K^{\ts}$ which do not satisfy \eref{eq:finaleq0}, we will have a residual for each time step $s+1$. Hence, we have -\begin{equation} - \boxed{ - M_{IJ}T_{J}^{s+1} + \alpha \Delta t L_I^{s+1} - M_{IJ}T_{J}^s + \Delta t \pl 1-\alpha \pr L_I^s - \Delta t Q_I = R_I^{s+1} -} - \label{eq:finaleq} -\end{equation} -is the residual at the $s+1$$^{th}$ time step. After assembling, we need to find the $T_{J}^{s+1}$ that essentially makes the assembled $R_I^{s+1}$ go to zero. Numerically, we will require some norm of the residual to go to zero. Note that the residual we developed is local to the element. This will have to be assembled in the usual manner before solving. \textit{We again note that $T_K^{s}$ are all known from the previous time step $s$. } - -\section{The Jacobian} -The solution (i.e. to find $T_K^{\ts}$) which make the residual to go to zero, can be found using the Newton Raphson's technique. For this we need to calculate the \textbf{Jacobian} or the \textbf{Tangent stiffness} matrix which involves calculating - $\pd{R_I^{s+1}}{T_{Q}^{s+1}}$. From the final equation \eref{eq:finaleq}, we get, -\begin{eqnarray} -J_{IQ}^{\ts}=\pd{R_I^{s+1}}{T_{Q}^{s+1}} = M_{IJ}\delta_{JQ} + \alpha \Delta t \pd{L_I^{s+1}}{T_Q}\\ -J_{IQ}^{\ts}=\pd{R_I^{s+1}}{T_{Q}^{s+1}} = M_{IQ} + \alpha \Delta t \pd{L_I^{s+1}}{T_Q} \label{jacob} -\end{eqnarray} -This differentiation, is very simple for this problem. But it is not difficult to imagine situations, where this is terribly complicated and algebraically involved. Again, the terms which have a superfix $s$ depend on values of $T_H$ known from the previous step $s$ and hence are known values and not variables; we do not need to differentiate the residual with respect to these known quantities. We know that -\begin{equation} - \int_{\Omega_e} \pl \psi_{I,i} \pl k (\psi_PT_P) \psi_{J,i}T_J \pr - \psi_I f \pr dx dy = L_I^{s+1} -\end{equation} -So, -\begin{equation} - \pd{L_I^{s+1}}{T_Q} = \int_{\Omega_e} \pl \psi_{I,i}\psi_{J,i} \pd{k}{T_Q} T_J + \psi_{I,i}\psi_{J,i} k \delta_{JQ} \pr dx dy -\end{equation} -if we now assume $k = a+bT+cT^2$, with $a,b$ and $c$ being constants, we have, -\begin{equation} - k = a + b (\psi_R T_R) + c \pl \psi_Y T_Y \pr^2 -\end{equation} -then -\begin{eqnarray} - \pd{k}{T_Q} = b\psi_R\delta_{RQ} + 2c (\psi_Y T_Y) \psi_Y\delta_{YQ} \\ - \pd{k}{T_Q} = b\psi_Q + 2c (\psi_Y T_Y) \psi_Q -\end{eqnarray} -and hence we get -\begin{equation} - \pd{L_I^{s+1}}{T_Q} = \int_{\Omega_e} \pl \psi_{I,i} (\psi_{J,i}T_J) (b\psi_Q + 2c (\psi_Y T_Y) \psi_Q) + \psi_{I,i}\psi_{Q,i}(a + b (\psi_R T_R) + c \pl \psi_Y T_Y)^2 \pr \pr dx dy -\end{equation} -Then, $J_{IQ}^{\ts}$ can be written using \eref{jacob}. It is this calculation of the Jacobian which is often complicated and we intend to automate it using the \textit{Automatic differentiation} facility. These are implemented in the functions \te{compute_residual()} and \te{copmute_jacobian()}. These functions are described in detail next. -\section{The functions \lstinline!compute_residual()! and \lstinline!compute_jacobian()!} -The functions \te{compute_residual()} and \te{compute_jacobian()} compute the residuals and the jacobian, respectively. These are done with automatic differentiation from TRILINOS. For dealii tutorials on these, please see step-71 and step-72. The residual is computed in such a way that the a single line will actually differentiate it and calculate the jacobian given by \eref{jacob}. While the difference between \te{compute_residual()} and \te{compute_jacobian()} is minimal and could have been incorporated in a single function, we separate them in this implementation, because we want to use the nonlinear solver provided by TRILINOS (NOX) to actually perform the Newton Raphson iteration. For further details on this (Nonlinear solvers) aspect, please see step-77 of the dealii tutorial. We describe some aspects of the \te{compute_residual()} and \te{copmute_jacobian()} here. - -\subsection{The \te{compute_residual()} function} -The \te{copmute_residual()} function actually takes in two arguments. \te{evaluation_point} and \te{residual}, both by reference. We will discuss this in some detail as we elaborate the function. - -The following lines -%---------------------------------------------------------------------------------- -\lstinputlisting[language=C++, linerange={1-36}]{../source/compute_residual.cc} -should be fairly straightforward. The variable \te{cell_rhs} holds the local (for an element or a cell) residual. Note that we use the \te{FEFaceValues<2>} as we will need to apply the Neuman boundary conditions on the right hand side of the domain. - -The lines below define the \textbf{Number type} we want to use to define our variables so that they are suitable for automatic differentiation. Here, we are using variables that will allow \textit{automatic differentiation}\footnote{The other option is symbolic differentiation using SYMENGINE package. Here we are using the automatic differentiation using SACADO.}. -%---------------------------------------------------------------------------------- -\lstinputlisting[language=C++, linerange={42-43}]{../source/compute_residual.cc} - -The following lines are once again straightforward to understand. -\lstinputlisting[language=C++, linerange={50-58}]{../source/compute_residual.cc} -We only comment that the variable \te{consol} holds the values of the variable \te{converged_solution} at the gauss points of the current, active cell. Similarly, the variable \te{consol_grad} holds the values of the gradients of the \te{converged_solution} at the gauss points. The variable \te{converged_solution} refers is solution at the previous time step (time step $s$). - - -In following lines the variables \te{old_solution} and \te{old_solution_grad} are -crucial. They define the variables of the appropriate number type using which we will define our residual, so that it is suitable for the automatic differentiation. \te{old_solution} is a solution variable of the appropriate number type with respect to which we need to define the residual and then differentiate it to get the Jacobian. For ease of understanding, one may consider this variable to be a special variable, that will look like this -\begin{equation} - \text{\te{old_solution}} = - \begin{bmatrix} - \psi_i (\zeta_1,\eta_1) d_i\\ - \psi_i (\zeta_2,\eta_2) d_i\\ - \psi_i (\zeta_3,\eta_3) d_i\\ - \psi_i (\zeta_4,\eta_4) d_i - \label{eq:sum} - \end{bmatrix} -\end{equation} -where $\psi_i$ are the shape functions and $d_i$ is some way of representing the primary variable in the problem (in this case the temperature). Could be viewed as something like a symbol, which could be used to define functions and then differentiated later. \te{old_solution} holds the values in the corresponding gauss points and the repeated indices in \eref{eq:sum} indicate summation over the nodes. Similarly, the variable \te{old_solution_grad} holds the corresponding gradients. -\lstinputlisting[language=C++, linerange={59-68}]{../source/compute_residual.cc} - -The lines -\lstinputlisting[language=C++, linerange={69-73}]{../source/compute_residual.cc} -are straightforward to understand. The vector \te{residual} will hold the \textbf{numerical} value of the residual and is hence initialized to be zero. - -The lines -\lstinputlisting[language=C++, linerange={74-80}]{../source/compute_residual.cc} -should be clear. In the following line we inform how many independent and dependent variables we will have. -\lstinputlisting[language=C++, linerange={81}]{../source/compute_residual.cc} -then, in the following line we tell the system, at which point it should numerically evaluate the residual. This is the point given by \te{evaluation_point}, which will contain the solution that is getting iterated to find the actual solution for each time step. -\lstinputlisting[language=C++, linerange={86}]{../source/compute_residual.cc} - -Now the lines -\lstinputlisting[language=C++, linerange={88-92}]{../source/compute_residual.cc} -the actual storage of the values of the shape functions and its derivatives at the gauss points is taking place in the variables \te{old_solution} and \te{old_solution_grad}. - -In the lines -\lstinputlisting[language=C++, linerange={97-98}]{../source/compute_residual.cc} -the values of the \te{converged_solution} and its gradients get stored in the variables \te{consol} and \te{con_sol_grad}. - -Then in the lines -\lstinputlisting[language=C++, linerange={104-120}]{../source/compute_residual.cc} -the actual evaluation of the residual is taking place. Note that each intermediate variable used is of the \te{ADNumberType}. Further, the residual is stored in a variable called \te{residual_ad[i]} because it is also of the special type that will be suitable for differentiation. Note that the symbols in the residual equation in this document and the variables used in the code, map in the following manner. The residual is pasted below again -\begin{equation} - \boxed{ - M_{IJ}T_{J}^{s+1} + \alpha \Delta t L_I^{s+1} - M_{IJ}T_{J}^s + \Delta t \pl 1-\alpha \pr L_I^s - \Delta t Q_I = R_I^{s+1} - } - \label{eq:finaleq2} -\end{equation} -\begin{itemize} - \item $M_{IJ}T_{J}^{s+1}$ $\rightarrow$ \te{MijTjcurr} - \item $M_{IJ}T_{J}^{s}$ $\rightarrow$ \te{MijTjprev} - \item $L_I^{s+1}$ $\rightarrow$ \te{Licurr} - \item $L_I^{s}$ $\rightarrow$ \te{Liprev} - \item $R_I^{s+1}$ $\rightarrow$ \te{residual_ad[i]} -\end{itemize} -then the following lines implement the Neumann boundary conditions -\lstinputlisting[language=C++, linerange={124-140}]{../source/compute_residual.cc} - -Now -\lstinputlisting[language=C++, linerange={144}]{../source/compute_residual.cc} -actually tells the \te{ad_helper} that this is the residual it should use, and the line -\lstinputlisting[language=C++, linerange={150}]{../source/compute_residual.cc} -evaluates the residual at the appropriate point, which is the \te{evaluation_point}. Then -\lstinputlisting[language=C++, linerange={151-156}]{../source/compute_residual.cc} are used to calculate the global \te{residual} from the local \te{cell_rhs}. Now \te{residual} is numerical and not of the special number type. Finally, lines -\lstinputlisting[language=C++, linerange={158-162}]{../source/compute_residual.cc} -are routine, in particular, it makes the global degrees of freedom, where the dirichlet boundary conditions are applied to be zero. -\subsection{The \te{compute_jacobian()} function} -The \te{compute_jacobian()} only takes in the \te{evaluation_point} as an input and is identical to the -\te{compute_residual()} except for some minor modification. In particular the line -\lstinputlisting[language=C++, linerange={109}]{../source/compute_jacobian.cc} -computes the actual jacobian by performing the automatic differentiation and evaluates it at the \te{evaluation_point}. The remaining lines are routine and should be straightforward. - -\section{General ideas involved in solving coupled nonlinear equations using Newton Raphson's technique} -Consider that we have the equations -\begin{eqnarray} - f_i^1(p_j^{\ts}, T_j^{\ts})=0 \label{eq13}\\ - f_i^2(p_j^{\ts}, T_j^{\ts})=0 \label{eq14} - \label{eqn} -\end{eqnarray} -are coupled non-linear algebraic equations in the variables $p_j$ and $T_j$. In the problem we are trying to solve, we only have one unknown, but this set of derivation clarifies what one should do for coupled problems as well. - -We need to solve the set \eref{eqn} using Newton-Raphson's technique. Now suppose we have 4 noded linear element, with values as $p_1^{\ts},p_2^{\ts},p_3^{\ts},p_4^{\ts}$ and $T_1^{\ts},T_2^{\ts},T_3^{\ts},T_4^{\ts}$, then actually $f_i^1$ and $f_i^2$ are -\begin{eqnarray} - f_i^1(p_1^{\ts},p_2^{\ts},p_3^{\ts},p_4^{\ts},T_1^{\ts},T_2^{\ts},T_3^{\ts},T_4^{\ts})=0 \\ - f_i^2(p_1^{\ts},p_2^{\ts},p_3^{\ts},p_4^{\ts},T_1^{\ts},T_2^{\ts},T_3^{\ts},T_4^{\ts})=0 -\end{eqnarray} -Also, note that $i$ goes from 1 to 4 for $f^1$ and from 1 to 4 for $f^2$. For the Newton Raphson's technique, we have to find the Jacobian. This is written explicitly below for clear understanding. -\begin{equation} - J_{pq}=- - \begin{bmatrix} - \JPP{1}{1} & \JPP{1}{2} & \JPP{1}{3} & \JPP{1}{4} & \JPT{1}{1} & \JPT{1}{2}&\JPT{1}{3}&\JPT{1}{4}\\ - - \JPP{2}{1} & \JPP{2}{2} & \JPP{2}{3} & \JPP{2}{4} & \JPT{2}{1} & \JPT{2}{2}&\JPT{2}{3}&\JPT{2}{4}\\ - - \JPP{3}{1} & \JPP{3}{2} & \JPP{3}{3} & \JPP{3}{4} & \JPT{3}{1} & \JPT{3}{2}&\JPT{3}{3}&\JPT{3}{4}\\ - - \JPP{4}{1} & \JPP{4}{2} & \JPP{4}{3} & \JPP{4}{4} & \JPT{4}{1} & \JPT{4}{2}&\JPT{4}{3}&\JPT{4}{4}\\ - %================================================================ - \JTP{1}{1} & \JTP{1}{2} & \JTP{1}{3} & \JTP{1}{4} & \JTT{1}{1} & \JTT{1}{2}&\JTT{1}{3}&\JTT{1}{4}\\ - - \JTP{2}{1} & \JTP{2}{2} & \JTP{2}{3} & \JTP{2}{4} & \JTT{2}{1} & \JTT{2}{2}&\JTT{2}{3}&\JTT{2}{4}\\ - - \JTP{3}{1} & \JTP{3}{2} & \JTP{3}{3} & \JTP{3}{4} & \JTT{3}{1} & \JTT{3}{2}&\JTT{3}{3}&\JTT{3}{4}\\ - - \JTP{4}{1} & \JTP{4}{2} & \JTP{4}{3} & \JTP{4}{4} & \JTT{4}{1} & \JTT{4}{2}&\JTT{4}{3}&\JTT{4}{4} - \end{bmatrix} -\end{equation} -We use the Jacobian evaluated at the values of the variables ($p_j^{\ts}$, $T_j^{\ts}$) at the \textit{previous iteration} (denoted by the variable \lstinline{present_solution} in the code) to obtain a new value of the required vector. That is -\begin{equation} - \begin{bmatrix} - p_1^{s+1}\\ - p_2^{s+1}\\ - p_3^{s+1}\\ - p_4^{s+1}\\ - T_1^{s+1}\\ - T_2^{s+1}\\ - T_3^{s+1}\\ - T_4^{s+1} - \end{bmatrix}^{k+1} - = - \begin{bmatrix} - p_1^{s+1}\\ - p_2^{s+1}\\ - p_3^{s+1}\\ - p_4^{s+1}\\ - T_1^{s+1}\\ - T_2^{s+1}\\ - T_3^{s+1}\\ - T_4^{s+1} - \end{bmatrix}^{k}- - \begin{bmatrix} - \JPP{1}{1} & \JPP{1}{2} & \JPP{1}{3} & \JPP{1}{4} & \JPT{1}{1} & \JPT{1}{2}&\JPT{1}{3}&\JPT{1}{4}\\ - - \JPP{2}{1} & \JPP{2}{2} & \JPP{2}{3} & \JPP{2}{4} & \JPT{2}{1} & \JPT{2}{2}&\JPT{2}{3}&\JPT{2}{4}\\ - - \JPP{3}{1} & \JPP{3}{2} & \JPP{3}{3} & \JPP{3}{4} & \JPT{3}{1} & \JPT{3}{2}&\JPT{3}{3}&\JPT{3}{4}\\ - - \JPP{4}{1} & \JPP{4}{2} & \JPP{4}{3} & \JPP{4}{4} & \JPT{4}{1} & \JPT{4}{2}&\JPT{4}{3}&\JPT{4}{4}\\ - %================================================================ - \JTP{1}{1} & \JTP{1}{2} & \JTP{1}{3} & \JTP{1}{4} & \JTT{1}{1} & \JTT{1}{2}&\JTT{1}{3}&\JTT{1}{4}\\ - - \JTP{2}{1} & \JTP{2}{2} & \JTP{2}{3} & \JTP{2}{4} & \JTT{2}{1} & \JTT{2}{2}&\JTT{2}{3}&\JTT{2}{4}\\ - - \JTP{3}{1} & \JTP{3}{2} & \JTP{3}{3} & \JTP{3}{4} & \JTT{3}{1} & \JTT{3}{2}&\JTT{3}{3}&\JTT{3}{4}\\ - - \JTP{4}{1} & \JTP{4}{2} & \JTP{4}{3} & \JTP{4}{4} & \JTT{4}{1} & \JTT{4}{2}&\JTT{4}{3}&\JTT{4}{4} - \end{bmatrix}^{-1} - \begin{bmatrix} - f_1^{1}\\ - f_2^{1}\\ - f_3^{1}\\ - f_4^{1}\\ - f_1^{2}\\ - f_2^{2}\\ - f_3^{2}\\ - f_4^{2} - \end{bmatrix}^{k} -\end{equation} -The $J_{pq}$ is evaluated with values obtained at the $k^{th}$ iteration. To avoid taking the inverse, we solve the following, for the changes $\Delta p$ and $\Delta T$. -\begin{equation} - \begin{bmatrix} - \JPP{1}{1} & \JPP{1}{2} & \JPP{1}{3} & \JPP{1}{4} & \JPT{1}{1} & \JPT{1}{2}&\JPT{1}{3}&\JPT{1}{4}\\ - - \JPP{2}{1} & \JPP{2}{2} & \JPP{2}{3} & \JPP{2}{4} & \JPT{2}{1} & \JPT{2}{2}&\JPT{2}{3}&\JPT{2}{4}\\ - - \JPP{3}{1} & \JPP{3}{2} & \JPP{3}{3} & \JPP{3}{4} & \JPT{3}{1} & \JPT{3}{2}&\JPT{3}{3}&\JPT{3}{4}\\ - - \JPP{4}{1} & \JPP{4}{2} & \JPP{4}{3} & \JPP{4}{4} & \JPT{4}{1} & \JPT{4}{2}&\JPT{4}{3}&\JPT{4}{4}\\ - %================================================================ - \JTP{1}{1} & \JTP{1}{2} & \JTP{1}{3} & \JTP{1}{4} & \JTT{1}{1} & \JTT{1}{2}&\JTT{1}{3}&\JTT{1}{4}\\ - - \JTP{2}{1} & \JTP{2}{2} & \JTP{2}{3} & \JTP{2}{4} & \JTT{2}{1} & \JTT{2}{2}&\JTT{2}{3}&\JTT{2}{4}\\ - - \JTP{3}{1} & \JTP{3}{2} & \JTP{3}{3} & \JTP{3}{4} & \JTT{3}{1} & \JTT{3}{2}&\JTT{3}{3}&\JTT{3}{4}\\ - - \JTP{4}{1} & \JTP{4}{2} & \JTP{4}{3} & \JTP{4}{4} & \JTT{4}{1} & \JTT{4}{2}&\JTT{4}{3}&\JTT{4}{4} - \end{bmatrix} - \begin{bmatrix} - \Delta p_1^{s+1}\\ - \Delta p_2^{s+1}\\ - \Delta p_3^{s+1}\\ - \Delta p_4^{s+1}\\ - \Delta T_1^{s+1}\\ - \Delta T_2^{s+1}\\ - \Delta T_3^{s+1}\\ - \Delta T_4^{s+1} - \end{bmatrix}= - \begin{bmatrix} - f_1^{1}\\ - f_2^{1}\\ - f_3^{1}\\ - f_4^{1}\\ - f_1^{2}\\ - f_2^{2}\\ - f_3^{2}\\ - f_4^{2} - \end{bmatrix}^{k} -\end{equation} -This is written as -\begin{equation} - J_{pq}\{\Delta solution\}^{k+1}=\{function\}^{k} -\end{equation} -Then add this to the value of -\begin{equation} - \begin{bmatrix} - p_1^{s+1}\\ - p_2^{s+1}\\ - p_3^{s+1}\\ - p_4^{s+1}\\ - T_1^{s+1}\\ - T_2^{s+1}\\ - T_3^{s+1}\\ - T_4^{s+1} - \end{bmatrix}^{k} -\end{equation} -We keep doing this, until the $L_2$-error of the vector is small. (Some criterion is used). -\begin{equation} - ||(function)^k||<\delta^0 ||function^0|| -\end{equation} -or -\begin{equation} - ||(function)^k||