deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
\(\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
245#include "allheaders.h"
246#include "nonlinear_heat.h"
253void nonlinear_heat::compute_residual(const Vector<double> & evaluation_point, Vector<double> & residual)
254{
260 fe.degree + 1);
261 const QGauss<1> face_quadrature_formula(fe.degree+1); //Define quadrature for integration over faces */
262 FEValues<2> fe_values(fe,
270 const unsigned int dofs_per_cell = fe.dofs_per_cell;
271 const unsigned int n_q_points = quadrature_formula.size();
272 const unsigned int n_q_f_points = face_quadrature_formula.size();
277 Vector<double> cell_rhs(dofs_per_cell);
279@endcode
280should 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.
281
283
284@code{.sh}
286 using ADNumberType = typename ADHelper::ad_type;
287
288@endcode
289
291@code{.sh}
293 0);
294 std::vector<types::global_dof_index> local_dof_indices(
295 dofs_per_cell);
300 std::vector<double> consol(n_q_points); /* Converged solution at the Gauss points from the previous time step*/
301
302@endcode
303We 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$).
304
305
307crucial. 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
308@f[
310 \begin{bmatrix}
316@f]
318@code{.sh}
319 std::vector<ADNumberType> old_solution(
320 n_q_points); /* Current solution at the Gauss points at this iteration for the current time*/
324 std::vector<Tensor<1, 2>> consol_grad(
325 n_q_points); /* Converged gradients of the solutions at the Gauss points from the previous time step */
326 std::vector<Tensor<1, 2,ADNumberType>> old_solution_grad(
327 n_q_points);
328
329
330@endcode
331
332The lines
333@code{.sh}
334 ComponentMask t_mask = fe.component_mask(t);
338 residual = 0.0;
339
340@endcode
341are straightforward to understand. The vector `residual` will hold the *numerical* value of the residual and is hence initialized to be zero.
342
343The lines
344@code{.sh}
345 for (const auto &cell: dof_handler.active_cell_iterators())
346 {
347 cell_rhs = 0;
348 fe_values.reinit(cell);
349 cell->get_dof_indices(local_dof_indices);
350 const unsigned int n_independent_variables = local_dof_indices.size();
351 const unsigned int n_dependent_variables = dofs_per_cell;
352
353@endcode
355@code{.sh}
356 ADHelper ad_helper(n_independent_variables, n_dependent_variables);
357
358@endcode
359then, 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.
360@code{.sh}
361 ad_helper.register_dof_values(evaluation_point,local_dof_indices);
362
363@endcode
364
365Now the lines
366@code{.sh}
367 const std::vector<ADNumberType> &dof_values_ad = ad_helper.get_sensitive_dof_values();
368 fe_values[t].get_function_values_from_local_dof_values(dof_values_ad,
370 fe_values[t].get_function_gradients_from_local_dof_values(dof_values_ad,
372
373@endcode
375
376In the lines
377@code{.sh}
378 fe_values[t].get_function_values(converged_solution, consol);
379 fe_values[t].get_function_gradients(converged_solution, consol_grad);
380
381@endcode
383
384Then in the lines
385@code{.sh}
386 for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
387 {
388 for (unsigned int i = 0; i < dofs_per_cell; ++i)
389 {
393 ADNumberType MijTjcurr = Cp*rho*fe_values[t].value(i,q_index)*old_solution[q_index];
394 ADNumberType MijTjprev = Cp*rho*fe_values[t].value(i,q_index)*consol[q_index];
395 ADNumberType k_curr = a + b*old_solution[q_index] + c*std::pow(old_solution[q_index],2);
396 ADNumberType k_prev = a + b*consol[q_index] + c*std::pow(consol[q_index],2);
397 ADNumberType Licurr = alpha * delta_t * (fe_values[t].gradient(i,q_index)*k_curr*old_solution_grad[q_index]);
398 ADNumberType Liprev = (1-alpha) * delta_t * (fe_values[t].gradient(i,q_index)*k_prev*consol_grad[q_index]);
399 residual_ad[i] += (MijTjcurr+Licurr-MijTjprev+Liprev)*fe_values.JxW(q_index);
400 }
401 }
432 residual_ad[i]+= -delta_t*(-10)*fe_face_values[t].value(i,q_point)*fe_face_values.JxW(q_point);
433
434 }
435 }
436 }
437
438@endcode
439
440Now
441@code{.sh}
442 ad_helper.register_residual_vector(residual_ad);
443
444@endcode
445actually tells the `ad_helper` that this is the residual it should use, and the line
446@code{.sh}
447 ad_helper.compute_residual(cell_rhs);
448
449@endcode
451@code{.sh}
452 cell->get_dof_indices(local_dof_indices);
453
454 for (unsigned int i =0;i < dofs_per_cell; ++i)
455 residual(local_dof_indices[i])+= cell_rhs(i);
456
457 }
458
459@endcode
460are used to calculate the global `residual` from the local `cell_rhs`. Now `residual` is numerical and not of the special number type. Finally, lines
461@code{.sh}
462 for(const types::global_dof_index i: DoFTools::extract_boundary_dofs(dof_handler,t_mask,{1}))
463 residual(i) = 0;
464
465 std::cout << " The Norm is :: = " << residual.l2_norm() << std::endl;
466}
467
468@endcode
470
471
472#### The compute_jacobian() function
473
474The `compute_jacobian()` only takes in the `evaluation_point` as an input and is identical to the
475`compute_residual()` except for some minor modification. In particular the line
476@code{.sh}
477 ad_helper.compute_linearization(cell_matrix);
478
479@endcode
481
482
483### General ideas involved in solving coupled nonlinear equations using Newton Raphson's technique
485@f{eqnarray}{
486 f_i^1(p_j^{{s+1}}, T_j^{{s+1}})=0 \\
487 f_i^2(p_j^{{s+1}}, T_j^{{s+1}})=0
488@f}
490
491We 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
492@f{eqnarray}{
493 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 \\
494 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
495@f}
496Also, 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.
497@f[
498 J_{pq}=-
499 \begin{bmatrix}
500 \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}\\
501
502 \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}\\
503
504 \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}\\
505
506 \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}\\
507 %================================================================
508 \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}\\
509
510 \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}\\
511
512 \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}\\
513
514 \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}
515 \end{bmatrix}
516@f]
518@f[
519 \begin{bmatrix}
520 p_1^{s+1}\\
521 p_2^{s+1}\\
522 p_3^{s+1}\\
523 p_4^{s+1}\\
524 T_1^{s+1}\\
525 T_2^{s+1}\\
526 T_3^{s+1}\\
527 T_4^{s+1}
528 \end{bmatrix}^{k+1}
529 =
530 \begin{bmatrix}
531 p_1^{s+1}\\
532 p_2^{s+1}\\
533 p_3^{s+1}\\
534 p_4^{s+1}\\
535 T_1^{s+1}\\
536 T_2^{s+1}\\
537 T_3^{s+1}\\
538 T_4^{s+1}
539 \end{bmatrix}^{k}-
540 \begin{bmatrix}
541 \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}\\
542
543 \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}\\
544
545 \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}\\
546
547 \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}\\
548 %================================================================
549 \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}\\
550
551 \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}\\
552
553 \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}\\
554
555 \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}
556 \end{bmatrix}^{-1}
557 \begin{bmatrix}
558 f_1^{1}\\
559 f_2^{1}\\
560 f_3^{1}\\
561 f_4^{1}\\
562 f_1^{2}\\
563 f_2^{2}\\
564 f_3^{2}\\
565 f_4^{2}
566 \end{bmatrix}^{k}
567@f]
569@f[
570 \begin{bmatrix}
571 \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}\\
572
573 \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}\\
574
575 \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}\\
576
577 \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}\\
578 %================================================================
579 \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}\\
580
581 \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}\\
582
583 \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}\\
584
585 \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}
586 \end{bmatrix}
587 \begin{bmatrix}
588 \Delta p_1^{s+1}\\
589 \Delta p_2^{s+1}\\
590 \Delta p_3^{s+1}\\
591 \Delta p_4^{s+1}\\
592 \Delta T_1^{s+1}\\
593 \Delta T_2^{s+1}\\
594 \Delta T_3^{s+1}\\
595 \Delta T_4^{s+1}
596 \end{bmatrix}=
597 \begin{bmatrix}
598 f_1^{1}\\
599 f_2^{1}\\
600 f_3^{1}\\
601 f_4^{1}\\
602 f_1^{2}\\
603 f_2^{2}\\
604 f_3^{2}\\
605 f_4^{2}
606 \end{bmatrix}^{k}
607@f]
609@f[
610 J_{pq}\{\Delta solution\}^{k+1}=\{function\}^{k}
611@f]
612Then add this to the value of
613@f[
614 \begin{bmatrix}
615 p_1^{s+1}\\
616 p_2^{s+1}\\
617 p_3^{s+1}\\
618 p_4^{s+1}\\
619 T_1^{s+1}\\
620 T_2^{s+1}\\
621 T_3^{s+1}\\
622 T_4^{s+1}
623 \end{bmatrix}^{k}
624@f]
626@f[
627 ||(function)^k||<\delta^0 ||function^0||
628@f]
629or
630@f[
631 ||(function)^k||<tol
632@f]
634
635
636#### The run() function
637
638The function which implements the nonlinear solver is the `run()` function for every time step. In particular the lines
639@code{.sh}
640 typename TrilinosWrappers::NOXSolver<Vector<double>>::AdditionalData additional_data;
641 additional_data.abs_tol = target_tolerance;
642 additional_data.max_iter = 100;
644 additional_data);
645
649 nonlinear_solver.residual =
651 Vector<double> &residual) {
652 compute_residual(evaluation_point, residual);
653 };
654
659 nonlinear_solver.setup_jacobian =
660 [&](const Vector<double> &current_u) {
661 compute_jacobian(current_u);
662 };
667 nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
668 Vector<double> &dst,
669 const double tolerance) {
670 solve(rhs, dst, tolerance);
671 };
677 }
678
679@endcode
680Here, 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.
681
682The lines
683@code{.sh}
684 typename TrilinosWrappers::NOXSolver<Vector<double>>::AdditionalData additional_data;
685 additional_data.abs_tol = target_tolerance;
686 additional_data.max_iter = 100;
688 additional_data);
689
690@endcode
691setup the solver and give any basic data needed, such as the tolerance to converge or the maximum number of iterations it can try.
692
693The lines
694@code{.sh}
695 nonlinear_solver.residual =
697 Vector<double> &residual) {
698 compute_residual(evaluation_point, residual);
699 };
700
701@endcode
702essentially, defines the function from where the residual will be calculated. Notice here that the `compute_residual()` function is called with `evaluation_point` and `residual`.
703
704Similarly, the lines,
705@code{.sh}
706 nonlinear_solver.setup_jacobian =
707 [&](const Vector<double> &current_u) {
708 compute_jacobian(current_u);
709 };
710
711@endcode
712define 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.
713
714The lines
715@code{.sh}
716 nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
717 Vector<double> &dst,
718 const double tolerance) {
719 solve(rhs, dst, tolerance);
720 };
721
722@endcode
723tells the solution needs to be done with the jacobian using the function `solve()`.
724
725Then the line
726@code{.sh}
727 nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
728 Vector<double> &dst,
729 const double tolerance) {
730 solve(rhs, dst, tolerance);
731 };
732
733@endcode
735
736For 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
737@code{.sh}
739
740@endcode.
742@code{.sh}
744
745@endcode
746
747
748### Results
749
750The 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.
751
753
754Contour plot of the temperature at the final step
755
757
758Evolution of temperature at a point close to the right edge @f$\approx (0.49, 0.12)@f$
759
760In closing, we give some ideas as to how the residuals and the jacobians are actually calculated in a finite element setting.
761
762
763### Evaluating the function at previous iterations
764
766@f[
767 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}
768@f]
771@f[
772\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}
773@f]
775@f[
776 \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}
777@f]
779@f[
780L_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
781@f]
782where @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.
783
785@f[
786 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}
787@f]
789@f[
790 J_{IQ}^e = \int_{\Omega_e}\rho C_p \psi_I \psi_J dx dy +
791 \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
792@f]
794
795
796
797
798## The dealii steps
799To understand this application, @ref step_71 "step-71", @ref step_72 "step-72" and @ref step_77 "step-77" are needed.
800
801
802<a name="ann-mesh/README.md"></a>
805@code{.sh}
806/*
807 * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
808 * Copyright (C) 2024 by Narasimhan Swaminathan
809 *
810 * This file is part of the deal.II code gallery.
811 *
812 * -----------------------------------------------------------------------------
813 */
814@endcode
815
816
817<a name="ann-include/allheaders.h"></a>
819 *
820 *
821 *
822 *
823 * @code
824 *   /* -----------------------------------------------------------------------------
825 *   *
826 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
827 *   * Copyright (C) 2024 by Narasimhan Swaminathan
828 *   *
829 *   * This file is part of the deal.II code gallery.
830 *   *
831 *   * -----------------------------------------------------------------------------
832 *   */
833 *  
836 * @endcode
837 *
839 * are common from previous steps.
840 *
841 * @code
842 *   #include <deal.II/base/quadrature_lib.h>
843 *   #include <deal.II/base/function_lib.h>
845 *   #include <deal.II/base/function.h>
846 *  
847 *   #include <deal.II/dofs/dof_handler.h>
848 *   #include <deal.II/dofs/dof_tools.h>
849 *  
850 *   #include <deal.II/fe/fe_values.h>
851 *   #include <deal.II/fe/fe.h>
852 *  
853 *   #include <deal.II/grid/tria.h>
854 *   #include <deal.II/grid/grid_generator.h>
855 *  
857 *   #include <deal.II/lac/full_matrix.h>
858 *   #include <deal.II/lac/precondition.h>
859 *   #include <deal.II/lac/solver_cg.h>
861 *   #include <deal.II/lac/vector.h>
862 *   #include <deal.II/grid/grid_tools.h>
863 *  
864 *   #include <deal.II/numerics/data_out.h>
865 *   #include <deal.II/numerics/vector_tools.h>
866 *   #include <deal.II/fe/fe_q.h>
867 *   #include <deal.II/grid/grid_out.h>
868 *   #include <deal.II/grid/grid_in.h>
869 *  
870 *   #include <deal.II/fe/fe_system.h>
871 *  
872 *   #include <deal.II/grid/tria_accessor.h>
873 *   #include <deal.II/grid/tria_iterator.h>
874 *   #include <deal.II/dofs/dof_handler.h>
875 *   #include <deal.II/dofs/dof_accessor.h>
876 *   #include <deal.II/dofs/dof_tools.h>
877 *   #include <deal.II/fe/fe_values.h>
878 *   #include <deal.II/fe/fe.h>
879 *   #include <deal.II/numerics/matrix_tools.h>
880 *   #include <deal.II/numerics/vector_tools.h>
881 *   #include <deal.II/lac/vector.h>
882 *   #include <deal.II/base/quadrature.h>
883 *   #include <deal.II/distributed/tria.h>
884 *   #include <deal.II/distributed/grid_refinement.h>
885 *  
887 *   #include <deal.II/base/timer.h>
888 *   #include <deal.II/base/utilities.h>
889 *  
890 *   #include <deal.II/base/exceptions.h>
891 *   #include <deal.II/base/geometric_utilities.h>
892 *   #include <deal.II/base/conditional_ostream.h>
893 *   #include <deal.II/base/mpi.h>
894 *  
895 *   #include <deal.II/grid/grid_tools.h>
896 *   #include <deal.II/dofs/dof_renumbering.h>
897 *   #include <deal.II/numerics/solution_transfer.h>
898 *   #include <deal.II/base/index_set.h>
901 * @endcode
902 *
905 *
906 * @code
908 * @endcode
909 *
910 * We use Trilinos wrappers based NOX to solve the non-linear equations.
911 *
912 * @code
913 *   #include <deal.II/trilinos/nox.h>
914 *  
915 *   #include <math.h>
918 *  
919 *   using namespace dealii;
920 *   #endif //__ALLHEADERS_H_INCLUDED__
921 * @endcode
922
923
924<a name="ann-include/nonlinear_heat.h"></a>
926 *
927 *
928 *
929 *
930 * @code
931 *   /* -----------------------------------------------------------------------------
932 *   *
933 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
934 *   * Copyright (C) 2024 by Narasimhan Swaminathan
935 *   *
936 *   * This file is part of the deal.II code gallery.
937 *   *
938 *   * -----------------------------------------------------------------------------
939 *   */
940 *  
943 *   #include "allheaders.h"
944 *  
947 *   class nonlinear_heat
948 *   {
949 *   public:
950 *  
953 *   nonlinear_heat();
955 *  
958 *   void run();
959 *  
962 *   const double delta_t;
963 *  
966 *   const double alpha;
967 *  
970 *   const double tot_time;
971 *  
972 *  
975 *   const double a;
976 *   const double b;
977 *   const double c;
978 *  
981 *   const double Cp;
982 *  
985 *   const double rho;
986 *  
989 *   double time;
990 *   private:
991 *  
995 *   void compute_jacobian(const Vector<double> &evaluation_point);
996 *  
1002 *   void compute_residual(const Vector<double> &evaluation_point, Vector<double> & residual);
1003 *  
1007 *   void setup_system(unsigned int time_step);
1008 *  
1014 *   void solve(const Vector<double> &rhs, Vector<double> & solution, const double tolerance);
1015 *  
1019 *   void output_results(unsigned int prn) const;
1020 *  
1024 *   void set_boundary_conditions(double time);
1025 *  
1027 *   DoFHandler<2> dof_handler;
1028 *   FESystem<2> fe;
1029 *   SparsityPattern sparsity_pattern;
1030 *   SparseMatrix<double> system_matrix; /*Matrix holding the global Jacobian*/
1031 *  
1035 *   std::unique_ptr<SparseDirectUMFPACK> matrix_factorization;
1036 *  
1040 *   Vector<double> converged_solution;/* Converged solution in the previous time step */
1041 *  
1045 *   Vector<double> present_solution;/* Converged solution in the previous time step */
1046 *   };
1047 *  
1048 *  
1052 *   class Initialcondition : public Function<2>
1053 *   {
1054 *   public:
1055 *   Initialcondition(): Function<2>(1)
1056 *   {}
1057 * @endcode
1058 *
1060 *
1061 * @code
1062 *   virtual double value(const Point<2> &p,
1063 *   const unsigned int component =0) const override;
1064 *   };
1065 *  
1066 *  
1070 *   class Boundary_values_left:public Function<2>
1071 *   {
1072 *   public:
1074 *   {}
1075 *   virtual double value(const Point<2> & p,const unsigned int component = 0) const override;
1076 *   };
1077 *   #endif //__MAIN_ALL_HEADER_H_INCLUDED__
1078 * @endcode
1079
1080
1081<a name="ann-nonlinear_heat.cc"></a>
1083 *
1084 *
1085 *
1086 *
1087 * @code
1088 *   /* -----------------------------------------------------------------------------
1089 *   *
1090 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1091 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1092 *   *
1093 *   * This file is part of the deal.II code gallery.
1094 *   *
1095 *   * -----------------------------------------------------------------------------
1096 *   */
1097 *  
1098 *   #include "allheaders.h"
1099 *   #include "nonlinear_heat.h"
1100 *   int main()
1101 *   {
1102 *   try
1103 *   {
1105 *   nlheat.run();
1106 *   }
1107 *   catch (std::exception &exc)
1108 *   {
1109 *   std::cerr << std::endl
1110 *   << std::endl
1111 *   << "----------------------------------------------------"
1112 *   << std::endl;
1113 *   std::cerr << "Exception on processing: " << std::endl
1114 *   << exc.what() << std::endl
1115 *   << "Aborting!" << std::endl
1116 *   << "----------------------------------------------------"
1117 *   << std::endl;
1118 *  
1119 *   return 1;
1120 *   }
1121 *   catch (...)
1122 *   {
1123 *   std::cerr << std::endl
1124 *   << std::endl
1125 *   << "----------------------------------------------------"
1126 *   << std::endl;
1127 *   std::cerr << "Unknown exception!" << std::endl
1128 *   << "Aborting!" << std::endl
1129 *   << "----------------------------------------------------"
1130 *   << std::endl;
1131 *   return 1;
1132 *   }
1133 *  
1134 *   return 0;
1135 *   }
1136 * @endcode
1137
1138
1139<a name="ann-source/boundary_values.cc"></a>
1141 *
1142 *
1143 *
1144 *
1145 * @code
1146 *   /* -----------------------------------------------------------------------------
1147 *   *
1148 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1149 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1150 *   *
1151 *   * This file is part of the deal.II code gallery.
1152 *   *
1153 *   * -----------------------------------------------------------------------------
1154 *   */
1155 *  
1156 *   #include "nonlinear_heat.h"
1157 *  
1163 *   double Boundary_values_left::value(const Point<2> & /*p*/, const unsigned int /*comp*/) const
1164 *   {
1165 *   return 100;
1166 *  
1167 *  
1172 * @endcode
1173 *
1175 * double total_time = nlheat.tot_time;
1176 * return this->get_time() * 100.0/total_time;
1177 *
1178 * @code
1179 *   }
1180 * @endcode
1181
1182
1183<a name="ann-source/compute_jacobian.cc"></a>
1184<h1>Annotated version of source/compute_jacobian.cc</h1>
1185 *
1186 *
1187 *
1188 *
1189 * @code
1190 *   /* -----------------------------------------------------------------------------
1191 *   *
1192 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1193 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1194 *   *
1195 *   * This file is part of the deal.II code gallery.
1196 *   *
1197 *   * -----------------------------------------------------------------------------
1198 *   */
1199 *  
1200 *   #include "allheaders.h"
1201 *   #include "nonlinear_heat.h"
1202 *  
1208 *   void nonlinear_heat::compute_jacobian(const Vector<double> &evaluation_point)
1209 *   {
1211 *   fe.degree + 1);
1212 *  
1213 *   const QGauss<1> face_quadrature_formula(fe.degree+1);
1214 *  
1215 *   FEValues<2> fe_values(fe,
1220 *  
1223 *  
1224 *  
1225 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
1226 *   const unsigned int n_q_points = quadrature_formula.size();
1227 *   const unsigned int n_q_f_points = face_quadrature_formula.size();
1228 *   Vector<double> cell_rhs(dofs_per_cell);
1229 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1231 *   using ADNumberType = typename ADHelper::ad_type;
1232 *   const FEValuesExtractors::Scalar t(0);
1233 *  
1236 *   system_matrix = 0.0;
1237 *   std::vector<types::global_dof_index> local_dof_indices(
1238 *   dofs_per_cell);
1239 *   /*==================================================================*/
1240 *   std::vector<double> consol(n_q_points);
1241 *   std::vector<ADNumberType> old_solution(
1242 *   n_q_points);
1243 *  
1244 *   std::vector<Tensor<1, 2>> consol_grad(
1245 *   n_q_points);
1246 *   std::vector<Tensor<1, 2,ADNumberType>> old_solution_grad(
1247 *   n_q_points);
1248 *  
1249 *   for (const auto &cell: dof_handler.active_cell_iterators())
1250 *   {
1251 *   cell_rhs = 0;
1252 *   cell_matrix = 0;
1253 *   fe_values.reinit(cell);
1254 *   cell->get_dof_indices(local_dof_indices);
1255 *   const unsigned int n_independent_variables = local_dof_indices.size();
1256 *   const unsigned int n_dependent_variables = dofs_per_cell;
1257 *   ADHelper ad_helper(n_independent_variables, n_dependent_variables);
1258 *   ad_helper.register_dof_values(evaluation_point,local_dof_indices);
1259 *   const std::vector<ADNumberType> &dof_values_ad = ad_helper.get_sensitive_dof_values();
1260 *   fe_values[t].get_function_values_from_local_dof_values(dof_values_ad,
1261 *   old_solution);
1262 *   fe_values[t].get_function_gradients_from_local_dof_values(dof_values_ad,
1264 *  
1265 *   fe_values[t].get_function_values(converged_solution, consol);
1266 *   fe_values[t].get_function_gradients(converged_solution, consol_grad);
1267 *   std::vector<ADNumberType> residual_ad(n_dependent_variables,
1268 *   ADNumberType(0.0));
1269 *   for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
1270 *   {
1271 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1272 *   {
1273 *  
1274 *   ADNumberType MijTjcurr = Cp*rho*fe_values[t].value(i,q_index)*old_solution[q_index];
1275 *   ADNumberType MijTjprev = Cp*rho*fe_values[t].value(i,q_index)*consol[q_index];
1276 *   ADNumberType k_curr = a + b*old_solution[q_index] + c*std::pow(old_solution[q_index],2);
1277 *   ADNumberType k_prev = a + b*consol[q_index] + c*std::pow(consol[q_index],2);
1278 *   ADNumberType Licurr = alpha * delta_t * (fe_values[t].gradient(i,q_index)*k_curr*old_solution_grad[q_index]);
1279 *   ADNumberType Liprev = (1-alpha) * delta_t * (fe_values[t].gradient(i,q_index)*k_prev*consol_grad[q_index]);
1280 *   residual_ad[i] += (MijTjcurr+Licurr-MijTjprev+Liprev)*fe_values.JxW(q_index);
1281 *   }
1282 *   }
1283 *  
1284 *  
1285 *   for (unsigned int face_number = 0;face_number<GeometryInfo<2>::faces_per_cell; ++face_number)
1286 *   {
1287 *  
1288 *   if (cell->face(face_number)->boundary_id() == 3)
1289 *   {
1290 *   fe_face_values.reinit(cell, face_number);
1291 *   for (unsigned int q_point=0;q_point<n_q_f_points;++q_point)
1292 *   {
1293 *   for (unsigned int i =0;i<dofs_per_cell;++i)
1294 *   residual_ad[i]+= -delta_t*(-10)*fe_face_values[t].value(i,q_point)*fe_face_values.JxW(q_point);
1295 *  
1296 *   }
1297 *   }
1298 *   }
1299 *  
1300 *   ad_helper.register_residual_vector(residual_ad);
1301 *   ad_helper.compute_residual(cell_rhs);
1302 *  
1305 *   ad_helper.compute_linearization(cell_matrix);
1306 *   cell->get_dof_indices(local_dof_indices);
1307 *  
1310 *  
1311 *   for (unsigned int i =0;i < dofs_per_cell; ++i){
1312 *   for(unsigned int j = 0;j < dofs_per_cell;++j){
1313 *   system_matrix.add(local_dof_indices[i],local_dof_indices[j],cell_matrix(i,j));
1314 *   }
1315 *   }
1316 *   }
1317 *  
1321 *   std::map<types::global_dof_index, double> boundary_values;
1323 *   1,
1326 *  
1327 *   Vector<double> dummy_solution(dof_handler.n_dofs());
1328 *   Vector<double> dummy_rhs(dof_handler.n_dofs());
1330 *   system_matrix,
1332 *   dummy_rhs);
1333 *  
1334 *   {
1335 *   std::cout << " Factorizing Jacobian matrix" << std::endl;
1336 *   matrix_factorization = std::make_unique<SparseDirectUMFPACK>();
1337 *   matrix_factorization->factorize(system_matrix);
1338 *   }
1339 *   }
1340 * @endcode
1341
1342
1343<a name="ann-source/compute_residual.cc"></a>
1344<h1>Annotated version of source/compute_residual.cc</h1>
1345 *
1346 *
1347 *
1348 *
1349 * @code
1350 *   /* -----------------------------------------------------------------------------
1351 *   *
1352 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1353 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1354 *   *
1355 *   * This file is part of the deal.II code gallery.
1356 *   *
1357 *   * -----------------------------------------------------------------------------
1358 *   */
1359 *  
1360 *   #include "allheaders.h"
1361 *   #include "nonlinear_heat.h"
1362 *  
1368 *   void nonlinear_heat::compute_residual(const Vector<double> & evaluation_point, Vector<double> & residual)
1369 *   {
1370 *  
1375 *   fe.degree + 1);
1376 *   const QGauss<1> face_quadrature_formula(fe.degree+1); //Define quadrature for integration over faces */
1377 *   FEValues<2> fe_values(fe,
1382 *  
1385 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
1386 *   const unsigned int n_q_points = quadrature_formula.size();
1387 *   const unsigned int n_q_f_points = face_quadrature_formula.size();
1388 *  
1392 *   Vector<double> cell_rhs(dofs_per_cell);
1393 *  
1394 *  
1399 *   using ADNumberType = typename ADHelper::ad_type;
1400 *  
1407 *   0);
1408 *   std::vector<types::global_dof_index> local_dof_indices(
1409 *   dofs_per_cell);
1410 *  
1413 *  
1414 *   std::vector<double> consol(n_q_points); /* Converged solution at the Gauss points from the previous time step*/
1415 *   std::vector<ADNumberType> old_solution(
1416 *   n_q_points); /* Current solution at the Gauss points at this iteration for the current time*/
1417 *  
1420 *   std::vector<Tensor<1, 2>> consol_grad(
1421 *   n_q_points); /* Converged gradients of the solutions at the Gauss points from the previous time step */
1422 *   std::vector<Tensor<1, 2,ADNumberType>> old_solution_grad(
1423 *   n_q_points);
1424 *  
1425 *   ComponentMask t_mask = fe.component_mask(t);
1426 *  
1429 *   residual = 0.0;
1430 *   for (const auto &cell: dof_handler.active_cell_iterators())
1431 *   {
1432 *   cell_rhs = 0;
1433 *   fe_values.reinit(cell);
1434 *   cell->get_dof_indices(local_dof_indices);
1435 *   const unsigned int n_independent_variables = local_dof_indices.size();
1436 *   const unsigned int n_dependent_variables = dofs_per_cell;
1437 *   ADHelper ad_helper(n_independent_variables, n_dependent_variables);
1438 *  
1442 *   ad_helper.register_dof_values(evaluation_point,local_dof_indices);
1443 *  
1444 *   const std::vector<ADNumberType> &dof_values_ad = ad_helper.get_sensitive_dof_values();
1445 *   fe_values[t].get_function_values_from_local_dof_values(dof_values_ad,
1446 *   old_solution);
1447 *   fe_values[t].get_function_gradients_from_local_dof_values(dof_values_ad,
1449 *  
1453 *   fe_values[t].get_function_values(converged_solution, consol);
1454 *   fe_values[t].get_function_gradients(converged_solution, consol_grad);
1455 *  
1458 *   std::vector<ADNumberType> residual_ad(n_dependent_variables,
1459 *   ADNumberType(0.0));
1460 *   for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
1461 *   {
1462 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1463 *   {
1464 *  
1467 *   ADNumberType MijTjcurr = Cp*rho*fe_values[t].value(i,q_index)*old_solution[q_index];
1468 *   ADNumberType MijTjprev = Cp*rho*fe_values[t].value(i,q_index)*consol[q_index];
1469 *   ADNumberType k_curr = a + b*old_solution[q_index] + c*std::pow(old_solution[q_index],2);
1470 *   ADNumberType k_prev = a + b*consol[q_index] + c*std::pow(consol[q_index],2);
1471 *   ADNumberType Licurr = alpha * delta_t * (fe_values[t].gradient(i,q_index)*k_curr*old_solution_grad[q_index]);
1472 *   ADNumberType Liprev = (1-alpha) * delta_t * (fe_values[t].gradient(i,q_index)*k_prev*consol_grad[q_index]);
1473 *   residual_ad[i] += (MijTjcurr+Licurr-MijTjprev+Liprev)*fe_values.JxW(q_index);
1474 *   }
1475 *   }
1476 *  
1479 *  
1480 *   for (unsigned int face_number = 0;face_number<GeometryInfo<2>::faces_per_cell; ++face_number)
1481 *   {
1482 *  
1483 *   if (cell->face(face_number)->boundary_id() == 3)
1484 *   {
1485 *   fe_face_values.reinit(cell, face_number);
1486 *   for (unsigned int q_point=0;q_point<n_q_f_points;++q_point)
1487 *   {
1488 *   for (unsigned int i =0;i<dofs_per_cell;++i)
1489 *  
1492 *   residual_ad[i]+= -delta_t*(-10)*fe_face_values[t].value(i,q_point)*fe_face_values.JxW(q_point);
1493 *  
1494 *   }
1495 *   }
1496 *   }
1497 *  
1500 *   ad_helper.register_residual_vector(residual_ad);
1501 *  
1506 *   ad_helper.compute_residual(cell_rhs);
1507 *   cell->get_dof_indices(local_dof_indices);
1508 *  
1509 *   for (unsigned int i =0;i < dofs_per_cell; ++i)
1510 *   residual(local_dof_indices[i])+= cell_rhs(i);
1511 *  
1512 *   }
1513 *  
1514 *   for(const types::global_dof_index i: DoFTools::extract_boundary_dofs(dof_handler,t_mask,{1}))
1515 *   residual(i) = 0;
1516 *  
1517 *   std::cout << " The Norm is :: = " << residual.l2_norm() << std::endl;
1518 *   }
1519 * @endcode
1520
1521
1522<a name="ann-source/initial_conditions.cc"></a>
1524 *
1525 *
1526 *
1527 *
1528 * @code
1529 *   /* -----------------------------------------------------------------------------
1530 *   *
1531 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1532 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1533 *   *
1534 *   * This file is part of the deal.II code gallery.
1535 *   *
1536 *   * -----------------------------------------------------------------------------
1537 *   */
1538 *  
1539 *   #include "nonlinear_heat.h"
1540 *  
1541 *  
1547 *   double Initialcondition::value(const Point<2> & /*p*/, const unsigned int /*comp*/) const
1548 *   {
1549 *  
1552 *   return 0.0;
1553 *   }
1554 * @endcode
1555
1556
1557<a name="ann-source/nonlinear_heat_cons_des.cc"></a>
1559 *
1560 *
1561 *
1562 *
1563 * @code
1564 *   /* -----------------------------------------------------------------------------
1565 *   *
1566 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1567 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1568 *   *
1569 *   * This file is part of the deal.II code gallery.
1570 *   *
1571 *   * -----------------------------------------------------------------------------
1572 *   */
1573 *  
1574 *   #include "nonlinear_heat.h"
1575 *  
1576 *  
1581 *   :delta_t(0.1),
1582 *   alpha(0.5),
1583 *   tot_time(5),
1584 *   a(0.3),
1585 *   b(0.003),
1586 *   c(0),
1587 *   Cp(1),
1588 *   rho(1),
1589 *   dof_handler(triangulation),
1590 *   fe(FE_Q<2>(1), 1)
1591 *   {}
1592 *  
1593 *  
1597 *   {
1598 *   dof_handler.clear();
1599 *   }
1600 *  
1601 *  
1602 * @endcode
1603
1604
1605<a name="ann-source/output_results.cc"></a>
1607 *
1608 *
1609 *
1610 *
1611 * @code
1612 *   /* -----------------------------------------------------------------------------
1613 *   *
1614 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1615 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1616 *   *
1617 *   * This file is part of the deal.II code gallery.
1618 *   *
1619 *   * -----------------------------------------------------------------------------
1620 *   */
1621 *  
1622 *   #include "nonlinear_heat.h"
1623 *  
1624 *  
1628 *   void nonlinear_heat::output_results(unsigned int prn) const
1629 *   {
1630 *   DataOut<2> data_out;
1631 *   data_out.attach_dof_handler(dof_handler);
1632 *   std::vector<std::string> solution_names;
1633 *   solution_names.emplace_back ("Temperature");
1634 *   data_out.add_data_vector(converged_solution, solution_names);
1635 *   data_out.build_patches();
1636 *   const std::string filename =
1637 *   "output/solution-" + Utilities::int_to_string(prn , 3) + ".vtu";
1638 *   std::ofstream output(filename);
1639 *   data_out.write_vtu(output);
1640 *   }
1641 * @endcode
1642
1643
1644<a name="ann-source/set_boundary_conditions.cc"></a>
1646 *
1647 *
1648 *
1649 *
1650 * @code
1651 *   /* -----------------------------------------------------------------------------
1652 *   *
1653 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1654 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1655 *   *
1656 *   * This file is part of the deal.II code gallery.
1657 *   *
1658 *   * -----------------------------------------------------------------------------
1659 *   */
1660 *  
1661 *   #include "nonlinear_heat.h"
1662 *   #include "allheaders.h"
1663 *  
1664 *  
1669 *   {
1671 *   bl_left.set_time(time);
1672 *   std::map<types::global_dof_index,double> boundary_values;
1674 *   1,
1675 *   bl_left,
1677 *  
1678 *   for (auto &boundary_value: boundary_values)
1680 *   }
1681 * @endcode
1682
1683
1684<a name="ann-source/setup_system.cc"></a>
1686 *
1687 *
1688 *
1689 *
1690 * @code
1691 *   /* -----------------------------------------------------------------------------
1692 *   *
1693 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1694 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1695 *   *
1696 *   * This file is part of the deal.II code gallery.
1697 *   *
1698 *   * -----------------------------------------------------------------------------
1699 *   */
1700 *  
1701 *   #include "allheaders.h"
1702 *   #include "nonlinear_heat.h"
1703 *  
1704 *  
1708 *   void nonlinear_heat::setup_system(unsigned int time_step)
1709 *   {
1710 *   if (time_step ==0) {
1711 *   dof_handler.distribute_dofs(fe);
1712 *   converged_solution.reinit(dof_handler.n_dofs());
1713 *   present_solution.reinit(dof_handler.n_dofs());
1714 *   }
1715 *  
1716 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
1717 *   DoFTools::make_sparsity_pattern(dof_handler, dsp);
1718 *   sparsity_pattern.copy_from(dsp);
1719 *   system_matrix.reinit(sparsity_pattern);
1720 *   matrix_factorization.reset();
1721 *   }
1722 * @endcode
1723
1724
1725<a name="ann-source/solve_and_run.cc"></a>
1727 *
1728 *
1729 *
1730 *
1731 * @code
1732 *   /* -----------------------------------------------------------------------------
1733 *   *
1734 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1735 *   * Copyright (C) 2024 by Narasimhan Swaminathan
1736 *   *
1737 *   * This file is part of the deal.II code gallery.
1738 *   *
1739 *   * -----------------------------------------------------------------------------
1740 *   */
1741 *  
1742 *   #include "nonlinear_heat.h"
1743 *  
1744 *  
1749 *   void nonlinear_heat::solve(const Vector<double> & rhs, Vector<double> & solution, const double /*tolerance*/)
1750 *   {
1751 *   std::cout << " Solving linear system" << std::endl;
1752 *   matrix_factorization->vmult(solution, rhs);
1753 *   }
1754 *  
1755 *   void nonlinear_heat::run()
1756 *   {
1758 *   gridin.attach_triangulation(triangulation);
1759 *   std::ifstream f("mesh/mesh.msh");
1760 *   gridin.read_msh(f);
1761 *   triangulation.refine_global(1);
1762 *  
1763 *   double time = 0;
1764 *   unsigned int timestep_number = 0;
1765 *   unsigned int prn =0;
1766 *  
1767 *  
1768 *   while (time <=tot_time)
1769 *   {
1770 *   if(time ==0)
1771 *   {
1773 *  
1774 *   VectorTools::interpolate(dof_handler,
1775 *   Initialcondition(),
1777 *  
1778 *   VectorTools::interpolate(dof_handler,
1779 *   Initialcondition(),
1781 *   }
1782 *   else
1783 *   {
1784 *  
1788 *   }
1789 *   std::cout<<">>>>> Time now is: "<<time <<std::endl;
1790 *   std::cout<<">>>>> Time step is:"<<timestep_number<<std::endl;
1792 *   {
1793 *  
1794 *   const double target_tolerance = 1e-3;
1795 *  
1799 *   typename TrilinosWrappers::NOXSolver<Vector<double>>::AdditionalData additional_data;
1800 *   additional_data.abs_tol = target_tolerance;
1801 *   additional_data.max_iter = 100;
1803 *   additional_data);
1804 *  
1805 *  
1808 *   nonlinear_solver.residual =
1809 *   [&](const Vector<double> &evaluation_point,
1810 *   Vector<double> &residual) {
1811 *   compute_residual(evaluation_point, residual);
1812 *   };
1813 *  
1814 *  
1817 *  
1818 *   nonlinear_solver.setup_jacobian =
1819 *   [&](const Vector<double> &current_u) {
1820 *   compute_jacobian(current_u);
1821 *   };
1822 *  
1825 *  
1826 *   nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs,
1827 *   Vector<double> &dst,
1828 *   const double tolerance) {
1829 *   solve(rhs, dst, tolerance);
1830 *   };
1831 *  
1836 *   }
1837 *  
1841 *  
1842 *   if(timestep_number % 1 == 0) {
1843 *  
1845 *   prn = prn +1;
1846 *   }
1847 *   timestep_number++;
1848 *   time=time+delta_t;
1849 *  
1850 *   }
1851 *   }
1852 * @endcode
1853
1854
1855*/
Definition point.h:113
constexpr void clear()
friend class Tensor
Definition tensor.h:865
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
Definition tensor.h:851
Point< 2 > second
Definition grid_out.cc:4630
Point< 2 > first
Definition grid_out.cc:4629
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.
std::vector< index_type > data
Definition mpi.cc:740
const Event initial
Definition event.cc:68
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:616
constexpr char T
constexpr types::blas_int zero
constexpr types::blas_int one
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:193
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
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:474
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={}, const unsigned int level=numbers::invalid_unsigned_int)
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