nonlinear-heat_transfer_with_AD_NOX.h

Go to the documentation of this file.

 
#include "allheaders.h"
#include "nonlinear_heat.h"
void nonlinear_heat::compute_residual(const Vector<double> & evaluation_point, Vector<double> & residual)
{
    const QGauss<2> quadrature_formula(
            fe.degree + 1);
    const QGauss<1> face_quadrature_formula(fe.degree+1); //Define quadrature for integration over faces */
    FEValues<2> fe_values(fe,
                          quadrature_formula,
                          update_values | update_gradients |
                          update_quadrature_points |
                          update_JxW_values); 
    FEFaceValues<2> fe_face_values(fe,face_quadrature_formula,update_values|update_quadrature_points|
                                                              update_normal_vectors | update_JxW_values);
    const unsigned int dofs_per_cell = fe.dofs_per_cell; 
    const unsigned int n_q_points = quadrature_formula.size();
    const unsigned int n_q_f_points = face_quadrature_formula.size();
    Vector<double> cell_rhs(dofs_per_cell); 
@endcode
should be fairly straightforward. The variable `cell_rhs` holds the local (for an element or a cell) residual. Note that we use the `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 *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 /automatic differentiation/. (The other option is symbolic differentiation using SYMENGINE package. Here we are using the automatic differentiation using SACADO.)
 
@code{.sh}
    using ADHelper = Differentiation::AD::ResidualLinearization<Differentiation::AD::NumberTypes::sacado_dfad,double>;
    using ADNumberType = typename ADHelper::ad_type;
 
@endcode
 
The following lines are once again straightforward to understand. 
@code{.sh}
    const FEValuesExtractors::Scalar t(
            0);
    std::vector<types::global_dof_index> local_dof_indices(
            dofs_per_cell); 
    std::vector<double> consol(n_q_points); /* Converged solution at the Gauss points from the previous time step*/
 
@endcode
We only comment that the variable `consol` holds the values of the variable `converged_solution` at the gauss points of the current, active cell. Similarly, the variable `consol_grad` holds the values of the gradients of the `converged_solution` at the gauss points. The variable `converged_solution` refers is solution at the previous time step (time step @f$s@f$). 
 
 
In following lines the variables `old_solution` and `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.  `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
@f[
        \text{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
        \end{bmatrix}
@f] 
where @f$\psi_i@f$ are the shape functions and @f$d_i@f$ 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.  `old_solution` holds the values in the corresponding gauss points and the repeated indices in the previous equation indicate summation over the nodes. Similarly, the variable `old_solution_grad` holds the corresponding gradients. 
@code{.sh}
    std::vector<ADNumberType> old_solution(
            n_q_points); /* Current solution at the Gauss points at this iteration for the current time*/
    std::vector<Tensor<1, 2>> consol_grad(
            n_q_points); /* Converged gradients of the solutions at the Gauss points from the previous time step */
    std::vector<Tensor<1, 2,ADNumberType>> old_solution_grad(
            n_q_points);
 
 
@endcode
 
The lines
@code{.sh}
    ComponentMask t_mask = fe.component_mask(t);
    residual = 0.0;
 
@endcode
are straightforward to understand. The vector `residual` will hold the *numerical* value of the residual and is hence initialized to be zero. 
 
The lines
@code{.sh}
    for (const auto &cell: dof_handler.active_cell_iterators())
    {
        cell_rhs = 0;
        fe_values.reinit(cell);
        cell->get_dof_indices(local_dof_indices);
        const unsigned int n_independent_variables = local_dof_indices.size();
        const unsigned int n_dependent_variables = dofs_per_cell;
 
@endcode
should be clear. In the following line we inform how many independent and dependent variables we will have. 
@code{.sh}
        ADHelper ad_helper(n_independent_variables, n_dependent_variables);
 
@endcode
then, in the following line we tell the system, at which point it should numerically evaluate the residual. This is the point given by `evaluation_point`, which will contain the solution that is getting iterated to find the actual solution for each time step.
@code{.sh}
        ad_helper.register_dof_values(evaluation_point,local_dof_indices);
 
@endcode
 
Now the lines
@code{.sh}
        const std::vector<ADNumberType> &dof_values_ad = ad_helper.get_sensitive_dof_values();
        fe_values[t].get_function_values_from_local_dof_values(dof_values_ad,
                                         old_solution);
        fe_values[t].get_function_gradients_from_local_dof_values(dof_values_ad,
                                            old_solution_grad);
 
@endcode
the actual storage of the values of the shape functions and its derivatives at the gauss points is taking place in the variables `old_solution` and `old_solution_grad`. 
 
In the lines 
@code{.sh}
        fe_values[t].get_function_values(converged_solution, consol);
        fe_values[t].get_function_gradients(converged_solution, consol_grad);
 
@endcode
the values of the `converged_solution` and its gradients get stored in the variables `consol` and `con_sol_grad`.
 
Then in the lines
@code{.sh}
        for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
        {
            for (unsigned int i = 0; i < dofs_per_cell; ++i)
            {
                ADNumberType MijTjcurr = Cp*rho*fe_values[t].value(i,q_index)*old_solution[q_index];
                ADNumberType MijTjprev = Cp*rho*fe_values[t].value(i,q_index)*consol[q_index];
                ADNumberType k_curr = a + b*old_solution[q_index] + c*std::pow(old_solution[q_index],2);
                ADNumberType k_prev = a + b*consol[q_index] + c*std::pow(consol[q_index],2);
                ADNumberType Licurr =  alpha * delta_t *  (fe_values[t].gradient(i,q_index)*k_curr*old_solution_grad[q_index]);
                ADNumberType Liprev =  (1-alpha) * delta_t *  (fe_values[t].gradient(i,q_index)*k_prev*consol_grad[q_index]);
                residual_ad[i] +=  (MijTjcurr+Licurr-MijTjprev+Liprev)*fe_values.JxW(q_index);
            }
        }
                        residual_ad[i]+= -delta_t*(-10)*fe_face_values[t].value(i,q_point)*fe_face_values.JxW(q_point);
 
                }
            }
        }
 
@endcode
 
Now 
@code{.sh}
        ad_helper.register_residual_vector(residual_ad);
 
@endcode
actually tells the `ad_helper` that this is the residual it should use, and the line
@code{.sh}
        ad_helper.compute_residual(cell_rhs);
 
@endcode
evaluates the residual at the appropriate point, which is the `evaluation_point`. Then
@code{.sh}
        cell->get_dof_indices(local_dof_indices);
 
        for (unsigned int i =0;i < dofs_per_cell; ++i)
            residual(local_dof_indices[i])+= cell_rhs(i);
 
    }
 
@endcode
are used to calculate the global `residual` from the local `cell_rhs`. Now `residual` is numerical and not of the special number type. Finally, lines
@code{.sh}
    for(const types::global_dof_index i: DoFTools::extract_boundary_dofs(dof_handler,t_mask,{1}))
        residual(i) = 0;
 
    std::cout << " The Norm is :: = " << residual.l2_norm() << std::endl;
}
 
@endcode 
are routine, in particular, it makes the global degrees of freedom, where the dirichlet boundary conditions are applied to be zero.
 
 
#### The compute_jacobian() function
 
The `compute_jacobian()` only takes in the `evaluation_point` as an input and is identical to the 
`compute_residual()` except for some minor modification. In particular the line
@code{.sh}
        ad_helper.compute_linearization(cell_matrix);
 
@endcode 
computes the actual jacobian by performing the automatic differentiation and evaluates it at the `evaluation_point`. The remaining lines are routine and should be straightforward. 
 
 
### General ideas involved in solving coupled nonlinear equations using Newton Raphson's technique
Consider that we have the equations 
@f{eqnarray}{
        f_i^1(p_j^{{s+1}}, T_j^{{s+1}})=0 \\
        f_i^2(p_j^{{s+1}}, T_j^{{s+1}})=0 
@f}
are coupled non-linear algebraic equations in the variables @f$p_j@f$ and @f$T_j@f$. 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 equations immediately above using Newton-Raphson's technique. Now suppose we have 4 noded linear element, with values as @f$p_1^{{s+1}},p_2^{{s+1}},p_3^{{s+1}},p_4^{{s+1}}@f$ and @f$T_1^{{s+1}},T_2^{{s+1}},T_3^{{s+1}},T_4^{{s+1}}@f$, then actually @f$f_i^1@f$ and @f$f_i^2@f$ are 
@f{eqnarray}{
        f_i^1(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}})=0 \\
        f_i^2(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}})=0
@f}
Also, note that @f$i@f$ goes from 1 to 4 for @f$f^1@f$ and from 1 to 4 for @f$f^2@f$. For the Newton Raphson's technique, we have to find the Jacobian. This is written explicitly below for clear understanding. 
@f[
        J_{pq}=-
        \begin{bmatrix}
                \frac{\partial f^1_1}{\partial p^{s+1}_1} & \frac{\partial f^1_1}{\partial p^{s+1}_2} & \frac{\partial f^1_1}{\partial p^{s+1}_3} & \frac{\partial f^1_1}{\partial p^{s+1}_4} & \frac{\partial f^1_1}{\partial T^{s+1}_1} & \frac{\partial f^1_1}{\partial T^{s+1}_2}&\frac{\partial f^1_1}{\partial T^{s+1}_3}&\frac{\partial f^1_1}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^1_2}{\partial p^{s+1}_1} & \frac{\partial f^1_2}{\partial p^{s+1}_2} & \frac{\partial f^1_2}{\partial p^{s+1}_3} & \frac{\partial f^1_2}{\partial p^{s+1}_4} & \frac{\partial f^1_2}{\partial T^{s+1}_1} & \frac{\partial f^1_2}{\partial T^{s+1}_2}&\frac{\partial f^1_2}{\partial T^{s+1}_3}&\frac{\partial f^1_2}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^1_3}{\partial p^{s+1}_1} & \frac{\partial f^1_3}{\partial p^{s+1}_2} & \frac{\partial f^1_3}{\partial p^{s+1}_3} & \frac{\partial f^1_3}{\partial p^{s+1}_4} & \frac{\partial f^1_3}{\partial T^{s+1}_1} & \frac{\partial f^1_3}{\partial T^{s+1}_2}&\frac{\partial f^1_3}{\partial T^{s+1}_3}&\frac{\partial f^1_3}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^1_4}{\partial p^{s+1}_1} & \frac{\partial f^1_4}{\partial p^{s+1}_2} & \frac{\partial f^1_4}{\partial p^{s+1}_3} & \frac{\partial f^1_4}{\partial p^{s+1}_4} & \frac{\partial f^1_4}{\partial T^{s+1}_1} & \frac{\partial f^1_4}{\partial T^{s+1}_2}&\frac{\partial f^1_4}{\partial T^{s+1}_3}&\frac{\partial f^1_4}{\partial T^{s+1}_4}\\
                %================================================================
                \frac{\partial f^2_1}{\partial p^{s+1}_1} & \frac{\partial f^2_1}{\partial p^{s+1}_2} & \frac{\partial f^2_1}{\partial p^{s+1}_3} & \frac{\partial f^2_1}{\partial p^{s+1}_4} & \frac{\partial f^2_1}{\partial T^{s+1}_1} & \frac{\partial f^2_1}{\partial T^{s+1}_2}&\frac{\partial f^2_1}{\partial T^{s+1}_3}&\frac{\partial f^2_1}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^2_2}{\partial p^{s+1}_1} & \frac{\partial f^2_2}{\partial p^{s+1}_2} & \frac{\partial f^2_2}{\partial p^{s+1}_3} & \frac{\partial f^2_2}{\partial p^{s+1}_4} & \frac{\partial f^2_2}{\partial T^{s+1}_1} & \frac{\partial f^2_2}{\partial T^{s+1}_2}&\frac{\partial f^2_2}{\partial T^{s+1}_3}&\frac{\partial f^2_2}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^2_3}{\partial p^{s+1}_1} & \frac{\partial f^2_3}{\partial p^{s+1}_2} & \frac{\partial f^2_3}{\partial p^{s+1}_3} & \frac{\partial f^2_3}{\partial p^{s+1}_4} & \frac{\partial f^2_3}{\partial T^{s+1}_1} & \frac{\partial f^2_3}{\partial T^{s+1}_2}&\frac{\partial f^2_3}{\partial T^{s+1}_3}&\frac{\partial f^2_3}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^2_4}{\partial p^{s+1}_1} & \frac{\partial f^2_4}{\partial p^{s+1}_2} & \frac{\partial f^2_4}{\partial p^{s+1}_3} & \frac{\partial f^2_4}{\partial p^{s+1}_4} & \frac{\partial f^2_4}{\partial T^{s+1}_1} & \frac{\partial f^2_4}{\partial T^{s+1}_2}&\frac{\partial f^2_4}{\partial T^{s+1}_3}&\frac{\partial f^2_4}{\partial T^{s+1}_4}
        \end{bmatrix}
@f]
We use the Jacobian evaluated at the values of the variables (@f$p_j^{{s+1}}@f$, @f$T_j^{{s+1}}@f$) at the /previous iteration/ (denoted by the variable `present_solution` in the code) to obtain a new value of the required vector. That is
@f[
        \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}
                \frac{\partial f^1_1}{\partial p^{s+1}_1} & \frac{\partial f^1_1}{\partial p^{s+1}_2} & \frac{\partial f^1_1}{\partial p^{s+1}_3} & \frac{\partial f^1_1}{\partial p^{s+1}_4} & \frac{\partial f^1_1}{\partial T^{s+1}_1} & \frac{\partial f^1_1}{\partial T^{s+1}_2}&\frac{\partial f^1_1}{\partial T^{s+1}_3}&\frac{\partial f^1_1}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^1_2}{\partial p^{s+1}_1} & \frac{\partial f^1_2}{\partial p^{s+1}_2} & \frac{\partial f^1_2}{\partial p^{s+1}_3} & \frac{\partial f^1_2}{\partial p^{s+1}_4} & \frac{\partial f^1_2}{\partial T^{s+1}_1} & \frac{\partial f^1_2}{\partial T^{s+1}_2}&\frac{\partial f^1_2}{\partial T^{s+1}_3}&\frac{\partial f^1_2}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^1_3}{\partial p^{s+1}_1} & \frac{\partial f^1_3}{\partial p^{s+1}_2} & \frac{\partial f^1_3}{\partial p^{s+1}_3} & \frac{\partial f^1_3}{\partial p^{s+1}_4} & \frac{\partial f^1_3}{\partial T^{s+1}_1} & \frac{\partial f^1_3}{\partial T^{s+1}_2}&\frac{\partial f^1_3}{\partial T^{s+1}_3}&\frac{\partial f^1_3}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^1_4}{\partial p^{s+1}_1} & \frac{\partial f^1_4}{\partial p^{s+1}_2} & \frac{\partial f^1_4}{\partial p^{s+1}_3} & \frac{\partial f^1_4}{\partial p^{s+1}_4} & \frac{\partial f^1_4}{\partial T^{s+1}_1} & \frac{\partial f^1_4}{\partial T^{s+1}_2}&\frac{\partial f^1_4}{\partial T^{s+1}_3}&\frac{\partial f^1_4}{\partial T^{s+1}_4}\\
                %================================================================
                \frac{\partial f^2_1}{\partial p^{s+1}_1} & \frac{\partial f^2_1}{\partial p^{s+1}_2} & \frac{\partial f^2_1}{\partial p^{s+1}_3} & \frac{\partial f^2_1}{\partial p^{s+1}_4} & \frac{\partial f^2_1}{\partial T^{s+1}_1} & \frac{\partial f^2_1}{\partial T^{s+1}_2}&\frac{\partial f^2_1}{\partial T^{s+1}_3}&\frac{\partial f^2_1}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^2_2}{\partial p^{s+1}_1} & \frac{\partial f^2_2}{\partial p^{s+1}_2} & \frac{\partial f^2_2}{\partial p^{s+1}_3} & \frac{\partial f^2_2}{\partial p^{s+1}_4} & \frac{\partial f^2_2}{\partial T^{s+1}_1} & \frac{\partial f^2_2}{\partial T^{s+1}_2}&\frac{\partial f^2_2}{\partial T^{s+1}_3}&\frac{\partial f^2_2}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^2_3}{\partial p^{s+1}_1} & \frac{\partial f^2_3}{\partial p^{s+1}_2} & \frac{\partial f^2_3}{\partial p^{s+1}_3} & \frac{\partial f^2_3}{\partial p^{s+1}_4} & \frac{\partial f^2_3}{\partial T^{s+1}_1} & \frac{\partial f^2_3}{\partial T^{s+1}_2}&\frac{\partial f^2_3}{\partial T^{s+1}_3}&\frac{\partial f^2_3}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^2_4}{\partial p^{s+1}_1} & \frac{\partial f^2_4}{\partial p^{s+1}_2} & \frac{\partial f^2_4}{\partial p^{s+1}_3} & \frac{\partial f^2_4}{\partial p^{s+1}_4} & \frac{\partial f^2_4}{\partial T^{s+1}_1} & \frac{\partial f^2_4}{\partial T^{s+1}_2}&\frac{\partial f^2_4}{\partial T^{s+1}_3}&\frac{\partial f^2_4}{\partial T^{s+1}_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}
@f]
The @f$J_{pq}@f$ is evaluated with values obtained at the @f$k^{th}@f$ iteration. To avoid taking the inverse, we solve the following, for the changes @f$\Delta p@f$ and @f$\Delta T@f$.
@f[
        \begin{bmatrix}
                \frac{\partial f^1_1}{\partial p^{s+1}_1} & \frac{\partial f^1_1}{\partial p^{s+1}_2} & \frac{\partial f^1_1}{\partial p^{s+1}_3} & \frac{\partial f^1_1}{\partial p^{s+1}_4} & \frac{\partial f^1_1}{\partial T^{s+1}_1} & \frac{\partial f^1_1}{\partial T^{s+1}_2}&\frac{\partial f^1_1}{\partial T^{s+1}_3}&\frac{\partial f^1_1}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^1_2}{\partial p^{s+1}_1} & \frac{\partial f^1_2}{\partial p^{s+1}_2} & \frac{\partial f^1_2}{\partial p^{s+1}_3} & \frac{\partial f^1_2}{\partial p^{s+1}_4} & \frac{\partial f^1_2}{\partial T^{s+1}_1} & \frac{\partial f^1_2}{\partial T^{s+1}_2}&\frac{\partial f^1_2}{\partial T^{s+1}_3}&\frac{\partial f^1_2}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^1_3}{\partial p^{s+1}_1} & \frac{\partial f^1_3}{\partial p^{s+1}_2} & \frac{\partial f^1_3}{\partial p^{s+1}_3} & \frac{\partial f^1_3}{\partial p^{s+1}_4} & \frac{\partial f^1_3}{\partial T^{s+1}_1} & \frac{\partial f^1_3}{\partial T^{s+1}_2}&\frac{\partial f^1_3}{\partial T^{s+1}_3}&\frac{\partial f^1_3}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^1_4}{\partial p^{s+1}_1} & \frac{\partial f^1_4}{\partial p^{s+1}_2} & \frac{\partial f^1_4}{\partial p^{s+1}_3} & \frac{\partial f^1_4}{\partial p^{s+1}_4} & \frac{\partial f^1_4}{\partial T^{s+1}_1} & \frac{\partial f^1_4}{\partial T^{s+1}_2}&\frac{\partial f^1_4}{\partial T^{s+1}_3}&\frac{\partial f^1_4}{\partial T^{s+1}_4}\\
                %================================================================
                \frac{\partial f^2_1}{\partial p^{s+1}_1} & \frac{\partial f^2_1}{\partial p^{s+1}_2} & \frac{\partial f^2_1}{\partial p^{s+1}_3} & \frac{\partial f^2_1}{\partial p^{s+1}_4} & \frac{\partial f^2_1}{\partial T^{s+1}_1} & \frac{\partial f^2_1}{\partial T^{s+1}_2}&\frac{\partial f^2_1}{\partial T^{s+1}_3}&\frac{\partial f^2_1}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^2_2}{\partial p^{s+1}_1} & \frac{\partial f^2_2}{\partial p^{s+1}_2} & \frac{\partial f^2_2}{\partial p^{s+1}_3} & \frac{\partial f^2_2}{\partial p^{s+1}_4} & \frac{\partial f^2_2}{\partial T^{s+1}_1} & \frac{\partial f^2_2}{\partial T^{s+1}_2}&\frac{\partial f^2_2}{\partial T^{s+1}_3}&\frac{\partial f^2_2}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^2_3}{\partial p^{s+1}_1} & \frac{\partial f^2_3}{\partial p^{s+1}_2} & \frac{\partial f^2_3}{\partial p^{s+1}_3} & \frac{\partial f^2_3}{\partial p^{s+1}_4} & \frac{\partial f^2_3}{\partial T^{s+1}_1} & \frac{\partial f^2_3}{\partial T^{s+1}_2}&\frac{\partial f^2_3}{\partial T^{s+1}_3}&\frac{\partial f^2_3}{\partial T^{s+1}_4}\\
                
                \frac{\partial f^2_4}{\partial p^{s+1}_1} & \frac{\partial f^2_4}{\partial p^{s+1}_2} & \frac{\partial f^2_4}{\partial p^{s+1}_3} & \frac{\partial f^2_4}{\partial p^{s+1}_4} & \frac{\partial f^2_4}{\partial T^{s+1}_1} & \frac{\partial f^2_4}{\partial T^{s+1}_2}&\frac{\partial f^2_4}{\partial T^{s+1}_3}&\frac{\partial f^2_4}{\partial T^{s+1}_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}
@f]
This is written as
@f[
        J_{pq}\{\Delta solution\}^{k+1}=\{function\}^{k}
@f]
Then add this to the value of 
@f[
        \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}
@f]
We keep doing this, until the @f$L_2@f$-error of the vector is small. Determining whether the norm is small can be done using either
@f[
        ||(function)^k||<\delta^0 ||function^0||
@f]
or
@f[
        ||(function)^k||<tol
@f]
where @f$tol@f$ is some small number. In the code, we use the second of these two choices. The final converged solution in the one that satisfies the set of equations and is called by the variable `converged_solution` in the code. This iterative process is carried out for every time step. The variable `present_solution` is given the value of `converged_solution` at the beginning of the iteration of a given time step, so that it serves as a better initial guess for the nonlinear set of equations for the next time step.
 
 
#### The run() function
 
The function which implements the nonlinear solver is the `run()` function for every time step. In particular the lines
@code{.sh}
            typename TrilinosWrappers::NOXSolver<Vector<double>>::AdditionalData additional_data;
            additional_data.abs_tol = target_tolerance;
            additional_data.max_iter = 100;
            TrilinosWrappers::NOXSolver<Vector<double>> nonlinear_solver(
                    additional_data);
 
            nonlinear_solver.residual =
                    [&](const Vector<double> &evaluation_point,
                        Vector<double> &residual) {
                        compute_residual(evaluation_point, residual);
                    };
 
            nonlinear_solver.setup_jacobian =
                    [&](const Vector<double> &current_u) {
                        compute_jacobian(current_u);
                    };
            nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
                                                       Vector<double> &dst,
                                                       const double tolerance) {
                solve(rhs, dst, tolerance);
            };
            nonlinear_solver.solve(present_solution);
        }
 
@endcode
Here, one needs to give how the residual is to be calculated and at what point, how the jacobian is to be calculated and at what point and finally, how to solve the linear system occurring during each iteration. These are given as `lambda` functions. 
 
The lines
@code{.sh}
            typename TrilinosWrappers::NOXSolver<Vector<double>>::AdditionalData additional_data;
            additional_data.abs_tol = target_tolerance;
            additional_data.max_iter = 100;
            TrilinosWrappers::NOXSolver<Vector<double>> nonlinear_solver(
                    additional_data);
 
@endcode
setup the solver and give any basic data needed, such as the tolerance to converge or the maximum number of iterations it can try. 
 
The lines
@code{.sh}
            nonlinear_solver.residual =
                    [&](const Vector<double> &evaluation_point,
                        Vector<double> &residual) {
                        compute_residual(evaluation_point, residual);
                    };
 
@endcode
essentially, defines the function from where the residual will be calculated. Notice here that the `compute_residual()` function is called with `evaluation_point` and `residual`. 
 
Similarly, the lines,
@code{.sh}
            nonlinear_solver.setup_jacobian =
                    [&](const Vector<double> &current_u) {
                        compute_jacobian(current_u);
                    };
 
@endcode
define the function from where the jacobian will be calculated. Note that it takes in a variable `current_u` which will be the `evaluation_point` for the computing the jacobian. Thus, the actual points for evaluating the jacobian and the residual need not be the same. 
 
The lines
@code{.sh}
            nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
                                                       Vector<double> &dst,
                                                       const double tolerance) {
                solve(rhs, dst, tolerance);
            };
 
@endcode
tells the solution needs to be done with the jacobian using the function `solve()`. 
 
Then the line
@code{.sh}
            nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
                                                       Vector<double> &dst,
                                                       const double tolerance) {
                solve(rhs, dst, tolerance);
            };
 
@endcode
actually calls the nonlinear solver.
 
For details as to how these lambda functions work together, please see @ref step_77 "step-77" and also the file `tests/trilinos/step-77-with-nox.cc`. Overall, the variable `present_solution` is presented as an initial guess to the non-linear solver, which performs the iteration and gives the final converged solution in the same variable. Hence, this will be the converged solution for the current step and this assignment happens in the following line
@code{.sh}
        converged_solution = present_solution;
 
@endcode.
The variable `present_solution` is assigned the value of `converged_solution` from the previous time step to serve as an initial guess. This is done in line 45. 
@code{.sh}
            present_solution = converged_solution;
 
@endcode
 
 
### Results
 
The results are essentially the time evolution of the temperature throughout the domain. The first of the pictures below shows the temperature distribution at the final step, i.e. at time @f$t=5@f$. This should be very similar to the figure at the bottom on the page [here](https://www.mathworks.com/help/pde/ug/heat-transfer-problem-with-temperature-dependent-properties.html). We also plot the time evolution of the temperature at a point close to the right edge of the domain indicated by the small magenta dot (close to @f$(0.49, 0.12)@f$) in the second of the pictures below. This is also similar to the second figure at the [bottom of this page](https://www.mathworks.com/help/pde/ug/heat-transfer-problem-with-temperature-dependent-properties.html). There could be minor differences due to the choice of the point. Further, note that, we have plotted in the second of the pictures below the temperature as a function of time steps instead of time. Since the @f$\Delta t@f$ chosen is 0.1, 50 steps maps to @f$t=5@f$ as in the link.
 
![image](../code-gallery/nonlinear-heat_transfer_with_AD_NOX/doc/Images/contour.png)
 
Contour plot of the temperature at the final step
 
![image](../code-gallery/nonlinear-heat_transfer_with_AD_NOX/doc/Images/Temp_evol.png)
 
Evolution of temperature at a point close to the right edge @f$\approx (0.49, 0.12)@f$
 
In closing, we give some ideas as to how the residuals and the jacobians are actually calculated in a finite element setting. 
 
 
### Evaluating the function at previous iterations
 
For solving the linear equation every iteration, we have to evaluate the right hand side functions (@f$f_i^1@f$ and @f$f_i^2@f$) at the values  @f$p@f$ and @f$T@f$ had at their previous iteration. For instance, in the current problem consider the final equation above (we only have only one function @f$f_I@f$ in this case), which is 
@f[
                M_{IJ}T_{J}^{s+1} +     \alpha \Delta t L_I^{s+1}  - M_{IJ}T_{J}^s + \Delta t \left( 1-\alpha \right) L_I^s - \Delta t Q_I = R_I^{s+1}
@f]
Each term is written in the following manner. We write it for only one term to illustrate the details. 
Consider @f$\alpha \Delta t L_I^{s+1} @f$
@f[
\alpha \Delta t \int_{\Omega_e} \left( \psi_{I,i} \left( k (\psi_PT_P) \psi_{J,i}T_J \right) - \psi_I f \right)  dx dy = L_I^{s+1}
@f]
to be specific
@f[
        \alpha \Delta t \int_{\Omega_e} \left( \psi_{I,i} \left( k (\psi_PT_P^{s+1}) \psi_{J,i}T_J^{s+1} \right) - \psi_I f \right)  dx dy = L_I^{s+1}
@f]
again, the term @f$\psi_{J,i} T_J^{{s+1}}@f$ is nothing but the value of gradient of @f$T^{{s+1}}@f$ at whatever point the gradient of @f$\psi_{J}@f$ is evaluated in. We write this gradient as @f$\nabla T^{{s+1}}(x,y)@f$. Similarly, @f$\psi_PT_P^{{s+1}}@f$ is the value of the temperature at the point where @f$\psi_P@f$ is evaluated and we call it @f$T^{{s+1}}@f$. While evaluating these terms to calculate the residual, we will use Gauss quadrature and this evaluation of this integral would become something like
@f[
L_I^{{s+1}} =   \Delta t \alpha \sum_{h} \nabla \psi_I (\xi_h,\eta_h) \cdot k(T^{{s+1}}) \nabla  T^{{s+1}}(\xi_h,\eta_h) w_hJ_h
@f]
where @f$h@f$ runs over the total number of quadrature points in the finite element and @f$w_hJ_h@f$ takes care of appropriate weights needed and other corrections needed to due to transforming the finite element onto a master element (See any finite element book for details). All other terms can also be calculated in a similar manner to obtain the per-element residual from the first equation in this section. This has to assembled over all finite elements to get the global residual.
 
Similar evaluations have to be done to obtain the Jacobian as well. In the current case, the per element Jacobian is given by
@f[
        J^e_{IQ} = \frac{\partial R_I^{s+1}}{\partial T_{Q}^{s+1}} = M_{IQ} + \alpha \Delta t \frac{\partial L_I^{s+1}}{\partial T_Q}
@f]
which is written as
@f[
        J_{IQ}^e = \int_{\Omega_e}\rho C_p \psi_I  \psi_J  dx dy  +
        \alpha \Delta t \int_{\Omega_e} \left( \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 \left( \psi_Y T_Y) \right)^2 \right)  dx dy 
@f]
which are evaluated using gauss quadrature by summing the functions after evaluation at the gauss points. Once again, the terms @f$\psi_{J,i}T_J@f$ and @f$\psi_R T_R@f$ are the values of the gradients of the temperature and the temperature itself at the Gauss points. 
 
 
 
 
## The dealii steps
To understand this application, @ref step_71 "step-71", @ref step_72 "step-72" and @ref step_77 "step-77" are needed. 
 
 
<a name="ann-mesh/README.md"></a>
<h1>Annotated version of mesh/README.md</h1>
The following copyright notice applies to the mesh files in this directory:
@code{.sh}
/*
 * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 * Copyright (C) 2024 by Narasimhan Swaminathan
 *
 * This file is part of the deal.II code gallery.
 *
 * -----------------------------------------------------------------------------
 */
@endcode
 
 
<a name="ann-include/allheaders.h"></a>
<h1>Annotated version of include/allheaders.h</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Narasimhan Swaminathan
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #ifndef __ALLHEADERS_H_INCLUDED__
 *   #define __ALLHEADERS_H_INCLUDED__
 * @endcode
 * 
 * These are some of the header files for this programme. Most of them
 * are common from previous steps.
 * 
 * @code
 *   #include <deal.II/base/quadrature_lib.h>
 *   #include <deal.II/base/function_lib.h>
 *   #include <deal.II/lac/generic_linear_algebra.h>
 *   #include <deal.II/base/function.h>
 *   
 *   #include <deal.II/dofs/dof_handler.h>
 *   #include <deal.II/dofs/dof_tools.h>
 *   
 *   #include <deal.II/fe/fe_values.h>
 *   #include <deal.II/fe/fe.h>
 *   
 *   #include <deal.II/grid/tria.h>
 *   #include <deal.II/grid/grid_generator.h>
 *   
 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
 *   #include <deal.II/lac/full_matrix.h>
 *   #include <deal.II/lac/precondition.h>
 *   #include <deal.II/lac/solver_cg.h>
 *   #include <deal.II/lac/sparse_matrix.h>
 *   #include <deal.II/lac/vector.h>
 *   #include <deal.II/grid/grid_tools.h>
 *   
 *   #include <deal.II/numerics/data_out.h>
 *   #include <deal.II/numerics/vector_tools.h>
 *   #include <deal.II/fe/fe_q.h>
 *   #include <deal.II/grid/grid_out.h>
 *   #include <deal.II/grid/grid_in.h>
 *   
 *   #include <deal.II/fe/fe_system.h>
 *   
 *   #include <deal.II/grid/tria_accessor.h>
 *   #include <deal.II/grid/tria_iterator.h>
 *   #include <deal.II/dofs/dof_handler.h>
 *   #include <deal.II/dofs/dof_accessor.h>
 *   #include <deal.II/dofs/dof_tools.h>
 *   #include <deal.II/fe/fe_values.h>
 *   #include <deal.II/fe/fe.h>
 *   #include <deal.II/numerics/matrix_tools.h>
 *   #include <deal.II/numerics/vector_tools.h>
 *   #include <deal.II/lac/vector.h>
 *   #include <deal.II/base/quadrature.h>
 *   #include <deal.II/distributed/tria.h>
 *   #include <deal.II/distributed/grid_refinement.h>
 *   
 *   #include <deal.II/lac/sparse_direct.h>
 *   #include <deal.II/base/timer.h>
 *   #include <deal.II/base/utilities.h>
 *   
 *   #include <deal.II/base/exceptions.h>
 *   #include <deal.II/base/geometric_utilities.h>
 *   #include <deal.II/base/conditional_ostream.h>
 *   #include <deal.II/base/mpi.h>
 *   
 *   #include <deal.II/grid/grid_tools.h>
 *   #include <deal.II/dofs/dof_renumbering.h>
 *   #include <deal.II/numerics/solution_transfer.h>
 *   #include <deal.II/base/index_set.h>
 *   #include <deal.II/lac/sparsity_tools.h>
 *   #include <deal.II/fe/fe_values_extractors.h>
 * @endcode
 * 
 * Automatic differentiation is carried out with TRILINOS Automatic differentiation scheme/
 * This is invoked with the following include.
 * 
 * @code
 *   #include <deal.II/differentiation/ad.h>
 * @endcode
 * 
 * We use Trilinos wrappers based NOX to solve the non-linear equations.
 * 
 * @code
 *   #include <deal.II/trilinos/nox.h>
 *   
 *   #include <math.h>
 *   #include <fstream>
 *   #include <iostream>
 *   
 *   using namespace dealii;
 *   #endif //__ALLHEADERS_H_INCLUDED__
 * @endcode
 
 
<a name="ann-include/nonlinear_heat.h"></a>
<h1>Annotated version of include/nonlinear_heat.h</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Narasimhan Swaminathan
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #ifndef __MAIN_ALL_HEADER_H_INCLUDED__
 *   #define __MAIN_ALL_HEADER_H_INCLUDED__
 *   #include "allheaders.h"
 *   
 *   class nonlinear_heat
 *   {
 *   public:
 *       
 *       nonlinear_heat();
 *       ~nonlinear_heat();
 *       
 *       void run(); 
 *       
 *       const double delta_t;
 *       
 *       const double alpha;
 *       
 *       const double tot_time;
 *   
 *       
 *       const double a;
 *       const double b;
 *       const double c;
 *       
 *       const double Cp;
 *       
 *       const double rho;
 *       
 *       double time;
 *   private:
 *       
 *       void compute_jacobian(const Vector<double> &evaluation_point);
 *       
 *       void compute_residual(const Vector<double> &evaluation_point, Vector<double> & residual);
 *       
 *       void setup_system(unsigned int time_step);
 *       
 *       void solve(const Vector<double> &rhs, Vector<double> & solution, const double tolerance);
 *       
 *       void output_results(unsigned int prn) const;
 *       
 *       void set_boundary_conditions(double time);
 *   
 *       Triangulation<2> triangulation;
 *       DoFHandler<2>         dof_handler; 
 *       FESystem<2>       fe;
 *       SparsityPattern     sparsity_pattern;
 *       SparseMatrix<double> system_matrix; /*Matrix holding the global Jacobian*/
 *       
 *       std::unique_ptr<SparseDirectUMFPACK> matrix_factorization;
 *       
 *       Vector<double> converged_solution;/* Converged solution in the previous time step */
 *       
 *       Vector<double> present_solution;/* Converged solution in the previous time step */
 *   };
 *   
 *   
 *   class Initialcondition : public Function<2>
 *   {
 *   public:
 *       Initialcondition(): Function<2>(1)
 *      {}
 * @endcode
 * 
 * Returns the initial values.
 * 
 * @code
 *       virtual double value(const Point<2> &p,
 *                                 const unsigned int component =0) const override;
 *   };
 *   
 *   
 *   class Boundary_values_left:public Function<2>
 *   {
 *   public:
 *       Boundary_values_left(): Function<2>(1)
 *       {}
 *       virtual double value(const Point<2> & p,const unsigned int component = 0) const override;
 *   }; 
 *   #endif //__MAIN_ALL_HEADER_H_INCLUDED__
 * @endcode
 
 
<a name="ann-nonlinear_heat.cc"></a>
<h1>Annotated version of nonlinear_heat.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Narasimhan Swaminathan
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "allheaders.h"
 *   #include "nonlinear_heat.h"
 *   int main()
 *   {
 *     try
 *       {
 *         nonlinear_heat nlheat; 
 *         nlheat.run(); 
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *   
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *   
 *     return 0;
 *   }
 * @endcode
 
 
<a name="ann-source/boundary_values.cc"></a>
<h1>Annotated version of source/boundary_values.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Narasimhan Swaminathan
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "nonlinear_heat.h"
 *   
 *   double Boundary_values_left::value(const Point<2> & /*p*/, const unsigned int /*comp*/) const
 *   {
 *       return 100;
 *   
 *       
 * @endcode
 * 
 * nonlinear_heat nlheat;
 * double total_time = nlheat.tot_time;
 * return this->get_time() * 100.0/total_time;
 * 
 * @code
 *   }
 * @endcode
 
 
<a name="ann-source/compute_jacobian.cc"></a>
<h1>Annotated version of source/compute_jacobian.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Narasimhan Swaminathan
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "allheaders.h"
 *   #include "nonlinear_heat.h"
 *   
 *   void nonlinear_heat::compute_jacobian(const Vector<double> &evaluation_point)
 *   {
 *       const QGauss<2> quadrature_formula(
 *               fe.degree + 1);
 *   
 *       const QGauss<1> face_quadrature_formula(fe.degree+1);
 *   
 *       FEValues<2> fe_values(fe,
 *                             quadrature_formula,
 *                             update_values | update_gradients |
 *                             update_quadrature_points |
 *                             update_JxW_values);
 *   
 *       FEFaceValues<2> fe_face_values(fe,face_quadrature_formula,update_values|update_quadrature_points|
 *                                                                 update_normal_vectors | update_JxW_values);
 *   
 *   
 *       const unsigned int dofs_per_cell = fe.dofs_per_cell;
 *       const unsigned int n_q_points = quadrature_formula.size();
 *       const unsigned int n_q_f_points = face_quadrature_formula.size();
 *       Vector<double> cell_rhs(dofs_per_cell);
 *       FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell); 
 *       using ADHelper = Differentiation::AD::ResidualLinearization<Differentiation::AD::NumberTypes::sacado_dfad,double>;
 *       using ADNumberType = typename ADHelper::ad_type;
 *       const FEValuesExtractors::Scalar t(0);
 *       
 *       system_matrix = 0.0;
 *       std::vector<types::global_dof_index> local_dof_indices(
 *               dofs_per_cell);
 *       /*==================================================================*/
 *       std::vector<double> consol(n_q_points);
 *       std::vector<ADNumberType> old_solution(
 *               n_q_points);
 *   
 *       std::vector<Tensor<1, 2>> consol_grad(
 *               n_q_points);
 *       std::vector<Tensor<1, 2,ADNumberType>> old_solution_grad(
 *               n_q_points);
 *   
 *       for (const auto &cell: dof_handler.active_cell_iterators())
 *       {
 *           cell_rhs = 0;
 *           cell_matrix = 0;
 *           fe_values.reinit(cell);
 *           cell->get_dof_indices(local_dof_indices);
 *           const unsigned int n_independent_variables = local_dof_indices.size();
 *           const unsigned int n_dependent_variables = dofs_per_cell;
 *           ADHelper ad_helper(n_independent_variables, n_dependent_variables);
 *           ad_helper.register_dof_values(evaluation_point,local_dof_indices);
 *           const std::vector<ADNumberType> &dof_values_ad = ad_helper.get_sensitive_dof_values();
 *           fe_values[t].get_function_values_from_local_dof_values(dof_values_ad,
 *                                            old_solution);
 *           fe_values[t].get_function_gradients_from_local_dof_values(dof_values_ad,
 *                                               old_solution_grad);
 *   
 *           fe_values[t].get_function_values(converged_solution, consol);
 *           fe_values[t].get_function_gradients(converged_solution, consol_grad);
 *           std::vector<ADNumberType> residual_ad(n_dependent_variables,
 *                                                 ADNumberType(0.0));
 *           for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
 *           {
 *               for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *               {
 *   
 *                   ADNumberType MijTjcurr = Cp*rho*fe_values[t].value(i,q_index)*old_solution[q_index];
 *                   ADNumberType MijTjprev = Cp*rho*fe_values[t].value(i,q_index)*consol[q_index];
 *                   ADNumberType k_curr = a + b*old_solution[q_index] + c*std::pow(old_solution[q_index],2);
 *                   ADNumberType k_prev = a + b*consol[q_index] + c*std::pow(consol[q_index],2);
 *                   ADNumberType Licurr =  alpha * delta_t *  (fe_values[t].gradient(i,q_index)*k_curr*old_solution_grad[q_index]);
 *                   ADNumberType Liprev =  (1-alpha) * delta_t *  (fe_values[t].gradient(i,q_index)*k_prev*consol_grad[q_index]);
 *                   residual_ad[i] +=  (MijTjcurr+Licurr-MijTjprev+Liprev)*fe_values.JxW(q_index);
 *               }
 *           }
 *   
 *   
 *           for (unsigned int face_number = 0;face_number<GeometryInfo<2>::faces_per_cell; ++face_number)
 *           {
 *   
 *               if (cell->face(face_number)->boundary_id() == 3)
 *               {
 *                   fe_face_values.reinit(cell, face_number);
 *                   for (unsigned int q_point=0;q_point<n_q_f_points;++q_point)
 *                   {
 *                       for (unsigned int i =0;i<dofs_per_cell;++i)
 *                           residual_ad[i]+= -delta_t*(-10)*fe_face_values[t].value(i,q_point)*fe_face_values.JxW(q_point);
 *   
 *                   }
 *               }
 *           }
 *   
 *           ad_helper.register_residual_vector(residual_ad);
 *           ad_helper.compute_residual(cell_rhs);
 *           
 *           ad_helper.compute_linearization(cell_matrix);
 *           cell->get_dof_indices(local_dof_indices);
 *           
 *   
 *           for (unsigned int i =0;i < dofs_per_cell; ++i){
 *               for(unsigned int j = 0;j < dofs_per_cell;++j){
 *                   system_matrix.add(local_dof_indices[i],local_dof_indices[j],cell_matrix(i,j));
 *               }
 *           }
 *       }
 *       
 *       std::map<types::global_dof_index, double> boundary_values;
 *       VectorTools::interpolate_boundary_values(dof_handler,
 *                                                1,
 *                                                Functions::ZeroFunction<2>(1),
 *                                                boundary_values);
 *   
 *       Vector<double> dummy_solution(dof_handler.n_dofs());
 *       Vector<double> dummy_rhs(dof_handler.n_dofs());
 *       MatrixTools::apply_boundary_values(boundary_values,
 *                                          system_matrix,
 *                                          dummy_solution,
 *                                          dummy_rhs);
 *   
 *       {
 *           std::cout << "  Factorizing Jacobian matrix" << std::endl;
 *           matrix_factorization = std::make_unique<SparseDirectUMFPACK>();
 *           matrix_factorization->factorize(system_matrix);
 *       }
 *   }
 * @endcode
 
 
<a name="ann-source/compute_residual.cc"></a>
<h1>Annotated version of source/compute_residual.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Narasimhan Swaminathan
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "allheaders.h"
 *   #include "nonlinear_heat.h"
 *   
 *   void nonlinear_heat::compute_residual(const Vector<double> & evaluation_point, Vector<double> & residual)
 *   {
 *       
 *       const QGauss<2> quadrature_formula(
 *               fe.degree + 1);
 *       const QGauss<1> face_quadrature_formula(fe.degree+1); //Define quadrature for integration over faces */
 *       FEValues<2> fe_values(fe,
 *                             quadrature_formula,
 *                             update_values | update_gradients |
 *                             update_quadrature_points |
 *                             update_JxW_values); 
 *   
 *       FEFaceValues<2> fe_face_values(fe,face_quadrature_formula,update_values|update_quadrature_points|
 *                                                                 update_normal_vectors | update_JxW_values);
 *       const unsigned int dofs_per_cell = fe.dofs_per_cell; 
 *       const unsigned int n_q_points = quadrature_formula.size();
 *       const unsigned int n_q_f_points = face_quadrature_formula.size();
 *       
 *       Vector<double> cell_rhs(dofs_per_cell); 
 *   
 *       
 *       using ADHelper = Differentiation::AD::ResidualLinearization<Differentiation::AD::NumberTypes::sacado_dfad,double>;
 *       using ADNumberType = typename ADHelper::ad_type;
 *       
 *       const FEValuesExtractors::Scalar t(
 *               0);
 *       std::vector<types::global_dof_index> local_dof_indices(
 *               dofs_per_cell); 
 *       
 *   
 *       std::vector<double> consol(n_q_points); /* Converged solution at the Gauss points from the previous time step*/
 *       std::vector<ADNumberType> old_solution(
 *               n_q_points); /* Current solution at the Gauss points at this iteration for the current time*/
 *       
 *       std::vector<Tensor<1, 2>> consol_grad(
 *               n_q_points); /* Converged gradients of the solutions at the Gauss points from the previous time step */
 *       std::vector<Tensor<1, 2,ADNumberType>> old_solution_grad(
 *               n_q_points);
 *   
 *       ComponentMask t_mask = fe.component_mask(t);
 *       
 *       residual = 0.0;
 *       for (const auto &cell: dof_handler.active_cell_iterators())
 *       {
 *           cell_rhs = 0;
 *           fe_values.reinit(cell);
 *           cell->get_dof_indices(local_dof_indices);
 *           const unsigned int n_independent_variables = local_dof_indices.size();
 *           const unsigned int n_dependent_variables = dofs_per_cell;
 *           ADHelper ad_helper(n_independent_variables, n_dependent_variables);
 *           
 *           ad_helper.register_dof_values(evaluation_point,local_dof_indices);
 *   
 *           const std::vector<ADNumberType> &dof_values_ad = ad_helper.get_sensitive_dof_values();
 *           fe_values[t].get_function_values_from_local_dof_values(dof_values_ad,
 *                                            old_solution);
 *           fe_values[t].get_function_gradients_from_local_dof_values(dof_values_ad,
 *                                               old_solution_grad);
 *           
 *           fe_values[t].get_function_values(converged_solution, consol);
 *           fe_values[t].get_function_gradients(converged_solution, consol_grad);
 *           
 *           std::vector<ADNumberType> residual_ad(n_dependent_variables,
 *                                                 ADNumberType(0.0));
 *           for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
 *           {
 *               for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *               {
 *                   
 *                   ADNumberType MijTjcurr = Cp*rho*fe_values[t].value(i,q_index)*old_solution[q_index];
 *                   ADNumberType MijTjprev = Cp*rho*fe_values[t].value(i,q_index)*consol[q_index];
 *                   ADNumberType k_curr = a + b*old_solution[q_index] + c*std::pow(old_solution[q_index],2);
 *                   ADNumberType k_prev = a + b*consol[q_index] + c*std::pow(consol[q_index],2);
 *                   ADNumberType Licurr =  alpha * delta_t *  (fe_values[t].gradient(i,q_index)*k_curr*old_solution_grad[q_index]);
 *                   ADNumberType Liprev =  (1-alpha) * delta_t *  (fe_values[t].gradient(i,q_index)*k_prev*consol_grad[q_index]);
 *                   residual_ad[i] +=  (MijTjcurr+Licurr-MijTjprev+Liprev)*fe_values.JxW(q_index);
 *               }
 *           }
 *               
 *   
 *           for (unsigned int face_number = 0;face_number<GeometryInfo<2>::faces_per_cell; ++face_number)
 *           {
 *   
 *               if (cell->face(face_number)->boundary_id() == 3)
 *               {
 *                   fe_face_values.reinit(cell, face_number);
 *                   for (unsigned int q_point=0;q_point<n_q_f_points;++q_point)
 *                   {
 *                       for (unsigned int i =0;i<dofs_per_cell;++i)
 *                           
 *                           residual_ad[i]+= -delta_t*(-10)*fe_face_values[t].value(i,q_point)*fe_face_values.JxW(q_point);
 *   
 *                   }
 *               }
 *           }
 *           
 *           ad_helper.register_residual_vector(residual_ad);
 *           
 *           ad_helper.compute_residual(cell_rhs);
 *           cell->get_dof_indices(local_dof_indices);
 *   
 *           for (unsigned int i =0;i < dofs_per_cell; ++i)
 *               residual(local_dof_indices[i])+= cell_rhs(i);
 *   
 *       }
 *   
 *       for(const types::global_dof_index i: DoFTools::extract_boundary_dofs(dof_handler,t_mask,{1}))
 *           residual(i) = 0;
 *   
 *       std::cout << " The Norm is :: = " << residual.l2_norm() << std::endl;
 *   }
 * @endcode
 
 
<a name="ann-source/initial_conditions.cc"></a>
<h1>Annotated version of source/initial_conditions.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Narasimhan Swaminathan
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "nonlinear_heat.h"
 *   
 *   
 *   double Initialcondition::value(const Point<2> & /*p*/, const unsigned int /*comp*/) const
 *   {
 *       
 *       return 0.0;
 *   }
 * @endcode
 
 
<a name="ann-source/nonlinear_heat_cons_des.cc"></a>
<h1>Annotated version of source/nonlinear_heat_cons_des.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Narasimhan Swaminathan
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "nonlinear_heat.h"
 *   
 *   
 *   nonlinear_heat::nonlinear_heat ()
 *       :delta_t(0.1),
 *       alpha(0.5),
 *       tot_time(5),
 *       a(0.3),
 *       b(0.003),
 *       c(0),
 *       Cp(1),
 *       rho(1),
 *       dof_handler(triangulation),
 *       fe(FE_Q<2>(1), 1)
 *   {}
 *   
 *   
 *   nonlinear_heat::~nonlinear_heat()
 *   {
 *       dof_handler.clear();
 *   }
 *       
 *   
 * @endcode
 
 
<a name="ann-source/output_results.cc"></a>
<h1>Annotated version of source/output_results.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Narasimhan Swaminathan
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "nonlinear_heat.h"
 *   
 *   
 *   void nonlinear_heat::output_results(unsigned int prn) const
 *   {
 *       DataOut<2> data_out;
 *       data_out.attach_dof_handler(dof_handler);
 *       std::vector<std::string> solution_names;
 *       solution_names.emplace_back ("Temperature");
 *       data_out.add_data_vector(converged_solution, solution_names);
 *       data_out.build_patches();
 *       const std::string filename =
 *               "output/solution-" + Utilities::int_to_string(prn , 3) + ".vtu";
 *       std::ofstream output(filename);
 *       data_out.write_vtu(output);
 *   }
 * @endcode
 
 
<a name="ann-source/set_boundary_conditions.cc"></a>
<h1>Annotated version of source/set_boundary_conditions.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Narasimhan Swaminathan
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "nonlinear_heat.h"
 *   #include "allheaders.h"
 *   
 *   
 *   void nonlinear_heat::set_boundary_conditions(double time)
 *   {
 *       Boundary_values_left bl_left;
 *       bl_left.set_time(time);
 *       std::map<types::global_dof_index,double> boundary_values;
 *       VectorTools::interpolate_boundary_values(dof_handler,
 *               1,
 *               bl_left,
 *               boundary_values);
 *   
 *       for (auto &boundary_value: boundary_values)
 *           present_solution(boundary_value.first) = boundary_value.second;
 *   }
 * @endcode
 
 
<a name="ann-source/setup_system.cc"></a>
<h1>Annotated version of source/setup_system.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Narasimhan Swaminathan
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "allheaders.h"
 *   #include "nonlinear_heat.h"
 *   
 *   
 *   void nonlinear_heat::setup_system(unsigned int time_step)
 *   {
 *        if (time_step ==0) {
 *           dof_handler.distribute_dofs(fe);
 *            converged_solution.reinit(dof_handler.n_dofs());
 *            present_solution.reinit(dof_handler.n_dofs());
 *        }
 *   
 *       DynamicSparsityPattern dsp(dof_handler.n_dofs());
 *       DoFTools::make_sparsity_pattern(dof_handler, dsp);
 *       sparsity_pattern.copy_from(dsp);
 *       system_matrix.reinit(sparsity_pattern);
 *       matrix_factorization.reset();
 *   }
 * @endcode
 
 
<a name="ann-source/solve_and_run.cc"></a>
<h1>Annotated version of source/solve_and_run.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Narasimhan Swaminathan
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "nonlinear_heat.h"
 *   
 *   
 *   void nonlinear_heat::solve(const Vector<double> & rhs, Vector<double> & solution, const double /*tolerance*/)
 *   {
 *       std::cout << "  Solving linear system" << std::endl;
 *       matrix_factorization->vmult(solution, rhs);
 *   }
 *   
 *   void nonlinear_heat::run()
 *   {
 *       GridIn<2> gridin;
 *       gridin.attach_triangulation(triangulation);
 *       std::ifstream f("mesh/mesh.msh");
 *       gridin.read_msh(f);
 *       triangulation.refine_global(1);
 *   
 *       double time = 0;
 *       unsigned int timestep_number = 0;
 *       unsigned int prn =0;
 *   
 *   
 *       while (time <=tot_time)
 *       {
 *           if(time ==0)
 *           {
 *               setup_system(timestep_number);
 *   
 *               VectorTools::interpolate(dof_handler,
 *                                        Initialcondition(),
 *                                        present_solution);
 *   
 *               VectorTools::interpolate(dof_handler,
 *                                        Initialcondition(),
 *                                        converged_solution);
 *           }
 *           else
 *           {
 *               
 *               present_solution = converged_solution;
 *           }
 *           std::cout<<">>>>> Time now  is: "<<time <<std::endl;
 *           std::cout<<">>>>> Time step is:"<<timestep_number<<std::endl;
 *           set_boundary_conditions(time);
 *           {
 *   
 *               const double target_tolerance = 1e-3;
 *               
 *               typename TrilinosWrappers::NOXSolver<Vector<double>>::AdditionalData additional_data;
 *               additional_data.abs_tol = target_tolerance;
 *               additional_data.max_iter = 100;
 *               TrilinosWrappers::NOXSolver<Vector<double>> nonlinear_solver(
 *                       additional_data);
 *   
 *               
 *               nonlinear_solver.residual =
 *                       [&](const Vector<double> &evaluation_point,
 *                           Vector<double> &residual) {
 *                           compute_residual(evaluation_point, residual);
 *                       };
 *   
 *               
 *   
 *               nonlinear_solver.setup_jacobian =
 *                       [&](const Vector<double> &current_u) {
 *                           compute_jacobian(current_u);
 *                       };
 *               
 *   
 *               nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
 *                                                          Vector<double> &dst,
 *                                                          const double tolerance) {
 *                   solve(rhs, dst, tolerance);
 *               };
 *               
 *               nonlinear_solver.solve(present_solution);
 *           }
 *           
 *           converged_solution = present_solution;
 *   
 *           if(timestep_number % 1 == 0) {
 *   
 *               output_results(prn);
 *               prn = prn +1;
 *           }
 *           timestep_number++;
 *           time=time+delta_t;
 *   
 *       }
 *   }
 * @endcode
 
 
*/