Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
nonlinear-heat_transfer_with_AD_NOX.h
Go to the documentation of this file.
1
246#include "allheaders.h"
247#include "nonlinear_heat.h"
254void nonlinear_heat::compute_residual(const Vector<double> & evaluation_point, Vector<double> & residual)
255{
260 const QGauss<2> quadrature_formula(
261 fe.degree + 1);
262 const QGauss<1> face_quadrature_formula(fe.degree+1); //Define quadrature for integration over faces */
263 FEValues<2> fe_values(fe,
264 quadrature_formula,
269 FEFaceValues<2> fe_face_values(fe,face_quadrature_formula,update_values|update_quadrature_points|
271 const unsigned int dofs_per_cell = fe.dofs_per_cell;
272 const unsigned int n_q_points = quadrature_formula.size();
273 const unsigned int n_q_f_points = face_quadrature_formula.size();
278 Vector<double> cell_rhs(dofs_per_cell);
280@endcode
281should 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.
282
283The 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.)
284
285@code{.sh}
287 using ADNumberType = typename ADHelper::ad_type;
288
289@endcode
290
291The following lines are once again straightforward to understand.
292@code{.sh}
294 0);
295 std::vector<types::global_dof_index> local_dof_indices(
296 dofs_per_cell);
301 std::vector<double> consol(n_q_points); /* Converged solution at the Gauss points from the previous time step*/
302
303@endcode
304We 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$).
305
306
307In following lines the variables `old_solution` and `old_solution_grad` are
308crucial. 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
309@f[
310 \text{old_solution} =
311 \begin{bmatrix}
312 \psi_i (\zeta_1,\eta_1) d_i\\
313 \psi_i (\zeta_2,\eta_2) d_i\\
314 \psi_i (\zeta_3,\eta_3) d_i\\
315 \psi_i (\zeta_4,\eta_4) d_i
316 \end{bmatrix}
317@f]
318where @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.
319@code{.sh}
320 std::vector<ADNumberType> old_solution(
321 n_q_points); /* Current solution at the Gauss points at this iteration for the current time*/
325 std::vector<Tensor<1, 2>> consol_grad(
326 n_q_points); /* Converged gradients of the solutions at the Gauss points from the previous time step */
327 std::vector<Tensor<1, 2,ADNumberType>> old_solution_grad(
328 n_q_points);
329
330
331@endcode
332
333The lines
334@code{.sh}
335 ComponentMask t_mask = fe.component_mask(t);
339 residual = 0.0;
340
341@endcode
342are straightforward to understand. The vector `residual` will hold the *numerical* value of the residual and is hence initialized to be zero.
343
344The lines
345@code{.sh}
346 for (const auto &cell: dof_handler.active_cell_iterators())
347 {
348 cell_rhs = 0;
349 fe_values.reinit(cell);
350 cell->get_dof_indices(local_dof_indices);
351 const unsigned int n_independent_variables = local_dof_indices.size();
352 const unsigned int n_dependent_variables = dofs_per_cell;
353
354@endcode
355should be clear. In the following line we inform how many independent and dependent variables we will have.
356@code{.sh}
357 ADHelper ad_helper(n_independent_variables, n_dependent_variables);
358
359@endcode
360then, 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.
361@code{.sh}
362 ad_helper.register_dof_values(evaluation_point,local_dof_indices);
363
364@endcode
365
366Now the lines
367@code{.sh}
368 const std::vector<ADNumberType> &dof_values_ad = ad_helper.get_sensitive_dof_values();
369 fe_values[t].get_function_values_from_local_dof_values(dof_values_ad,
370 old_solution);
371 fe_values[t].get_function_gradients_from_local_dof_values(dof_values_ad,
372 old_solution_grad);
373
374@endcode
375the 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`.
376
377In the lines
378@code{.sh}
379 fe_values[t].get_function_values(converged_solution, consol);
380 fe_values[t].get_function_gradients(converged_solution, consol_grad);
381
382@endcode
383the values of the `converged_solution` and its gradients get stored in the variables `consol` and `con_sol_grad`.
384
385Then in the lines
386@code{.sh}
387 for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
388 {
389 for (unsigned int i = 0; i < dofs_per_cell; ++i)
390 {
394 ADNumberType MijTjcurr = Cp*rho*fe_values[t].value(i,q_index)*old_solution[q_index];
395 ADNumberType MijTjprev = Cp*rho*fe_values[t].value(i,q_index)*consol[q_index];
396 ADNumberType k_curr = a + b*old_solution[q_index] + c*std::pow(old_solution[q_index],2);
397 ADNumberType k_prev = a + b*consol[q_index] + c*std::pow(consol[q_index],2);
398 ADNumberType Licurr = alpha * delta_t * (fe_values[t].gradient(i,q_index)*k_curr*old_solution_grad[q_index]);
399 ADNumberType Liprev = (1-alpha) * delta_t * (fe_values[t].gradient(i,q_index)*k_prev*consol_grad[q_index]);
400 residual_ad[i] += (MijTjcurr+Licurr-MijTjprev+Liprev)*fe_values.JxW(q_index);
401 }
402 }
433 residual_ad[i]+= -delta_t*(-10)*fe_face_values[t].value(i,q_point)*fe_face_values.JxW(q_point);
434
435 }
436 }
437 }
438
439@endcode
440
441Now
442@code{.sh}
443 ad_helper.register_residual_vector(residual_ad);
444
445@endcode
446actually tells the `ad_helper` that this is the residual it should use, and the line
447@code{.sh}
448 ad_helper.compute_residual(cell_rhs);
449
450@endcode
451evaluates the residual at the appropriate point, which is the `evaluation_point`. Then
452@code{.sh}
453 cell->get_dof_indices(local_dof_indices);
454
455 for (unsigned int i =0;i < dofs_per_cell; ++i)
456 residual(local_dof_indices[i])+= cell_rhs(i);
457
458 }
459
460@endcode
461are used to calculate the global `residual` from the local `cell_rhs`. Now `residual` is numerical and not of the special number type. Finally, lines
462@code{.sh}
463 for(const types::global_dof_index i: DoFTools::extract_boundary_dofs(dof_handler,t_mask,{1}))
464 residual(i) = 0;
465
466 std::cout << " The Norm is :: = " << residual.l2_norm() << std::endl;
467}
468
469@endcode
470are routine, in particular, it makes the global degrees of freedom, where the dirichlet boundary conditions are applied to be zero.
471
472
473#### The compute_jacobian() function
474
475The `compute_jacobian()` only takes in the `evaluation_point` as an input and is identical to the
476`compute_residual()` except for some minor modification. In particular the line
477@code{.sh}
478 ad_helper.compute_linearization(cell_matrix);
479
480@endcode
481computes the actual jacobian by performing the automatic differentiation and evaluates it at the `evaluation_point`. The remaining lines are routine and should be straightforward.
482
483
484### General ideas involved in solving coupled nonlinear equations using Newton Raphson's technique
485Consider that we have the equations
486@f{eqnarray}{
487 f_i^1(p_j^{{s+1}}, T_j^{{s+1}})=0 \\
488 f_i^2(p_j^{{s+1}}, T_j^{{s+1}})=0
489@f}
490are 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.
491
492We 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
493@f{eqnarray}{
494 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 \\
495 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
496@f}
497Also, 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.
498@f[
499 J_{pq}=-
500 \begin{bmatrix}
501 \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}\\
502
503 \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}\\
504
505 \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}\\
506
507 \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}\\
508 %================================================================
509 \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}\\
510
511 \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}\\
512
513 \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}\\
514
515 \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}
516 \end{bmatrix}
517@f]
518We 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
519@f[
520 \begin{bmatrix}
521 p_1^{s+1}\\
522 p_2^{s+1}\\
523 p_3^{s+1}\\
524 p_4^{s+1}\\
525 T_1^{s+1}\\
526 T_2^{s+1}\\
527 T_3^{s+1}\\
528 T_4^{s+1}
529 \end{bmatrix}^{k+1}
530 =
531 \begin{bmatrix}
532 p_1^{s+1}\\
533 p_2^{s+1}\\
534 p_3^{s+1}\\
535 p_4^{s+1}\\
536 T_1^{s+1}\\
537 T_2^{s+1}\\
538 T_3^{s+1}\\
539 T_4^{s+1}
540 \end{bmatrix}^{k}-
541 \begin{bmatrix}
542 \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}\\
543
544 \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}\\
545
546 \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}\\
547
548 \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}\\
549 %================================================================
550 \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}\\
551
552 \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}\\
553
554 \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}\\
555
556 \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}
557 \end{bmatrix}^{-1}
558 \begin{bmatrix}
559 f_1^{1}\\
560 f_2^{1}\\
561 f_3^{1}\\
562 f_4^{1}\\
563 f_1^{2}\\
564 f_2^{2}\\
565 f_3^{2}\\
566 f_4^{2}
567 \end{bmatrix}^{k}
568@f]
569The @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$.
570@f[
571 \begin{bmatrix}
572 \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}\\
573
574 \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}\\
575
576 \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}\\
577
578 \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}\\
579 %================================================================
580 \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}\\
581
582 \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}\\
583
584 \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}\\
585
586 \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}
587 \end{bmatrix}
588 \begin{bmatrix}
589 \Delta p_1^{s+1}\\
590 \Delta p_2^{s+1}\\
591 \Delta p_3^{s+1}\\
592 \Delta p_4^{s+1}\\
593 \Delta T_1^{s+1}\\
594 \Delta T_2^{s+1}\\
595 \Delta T_3^{s+1}\\
596 \Delta T_4^{s+1}
597 \end{bmatrix}=
598 \begin{bmatrix}
599 f_1^{1}\\
600 f_2^{1}\\
601 f_3^{1}\\
602 f_4^{1}\\
603 f_1^{2}\\
604 f_2^{2}\\
605 f_3^{2}\\
606 f_4^{2}
607 \end{bmatrix}^{k}
608@f]
609This is written as
610@f[
611 J_{pq}\{\Delta solution\}^{k+1}=\{function\}^{k}
612@f]
613Then add this to the value of
614@f[
615 \begin{bmatrix}
616 p_1^{s+1}\\
617 p_2^{s+1}\\
618 p_3^{s+1}\\
619 p_4^{s+1}\\
620 T_1^{s+1}\\
621 T_2^{s+1}\\
622 T_3^{s+1}\\
623 T_4^{s+1}
624 \end{bmatrix}^{k}
625@f]
626We 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
627@f[
628 ||(function)^k||<\delta^0 ||function^0||
629@f]
630or
631@f[
632 ||(function)^k||<tol
633@f]
634where @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.
635
636
637#### The run() function
638
639The function which implements the nonlinear solver is the `run()` function for every time step. In particular the lines
640@code{.sh}
641 typename TrilinosWrappers::NOXSolver<Vector<double>>::AdditionalData additional_data;
642 additional_data.abs_tol = target_tolerance;
643 additional_data.max_iter = 100;
645 additional_data);
646
650 nonlinear_solver.residual =
651 [&](const Vector<double> &evaluation_point,
652 Vector<double> &residual) {
653 compute_residual(evaluation_point, residual);
654 };
655
660 nonlinear_solver.setup_jacobian =
661 [&](const Vector<double> &current_u) {
662 compute_jacobian(current_u);
663 };
668 nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
669 Vector<double> &dst,
670 const double tolerance) {
671 solve(rhs, dst, tolerance);
672 };
677 nonlinear_solver.solve(present_solution);
678 }
679
680@endcode
681Here, 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.
682
683The lines
684@code{.sh}
685 typename TrilinosWrappers::NOXSolver<Vector<double>>::AdditionalData additional_data;
686 additional_data.abs_tol = target_tolerance;
687 additional_data.max_iter = 100;
689 additional_data);
690
691@endcode
692setup the solver and give any basic data needed, such as the tolerance to converge or the maximum number of iterations it can try.
693
694The lines
695@code{.sh}
696 nonlinear_solver.residual =
697 [&](const Vector<double> &evaluation_point,
698 Vector<double> &residual) {
699 compute_residual(evaluation_point, residual);
700 };
701
702@endcode
703essentially, defines the function from where the residual will be calculated. Notice here that the `compute_residual()` function is called with `evaluation_point` and `residual`.
704
705Similarly, the lines,
706@code{.sh}
707 nonlinear_solver.setup_jacobian =
708 [&](const Vector<double> &current_u) {
709 compute_jacobian(current_u);
710 };
711
712@endcode
713define 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.
714
715The lines
716@code{.sh}
717 nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
718 Vector<double> &dst,
719 const double tolerance) {
720 solve(rhs, dst, tolerance);
721 };
722
723@endcode
724tells the solution needs to be done with the jacobian using the function `solve()`.
725
726Then the line
727@code{.sh}
728 nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
729 Vector<double> &dst,
730 const double tolerance) {
731 solve(rhs, dst, tolerance);
732 };
733
734@endcode
735actually calls the nonlinear solver.
736
737For 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
738@code{.sh}
739 converged_solution = present_solution;
740
741@endcode.
742The 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.
743@code{.sh}
744 present_solution = converged_solution;
745
746@endcode
747
748
749### Results
750
751The 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.
752
753![image](../code-gallery/nonlinear-heat_transfer_with_AD_NOX/doc/Images/contour.png)
754
755Contour plot of the temperature at the final step
756
757![image](../code-gallery/nonlinear-heat_transfer_with_AD_NOX/doc/Images/Temp_evol.png)
758
759Evolution of temperature at a point close to the right edge @f$\approx (0.49, 0.12)@f$
760
761In closing, we give some ideas as to how the residuals and the jacobians are actually calculated in a finite element setting.
762
763
764### Evaluating the function at previous iterations
765
766For 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
767@f[
768 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}
769@f]
770Each term is written in the following manner. We write it for only one term to illustrate the details.
771Consider @f$\alpha \Delta t L_I^{s+1} @f$
772@f[
773\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}
774@f]
775to be specific
776@f[
777 \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}
778@f]
779again, 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
780@f[
781L_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
782@f]
783where @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.
784
785Similar evaluations have to be done to obtain the Jacobian as well. In the current case, the per element Jacobian is given by
786@f[
787 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}
788@f]
789which is written as
790@f[
791 J_{IQ}^e = \int_{\Omega_e}\rho C_p \psi_I \psi_J dx dy +
792 \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
793@f]
794which 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.
795
796
797
798
799## The dealii steps
800To understand this application, @ref step_71 "step-71", @ref step_72 "step-72" and @ref step_77 "step-77" are needed.
801
802
803<a name="ann-mesh/README.md"></a>
804<h1>Annotated version of mesh/README.md</h1>
805The following copyright notice applies to the mesh files in this directory:
806@code{.sh}
807/*
808 * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
809 * Copyright (C) 2024 by Narasimhan Swaminathan
810 *
811 * This file is part of the deal.II code gallery.
812 *
813 * -----------------------------------------------------------------------------
814 */
815@endcode
816
817
818<a name="ann-include/allheaders.h"></a>
819<h1>Annotated version of include/allheaders.h</h1>
820 *
821 *
822 *
823 *
824 * @code
825 *   /* -----------------------------------------------------------------------------
826 *   *
827 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
828 *   * Copyright (C) 2024 by Narasimhan Swaminathan
829 *   *
830 *   * This file is part of the deal.II code gallery.
831 *   *
832 *   * -----------------------------------------------------------------------------
833 *   */
834 *  
835 *   #ifndef __ALLHEADERS_H_INCLUDED__
836 *   #define __ALLHEADERS_H_INCLUDED__
837 * @endcode
838 *
839 * These are some of the header files for this programme. Most of them
840 * are common from previous steps.
841 *
842 * @code
843 *   #include <deal.II/base/quadrature_lib.h>
844 *   #include <deal.II/base/function_lib.h>
845 *   #include <deal.II/lac/generic_linear_algebra.h>
846 *   #include <deal.II/base/function.h>
847 *  
848 *   #include <deal.II/dofs/dof_handler.h>
849 *   #include <deal.II/dofs/dof_tools.h>
850 *  
851 *   #include <deal.II/fe/fe_values.h>
852 *   #include <deal.II/fe/fe.h>
853 *  
854 *   #include <deal.II/grid/tria.h>
855 *   #include <deal.II/grid/grid_generator.h>
856 *  
857 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
858 *   #include <deal.II/lac/full_matrix.h>
859 *   #include <deal.II/lac/precondition.h>
860 *   #include <deal.II/lac/solver_cg.h>
861 *   #include <deal.II/lac/sparse_matrix.h>
862 *   #include <deal.II/lac/vector.h>
863 *   #include <deal.II/grid/grid_tools.h>
864 *  
865 *   #include <deal.II/numerics/data_out.h>
866 *   #include <deal.II/numerics/vector_tools.h>
867 *   #include <deal.II/fe/fe_q.h>
868 *   #include <deal.II/grid/grid_out.h>
869 *   #include <deal.II/grid/grid_in.h>
870 *  
871 *   #include <deal.II/fe/fe_system.h>
872 *  
873 *   #include <deal.II/grid/tria_accessor.h>
874 *   #include <deal.II/grid/tria_iterator.h>
875 *   #include <deal.II/dofs/dof_handler.h>
876 *   #include <deal.II/dofs/dof_accessor.h>
877 *   #include <deal.II/dofs/dof_tools.h>
878 *   #include <deal.II/fe/fe_values.h>
879 *   #include <deal.II/fe/fe.h>
880 *   #include <deal.II/numerics/matrix_tools.h>
881 *   #include <deal.II/numerics/vector_tools.h>
882 *   #include <deal.II/lac/vector.h>
883 *   #include <deal.II/base/quadrature.h>
884 *   #include <deal.II/distributed/tria.h>
885 *   #include <deal.II/distributed/grid_refinement.h>
886 *  
887 *   #include <deal.II/lac/sparse_direct.h>
888 *   #include <deal.II/base/timer.h>
889 *   #include <deal.II/base/utilities.h>
890 *  
891 *   #include <deal.II/base/exceptions.h>
892 *   #include <deal.II/base/geometric_utilities.h>
893 *   #include <deal.II/base/conditional_ostream.h>
894 *   #include <deal.II/base/mpi.h>
895 *  
896 *   #include <deal.II/grid/grid_tools.h>
897 *   #include <deal.II/dofs/dof_renumbering.h>
898 *   #include <deal.II/numerics/solution_transfer.h>
899 *   #include <deal.II/base/index_set.h>
900 *   #include <deal.II/lac/sparsity_tools.h>
901 *   #include <deal.II/fe/fe_values_extractors.h>
902 * @endcode
903 *
904 * Automatic differentiation is carried out with TRILINOS Automatic differentiation scheme/
905 * This is invoked with the following include.
906 *
907 * @code
908 *   #include <deal.II/differentiation/ad.h>
909 * @endcode
910 *
911 * We use Trilinos wrappers based NOX to solve the non-linear equations.
912 *
913 * @code
914 *   #include <deal.II/trilinos/nox.h>
915 *  
916 *   #include <math.h>
917 *   #include <fstream>
918 *   #include <iostream>
919 *  
920 *   using namespace dealii;
921 *   #endif //__ALLHEADERS_H_INCLUDED__
922 * @endcode
923
924
925<a name="ann-include/nonlinear_heat.h"></a>
926<h1>Annotated version of include/nonlinear_heat.h</h1>
927 *
928 *
929 *
930 *
931 * @code
932 *   /* -----------------------------------------------------------------------------
933 *   *
934 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
935 *   * Copyright (C) 2024 by Narasimhan Swaminathan
936 *   *
937 *   * This file is part of the deal.II code gallery.
938 *   *
939 *   * -----------------------------------------------------------------------------
940 *   */
941 *  
942 *   #ifndef __MAIN_ALL_HEADER_H_INCLUDED__
943 *   #define __MAIN_ALL_HEADER_H_INCLUDED__
944 *   #include "allheaders.h"
945 *  
948 *   class nonlinear_heat
949 *   {
950 *   public:
951 *  
954 *   nonlinear_heat();
955 *   ~nonlinear_heat();
956 *  
959 *   void run();
960 *  
963 *   const double delta_t;
964 *  
967 *   const double alpha;
968 *  
971 *   const double tot_time;
972 *  
973 *  
976 *   const double a;
977 *   const double b;
978 *   const double c;
979 *  
982 *   const double Cp;
983 *  
986 *   const double rho;
987 *  
990 *   double time;
991 *   private:
992 *  
996 *   void compute_jacobian(const Vector<double> &evaluation_point);
997 *  
1003 *   void compute_residual(const Vector<double> &evaluation_point, Vector<double> & residual);
1004 *  
1008 *   void setup_system(unsigned int time_step);
1009 *  
1015 *   void solve(const Vector<double> &rhs, Vector<double> & solution, const double tolerance);
1016 *  
1020 *   void output_results(unsigned int prn) const;
1021 *  
1025 *   void set_boundary_conditions(double time);
1026 *  
1028 *   DoFHandler<2> dof_handler;
1029 *   FESystem<2> fe;
1030 *   SparsityPattern sparsity_pattern;
1031 *   SparseMatrix<double> system_matrix; /*Matrix holding the global Jacobian*/
1032 *  
1036 *   std::unique_ptr<SparseDirectUMFPACK> matrix_factorization;
1037 *  
1041 *   Vector<double> converged_solution;/* Converged solution in the previous time step */
1042 *  
1046 *   Vector<double> present_solution;/* Converged solution in the previous time step */
1047 *   };
1048 *  
1049 *  
1053 *   class Initialcondition : public Function<2>
1054 *   {
1055 *   public:
1056 *   Initialcondition(): Function<2>(1)
1057 *   {}
1058 * @endcode
1059 *
1060 * Returns the initial values.
1061 *
1062 * @code
1063 *   virtual double value(const Point<2> &p,
1064 *   const unsigned int component =0) const override;
1065 *   };
1066 *  
1067 *  
1071 *   class Boundary_values_left:public Function<2>
1072 *   {
1073 *   public:
1074 *   Boundary_values_left(): Function<2>(1)
1075 *   {}
1076 *   virtual double value(const Point<2> & p,const unsigned int component = 0) const override;
1077 *   };
1078 *   #endif //__MAIN_ALL_HEADER_H_INCLUDED__
1079 * @endcode
1080
1081
1082<a name="ann-nonlinear_heat.cc"></a>
1083<h1>Annotated version of nonlinear_heat.cc</h1>
1084 *
1085 *
1086 *
1087 *
1088 * @code
1089 *   /* -----------------------------------------------------------------------------
1090 *   *
1091 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1092 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1093 *   *
1094 *   * This file is part of the deal.II code gallery.
1095 *   *
1096 *   * -----------------------------------------------------------------------------
1097 *   */
1098 *  
1099 *   #include "allheaders.h"
1100 *   #include "nonlinear_heat.h"
1101 *   int main()
1102 *   {
1103 *   try
1104 *   {
1105 *   nonlinear_heat nlheat;
1106 *   nlheat.run();
1107 *   }
1108 *   catch (std::exception &exc)
1109 *   {
1110 *   std::cerr << std::endl
1111 *   << std::endl
1112 *   << "----------------------------------------------------"
1113 *   << std::endl;
1114 *   std::cerr << "Exception on processing: " << std::endl
1115 *   << exc.what() << std::endl
1116 *   << "Aborting!" << std::endl
1117 *   << "----------------------------------------------------"
1118 *   << std::endl;
1119 *  
1120 *   return 1;
1121 *   }
1122 *   catch (...)
1123 *   {
1124 *   std::cerr << std::endl
1125 *   << std::endl
1126 *   << "----------------------------------------------------"
1127 *   << std::endl;
1128 *   std::cerr << "Unknown exception!" << std::endl
1129 *   << "Aborting!" << std::endl
1130 *   << "----------------------------------------------------"
1131 *   << std::endl;
1132 *   return 1;
1133 *   }
1134 *  
1135 *   return 0;
1136 *   }
1137 * @endcode
1138
1139
1140<a name="ann-source/boundary_values.cc"></a>
1141<h1>Annotated version of source/boundary_values.cc</h1>
1142 *
1143 *
1144 *
1145 *
1146 * @code
1147 *   /* -----------------------------------------------------------------------------
1148 *   *
1149 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1150 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1151 *   *
1152 *   * This file is part of the deal.II code gallery.
1153 *   *
1154 *   * -----------------------------------------------------------------------------
1155 *   */
1156 *  
1157 *   #include "nonlinear_heat.h"
1158 *  
1164 *   double Boundary_values_left::value(const Point<2> & /*p*/, const unsigned int /*comp*/) const
1165 *   {
1166 *   return 100;
1167 *  
1168 *  
1173 * @endcode
1174 *
1175 * nonlinear_heat nlheat;
1176 * double total_time = nlheat.tot_time;
1177 * return this->get_time() * 100.0/total_time;
1178 *
1179 * @code
1180 *   }
1181 * @endcode
1182
1183
1184<a name="ann-source/compute_jacobian.cc"></a>
1185<h1>Annotated version of source/compute_jacobian.cc</h1>
1186 *
1187 *
1188 *
1189 *
1190 * @code
1191 *   /* -----------------------------------------------------------------------------
1192 *   *
1193 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1194 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1195 *   *
1196 *   * This file is part of the deal.II code gallery.
1197 *   *
1198 *   * -----------------------------------------------------------------------------
1199 *   */
1200 *  
1201 *   #include "allheaders.h"
1202 *   #include "nonlinear_heat.h"
1203 *  
1209 *   void nonlinear_heat::compute_jacobian(const Vector<double> &evaluation_point)
1210 *   {
1211 *   const QGauss<2> quadrature_formula(
1212 *   fe.degree + 1);
1213 *  
1214 *   const QGauss<1> face_quadrature_formula(fe.degree+1);
1215 *  
1216 *   FEValues<2> fe_values(fe,
1217 *   quadrature_formula,
1220 *   update_JxW_values);
1221 *  
1222 *   FEFaceValues<2> fe_face_values(fe,face_quadrature_formula,update_values|update_quadrature_points|
1224 *  
1225 *  
1226 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
1227 *   const unsigned int n_q_points = quadrature_formula.size();
1228 *   const unsigned int n_q_f_points = face_quadrature_formula.size();
1229 *   Vector<double> cell_rhs(dofs_per_cell);
1230 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1232 *   using ADNumberType = typename ADHelper::ad_type;
1233 *   const FEValuesExtractors::Scalar t(0);
1234 *  
1237 *   system_matrix = 0.0;
1238 *   std::vector<types::global_dof_index> local_dof_indices(
1239 *   dofs_per_cell);
1240 *   /*==================================================================*/
1241 *   std::vector<double> consol(n_q_points);
1242 *   std::vector<ADNumberType> old_solution(
1243 *   n_q_points);
1244 *  
1245 *   std::vector<Tensor<1, 2>> consol_grad(
1246 *   n_q_points);
1247 *   std::vector<Tensor<1, 2,ADNumberType>> old_solution_grad(
1248 *   n_q_points);
1249 *  
1250 *   for (const auto &cell: dof_handler.active_cell_iterators())
1251 *   {
1252 *   cell_rhs = 0;
1253 *   cell_matrix = 0;
1254 *   fe_values.reinit(cell);
1255 *   cell->get_dof_indices(local_dof_indices);
1256 *   const unsigned int n_independent_variables = local_dof_indices.size();
1257 *   const unsigned int n_dependent_variables = dofs_per_cell;
1258 *   ADHelper ad_helper(n_independent_variables, n_dependent_variables);
1259 *   ad_helper.register_dof_values(evaluation_point,local_dof_indices);
1260 *   const std::vector<ADNumberType> &dof_values_ad = ad_helper.get_sensitive_dof_values();
1261 *   fe_values[t].get_function_values_from_local_dof_values(dof_values_ad,
1262 *   old_solution);
1263 *   fe_values[t].get_function_gradients_from_local_dof_values(dof_values_ad,
1264 *   old_solution_grad);
1265 *  
1266 *   fe_values[t].get_function_values(converged_solution, consol);
1267 *   fe_values[t].get_function_gradients(converged_solution, consol_grad);
1268 *   std::vector<ADNumberType> residual_ad(n_dependent_variables,
1269 *   ADNumberType(0.0));
1270 *   for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
1271 *   {
1272 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1273 *   {
1274 *  
1275 *   ADNumberType MijTjcurr = Cp*rho*fe_values[t].value(i,q_index)*old_solution[q_index];
1276 *   ADNumberType MijTjprev = Cp*rho*fe_values[t].value(i,q_index)*consol[q_index];
1277 *   ADNumberType k_curr = a + b*old_solution[q_index] + c*std::pow(old_solution[q_index],2);
1278 *   ADNumberType k_prev = a + b*consol[q_index] + c*std::pow(consol[q_index],2);
1279 *   ADNumberType Licurr = alpha * delta_t * (fe_values[t].gradient(i,q_index)*k_curr*old_solution_grad[q_index]);
1280 *   ADNumberType Liprev = (1-alpha) * delta_t * (fe_values[t].gradient(i,q_index)*k_prev*consol_grad[q_index]);
1281 *   residual_ad[i] += (MijTjcurr+Licurr-MijTjprev+Liprev)*fe_values.JxW(q_index);
1282 *   }
1283 *   }
1284 *  
1285 *  
1286 *   for (unsigned int face_number = 0;face_number<GeometryInfo<2>::faces_per_cell; ++face_number)
1287 *   {
1288 *  
1289 *   if (cell->face(face_number)->boundary_id() == 3)
1290 *   {
1291 *   fe_face_values.reinit(cell, face_number);
1292 *   for (unsigned int q_point=0;q_point<n_q_f_points;++q_point)
1293 *   {
1294 *   for (unsigned int i =0;i<dofs_per_cell;++i)
1295 *   residual_ad[i]+= -delta_t*(-10)*fe_face_values[t].value(i,q_point)*fe_face_values.JxW(q_point);
1296 *  
1297 *   }
1298 *   }
1299 *   }
1300 *  
1301 *   ad_helper.register_residual_vector(residual_ad);
1302 *   ad_helper.compute_residual(cell_rhs);
1303 *  
1306 *   ad_helper.compute_linearization(cell_matrix);
1307 *   cell->get_dof_indices(local_dof_indices);
1308 *  
1311 *  
1312 *   for (unsigned int i =0;i < dofs_per_cell; ++i){
1313 *   for(unsigned int j = 0;j < dofs_per_cell;++j){
1314 *   system_matrix.add(local_dof_indices[i],local_dof_indices[j],cell_matrix(i,j));
1315 *   }
1316 *   }
1317 *   }
1318 *  
1322 *   std::map<types::global_dof_index, double> boundary_values;
1324 *   1,
1326 *   boundary_values);
1327 *  
1328 *   Vector<double> dummy_solution(dof_handler.n_dofs());
1329 *   Vector<double> dummy_rhs(dof_handler.n_dofs());
1330 *   MatrixTools::apply_boundary_values(boundary_values,
1331 *   system_matrix,
1332 *   dummy_solution,
1333 *   dummy_rhs);
1334 *  
1335 *   {
1336 *   std::cout << " Factorizing Jacobian matrix" << std::endl;
1337 *   matrix_factorization = std::make_unique<SparseDirectUMFPACK>();
1338 *   matrix_factorization->factorize(system_matrix);
1339 *   }
1340 *   }
1341 * @endcode
1342
1343
1344<a name="ann-source/compute_residual.cc"></a>
1345<h1>Annotated version of source/compute_residual.cc</h1>
1346 *
1347 *
1348 *
1349 *
1350 * @code
1351 *   /* -----------------------------------------------------------------------------
1352 *   *
1353 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1354 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1355 *   *
1356 *   * This file is part of the deal.II code gallery.
1357 *   *
1358 *   * -----------------------------------------------------------------------------
1359 *   */
1360 *  
1361 *   #include "allheaders.h"
1362 *   #include "nonlinear_heat.h"
1363 *  
1369 *   void nonlinear_heat::compute_residual(const Vector<double> & evaluation_point, Vector<double> & residual)
1370 *   {
1371 *  
1375 *   const QGauss<2> quadrature_formula(
1376 *   fe.degree + 1);
1377 *   const QGauss<1> face_quadrature_formula(fe.degree+1); //Define quadrature for integration over faces */
1378 *   FEValues<2> fe_values(fe,
1379 *   quadrature_formula,
1382 *   update_JxW_values);
1383 *  
1384 *   FEFaceValues<2> fe_face_values(fe,face_quadrature_formula,update_values|update_quadrature_points|
1386 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
1387 *   const unsigned int n_q_points = quadrature_formula.size();
1388 *   const unsigned int n_q_f_points = face_quadrature_formula.size();
1389 *  
1393 *   Vector<double> cell_rhs(dofs_per_cell);
1394 *  
1395 *  
1400 *   using ADNumberType = typename ADHelper::ad_type;
1401 *  
1407 *   const FEValuesExtractors::Scalar t(
1408 *   0);
1409 *   std::vector<types::global_dof_index> local_dof_indices(
1410 *   dofs_per_cell);
1411 *  
1414 *  
1415 *   std::vector<double> consol(n_q_points); /* Converged solution at the Gauss points from the previous time step*/
1416 *   std::vector<ADNumberType> old_solution(
1417 *   n_q_points); /* Current solution at the Gauss points at this iteration for the current time*/
1418 *  
1421 *   std::vector<Tensor<1, 2>> consol_grad(
1422 *   n_q_points); /* Converged gradients of the solutions at the Gauss points from the previous time step */
1423 *   std::vector<Tensor<1, 2,ADNumberType>> old_solution_grad(
1424 *   n_q_points);
1425 *  
1426 *   ComponentMask t_mask = fe.component_mask(t);
1427 *  
1430 *   residual = 0.0;
1431 *   for (const auto &cell: dof_handler.active_cell_iterators())
1432 *   {
1433 *   cell_rhs = 0;
1434 *   fe_values.reinit(cell);
1435 *   cell->get_dof_indices(local_dof_indices);
1436 *   const unsigned int n_independent_variables = local_dof_indices.size();
1437 *   const unsigned int n_dependent_variables = dofs_per_cell;
1438 *   ADHelper ad_helper(n_independent_variables, n_dependent_variables);
1439 *  
1443 *   ad_helper.register_dof_values(evaluation_point,local_dof_indices);
1444 *  
1445 *   const std::vector<ADNumberType> &dof_values_ad = ad_helper.get_sensitive_dof_values();
1446 *   fe_values[t].get_function_values_from_local_dof_values(dof_values_ad,
1447 *   old_solution);
1448 *   fe_values[t].get_function_gradients_from_local_dof_values(dof_values_ad,
1449 *   old_solution_grad);
1450 *  
1454 *   fe_values[t].get_function_values(converged_solution, consol);
1455 *   fe_values[t].get_function_gradients(converged_solution, consol_grad);
1456 *  
1459 *   std::vector<ADNumberType> residual_ad(n_dependent_variables,
1460 *   ADNumberType(0.0));
1461 *   for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
1462 *   {
1463 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1464 *   {
1465 *  
1468 *   ADNumberType MijTjcurr = Cp*rho*fe_values[t].value(i,q_index)*old_solution[q_index];
1469 *   ADNumberType MijTjprev = Cp*rho*fe_values[t].value(i,q_index)*consol[q_index];
1470 *   ADNumberType k_curr = a + b*old_solution[q_index] + c*std::pow(old_solution[q_index],2);
1471 *   ADNumberType k_prev = a + b*consol[q_index] + c*std::pow(consol[q_index],2);
1472 *   ADNumberType Licurr = alpha * delta_t * (fe_values[t].gradient(i,q_index)*k_curr*old_solution_grad[q_index]);
1473 *   ADNumberType Liprev = (1-alpha) * delta_t * (fe_values[t].gradient(i,q_index)*k_prev*consol_grad[q_index]);
1474 *   residual_ad[i] += (MijTjcurr+Licurr-MijTjprev+Liprev)*fe_values.JxW(q_index);
1475 *   }
1476 *   }
1477 *  
1480 *  
1481 *   for (unsigned int face_number = 0;face_number<GeometryInfo<2>::faces_per_cell; ++face_number)
1482 *   {
1483 *  
1484 *   if (cell->face(face_number)->boundary_id() == 3)
1485 *   {
1486 *   fe_face_values.reinit(cell, face_number);
1487 *   for (unsigned int q_point=0;q_point<n_q_f_points;++q_point)
1488 *   {
1489 *   for (unsigned int i =0;i<dofs_per_cell;++i)
1490 *  
1493 *   residual_ad[i]+= -delta_t*(-10)*fe_face_values[t].value(i,q_point)*fe_face_values.JxW(q_point);
1494 *  
1495 *   }
1496 *   }
1497 *   }
1498 *  
1501 *   ad_helper.register_residual_vector(residual_ad);
1502 *  
1507 *   ad_helper.compute_residual(cell_rhs);
1508 *   cell->get_dof_indices(local_dof_indices);
1509 *  
1510 *   for (unsigned int i =0;i < dofs_per_cell; ++i)
1511 *   residual(local_dof_indices[i])+= cell_rhs(i);
1512 *  
1513 *   }
1514 *  
1515 *   for(const types::global_dof_index i: DoFTools::extract_boundary_dofs(dof_handler,t_mask,{1}))
1516 *   residual(i) = 0;
1517 *  
1518 *   std::cout << " The Norm is :: = " << residual.l2_norm() << std::endl;
1519 *   }
1520 * @endcode
1521
1522
1523<a name="ann-source/initial_conditions.cc"></a>
1524<h1>Annotated version of source/initial_conditions.cc</h1>
1525 *
1526 *
1527 *
1528 *
1529 * @code
1530 *   /* -----------------------------------------------------------------------------
1531 *   *
1532 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1533 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1534 *   *
1535 *   * This file is part of the deal.II code gallery.
1536 *   *
1537 *   * -----------------------------------------------------------------------------
1538 *   */
1539 *  
1540 *   #include "nonlinear_heat.h"
1541 *  
1542 *  
1548 *   double Initialcondition::value(const Point<2> & /*p*/, const unsigned int /*comp*/) const
1549 *   {
1550 *  
1553 *   return 0.0;
1554 *   }
1555 * @endcode
1556
1557
1558<a name="ann-source/nonlinear_heat_cons_des.cc"></a>
1559<h1>Annotated version of source/nonlinear_heat_cons_des.cc</h1>
1560 *
1561 *
1562 *
1563 *
1564 * @code
1565 *   /* -----------------------------------------------------------------------------
1566 *   *
1567 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1568 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1569 *   *
1570 *   * This file is part of the deal.II code gallery.
1571 *   *
1572 *   * -----------------------------------------------------------------------------
1573 *   */
1574 *  
1575 *   #include "nonlinear_heat.h"
1576 *  
1577 *  
1581 *   nonlinear_heat::nonlinear_heat ()
1582 *   :delta_t(0.1),
1583 *   alpha(0.5),
1584 *   tot_time(5),
1585 *   a(0.3),
1586 *   b(0.003),
1587 *   c(0),
1588 *   Cp(1),
1589 *   rho(1),
1590 *   dof_handler(triangulation),
1591 *   fe(FE_Q<2>(1), 1)
1592 *   {}
1593 *  
1594 *  
1597 *   nonlinear_heat::~nonlinear_heat()
1598 *   {
1599 *   dof_handler.clear();
1600 *   }
1601 *  
1602 *  
1603 * @endcode
1604
1605
1606<a name="ann-source/output_results.cc"></a>
1607<h1>Annotated version of source/output_results.cc</h1>
1608 *
1609 *
1610 *
1611 *
1612 * @code
1613 *   /* -----------------------------------------------------------------------------
1614 *   *
1615 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1616 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1617 *   *
1618 *   * This file is part of the deal.II code gallery.
1619 *   *
1620 *   * -----------------------------------------------------------------------------
1621 *   */
1622 *  
1623 *   #include "nonlinear_heat.h"
1624 *  
1625 *  
1629 *   void nonlinear_heat::output_results(unsigned int prn) const
1630 *   {
1631 *   DataOut<2> data_out;
1632 *   data_out.attach_dof_handler(dof_handler);
1633 *   std::vector<std::string> solution_names;
1634 *   solution_names.emplace_back ("Temperature");
1635 *   data_out.add_data_vector(converged_solution, solution_names);
1636 *   data_out.build_patches();
1637 *   const std::string filename =
1638 *   "output/solution-" + Utilities::int_to_string(prn , 3) + ".vtu";
1639 *   std::ofstream output(filename);
1640 *   data_out.write_vtu(output);
1641 *   }
1642 * @endcode
1643
1644
1645<a name="ann-source/set_boundary_conditions.cc"></a>
1646<h1>Annotated version of source/set_boundary_conditions.cc</h1>
1647 *
1648 *
1649 *
1650 *
1651 * @code
1652 *   /* -----------------------------------------------------------------------------
1653 *   *
1654 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1655 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1656 *   *
1657 *   * This file is part of the deal.II code gallery.
1658 *   *
1659 *   * -----------------------------------------------------------------------------
1660 *   */
1661 *  
1662 *   #include "nonlinear_heat.h"
1663 *   #include "allheaders.h"
1664 *  
1665 *  
1669 *   void nonlinear_heat::set_boundary_conditions(double time)
1670 *   {
1671 *   Boundary_values_left bl_left;
1672 *   bl_left.set_time(time);
1673 *   std::map<types::global_dof_index,double> boundary_values;
1675 *   1,
1676 *   bl_left,
1677 *   boundary_values);
1678 *  
1679 *   for (auto &boundary_value: boundary_values)
1680 *   present_solution(boundary_value.first) = boundary_value.second;
1681 *   }
1682 * @endcode
1683
1684
1685<a name="ann-source/setup_system.cc"></a>
1686<h1>Annotated version of source/setup_system.cc</h1>
1687 *
1688 *
1689 *
1690 *
1691 * @code
1692 *   /* -----------------------------------------------------------------------------
1693 *   *
1694 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1695 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1696 *   *
1697 *   * This file is part of the deal.II code gallery.
1698 *   *
1699 *   * -----------------------------------------------------------------------------
1700 *   */
1701 *  
1702 *   #include "allheaders.h"
1703 *   #include "nonlinear_heat.h"
1704 *  
1705 *  
1709 *   void nonlinear_heat::setup_system(unsigned int time_step)
1710 *   {
1711 *   if (time_step ==0) {
1712 *   dof_handler.distribute_dofs(fe);
1713 *   converged_solution.reinit(dof_handler.n_dofs());
1714 *   present_solution.reinit(dof_handler.n_dofs());
1715 *   }
1716 *  
1717 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
1718 *   DoFTools::make_sparsity_pattern(dof_handler, dsp);
1719 *   sparsity_pattern.copy_from(dsp);
1720 *   system_matrix.reinit(sparsity_pattern);
1721 *   matrix_factorization.reset();
1722 *   }
1723 * @endcode
1724
1725
1726<a name="ann-source/solve_and_run.cc"></a>
1727<h1>Annotated version of source/solve_and_run.cc</h1>
1728 *
1729 *
1730 *
1731 *
1732 * @code
1733 *   /* -----------------------------------------------------------------------------
1734 *   *
1735 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1736 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1737 *   *
1738 *   * This file is part of the deal.II code gallery.
1739 *   *
1740 *   * -----------------------------------------------------------------------------
1741 *   */
1742 *  
1743 *   #include "nonlinear_heat.h"
1744 *  
1745 *  
1750 *   void nonlinear_heat::solve(const Vector<double> & rhs, Vector<double> & solution, const double /*tolerance*/)
1751 *   {
1752 *   std::cout << " Solving linear system" << std::endl;
1753 *   matrix_factorization->vmult(solution, rhs);
1754 *   }
1755 *  
1756 *   void nonlinear_heat::run()
1757 *   {
1758 *   GridIn<2> gridin;
1760 *   std::ifstream f("mesh/mesh.msh");
1761 *   gridin.read_msh(f);
1763 *  
1764 *   double time = 0;
1765 *   unsigned int timestep_number = 0;
1766 *   unsigned int prn =0;
1767 *  
1768 *  
1769 *   while (time <=tot_time)
1770 *   {
1771 *   if(time ==0)
1772 *   {
1773 *   setup_system(timestep_number);
1774 *  
1775 *   VectorTools::interpolate(dof_handler,
1776 *   Initialcondition(),
1777 *   present_solution);
1778 *  
1779 *   VectorTools::interpolate(dof_handler,
1780 *   Initialcondition(),
1781 *   converged_solution);
1782 *   }
1783 *   else
1784 *   {
1785 *  
1788 *   present_solution = converged_solution;
1789 *   }
1790 *   std::cout<<">>>>> Time now is: "<<time <<std::endl;
1791 *   std::cout<<">>>>> Time step is:"<<timestep_number<<std::endl;
1792 *   set_boundary_conditions(time);
1793 *   {
1794 *  
1795 *   const double target_tolerance = 1e-3;
1796 *  
1800 *   typename TrilinosWrappers::NOXSolver<Vector<double>>::AdditionalData additional_data;
1801 *   additional_data.abs_tol = target_tolerance;
1802 *   additional_data.max_iter = 100;
1803 *   TrilinosWrappers::NOXSolver<Vector<double>> nonlinear_solver(
1804 *   additional_data);
1805 *  
1806 *  
1809 *   nonlinear_solver.residual =
1810 *   [&](const Vector<double> &evaluation_point,
1811 *   Vector<double> &residual) {
1812 *   compute_residual(evaluation_point, residual);
1813 *   };
1814 *  
1815 *  
1818 *  
1819 *   nonlinear_solver.setup_jacobian =
1820 *   [&](const Vector<double> &current_u) {
1821 *   compute_jacobian(current_u);
1822 *   };
1823 *  
1826 *  
1827 *   nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
1828 *   Vector<double> &dst,
1829 *   const double tolerance) {
1830 *   solve(rhs, dst, tolerance);
1831 *   };
1832 *  
1836 *   nonlinear_solver.solve(present_solution);
1837 *   }
1838 *  
1841 *   converged_solution = present_solution;
1842 *  
1843 *   if(timestep_number % 1 == 0) {
1844 *  
1845 *   output_results(prn);
1846 *   prn = prn +1;
1847 *   }
1848 *   timestep_number++;
1849 *   time=time+delta_t;
1850 *  
1851 *   }
1852 *   }
1853 * @endcode
1854
1855
1856*/
std::vector< bool > component_mask
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
Definition fe_q.h:554
void attach_triangulation(Triangulation< dim, spacedim > &tria)
Definition grid_in.cc:153
Definition point.h:111
void refine_global(const unsigned int times=1)
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
__global__ void set(Number *val, const Number s, const size_type N)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
const Event initial
Definition event.cc:64
Expression differentiate(const Expression &f, const Expression &x)
IndexSet extract_boundary_dofs(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask={}, const std::set< types::boundary_id > &boundary_ids={})
Definition dof_tools.cc:614
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471
void apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
std::string get_time()
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={})
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
int(& functions)(const void *v1, const void *v2)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation